# Recurrence relation for matrices

I have matrices ($$S_0$$ thought $$S_N$$) and I have a recurrence relation that link successive matrices together. $$S_i + S_i(aS_{i-1}^{-1})S_i=C_i+aS_{i+1}$$

We can assume for this problem that $$S_0=S_N=I$$, the identity matrix. The $$C_i$$ in this equation are all known beforehand. All matrices $$S_i$$ are definite positive and $$a>0$$. The matrices $$C_i$$ are positive semi-definite, with some eigenvalues equal to 0 but at least one eigenvalue is positive.

If I look at the problem with $$1\times1$$ matrices, the analytic solution for $$S_t$$ in terms of $$S_{t-1}$$ and $$S_{t+1}$$ is quite straightforward to write and there is always exactly one positive root. I would like to see if I can find an iterative scheme to compute the matrices $$S_i$$, even if it requires many iterations and it takes forever...

Edit: The simplified question is to compute $$X$$ from $$X+XAX=B$$, where $$A$$ and $$B$$ are known. All matrices here are symmetrical and definite-positive.

Edit 2: Looking at the equation $$X+XAX=B$$, it is very similar to a quadratic equation. So I thought that maybe we can solve it in a similar way.

We can factor the equation as: \begin{align}XAX+X&=B\\ \left(A^\frac12X+\frac12A^{-\frac12}\right)^T\left(A^\frac12X+\frac12A^{-\frac12}\right)&=\frac14A^{-1}+B\\ &=\frac{A^{-1}}4\left(I+4AB\right)\end{align} The matrix $$A$$ is positive-definite so $$A^\frac12$$, $$A^{-\frac12}$$, and $$A^{-1}$$ all exist and also are positive-definite. That means that $$I+4AB$$ is also positive-definite. So we can "extract the square root" and rewrite: \begin{align}\left(A^\frac12X+\frac12A^{-\frac12}\right)&=\frac{A^{-\frac12}}2\left(I+4AB\right)^\frac12\\ X+\frac12A^{-1}&=\frac{A^{-1}}2\left(I+4AB\right)^\frac12\\ X&=\frac{A^{-1}}2\left[-I+\left(I+4AB\right)^\frac12\right] \end{align}

This is basically the same approach as the one used to solve a quadratic equation, except that I solve a matrix equation. Also, it's clear that the positive root only is used, so we can retrieve a matrix $$X$$ that is also definite-positive.

Updated question: Is this approach valid? I know I shouldn't do an updated question but I will copy it in an answer only if it's valid...

• You can check if the approach is valid by sticking the solution into the equation and checking if it's satisfied. It seems to be valid ... so this seems to simpler than solving generic Riccati eq -- colab.research.google.com/drive/… Commented Jun 29, 2022 at 21:36
• @YaroslavBulatov Indeed, it's probably simpler this way. Will check and compare. Thanks for the advice.
– PC1
Commented Jun 30, 2022 at 0:22
• @YaroslavBulatov I tried the ScyPi solver for CARE equations (scipy.linalg.solve_continuous_are) and it seems about 8x slower than the matrix version in my post. The results are very similar.
– PC1
Commented Jun 30, 2022 at 0:59
• @YaroslavBulatov Interesting, thanks for sharing the experiments. I'll keep my comment there though, just to point out that this formula might need a different proof that it works. :) Commented Jun 30, 2022 at 16:08
• @FedericoPoloni I would prefer to solve for al $S_i$ matrices at once but I have trouble to solve simultaneously. So I just iterate, computing $S'_i$ from $S_i$ over and over... Not efficient but it seems to converge to some sensible result. If there are more practical ways to approach the problem then I would be happy to see but I didn't find much information yet.
– PC1
Commented Jun 30, 2022 at 16:48

Another different way to obtain an equivalent formula: let $$Y = A^{1/2}XA^{1/2}$$. Then, multiplying your equation from both sides by $$A^{1/2}$$, we have $$Y^2+Y = A^{1/2}BA^{1/2} = C$$.

The matrix $$C$$ is symmetric and positive definite, so we can diagonalize it $$C = U \operatorname{diag}(\lambda_1,\dots,\lambda_n)U^*$$ and obtain the positive definite solution as $$Y = U\operatorname{diag}(\mu_1,\dots,\mu_n)n)U^*$$, where $$\mu_i>0$$ solves the scalar equation $$\mu_i^2+\mu_i = \lambda_i$$. Then $$X = A^{-1/2}YA^{-1/2}$$.

• Indeed, this is an excellent way to get the result. I have the feeling that it might be more stable than my own try as there is no matrix inversion required.
– PC1
Commented Jul 1, 2022 at 18:29
• There is still a matrix inversion to get $X$ from $Y$, unfortunately. Commented Jul 2, 2022 at 7:29
• I think that I could also do a Cholesky decomposition. Let's say $A=LL^T$ so $Y=L^TXL$, the rest of the argument is valid.
– PC1
Commented Jul 9, 2022 at 3:42
• @PC1 Yes, that works too! Commented Jul 9, 2022 at 16:52

This is an instance of the Riccati equation, which can be solved using the Matrix Sign function. Relevant section from Higham's "Functions of Matrices" book: