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

[josef.pktd at gmail.com]

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
>