by **Julien Langou** » Wed Oct 21, 2009 4:30 pm

Here are three emails we had with Peter Strobach. The emails are in

chronological order. (So the first email is the orginal email from Peter,

then you will read my answer, finally his answer.)

To put it short, Peter is reporting than using close form formula to solve

degree 4 polynomials ("quartic polynomials") is faster and more stable

than DGEEV + companion matrix.

>> From: "Peter Strobach"

>> Date: October 21, 2009 3:52:31 AM MDT

>> Subject: Good-Quality Solver for Quartic Equation

>>

>> Dear Julie,

>>

>> I have seen your comment on the above topic "good-quality solver for

>> quartic equation" in the Lapack/Scalapack forum user discussion.

>>

>> You wrote on Aug. 08, 2006: "Lapack does not have any 'direct' method

>> for calculating the roots of a real, quartic equation, but ...''

>>

>> and then you recommend the classical companion matrix method as an

>> alternative to the closed-form solution.

>>

>> Now it appears that I have a larger optimization problem where I have

>> to solve very many quartic equations, literally tens of thousands of

>> them, in one of these larger problems.

>>

>> Therefore, for me, it is a pressing issue, and the solution of the

>> quartic, in my case, is time and accuracy critical, as this

>> dominates runtimes and overall accuracy of my larger optimization

>> program.

>>

>> In 2008, I began writing a code (in double precision Fortran) for

>> the closed-form quartic solver. Meanwhile, this program is well

>> developed.

>>

>> Now, as I have seen your comment in this forum, I had the desire

>> to run an experiment, how fast and accurate a companion matrix

>> type (iterative) quartic solver based on unsymmetric eigenvalue

>> computation might be, in direct comparison with my closed-form

>> quartic solver.

>>

>> Thus, I readily implemented the companion type solver and used

>> the Lapack unsymmetric eigenvalue routine DGEEV.f (computation

>> of eigenvectors disabled, of course).

>>

>> The experiment was carried out on a Acer TravelMate 6592G, my favourite

>> simulation computer, with processor Intel T9300 dual core, 2.5 GHz,

>> 800 MHz FSB, 6 MB L2 cache. The compilier was the latest version

>> of the Intel vectorizing Fortran compiler.

>>

>> My test program comprises a random generation of the quartic

>> coefficients and then the call of the solver routines.

>>

>> The test amounts the solution of 100.000 (one hundred thousand)

>> random quartics. The runtime was accumulated for the solver subroutines

>> using the cpu_time command in Fortran.

>>

>> The following overall cpu_time results were obtained this morning

>> for solving 100.000 random quartics on this TravelMate 6592G:

>>

>> - Closed-form quartic solver (Strobach): overall cpu_time = 0.25 sec

>> - Companion/Lapack eigenvalue solver : overall cpu_time = 2.78 sec

>>

>> You can see - quite striking differences: The closed form quartic

>> solver is more than 10 times faster than the Companion/Eigenvalue

>> based routine.

>>

>> And also the accuray of the solution was overall significantly better

>> for the closed-form solver when compared with the eigenvalue based

>> solver, where the overall accuracy was unsatisfactory.

>>

>> For me, this is very important, because the quality of my overall

>> optimization routine depends crucially on these two factors:

>> Runtime and Accuracy of the underlying quartic solver.

>>

>> Best regards,

>> Peter Strobach, Ph.D.

> Datum: Wed, 21 Oct 2009 08:58:05 -0600 (MDT)

> Von: Julien Langou

> An: Peter Strobach

> CC: Julie Langou

> Betreff: Re: Fwd: Good-Quality Solver for Quartic Equation

>

> Dear Peter,

>

> Thanks for your comment. Yes, there is close form solution for any

> polynomial of degree 4. (That stops here though! You wont do this for

> degree 5 or higher!) In the nonsymmetric eigenvalue solver, LAPACK stops

> at degree 2: when we have a 2-by-2 block diagonal, we get the eigenvalues

> using the discriminant formula. We could theoretically have stopped before

> at degree 4 as you mentionned or later at degree 1. Odd degree (1 and 3)

> being not particularly welcome since we want to leave complex conjugate

> pair in pair. (To avoid complex arithmetic.)

>

> It is good to know that the close formula for the quartic equation is

> faster than LAPACK on the companion matrix. Thanks for reporting your

> experiment to us. Yes a factor of 10 is relevant. The stability report is

> even more interesting. At least unknown to me. I am not sure what we will

> do with this. The software implication seems quite heavy. But it is good

> to know for sure.

>

> For information, current research on 'numerical linear algebra' and

> 'polynomial root finding' is on: how to develop algorithms that exploit

> the structure of the companion matrix. Exploiting the structure of the

> companion should enable to have faster and more stable root finding

> algorithms (for general n).

>

> May I forward your to the lapackers community? and post it on our forum as

> well? I believe what you are saying is relevant for our community.

>

> Best wishes,

> Julien and Julie.

>

Date: Wed, 21 Oct 2009 09:31:21 -0600

To: "Langou, Julien"

Subject: Re: Fwd: Good-Quality Solver for Quartic Equation

Dear Julien and Julie,

there is no problem in forwarding this report to the lapackers

community. I have seen there are relatively many views on this

subject matter.

I can see your argumentation below. With 2 x 2 blocks you can

stay within the real Schur form, avoiding complex arithmetic.

Avoiding complex arithmetic: Sometimes this is nice and sometimes

this is not wise ...

(for instance, my closed-form solver uses double complex in a

versatile way at some places ...)

Anyway, an old wisdom counts: Don't do anything iteratively

before you haven't done it in closed form (if you can)!

It was only a subproblem in a larger problem which is joint

matrix diagonalization: Finding the best generalized eigenvectors

for a set of more than two matrices, which is a challenging

nonlinear optimization problem that can be handled very nicely

using the unit-determinant constraint.

Yes, this is enlightening: To discover how restrictive it is to think

in the world of the orthonormality constraint, and how much freedom

one gains as soon as one arrives in the world of the unit-determinant

constraint!

All the best,

Peter.