MATLAB's EIGS vs PRIMME_EIGS

[Mbakhov]

Dear developers,

I'm planning to use PRIMME to calculate 10 interior eigenvalues and eigenvectors of large sparse matrices (~ 810^5 x 810^5) in Matlab. However, I found that PRIMME is ~ 15 times slower in moderately small matrices ( 3K x 3K) with respect to standard eigs function; Eigs is performing the task in ~ 0.05 sec, while for GD_Olsen_plusK it takes 14 seconds (in 8 cores, 3Ghz each).

So, is it possible to speed up the calculation process?

load H.mat;

[~,Emin]=eigs(H,1,'sa');   
[~,Emax]=eigs(H,1,'la');

 e=(Emax-Emin)/2+Emin;
 
tic
  opts = struct();
D = primme_eigs(H,10,e,opts,'JD_Olsen_plusK' );  
toc
%%
tic
[XX EE]=eigs(H,10,e);
toc

Thank you,

Best regards,

Murod

[Mbakhov]

H.zip

[stathopoulos]

Dear Murod, The reason eigs() is faster is that, for interior eigenvalues, it factorizes (A-e I) and runs Lanczos on (A-eI)^{-1}. Then it converges in a few iterations. However, the factorization is O(N^3) process and while very fast for 3K becomes quickly a bottleneck for larger sizes. PRIMME uses only matrix vector multiplications (and a preconditioner if available). Therefore convergence will be much slower (especially for interior eigenvalues) but every iteration is much faster (O(Matvec)). That’s why these methods are best suited for large sparse matrices. You may want to investigate where the cross over point is on your matrix size and your machine where factorizing (A-e I) becomes too expensive. If you let us know of the specifics of your problem, we may be able to provide more guidance. Andreas

…

On Jul 28, 2020, at 10:38 AM, Mbakhov ***@***.***> wrote: Dear developers, I'm planning to use PRIMME to calculate 10 interior eigenvalues and eigenvectors of large sparse matrices (~ 810^5 x 810^5) in Matlab. However, I found that PRIMME is ~ 15 times slower in moderately small matrices ( 3K x 3K) with respect to standard eigs function; Eigs is performing the task in ~ 0.05 sec, while for GD_Olsen_plusK it takes 14 seconds (in 8 cores, 3Ghz each). So, is it possible to speed up the calculation process? load H.mat; [~,Emin]=eigs(H,1,'sa'); [~,Emax]=eigs(H,1,'la'); e=(Emax-Emin)/2+Emin; tic opts = struct(); D = primme_eigs(H,10,e,opts,'JD_Olsen_plusK' ); toc %% tic [XX EE]=eigs(H,10,e); toc Thank you, Best regards, Murod — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#35>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADDKVJI6TFICOTFGJV3T5LDR53PGNANCNFSM4PKS5GJQ>.

----------------------------------------------------------------------------- Andreas Stathopoulos Professor William & Mary E-mail: andreas@cs.wm.edu Department of Computer Science Phone: 757-221-3483 P.O. Box 8795 Fax: 757-221-1717 Williamsburg, VA 23187-8795 http://www.cs.wm.edu/~andreas

[eromero-vlc]

Thanks for providing a matrix, that always helps! As Andreas said, eigs relays on doing the LU factorization of (H-e*eye(size(H))). For cases where finding interior eigenvalues takes many iterations, doing LU and using eigs may be the fastest choice if you can afford to do the factorization.

I played a little bit with the matrix that you provided. It seems that using refined extraction helps. Also, for many applications, it is sufficient to set the accuracy of the solver to 3-4 order of magnitudes below the average distance of the target eigenvalues. For H, the minimum distance between nearby eigenvalues is 1e-4 around e, and the norm of H is around 8.4, so a tol=avg_eig_dist/norm(H)*1e-3=1e-8 should be ok.

l=primme_eigs(H,10,e,struct('projection','primme_proj_refined','tol',1e-8),'DEFAULT_MIN_TIME')

[Mbakhov]

Dear Andreas,

Thanks for your quick respond.

I understood the point.

The matrices which I consider are real, hermitian and sparse. The dimension is very large ~ 850 000 x 850 000. Please, find the link to the example matrix.

I need :
10-20 eigenvalues and eigenvectors from the mid-spectrum.

So, my question is what would be the expected cpu memory requirements and number of cores to solve the problem in a reasonable timescale ? and whether it is practical to use PRIMME for such large matrices?

Thank you,

Best regards,

Murod

[Mbakhov]

Thanks for providing a matrix, that always helps! As Andreas said, eigs relays on doing the LU factorization of (H-e*eye(size(H))). For cases where finding interior eigenvalues takes many iterations, doing LU and using eigs may be the fastest choice if you can afford to do the factorization.

I played a little bit with the matrix that you provided. It seems that using refined extraction helps. Also, for many applications, it is sufficient to set the accuracy of the solver to 3-4 order of magnitudes below the average distance of the target eigenvalues. For H, the minimum distance between nearby eigenvalues is 1e-4 around e, and the norm of H is around 8.4, so a tol=avg_eig_dist/norm(H)*1e-3=1e-8 should be ok.

l=primme_eigs(H,10,e,struct('projection','primme_proj_refined','tol',1e-8),'DEFAULT_MIN_TIME')

Thank you!
Indeed, it helps. However, as I addressed in the previous post, it will take enormous time to perform the same for the matrix of ~ 850 000 x 850 000.

Best regards,

