SVD fails when computing singular values only

Issue #409 wontfix
Nils Vu created an issue

blaze::svd fails for a 3x3 matrix when computing only the singular values, but works when computing left and right singular vectors. For concreteness, this is the matrix:

H =
( 5.8170601797986201120e+00 1.1098039266568445527e+01 -9.8623087318593397745e-17 )
( 1.1098039266568445527e+01 4.4929100274260804326e+01 2.7160992219672966286e+01 )
( 0.0000000000000000000e+00 2.7160992219672962733e+01 9.1466495154621867414e+01 )

This snippet fails:

blaze::DynamicVector<double, blaze::columnVector> s;
blaze::svd(H, s);

This snippet works:

blaze::DynamicMatrix<double, blaze::rowMajor> U;
blaze::DynamicVector<double, blaze::columnVector> s;
blaze::DynamicMatrix<double, blaze::rowMajor> V;
blaze::svd(H, U, s, V);

For some reason, this also works:

blaze::DynamicVector<double, blaze::columnVector> s;
const auto s = blaze::svd(H);
printf("s = %s", s);

I tested this with Blaze 3.8. Please let me know if I can provide more information.

Comments (8)

  1. Klaus Iglberger

    Hi Nils!

    Thanks for reporting this potential defect. Thanks also for your minimum example. Unfortunately we’re not able to reproduce the problem. We used the following code to compare the three different approaches:

    #include <iostream>
    #include <cstdlib>
    #include <blaze/Blaze.h>
    
    int main()
    {
       blaze::DynamicMatrix<double, blaze::rowMajor> U;
       blaze::DynamicVector<double, blaze::columnVector> s1, s2, s3;
       blaze::DynamicMatrix<double, blaze::rowMajor> V;
    
       blaze::DynamicMatrix<double, blaze::rowMajor> H
          { { 5.8170601797986201120e+00, 1.1098039266568445527e+01, -9.8623087318593397745e-17 }
          , { 1.1098039266568445527e+01, 4.4929100274260804326e+01,  2.7160992219672966286e+01 }
          , { 0.0000000000000000000e+00, 2.7160992219672962733e+01,  9.1466495154621867414e+01 } };
    
       blaze::svd(H, U, s1, V);
       blaze::svd(H, s2);
       s3 = svd(H);
    
       std::cout << "\n"
                 << "Result 1:\n"
                 << "  s1 = " << trans(s1) << "\n"
                 << "Result 2:\n"
                 << "  s2 = " << trans(s2) << "\n"
                 << "Result 3:\n"
                 << "  s3 = " << trans(s3) << "\n"
                 << "\n";
    
       return EXIT_SUCCESS;
    }
    

    The output:

    Result 1:
      s1 = ( 104.185 35.7887 2.23896 )
    Result 2:
      s2 = ( 104.185 35.7887 2.23896 )
    Result 3:
      s3 = ( 104.185 35.7887 2.23896 )
    

    In order to understand what causes your problem, we unfortunately need more information.

    First, could you please give us the information which LAPACK library you’re working with?

    Second, as you can see in <blaze/math/dense/SVD.h>, the first two calls end up in a call to the LAPACK function gesdd() (see the wiki; the third one also does, but with a detour). Instead of calling the svd() functions, could you please directly use the gesdd() functions? In order to see if the problem is within Blaze or the LAPACK implementation, it would indeed help if you could use the second gesdd() function, which is merely a wrapper around the according LAPACK kernel (see <blaze/math/lapack/clapack/gesdd.h>). Unfortunately this is much more inconvenient (14 arguments!), but if you use column-major matrices instead of row-major matrices you shouldn’t have any problems to make it work.

    Hopefully this will help to understand the problem. Thanks again, we appreciate your help,

    Best regards,

    Klaus!

  2. Nils Vu reporter

    Hi Klaus, thanks for checking this issue! I’ll test calling into LAPACK directly and report back here.

  3. Klaus Iglberger

    Hi Nils!

    Is it correct to assume that the issue is no longer relevant and that the problem has been solved?

    Best regards,

    Klaus!

  4. Nils Vu reporter

    Hi Klaus, I have since worked around the issue, but have not yet investigated it any further.

  5. Klaus Iglberger

    Hi Nils!

    I will close the issue as there is no clear indication of an error. In case you encounter new "evidence", please reopen this issue.

    Best regards,

    Klaus!

  6. Log in to comment