Hi, I needed to make a small computation and tried Eigen for the first time, but got a wrong result for various solvers (see the code below). In Octave I got the correct solution for the same input matrices, as shown below.
Can you tell me how can I make Eigen come up with the correct solution? I guess it has something to do with the computational precision or something that needs to be somehow set in Eigen. Thanks in advance.

Code: Select all

#include <Eigen/Dense>
typedef std::complex<double> Complex;
...
MatrixXcd r(6, 6);
    r << Complex(0.16,15.702),  Complex(-0.16,-15.702), Complex(0,0), Complex(0,15.34), Complex(0,-15.34), Complex(0,0),
         Complex(0,0), Complex(0.16,15.702), Complex(-0.16,-15.702), Complex(0,0), Complex(0,15.34), Complex(0,-15.34),
         Complex(-0.16,-15.702), Complex(0,0), Complex(0.16,15.702), Complex(0,-15.34), Complex(0,0), Complex(0,15.34),
         Complex(0,15.34), Complex(8.68828,0), Complex(-8.68828,0), Complex(0.078,15.853), Complex(8.97883,0), Complex(-8.97883,0),
         Complex(-8.68828,0), Complex(0,15.34), Complex(8.68828,0), Complex(-8.97883,0), Complex(0.078,15.853), Complex(8.97883,0),
         Complex(8.68828,0), Complex(-8.68828,0), Complex(0,15.34), Complex(8.97883,0), Complex(-8.97883,0), Complex(0.078,15.853);

VectorXcd u(6);

    u << Complex(330,190.526),
         Complex(0,-381.051),
         Complex(-330,190.526),
         Complex(0,0),
         Complex(0,0),
         Complex(0,0);

    //Eigen::PartialPivLU<MatrixXcd> t(r);
    //Eigen::FullPivLU<MatrixXcd> t(r);
    Eigen::FullPivHouseholderQR<MatrixXcd> t(r);
    //Eigen::ColPivHouseholderQR<MatrixXcd> t(r);
    //Eigen::LLT<MatrixXcd> t(r);
    //Eigen::LDLT<MatrixXcd> t(r);

    VectorXcd i(6);
    i = t.solve(u);
    std::cout << i << std::endl << std::endl << i.array().abs() << std::endl << std::endl;

    double relative_error = (r*i - u).norm() / u.norm(); // norm() is L2 norm
    cout << "The relative error is:\n" << relative_error << endl;
  ...


Output:
(64.1162,-71.0894)
(-27.043,-78.3232)
(12.2721,4.23987)
(-64.7021,56.1353)
(16.2636,84.1014)
(0,0)

95.7319
82.8604
12.9838
85.6594
85.6595
0

The relative error is:
8.74773e-07

Here is the Octave code:

Code: Select all

Complex = @(x,y) (x + y*i);

r = [Complex(0.16,15.702), Complex(-0.16,-15.702), Complex(0,0), Complex(0,15.34), Complex(0,-15.34), Complex(0,0);
     Complex(0,0), Complex(0.16,15.702), Complex(-0.16,-15.702), Complex(0,0), Complex(0,15.34), Complex(0,-15.34);
     Complex(-0.16,-15.702), Complex(0,0), Complex(0.16,15.702), Complex(0,-15.34), Complex(0,0), Complex(0,15.34);
     Complex(0,15.34), Complex(8.68828,0), Complex(-8.68828,0), Complex(0.078,15.853), Complex(8.97883,0), Complex(-8.97883,0);
     Complex(-8.68828,0), Complex(0,15.34), Complex(8.68828,0), Complex(-8.97883,0), Complex(0.078,15.853), Complex(8.97883,0);
     Complex(8.68828,0), Complex(-8.68828,0), Complex(0,15.34), Complex(8.97883,0), Complex(-8.97883,0), Complex(0.078,15.853)];

u = [Complex(330,190.526);
     Complex(0,-381.051);
     Complex(-330,190.526);
     Complex(0,0);
     Complex(0,0);
     Complex(0,0)];

rank = rank(r)

i = r \ u

abs(i)

relative_error = norm(r*i - u) / norm(u)


Octave output (the solution below is definitely correct):
rank = 5
warning: matrix singular to machine precision, rcond = 3.70599e-19
i =

47.6678 - 22.6985i
-43.4914 - 29.9323i
-4.1764 + 52.6308i
-48.5559 + 9.3897i
32.4098 + 37.3558i
16.1461 - 46.7455i

ans =

52.796
52.796
52.796
49.455
49.456
49.455

relative_error = 8.7477e-07

Your matrix is singular and so your problem accepts a sub-space of solutions. As you can see, both the solution provided by Eigen and Octave have the same relative error. If you are looking for the unique least-norm solution, then you currently have to use the SVD solver.

ggael wrote:Your matrix is singular and so your problem accepts a sub-space of solutions. As you can see, both the solution provided by Eigen and Octave have the same relative error. If you are looking for the unique least-norm solution, then you currently have to use the SVD solver.

You are right. But I did try SVD too, and got an incorrect solution (e+11, etc.) with various preconditioners. So, does it mean that Eigen cannot offer a method that would solve my problem correctly, like Octave's "\" did? Are there any other "tweaks" for SVD other than those that I tried in the code below?

Code: Select all

    VectorXcd i(6);
    JacobiSVD<MatrixXcd, ColPivHouseholderQRPreconditioner> t(r, ComputeThinU | ComputeThinV);
    // JacobiSVD<MatrixXcd, FullPivHouseholderQRPreconditioner> t(r, ComputeFullU | ComputeFullV);
    // JacobiSVD<MatrixXcd, HouseholderQRPreconditioner> t(r, ComputeThinU | ComputeThinV);
    // JacobiSVD<MatrixXcd, NoQRPreconditioner> t(r, ComputeThinU | ComputeThinV);
    i = t.solve(u);

    std::cout << "The least-squares solution is:\n" << i << std::endl << i.array().abs() << std::endl;
    double relative_error = (r*i - u).norm() / u.norm(); // norm() is L2 norm
    cout << "The relative error is:\n" << relative_error << endl;
    std::cout << "reconstructin:\n" << r*i << std::endl << std::endl;


Output:
The least-squares solution is:
(-1.95492e+11,8.7412e+10)
(-1.95492e+11,8.7412e+10)
(-1.95492e+11,8.7412e+10)
(1.89578e+11,-8.36507e+10)
(1.89578e+11,-8.36507e+10)
(1.89578e+11,-8.36507e+10)
2.14145e+11
2.14145e+11
2.14145e+11
2.07213e+11
2.07213e+11
2.07213e+11
The relative error is:
4.35e-06
reconstructin:
(330,190.524)
(-0.000517249,-381.052)
(-330,190.527)
(0.00030756,0.000229239)
(0.000424862,-0.000717521)
(0.000412107,0.000659466)

oh, it seems this has been fixed a few months ago. Need to make a 3.2.2 release. In the meantime you can use:

https://bitbucket.org/eigen/eigen/get/3.2.tar.bz2