Re: [eigen] Unsupported Polynomial solver

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Sorry for late reply, but thanks for the patches.

Regarding the eigen-solver choice issue, it's possible to automatically use the right one based on NumTraits<Scalar>::IsComplex and conditional<>. 

We would also need to extended the unit test to cover complexes.

gael

On Tue, Nov 8, 2016 at 3:41 PM, Oleg Shirokobrod <oleg.shirokobrod@xxxxxxxxx> wrote:

Dear list,

Polynomial solver currently finds roots only for real polynomial coefficients. Having Eigen::ComplexEigenSolver lets Polynomial module to compute roots for complex polynomial coefficients. I have made necessary modifications. I have attached corresponding patch files. Unfortunately module cannot deduce type of the solver from partial specialization of corresponding matrix type, such that

template<typename RealScalar> class EigenSolver {};

template<typenameRealScalar> class EigenSolver<Eigen::Matrix<RealScalar,_Deg,_Deg> > {// EigenSolver };

template<typename  RealScalar> class EigenSolver<Eigen::Matrix<std::complex<RealScalar>,_Deg,_Deg > > { //ComplexEigenSolver};

So I have to add additional template parameter EigenSolverType with default value for real coefficients:

template< typename _Scalar, int _Deg , typename EigenSolverType = EigenSolver<Matrix<_Scalar,_Deg,_Deg> > >
class PolynomialSolver : public PolynomialSolverBase<_Scalar,_Deg>{};

I have to replace in number of places in the file Companion.h (mainly in functions balanced() and balance()) Scalar with RealScalar, where variables are really real.

With this code I have run test against Matlab and it gave similar results.

Test
VectorXcd roots = VectorXcd::Random(4);
VectorXcd polynomialCoefficients;
roots_to_monicPolynomial(roots, polynomialCoefficients);

cout << "roots : " << endl;
cout << setprecision(14) << roots << endl;
cout << "polynomialCoefficients : " << endl;
cout << setprecision(14)  << polynomialCoefficients << endl;

PolynomialSolver<std::complex<double>, Dynamic, ComplexEigenSolver<MatrixXcd> > psolve(polynomialCoefficients);

VectorXcd computedRoots = psolve.roots();

cout << "computedRoots : " << endl;
cout << setprecision(14) << computedRoots << endl;

for(int index = 0; index < computedRoots.size(); ++index)
{
    cout << setprecision(14) << "polynom at  computedRoots * 1.0 e-11: " << poly_eval(polynomialCoefficients, computedRoots(index))*1.0e11 << endl;
}

Output:
roots :
(0.12717062898648,-0.99749748222297)
(0.61748100222785,-0.61339152195807)
(-0.04025391399884,0.17001861629078)
(0.79192480239265,-0.29941709646901)

polynomialCoefficients :
(0.091649023983597,-0.091441416775918)
(0.24020598856634,0.37401934653925)
(-0.16301627948124,-1.8544616197629)
(-1.4963225196081,1.7402874843593)
(1,0)

computedRoots :
(-0.04025391399884,0.17001861629078)
(0.79192480239265,-0.29941709646901)
(0.61748100222785,-0.61339152195807)
(0.12717062898648,-0.99749748222297)

polynom at  computedRoots * 1.0 e-11: (8.3266726846887e-006,-1.2490009027033e-005)
polynom at  computedRoots * 1.0 e-11: (4.7184478546569e-005,-1.6653345369377e-005)
polynom at  computedRoots * 1.0 e-11: (2.2204460492503e-005,-1.3877787807814e-005)
polynom at  computedRoots * 1.0 e-11: (1.5163809063462e-005,-2.7286889655471e-005)

Matlab for the same coefficients gives following roots (remember that matlab array of coefficients and Eigen one are reversed to each other)

computedRoots :
0.127170628986480 - 0.997497482222969i
0.617481002227849 - 0.613391521958071i
0.791924802392650 - 0.299417096469009i
-0.040253913998840 + 0.170018616290780i

Best regards,

Oleg Shirokobrod