Re: [eigen] Unsupported Polynomial solver

[ Thread Index | Date Index | More lists.tuxfamily.org/eigen Archives ]


Hi,

the modification happened in the default (i.e., "development") branch, after the initial 3.3 release. Usually, we don't backport added features, but only bug-fixes to the stable branch. However, there is no strict rule about this, especially if this is not too complicated.

So no objections from me to backport this.

Christoph


On 21/03/2019 07.20, Oleg Shirokobrod wrote:
Hi Gael,

I mean that this modification was made in some branch but is not included
in any release 3.3.0 - 3.3.7. Would it be possible to make this
modification in the next 3.3.8 release?

On Thu, Mar 21, 2019 at 12:35 AM Gael Guennebaud <gael.guennebaud@xxxxxxxxx>
wrote:

Hi Oleg,

do you mean a backport to 3.3?

gael

On Wed, Mar 20, 2019 at 3:51 PM Oleg Shirokobrod <
oleg.shirokobrod@xxxxxxxxx> wrote:

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







--
 Dr.-Ing. Christoph Hertzberg

 Besuchsadresse der Nebengeschäftsstelle:
 DFKI GmbH
 Robotics Innovation Center
 Robert-Hooke-Straße 5
 28359 Bremen, Germany

 Postadresse der Hauptgeschäftsstelle Standort Bremen:
 DFKI GmbH
 Robotics Innovation Center
 Robert-Hooke-Straße 1
 28359 Bremen, Germany

 Tel.:     +49 421 178 45-4021
 Zentrale: +49 421 178 45-0
 E-Mail:   christoph.hertzberg@xxxxxxx

 Weitere Informationen: http://www.dfki.de/robotik
  -------------------------------------------------------------
  Deutsches Forschungszentrum für Künstliche Intelligenz GmbH
  Trippstadter Strasse 122, D-67663 Kaiserslautern, Germany

  Geschäftsführung:
  Prof. Dr. Jana Koehler (Vorsitzende)
  Dr. Walter Olthoff

  Vorsitzender des Aufsichtsrats:
  Prof. Dr. h.c. Hans A. Aukes
  Amtsgericht Kaiserslautern, HRB 2313
  -------------------------------------------------------------




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