Nico Schlömer | 13 Oct 21:42 2010
Picon

Re: scipy.sparse.linalg.lobpcg broken?

Aaaand another thing:
The matrix I'm looking at is in fact complex-valued, but Hermitian,
for which LOBPCG should work. However:

========================== *snip* ==========================
/opt/scipy/0.8.0/lib/python/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py:457:
ComplexWarning: Casting complex values to real discards the imaginary
part
  eigBlockVector = np.asarray( eigBlockVector[:,ii].astype( np.float64 ) )
========================== *snap* ==========================

And the computed result is not correct. :/

--Nico

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

Gmane