[eigen] Polynomial solver, eigenvalues of companion matrix and balancing

[Manuel Yguel <manuel.yguel@xxxxxxxxx>]

Hello,
I have written a class to compute the roots of a real polynomial.
I build the companion matrix and compute its eigenvalues.
You can see the code in the attached file (not a patch yet ... see the
following).
I have some problems with that method:

A] When testing: I encounter almost systematic failure for polynomials
with deg greater than 7 (see test file attached),
therefore I wonder:

1) Do I use the right eigensolver ?
2) Does the eigensolver balance the input matrix ?

I have written some code to balance a companion matrix explicitly (and
I am testing this stuff at the moment) but before investing more time
in that direction I need to know if doing the balancing is redundant.


B] A somehow related question is: the eigensolver make a copy of the
matrix however, the companion matrix is sparse and for high degree
polynomial this copy could be avoided.
I have thought about writing the companion matrix class as a
cwiseNullaryOp. Do I follow the right thread?


- Best regards,

Manuel

P.S. I am also writing a solver with a Bézier bissection, it is
claimed to perform faster.

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Manuel Yguel <manuel.yguel@xxxxxxxxx>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#ifndef EIGEN_POLYNOMIAL_SOLVER_H
#define EIGEN_POLYNOMIAL_SOLVER_H

// This file requires the user to include
// * Eigen/Core
// * Eigen/Eigenvalues

#ifndef EIGEN_PARSED_BY_DOXYGEN

template<int Size>
struct ei_decrement_if_fixed_size
{
  enum {
    ret = (Size == Dynamic) ? Dynamic : Size-1 };
};

#endif

/** \Polynomials_module \ingroup Polynomials_Module
  *
  * \class PolynomialSolver
  *
  * \brief A polynomial solver
  *
  * Computes the complex roots of a real polynomial.
  *
  * \param _Scalar the scalar type, i.e., the type of the polynomial coefficients
  * \param _Deg the degree of the polynomial, can be a compile time value or Dynamic.
  *             Notice that the number of polynomial coefficients is _Deg+1.
  *
  * This class implements a polynomial solver and provide convenient methods such as
  * - real roots,
  * - greatest, smallest complex roots,
  * - real roots with greatest, smallest absolute real value.
  * - greatest, smallest real roots.
  */
template< typename _Scalar, int _Deg >
class PolynomialSolver
{
  public:
    EIGEN_MAKE_ALIGNED_OPERATOR_NEW_IF_VECTORIZABLE_FIXED_SIZE(_Scalar,_Deg==Dynamic ? Dynamic : _Deg)

    typedef _Scalar                                  Scalar;
    typedef typename NumTraits<Scalar>::Real         RealScalar;
    typedef Matrix<Scalar,_Deg,_Deg>                 CompanionMatrixType;
    typedef EigenSolver<CompanionMatrixType>         EigenSolverType;
    typedef typename EigenSolverType::Complex        RootType;
    typedef typename EigenSolverType::EigenvalueType RootsType;

  public:
    /** Computes the complex roots of a new polynomial. */
    template< typename OtherPolynomial >
    void compute( const OtherPolynomial& poly )
    {
      typedef ei_decrement_if_fixed_size<_Deg>          Deg_1;

      typedef Matrix< Scalar, _Deg, Deg_1::ret >        LeftBlock;
      typedef Matrix< Scalar, Deg_1::ret, Deg_1::ret >  BottomLeftBlock;
      typedef Matrix< Scalar, 1, Deg_1::ret >           LeftBlockFirstRow;

      const int deg = poly.size()-1;
      const int deg_1 = deg-1;
      CompanionMatrixType companion(deg,deg);
      companion <<
        ( LeftBlock(deg,deg_1) << LeftBlockFirstRow::Zero(1,deg_1),
                                 BottomLeftBlock::Identity(deg_1,deg_1) ).finished(),
        -1.0/poly[deg]*poly.start(deg);
      m_eigenSolver.compute( companion );
      m_roots = m_eigenSolver.eigenvalues();
    }

  public:
    template< typename OtherPolynomial >
    inline PolynomialSolver( const OtherPolynomial& poly ){
      compute( poly ); }

    inline PolynomialSolver(){}

  public:
    /** \returns the complex roots of the polynomial */
    inline const RootsType& roots() const { return m_roots; }

