Re: [eigen] Unsupported Polynomial solver

[Oleg Shirokobrod <oleg.shirokobrod@xxxxxxxxx>]

Hi Gael,

You finally modified PolynomialSolver.h code for finding roots of polynomial with complex coefficients and wrote unit test. This enhancement is very important for signal processing application. This modification of PolynomialSolver.h and Companion.h are rather trivial and localized. So it is safe to modify existing code. Would it be possible to include them to the next release?

Best regards,

Oleg

On Wed, Nov 23, 2016 at 5:09 PM Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

Done:

https://bitbucket.org/eigen/eigen/commits/3842186ae759/
https://bitbucket.org/eigen/eigen/commits/eff1bfef71a2/
https://bitbucket.org/eigen/eigen/commits/df6c723906fc/

gael

On Wed, Nov 23, 2016 at 3:39 PM, Gael Guennebaud <gael.guennebaud@xxxxxxxxx> wrote:

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