[SciPy-dev] integration of symeig in SciPy

[Nils Wagner]

On Fri, 24 Oct 2008 08:04:46 -0400
  alex griffing <argriffi at ncsu.edu> wrote:
> Tiziano Zito wrote:
>> Dear SciPy devs, 
>>
>> there has been some discussion about integrating symeig 
>>in SciPy.
>> symeig <http://mdp-toolkit.sourceforge.net/symeig.html> 
>>is a Python
>> wrapper for the LAPACK functions to solve the standard 
>>and
>> generalized eigenvalue problems for symmetric 
>>(hermitian) positive
>> definite matrices. It is superior to the scipy 
>>"linalg.eigh" wrapper
>> because it wraps *all* relevant LAPACK routines, 
>>enabling for
>> example to just return a subset of all 
>>egenvectors/eigenvalues,
>> which is a killer feature if you are dealing with 
>>huge-nonsparse
>> matrices. Some of these LAPACK routines are available in
>> scipy.lib.lapack.flapack and can be accessed through
>> scipy.lib.lapack.get_lapack_funcs . Some of them are 
>>still missing
>> in scipy, while symeig offers a unified interface to all 
>>relevant
>> LAPACK functions.  symeig has been dowloaded  more than 
>>1300 times
>> since its first public appearance on sourceforge (under 
>>the
>> mdp-toolkit package) in 2004. It features a complete 
>>unit test and a
>> telling doc-string. I am one of the authors of symeig 
>>and have been
>> contacted by some scipy devs to discuss the integration 
>>of symeig in
>> scipy. I am more than willing to help doing this. The 
>>most difficult
>> part (the LGPL license) has been addressed already: I've 
>>re-issued
>> symeig under a BSD license. Next step would be to adapt 
>>symeig code
>> and tests to the scipy framework. Some questions to you:
>>
>> - where would you like to have symeig? (I would suggest 
>>scipy.linalg)
>>
>> - should we change the signature to match or resamble 
>>that of other
>>   eigenproblem-related functions in scipy?
>>
>> - is there a responsible scipy dev for the linalg 
>>package I can bother
>>   with any questions I may have in the process?
>>
>> - do I get svn commit rights or should I work on a local 
>>copy and
>>   send a patch?
>>
>> let me know, 
>>
>> tiziano
>>   
> 
> This sounds like it would be very useful for scipy 
>although I am just a 
> user not a dev.  Could I use symeig to extract only the 
>eigenvector 
> corresponding to the largest eigenvalue, given a 
>non-sparse symmetric 
> real matrix whose eigenvalues are all positive except 
>for exactly one 
> that is zero?  The documentation for symeig says that it 
>works for 
> "symmetric (hermitian) definite positive matrices" but I 
>guess mine 
> would be a "symmetric (hermitian) semi-definite positive 
>matrix".  If 
> symeig wouldn't work, is there another way I could do 
>this without using 
> eigh to get the full eigendecomposition?
> 
> Thanks,
> Alex
> _______________________________________________
> Scipy-dev mailing list
> Scipy-dev at scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-dev
  
Here is an example where K is positive semidefinite.

from symeig import symeig
from numpy import diag, ones, r_

n = 10
M = 2.0*diag(ones(n))
#
# Symmetric positive-semidefinite
#
K = 
diag(r_[1.0,2*ones(n-2),1])-diag(ones(n-1),1)-diag(ones(n-1),-1)
K[0,0]   = 1.0
K[-1,-1] = 1.0
print M
print K

lo = 1
hi = 1
w,Z = symeig(K,M,range=(lo,hi))
print
print 'Smallest eigenvalue',w
print
print 'Rigid body mode'
print
print Z

lo = n
hi = n
w,Z = symeig(K,M,range=(lo,hi))
print
print 'Largest eigenvalue',w
print
print 'Corresponding eigenvector'
print
print Z

Nils