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

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

On Mon, Jan 18, 2010 at 6:51 PM, Benoit Jacob <jacob.benoit.1@xxxxxxxxx> wrote:
> 2010/1/18 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 ?
>
> If you only want to support real polynomials, then yes.
>
>> 2) Does the eigensolver balance the input matrix ?
>
> Not as far as I can see.
>
>> 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.
>
> Ask Gael to be sure, but I don't think that we have this right now.

Indeed, there is no balancing. BTW, this eigen solver is still the one
coming from Jama, so perhaps it would be interesting to try with
another eigensolver, just to see if that one is really bad or not. For
instance if you have Linux, then you can very quickly try with
octave's eigensolver using the following function (you need to link to
liboctave.so):

#include <octave/config.h>
#include <octave/EIG.h>
#include <octave/dbleDET.h>

template<typename Matrix, typename Eigenvalues>
int eigenvalue_by_octave(const Matrix& m, Eigenvalues& evals)
{
  typedef typename Eigenvalues::Scalar EigenvalueType;
  int size = m.rows();

  // prepare the matrix for octave
  Matrix m_copy(size, size);
  for(int i=0 ; i<size ; ++i)
    for(int j=0 ; j<size ; ++j)
      m_copy.elem(i,j) = m(i,j);

  int info;
  EIG eig(m_copy, info, true);

  // get the eigen values
  for(int i=0 ; i<size ; ++i)
    evals(i) = eig.eigenvalues().elem(i);

  return (info==0 ? 1 : 0);
}

>
>> 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?
>
> In all these decomposition algorithms, the work matrix is initialized
> at the start by copying the input matrix into it. This is how these
> algorithms work. If you want to preserve sparsity, you need a
> completely different algorithm: this one won't preserve sparsity at
> all. But I am not sure that companion matrices offer real
> opportunities to take advantage of sparsity here.

Actually, here we are missing an "inplace" decomposition. Furthermore,
the eigensolver starts by doing a Hessenberg decomposition, and since
the companion matrix is already Hessenberg, we could skip this very
expensive part.

gael

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