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

[Benoit Jacob <jacob.benoit.1@xxxxxxxxx>]

2010/1/18 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.

Yes, this is what Gael is pointing out when he says that the companion
matrix is already hessenberg. We need to provide an API for
eigenvalues of a Hessenberg matrix. In the same vein, right now i'm
adding SVD of a bidiagonal matrix.

Benoit

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