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

Nico Schlömer nico.schloemer@gmail....
Wed Oct 13 11:57:57 CDT 2010


Possibly the devs know -- how to get back to them?

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


More information about the SciPy-User mailing list