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

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

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 结合使用，除非你知道自己在做什么！

以下应该有效。我用相应的 Eigen 功能替换了所有循环，并通过直接使用对角矩阵进行乘积来避免临时循环。或者，您可以使用 replicate()d 矩阵执行 cwiseProduct。

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;

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

Python Numpy methods correspond to C++ Eigen make crash

c++

numpy

python-3.x

eigen