# What is LyapunovSolve::meig message about?

According to Wikipedia, a Sylvester equation $AX+XB=C$ has a unique solution for $X$ exactly when there are no common eigenvalues of $A$ and $-B$. But the following code gives me a message I don't quite understand:

A = RandomReal[{-10, 10}, {3, 3}];
B = Table[5 KroneckerDelta[i, j], {i, 3}, {j, 3}];
c = IdentityMatrix[3];
LyapunovSolve[A, B, c];


LyapunovSolve::meig: Matrix with multiple eigenvalues encountered. The solution may be inaccurate or may not exist.

Here A and -B will have no common eigenvalue. So I suppose the message refers to something other. I've tried looking for mentioning of eigenvalues in the documentation for LyapunovSolve (where LyapunovSolve::meig points unhelpfully anyway), but there the stability for $ax+xa^\dagger=c$ is considered only.

I've checked that the message vanishes if I break degeneracy of spectrum of B by adding RandomReal[] to its elements.

So what does this message refer to? Are there some special stability problems with this equation, not covered in the documentation?

• "The solution may be inaccurate" - and this is the case. If X = LyapunovSolve[A, B, c], then A.X + X.B will give an identity matrix with not exact zeros, but sth like 10^-16. Contrary, if you simply take A = RandomInteger[....], there will be no message. Rationalizing A helps. – corey979 Dec 7 '16 at 20:48
• I think it's complaining about the repeated eigenvalues of B. I don't know the underlying algorithm, but it must be a case in which the algorithm is sometimes numerically unstable. It's a warning that you should check the answer (which seems pretty good). – Michael E2 Dec 7 '16 at 22:22
• @corey979 I disagree. Consider LyapunovSolve[N@DiagonalMatrix[{1, 2, 3}], DiagonalMatrix[{4, 4, 5}], IdentityMatrix[3]]. It's complaining about a potential numerics issue. OTOH, an exact A, B, and c invokes an exact solver, which has no numerics issues. – Michael E2 Dec 7 '16 at 23:12
• @MichaelE2 That essentially illustrates my point: it complains when given approximate numbers. Changing to exact (via Rationalize, e.g.) allows to give an exact answer, so it does not complain. Or sth. – corey979 Dec 7 '16 at 23:16
• @corey979 It does not always complain when given approximate numbers. It complains when at least one matrix has inexact numbers and B has repeated eigenvalues. My point is similar (or complementary) to yours: M is warning, I surmise, that a numerical routine is being used that is known to be unstable on matrices with repeated eigenvalues. – Michael E2 Dec 7 '16 at 23:38

Classically, a Sylvester/Lyapunov equation is solved via the Bartels-Stewart algorithm. In your given example, however, the algorithm (suitably simplified in this case, since your $\mathbf B$ is diagonal) works without a hitch, so I suspect something else is being done behind the scenes:

BlockRandom[SeedRandom[42]; A = RandomReal[{-10, 10}, {3, 3}]];
B = 5 IdentityMatrix[3];
c = IdentityMatrix[3];

{q, t} = SchurDecomposition[A];
q.Transpose[MapThread[LinearSolve[t + #1 IdentityMatrix[3], #2] &,
{Diagonal[B], Transpose[c].q}]]
{{0.0828636, -0.113095, -0.148619}, {-0.098646, 0.0380146, -0.0940719},
{-0.161339, -0.157174, -0.103911}}


which is virtually the same as the solution obtained from LyapunovSolve[].

If one wants a further check, one can do the (usually wasteful) reformulation of the Sylvester equation into a linear system involving Kronecker products:

Transpose[Partition[LinearSolve[KroneckerProduct[IdentityMatrix[Length[B]], A] +
KroneckerProduct[Transpose[B], IdentityMatrix[Length[A]]],
Flatten[Transpose[c]]], 3]]


which yields the same matrix up to roundoff error.

A = RandomReal[{-10, 10}, {3, 3}];
B = Table[5 KroneckerDelta[i, j], {i, 3}, {j, 3}];
c = IdentityMatrix[3];
X = LyapunovSolve[A, B, c]


indeed gives the cited error, but

MatrixForm[A.X + X.B]


explains that the issue is that "the solution may be inaccurate". Contrarily, changing RandomReal to RandomInteger leads to no errors/warnings (showing that the repeated eigenvalues of B are not a problem here - the relevant theorem says that the eigenvalues of A and -B should be different):

A = RandomInteger[{-10, 10}, {3, 3}];
B = Table[5 KroneckerDelta[i, j], {i, 3}, {j, 3}];
c = IdentityMatrix[3];
X = LyapunovSolve[A, B, c]

MatrixForm[A.X + X.B]


Going back to RandomReal, if one Rationalizes the matrix A:

A = Rationalize[#, 10^-16] & @ RandomReal[{-10, 10}, {3, 3}];
B = Table[5 KroneckerDelta[i, j], {i, 3}, {j, 3}];
c = IdentityMatrix[3];
X = LyapunovSolve[A, B, c]


it gives no warning, and

MatrixForm[A.X + X.B]


is also

For comparison:

MatrixForm[N @ A.X + X.B]


Hence, LyapunovSolve seems to have some issues with the underlying numerics. With exact values, it gives an exact answer and does not complain.

• Note that A.X + X.B is accurate to an error of 8 \$MachineEpsilon * Norm[c], i.e. to machine precision less 3 bits error, which is not at all bad. – Michael E2 Dec 7 '16 at 23:14