[SciPy-Dev] ARPACK fixes before 0.9

[David Cournapeau]

On Sun, Nov 28, 2010 at 1:33 AM, Pauli Virtanen <pav at iki.fi> wrote:
> Hi,
>
> The Arpack interface in Scipy needs some fixes, which would be useful to
> get in for Scipy 0.9. I have them in a branch here:
>
>        https://github.com/pv/scipy-work/tree/bug/1313-arpack
>
> Comments would be appreciated (esp. from David who wrote the original
> interface). I'm going to merge them probably the next W/E unless
> objections arise.
>
> Bugfixes:
>
> - Raise ArpackNoConvergence when ARPACK iteration does not converge.
>  Previously, there was no way to know if convergence was obtained.
>
> - Fix a bug in return value extraction from dneupd, which resulted
>  to invalid eigenvectors/eigenvalues being returned on non-convergence.
>
> - Allow complex matrices in the sparse SVD. It's a naive approach,
>  but probably still better than nothing.

I think so too, that's why I put SVD in the first place. When I was
looking at sparse svd implementations, the best solution looked like
adding support for PROPACK (which is much faster than ARPACK, up to 10
times on some middle-size problems I tried), but there is no license
from the author (and any known attempt to get clarification has been
ignored AFAIK).

The PROPACK method has been fairly well documented, though, with both
fortran and (early) matlab implementation, so it is doable. The
original code cannot be used as is anyway because of some hardcoded
limitations.

Besides speed, the main issue with solving A^h A is precision, I was
thinking about adding support for solving [A 0; 0 A^H] which is more
precise (but much slower), but never got around it. Should be fairly
easy to add.

>
> Other changes:
>
> - Rename routines eigen* -> eigs* to avoid shadowing the "eigen" module.
>  "eigs" is also the name for the equivalent routines in Octave/Matlab.
>
> - Remove the ``speigs`` ARPACK interface. It does not seem to make sense
>  two have two different interfaces to the same library.

agreed - I only avoid removing anything at that time because I did not
know what was the status of the interface (and needed svd quickly).

For the parts I am familiar with, the changes made sense - but as Aric
mentioned, I have only contributed the svd part. There are still some
hard (causing crash) bugs in the arpack submodule, but I suspect it is
a bug in arpack itself: when I tried to debug it, it was out of bound
array indexing in some cases, which did not seem to be caused by
invalid input. I think I would rather spend time implementing better
algorithms than debugging unstrucuted, goto-ridden fortran code :)

cheers,

David