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

[Manuel Yguel <manuel.yguel@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.
thanks, that what I guessed too.
>
>> 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.
>
>> 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.
Sorry I was not clear enough, so let me explain it again:
a companion matrix (balanced) needs 2d-1 coefficients for a polynomial
of degree d.
There is a big memory gain, by not building an entire matrix then
copying it, but just initializing the matrix inside the eigensolver
algorithm by the 2d-1 coefficients.
After, it is up to the eigensolver to do whatever is needed with its matrix..
I just do not want to build an entire dxd matrix that will just be
copied afterward.

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