[Numpy-discussion] SVD errors
M Trumpis
mtrumpis@berkeley....
Mon Feb 9 16:25:54 CST 2009
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
>
More information about the Numpy-discussion
mailing list