I have a symmetric matrix A. I'm not able to compute all the eigenvectors accurately, and I believe it is due to the last few eigenvalues for A being really, really minuscule.

Is there a method / library that will allow me to get really accurate eigenvalues/vectors, even if it's incredibly slow? (in R or Python)

(One suggestion was to iteratively eliminate the biggest eigenvalue/vector using "projectors", but I'm not sure how to go about this / if it would work)


Some context: The matrices I'm considering are basically Fisher Information matrices produced by taking the Hessian of the outputs of many-parameter models, each with significant internal parameter "redundancies", which you can uncover from an eigendecomposition. (It's not unusual to get really high condition numbers in these situations, as described here: https://www.lassp.cornell.edu/sethna/Sloppy/WhatAreSloppyModels.html.)

Here's some R code with an example of a matrix I'm interested in.

#A matrix we're interested in:

m <- matrix(c(
0.002307414, -1.318435973, 8.63e-11, -3.33e-11, 0.486238053, 0.41182779, -0.169525454, 0.624275558,
-1.318435973, 753.3425186, -4.93e-08, 1.90e-08, -277.8320732, -235.3147146, 96.86532753, -356.7054684,
8.63e-11, -4.93e-08, 3.23e-18, -1.24e-18, 1.82e-08, 1.54e-08, -6.34e-09, 2.34e-08,
-3.33e-11, 1.90e-08, -1.24e-18, 4.79e-19, -7.01e-09, -5.94e-09, 2.44e-09, -9.00e-09,
0.486238053, -277.8320732, 1.82e-08, -7.01e-09, 102.4642297, 86.78386443, -35.72384951, 131.5526701,
0.41182779, -235.3147146, 1.54e-08, -5.94e-09, 86.78386443, 73.50310589, -30.25693671, 111.4208258,
-0.169525454, 96.86532753, -6.34e-09, 2.44e-09, -35.72384951, -30.25693671, 12.45501408, -45.8654479,
0.624275558, -356.7054684, 2.34e-08, -9.00e-09, 131.5526701, 111.4208258, -45.8654479, 168.8989909
), nrow=8)

print('matrix m:')

#Get eigenvalues/vectors
eig <- eigen(m)

#Check - does Matrix*eigenvector == eigenvalue*eigenvector? (expected result is c(1,1,1,1,1,1,1,1)
print('Eigenvector 1')
print((m %*% eig$vectors[,1]) / (eig$values[1] * eig$vectors[,1]))   #Fine - yields c(1,1,1,1,1,1,1,1)
print('Eigenvector 2')
print((m %*% eig$vectors[,2]) / (eig$values[2] * eig$vectors[,2]))   #Mostly ok - all approximately 1
print('Eigenvector 3 (bad!)')
print((m %*% eig$vectors[,3]) / (eig$values[3] * eig$vectors[,3]))   #---Bad (elements vary a lot)!
print('Eigenvector 4 (really bad!)')
print((m %*% eig$vectors[,4]) / (eig$values[4] * eig$vectors[,4]))   #---Really bad!
print('Eigenvector 5')
print((m %*% eig$vectors[,5]) / (eig$values[5] * eig$vectors[,5]))   #Mostly ok
print('Eigenvector 6')
print((m %*% eig$vectors[,6]) / (eig$values[6] * eig$vectors[,6]))   #Mostly ok
print('Eigenvector 7')
print((m %*% eig$vectors[,7]) / (eig$values[7] * eig$vectors[,7]))   #Mostly ok
print('Eigenvector 8')
print((m %*% eig$vectors[,8]) / (eig$values[8] * eig$vectors[,8]))   #Mostly ok
  • 6
    $\begingroup$ What's the reason why you want to compute eigenvectors? (I ask because this may lead to different solutions) $\endgroup$ Commented Aug 7, 2021 at 12:24
  • 1
    $\begingroup$ Rough explaination: I take the hessian of a function, and find its eigenvectors. The eigenvectors with the smallest eigenvalues indicate how to change the variables in a way that least affects the function output. $\endgroup$
    – Mich55
    Commented Aug 7, 2021 at 12:52
  • 1
    $\begingroup$ In R, the package Rmpfr help to get really accurate number. But I don't know if it can be use to compute eigen value. Take a look $\endgroup$ Commented Aug 7, 2021 at 19:28
  • 3
    $\begingroup$ The matrix condition number for your example is on the order of 10^17!!! That poses extraordinary numerical problems. Please say more about how you are getting such an ill-conditioned Hessian in your application. I suspect that your real problem is at that earlier step, in formulating the function or using it to evaluate a data set. $\endgroup$
    – EdM
    Commented Aug 11, 2021 at 16:00
  • 2
    $\begingroup$ I don't know that "often have a huge condition number" is correct, at least in usual statistical practice. That's why, in large-data linear models, so much attention is paid to avoiding massive multicollinearity. Sounds like you have the equivalent problem here. Please edit your question to say more about where this type of matrix is coming from, as I don't think that your solution will be strictly computational. $\endgroup$
    – EdM
    Commented Aug 11, 2021 at 16:15

3 Answers 3


The problem is due to "leakage" from the large eigenvectors. I will present a brief analysis and then offer a solution, with code, followed by some remarks about the nature and limitations of the solution.


Let the eigenvalues of a real symmetric matrix $\mathbb M$ be $\lambda_1,\lambda_2,\ldots,\lambda_d$ ordered by decreasing absolute value so that $|\lambda_1|\ge |\lambda_2|\ge \cdots \ge |\lambda_d|.$ The objective, as stated in a comment to the question, is to find vectors $x$ so that $\mathbb{M}x$ is relatively small. Because the mapping $x\to\mathbb{M}x$ is linear, we may focus our study of such $x$ on those of unit length.

The Spectral Theorem guarantees the existence of an orthonormal basis $\{e_1,e_2,\ldots, e_d\}$ of corresponding eigenvectors in which $x$ may be expressed as the linear combination

$$x = \xi_1 e_1 + \xi_2 e_2 + \cdots + \xi_d e_d.$$

The orthonormality assures us that

$$1 = |x|^2 = \xi_1^2 + \xi_2^2 + \cdots + \xi_d^2.$$

Applying $\mathbb M$ yields

$$\mathbb{M}x = \xi_1 \mathbb{M} e_1 + \xi_2 \mathbb{M}e_2 + \cdots \xi_d \mathbb{M}e_d = \xi_1\lambda_1 e_1 + \xi_2\lambda_2 e_2 + \cdots + \xi_d\lambda_d e_d.$$

Suppose there is some slight error in the computation of the eigenvectors $e_i,$ so that where we think we are working with the eigenvector with eigenvalue $\lambda_i,$ we really are working with a perturbation $e_i^\prime$ where

$$e_i^\prime - e_i = \alpha_{i1}e_1 + \alpha_{i2}e_2 + \cdots + \alpha_{id}e_d.$$

These aren't quite eigenvectors anymore, as we may check by applying $\mathbb M$ to them:

$$\mathbb{M}e_i^\prime = \mathbb{M}e_i + \alpha_{i1}\mathbb{M}e_1 + \cdots + \alpha_{id}\mathbb{M}e_d = \lambda_i e_i + \alpha_{i1}\lambda_1 e_1 + \cdots + \alpha_{id}\lambda_d e_d.$$

Notice, in particular, that when $\lambda_i$ is tiny compared to $\lambda_1,$ the appearance of the multiple $\alpha_{i1}\lambda_1$ of $e_1$ can hugely alter the result, even when $\alpha_{i1}$ is tiny, due to that multiplication by $\lambda_1.$ This is the crux of the matter.

The problem propagates when applying $\mathbb{M}$ to the linear combination

$$x^\prime = \xi_1 e_1^\prime + \xi_2 e_2^\prime + \cdots + \xi_d e_d^\prime,$$

which we think is $x$ (because we think the $e_i^\prime$ are the $e_i$). Applying $\mathbb{M}$ separately to each of the $e_i^\prime$ on the right "leaks" potentially large multiples of $e_1$ into the result.


Fortunately there's a simple solution: remove the unexpected eigenvectors from the result. When (say) the first $k$ coefficients of $x$ are zero, $\xi_1=\xi_2=\cdots=\xi_k,$ then *there should not be any multiples of $e_1$ through $e_k$ in $\mathbb{M} x.$ We can remove them by projecting the result of $\mathbb{M}x$ onto the space generated by $e_{k+1},\ldots, e_d.$

This doesn't quite get us the correct value of $\mathbb{M}x,$ because of the "leakage" among the smaller eigenvectors. For instance, applying $\mathbb{M}$ to $x=e_d$ yields a linear combination of all the $e_i.$ Projecting out the first $k$ eigenvectors leaves us still with a linear combination

$$\lambda_{k+1}\alpha_{d,k+1} e_{k+1} + \lambda_{k+2}\alpha_{d,k+2} e_{k+2} + \cdots + \lambda_{d}(1+\alpha_{d,d}) e_{d}.$$

However, if (a) $|\lambda_{k+1}| \approx |\lambda_{k+2}| \approx \cdots \approx |\lambda_d|$ are of comparable orders of magnitude and (b) the $\alpha_{d,*}$ coefficients are all relatively small, the result will still be reasonably accurate.


Let's see how this works with the matrix $\mathbb M$ of the question. Here $k=8$ and $\lambda_1\approx 1111$ and $\lambda_2 \approx 5\times 10^{-8}$ are much greater than any of the other eigenvectors. The appearance of tiny negative eigenvalues in a setting where (one assumes) non-negative values are expected is further evidence that the last six (really seven) eigenvalues might just be noisy versions of zero.

The R code below creates random vectors $x^\prime$ in the space generated by the last $8-2=6$ eigenvectors. It applies $\mathbb M$ to them in two ways: directly as $y=\mathbb{M}x^\prime$; and with the "projection method" in which $y$ is regressed against the first two eigenvectors. It returns the norms of those two results. We hope that the norms will be really tiny, because $\mathbb M$ shouldn't change the norm by more than $\max\{|\lambda_3|,\ldots, |\lambda_8|\}\approx 6.25\times 10^{-8}.$ To see what it really does, I plotted histograms of the two sets of results (left and middle). To compare the two approaches I also drew their scatterplot at the right.


The direct method (left) is awful: a considerable amount of $e_1$ has leaked into the eigenvectors, causing the norms to be inflated by approximately ten(!) orders of magnitude. After projecting the first two eigenvectors out, the norms are always of the expected magnitude. The scatterplot shows no relationship among the results. (In some cases, for other matrices $\mathbb M,$ there are some intriguing relationships induced by the magnitudes of the error coefficients $\alpha_{ij}.$)


This solution guarantees that the result of computing $\mathbb{M}x$ will be orthogonal to the largest eigenspaces and will, with an accuracy relative to the magnitudes of the largest eigenvalues, reflect how much $\mathbb M$ affects the lengths. That's probably all that's needed in the intended application.

There's no prospect of doing any better, either. Many of the entries in matrices like $\mathbb M$ were already computed using floating point arithmetic and therefore cannot be considered any more precise than the results of those operations. The large condition numbers of matrices like $\mathbb M$ guarantee that even using infinitely precise arithmetic, the effects of either (a) changing any entry in $\mathbb M$ in its least significant bit or (b) changing any coefficient of $x$ in its least significant bit will hugely perturb the result.

The svd (singular value decomposition) function in R does a better job of computing the eigenvalues and vectors. This doesn't improve the direct calculation of $\mathbb M$ appreciably, though. If you would like to compare these methods, replace the line eig <- eigen(m) by eig <- with(svd(m), list(values = d, vectors = u)).


This picks up after the definition of m in the question and the calculation of its eigenvectors using eigen.

lambda <- eig$values
i <- lambda >= max(lambda) * 1e-12 # Indexes of large eigenvalues.
# Generate unit movements in the small directions.
nu <- function(x) { # Standardize to unit length if possible
  a <- sum(x^2)
  if(a==0) x else x/sqrt(a)
sim <- replicate(500, {
  # Naive (uncorrected) algorithm.
  dx <- eig$vectors %*% rnorm(length(lambda)) * ifelse(i, 0, 1)
  dy <- m %*% nu(dx)
  # Corrected algorithm with projection.
  dx <- residuals(lm(dx ~ eig$vectors[,i]))
  dz <- m %*% nu(dx)
  c(sqrt(sum(dy^2)), sqrt(sum(dz^2)))
# Compare.
hist(sim[1,], main="Magnitudes of resultants: direct method", cex.main=1)
hist(sim[2,], main="Magnitudes of resultants: projection  method", cex.main=1)
plot(t(sim), xlab="Direct method", ylab="Projection method", main="Comparison", cex.main=1)


Whuber has already done a better job of that than I could of answering the question, but I have some figures I generated for other work may help illustrate some of the issues. A similar version is in the supporting information associated with this paper on visual exploration of PCA for spectroscopists (open access)

Whuber indicated that SVD is good for accurate calculation of eigenvalues and vectors, and from my reading that is a widely held opinion, this arises from the use of matrix decomposition on the covariance matrix, which has the consequence that all eigenvectors are solved in conjunction and computation errors are spread over them all. Sequential/iterative methods such as the power method (sequential fitting of covariance matrix) or NIPALs (sequential fitting of residual) in contrast accrue errors in each eigenvector calculation and these can build up to the point they over power very low variance eigenvectors.

In that paper I had implemented PCA using NIPALs and as a sanity check I published data in the supporting information comparing the results to an established python function that used SVD. The first PC was very similar for both algorithms:

Plot overlaying the first PC calculated using the NIPALS algorithm and that calculated by the SKLearn SVD based PCA

The difference between them was miniscule ($10^{-15}$), on the order of computational accuracy: difference between PC1 calculated by NIPALs and SVD

The difference was dependent on the specified tolerance for the NIPALs algorithm. Next we see a plot of the SVD rank against the NIPALs rank for the SVD-eigenvector that the NIPALs eigenvector most closely correlated with:

SVD rank vs NIPALs rank for correlation matched eigenvectors

At first the rank is perfectly correlated, until we hit NIPALS PC24 where things get wobbly. NIPALS-24 is suddenly best correlated with SVD-1! Above NIPALS-24 the correlations get really messed up. It is not limited to rank, but actual eigenvector shape becomes much less consistent (even if we match by optimum correlation), as seen if we plot the correlation between SVD/NIPALs vs NIPALs rank:

correlation between SVD/NIPALs eigenvector correlation vs NIPALs rank

This shows that information in the lower rank PCs is getting muddled up in the NIPALS wrt to SVD. If we look at SVD-1 (blue) and NIPALS-24 (orange) and the difference between them (green) we can seen that the two are indeed very similar and the difference is very noisy.

Comparing SVD-1 and NIPALS-24 and their difference

As you mention in the OP and Whuber mentions, this is usually handled by projection which also goes by the name of renormalisation or re-orthogonalisation . The process starts by applying the $PC_{1-i-1}$ eigenvectors to the residual and reconstructing (inner product of scores and eigenvectors), then recalculating the residual. If we done this after NIPALS-23 then the residual would have popped up as mostly NIPALS-1. Reconstructing the data at this point would have created a residual without the NIPALS-1 ghost and therefore would not have caused the same level of problems.

As Whuber indicates it is not possible to achieve infinitesimal accuracy, so it is a matter of what limits of precision do you need.


If your eigenvalues are "sloppy" then the eigenvectors will necessarily be "sloppy" too. With a condition number on the order of $10^{17}$ as in your example, I'd be wary of any results you got regardless of what numerical approach you use.

The solution is to avoid this type of problem before it starts.

You don't provide details about the "multi-parameter models" that you're trying to fit, but something seems to be terribly wrong. For an 8 x 8 Hessian matrix, there is only 1 eigenvalue of any substantial magnitude.

> eig$values
[1]  1.110666e+03  5.014476e-08  7.523231e-14  1.113187e-14 -4.122444e-10
[6] -2.503845e-09 -1.832745e-08 -6.251041e-08

I suspect that there is substantial multicollinearity in your data, although other problems like an over-specified model might also be in play.

My first reaction is to use principal components of your data as predictors instead of the raw data. If multicollinearity is to blame for your problem, using just a few principal components could turn that weakness into a strength. Think also about whether the structure your model is contributing.

  • $\begingroup$ Thanks for the collinearity and PCA tips. I'm still convinced that some of my models are genuinely producing these matrices so can I ask, is the short answer: no, for matrices with condition numbers this high, it isn't really possible to compute accurate eigenvalues/vectors? As I mentioned to @whuber in a comment, it would be very useful to even estimate these, if possible. Thanks for your answer! $\endgroup$
    – Mich55
    Commented Aug 11, 2021 at 20:45

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.