[Numpy-discussion] SVD errors

[M Trumpis]

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 at iki.fi> wrote:
> Mon, 02 Feb 2009 18:27:05 -0600, Robert Kern wrote:
>> On Mon, Feb 2, 2009 at 18:21,  <mtrumpis at 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 at scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion
>