Home
Reading
Searching
Subscribe
Sponsors
Statistics
Posting
Contact
Spam
Lists
Links
About
Hosting
Filtering
Features
Download
Marketing
Archives
FAQ
Blog


Aaaand another thing:
The matrix I'm looking at is in fact complexvalued, 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
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
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 nonzero 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, wrote:
>>>>
>>>> On Wed, Oct 13, 2010 at 10:58 AM, Nico Schlömer
>>>> 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
>>>>> 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+1e6*np.random.randn(10)), np.random.randn(10,1),
>>>> verbosityLevel=10)
>>>> lobpcg( sparse.csr_matrix(np.diag(1+1e4*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.92296269e16]
>>>> (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, wrote:
>>>>>>>
>>>>>>> On Wed, Oct 13, 2010 at 10:21 AM, Nico Schlömer
>>>>>>> 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
>>>>>>>> lobpcg( A, X )
>>>>>>>> File
>>>>>>>>
"/usr/lib64/python2.6/sitepackages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py",
>>>>>>>> line 304, in lobpcg
>>>>>>>> blockVectorX, blockVectorBX = b_orthonormalize( B,
blockVectorX )
>>>>>>>> File
>>>>>>>>
"/usr/lib64/python2.6/sitepackages/scipy/sparse/linalg/eigen/lobpcg/lobpcg.py",
>>>>>>>> line 130, in b_orthonormalize
>>>>>>>> gramVBV = sla.cholesky( gramVBV )
>>>>>>>> File
>>>>>>>>
"/usr/lib64/python2.6/sitepackages/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/sitepackages/scipy/linalg/decomp_cholesky.py",
>>>>>>>> line 24, in _cholesky
>>>>>>>> raise LinAlgError("%dth leading minor not positive definite"
%
>>>>>>>> info)
>>>>>>>> numpy.linalg.linalg.LinAlgError: 1th 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
>>>>>>> nonzero, it seems to work.
>>>>>>>
>>>>>>> Josef
>>>>>>>
>>>>>>>>
>>>>>>>> Cheers,
>>>>>>>> Nico
>>>>>>>> _______________________________________________
>>>>>>>> SciPyUser mailing list
>>>>>>>> [email protected]
>>>>>>>> http://mail.scipy.org/mailman/listinfo/scipyuser
>>>>>>>>
>>>>>>> _______________________________________________
>>>>>>> SciPyUser mailing list
>>>>>>> [email protected]
>>>>>>> http://mail.scipy.org/mailman/listinfo/scipyuser
>>>>>>>
>>>>>>
>>>>> _______________________________________________
>>>>> SciPyUser mailing list
>>>>> [email protected]
>>>>> http://mail.scipy.org/mailman/listinfo/scipyuser
>>>>>
>>>> _______________________________________________
>>>> SciPyUser mailing list
>>>> [email protected]
>>>> http://mail.scipy.org/mailman/listinfo/scipyuser
>>>>
>>> _______________________________________________
>>> SciPyUser mailing list
>>> [email protected]
>>> http://mail.scipy.org/mailman/listinfo/scipyuser
>>>
>>
>> _______________________________________________
>> SciPyUser mailing list
>> [email protected]
>> http://mail.scipy.org/mailman/listinfo/scipyuser
>>
>>
>

