Sorry, another sparse question: Is there anything is scipy which can find the eigenvalues and eigenvectors of a sparse matrix? -- David J. Grant Scientific Officer Intellectual Property D-Wave Systems Inc. tel: 604.732.6604 fax: 604.732.6614 ************************** CONFIDENTIAL COMMUNICATION This electronic transmission, and any documents attached hereto, is confidential. The information is intended only for use by the recipient named above. If you have received this electronic message in error, please notify the sender and delete the electronic message. Any disclosure, copying, distribution, or use of the contents of information received in error is strictly prohibited.
David Grant wrote:
Sorry, another sparse question: Is there anything is scipy which can find the eigenvalues and eigenvectors of a sparse matrix?
Unfortunately not yet. I have added an entry in the Wiki page (Support for ARPACK in scipy) http://www.scipy.org/wikis/featurerequests/PrioritiesAndWishlist There are some external tools http://jrfonseca.dyndns.org/work/phd/ http://sourceforge.net/projects/pysparse/ http://people.web.psi.ch/geus/pyfemax/ Nils
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
Nils Wagner wrote:
David Grant wrote:
Sorry, another sparse question: Is there anything is scipy which can find the eigenvalues and eigenvectors of a sparse matrix?
Unfortunately not yet. I have added an entry in the Wiki page (Support for ARPACK in scipy)
http://www.scipy.org/wikis/featurerequests/PrioritiesAndWishlist
There are some external tools
http://jrfonseca.dyndns.org/work/phd/ http://sourceforge.net/projects/pysparse/ http://people.web.psi.ch/geus/pyfemax/
I just looked at pysparse again. It looks like complex numbers cannot be stored. Am I correct? David
Nils Wagner wrote:
David Grant wrote:
Sorry, another sparse question: Is there anything is scipy which can find the eigenvalues and eigenvectors of a sparse matrix?
Unfortunately not yet. I have added an entry in the Wiki page (Support for ARPACK in scipy)
http://www.scipy.org/wikis/featurerequests/PrioritiesAndWishlist
There are some external tools
Does anyone know how much work it would be to get ARPACK support in scipy? Even just minimal support. This is totally new ground for me, so I'd like to know if I'm biting off more than I can chew. -- David J. Grant http://www.davidgrant.ca:81
Does anyone know how much work it would be to get ARPACK support in scipy? Even just minimal support. This is totally new ground for me, so I'd like to know if I'm biting off more than I can chew.
Not a lot for simple support if you use f2py. My suggestion would be to write and f2py interface for the desired routines and we can build up from there. -Travis
Travis Oliphant wrote:
Does anyone know how much work it would be to get ARPACK support in scipy? Even just minimal support. This is totally new ground for me, so I'd like to know if I'm biting off more than I can chew.
Not a lot for simple support if you use f2py.
My suggestion would be to write an f2py interface for the desired routines and we can build up from there.
I assume your suggestion does not involve using this at all: http://jrfonseca.dyndns.org/work/phd/ ? I believe Josa Fonseca was going to use f2py, but he didn't. Here is from his web page: "These bindings weren't generated with any automatic binding generation tool. Even though I initially tried both PyFortran and f2py, both showed to be inappropriate to handle the specificity of the ARPACK API. ARPACK uses a /reverse communication interface/ where basically the API sucessively returns to caller which must take some update steps, and re-call the API with most arguments untouched. The intelligent (and silent) argument conversions made by the above tools made very difficult to implement and debug the most simple example. Also, for large-scale problems we wouldn't want any kind of array conversion/transposing happening behind the scenes as that would completely kill performance." "Therefore the bindings are faithfull to the original FORTRAN API and are not very friendly as is. Nevertheless a Python wrapper to these calls can easily be made, where all kind of type conversions and dummy-safe actions can be made." So who is right? David
On Thu, 28 Oct 2004, David Grant wrote:
Travis Oliphant wrote:
Does anyone know how much work it would be to get ARPACK support in scipy? Even just minimal support. This is totally new ground for me, so I'd like to know if I'm biting off more than I can chew.
Not a lot for simple support if you use f2py.
My suggestion would be to write an f2py interface for the desired routines and we can build up from there.
I assume your suggestion does not involve using this at all: http://jrfonseca.dyndns.org/work/phd/ ? I believe Josa Fonseca was going to use f2py, but he didn't.
Here is from his web page:
"These bindings weren't generated with any automatic binding generation tool. Even though I initially tried both PyFortran and f2py, both showed to be inappropriate to handle the specificity of the ARPACK API. ARPACK uses a /reverse communication interface/ where basically the API sucessively returns to caller which must take some update steps, and re-call the API with most arguments untouched. The intelligent (and silent) argument conversions made by the above tools made very difficult to implement and debug the most simple example. Also, for large-scale problems we wouldn't want any kind of array conversion/transposing happening behind the scenes as that would completely kill performance."
I cannot agree with these conclusions, at least when using f2py. The array conversion can be easily avoided behind the scenes (because the main loop is in Fortran) and transposing is not even an issue as user-specific parts of ARPACK API contains only 1-dimensional arrays.
"Therefore the bindings are faithfull to the original FORTRAN API and are not very friendly as is. Nevertheless a Python wrapper to these calls can easily be made, where all kind of type conversions and dummy-safe actions can be made."
So who is right?
I looked into ARPACK API and I think f2py could be effectively used to wrap ARPACK. For example, when considering SIMPLE/dssimp.f case then I would write a Fortran subroutine that contains dsaupd main-loop and call to dseupd (the subroutine av would be provided as a callback argument) and then wrap this subroutine with f2py. Pearu
Pearu Peterson wrote:
On Thu, 28 Oct 2004, David Grant wrote:
Travis Oliphant wrote:
Does anyone know how much work it would be to get ARPACK support in scipy? Even just minimal support. This is totally new ground for me, so I'd like to know if I'm biting off more than I can chew.
Not a lot for simple support if you use f2py.
My suggestion would be to write an f2py interface for the desired routines and we can build up from there.
I assume your suggestion does not involve using this at all: http://jrfonseca.dyndns.org/work/phd/ ? I believe Josa Fonseca was going to use f2py, but he didn't.
Here is from his web page:
"These bindings weren't generated with any automatic binding generation tool. Even though I initially tried both PyFortran and f2py, both showed to be inappropriate to handle the specificity of the ARPACK API. ARPACK uses a /reverse communication interface/ where basically the API sucessively returns to caller which must take some update steps, and re-call the API with most arguments untouched. The intelligent (and silent) argument conversions made by the above tools made very difficult to implement and debug the most simple example. Also, for large-scale problems we wouldn't want any kind of array conversion/transposing happening behind the scenes as that would completely kill performance."
I cannot agree with these conclusions, at least when using f2py. The array conversion can be easily avoided behind the scenes (because the main loop is in Fortran) and transposing is not even an issue as user-specific parts of ARPACK API contains only 1-dimensional arrays.
I am in agreement with Pearu. Look at the iterative package in scipy.linalg. The underlying Fortran code uses reverse communication and this was managed with f2py-generated wrappers.
"Therefore the bindings are faithfull to the original FORTRAN API and are not very friendly as is. Nevertheless a Python wrapper to these calls can easily be made, where all kind of type conversions and dummy-safe actions can be made."
So who is right?
I looked into ARPACK API and I think f2py could be effectively used to wrap ARPACK. For example, when considering SIMPLE/dssimp.f case then I would write a Fortran subroutine that contains dsaupd main-loop and call to dseupd (the subroutine av would be provided as a callback argument) and then wrap this subroutine with f2py.
My suggestion is to do something very similar as is done in scipy/linalg/iterative.py wherein the object passed to the Python-level code simply requires a matvec method (like all sparse arrays have), or is a normal 2-d array. Then, call python-wrapped dsaupd in a while loop which does the required thing after returning from every dsaupd call. This is the most flexible interface and deserves a place in scipy. Pearu's suggestion is also useful and may be a bit faster (though I'm not sure it would be noticeably faster). -Travis
On Thu, 28 Oct 2004, Travis Oliphant wrote:
I looked into ARPACK API and I think f2py could be effectively used to wrap ARPACK. For example, when considering SIMPLE/dssimp.f case then I would write a Fortran subroutine that contains dsaupd main-loop and call to dseupd (the subroutine av would be provided as a callback argument) and then wrap this subroutine with f2py.
My suggestion is to do something very similar as is done in scipy/linalg/iterative.py wherein the object passed to the Python-level code simply requires a matvec method (like all sparse arrays have), or is a normal 2-d array. Then, call python-wrapped dsaupd in a while loop which does the required thing after returning from every dsaupd call. This is the most flexible interface and deserves a place in scipy.
Pearu's suggestion is also useful and may be a bit faster (though I'm not sure it would be noticeably faster).
My suggestion was not to make things faster (though it might turn out to be slightly so) but more robustly to share common data between dsaupd and dseupd calls: the second part of dseupd arguments should be exactly the same as dsaupd plus the input arguments determining the task must be the same. Though wrapping both dsaupd and dseupd to Python gives more flexibilty but I am not sure that will be ever needed nor would avoid incorrect usage of the ARPACK wrappers. Just my 2 Estonian cents, Pearu
Does anyone know how much work it would be to get ARPACK support in scipy? Even just minimal support. This is totally new ground for me, so I'd like to know if I'm biting off more than I can chew.
After looking over ARPACK, it looks to me that you could get a functional ARPACK interface in 1-2 days, using f2py. Then, you could look at the linalg/iterative methods that interface with reverse-communication fortran code to get an idea how you might write simpler Python wrappers to the raw wrappers produced by f2py. I don't think this would be hard at all. It would be a good way for you to learn f2py, if you are not already familiar with it. Also, you might notice how iterative and sparse f2py interfaces use a (very simple) templating mechanism (in scipy_distutils/from_template.py) so that interfaces for all four precisions are maintained with a single code-base. But, if you produce an interface for the precision you care about, I could handle converting it to a general-purpose interface for all precisions. Is the ARPACK code distributible in Python? -Travis
Travis Oliphant wrote:
Does anyone know how much work it would be to get ARPACK support in scipy? Even just minimal support. This is totally new ground for me, so I'd like to know if I'm biting off more than I can chew.
After looking over ARPACK, it looks to me that you could get a functional ARPACK interface in 1-2 days, using f2py. Then, you could look at the linalg/iterative methods that interface with reverse-communication fortran code to get an idea how you might write simpler Python wrappers to the raw wrappers produced by f2py.
So use f2py. Then look at the stuff here: *http://tinyurl.com/3uzmy* to figure out how to get reverse-communication to work, then write higher-level code which talks to the python code generated by f2py?
I don't think this would be hard at all. It would be a good way for you to learn f2py, if you are not already familiar with it. Also, you might notice how iterative and sparse f2py interfaces use a (very simple) templating mechanism (in scipy_distutils/from_template.py) so that interfaces for all four precisions are maintained with a single code-base.
But, if you produce an interface for the precision you care about, I could handle converting it to a general-purpose interface for all precisions.
ok
Is the ARPACK code distributible in Python?
ARPACK is "freely distributable" or "public domain". I've seen both of these used to describe ARPACK's license or lack-therof. -- David J. Grant Scientific Officer Intellectual Property D-Wave Systems Inc. tel: 604.732.6604 fax: 604.732.6614 ************************** CONFIDENTIAL COMMUNICATION This electronic transmission, and any documents attached hereto, is confidential. The information is intended only for use by the recipient named above. If you have received this electronic message in error, please notify the sender and delete the electronic message. Any disclosure, copying, distribution, or use of the contents of information received in error is strictly prohibited.
David Grant wrote:
Sorry, another sparse question: Is there anything is scipy which can find the eigenvalues and eigenvectors of a sparse matrix?
Having read all the followups, it seems to me that I must be asking a dumb question -- but won't linalg.svd do this? Begging for education, really... Nick
Nick Arnett wrote:
David Grant wrote:
Sorry, another sparse question: Is there anything is scipy which can find the eigenvalues and eigenvectors of a sparse matrix?
Having read all the followups, it seems to me that I must be asking a dumb question -- but won't linalg.svd do this?
Begging for education, really...
An SVD is not the same as an eigen decomposition. http://mathworld.wolfram.com/EigenDecomposition.html http://mathworld.wolfram.com/SingularValueDecomposition.html In any case, all of the linalg.* functions only operate on dense arrays, not sparse matrices. -- Robert Kern rkern@ucsd.edu "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
Robert Kern wrote:
An SVD is not the same as an eigen decomposition.
http://mathworld.wolfram.com/EigenDecomposition.html http://mathworld.wolfram.com/SingularValueDecomposition.html
That dawned on me not too long after posting. Fuzzy brain this week, a birth and a death in the family a few days apart.
In any case, all of the linalg.* functions only operate on dense arrays, not sparse matrices.
Ummm, now I'm thinking I'm in over my head. SVD is a matrix operation isn't it? I was imagining that an array of mostly zeros would be the equivalent of a sparse matrix, so that SVD would apply to it for what I'm doing, which is related to latent semantic analysis. Am I way off-base here? My understanding of Scipy isn't great, compounded by a less-than-great grasp of the mathematics. Learning as go, in large part. Any help gratefully accepted. Nick
Nick Arnett wrote:
Robert Kern wrote:
An SVD is not the same as an eigen decomposition.
http://mathworld.wolfram.com/EigenDecomposition.html http://mathworld.wolfram.com/SingularValueDecomposition.html
That dawned on me not too long after posting. Fuzzy brain this week, a birth and a death in the family a few days apart.
In any case, all of the linalg.* functions only operate on dense arrays, not sparse matrices.
Ummm, now I'm thinking I'm in over my head. SVD is a matrix operation isn't it? I was imagining that an array of mostly zeros would be the equivalent of a sparse matrix, so that SVD would apply to it for what I'm doing, which is related to latent semantic analysis.
The SVD of A is constructed using the eigen decompostion of A A^H and A^H A which have the same real-valued eigenvalues (though one may have some zero-valued eigenvalues the other doesn't) and unitary eigenvector matrices: The decompostion provides A = U D V^H where U U^H =I and V V^H = I so that AA^H = U D^2 U^H and A^H A = V D^2 V^H showing that U are the eigenvectors of A A^H corresponding to non-zero eigenvalues V are the eigenvectors of A^H A corresponding to non-zero eigenvalues D^2 are the non-zero eigenvalues of both So, while the SVD is not an eigenvector decomposition it is related to one. -Travis
Travis Oliphant wrote:
So, while the SVD is not an eigenvector decomposition it is related to one.
Right... but I still don't understand the statement, "In any case, all of the linalg.* functions only operate on dense arrays, not sparse matrices." Why would linalg.svd not operate on a sparse matrix? I was working from the example (below) from your Scipy tutorial, in fact. Would the results not be meaningful if the matrix is sparse?
A = mat('[1 3 2; 1 2 3]') M,N = A.shape U,s,Vh = linalg.svd(A) Sig = mat(diagsvd(s,M,N)) U, Vh = mat(U), mat(Vh) print U Matrix([[-0.7071, -0.7071], [-0.7071, 0.7071]]) print Sig Matrix([[ 5.1962, 0. , 0. ], [ 0. , 1. , 0. ]]) print Vh Matrix([[-0.2722, -0.6804, -0.6804], [-0. , -0.7071, 0.7071], [-0.9623, 0.1925, 0.1925]]) print A Matrix([[1, 3, 2], [1, 2, 3]]) print U*Sig*Vh Matrix([[ 1., 3., 2.], [ 1., 2., 3.]])
Thanks! Nick
On Wed, 17 Nov 2004, Nick Arnett wrote:
Travis Oliphant wrote:
So, while the SVD is not an eigenvector decomposition it is related to one.
Right... but I still don't understand the statement, "In any case, all of the linalg.* functions only operate on dense arrays, not sparse matrices."
The statement means that the current *implementation* of linalg functions cannot operate on sparse matrices. Mathematically there is no difference between sparse or dense matrices. Sparse or dense are notions of matrix representation.
Why would linalg.svd not operate on a sparse matrix? I was working from the example (below) from your Scipy tutorial, in fact. Would the results not be meaningful if the matrix is sparse?
<snip> Of course they are. Just with the current implementation one must transform a sparse matrix to a full matrix (that may contain lots of zeros) before applying linalg functions. I believe that in future linalg functions will support sparse matrices as well. Pearu
Pearu Peterson wrote:
Of course they are. Just with the current implementation one must transform a sparse matrix to a full matrix (that may contain lots of zeros) before applying linalg functions.
Okay, now I get it -- I was imagining that there was some mathematical difference between dense and sparse matrices that I wasn't aware of... What we're talking about is the representation, if I understand correctly. And... whew. Nick
Nick Arnett wrote:
Travis Oliphant wrote:
So, while the SVD is not an eigenvector decomposition it is related to one.
Right... but I still don't understand the statement, "In any case, all of the linalg.* functions only operate on dense arrays, not sparse matrices."
Why would linalg.svd not operate on a sparse matrix? I was working from the example (below) from your Scipy tutorial, in fact. Would the results not be meaningful if the matrix is sparse?
Sparse matrices and dense matrices are very different objects and linalg.svd is a wrapper around LAPACK which only works on dense matrices. Getting all linear algebra operations working for sparse matrices is a very tall order and has not been done yet. -Travis
participants (8)
-
David Grant -
David Grant -
Nick Arnett -
Nick Arnett -
Nils Wagner -
Pearu Peterson -
Robert Kern -
Travis Oliphant