9
$\begingroup$

I have a Hermitian matrix, and I would like to get a list of orthogonal eigenvectors and corresponding eigenvalues. But when there are degenerate eigenvalues, sometimes Eigensystem/Eigenvectors returns eigenvectors that are not orthogonal.

Now I could fix this by hunting down every eigenspace, and getting an orthogonal basis for it, but I feel like there must be a more efficient and less messy way.

Example:

eVectors =  Eigenvectors[{{1,-2,0,-2,0,0},{-2,1,0,-2,0,0},{0,0,1,0,-2,-2},
{-2,-2,0,1,0,0},{0,0,-2,0,1,-2},{0,0,-2,0,-2,1}}];

Now to find out whether they are orthogonal:

eVectors.ConjugateTranspose[eVectors]

which returns some nonzero off-diagonal elements.

$\endgroup$
  • 2
    $\begingroup$ According to the Documentation Center page for Eigenvectors, "For approximate numerical matrices m, the eigenvectors are normalized." and "For exact or symbolic matrices m, the eigenvectors are not normalized.". Thus to fix this, apply N to your matrix, and the result will be orthonormal, as you desired. $\endgroup$ – DumpsterDoofus Jan 31 '15 at 15:25
  • 7
    $\begingroup$ The reason for this behavior is simple: for exact input, the resulting eigenvectors can be complex numeric expressions, and calculating norms and other such arithmetic used in the Gram-Schmidt procedure becomes frightening from a computation-time perspective. As a result, Mathematica does not normalize symbolic eigenvectors because doing so could be catastrophic. In contrast, with floating-point results, such performance considerations do not apply, and so by default it orthonormalizes them. $\endgroup$ – DumpsterDoofus Jan 31 '15 at 15:28
  • 1
    $\begingroup$ @DumpsterDoofus Could you post it as an answer? Normalization and orthogonalization are different procedures so I don't think it is RTFM. May be a duplicate, but I can't find such Q&A. $\endgroup$ – ybeltukov Jan 31 '15 at 18:11
  • $\begingroup$ That's useful @DumpsterDoofus, thanks. I would like to be able to do this exactly in some cases, with a (real) parameter or two in the matrix, so I'm still interested in how best to go about it efficiently. $\endgroup$ – Holographer Jan 31 '15 at 18:15
  • $\begingroup$ @DumpsterDoofus I agree with ybeltukov, you should post that as an answer. Not only because it seems right but also because it shows a good way to reason about the inner workings of a black box. $\endgroup$ – Dr. belisarius Jan 31 '15 at 19:02
12
$\begingroup$

While the documentation does not specifically say that symbolic Hermitian matrices are not necessarily given orthonormal eigenbases, it does say

For approximate numerical matrices m, the eigenvectors are normalized.

For exact or symbolic matrices m, the eigenvectors are not normalized.

From this, it's reasonable to guess that if Mathematica isn't smart enough to normalize symbolic results, then it's also probably not smart enough to orthonormalize them, either.

Why is Mathematica not smart enough to do this? For exact input, the resulting eigenvectors can be complex numeric expressions, and calculating norms and other such arithmetic used in the Gram-Schmidt procedure becomes frightening from a computation-time perspective. As a result, Mathematica does not normalize symbolic eigenvectors because doing so could be catastrophic.

In contrast, with floating-point results, such performance considerations do not apply, and so by default it orthonormalizes them. Thus, if numeric results are acceptable, you can simply apply N to your array, and the results will automatically be orthonormal.

Otherwise, if numeric results are not acceptable, then you will have to manually orthogonalize the individual subspaces. Here is one method:

A = {{1, -2, 0, -2, 0, 0}, {-2, 1, 0, -2, 0, 0}, {0, 0, 1, 
    0, -2, -2}, {-2, -2, 0, 1, 0, 0}, {0, 0, -2, 0, 1, -2}, {0, 0, -2,
     0, -2, 1}};
{val, vec} = Eigensystem[A];
o = Flatten[Orthogonalize /@ FindClusters[N@val -> vec], 1];