Murod

[stathopoulos]

Murod, If we assume your spectrum has similar separations at the large size as at the small size, the number of iterations will be similar for your 850K matrix, so the time complexity goes up linearly with the size. This is the best case scenario. Typically you can expect denser separations as size increases so you need to adjust the tolerances accordingly (and the time will go up). Let me say that we often run matrices of size more than a million, so I suggest you give it a try. Andreas

…

On Jul 28, 2020, at 12:36 PM, Mbakhov ***@***.***> wrote: Thanks for providing a matrix, that always helps! As Andreas said, eigs relays on doing the LU factorization of (H-e*eye(size(H))). For cases where finding interior eigenvalues takes many iterations, doing LU and using eigs may be the fastest choice if you can afford to do the factorization. I played a little bit with the matrix that you provided. It seems that using refined extraction helps. Also, for many applications, it is sufficient to set the accuracy of the solver to 3-4 order of magnitudes below the average distance of the target eigenvalues. For H, the minimum distance between nearby eigenvalues is 1e-4 around e, and the norm of H is around 8.4, so a tol=avg_eig_dist/norm(H)*1e-3=1e-8 should be ok. l=primme_eigs(H,10,e,struct('projection','primme_proj_refined','tol',1e-8),'DEFAULT_MIN_TIME') Thank you! Indeed, it helps. However, as I addressed in the previous post, it will take enormous time to perform the same for the matrix of ~ 850 000 x 850 000. Best regards, Murod — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#35 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADDKVJK2FYW2OYBB5MUTXU3R53477ANCNFSM4PKS5GJQ>.

----------------------------------------------------------------------------- Andreas Stathopoulos Professor William & Mary E-mail: andreas@cs.wm.edu Department of Computer Science Phone: 757-221-3483 P.O. Box 8795 Fax: 757-221-1717 Williamsburg, VA 23187-8795 http://www.cs.wm.edu/~andreas

[tanlin2013]

Hi,

A related question concerning about extracting interior eigenvalues and eigenvectors from a large (Hermition) matrix. I found a rather old paper about this kind of problem. But as an end user, I'm not really certain the efficiency between this modified Rayleigh-Ritz method and Jacobi-Davison-like algorithms. The number of interrior eigenpairs I need is in general around few hundreds, while the size of matrix is about (10^7 * 10^7) also.

I'm still on early planing to see whether this problem is feasible with numerical methods. Hopefully there're some answers in this kind of problem already. Thanks!

best,
Tan

[stathopoulos]

Hello, The paper you refer to introduces what is now known as Harmonic Ritz eigenvalues. This not a separate iterative method but a different way to extract information from the projection space we build (with Krylov methods, Jacobi-Davidson methods, or any iterative method you have). It provides monotonic convergence of the values towards a given shift. There are some cases where it may avoid missing eigenvalues, especially if you use a method that restarts very frequently (i.e., a small basis size). However, in my experience, it is usually slower than Rayleigh-Ritz. In your case, needing 100s of eigenvalues in the middle of the spectrum implies that you will need a large basis size and you will have long run times which tends to recover some eigenvalues if missed. Either way, PRIMME implements Harmonic Ritz value extraction, along with Rayleigh-Ritz, and Refined-Ritz. You can experiment with these three and see which one performs best for your problem. Best, Andreas

…

On Sep 14, 2021, at 5:21 AM, Tan Tao-Lin ***@***.*** ***@***.***>> wrote: Hi, A related question concerning about extracting interior eigenvalues and eigenvectors from a large (Hermition) matrix. I found a rather old paper <https://www.sciencedirect.com/science/article/pii/0024379591903816> about this kind of problem. But as an end user, I'm not really certain the efficiency between this modified Rayleigh-Ritz method and Jacobi-Davison-like algorithms. The number of interrior eigenpairs I need is in general around few hundreds, while the size of matrix is about (10^7 * 10^7) also. I'm still on early planing to see whether this problem is feasible with numerical methods. Hopefully there're some answers in this kind of problem. Thanks! best, Tan — You are receiving this because you commented. Reply to this email directly, view it on GitHub <#35 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ADDKVJP4OH4EAZ35D3LCQB3UB4H2LANCNFSM4PKS5GJQ>. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>.

[tanlin2013]

Thank you Andreas! I will give it a try.

Then how exactly I could switch to Rayleigh-Ritz in PRIMME? I do find it in the documentation. But if there some some examples that will be helpful. Thanks!

p.s. I'm using python actually.

Best,
Tan

[eromero-vlc]

PRIMME by default uses Rayleigh-Ritz. You can change the extraction method with the option projection_projection. The documentation of the option is here as you pointed out. Refined and harmonic only works when finding the eigenvalues closest to a target, which is passed in the option which. These are some examples:

>>> # Rayleigh-Ritz extraction (default)
>>> evals, evecs,stats = primme.eigsh(A, 3, tol=1e-6, which=50.1, projection_projection='primme_proj_RR', return_stats=True); stats['numMatvecs']
212
>>> # Refine extraction
>>> evals, evecs,stats = primme.eigsh(A, 3, tol=1e-6, which=50.1, projection_projection='primme_proj_refined', return_stats=True); stats['numMatvecs']
235
>>> # Harmonic-Ritz extraction
>>> evals, evecs,stats = primme.eigsh(A, 3, tol=1e-6, which=50.1, projection_projection='primme_proj_harmonic', return_stats=True); stats['numMatvecs']
191