I am writing a wrapper to Eigen for my personal use and I encountered the following weird behavior:

Code: Select all

void get_QR(MatrixXd A, MatrixXd& Q, MatrixXd& R) {
        HouseholderQR<MatrixXd> qr(A);
        Q  =   qr.householderQ()*(MatrixXd::Identity(A.rows(),A.cols()));
        R  =   qr.matrixQR().block(0,0,A.cols(),A.cols()).triangularView<Upper>();
}

Code: Select all

void get_QR(double* A, int m, int n, double*& Q, double*& R) {
        //      Maps the double to MatrixXd.
        Map<MatrixXd> A_E(A, m, n);
       
        //      Obtains the QR of A_E.
        MatrixXd Q_E, R_E;
        get_QR(A_E, Q_E, R_E);
       
        //      Maps the MatrixXd to double.
        Q       =       Q_E.data();
        R       =       R_E.data();
}

Code: Select all

        srand(time(NULL));

        int m   =       atoi(argv[1]);
        int n   =       atoi(argv[2]);

//      Check the double version.
        double* A       =       new double[m*n];
        double* Q       =       new double[m*n];
        double* R       =       new double[n*n];

        double RANDMAX  =       double(RAND_MAX);
        //      Initialize A as a random matrix.
        for (int index=0; index<m*n; ++index) {
                A[index]        =       rand()/RANDMAX;
        }

        get_QR(A, m, n, Q, R);
        std::cout << Q[0] << std::endl;

//      Check the MatrixXd version.
        Map<MatrixXd> A_E(A, m, n);
        MatrixXd Q_E(m,n), R_E(n,n);
        get_QR(A_E, Q_E, R_E);
        std::cout << Q[0] << std::endl;

I get different values of Q[0]. For instance, I get "-0.421857" and "-1.49563".

Thanks

This question has been answered here: http://stackoverflow.com/questions/2210 ... sing-eigen

ggael, I added couple more questions in the comments over there. I would greatly appreciate if you could provide an answer for them as well.

@ggael If I implement what you have in the post on stack overflow, I then get the following:

If I call the function get_QR(Ref<const MatrixXd> A, Ref<MatrixXd> Q, Ref<MatrixXd> R), from outside, where I declare Q and R as MatrixXd Q, MatrixXd R, I get the following error: Assertion failed: (rows() == other.rows() && cols() == other.cols()), function lazyAssign, file /usr/local/include/Eigen/src/Core/Assign.h, line 498.

Works for me:

Code: Select all

#include <Eigen/Dense>
using namespace Eigen;

void get_QR(Ref<const MatrixXd> A, Ref<MatrixXd> Q, Ref<MatrixXd> R) {
    HouseholderQR<MatrixXd> qr(A);
    Q = qr.householderQ()*(MatrixXd::Identity(A.rows(),A.cols()));
    R = qr.matrixQR().block(0,0,A.cols(),A.cols()).triangularView<Upper>();
}

void get_QR(const double* A, int m, int n, double* Q, double* R) {
    Map<const MatrixXd> A_E(A, m, n);
    Map<MatrixXd> Q_E(Q, m, n);
    Map<MatrixXd> R_E(R, n, n);
    get_QR(A_E, Q_E, R_E);
}

int main() {
  MatrixXd A = MatrixXd::Random(10,5);
  MatrixXd Q(10,5);
  MatrixXd R(5,5);
  get_QR(A, Q, R);
  get_QR(A.data(), 10, 5, Q.data(), R.data());
}


Thanks for replying. However, from the main, if I call as follows:

Code: Select all

int main() {
  MatrixXd A = MatrixXd::Random(10,5);
  MatrixXd Q;
  MatrixXd R;
  get_QR(A, Q, R);
  get_QR(A.data(), 10, 5, Q.data(), R.data());
}

without explicitly specifying the size of the matrices $Q$ and $R$, it fails giving the error I mentioned.

oh, I see. right. Ref objects cannot be resized. If you need both features, then add one more overload with Matrix arguments which resize and then call the Ref<> version.

ggael wrote:oh, I see. right. Ref objects cannot be resized. If you need both features, then add one more overload with Matrix arguments which resize and then call the Ref<> version.

Thanks ggael. Could you be more explicit (ideally with an example) as to what you mean by

add one more overload with Matrix arguments which resize and then call the Ref<> version.

Also, if I were to pass the matrix to the function with prototype as follows:

Code: Select all

void get_QR(MatrixXd A, MatrixXd& Q, MatrixXd& R)

I do not need to explicitly pass the size of Q and R. The size is need only when the function prototype is

Code: Select all

void get_QR(Ref<const MatrixXd> A, Ref<MatrixXd> Q, Ref<MatrixXd> R)

The reason why I am interested in dynamic sizes is that if I want to do an SVD with a desired tolerance, the size of my matrices, U, S and V are not known apriori. I would like the function to figure out the sizes and then allocate appropriate memory.

Thanks for your help.

then one version taking Matrix objects should be enough because you'll never be able to resize Map objects anyway.

ggael wrote:then one version taking Matrix objects should be enough because you'll never be able to resize Map objects anyway.

I do not understand your statement.