It is orthonormal, and it diagonalizes A:

o.o\[Transpose] == IdentityMatrix[6]
Simplify[o.A.o\[Transpose]] == DiagonalMatrix[val]
(*True*)
(*True*)

While it works for your example matrix, it may fail for more complicated cases, so check your results.

As a further example of the catastrophic issues involved with orthogonalizing symbolic vectors, try executing the following:

Orthogonalize[Eigenvectors[#.#\[Transpose] &@RandomInteger[{-5, 5}, {6, 6}]]]

The calculation just goes on and on, because the eigenvectors are comprised of giant Root objects. So if symbolic results are what you need, you may run into trouble.

$\endgroup$
  • 4
    $\begingroup$ "it's reasonable to guess that if Mathematica isn't smart enough to normalize symbolic results, then it's also probably not smart enough..."?? I beg to differ. In point of fact, Mathematica is smart enough to NOT normalize those results. It would be a bad idea. Orthonormalizing them would compound the badness considerably. $\endgroup$ – Daniel Lichtblau Jan 31 '15 at 21:27
  • 1
    $\begingroup$ @DanielLichtblau: Yup, that was the main point of my answer (ie, that orthonormalizing symbolic results would be a disaster). I thought I made this point clear in the sentences that came afterwards, but maybe I should reword it? $\endgroup$ – DumpsterDoofus Jan 31 '15 at 21:36
  • 1
    $\begingroup$ Your comment is sufficient, I think. (If someone sees mine, they'll see yours. And if neither is seen, it's not an issue.) $\endgroup$ – Daniel Lichtblau Jan 31 '15 at 21:38
4
$\begingroup$

Here is a method that works when eigenvalues do not involve Root objects. Perturb symmetrically, and in such a way that equal eigenvalues become unequal (or enough do that we can get an orthogonal set of eigenvectors). Then take the limit as the perturbation goes to zero.

Here I add e to the (1,3) and (3,1) positions.

pmat = {{1, -2, e, -2, 0, 0}, {-2, 1, 0, -2, 0, 0}, {e, 0, 1, 0, -2, -2},
  {-2, -2, 0, 1, 0, 0}, {0, 0, -2, 0, 1, -2}, {0, 0, -2, 0, -2, 1}};

{pvals, pvecs} = Eigensystem[pmat];
vecs = Limit[pvecs, e -> 0]

(* Out[3691]= {{0, 0, 0, 0, -1, 1}, {0, -1, 0, 1, 0, 0}, {-1, -1, 1, -1, 1, 1},
  {2, -1, -2, -1, 1, 1}, {1, 1, 1, 1, 1, 1}, {-2, 1, -2, 1, 1, 1}} *)

Check orthogonality:

In[3692]:= vecs.Transpose[vecs]

(* Out[3692]= {{2, 0, 0, 0, 0, 0}, {0, 2, 0, 0, 0, 0}, {0, 0, 6, 0, 0, 0},
  {0, 0, 0, 12, 0, 0}, {0, 0, 0, 0, 6, 0}, {0, 0, 0, 0, 0, 12}} *)
$\endgroup$
1
$\begingroup$

The reason that that the orthogonalization fails is not because there is any problem performing Graham Schmidt. It is because no general symbolic solutions to polynomials beyond 5th order are known and computing eigenvalues of a 5x5 matrix is equivalent to solving a 5th order polynomial. (http://mathworld.wolfram.com/Polynomial.html ) .

The original problem had trouble getting orthogonal symbolic eigenvectors for the matrix

M = {{1, -2, 0, -2, 0, 0}, {-2, 1, 0, -2, 0, 0}, {0, 0, 1, 0, -2, -2}, {-2, -2, 0, 1, 0, 0}, {0, 0, -2, 0, 1, -2}, {0, 0, -2, 0, -2, 1}}

In spite of the fact that this is a 6x6 matrix, it's eigenvalues can be computed exactly and they are plus and minus 3.

 ev=Eigenvalues[M]

with that knowledge in hand it is trivial to compute an orthogonal basis for the degenerate eigenvalues. First note that that a basis can be computed by asking for the null space of M-lambda I for each lambda in the spectrum of M. In the case at hand the following two sets of eigenvectors span the two subspaces. The eigenvectors in one set are orthogonal to those in the other set, as they must be.

 evp = NullSpace[(M - 3 IdentityMatrix[6])]
 evm = NullSpace[(M + 3 IdentityMatrix[6])]

evp[[1]].evm[[1]]

Orthogonalization of the degenerate subspaces proceeds without difficulty as can be seen from the following.

orthogplus=Orthogonalize[evp]
orthogminux=Orthogonalize[evm]

all work quite cleanly.

If you can get the exact eigenvalues you will be able to obtain the exact eigenvectors and you'll be able to orthogonalize those that correspond to degenerate eigenvalues without difficulty. A general procedure for obtaining a symbolic orthogonal basis for any matrix M is:

M = {{1, -2, 0, -2, 0, 0}, {-2, 1, 0, -2, 0, 0}, {0, 0, 1, 
0, -2, -2}, {-2, -2, 0, 1, 0, 0}, {0, 0, -2, 0, 1, -2}, {0, 0, -2,
 0, -2, 1}};

{evals,evecs}=Eigensystem[M];

If the spectrum (evals) is not exact, STOP. There is no symbolic eigenbasis because it is not possible to compute the eigenvalues symbolically. I don't know how to tell Mathematica to tell if a spectrum is exact or not but it will be apparent by inspection. If the eigenvalues are exact they will either be rational or will involve symbolic expressions containing radicals.

orthogaleigenvectors=Orthogonalize[evecs];
$\endgroup$
  • $\begingroup$ Umm, I don't see what might be the impediment for 5x5 and beyond. The eigenvectors might be expressed in terms of Root objects. But that won't stop one from using Gram-Schmidt. It might make it quite slow perhaps, but that's a computational complexity issue and not a mathematical obstruction per se. $\endgroup$ – Daniel Lichtblau Feb 2 '15 at 18:11
  • $\begingroup$ There is a fundamental difference between the solutions to a 4th order and lower polynomials and 5th order and higher. The solution to 4th order and lower polynomials can be calculated using the operations of root extraction, addition, subtraction, multiplication, and division. It was proven by Abel in the early 19th century that 5th order and higher equations have no symbolic solution in terms of those operations. Later Hermite and others found solutions to the 5th and higher order polynomials in terms of special functions but as far as I know Mathematica cannot implement those solutions. $\endgroup$ – JEP Feb 3 '15 at 2:24
  • 1
    $\begingroup$ No. I'm claiming that Gram-Schmidt will work regardless of whether fed vectors comprised of radicals or (as I explicitly noted) Root objects. So this is not a case where "orthogonalization fails". $\endgroup$ – Daniel Lichtblau Feb 3 '15 at 16:18
  • 1
    $\begingroup$ No, it's slow-to-forever because of intermediate swell. This is going to happen whether one has Root objects or radicals. For example, try Orthogonalize on the eigenvectors from a random 4x4 matrix of integers. You will find it to be slow with either setting of Quartics e.g. {vals, vecs} = Eigensystem[mat, Quartics -> True]; Timing[v2 = Orthogonalize[vecs];] (were True to work much better, it would be the default setting.) This has little to do with the form of the root solutions per se. It has to do with the computational complexity in the algebraic manipulations thereof. $\endgroup$ – Daniel Lichtblau Feb 3 '15 at 23:48
  • 1
    $\begingroup$ I think I'm going to defer at this point to your expertise. My own understanding of the underlying math and computation is mostly anecdotal. From working with this stuff day in and out since the Bush 41 years. $\endgroup$ – Daniel Lichtblau Feb 4 '15 at 16:09

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.