使用特征时的奇怪行为

Weird Behavior when using Eigen

本文关键字:特征      更新时间:2023-10-16

我正在为Eigen编写一个包装器供我个人使用,我遇到了以下奇怪的行为:

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>();
}
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();
}

以下是测试:

int main(int argc, char* argv[]) {
    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;
    double* R;
    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;
}

我得到不同的 Q[0] 值。例如,我得到"-0.421857"和"-1.49563"。

谢谢

乔治的答案是正确的,但存在不必要的副本。更好的解决方案包括映射 Q 和 R:

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(Q, n, n);
    HouseholderQR<MatrixXd> qr(A_E);
    Q_E = qr.householderQ()*(MatrixXd::Identity(m,n));
    R_E = qr.matrixQR().block(0,0,n,n).triangularView<Upper>();
}

为了能够重用采用 Eigen 对象的 get_QR 函数,请使用 Ref<MatrixXd> 而不是 MatrixXd

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);
}

Ref<MatrixXd>可以包装任何类似于MatrixXd的特征对象,而无需任何副本。这包括MatrixXd本身,以及MapBlock表达式。这需要特征 3.2。

我认为这与Eigen没有任何关系。

看起来您正在将指针 Q 分配给属于局部变量的内存位置Q_E

        Q       =       Q_E.data();

泄漏以前的内存分配

double* Q       =       new double[m*n];

并且在get_QR()函数之外毫无意义或未定义。

你应该改用memcpy:

memcpy(Q, Q_E.data(), m*n*sizeof(double));