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

Robert Cimrman cimrman3 at ntc.zcu.cz
Wed Oct 13 13:09:35 EDT 2010 

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 at gmail.com> wrote:
>> On Wed, Oct 13, 2010 at 10:58 AM, Nico Schlömer
>> <nico.schloemer at 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 at 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 at gmail.com> wrote:
>>>>> On Wed, Oct 13, 2010 at 10:21 AM, Nico Schlömer
>>>>> <nico.schloemer at 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 at scipy.org
>>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>>
>>>>> _______________________________________________
>>>>> SciPy-User mailing list
>>>>> SciPy-User at scipy.org
>>>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>>>
>>>>
>>> _______________________________________________
>>> SciPy-User mailing list
>>> SciPy-User at scipy.org
>>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>>
>> _______________________________________________
>> SciPy-User mailing list
>> SciPy-User at scipy.org
>> http://mail.scipy.org/mailman/listinfo/scipy-user
>>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>

More information about the SciPy-User mailing list