Invalid Results Using Linear System Solvers Eigen 3
I am porting my MATLAB code to C ++ from Eigen 3 and I decided to build my linear solver rather than call it from a matrix object so I can reuse it. Unfortunately, it does not give the expected results. After several tests, I traced the issue with what looks like a linear solver object, as evidenced by the following relevant snippet of my code:
MatrixXcd M(6,6), Y(6,6), dvalor(6,6);
// Initialization of the matrices...
FullPivHouseholderQR<MatrixXcd> mdivide;
auto result = mldivide.compute(Y).solve(M).eval() *
dvalor * mldivide.compute(M).solve(Y).eval();
std::cout << "result:" << std::endl
<< result << std::endl << std::endl;
std::cout << "using fullPivHouseholderQr:" << std::endl
<< Y.fullPivHouseholderQr().solve(M) * dvalor *
M.fullPivHouseholderQr().solve(Y) << std::endl << std::endl;
std::cout << "using jacobi:" << std::endl
<< Y.jacobiSvd(ComputeThinU | ComputeThinV).solve(M) * dvalor *
M.jacobiSvd(ComputeThinU | ComputeThinV).solve(Y) << std::endl << std::endl;
std::cout << "using inverse:" << std::endl
<< Y.inverse() * M * dvalor * M.inverse() * Y << std::endl << std::endl;
This is the result:
result:
(0.564196,0.00606298) (-0.15282,-0.00179293) (-0.16564,-0.000220726) (-0.16564,-0.000220726) (-0.15282,-0.00179293) (-0.179235,0.00541725)
(2.73563e+184,-8.89017e+185) (1,1.12229e-10) (-6.04297e-11,5.62259e-12) (-6.02256e-11,2.60221e-12) (-5.33448e-11,-3.05694e-11) (-6.55427e-11,8.17723e-12)
(2.11057e+184,-8.81209e+185) (-5.53615e-11,4.54731e-12) (1,1.11384e-10) (-5.99088e-11,9.29495e-12) (-5.51744e-11,1.52693e-12) (-6.21597e-11,-3.11176e-11)
(2.11057e+184,-8.81209e+185) (-0.171319,-0.00113681) (0.134093,-2.05541e-05) (0.865907,2.05542e-05) (0.171319,0.00113681) (-0.140389,-0.000216684)
(2.73563e+184,-8.89017e+185) (0.218904,0.00143522) (-0.171337,3.98397e-05) (0.171337,-3.98397e-05) (0.781096,-0.00143522) (0.179383,0.000262654)
(1.4668e+184,-8.78677e+185) (-5.52732e-11,5.99039e-12) (-5.69133e-11,-3.2223e-11) (-5.97328e-11,9.63533e-12) (-5.51767e-11,4.01574e-12) (1,1.1063e-10)
using fullPivHouseholderQr:
(1,1.1063e-10) (-5.51766e-11,4.01578e-12) (-5.97331e-11,9.63535e-12) (-5.6913e-11,-3.22229e-11) (-5.52733e-11,5.99029e-12) (-6.46967e-11,1.0973e-11)
(-6.54343e-11,6.20279e-12) (1,1.12229e-10) (-6.04297e-11,5.62259e-12) (-6.02256e-11,2.60221e-12) (-5.33448e-11,-3.05694e-11) (-6.55427e-11,8.17723e-12)
(-6.48883e-11,1.07408e-11) (-5.53615e-11,4.54731e-12) (1,1.11384e-10) (-5.99088e-11,9.29495e-12) (-5.51744e-11,1.52693e-12) (-6.21597e-11,-3.11176e-11)
(-6.21592e-11,-3.11175e-11) (-5.51743e-11,1.52697e-12) (-5.99093e-11,9.29482e-12) (1,1.11385e-10) (-5.53614e-11,4.54728e-12) (-6.48881e-11,1.07408e-11)
(-6.55427e-11,8.1773e-12) (-5.33444e-11,-3.05693e-11) (-6.02261e-11,2.60223e-12) (-6.04293e-11,5.62257e-12) (1,1.12229e-10) (-6.54344e-11,6.20272e-12)
(-6.46967e-11,1.0973e-11) (-5.52732e-11,5.99039e-12) (-5.69133e-11,-3.2223e-11) (-5.97328e-11,9.63533e-12) (-5.51767e-11,4.01574e-12) (1,1.1063e-10)
using jacobi:
(1,1.1063e-10) (-5.51761e-11,4.01562e-12) (-5.9733e-11,9.63488e-12) (-5.69134e-11,-3.22232e-11) (-5.5273e-11,5.9903e-12) (-6.46969e-11,1.09735e-11)
(-6.54346e-11,6.20257e-12) (1,1.12228e-10) (-6.04303e-11,5.62232e-12) (-6.0225e-11,2.60224e-12) (-5.33448e-11,-3.05694e-11) (-6.55414e-11,8.17752e-12)
(-6.48874e-11,1.07408e-11) (-5.53617e-11,4.54725e-12) (1,1.11384e-10) (-5.99097e-11,9.29465e-12) (-5.51742e-11,1.52688e-12) (-6.21596e-11,-3.11175e-11)
(-6.21608e-11,-3.11175e-11) (-5.5174e-11,1.52688e-12) (-5.99104e-11,9.29433e-12) (1,1.11384e-10) (-5.53612e-11,4.54728e-12) (-6.48873e-11,1.07413e-11)
(-6.55423e-11,8.17718e-12) (-5.33459e-11,-3.05695e-11) (-6.02258e-11,2.60206e-12) (-6.04294e-11,5.6224e-12) (1,1.12229e-10) (-6.5434e-11,6.20302e-12)
(-6.46966e-11,1.09729e-11) (-5.5273e-11,5.99037e-12) (-5.69134e-11,-3.22233e-11) (-5.97333e-11,9.63543e-12) (-5.51767e-11,4.01582e-12) (1,1.1063e-10)
using inverse:
(1,1.1063e-10) (-5.51768e-11,4.01574e-12) (-5.97329e-11,9.63535e-12) (-5.6914e-11,-3.2223e-11) (-5.52733e-11,5.99031e-12) (-6.46967e-11,1.0973e-11)
(-6.54338e-11,6.20274e-12) (1,1.12229e-10) (-6.04293e-11,5.6225e-12) (-6.02264e-11,2.6022e-12) (-5.33451e-11,-3.05693e-11) (-6.55425e-11,8.1773e-12)
(-6.48877e-11,1.07408e-11) (-5.53615e-11,4.54726e-12) (1,1.11384e-10) (-5.99097e-11,9.29488e-12) (-5.51743e-11,1.52697e-12) (-6.21592e-11,-3.11176e-11)
(-6.21587e-11,-3.11176e-11) (-5.51744e-11,1.52696e-12) (-5.99093e-11,9.29489e-12) (1,1.11384e-10) (-5.53616e-11,4.54727e-12) (-6.4888e-11,1.07408e-11)
(-6.55419e-11,8.17731e-12) (-5.33448e-11,-3.05693e-11) (-6.02257e-11,2.6022e-12) (-6.04301e-11,5.6225e-12) (1,1.12229e-10) (-6.54343e-11,6.20273e-12)
(-6.46964e-11,1.0973e-11) (-5.52732e-11,5.99031e-12) (-5.69131e-11,-3.2223e-11) (-5.97336e-11,9.63535e-12) (-5.51767e-11,4.01574e-12) (1,1.1063e-10)
The first iteration of the algorithm in MATLAB creates an almost identical matrix, which can be seen in C ++ by accessing matrices and calling the solver from the matrix object.
But the result is completely wrong when using the solver object. I have also tried other solvers from my own documentation page with the same result: https://eigen.tuxfamily.org/dox/group__TutorialLinearAlgebra.html .
Am I missing any step or will it call the solver from the matrix, else will it initialize the solver?
The native version is 3.3.4 and I am not using compiler optimizers yet, just C ++ 11:
g++ -I. -std=c++11 main.cpp
The size of the matrices is constant and it will run inside a large loop. Thus, the more memory allocation I can avoid by pulling out the solver, the better.
Thank.
Be careful with auto
, see common pitfalls , mainly in your case result
is an expression Product<>
referencing dead objects, Replace auto
with MatrixXcd
and solve the problem.
In the first example, you are using mdivide
twice in the same expression. .compute()
changes the internal state.
Using two objects mdivide
, one for Y
and one for, M
may give you better results.