Complex Matrix matrix multiplication using Eigen

Question

I am doing complex matrix matrix product as shown in the code here. The Eigen * operator doesn't seem to give the right result. To verify I have also written another routine to find the product. From this it was seen that the imaginary part of the result was fine but the real part is inaccurate. Can someone let me know why?

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;

which gave the result:

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

I have observed that this problem didn't occur when this was the only piece of code in the project. But I have this as part of a bigger project, where it seems to fail.

How large is `W.real().norm()`? It's not hard to construct matrices `U2` and `V2` such that the real part of the product (or even the entire product) gets very small, thus the relative error can get huge. This may just have happened by chance in your case. Also, you did not really show whether `W2` or `Q` is "more accurate". — chtz, Dec 17 '19 at 07:04
Also, just of curiosity, how large are `N1`, `N2`, `computed_rank`? Are you able to produce a [mre] (with fast-math en-/disabled)? — chtz, Dec 17 '19 at 07:08
@chtz Please see my answer below for a minimum reproducible example. — Adhvaitha, Dec 17 '19 at 08:09

score 2 · Accepted Answer · answered Dec 16 '19 at 14:22

2

It turns out that -ffast-math optimization flag of GNU is not compatible with all programs. I have used it while compiling the bigger project and not while doing the one that is in isolation. This is the reason it failed! For more details refer https://gcc.gnu.org/onlinedocs/gcc/Optimize-Options.html.

answered Dec 16 '19 at 14:22

Vaishnavi Gujjula ma16d301

41
3

You may suggest this in the Eigen forums: https://forum.kde.org/viewforum.php?f=74 – Adhvaitha Dec 16 '19 at 18:30
1

A while ago, when I had a similar problem, I found [this answer](https://stackoverflow.com/a/38980445/4770166) on the topic particularly helpful. – RHertel Dec 17 '19 at 08:12

score 1 · Answer 2 · answered Dec 17 '19 at 08:08

Below is a minimum reproducible example: Filename is 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;
}

Compiling the above with -ffast-math flag produces the following output.

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

./a.out

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

Compiling the above without -ffast-math flag produces the following output.

g++-9 -O4 test.cpp

./a.out

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

Complex Matrix matrix multiplication using Eigen

2 Answers2