# spectral decomposition of matrix with repeated eigenvalues

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?

• Looks like a bug to me. According to the documentation, "All the nonzero eigenvectors given are independent." However, vecs[] . vecs[] gives $6 \sqrt{2}$, not zero. So something is wrong here. The problem disappears when using a numerical matrix: {vals, vecs} = Eigensystem[N[m]]. Aug 2, 2021 at 7:22
• @Roman Independent != orthogonal. This matrix has multiplicity in the eigenvalues. And eigenvectors corresponding to the same eigenvalue need not be orthogonal. Aug 2, 2021 at 13:38
• @DanielLichtblau thanks, I missed that distinction. Aug 2, 2021 at 14:17
• @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. Aug 2, 2021 at 16:27
• You can use Orthogonalize on the eigenvectors. Aug 2, 2021 at 22:56

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    *)

• 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. Mar 24, 2022 at 0:06
• @benrg For symmetric matrices, eigenvectors associated with different eigenvectors are always orthogonal. Consequently there will never be mixing of eigenvectors of different eigenvalues. Mar 24, 2022 at 6:09
• 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. Mar 24, 2022 at 6:18
• @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! Oct 26, 2022 at 7:35