Let us proceed systematically:
- numerical precision of data (you said from medical imaging)
- number of operations required for standard methods (as from libraries)
- 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.
So, what are we left with: possible out-of-core methods? I can only think of an alternative: an out-of-core method that would be based on matri-vector products and would require reloading portions the matrix every time you need to compute a product. This require considerable expertise in maximising throughput, but please read on.
I have used to good effect (but other methods are available too) a combination of Lanczos method to generate a suitable subspace of vectors (see further on please), followed by some subspace power iteration. Subspace: a number of vector (> 20) which approximate the vector space of your largest eigenvectors (larhest: largest eigenvalues). I will not go through details here.
So first go through a series of Krylov method steps (Lanczos? CG?) to get some starting points for say m > 20 vectors. Then use subspace iteration to get the largest eigenvectors from the "pool" given by the subspace. It does work, provided there is a sufficient gap between eigenvalues (see above), and also allows for degeneracy (i.e. ultiple eigenvalues).
Cost: ~10^13 bytes transfer per matrix-vector product (matrix of course out-of-core). It requires ~m10^7 memory (if m is 100, say order of some Gigabytes, pretty trivial). if L (L>m) Lanczos (CG) iterations were required each would need at least one matrix-vector product, the cost would be > ~L10^12 flops, about 100 Tflops; if the subspace iteration required S iterations, the cost would be ~ Sm10^12, if S were again 100, than about 10 Petaflops in total.
So, provided you had a system that gave 1TFlops, you could carry out the computation in 100 (Lacnzos) + 10000 (iteration) seconds, or just about a couple of hours. However, it would require ~(L+S) matrix loads. So it would need, if we can load the matrices at 1Gbyte/sec (L+S)*10^13/10^9 ~ 10^6 seconds. But, here comes a useful point: it can be parallelised: if you split the matrix into "chunks" each processor could load a chunk of the matrix and some vectors could be duplicated (nbut they hold many fewer elements) although vectors should be added and synchronized across processors - but they are many fewer than matrix elements! So suppose, we split the lot across 100 processors, and ignoring synch times, the whole process could be carried out in (10^4+10^6)/100 ~ 10^4 seconds, just a few hours. So, it would be feasible this wy, BUT.... it requires quite a bit of expertise to write the necessary bespoke code.
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