Python Numpy方法对应于C ++ Eigen make crash

Question

我有一个NumPy算术需要用Eigen转换C ++。

# D is a 20001x13 matrix stacked from Dva and Dvb, then multiply by w_f.
# dtype=complex<double>
D = np.column_stack((Dva, Dvb)) * w_f.reshape((20001, 1)) * np.ones((1, 13))
R = np.dot(D.conj().T, D)

这是我的C ++代码（最小测试）：

#include <Eigen/Core>
#include <Eigen/Dense>
#include <vector>
#include <complex>

using namespace std;

typedef complex<double> dcomplex;

void foo()
{
    vector<dcomplex> wf;
    wf.resize(20001);
    Eigen::Matrix<dcomplex, 20001, 13> *tmp_1 = new Eigen::Matrix<dcomplex, 20001, 13>;
    Eigen::Matrix<dcomplex, 20001, 13> *tmp_2 = new Eigen::Matrix<dcomplex, 20001, 13>;
    Eigen::Matrix<dcomplex, 20001, 7> *Dva = new Eigen::Matrix<dcomplex, 20001, 7>;
    Eigen::Matrix<dcomplex, 20001, 6> *Dvb = new Eigen::Matrix<dcomplex, 20001, 6>;

    for (int i = 0; i < 20001; i++){
        for (int j = 0; j < 7; j++)
            (*Dva)(i, j) = 0;
        for (int j = 0; j < 6; j++)
            (*Dvb)(i, j) = 0;
        for (int j = 0; j < 13; j++)
            (*tmp_2)(i, j) = wf[i];
    }
    *tmp_1 << *Dva, *Dvb;
    auto *D = &tmp_1->cwiseProduct(*tmp_2);

    auto R = (D->transpose() * (*D));
    R(0,0);
}

矩阵R的形状在Eigen中是13x13，与NumPy相同。但是变量R不能用C ++表示。

R.rows() == 13;  // true
R.cols() == 13;  // true
R(0, 0);  // or what ever makes it crash

引发异常“0xC00000FD：堆栈溢出”。

Answer 1

首先，你应该在C ++代码中使用new。大多数时候使用本地对象（或std::vector），如果需要，使用std::unique_ptr或std::shared_ptr等智能指针。

关于问题的特征部分，避免使用非常大（超过几KiB）的固定大小的矩阵。你可以修复一个维度，另一个Dynamic。最后，避免将auto与Eigen结合使用，除非你知道自己在做什么！

以下应该有效。我用相应的特征函数替换了所有的循环，并通过直接用对角矩阵做一个产品来避免临时性。或者，您可以使用cwiseProduct矩阵进行replicate()d。

typedef Eigen::Matrix<dcomplex, Eigen::Dynamic,1> VectorXcd;
typedef Eigen::Matrix<dcomplex, Eigen::Dynamic,13> MatrixX13cd;
typedef Eigen::Matrix<dcomplex, Eigen::Dynamic,7> MatrixX7cd;
typedef Eigen::Matrix<dcomplex, Eigen::Dynamic,6> MatrixX6cd;
typedef Eigen::Matrix<dcomplex, 13,13> Matrix13cd;


MatrixX7cd  Dva(20001,  7);
MatrixX6cd  Dvb(20001,  6);

Dva.setZero(); Dvb.setZero();
MatrixX13cd D(20001, 13);
D.leftCols(7).noalias()  = VectorXcd::Map(wf.data(), wf.size()).asDiagonal() * Dva;
D.rightCols(6).noalias() = VectorXcd::Map(wf.data(), wf.size()).asDiagonal() * Dvb;

Matrix13cd R = D.transpose() * D;