Non explicit error message in lobpcg

[glemaitre]

In scikit-learn, we are benchmarking the integration of lobpcg to compute some svd. On an ill-posed problem, I got a weird error about array being None instead of a LinAlgError.

The issue is coming from the following line:

scipy/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py

Lines 105 to 108 in 91f4394

	except LinAlgError:
	#raise ValueError('Cholesky has failed')
	blockVectorV = None
	blockVectorBV = None

The LinAlgError is catched. Then the arrays are set to None for some reason but this state will not allow to resolve any issue later and just failed with a more cryptic error (see traceback below).
I was wondering if it was really intended to not fail early?

Reproducing code example:

Here, the matrix X should be transposed to be used in lobpcg and it is the reason for having an ill-posed problem which will raise an error.

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import lobpcg, aslinearoperator, LinearOperator

rng = np.random.RandomState(42)
X = sparse.random(32561, 123, density=0.4)

MLO = aslinearoperator(X)

def _matvec(V):
    return MLO(MLO.H(V))

A = LinearOperator(
    dtype=X.dtype, shape=(X.shape[0], X.shape[0]), matvec=_matvec,
    matmat=_matvec
)

n_components = 60
Q = rng.normal(size=(X.shape[0], 60))

lobpcg(A, Q, maxiter=10, verbosityLevel=0)

Error message:

Traceback (most recent call last):
  File "tmp.py", line 25, in <module>
    lobpcg(A, Q, maxiter=10, verbosityLevel=0)
  File "/home/lemaitre/miniconda3/envs/scipy-dev/lib/python3.7/site-packages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py", line 489, in lobpcg
    activeBlockVectorAR = A(activeBlockVectorR)
  File "/home/lemaitre/miniconda3/envs/scipy-dev/lib/python3.7/site-packages/scipy/sparse/linalg/interface.py", line 387, in __call__
    return self*x
  File "/home/lemaitre/miniconda3/envs/scipy-dev/lib/python3.7/site-packages/scipy/sparse/linalg/interface.py", line 390, in __mul__
    return self.dot(x)
  File "/home/lemaitre/miniconda3/envs/scipy-dev/lib/python3.7/site-packages/scipy/sparse/linalg/interface.py", line 420, in dot
    % x)
ValueError: expected 1-d or 2-d array or matrix, got array(None, dtype=object)

I was expected something like:

LinAlgError: 49-th leading minor of the array is not positive definite

Scipy/Numpy/Python version information:

>>> import sys, scipy, numpy; print(scipy.__version__, numpy.__version__, sys.version_info)
1.4.0.dev0+f3d6121 1.17.3 sys.version_info(major=3, minor=7, micro=4, releaselevel='final', serial=0)

[lobpcg]

@glemaitre Yes, it was really intended in scipy/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py Lines 105 to 108 in 91f4394 not not fail early.

It appears that the error is coming from LinearOperator, evaluating
A(activeBlockVectorR) possibly because activeBlockVectorR is None. I'll try to find time to double-check the code logic - why it gets to line 489 in lobpcg if activeBlockVectorR is indeed None. - Thanks for reporting!

Having said that, let me comment that your script is not using lobpcg well for svd:
size(X) is 32561x123 so to do SVD(X) the easy way one should use the normal matrix 123x123, which is MLO.H(MLO(V)) where V.shape[0]=123 , not your MLO(MLO.H(V)) which is size 32561x32561, and has all (but 123) zero eigenvalues.

[glemaitre]

Having said that, let me comment that your script is not using lobpcg well for svd:
size(X) is 32561x123 so to do SVD(X) the easy way one should use the normal matrix 123x123, which is MLO.H(MLO(V)) where V.shape[0]=123 , not your MLO(MLO.H(V)) which is size 32561x32561, and has all (but 123) zero eigenvalues.

Yes that's why I am expecting an error to be raised.

[lobpcg]