  public:
    /** Clear and fills the back insertion sequence with the real roots of the polynomial
     * i.e. the real part of the complex roots that have an imaginary part which
     * absolute value is smaller than absImaginaryThreshold.
     * absImaginaryThreshold takes the dummy_precision associated
     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
     * */
    template<typename Stl_back_insertion_sequence>
    inline void realRoots( Stl_back_insertion_sequence& bi_seq,
        const RealScalar& absImaginaryThreshold = dummy_precision<Scalar>() ) const
    {
      bi_seq.clear();
      for( int i=0; i<_Deg; ++i )
      {
        if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold ){
          bi_seq.push_back( m_roots[i].real() ); }
      }
    }

  protected:
    template<typename squaredNormBinaryPredicate>
    inline const RootType& selectComplexRoot_withRespectToNorm( squaredNormBinaryPredicate& pred )
    {
      int res=0;
      RealScalar norm2 = m_roots[0].squaredNorm();
      for( int i=1; i<_Deg; ++i )
      {
        const RealScalar currNorm2 = m_roots[i].squaredNorm();
        if( pred( currNorm2, norm2 ) ){
          res=i; norm2=currNorm2; }
      }
      return m_roots[res];

    }

  public:
    /**
     * \returns the complex root with greatest norm.
     */
    inline const RootType& greatestRoot() const
    { return selectComplexRoot_withRespectToNorm( std::greater<Scalar>() ); }

    /**
     * \returns the complex root with smallest norm.
     */
    inline const RootType& smallestRoot() const
    { return selectComplexRoot_withRespectToNorm( std::less<Scalar>() ); }

  protected:
    template<typename squaredRealPartBinaryPredicate>
    inline const RealScalar& selectRealRoot_withRespectToAbsRealPart(
        squaredRealPartBinaryPredicate& pred,
        bool& hasArealRoot,
        const RealScalar& absImaginaryThreshold = dummy_precision<Scalar>() ) const
    {
      hasArealRoot = false;
      int res=0;
      RealScalar abs2;

      for( int i=0; i<_Deg; ++i )
      {
        if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold )
        {
          if( !hasArealRoot )
          {
            hasArealRoot = true;
            res = i;
            abs2 = m_roots[i].real() * m_roots[i].real();
          }
          else
          {
            const RealScalar currAbs2 = m_roots[i].real() * m_roots[i].real();
            if( pred( currAbs2, abs2 ) )
            {
              abs2 = currAbs2;
              res = i;
            }
          }
        }
        else
        {
          if( ei_abs( m_roots[i].imag() ) < ei_abs( m_roots[res].imag() ) ){
            res = i; }
        }
      }
      return m_roots[res].real();
    }


    template<typename RealPartBinaryPredicate>
    inline const RealScalar& selectRealRoot_withRespectToRealPart(
        RealPartBinaryPredicate& pred,
        bool& hasArealRoot,
        const RealScalar& absImaginaryThreshold = dummy_precision<Scalar>() ) const
    {
      hasArealRoot = false;
      int res=0;
      RealScalar val;

      for( int i=0; i<_Deg; ++i )
      {
        if( ei_abs( m_roots[i].imag() ) < absImaginaryThreshold )
        {
          if( !hasArealRoot )
          {
            hasArealRoot = true;
            res = i;
            val = m_roots[i].real();
          }
          else
          {
            const RealScalar curr = m_roots[i].real();
            if( pred( curr, val ) )
            {
              val = curr;
              res = i;
            }
          }
        }
        else
        {
          if( ei_abs( m_roots[i].imag() ) < ei_abs( m_roots[res].imag() ) ){
            res = i; }
        }
      }
      return m_roots[res].real();
    }

  public:
    /**
     * \returns the real root with greatest absolute magnitude.
     * A real root is defined as the real part of a complex root with absolute imaginary
     * part smallest than absImaginaryThreshold.
     * absImaginaryThreshold takes the dummy_precision associated
     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
     * If no real root is found the boolean hasArealRoot is set to false and the real part of
     * the root with smallest absolute imaginary part is returned instead.
     * @param[out] hasArealRoot: boolean true if a real root is found according to the
     * absImaginaryThreshold criterion, false otherwise.
     * @param[in] absImaginaryThreshold: threshold on the absolute imaginary part to decide
     * whether or not a root is real.
     */
    inline const RealScalar& absGreatestRealRoot(
        bool& hasArealRoot,
        const RealScalar& absImaginaryThreshold = dummy_precision<Scalar>() ) const
    { return selectRealRoot_withRespectToAbsRealPart( std::greater<Scalar>(), hasArealRoot, absImaginaryThreshold ); }


    /**
     * \returns the real root with smallest absolute magnitude.
     * A real root is defined as the real part of a complex root with absolute imaginary
     * part smallest than absImaginaryThreshold.
     * absImaginaryThreshold takes the dummy_precision associated
     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
     * If no real root is found the boolean hasArealRoot is set to false and the real part of
     * the root with smallest absolute imaginary part is returned instead.
     * @param[out] hasArealRoot: boolean true if a real root is found according to the
     * absImaginaryThreshold criterion, false otherwise.
     * @param[in] absImaginaryThreshold: threshold on the absolute imaginary part to decide
     * whether or not a root is real.
     */
    inline const RealScalar& absSmallestRealRoot(
        bool& hasArealRoot,
        const RealScalar& absImaginaryThreshold = dummy_precision<Scalar>() ) const
    { return selectRealRoot_withRespectToAbsRealPart( std::less<Scalar>(), hasArealRoot, absImaginaryThreshold ); }


    /**
     * \returns the real root with greatest value.
     * A real root is defined as the real part of a complex root with absolute imaginary
     * part smallest than absImaginaryThreshold.
     * absImaginaryThreshold takes the dummy_precision associated
     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
     * If no real root is found the boolean hasArealRoot is set to false and the real part of
     * the root with smallest absolute imaginary part is returned instead.
     * @param[out] hasArealRoot: boolean true if a real root is found according to the
     * absImaginaryThreshold criterion, false otherwise.
     * @param[in] absImaginaryThreshold: threshold on the absolute imaginary part to decide
     * whether or not a root is real.
     */
    inline const RealScalar& greatestRealRoot(
        bool& hasArealRoot,
        const RealScalar& absImaginaryThreshold = dummy_precision<Scalar>() ) const
    { return selectRealRoot_withRespectToRealPart( std::greater<Scalar>(), hasArealRoot, absImaginaryThreshold ); }


    /**
     * \returns the real root with smallest value.
     * A real root is defined as the real part of a complex root with absolute imaginary
     * part smallest than absImaginaryThreshold.
     * absImaginaryThreshold takes the dummy_precision associated
     * with the _Scalar template parameter of the PolynomialSolver class as the default value.
     * If no real root is found the boolean hasArealRoot is set to false and the real part of
     * the root with smallest absolute imaginary part is returned instead.
     * @param[out] hasArealRoot: boolean true if a real root is found according to the
     * absImaginaryThreshold criterion, false otherwise.
     * @param[in] absImaginaryThreshold: threshold on the absolute imaginary part to decide
     * whether or not a root is real.
     */
    inline const RealScalar& smallestRealRoot(
        bool& hasArealRoot,
        const RealScalar& absImaginaryThreshold = dummy_precision<Scalar>() ) const
    { return selectRealRoot_withRespectToRealPart( std::less<Scalar>(), hasArealRoot, absImaginaryThreshold ); }

  protected:
    EigenSolverType         m_eigenSolver;
    RootsType               m_roots;
};

