Re: [eigen] solve_complex_linear_system
• To: eigen@xxxxxxxxxxxxxxxxxxx
• Subject: Re: [eigen] solve_complex_linear_system
• From: Bill Greene <w.h.greene@xxxxxxxxx>
• Date: Fri, 28 Apr 2017 11:10:55 -0400
• Dkim-signature: v=1; a=rsa-sha256; c=relaxed/relaxed; d=gmail.com; s=20161025; h=mime-version:in-reply-to:references:from:date:message-id:subject:to; bh=FIOsi3CMoNvO9kwphNkyzG0YqXN8Er1edryj86rUogc=; b=cOYQ8M9OxZT6pryeYKMkbRvYLdGPoYQC9dYwheNquZ7hZAZIkcLCJQSCdbSfUQf0DP XiPa8zTq2RT1509iLuECuCpydgFhL6+VpqzPGIPp/ZpPHZtycwgCI16tx3oTiQkHBRjh sZivUJx01jZ2R462CWC+Zqveazu2YU70DNxv1yY34BUi4EnNAtf0r3H1C2ztOBVIrjdb djrE2gey5yCSKgoL0Z8heKsY6/T2n948sH9CeyjrbYVe1z50Lv1z4WkfAQlB+VX0H4T8 Rr08jy5sJDk8eVePLgHqeJRZLkHhMfhGRImpgeL9+JDADBK+FwGrLdv0gDsFs3t6dapF 8Rtw==

```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

```

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