scipy: LOBPCG convergence failure if guess provided

When using the lobpcg implementation in scipy.sparse.linalg with scipy version 1.3.0 on a preconditioned, normal eigenproblem, the algorithm aborts with the error message shown below if a particular guess is employed. This is surprising, because the guess is derived from an approximate physical model and a single Rayleigh-Ritz step in the subspace spanned by the guess vectors already provides a very good approximation to the desired eigenvalues.

Reproducing code example:

Can be found in mfherbst/issue_scipy_lobpcg_guess. For running the code the datafile is required. On execution the following output is produced on my machine:

scipy version:  1.3.0
Problem properties:
N:                  411
m:                  10
λref:               [0.21587252 0.7147382  0.72196884 0.72214299 0.72648404 0.73109012
 0.73639964 0.80308925 0.8303687  0.91917443]
λapprox:            [0.286794   0.75381381 0.75949557 0.76511668 0.76580448 0.77711184
 0.78372357 0.86710816 0.95792697 0.96882807]
A λ dist:           0.011713478159563806
P*A ritz dist:      0.046591195439390276

read_guess=True
ERROR:  9-th leading minor of the array is not positive definite

read_guess=False

As can be deduced from the depicted output, the first execution, which employs the guess, fails, whereas the second (employing a random guess) succeeds. λref are the eigenvalues on the full matrix from LAPACK, λapprox the Ritz values inside the subspace spanned by the guesses, λ dist and ritz dist are the absolute differences between the last eigenvalue / Ritz value of interest and the next value.

Error message:

Traceback (most recent call last):
  File "./test_lobpcg.py", line 47, in <module>
    test_lobpcg(read_guess=True)
  File "./test_lobpcg.py", line 39, in test_lobpcg
    λ, v = lobpcg(A, X, largest=False, M=P, maxiter=100, tol=tol)
  File "/home/mfh/.julia/python_virtualenv/lib/python3.7/site-packages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py", line 476, in lobpcg
    aux = _b_orthonormalize(B, activeBlockVectorP, retInvR=True)
  File "/home/mfh/.julia/python_virtualenv/lib/python3.7/site-packages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py", line 93, in _b_orthonormalize
    gramVBV = cholesky(gramVBV)
  File "/home/mfh/.julia/python_virtualenv/lib/python3.7/site-packages/scipy/linalg/decomp_cholesky.py", line 91, in cholesky
    check_finite=check_finite)
  File "/home/mfh/.julia/python_virtualenv/lib/python3.7/site-packages/scipy/linalg/decomp_cholesky.py", line 40, in _cholesky
    "definite" % info)
numpy.linalg.LinAlgError: 9-th leading minor of the array is not positive definite

Scipy/Numpy/Python version information:

1.3.0 1.16.2 sys.version_info(major=3, minor=7, micro=3, releaselevel='final', serial=0)

About this issue

Original URL
State: closed
Created 5 years ago
Comments: 19 (10 by maintainers)

Commits related to this issue

multiple stability updates multiple stability updates to get the code closer in sync with the latest version of https://github.com/lobpcg/blopex/blob/master/blopex_tools/matlab/lobpcg/lobpcg.m ... — committed to aburrell/scipy by lobpcg 5 years ago

Most upvoted comments

With the changes suggested in #10621 the issue in the provided example will be resolved.

mfherbst on Aug 19, 2019

@rc @mfherbst @joseeroman @antoine-levitt @mohamed82008

I think that I have found easy fixes:

https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py follows https://arxiv.org/abs/0705.2626 implemented in https://github.com/lobpcg/blopex/blob/master/blopex_abstract/krylov/lobpcg.c , which is a somewhat simplified and inexact version of my MATLAB implementation https://github.com/lobpcg/blopex/blob/master/blopex_tools/matlab/lobpcg/lobpcg.m

Comparing my latest MATLAB implementation https://github.com/lobpcg/blopex/blob/master/blopex_tools/matlab/lobpcg/lobpcg.m with this scipy code, I see 3 steps missing.

The most important one is MATLAB lines 601-608

    %Making active (preconditioned) residuals orthogonal to blockVectorX
    if isempty(operatorB)
        blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
            blockVectorX*(blockVectorX'*blockVectorR(:,activeMask));
    else
        blockVectorR(:,activeMask) = blockVectorR(:,activeMask) - ...
            blockVectorX*(blockVectorBX'*blockVectorR(:,activeMask));
    end

that should have been before line 460 of https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py but not there.

If the directions P are linearly depended at a given iteration, which is the present issue at the first iteration where P appears, the easiest solution is to drop the whole P altogether reducing to the steepest descent, but only at this iteration. In other words, each time check if Cholesky fails for orthonormalization of P and if so, drop P. I have tested in https://github.com/lobpcg/blopex/blob/master/blopex_tools/matlab/lobpcg/lobpcg.m that the present issue gets fixed without reducing the convergence speed, since the trouble with P actually only appears once. I’ll try to fix the python version accordingly and report here.
The last big omission in the BLOPEX-C and scipy versions is the use explicitGramFlag = 1; %suggested by Garrett Moran in MATLAB that forces explicitly computing several blocks in the Gram matrices for the Rayleigh-Ritz method, see starting at line 709 of https://github.com/lobpcg/blopex/blob/master/blopex_tools/matlab/lobpcg/lobpcg.m But this only kicks off at small tolerances so not as crucial as the other 2 omission.

lobpcg on Jun 12, 2019