Solving underdetermined Lyapunov equations?

I'm wondering if there's an efficient way to get a solution (ie, LeastSquaressolution) for Lyapunov equation $$AX+XA=C$$ with symmetric positive definite $$A$$ and $$C$$.

I want something that would work like LyapunovSolve, but would work for underconstrained problems, i.e., LyapunovSolve[A, A] should give me something whose spectrum looks like $$I$$.

I tried a naive approach which is to do Kronecker expansion followed by LeastSquares, which gives the desired result

kronExpand[x_] := Module[{ii},
ii = IdentityMatrix[First[Dimensions[x]]];
ii\[TensorProduct]x + Transpose[x]\[TensorProduct]ii
];
lyapLeastSquares[A_, B_] := Module[{d, X},
X = LeastSquares[kronExpand[A], vec[B]];
X = unvec[X, d];

(check this notebook for an end-to-end example)

However, this expansion is too large to be practical. IE, my matrices are on the order of 1000 which is fast using LyapunovSolve, but doing Kronecker expansion means that I have matrices on the order of 1M-by-1M. Any suggestions on how to make this feasible?

• Are you trying to solve AX+XA=A? – Mirko Aveta Sep 13 at 19:54
• Trying to solve AX+XA=C where C is allowed to be A (but can also be other values) – Yaroslav Bulatov Sep 14 at 2:23
• It seems in principle one could borrow the logic from Matlab Lyapunov solver and skip all divisions by zero (which result from system being underdetermined) -- github.com/msubhransu/matrix-sqrt/blob/master/lyap2.m – Yaroslav Bulatov Sep 14 at 2:25

I ended up with two versions below by reverse engineering Matlab implementation. They both seem to give me the same answer as the least squares solution when A, B are positive semi-positive (experiments notebook)

lyapSvd[A_, B_] := Module[{U, s, C, cutoff, sdiv, Y},
{U, S, V} = SingularValueDecomposition[A];
s = Diagonal[S];
C = U\[Transpose].B.V;
smat = Outer[Plus, s, s];
sinv = Map[If[# > 0, 1/#, #] &, smat, {2}];
Y = C*sinv;
X = U.Y.Transpose[V];
Print["errorSvd=", Norm[X.A + A.X - B]];
X
];

lyapSpectral[A_, B_] := Module[{U, s, C, cutoff, sdiv, Y},
U = Transpose@Eigenvectors[A];
s = Eigenvalues[A];
C = U\[Transpose].B.U;
cutoff = Max[s]*\$MachineEpsilon;
sdiv = Map[If[# > cutoff, 1/#, #] &, Outer[Plus, s, s], {2}];
Y = C*sdiv;
X = U.Y.Transpose[U];
Print["errorSpectral=", Norm[A.X + X.A - B]];
X
];

• You might want to include the definition of erank for completeness. – J. M. will be back soon Sep 15 at 9:33
• Oh that's just the trace divide by norm, was left over from experiments, fixed – Yaroslav Bulatov Sep 15 at 20:45