# 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. Aug 13, 2012 at 12:38
• @J.M. Not in Mma 8.0.0.0 Aug 13, 2012 at 14:19
• @Verde, huh. I'm using 8.0.1.0 myself. I wonder... Aug 13, 2012 at 14:58
• Please write down your code for using LeastSquares and LyapunovSolve Aug 13, 2012 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]. Aug 14, 2012 at 0:29

If $$A$$ were symmetric, you could do it by using eigenvector representation of Lyapunov's equation and by replacing 1/0 divisions with 0

(* normalized right eigenvectors. *)
normalize[vec_] := vec/Norm[vec];
evecs[mat_] := Transpose[normalize /@ Eigenvectors[mat]];
norm2[mat_] := Total[Flatten[mat*mat]];

lyapSpectral[A_, B_] := Module[{U, s, C, cutoff, sdiv, Y},
(* A must be PSD *)
U = evecs[A]; (* A=U.D.U' *)
s = Eigenvalues[A];
C = U\[Transpose] . B . U;
cutoff = Max[s]*10^6*\$MachineEpsilon;
sdiv = Map[If[# > cutoff, 1/#, #] &, Outer[Plus, s, s], {2}];
Y = C*sdiv;
X = U . Y . Transpose[U];
Print["errorSpectral=", N@ norm2[A . X + X . A - B]];
X
];
A = {{10, 17}, {17, 29}};
B = {{26, 2}, {2, 4}};
lyapSpectral[A, B];


There's a similar eigenvector formula in non-symmetric case, but unclear if there's a way to transfer the symmetric trick to this case.

(* https://mathematica.stackexchange.com/questions/206675/making-work-as-matrixpower-for-matrices, need to use Notation palette *)
Needs["Notation"]
distinguishMatrixPowerAndPower[A_, d_] :=
If[SquareMatrixQ[A], MatrixPower[A, d], Power[A, d]];
Notation[
NotationTemplateTag[A_^d_] \[DoubleLongLeftRightArrow]
NotationTemplateTag[distinguishMatrixPowerAndPower[A_, d_]]]

lyapGeneric[A_, B_, C_] :=
Module[{n, ones, Avals, Avecs, Bvals, Bvecs, T, U, DA, DB},
n = First@Dimensions[A];
Assert[Dimensions[A] == {n, n}];
Assert[Dimensions[B] == {n, n}];

normalize[x_] := x/Norm[x];
ones := Array[1 &, {n, n}];
{Avals, Avecs} = Eigensystem[A];
{Bvals, Bvecs} = Eigensystem[B];
T = Transpose[Avecs];
U = Transpose[Bvecs];
DA = DiagonalMatrix[Avals];
DB = DiagonalMatrix[Bvals];
(*Assert[T.DA.Inverse[T]\[DotEqual]A]
Assert[U.DB.Inverse[U]\[DotEqual]B]*)
X = T . (T^-1 . C . U/(DA . ones + ones . DB)) . U^-1;
Print["errorGeneric=", N@ norm2[A . X + X . B - C]];
]
A = {{-1, 4}, {3, 0}};
B = {{-4, 3}, {2, 6}};
C0 = {{-5, 0}, {-2, -6}};
lyapGeneric[A, B, C0];

`