3
$\begingroup$

I have a real symmetric matrix m (below) which I would like to decompose as $Q^TDQ$ with $Q$ orthogonal and $D$ diagonal. When I try to diagonalize it using Eigensystem as below, however, the result is block diagonal:

m - Transpose[m] // Norm

returns 0 but

{vals, vecs} = Eigensystem[m];
vecs . m . Transpose[vecs] // Simplify // MatrixForm

returns a block diagonal matrix. Why not a simple diagonal matrix? (The eigenvectors returned aren't unit length, but normalizing them doesn't change make the result diagonal.)

MWE:

m := {{-8, 1, 1, 1, 4, -2, 1, -2, 4}, {1, 4, -2, 4, 
    1, -2, -2, -2, -2}, {1, -2, 4, -2, -2, -2, 4, -2, 1}, {1, 4, -2, 
    4, 1, -2, -2, -2, -2}, {4, 1, -2, 1, -8, 1, -2, 1, 
    4}, {-2, -2, -2, -2, 1, 4, -2, 4, 1}, {1, -2, 4, -2, -2, -2, 
    4, -2, 1}, {-2, -2, -2, -2, 1, 4, -2, 4, 1}, {4, -2, 1, -2, 4, 1, 
    1, 1, -8}};
Norm[m - Transpose[m]]
{vals, vecs2} = Eigensystem[m];
Transpose[vecs] . DiagonalMatrix[vals] . vecs // Simplify // 
  MatrixForm;
vecs2 . m . Transpose[vecs2] // Simplify // MatrixForm

There are repeated eigenvalues, and I guess Mathematica may choose a nonorthogonal basis for the eigenspaces of dimension >1. (I say "may" since sometimes it does choose an orthogonal basis, e.g., running Eigensystem on the identity matrix). That's reasonable since Eigensystem doesn't apply just to symmetric matrices and an orthogonal eigenbasis may not exist on its input. If I use SingularValueDecomposition it returns different left and right eigenvectors ie $Q_1^TDQ_2$. Neither routine takes into account that the matrix is symmetric. Is there not a built-in routine for diagonalizing a symmetric matrix a la standard spectral decomposition?

Improve this question
$\endgroup$
5
  • $\begingroup$ Looks like a bug to me. According to the documentation, "All the nonzero eigenvectors given are independent." However, vecs[[1]] . vecs[[2]] gives $6 \sqrt{2}$, not zero. So something is wrong here. The problem disappears when using a numerical matrix: {vals, vecs} = Eigensystem[N[m]]. $\endgroup$
    – Roman
    Commented Aug 2, 2021 at 7:22
  • 2
    $\begingroup$ @Roman Independent != orthogonal. This matrix has multiplicity in the eigenvalues. And eigenvectors corresponding to the same eigenvalue need not be orthogonal. $\endgroup$
    – Daniel Lichtblau
    Commented Aug 2, 2021 at 13:38
  • $\begingroup$ @DanielLichtblau thanks, I missed that distinction. $\endgroup$
    – Roman
    Commented Aug 2, 2021 at 14:17
  • $\begingroup$ @DanielLichtblau But they can be chosen to be orthogonal for a symmetric matrix like this. So I guess my question is how to do a spectral decomposition in mathematica, which I've updated. $\endgroup$
    – sayda
    Commented Aug 2, 2021 at 16:27
  • $\begingroup$ You can use Orthogonalize on the eigenvectors. $\endgroup$
    – Daniel Lichtblau
    Commented Aug 2, 2021 at 22:56

1 Answer 1

Reset to default
5
$\begingroup$

Just ortho-normalize the eigenvectors with the Gram–Schmidt process:

m = {{-8,  1,  1,  1,  4, -2,  1, -2,  4},
     { 1,  4, -2,  4,  1, -2, -2, -2, -2},
     { 1, -2,  4, -2, -2, -2,  4, -2,  1},
     { 1,  4, -2,  4,  1, -2, -2, -2, -2},
     { 4,  1, -2,  1, -8,  1, -2,  1,  4},
     {-2, -2, -2, -2,  1,  4, -2,  4,  1},
     { 1, -2,  4, -2, -2, -2,  4, -2,  1},
     {-2, -2, -2, -2,  1,  4, -2,  4,  1},
     { 4, -2,  1, -2,  4,  1,  1,  1, -8}};

{vals, vecs} = Eigensystem[m];
nvecs = Orthogonalize[vecs, Method -> "GramSchmidt"];

nvecs . m . Transpose[nvecs] // FullSimplify // MatrixForm

$$ \left( \begin{array}{ccccccccc} -9 \sqrt{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -9 \sqrt{2} & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 9 \sqrt{2} & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 9 \sqrt{2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{array} \right) $$

Additionally, if you want simpler expressions for the ortho-normalized eigenvectors, RootReduce can help:

nvecs // RootReduce
(*    simplified ortho-normal eigenvectors    *)
Improve this answer
$\endgroup$
4
  • $\begingroup$ It's far from obvious that this will work. The output of Orthogonalize is only correct if it doesn't mix eigenvectors of different eigenvalues, but it has no access to the eigenvalues. I think it does work, but only because Orthogonalize uses Gram-Schmidt by default, and the computation is exact/symbolic. $\endgroup$
    – benrg
    Commented Mar 24, 2022 at 0:06
  • $\begingroup$ @benrg For symmetric matrices, eigenvectors associated with different eigenvectors are always orthogonal. Consequently there will never be mixing of eigenvectors of different eigenvalues. $\endgroup$
    – Roman
    Commented Mar 24, 2022 at 6:09
  • $\begingroup$ At a high level, Orthogonalize just returns an orthonormal basis for the space spanned by the inputs. If the input is an invertible square matrix, as it may be here, then the identity matrix is a valid return value. To prove it doesn't return the identity matrix, you have to know what algorithm it uses. $\endgroup$
    – benrg
    Commented Mar 24, 2022 at 6:18
  • $\begingroup$ @benrg you are correct, of course; I was being unclear about the implicit default method of orthogonalization. For clarity and future compatibility I've made this default explicit now. Thanks! $\endgroup$
    – Roman
    Commented Oct 26, 2022 at 7:35

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.