[SciPy-dev] Ticket #709 (was: Re: Next scipy release 0.7)

[Pauli Virtanen]

Tue, 30 Sep 2008 19:56:58 +0200, Nils Wagner wrote:
> IMHO,
> 
> http://projects.scipy.org/scipy/scipy/ticket/709
> 
> should be resolved, too.

Short summary about this one: it's about a special matrix pair for which 
LAPACK's generalized eigenvalue problem solver DGGEV returns bogus 
output, losing one eigenvalue from a complex eigenpair. (The eigenvalue 
problem should nonetheless be well-defined, AFAIK, this looks like a bug 
in LAPACK.)

The question here is what we should/can do:

1. Raise LinAlgError if we detect this condition.

2. Try to repair the returned data by filling in the other eigenvalue of
   the pair, as we know all complex eigenvalues come in pairs.

3. Try to return as much as we can make sense of, but put NaN in place of
   the missing eigenvalue.


Details:

Instead of expected complex eigenvalue pair (a_R +- i a_I)/beta, the 
LAPACK routine instead returns the values

	a_R  = whatever
	a_I  = [0., -2.18645248e-16j]
	beta = [0., 7.95804395e-17]

The correct result would have been something like

	a_R  = whatever
	a_I  = [+2.18645248e-16j, -2.18645248e-16j]
	beta = [7.95804395e-17, 7.95804395e-17]

but one eigenvalue of the pair was lost, apparently due to rounding etc.

Eigenvector data corresponding to zeros in a_I are interpreted 
differently than for positive entries, so the zero in the wrong place 
messed up the Scipy routine. We could detect this and raise an error (1), 
substitute in the missing eigenvalue based on the correctly calculated 
one (2), or leave the missing eigenvalue as corrupted and return only the 
correct one (3).

-- 
Pauli Virtanen