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

Nico Schlömer nico.schloemer@gmail....
Wed Oct 13 09:58:03 CDT 2010


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
>
>
>
> 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
>>
>


More information about the SciPy-User mailing list