c++ - 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相同。但是在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.

关于问题的 Eigen 部分，避免使用非常大(超过几 KiB)的固定大小的矩阵。您可以有一个维度是固定的，而另一个维度是动态的。最后，避免将 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;