使用 Eigen 的复杂矩阵矩阵乘法

Question

我正在做复杂的矩阵矩阵乘积，如这里的代码所示。 Eigen * 运算符似乎没有给出正确的结果。为了验证我还编写了另一个例程来查找产品。由此可见，结果的虚部是好的，但实部是不准确的。谁能告诉我为什么？

Eigen::MatrixXcd U2 = Eigen::MatrixXcd::Random(N1, computed_rank);
Eigen::MatrixXcd V2 = Eigen::MatrixXcd::Random(computed_rank, N2);
Eigen::MatrixXcd W2 = Mat::Zero(N1,N2);
for (size_t k = 0; k < computed_rank; k++) {
    W2 += U2.col(k)*V2.row(k);
}
Eigen::MatrixXcd Q = U2*V2;
Eigen::MatrixXcd E2 = Q-W2;
cout << "Err Total: " << E2.norm()/W2.norm() << endl;
cout << "Err Real: " << E2.real().norm()/W2.real().norm() << endl;
cout << "Err Imag: " << E2.imag().norm()/W2.imag().norm() << endl;

给出了结果：

Err Total: 0.84969
Err Real: 1.17859
Err Imag: 1.18274e-16

我观察到当这是项目中唯一的一段代码时并没有出现这个问题。但我把它作为一个更大项目的一部分，它似乎失败了。

Answer 1

事实证明，GNU 的-ffast-math 优化标志并不兼容所有程序。我在编译更大的项目时使用过它，而不是在做一个孤立的项目时使用过它。这就是它失败的原因！有关详细信息，请参阅 https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html.

Answer 2

下面是一个最小的可重现示例：文件名是 test.cpp

#include<iostream>
#include<Eigen/Dense>

using namespace std;

int main() {
    int N1 = 12;
    int N2 = 12;
    int computed_rank = 5;
    Eigen::MatrixXcd U2 = Eigen::MatrixXcd::Random(N1, computed_rank);
    Eigen::MatrixXcd V2 = Eigen::MatrixXcd::Random(computed_rank, N2);
    Eigen::MatrixXcd W2 = Eigen::MatrixXcd::Zero(N1,N2);
    for (size_t k = 0; k < computed_rank; k++) {
        W2 += U2.col(k)*V2.row(k);
    }
    Eigen::MatrixXcd Q = U2*V2;
    Eigen::MatrixXcd E2 = Q-W2;
    cout << "Err Total: " << E2.norm()/W2.norm() << endl;
    cout << "Err Real: " << E2.real().norm()/W2.real().norm() << endl;
    cout << "Err Imag: " << E2.imag().norm()/W2.imag().norm() << endl;
}

使用 -ffast-math 标志编译上述 会产生以下输出。

g++-9 -O4 -ffast-math test.cpp

./a.out

Err Total: 0.949251 Err Real: 1.41443 Err Imag: 1.35749e-16

在没有 -ffast-math 标志的情况下编译上述 会产生以下输出。

g++-9 -O4 test.cpp

./a.out

Err Total: 1.35029e-16 Err Real: 1.35181e-16 Err Imag: 1.34903e-16

使用 Eigen 的复杂矩阵矩阵乘法

Complex Matrix matrix multiplication using Eigen

c++

matrix

complex-numbers

matrix-multiplication

eigen