Having said that, let me comment that your script is not using lobpcg well for svd:
size(X) is 32561x123 so to do SVD(X) the easy way one should use the normal matrix 123x123, which is MLO.H(MLO(V)) where V.shape[0]=123 , not your MLO(MLO.H(V)) which is size 32561x32561, and has all (but 123) zero eigenvalues.

Yes that's why I am expecting an error to be raised.

Well, it does give an error :-)

However, lobpcg should still be expected to run on this example without errors, in my opinion, just given a symmetric matrix even with mostly zero eigenvalues.

[lobpcg]

OK, I get it and I have found a simple example with the same issue:

import numpy as np
from scipy import sparse
from scipy.sparse.linalg import lobpcg, aslinearoperator, LinearOperator

X = np.random.randn(100, 10)

def _matvec(V):
    return X @ (X.T @ V)

A = LinearOperator(
    dtype=X.dtype, shape=(X.shape[0], X.shape[0]), matvec=_matvec,
    matmat=_matvec)

Q = np.random.randn(X.shape[0], 4)
lobpcg(A, Q, maxiter=20, verbosityLevel=1)

always fails, while with Q = np.random.randn(X.shape[0], 3) it always works. The operator A of the eigenvalue problem is 100x100 and has only 10 nonzero eigenvalues. LOBPCG runs in the Krylov space generated by A that cannot thus have more than 10 dimensions. LOBPCG uses 3-times-block_size subspace to operate and runs out of space when block_size = 4, since 3-times-4 = 12 > 10, on Iteration 2. The iterative vectors become linearly dependent, so the orthogonalization fails, activating Lines 105 to 108 in scipy/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py :

 except LinAlgError: 
     #raise ValueError('Cholesky has failed') 
     blockVectorV = None 
     blockVectorBV = None 

Let me think how to handle it best... @glemaitre - thanks for finding and reporting such a case!

[FTB-B]

I am using scikit-learn spectral clustering for my clustering problem. I use the following configuration for the spectral clustering

clustering = sklearn.cluster.SpectralClustering(n_clusters = number_clusters , affinity="cosine",assign_labels="clusterQR",eigen_solver='lobpcg',n_jobs=psutil.cpu_count()).fit(embedding_matrix)

but I get the error

File "/data/fatemeh/mem2Vec/kym_meme/scikit-learn-clusterQR/sklearn/cluster/_spectral.py", line 559, in fit
    assign_labels=self.assign_labels)
  File "/data/fatemeh/mem2Vec/kym_meme/scikit-learn-clusterQR/sklearn/cluster/_spectral.py", line 301, in spectral_clustering
    eigen_tol=eigen_tol, drop_first=False)
  File "/data/fatemeh/mem2Vec/kym_meme/scikit-learn-clusterQR/sklearn/manifold/_spectral_embedding.py", line 339, in spectral_embedding
    largest=False, maxiter=2000)
  File "/home/ftahmas/venv/lib/python3.6/site-packages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py", line 489, in lobpcg
    activeBlockVectorAR = A(activeBlockVectorR)
  File "/home/ftahmas/venv/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 387, in __call__
    return self*x
  File "/home/ftahmas/venv/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 390, in __mul__
    return self.dot(x)
  File "/home/ftahmas/venv/lib/python3.6/site-packages/scipy/sparse/linalg/interface.py", line 420, in dot
    % x)
