Re: [eigen] solve_complex_linear_system

[Bill Greene <w.h.greene@xxxxxxxxx>]

It is unclear to me what your right-hand-side vector is.
You show:
b=[1 2+i 3];
But your matlab solution is correct only if what you mean is
b=[1;2-i;3];
i.e. you are taking the conjugate transpose of your b row vector.

If I take b to be
b=[1; 2+i; 3];

I get the solution to be

   1.18182 - 0.27273i
   0.54545 + 0.81818i
   1.90909 - 0.27273i

from Octave, Eigen::ConjugateGradient, and Eigen::BiCGSTAB. I am using
Eigen 3.3.3.

On Thu, Apr 27, 2017 at 3:43 PM, Moritz Dannhauer
<moritz.dannhauer@xxxxxxxxx> wrote:
> 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