#endif // EIGEN_POLYNOMIAL_SOLVER_H

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2010 Manuel Yguel <manuel.yguel@xxxxxxxxx>
//
// Eigen is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 3 of the License, or (at your option) any later version.
//
// Alternatively, you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation; either version 2 of
// the License, or (at your option) any later version.
//
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License and a copy of the GNU General Public License along with
// Eigen. If not, see <http://www.gnu.org/licenses/>.

#include "main.h"
#include <Eigen/Polynomials>
#include <iostream>

using namespace std;

template<int Size>
struct ei_increment_if_fixed_size
{
  enum {
    ret = (Size == Dynamic) ? Dynamic : Size+1
  };
};

template< typename T >
struct CurrComplex
{
  typedef typename std::complex<typename NumTraits<typename T::Scalar>::Real> C;
};

template <typename Polynomials>
typename CurrComplex<Polynomials>::C eval( const Polynomials& poly, const typename CurrComplex<Polynomials>::C& root )
{
  typename CurrComplex<Polynomials>::C val=poly[poly.size()-1];
  for( int i=poly.size()-2; i>=0; --i )
  {
    val *= root; val += poly[i];
  }
  return val;
}

template<typename _Scalar, int _Deg> void polynomialsolver(int deg)
{
  typedef ei_increment_if_fixed_size<_Deg>            Dim;
  typedef Matrix<_Scalar,Dim::ret,1>                  PolynomialType;
  typedef PolynomialSolver<_Scalar, _Deg >            PolynomialSolverType;
  typedef typename PolynomialSolverType::RootsType    RootsType;

  typedef Matrix<_Scalar,_Deg,1>                      EvalRootsType;

  PolynomialType pols = PolynomialType::Random(deg+1);
  PolynomialSolverType psolve( pols );
  const RootsType& roots( psolve.roots() );
  EvalRootsType evr( deg );
  for( int i=0; i<roots.size(); ++i )
  {
    evr[i] = std::abs( eval( pols, roots[i] ) );
  }
  cout << evr.transpose() << endl;
  VERIFY( evr.isZero( test_precision<_Scalar>() ) );
}

template<typename _Scalar> void polynomialsolver_scalar()
{
  CALL_SUBTEST_2( (polynomialsolver<_Scalar,2>(2)) );
  CALL_SUBTEST_3( (polynomialsolver<_Scalar,3>(3)) );
  CALL_SUBTEST_4( (polynomialsolver<_Scalar,4>(4)) );
  CALL_SUBTEST_5( (polynomialsolver<_Scalar,5>(5)) );
  CALL_SUBTEST_6( (polynomialsolver<_Scalar,6>(6)) );
  CALL_SUBTEST_7( (polynomialsolver<_Scalar,7>(7)) );
  CALL_SUBTEST_8( (polynomialsolver<_Scalar,17>(17)) );

  int deg = ei_random<int>(17,1000);
  CALL_SUBTEST_9( (polynomialsolver<_Scalar,Dynamic>(deg)) );
}

void test_polynomialsolver()
{
  for(int i = 0; i < g_repeat; i++)
  {
    polynomialsolver_scalar<double>();
    polynomialsolver_scalar<float>();
  }

  //CALL_SUBTEST_5( eigensolver_verify_assert<MatrixXf>() );
}