ValueError: expected 1-d or 2-d array or matrix, got array(None, `dtype=object)`

when I use affinity="rbf" it works without error!

any idea why?

[lobpcg]

@FTB-B I could investigate if you provide a reproducible example, please.

The issue should not be related to clusterqr, so please run with a different already available in scikit-learn function for labeling and submit a formal bug report with a ping to me.

[ogauthe]

I ran into the same problem as @glemaitre.

I have a (504, 504) np.float64 matrix A, I wish to compute a few singular values and singular vectors. scipy.linalg.svd and scipy.linalg.svds(solver='arpack') work well, however scipy.linalg.svds(solver='lobpcg') crashes with a very confusing error message. I can share the matrix if requested, github does not let me upload a npy file.

First, I think the line 105 raise ValueError('Cholesky has failed') should be uncommented to make the error message explicit (the error also comes from there).

Concerning the fix, my matrix has 78 non-zero singular values, but it crashes for k=5. Surprisingly, lobpcg crashes only for some values of k, worse when it does not crash the result is wrong. Comp

import numpy as np
import scipy.linalg as lg
import scipy.sparse.linalg as slg

A = np.load("data_crash.npy")
u0, s0, v0 = lg.svd(A)
x = lg.norm(A)
for k in range(1, 78):
    try:    
         u1, s1, v1 = slg.svds(A, k=k, solver="lobpcg")
         u2, s2, v2 = slg.svds(A, k=k, solver="arpack")
         print(k, "success")
         print("dense  truncation", lg.norm(u0[:,:k] * s0[:k] @ v0[:k] - A) / x) 
         print("arpack truncation", lg.norm(u2 * s2 @ v2 - A) / x) 
         print("lobpcg truncation", lg.norm(u1 * s1 @ v1 - A) / x) 
    except ValueError as err: 
         print(k, err)

gives me

1 success
dense  truncation 0.37410744042017613
arpack truncation 0.37410744042017613
lobpcg truncation 0.37410744042017613
2 success
dense  truncation 0.002749949291553053
arpack truncation 0.0027499492915530523
lobpcg truncation 1.2256989944723657
3 success
dense  truncation 0.0020247195761314294
arpack truncation 0.002024719576131428
lobpcg truncation 5.526570759513649
4 success
dense  truncation 0.0012373784439977066
arpack truncation 0.0012373784439977053
lobpcg truncation 0.06735836840347356
5 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
6 success
dense  truncation 0.0008082565817324968
arpack truncation 0.0008082565817324936
lobpcg truncation 0.0009858440032944039
7 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
8 success
dense  truncation 0.00044906990420738645
arpack truncation 0.0004490699042073811
lobpcg truncation 3.7862403491418792
9 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
10 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
11 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
12 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
13 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
14 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
15 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
16 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
17 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
18 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
19 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
20 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
21 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
22 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
23 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
24 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
25 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
26 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
27 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
28 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
29 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
30 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
31 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
32 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
33 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
34 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
35 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
36 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
37 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
38 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
39 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
40 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
41 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
42 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
43 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
44 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
45 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
46 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
47 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
48 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
49 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
50 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
51 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
52 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
53 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
54 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
55 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
56 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
57 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
58 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
59 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
60 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
61 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
62 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
63 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
64 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
65 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
66 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
67 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
68 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
69 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
70 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
71 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
72 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
73 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
74 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
75 expected 1-d or 2-d array or matrix, got array(None, dtype=object)
76 success
dense  truncation 8.826353786938982e-12
arpack truncation 2.7656568280904688e-06
lobpcg truncation 5.069589637576937
77 expected 1-d or 2-d array or matrix, got array(None, dtype=object)

[lobpcg]

@ogauthe this is a tricky matrix for lobpcg due to so many 0 eigenvalues; see #10974 (comment)

One simple fix could be just adding a small-magnitude random matrix to A to turn exact zeros into small nonzero values.

[ogauthe]

Thank you for your answer.

Beyond the error message when the function crashes, I think there is a real problem here. For k=2, the result is completely wrong, but this is silent. The matrix A.T @ A is symmetric positive by construction and quite small. Neither svds nor lobpcg documentation mention that the result may be wrong if too many eigenvalues are numerical zeros.

I can open a separate issue for this, a wrong result is much worse than a crash for me.

[lobpcg]

I have not paid attention to the values, sorry. Wrong values for small k is surely very odd, and unrelated to many zeros.

Please submit a separate bug report with a ping to me and upload the matrix to reproduce the problem.

[lobpcg]

When File "/home/ftahmas/venv/lib/python3.6/site-packages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py", line 489, in lobpcg activeBlockVectorAR = A(activeBlockVectorR) gives an error, the method cannot continue the iterations. But instead of crashing the code should exit gracefully outputting the current approximation.