SVD and QR for non-square matrices not matching Eigen

[cdeterman]

I am experimenting with the different solvers (which I know themselves experimental). Of particular interest is SVD. I decided to compare the results to the Eigen C++ library which as you know is a fantastic linear algebra library (hence why its integrated here). The results between ViennaCL and Eigen match when provided square matrices but deviate when provided non-square matrices. For example, the following outputs differ noticeably. The order of the columns is not really important but the actual numbers don't line up between the two libraries. Is this a known bug?

#include <Eigen/Core>
#include <Eigen/SVD>
#include <iostream>

#define VIENNACL_WITH_OPENCL 1
#define VIENNACL_WITH_EIGEN 1

#include "viennacl/vector.hpp"
#include "viennacl/matrix.hpp"
#include "viennacl/linalg/svd.hpp"

using namespace Eigen;
using std::cout;

int main()
{
    // initializers
    MatrixXd C;
    viennacl::matrix<double> vcl_C = viennacl::zero_matrix<double>(4,6);   
    viennacl::matrix<double> vcl_U = viennacl::zero_matrix<double>(4, 4);
    viennacl::matrix<double> vcl_V = viennacl::zero_matrix<double>(6, 6);

    // random data
    C.setRandom(4,6);

    // copy eigen to viennacl
    copy(C, vcl_C);

    // eigen SVD
    JacobiSVD<MatrixXd> svd( C, ComputeThinU | ComputeThinV);

    std::cout << "diag" << std::endl;
    std::cout << svd.singularValues() << std::endl; 

    std::cout << "U" << std::endl;
    std::cout << svd.matrixU() << std::endl;

    std::cout << "V" << std::endl;
    std::cout << svd.matrixV() << std::endl;

    // viennacl svd

    viennacl::linalg::svd(vcl_C, vcl_U, vcl_V);

    viennacl::vector_base<double> D(vcl_C.handle(), std::min(vcl_C.size1(), vcl_C.size2()), 0, vcl_C.internal_size2() + 1);

    std::cout << "diag" << std::endl;
    std::cout << D << std::endl;

    std::cout << "U" << std::endl;
    std::cout << vcl_U << std::endl;

    std::cout << "V" << std::endl;
    std::cout << vcl_V << std::endl;

    return 0;
}

[cdeterman]

Not meaning to be pushy but have you been able to have a look at this? I'm not necessarily asking for a major fix, I just want to know if this is expected behavior with the current version. Thanks

[karlrupp]

Sorry for the late response, I was on vacation. It looks like an issue on our front: ViennaCL's results are wrong as soon as the number of columns in C is larger than the number of rows. It is know a known bug, thank you for reporting.

(I suspect that the SVD code, which was written at an early stage, does not correctly account for the number of rows and columns in one of the worker kernels)

[cdeterman]

@karlrupp the same thing appears to be happening with QR. Would you prefer a separate issue created for that or should I change the title of this one to address both SVD and QR.

[karlrupp]

Please update the issue title and specify whether you are referring to the QR method for eigenvalues (which I assume) or to a QR factorization.

[cdeterman]

I am referring to the QR factorization. Below is a reproducible example showing the difference between the Eigen library and ViennaCL. It's primarily the R matrix that isn't filled correctly. The Q matrices appear to match correctly.

#include <Eigen/Core>
#include <Eigen/Dense>
#include <iostream>

#define VIENNACL_WITH_OPENCL 1
#define VIENNACL_WITH_EIGEN 1

#include "viennacl/vector.hpp"
#include "viennacl/matrix.hpp"
#include "viennacl/linalg/qr.hpp"

using namespace Eigen;
using std::cout;

int main()
{
    // initializers
    MatrixXd C;
    viennacl::matrix<double> vcl_C = viennacl::zero_matrix<double>(4,6);   
    viennacl::matrix<double> vcl_R = viennacl::zero_matrix<double>(4, 6); 
    viennacl::matrix<double> vcl_Q = viennacl::zero_matrix<double>(4, 4); 

    // random data
    C.setRandom(4,6);

    // copy eigen to viennacl
    copy(C, vcl_C);

    // eigen QR
    HouseholderQR<MatrixXd> qr(C);
    MatrixXd R = qr.matrixQR().triangularView<Upper>();
    MatrixXd Q = qr.householderQ();

    std::cout << "R" << std::endl;
    std::cout << R << std::endl;

    std::cout << "Q" << std::endl;
    std::cout << Q << std::endl;

    // viennacl QR

    std::vector<double> betas = viennacl::linalg::inplace_qr(vcl_C);

    viennacl::linalg::recoverQ(vcl_C, betas, vcl_Q, vcl_R);

    std::cout << "R" << std::endl;
    std::cout << vcl_R << std::endl;

    std::cout << "Q" << std::endl;
    std::cout << vcl_Q << std::endl;

    return 0;
}

[karlrupp]

Thanks for the example, this will help me fixing it.

[cdeterman]

@karlrupp I was wondering what the status of this issue is for the 1.8.0 release. If I recall, you said you expected 1.8.0 to be released by the end of March but I could be mistaken. Any updates you can provide would be appreciated. Thanks

[karlrupp]

@cdeterman It is still scheduled for 1.8.0, yes. A release of 1.8.0 in March is unlikely, but some time in April or May will do. Regardless, you will get notified when the issue is fixed and closed :-)

[cdeterman]

@karlrupp Sorry to bug again as I know you will notify when the issue is fixed here but I need to ask for a new estimated timeline. I have users asking about this functionality and I would like to know if I should be waiting for this fix or if I should explore some custom OpenCL kernels in the interim. I just don't want to pursue a bunch of work if this issue is intended to be fixed in the 'near future'. Thanks.

[karlrupp]

@cdeterman yeah, unfortunately my estimates were too optimistic. Take end of September or early October as a new estimate. Also, I'll look into this issue in a couple of days, hopefully I can spot the issue (I do know the QR factorization code, but I don't know all details of the SVD code)