Re: LAPACK problem affecting full SVD, xGESDD (Octave 4.4)

[Tim Mitchell]

----- On Jan 24, 2019, at 7:11 PM, Rik <address@hidden> wrote:

On 01/24/2019 09:00 AM, address@hidden wrote:

Subject:

LAPACK problem affecting full SVD, xGESDD (Octave 4.4)

From:

Tim Mitchell <address@hidden>

Date:

01/23/2019 08:00 AM

To:

address@hidden

List-Post:

<mailto:address@hidden>

Precedence:

list

MIME-Version:

1.0

Message-ID:

<address@hidden>

Content-Type:

multipart/alternative; boundary="=_f90dbfe4-44b8-450f-88d0-51623ea9d5ab"

Message:

1

Dear Octave developers,

Last night I found a bad bug when computing a full SVD using xGESDD.  On certain problems, the computed singular values can be vastly different when also requesting singular vectors.

I've already opened an issue on LAPACK's Github page so I will just link to that for more details:

https://github.com/Reference-LAPACK/lapack/issues/316

Best regards,

Tim Mitchell

Could you also file a bug report about this at bugs.octave.org?  This is an upstream issue, but when users have trouble with Octave's SVD they are likely to start searching Octave's bug tracker, rather than LAPACK's bug tracker.

The specific driver to use for SVD in Octave is controlled by the svd_driver() function for which I quote the help text:

 -- VAL = svd_driver ()
 -- OLD_VAL = svd_driver (NEW_VAL)
 -- svd_driver (NEW_VAL, "local")
     Query or set the underlying LAPACK driver used by 'svd'.

     Currently recognized values are "gesdd" and "gesvd".  The default
     is "gesdd".

     When called from inside a function with the "local" option, the
     variable is changed locally for the function and any subroutines it
     calls.  The original variable value is restored when exiting the
     function.

     Algorithm Notes: The LAPACK library provides two routines for
     calculating the full singular value decomposition (left and right
     singular matrices as well as singular values).  When calculating
     just the singular values the following discussion is not relevant.

     The default routine use by Octave is the newer 'gesdd' which is
     based on a Divide-and-Conquer algorithm that is 5X faster than the
     alternative 'gesvd', which is based on QR factorization.  However,
     the new algorithm can use significantly more memory.  For an MxN
     input matrix the memory usage is of order O(min(M,N) ^ 2), whereas
     the alternative is of order O(max(M,N)). In general, modern
     computers have abundant memory so Octave has chosen to prioritize
     speed.

     In addition, there have been instances in the past where some input
     matrices were not accurately decomposed by 'gesdd'.  This appears
     to have been resolved with modern versions of LAPACK.  However, if
     certainty is required the accuracy of the decomposition can always
     be tested after the fact with

          [U, S, V] = svd (X);
          norm (X - U*S*V', "fro")

     See also: svd.

It appears that the last paragraph, "in the past some input matrices were not accurately decomposed by 'gesdd'", has returned.

Assuming, that the bug is verified to be real on the LAPACK site, a fix is unlikely to be available in distributions anytime soon (time to code fix, time to make a new LAPACK release, time for distributions to incorporate new library, time for desktop users to upgrade OS).

Octave itself is about to produce a 5.0 release.  We might minimize problems going forward by changing the default SVD calculation routine to the, apparently more stable, 'gesvd'.

--Rik

Hi Rik,

Thanks for the quick response.  I just posted an Octave bug report as requested.  There I also elaborated a bit how bad these errors can really be and why

[U,S,V]=svd(X);
norm(X-U*S*V',"fro")

is actually not a sufficient test.

I guess discussions can be continued on the Octave bug board or LAPACK issue page as appropriate.  As more users become more aware of this issue, I'm hoping they will send me any examples that trigger this issue.  It would be great to compile a good test set of challenging examples for GESDD.

Best regards,

Tim