[eigen] solve_complex_linear_system
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: [eigen] solve_complex_linear_system
• From: Moritz Dannhauer <moritz.dannhauer@xxxxxxxxx>
• Date: Thu, 27 Apr 2017 15:43:24 -0400
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20161025; h=mime-version:from:date:message-id:subject:to; bh=w9ocLnHoc2EfYQLWcB2IgXHF2TX6QGx21r3/rjmYQTk=; b=trNJxlKBH7vUOn6KznqRevcEZS0a3b5icIXlO0rEM59eu80sEWmfOZm8r491urSXVx ITS9pSw5jS876IMvPgvFeu8TmqpJucCJ8eemvffBm+RIOm5bq05s7zrGpI8Tp5j2zkQI bxQjPOILQfVktQ5t/PrWY7MiLmkaP6IrPUDhY23Owl2WmGTUkLfjTA18pvaJcb6EzGpo 9yPDIUV72DvyAMCqyL/HsOQrIyzuEFeJMqj1XflOlCHBxGvHYN1DRDbV+fVARIS1Vzj3 brzy7pmqsyykgOdjdLUKf+VfF/uwpgES+bPNrv6MkGYOrPi1Lt0Lg6UoWtBKlqIcabeY P0Nw==

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_;
};

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 in that regard?

Thanks so much!
Best,
Moritz

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/