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

I am trying to do eigenvalue decomposition for a huge matrix larger than 788000×788000 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 largest eigenvalues and their 20 corresponding eigenvectors.

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.

• As posed I very strongly suspect the problem is not solvable with current computer hardware. Assuming double precision to store the matrix requires ~46263933 GBytes of memory. The number 1 in the top 500 list of supercomputers has 5087232 GBytes (top500.org/system/179807), roughly 10% of what you need. Even if you can find a scalable out of core solvers this is still a phenomenal amount of disk space. You are going to have to approximate the problem somehow to get any sort of "solution". – Ian Bush Apr 23 at 11:12
• I haven't mentioned time as the memory makes the problem impossible without approximation. But just to give you an idea for us finding all eigenpairs of a dense, real symmetric order 124640 matrix takes about 4 minutes on 16384 cores of a modern supercomputer. Assuming exact O(N^3) scaling to solve a problem of the order you are proposing would take ~1923 years. Of course you only want a small proportion of the eigenpairs, but I think the number is indicative of the difficulty of what you want to do! – Ian Bush Apr 23 at 11:44
• I think the stated matrix size is not right. As you say, you have $512\times 512\times 300$ pixels, but that just means that you have 78800000 pieces of information. You don't have a $78800000 \times 78800000$ matrix. – Wolfgang Bangerth Apr 23 at 16:34
• More specifically, if you had a matrix this size, you would not be able to store it on the largest supercomputer in the world, so I suspect that you do not actually have a matrix of this size :-) – Wolfgang Bangerth Apr 23 at 16:35
• One quick comment: Are you sure you need to do whatever operation on your medical image in its entire domain? I mean you are saying that your image is a 3D array of 512x512x300 and you want to do some eigendecomposition operation on it, but I doubt that you actually have something useful on the entire of this 3D domain. In fact, I'm saying that if you are looking for an specific region of your medical image, e.g. liver or whatever organ, tissue, or vessel, etc., you might be able to crop your 3D image to only contain those interesting regions and then do whatever operation you want. – Alone Programmer Apr 23 at 19:35

Take a look at the literature that does similar things for facial recognition -- search for the term "eigenface", for example.

The point to make in this context is that the information you are looking for does not actually require you to consider high-resolution images. You may have $$10,000\times 10,000$$ pixels, for which any non-trivial computations are all but infeasible in the time scales you probably have in mind. But you really only care about the large-scale features of things when you look at eigenvector decompositions (the "eigenface" decomposition) and so people typically scale down their images to $$128\times 128$$ or similarly low resolutions, and do computations on these. For these sizes, eigenvector decompositions, principal component analyses, etc., are possible because the matrices you get are relatively small.

The underlying principle has been known in other communities for decades; for example, in the PDE solver communities this is called "multilevel" decompositions, and in other communities it goes by the name "hierarchical decomposition". The idea is that you might have a lot of data, but the information that is in it can be well approximated with far fewer bits. If, for example, you are interested in X-ray images of the whole body, you generally don't need extremely high resolution if your goal is to compare between individuals (for which you might want to compute "eigenskeletons"). It would be different if you were interested in the specific branching pattern of the lung of a single individual, but in that case you care about one picture, not the statistical properties of hundreds of imagines that you consider when you are interested in PCA or eigenvector decompositions.

As noted above by Thijs Steel a randomized svd is a solution but the number 78800000 is out of our computers computation ability.So you can proceed to the rsvd algorithm by :

import numpy as np
n  = 788
mu = 0
sigma = 1
A  = np.random.normal(mu, sigma, (n,n))
Omega  = np.random.normal(mu, sigma, (n,n))
def rsvd(A, Omega):
Y = A @ Omega
Q, _ = np.linalg.qr(Y)
B = Q.T @ A
u_tilde, s, v = np.linalg.svd(B, full_matrices = 0)
u = Q @ u_tilde
return u, s, v
u,s,v = rsvd(A,Omega)

• Thank you very much for all of you. I will try to implement the methods you have mentioned. – Michael Cheng Apr 23 at 13:27

Let us proceed systematically:

1. numerical precision of data (you said from medical imaging)
2. number of operations required for standard methods (as from libraries)
3. possible out-of-core computation (i.e. not the whole matrix at all times in memory). In all cases, I am afraid, you would have to be prepared to suffer. Incidentally, out-of-core methods are very unlikely to be availabe in any standard library like hose you mentioned. So you would need to do much coding yourelf or get suitable advice.

If your matrix is really ~10^6 x 10^6 elements, and it comes from medical imaging it is important to know how many significant digits (if you prefer your accuracy/noise level) you have. Too few significant digits, then your question depends very much on how strong the principal components (is that all about, is it not?) are wrt the other eigenvalues. If they are not well separated, than you could not separate the largest 20 from the others without loss of information. I also remind you here that the error in the eigenvectors depends on the sepratation of the eigenvalues (i.e. their diofference in absoloute value). If you had few digits available and the separation between your top 20 eigenvalues and te reast is small, than the error would be overwhelming and is possible that any solution would be misleading. Also, suppose that your data are 4 decimal digits accurate, then, given the size of your problem, you would have a very high "density" of eigenvalues which would make choosing just the top 20 at best a lottery. ou would need many more eigenvalues than 20.

Standard methods for such a matrix, would require ~10^13 bytes (~10 Tbytes) in double precision (which you need: single precision would be useless), which is a lot, of course. This cannot be handled by standard libraries. So this avenue is close. Also, very importantly, standard methids rely on some orthogonal reduction before computing the eigenvalues. These methods require O(N^3) operations and, beacause the matrix is symmetric (ouch!), these reductions are not very parallelisable (one global synchronisation for each row/column). So, you would need ~10^18 flops (all other stages would require the same order of magnitude of operations), which is about 1,000 Petaflops. The fastest computer system is about 90 PFlops in total, On that system, because of the poor parallelism of the orthogonalisation, you would need to wait several days/weeks/months for a soltion, if feasible. On a syystem capable of 1 TFlop, yo would need to wait for 1,000,000 hours, i.e. the time between Charlemagne (~800 ad) and now. So, the standard approach tghrough libraries cannot be done. It is an impossibility and following this path would only lead to frustration.

To summarize:

Point 1 is absolutely essential. If the numerics do not support the computation, do not do it. Any results would be useless.

Forget about standard libraries. Standard methods could not be used.

It could be done-ish by out-of-core iterative methods, but it would be a pretty steep path.

The overriding thing is point 1 above.

BTW, a constant image has, of course, just one non-zero eigenvalue. An image with random pixels, would have >> 20 eigenvalues.

hope this is of some help.

All the best

It can be done. Back in the 1990s certain groups (think radar cross section computation) were working with 1,000,000 x 1,000,000 dense systems. Consistent with other replies, this was done iteratively and out of core.