6
$\begingroup$

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)

Thank you in advance.

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}]
$\endgroup$
7
  • $\begingroup$ At least for your example, LyapunovSolve[A, A // Transpose, G] works nicely. $\endgroup$ Commented Aug 13, 2012 at 12:38
  • $\begingroup$ @J.M. Not in Mma 8.0.0.0 $\endgroup$ Commented Aug 13, 2012 at 14:19
  • $\begingroup$ @Verde, huh. I'm using 8.0.1.0 myself. I wonder... $\endgroup$ Commented Aug 13, 2012 at 14:58
  • $\begingroup$ Please write down your code for using LeastSquares and LyapunovSolve $\endgroup$ Commented Aug 13, 2012 at 20:42
  • $\begingroup$ 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]. $\endgroup$ Commented Aug 14, 2012 at 0:29

1 Answer 1

1
$\begingroup$

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];

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.