Hello list.. I've run into two SVD errors over the last few days. Both errors are identical in numpy/scipy.
I've submitted a ticket for the 1st problem (numpy ticket #990). Summary is: some builds of the lapack_lite module linking against system LAPACK (not the bundled dlapack_lite.o, etc) give a "LinAlgError: SVD did not converge" exception on my matrix. This error does occur using Mac's Accelerate framework LAPACK, and a coworker's Ubuntu LAPACK version. It does not seem to happen using ATLAS LAPACK (nor using Octave/Matlab on said Ubuntu)
Just today I've come across a negative singular value cropping up in an SVD of a different matrix. This error does occur on my ATLAS LAPACK based numpy, as well as on the Ubuntu setup. And once again, it does not happen in Octave/Matlab.
I'm using numpy 1.3.0.dev6336 -- don't know what the Ubuntu box is running.
Here are some npy files for the two different cases:
https://cirl.berkeley.edu/twiki/pub/User/MikeTrumpis/noconverge_operator.npy https://cirl.berkeley.edu/twiki/pub/User/MikeTrumpis/negsval_operator.npy
Mike
On Mon, Feb 2, 2009 at 18:21, mtrumpis@berkeley.edu wrote:
Hello list.. I've run into two SVD errors over the last few days. Both errors are identical in numpy/scipy.
I've submitted a ticket for the 1st problem (numpy ticket #990). Summary is: some builds of the lapack_lite module linking against system LAPACK (not the bundled dlapack_lite.o, etc) give a "LinAlgError: SVD did not converge" exception on my matrix. This error does occur using Mac's Accelerate framework LAPACK, and a coworker's Ubuntu LAPACK version. It does not seem to happen using ATLAS LAPACK (nor using Octave/Matlab on said Ubuntu)
These are almost certainly issues with the particular implementations of LAPACK that you are using. I don't think there is anything we can do from numpy or scipy to change this.
I ran into this problem a year or so ago. I suspect my messages to the list are in the archives somewhere. It is a known problem and involves a hard-coded maximum number of iterations in the SVD code. The problem is on the LaPack side. You can go in and change it, but then you have to recompile everything and rebuild Numpy, etc. etc. Not sure how easy/hard this is. I avoided it.
What I found that worked for me (depends on your numerical situation) is to take the original matrix you are trying to decompose, say A, and examine, instead, the SVD of A^T A. Then the singular values of that matrix are the square of the singular values of A. This worked for me, but my original matrix was square. Maybe that helped. Don't know. It's worth a try.
-- Lou Pecora, my views are my own.
--- On Mon, 2/2/09, mtrumpis@berkeley.edu mtrumpis@berkeley.edu wrote:
From: mtrumpis@berkeley.edu mtrumpis@berkeley.edu Subject: [Numpy-discussion] SVD errors To: numpy-discussion@scipy.org Date: Monday, February 2, 2009, 7:21 PM Hello list.. I've run into two SVD errors over the last few days. Both errors are identical in numpy/scipy.
I've submitted a ticket for the 1st problem (numpy ticket #990). Summary is: some builds of the lapack_lite module linking against system LAPACK (not the bundled dlapack_lite.o, etc) give a "LinAlgError: SVD did not converge" exception on my matrix. This error does occur using Mac's Accelerate framework LAPACK, and a coworker's Ubuntu LAPACK version. It does not seem to happen using ATLAS LAPACK (nor using Octave/Matlab on said Ubuntu)
Just today I've come across a negative singular value cropping up in an SVD of a different matrix. This error does occur on my ATLAS LAPACK based numpy, as well as on the Ubuntu setup. And once again, it does not happen in Octave/Matlab.
I'm using numpy 1.3.0.dev6336 -- don't know what the Ubuntu box is running.
Here are some npy files for the two different cases:
https://cirl.berkeley.edu/twiki/pub/User/MikeTrumpis/noconverge_operator.npy https://cirl.berkeley.edu/twiki/pub/User/MikeTrumpis/negsval_operator.npy
Mike
Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Mon, 02 Feb 2009 18:27:05 -0600, Robert Kern wrote:
On Mon, Feb 2, 2009 at 18:21, mtrumpis@berkeley.edu wrote:
Hello list.. I've run into two SVD errors over the last few days. Both errors are identical in numpy/scipy.
I've submitted a ticket for the 1st problem (numpy ticket #990). Summary is: some builds of the lapack_lite module linking against system LAPACK (not the bundled dlapack_lite.o, etc) give a "LinAlgError: SVD did not converge" exception on my matrix. This error does occur using Mac's Accelerate framework LAPACK, and a coworker's Ubuntu LAPACK version. It does not seem to happen using ATLAS LAPACK (nor using Octave/Matlab on said Ubuntu)
These are almost certainly issues with the particular implementations of LAPACK that you are using. I don't think there is anything we can do from numpy or scipy to change this.
Yes, this is almost certainly a LAPACK problem. If in doubt, you can test it with the following F90 program (be sure to link it against the same LAPACK as Numpy). Save the matrices with 'np.savetxt("foo.txt", x.ravel ())' and run './test < foo.txt'.
---- program test implicit none integer, parameter :: N = 128 double precision :: A(N*N), S(N), U(N,N), Vh(N,N) double precision, allocatable :: WORK(:) double precision :: tmp integer :: IWORK(8*N) integer :: INFO = 0, LWORK
read(*,*) A A = reshape(transpose(reshape(A, (/N, N/))), (/ N*N /)) call dgesdd('A', N, N, A, N, S, U, N, Vh, N, tmp, -1, IWORK, INFO) LWORK = tmp if (info .ne. 0) stop 'lwork query failed' write(*,*) 'lwork:', lwork allocate(WORK(LWORK)) call dgesdd('A', N, N, A, N, S, U, N, Vh, N, WORK, LWORK, IWORK, INFO) write(*,*) 'info:', INFO write(*,*) 'min(S):', minval(S) if (INFO .ne. 0) then write(*,*) ' -> SVD failed to converge' end if if (minval(S) .lt. 0) then write(*,*) ' -> negative singular value' end if end program ----
I played around with a C translation of that test program, and found that dgesvd (but not dgesdd) happens to converge and return all non-negative singular values for both operators I was having trouble with. I'm also looking at the Octave code just now, and I think they're using dgesvd also. Any one know why numpy uses dgesdd? speed?
Mike
On Tue, Feb 3, 2009 at 12:46 AM, Pauli Virtanen pav@iki.fi wrote:
Mon, 02 Feb 2009 18:27:05 -0600, Robert Kern wrote:
On Mon, Feb 2, 2009 at 18:21, mtrumpis@berkeley.edu wrote:
Hello list.. I've run into two SVD errors over the last few days. Both errors are identical in numpy/scipy.
I've submitted a ticket for the 1st problem (numpy ticket #990). Summary is: some builds of the lapack_lite module linking against system LAPACK (not the bundled dlapack_lite.o, etc) give a "LinAlgError: SVD did not converge" exception on my matrix. This error does occur using Mac's Accelerate framework LAPACK, and a coworker's Ubuntu LAPACK version. It does not seem to happen using ATLAS LAPACK (nor using Octave/Matlab on said Ubuntu)
These are almost certainly issues with the particular implementations of LAPACK that you are using. I don't think there is anything we can do from numpy or scipy to change this.
Yes, this is almost certainly a LAPACK problem. If in doubt, you can test it with the following F90 program (be sure to link it against the same LAPACK as Numpy). Save the matrices with 'np.savetxt("foo.txt", x.ravel ())' and run './test < foo.txt'.
program test implicit none integer, parameter :: N = 128 double precision :: A(N*N), S(N), U(N,N), Vh(N,N) double precision, allocatable :: WORK(:) double precision :: tmp integer :: IWORK(8*N) integer :: INFO = 0, LWORK
read(*,*) A A = reshape(transpose(reshape(A, (/N, N/))), (/ N*N /)) call dgesdd('A', N, N, A, N, S, U, N, Vh, N, tmp, -1, IWORK, INFO) LWORK = tmp if (info .ne. 0) stop 'lwork query failed' write(*,*) 'lwork:', lwork allocate(WORK(LWORK)) call dgesdd('A', N, N, A, N, S, U, N, Vh, N, WORK, LWORK, IWORK, INFO) write(*,*) 'info:', INFO write(*,*) 'min(S):', minval(S) if (INFO .ne. 0) then write(*,*) ' -> SVD failed to converge' end if if (minval(S) .lt. 0) then write(*,*) ' -> negative singular value' end if end program
Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
On Mon, Feb 9, 2009 at 16:25, M Trumpis mtrumpis@berkeley.edu wrote:
I played around with a C translation of that test program, and found that dgesvd (but not dgesdd) happens to converge and return all non-negative singular values for both operators I was having trouble with. I'm also looking at the Octave code just now, and I think they're using dgesvd also. Any one know why numpy uses dgesdd? speed?
Yes.
Robert Kern <robert.kern <at> gmail.com> writes:
On Mon, Feb 9, 2009 at 16:25, M Trumpis <mtrumpis <at> berkeley.edu> wrote:
I played around with a C translation of that test program, and found that dgesvd (but not dgesdd) happens to converge and return all non-negative singular values for both operators I was having trouble with. I'm also looking at the Octave code just now, and I think they're using dgesvd also. Any one know why numpy uses dgesdd? speed?
Yes.
Hello Everyone,
I am trying to write my own svd function to use the "dgesvd" method from lapack, but my problem is that I cannot find the "dgesvd" method from this import "from numpy.linalg import lapack_lite".
Do I have to install lapack on my system(Ubuntu 13.04 64-bit)? I want to get the same results as MATLAB 2011b which is why I am trying to write a slightly different svd function than the function from "numpy.linalg.svd".
I would appreciate your help, Thank you. Amir
On Tue, Jun 4, 2013 at 5:18 PM, Amir Mohammadi 183.amir@gmail.com wrote:
Hello Everyone,
I am trying to write my own svd function to use the "dgesvd" method from lapack, but my problem is that I cannot find the "dgesvd" method from this import "from numpy.linalg import lapack_lite".
Do I have to install lapack on my system(Ubuntu 13.04 64-bit)? I want to get the same results as MATLAB 2011b which is why I am trying to write a slightly different svd function than the function from "numpy.linalg.svd".
numpy does not wrap this LAPACK function. You will have to install LAPACK yourself and write your own extension module that wraps the DGESVD subroutine.
scipy also wraps *gesdd, and it would probably be a simple matter of copy-paste-fix to get it to wrap *gesvd too, since the interfaces are essentially the same.
https://github.com/scipy/scipy/blob/master/scipy/lib/lapack/flapack_esv.pyf....
-- Robert Kern