Eigenvalue decomposition for a very huge matrix of medical images (such as the pixel physical coordinates of CT images)

Question

_{Crossposted on Computational Science SE}

---

I am trying to do eigenvalue decomposition for a huge matrix larger than $78.8 \cdot 10^4 \times 78.8 \cdot 10^4$ for medical image analysis. The matrix is not sparse and every element in the matrix has a real value. And, for example, I want to obtain the first $20$ eigenvectors corresponding to the first $20$ largest eigenvalues.

The computer is not able to do eigenvalue decomposition for the huge matrix and the memory overflows, although my computer configuration is very excellent. I write the computer codes with Python language and other related packages (such as NumPy, OpenCV, Matplotlib and so on). Is there any other Python library or related package that can do eigenvalue decomposition and solve the computation problem? Or, is there any other method that can solve this problem with Python?

I am in a difficult situation now, and hope someone can help me. Thank you so much.

So sorry, I wrote wrongly, the huge matrix is ​​also symmetric.

Answer 1

A Krylov based method would be to recommend as you can get away with matrix-vector operations.

Any method that tries to store and manipulate actual matrices will require too much memory.

You can read more about Krylov subspace methods for example at Wikipedia

Answer 2

The eigenvalues of an upper or lower triangular matrix are the diagonal entries of the matrix

So the problem now is how to convert your original matrix into a triangular one.

Because you matrix is huge, I suggest a home-made way.
The most basic method is to find k factors between row j-1 and row j that makes matrix[j,k]=0 for k<j. There are plenty of literature about this, so I skip explanations here.

The thing is that with this simple method you don't have to load the whole matrix into memory, but just the two rows you are working with.
But you must store the modified rows because they will be used in further rows modification.

This method is slow, and likely prone to wrong values due to roundings in calculations. To make results a bit less wrong you may calculate both upper and lower triangulations and get the average diagonal values.