1
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ "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. $\endgroup$ – corey979 Dec 7 '16 at 20:48
  • $\begingroup$ 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). $\endgroup$ – Michael E2 Dec 7 '16 at 22:22
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Dec 7 '16 at 23:12
  • $\begingroup$ @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. $\endgroup$ – corey979 Dec 7 '16 at 23:16
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Dec 7 '16 at 23:38
2
$\begingroup$

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.

$\endgroup$
1
$\begingroup$
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]

result

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]

result of exact case

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

identity matrix

For comparison:

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

numerical result

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

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – Michael E2 Dec 7 '16 at 23:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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