2
$\begingroup$

I use Eigen to compute the eigenvalues of symmetric matrices.

The problem is, that sometimes the matrices not nice at all numerically. Because of this, I get NaN among the eigenvalues. I have tested the matrix with maxima and WolframAlpha as well. Maxima would give me two complex and one real eigenvalues, while apparently WolframAlpha is able to compute all three real eigenvalues.

Is there a method or library, that could guarantee the eigenvalues? If not, what can be done in these cases?

For the curious, this is the matrix I've tested: $\begin{bmatrix} 1.4580680&-0.0028459986&0.10822548 \\\\ -0.0028459986&5.5612186 \times 10^{-6}&-0.00012824166 \\\\ 0.10822548&-0.00012824166&1.3172486 \end{bmatrix}$

  • 0
    If you are specifically interested in 3x3 matrices (it's mentioned in the question title, not in the body) the problem reduces to solving a cubic polynomial, for which there are 'exact' formulas. http://mathworld.wolfram.com/CubicFormula.html2011-11-18

3 Answers 3

0

The problem was twofold:

The first part is of some hard-to-put-a-finger-on-it compiler issue. On a different machine, Eigen yielded real eigenvalues for same matrix, as it should.

The second part is about numerical precision. On the other machine, I get real eigenvalues either using floats or doubles, but the numbers only match with the ones from the other answers, if doubles are used.

1

Using Pari/Gp I get three real eigenvalues $\small A = P \cdot D \cdot P^{-1} $ and with the eigenvectors normed I get by
$\small \begin{array} {lll} P = \text{mateigen}(A) \\ P = P \cdot \text{matdiagonal}(P[1,]^{-1.0}) & \text{ // normalize the matrix P} \\ Q= P^{-1} \\ D= Q \cdot A \cdot P \end{array}$

the tree matrices

$\small \begin{array} {r|rrr} & 1.00000000000 & 1.00000000000 & 1.00000000000 \\ P=& -0.00192221444287 & -0.00207348773146 & 511.089920195 \\ & 0.542410157237 & -1.84363064066 & -0.0324027372303 \\ \hline \\ \text{diag}(D)=& 1.51677607024 & 1.25854609012 & 0.000000000860027279267 \\ \hline \\ & 0.772670655253 & -0.00148523869311 & 0.419104411609 \\ Q=P^{-1}= & 0.227325516467 & -0.000471356669442 & -0.419104287562 \\ & 0.00000382828008387 & 0.00195659536255 & -0.000000124046753602 \end{array} $

If A is really symmetric, then there is also a solution having P orthogonal, meaning it is a pure rotation(matrix).

Using Pari/GP you can even get eigenvalues of matrices over the complex numbers, (with arbitrary precision) so I think it could serve as recommendation...(As you are asking for a lib (I think library for use in Fortran/C/et al programming languages) I think you can even use the API-library of Pari)

1

I get all the real eigenvalues using SymPy, and it gives me three real eigenvalues for the matrix you gave. Unfortunately though, it cannot determine the eigenvector for $\lambda = 1.51677607024351$, but it does determine all three real eigenvalues.

Specifically, I obtain that $\lambda \in \{ 8.60027278426881*10^{-10}, 1.25854609011506, 1.51677607024351 \}$.

SymPy uses the Berkowitz algorithm to compute eigenvalues, so I imagine it can find eigenvalues as accurately as the algorithm allows.

  • 1
    After perusing SymPy's docs, it seems to me that the method of Berkowitz involves computing the coefficients of the eigenpolynomial. While this might be fine for a small problem, it isn't an approach that can be recommended in general (for inexact arithmetic, that is; there's not much trouble computing coefficients when you can afford rational arithmetic). Besides, the symmetric eigenproblem has a very nice structure that can be leveraged, which the route of forming the characteristic polynomial ignores...2011-10-19