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?

  • $\begingroup$ Are you trying to solve AX+XA=A? $\endgroup$ Sep 13, 2019 at 19:54
  • $\begingroup$ Trying to solve AX+XA=C where C is allowed to be A (but can also be other values) $\endgroup$ Sep 14, 2019 at 2:23
  • $\begingroup$ 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 $\endgroup$ Sep 14, 2019 at 2:25

1 Answer 1


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

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

  • $\begingroup$ You might want to include the definition of erank for completeness. $\endgroup$ Sep 15, 2019 at 9:33
  • $\begingroup$ Oh that's just the trace divide by norm, was left over from experiments, fixed $\endgroup$ Sep 15, 2019 at 20:45

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.