# LeastSquare Solution for the Continuous Time Lyapunov Equation

I have been working with a problem which involves solving the continuous time Lyapunov equation

$$A R + R A^\top = G$$

for the symmetric positive definite matrix $R$. Here $A$ is real, invertible and sparse and $G$ is diagonal but singular.

I am solving this equation for various matrices $A$ using LyapunovSolve. However, in some situations I get the error message "LyapunovSolve::nuniq: The matrix equation has no unique solution"

I can always transform this into an equation for each entry using CoefficientArrays and then solve with LeastSquares. This however, I have found to be incredibly slower than LyapunovSolve. Notwithstanding, it does yield a solution whose error (measured via the norm $\|A R + R A^\top - G\|$) is quite satisfactory.

My question is the following: is it possible to obtain with LyapunovSolve a solution which, albeit not unique, is the best possible solution (in the least square sense)? Below is an example matrix for reference (hopefully, copying it into a notebook should work)

G = DiagonalMatrix[SparseArray[{11 -> 2., 20 -> 1.}, {20}]]

A = SparseArray[{{1, 11} -> -1., {2, 12} -> -1., {3, 13} -> -1.,
{4, 14} -> -1., {5, 15} -> -1., {6, 16} -> -1.,
{7, 17} -> -1., {8, 18} -> -1., {9, 19} -> -1.,
{10, 20} -> -1., {11, 1} -> 10.5, {11, 2} -> -1.,
{11, 3} -> -1., {11, 4} -> -1., {11, 5} -> -1.,
{11, 6} -> -1., {11, 7} -> -1., {11, 8} -> -1.,
{11, 9} -> -1., {11, 11} -> 1., {12, 1} -> -1.,
{12, 2} -> 7., {12, 3} -> -1., {12, 4} -> -1.,
{12, 5} -> -1., {12, 6} -> -1., {12, 8} -> -1.,
{12, 10} -> -1., {13, 1} -> -1., {13, 2} -> -1.,
{13, 3} -> 8., {13, 4} -> -1., {13, 5} -> -1.,
{13, 6} -> -1., {13, 7} -> -1., {13, 9} -> -1.,
{13, 10} -> -1., {14, 1} -> -1., {14, 2} -> -1.,
{14, 3} -> -1., {14, 4} -> 6., {14, 5} -> -1.,
{14, 7} -> -1., {14, 9} -> -1., {15, 1} -> -1.,
{15, 2} -> -1., {15, 3} -> -1., {15, 4} -> -1.,
{15, 5} -> 6., {15, 6} -> -1., {15, 8} -> -1.,
{16, 1} -> -1., {16, 2} -> -1., {16, 3} -> -1.,
{16, 5} -> -1., {16, 6} -> 7., {16, 7} -> -1.,
{16, 8} -> -1., {16, 9} -> -1., {17, 1} -> -1.,
{17, 3} -> -1., {17, 4} -> -1., {17, 6} -> -1.,
{17, 7} -> 5., {17, 10} -> -1., {18, 1} -> -1.,
{18, 2} -> -1., {18, 5} -> -1., {18, 6} -> -1.,
{18, 8} -> 4., {19, 1} -> -1., {19, 3} -> -1.,
{19, 4} -> -1., {19, 6} -> -1., {19, 9} -> 5.,
{19, 10} -> -1., {20, 2} -> -1., {20, 3} -> -1.,
{20, 7} -> -1., {20, 9} -> -1., {20, 10} -> 5.5,
{20, 20} -> 1.}, {20, 20}]

-
At least for your example, LyapunovSolve[A, A // Transpose, G] works nicely. – J. M. Aug 13 '12 at 12:38
@J.M. Not in Mma 8.0.0.0 – Dr. belisarius Aug 13 '12 at 14:19
@Verde, huh. I'm using 8.0.1.0 myself. I wonder... – J. M. Aug 13 '12 at 14:58
Please write down your code for using LeastSquares and LyapunovSolve – Dr. belisarius Aug 13 '12 at 20:42
For least squares what I do is write a symbolic symmetric matrix in the sorts of R = Array[r,{n,n}]/.r[i_,j_]/;j<i->r[j,i]. Then I do {b,B} = CoefficientArrays[Flatten[A.R + R.A[Transpose] + G], vars] where vars = DeleteDuplicates@Flatten@R. Finally, I do LeastSquares[B,b]. – Gabriel Landi Aug 14 '12 at 0:29