# [SciPy-User] scipy.sparse.linalg.lobpcg broken?

Nico Schlömer nico.schloemer@gmail....
Wed Oct 13 14:23:19 CDT 2010

```> It actually works for
> me:

Okay, nice. I tried to run the script on another machine too, where it
appears to work alright.
This exact same code fails on my own machine locally, meaning I'm
probably missing some crucial dependencies (of Scipy?). Using Scipy
0.8.0 on Gentoo.
Any idea what that could possibly be?

Cheers,
Nico

On Wed, Oct 13, 2010 at 7:09 PM, Robert Cimrman <cimrman3@ntc.zcu.cz> wrote:
> On Wed, 13 Oct 2010, Nico Schlömer wrote:
>>
>> Possibly the devs know -- how to get back to them?
>
> If I recall correctly, any non-zero X should be ok. It actually works for
> me:
>
> In [20]: X = np.random.rand( n,1 )
>
> In [21]: lobpcg( A, X )
> Out[21]:
> (array([ 1.]),
>  array([[ 0.24881513],
>       [ 0.35018185],
>       [ 0.47060887],
>       [ 0.0544507 ],
>       [ 0.2243659 ],
>       [ 0.22597527],
>       [ 0.55431323],
>       [ 0.40576857],
>       [ 0.05081581],
>       [ 0.12299473]]))
>
> In [23]: lobpcg( A, X )
> Out[23]:
> (array([ 1.]),
>  array([[ 0.31622777],
>       [ 0.31622777],
>       [ 0.31622777],
>       [ 0.31622777],
>       [ 0.31622777],
>       [ 0.31622777],
>       [ 0.31622777],
>       [ 0.31622777],
>       [ 0.31622777],
>       [ 0.31622777]]))
>
> The eigenvalue is correctly computed (lambda = 1). Since A is a unit matrix
> and we are solving A x = lambda x, any nonzero x is an eigenvector.
>
> r.
>
>> --Nico
>>
>>
>> On Wed, Oct 13, 2010 at 5:32 PM,  <josef.pktd@gmail.com> wrote:
>>>
>>> On Wed, Oct 13, 2010 at 10:58 AM, Nico Schlömer
>>> <nico.schloemer@gmail.com> wrote:
>>>>
>>>> In fact, for general
>>>>  X = np.random.rand( n,1 )
>>>> this thing always blows up -- look at the residual norms.
>>>>
>>>> Bug report?
>>>> Nico
>>>>
>>>>
>>>>
>>>> On Wed, Oct 13, 2010 at 4:50 PM, Nico Schlömer
>>>> <nico.schloemer@gmail.com> wrote:
>>>>>
>>>>> Really? I tried ones() instead, and got (with verbosity=10)
>>>>>
>>>>> ======================= *snip* =======================
>>>>> Solving generalized eigenvalue problem with preconditioning
>>>>>
>>>>> matrix size 10
>>>>> block size 1
>>>>>
>>>>> No constraints
>>>>>
>>>>>
>>>>> iteration 0
>>>>> [ True]
>>>>> current block size: 1
>>>>> eigenvalue: [ 100.]
>>>>> residual norms: [ 990.]
>>>>> iteration 1
>>>>> [ True]
>>>>> current block size: 1
>>>>> eigenvalue: [ 0.]
>>>>> residual norms: [  9.60596010e+12]
>>>>> iteration 2
>>>>> [ True]
>>>>> current block size: 1
>>>>> eigenvalue: [ 0.]
>>>>> residual norms: [  1.63581388e+65]
>>>>> iteration 3
>>>>> Warning: invalid value encountered in multiply
>>>>> [False]
>>>>> Warning: invalid value encountered in multiply
>>>>> final eigenvalue: [ 0.]
>>>>> final residual norms: [ nan]
>>>>> ======================= *snap* =======================
>>>>>
>>>>> We're still talking about the identity matrix, so I don't expect this
>>>>> breakdown to be inherent in the method.
>>>>>
>>>>> Cheers,
>>>>> Nico
>>>
>>> noisy identity matrix seems to work
>>> lobpcg( np.diag(1+1e-6*np.random.randn(10)), np.random.randn(10,1),
>>> verbosityLevel=10)
>>> lobpcg( sparse.csr_matrix(np.diag(1+1e-4*np.random.randn(10))),
>>> np.random.randn(10,1), verbosityLevel=10)
>>>
>>> I'm not sure what this means (since I'm no expert on this)
>>>
>>>>>> X = np.zeros( (n,1) )
>>>>>> X[-3:]=.1
>>>>>> lobpcg( A, X, verbosityLevel=10)
>>>
>>> Solving generalized eigenvalue problem with preconditioning
>>>
>>> matrix size 10
>>> block size 1
>>>
>>> No constraints
>>>
>>>
>>> iteration 0
>>> [False]
>>> final eigenvalue: [ 1.]
>>> final residual norms: [  1.92296269e-16]
>>> (array([ 1.]), array([[ 0.        ],
>>>       [ 0.        ],
>>>       [ 0.        ],
>>>       [ 0.        ],
>>>       [ 0.        ],
>>>       [ 0.        ],
>>>       [ 0.        ],
>>>       [ 0.57735027],
>>>       [ 0.57735027],
>>>       [ 0.57735027]]))
>>>
>>>
>>> final residual norm = 0.
>>>
>>>>>> X = np.ones( (n,1) )
>>>>>> lobpcg( A, X, verbosityLevel=10)
>>>
>>> Solving generalized eigenvalue problem with preconditioning
>>>
>>> matrix size 10
>>> block size 1
>>>
>>> No constraints
>>>
>>>
>>> iteration 0
>>> [False]
>>> final eigenvalue: [ 1.]
>>> final residual norms: [ 0.]
>>> (array([ 1.]), array([[ 0.31622777],
>>>       [ 0.31622777],
>>>       [ 0.31622777],
>>>       [ 0.31622777],
>>>       [ 0.31622777],
>>>       [ 0.31622777],
>>>       [ 0.31622777],
>>>       [ 0.31622777],
>>>       [ 0.31622777],
>>>       [ 0.31622777]]))
>>>
>>> I have no idea if there are some inherent problems with starting
>>> values and whether lobpcg is supposed to converge from any starting
>>> values.
>>>
>>> Josef
>>>
>>>
>>>>>
>>>>>
>>>>>
>>>>> On Wed, Oct 13, 2010 at 4:29 PM,  <josef.pktd@gmail.com> wrote:
>>>>>>
>>>>>> On Wed, Oct 13, 2010 at 10:21 AM, Nico Schlömer
>>>>>> <nico.schloemer@gmail.com> wrote:
>>>>>>>
>>>>>>> Hi,
>>>>>>>
>>>>>>> I thought I give lobpcg a shot, and tried
>>>>>>>
>>>>>>> ====================== *snip* ======================
>>>>>>> from scipy.sparse.linalg import lobpcg
>>>>>>> from scipy.sparse import identity
>>>>>>> import numpy as np
>>>>>>>
>>>>>>> n = 10
>>>>>>> X = np.zeros( (n,1) )
>>>>>>> A = identity( n )
>>>>>>> lobpcg( A, X )
>>>>>>> ====================== *snap* ======================
>>>>>>>
>>>>>>> On my machine, this yields
>>>>>>>
>>>>>>> ====================== *snip* ======================
>>>>>>> Traceback (most recent call last):
>>>>>>>  File "logpcg_test.py", line 8, in <module>
>>>>>>>    lobpcg( A, X )
>>>>>>>  File
>>>>>>> "/usr/lib64/python2.6/site-packages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py",
>>>>>>> line 304, in lobpcg
>>>>>>>    blockVectorX, blockVectorBX = b_orthonormalize( B, blockVectorX )
>>>>>>>  File
>>>>>>> "/usr/lib64/python2.6/site-packages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py",
>>>>>>> line 130, in b_orthonormalize
>>>>>>>    gramVBV = sla.cholesky( gramVBV )
>>>>>>>  File
>>>>>>> "/usr/lib64/python2.6/site-packages/scipy/linalg/decomp_cholesky.py",
>>>>>>> line 66, in cholesky
>>>>>>>    c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a,
>>>>>>> clean=True)
>>>>>>>  File
>>>>>>> "/usr/lib64/python2.6/site-packages/scipy/linalg/decomp_cholesky.py",
>>>>>>> line 24, in _cholesky
>>>>>>>    raise LinAlgError("%d-th leading minor not positive definite" %
>>>>>>> info)
>>>>>>> numpy.linalg.linalg.LinAlgError: 1-th leading minor not positive
>>>>>>> definite
>>>>>>> ====================== *snap* ======================
>>>>>>>
>>>>>>> Fail!
>>>>>>>
>>>>>>> Am I missing a library, or is that routine broken?
>>>>>>
>>>>>> It looks like a bug if X is all zeros. If at least 1 element of X is
>>>>>> non-zero, it seems to work.
>>>>>>
>>>>>> Josef
>>>>>>
>>>>>>>
>>>>>>> Cheers,
>>>>>>> Nico
>>>>>>> _______________________________________________
>>>>>>> SciPy-User mailing list
>>>>>>> SciPy-User@scipy.org
>>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>>
>>>>>> _______________________________________________
>>>>>> SciPy-User mailing list
>>>>>> SciPy-User@scipy.org
>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>
>>>>>
>>>> _______________________________________________
>>>> SciPy-User mailing list
>>>> SciPy-User@scipy.org
>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User@scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User@scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User@scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
```