Hi Eigen experts,
I am part of a developer team to evaluate Eigen's capabilities to solve a system of complex linear equation, such as A * x = b.
Our test case are small self-adjoint matrices such as (here as matlab code):
A=[-1 1-2i 0; 1+2i 1 -i; 0 i 2];
b=[1 2+i 3];
Eigen's solution is:
x_eigen=[0.9245 - 0.7358i 0.6792 + 0.6226i 0.7358 - 1.0755i];
Matlab's solution is correct though:
x_matlab=[0.4545 + 0.6364i 0.5455 - 0.4545i 1.7273 + 0.2727i];

We are using Eigen's Conjugate Gradient or BiConjugate Gradient algorithm to solve this, e.g.:

template <class ColumnMatrixType, template <typename> class SolverType>
class SolveLinearSystemAlgorithmEigenCGImpl
{
public:
SolveLinearSystemAlgorithmEigenCGImpl(SharedPointer<ColumnMatrixType> rhs, double tolerance, int maxIterations) :
tolerance_(tolerance), maxIterations_(maxIterations), rhs_(rhs) {}

using SolutionType = ColumnMatrixType;

template <class MatrixType>
typename ColumnMatrixType::EigenBase solveWithEigen(const MatrixType& lhs)
{
SolverType<typename MatrixType::EigenBase> solver;
solver.compute(lhs);

if (solver.info() != Eigen::Success)
BOOST_THROW_EXCEPTION(AlgorithmInputException()
<< LinearAlgebraErrorMessage("Eigen solver initialization was unsuccessful")
<< EigenComputationInfo(solver.info()));

solver.setTolerance(tolerance_);
solver.setMaxIterations(maxIterations_);
auto solution = solver.solve(*rhs_).eval();
tolerance_ = solver.error();
maxIterations_ = solver.iterations();
return solution;
}

double tolerance_;
int maxIterations_;
private:
SharedPointer<ColumnMatrixType> rhs_;
};
----cut----

The algorithm runs (e.g., 5000 iterations and target error of 0.00001) quickly but outputs a wrong solution (see above).
We checked obvious things (like that the inputs and outputs gets messed up somehow) indicating that we did a mistake - but that seems not to be the case.
Has anyone experienced that as well?
Is there anything we need to do to make it work (e.g., specific algorithm parameters, specific matrix properties we missed to check for)?
Are there limits what Eigen can do?

Thanks so much!
Best,
Moritz

Sorry for late reply,

there seems to be a problem with the way you setup the Eigen's code, here is a self-contained example showing that both Eigen::LDLT and Eigen::ConjugateGradient properly solved the problem, whereas the solution you gave for MatLab does not solve the initial problem:

Code: Select all

#include <iostream>
#include <Eigen/Dense>
#include <Eigen/IterativeLinearSolvers>
using namespace Eigen;
using namespace std;

int main()
{
  typedef std::complex<double> c;
  MatrixXcd A(3,3);
  // A=[-1 1-2i 0; 1+2i 1 -i; 0 i 2];
  A << -1,  c(1,-2), 0,
       c(1,2), 1, c(0,-1),
       0, c(0,1), 2;
  VectorXcd b(3);
  b << 1, c(2,1), 3;

  VectorXcd x1(3), x2(3), x3(3), xm(3);
  x1 << c(0.9245,-0.7358), c(0.6792,0.6226), c(0.7358, -1.0755);
  xm << c(0.4545,0.6364), c(0.5455, -0.4545), c(1.7273,0.2727);

  x2 = A.ldlt().solve(b);
  ConjugateGradient<MatrixXcd> cg(A);
  x3 = cg.solve(b);

  cout << "x1: " << x1.transpose() << "   -> " << (A*x1-b).norm() / b.norm() << endl;
  cout << "x2: " << x2.transpose() << "   -> " << (A*x2-b).norm() / b.norm() << endl;
  cout << "x3: " << x3.transpose() << "   -> " << (A*x3-b).norm() / b.norm() << endl;
  cout << "xm: " << xm.transpose() << "   -> " << (A*xm-b).norm() / b.norm() << endl;
}