How to obtain matrix from summation inverse equation I have a set of square matrices $\{A_i\}_{i \in \{1,..., n\}}$ and another square matrix of equal size $K$. Under the assumption that such a matrix exists and is unique, I want to find the unique $B$ that satisfies the following equation:
$$\sum_{i=1}^n(A_i+B)^{-1}=K$$
All matrices are of full rank and symmetric. Does anyone know of a way, algebraically or computationally, for me to work out the value of $B$?
 A: This does not seem a simple matrix equation to solve. Computationally, my first attempt would be with Newton's method, even if it takes $O(k^6)$ per iteration, where $k$ is the size of the matrices. The Jacobian of the map in the LHS is
$$
L_B f[H] = \sum_{i=1}^n (A_i+B)^{-1}H(A_i+B)^{-1},
$$
and to solve the equation $L_Bf[H] = Y$ for $H$ given $Y$ you need to convert it to a $k^2 \times k^2$ linear system (there are faster algorithms to solve this linear matrix equation for $n=2$, but I do not think there is anything better otherwise).
If you need to solve it for dimensions for which this is unfeasible, then I would try turning it into a fixed-point equation, for instance
$$
B = \left(K - \sum_{i=2}^n (A_i+B)^{-1}\right)^{-1} - A_1,
$$
and then hope that the iteration
$$
B_{m+1} = \left(K - \sum_{i=2}^n (A_i+B_m)^{-1}\right)^{-1} - A_1,
$$
converges.
The scalar version of this equation is a secular equation, but searching for this term I found nothing interesting for matrix arguments.
A: If $n\geq 2$, generically, you are unlikely to encounter such an equation with only one solution.
Let $d$ be the dimension of our symmetric matrices.
$\bullet$ Consider the case $n=2$. Then the algebraic resolution of the system leads to obtaining the roots of a real polynomial $P$ of degree $2^d$ (generically). Our system admits a unique solution iff $P$ admits only one real root. Since $degree(P)$ is even, this is not possible.
Then we must assume that $n\geq 3$.
$\bullet$ Note that a random polynomial of degree $p$ has
$O(\log(p))$ real zeros on average -when the coefficients are independent standard normals-
or $O(\sqrt{p})$ -when the variance varies with the index of the coefficient-.
If $A_1,A_2,K$ are randomly chosen, then $P$ can be roughly considered as random and the probability that $P$ admits a single real root is very small when $p$ is large.
$\bullet$ For example, if $d=2,n=3$, we obtain -generically- a polynomial of degree $11$. I did some tests and got $5,7,9$ or $11$ real roots.
Finally, the existence of a unique solution is a very special case.
