# Solving $A=B+AB$ without matrix inverse

I have a linear system of equations that can be expressed $$A=B+AB,$$ where $A$ and $B$ are real, symmetric matrices. I would like to solve for $A$ given $B$. At present, I solve for $A$ directly via $$A=B(I-B)^{-1},$$ with the inverse computed via $LU$ factorization. For most relevant inputs $B$, the inversion goes well, but of course, sometimes $(I-B)$ is simply ill-conditioned. Is there a way to solve for $A$ without explicitly computing a matrix inverse?

• In case of matrix being ill conditioned/singular, you can always add some random noise and try to solve it. Or else, there are matrix conditioning algorithms. Finally, you could take QR decomposition of SVD methods, which are found to be numerically better than LU. – Tolga Birdal Apr 25 '16 at 19:40

If $A$ and $B$ are real symmetric, then $A=B+AB$ if and only if the product $AB$ is also real symmetric. In turn, $AB=BA$ holds if any only if $A$ and $B$ share a common eigendecomposition. This latter statement gives a recipe for computing $A$ given $B$.

1. Given $B$, compute its eigendecomposition $B=V\Lambda V^T$, in which $V$ is the orthonormal set of of eigenvectors, and $\Lambda = \mathrm{diag}(\lambda_1,\dots,\lambda_n)$ is the diagonal matrix of real eigenvalues.
2. Compute the diagonal matrix $S= \mathrm{diag}(s_1,\dots,s_n)$ in which each $s_k$ is defined with respect to the corresponding $\lambda_k$ as in $$s_k=\frac{\lambda_k}{1-\lambda_k}\text{ for all }k\in\{1,\ldots,n\}.$$
3. Output $A=VSV^T$.

In fact, this is the most numerically stable procedure possible for solving the desired equation. The computed $A$ is exactly real symmetric by construction. Also, limiting all matrix-matrix operations to orthogonal transformations minimizes the numerical error that are propagated from the problem data to the computed solution.

Of course, we will none-the-less encounter ill-conditioning if the eigenvalues of $B$ are close to 1. But this latter phenomenon is physical, and reflects the underlying ill-conditioning of the equation to be solved. In this latter case, we may reformulate the problem and solve it in a least-squares sense, using essentially the same steps as the ones presented above.

• You should give some reasoning why this procedure is the "most numerically stable" for solving this equations. – Jan Apr 26 '16 at 10:42
• Thanks for the suggestion @Jan. I've added some further details. – Richard Zhang Apr 26 '16 at 14:36

Your equation is a Sylvester equation $$AB - A = -B,$$ for which there are numerous solution algorithms.

I would definitely recommend the Bartels-Stewart algorithm that bases on real Schur decompositions. These can be computed for general square matrices in a stable fashion and avoiding complex arithmetics; cf. Chapter 16.1 in the book Accuracy and Stability of Numerical Algorithms by N. Higham.

However, all algorithms assume that the equation is uniquely solvable, which in your case is equivalent to that $B$ does not have $1$ as eigenvalue. For this case, I have never come across any theory that may have helped to build a numerical algorithm on it. (And I have done quite some research on that)

• Note that the Schur decomposition of a real symmetric matrix is just its usual orthonormal eigendecomposition. So if you apply the Bartels-Stewart algorithm to this problem, you will simply end up with the procedure I had suggested. – Richard Zhang Apr 26 '16 at 13:50
• Alright, I didn't notice that there are only real symmetric matrices involved. – Jan Apr 27 '16 at 6:40