c++ - 新手问题 - 使用 Eigen 的复矩阵代数

Question

我是 C++ 新手，正在尝试将以下 MATLAB 代码转换为 C++:

x = [1,2,3];
y = [4,5,7];
px = 0.5;
py = 0.5;

m = -1:1:1;

v = [1+i,2+i,3+i;1+2i,2+2i,3+2i;1+3i,2+3i,3+3i];

Tr5 = real(exp(2i*pi/px*x(1)*m)*v*exp(2i*pi/py*m'*y(1)));

fprintf('Tr5 = %f\n',Tr5);

返回 Tr5=18。

这是我在 C++ 中的尝试，我在每一行都尽可能少地尝试查找错误:

#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>
#include <math.h>
#include <iterator>
#include <vector>
#include <algorithm>
#include <complex>
#include <iostream>

int main()
{

double x [3] = {1, 2, 3};
double y [3] = {4, 5, 7};
const std::complex<double> im(0.0, 1.0);
double px = 0.5;
double py = 0.5;

Eigen::MatrixXcd Tr1;
Eigen::MatrixXcd Tr2;
Eigen::MatrixXcd Tr3;
Eigen::MatrixXcd Tr4;

Eigen::MatrixXd m = Eigen::VectorXd::LinSpaced(2*1+1, -1, 1);
Eigen::MatrixXcd v4(3,3);
v4(0,0) = 1.0d + 1.0 * im;
v4(0,1) = 2.0d + 1.0 * im;
v4(0,2) = 3.0d + 1.0 * im;
v4(1,0) = 1.0d + 2.0 * im;
v4(1,1) = 2.0d + 2.0 * im;
v4(1,2) = 3.0d + 2.0 * im;
v4(2,0) = 1.0d + 3.0 * im;
v4(2,1) = 2.0d + 3.0 * im;
v4(2,2) = 3.0d + 3.0 * im;

Tr1 = 2.0*im*M_PI/px*x[1]*m;
Tr1 = Tr1.exp();

Tr2 = 2.0*im*M_PI/py*y[1]*m;
Tr2 = Tr2.exp();

Tr3 = Tr2.transpose();
Tr4 = Tr1 * v4 * Tr3;

Eigen::MatrixXcd Tr5 = Tr4.real(); // This should now be "18"

std::cout << "Tr5 = " << Tr5 << std::endl;

}

但是，当我尝试编译和运行时，出现错误:

/usr/local/include/unsupported/Eigen/src/MatrixFunctions/MatrixExponential.h:436: const Eigen::MatrixExponentialReturnValue<Derived> Eigen::MatrixBase<Derived>::exp() const [with Derived = Eigen::Matrix<std::complex<double>, -1, -1>]: Assertion `rows() == cols()' failed.

有人可以帮我处理代码吗？

非常感谢， 汤姆等待

Answer 1

您似乎没有使用系数指数，而是矩阵指数(包含在 中)矩阵指数(参见 https://en.wikipedia.org/wiki/Matrix_exponential)是方矩阵上的矩阵函数。您会看到该断言错误，因为您的矩阵“m”不是正方形。

您应该使用系数指数。在 Eigen 中，系数函数是为数组定义的，但它们也可用于矩阵，首先将其转换为数组:m.array()(参见 https://eigen.tuxfamily.org/dox/group__CoeffwiseMathFunctions.html)但是，您需要在执行矩阵乘法之前将数组转换回矩阵:leftSide.array().exp().matrix() * v * rightSide.array().exp().matrix()

最后，您可以使用逗号初始化语法(我相信这样更容易阅读)来填充您的矩阵 v。我将展示它作为示例，当然，您的初始化也很好。

示例(语法类似于 MATLAB):

#define _USE_MATH_DEFINES

#include "Eigen/Core"
#include <iostream>
#include <math.h>

std::complex<double> operator"" _i(long double x ) 
{ 
    return std::complex<double>(0.0, x);
} 

int main()
{
    double x [3] = {1, 2, 3};
    double y [3] = {4, 5, 7};

    double px = 0.5;
    double py = 0.5;

    Eigen::RowVectorXd m = Eigen::RowVectorXd::LinSpaced(3, -1, 1);

    Eigen::MatrixXcd v(3,3);
    v << (1.0 + 1.0_i), (2.0 + 1.0_i), (3.0 + 1.0_i),
         (1.0 + 2.0_i), (2.0 + 2.0_i), (3.0 + 2.0_i),
         (1.0 + 3.0_i), (2.0 + 3.0_i), (3.0 + 3.0_i);

    double Tr5 = ((2.0_i * M_PI / px * x[0] * m).array().exp().matrix() * v * (2.0_i * M_PI / py * y[0] * m).array().exp().matrix().transpose()).real()(0);

    std::cout << "Tr5 = " << Tr5 << std::endl;

    return 0;
}

这将产生您想要的输出:

Tr5 = 18