[SciPy-dev] [Fwd: Re: wrapping BLOPEX for SciPy]

[Robert Cimrman]

Hi all,

FYI: I have asked Andrew Knyazev, the author of BLOPEX (sparse)
eigenvalue solver, some questions related to making it usable from SciPy
(which I plan to do, unless somebody else did/does it already). Some
interesting points have appeared in the discussion, see below remarks of
Julien Langou, one of the LAPACK developers (the discussion leading to
those remarks is below below :).

r.


-------- Original Message --------
Subject: Re: wrapping BLOPEX for SciPy
Date: Wed, 21 Mar 2007 14:28:24 -0600 (MDT)
From: Julien Langou <langou at math.cudenver.edu>
To: Andrew Knyazev <andrew.knyazev at cudenver.edu>
CC: Robert Cimrman <cimrman3 at ntc.zcu.cz>


Hello,

Some points from the top of my heads for some SciPY developers.

* while Matlab eig correctly call SYEV (symmetric eigenvalue solver) or
GEEV (nonsymmetric eigenvalue solver) from LAPACK depending whether the
matrix is symmetric or not, Matlab eigs call ARPACK whether the matrix is
symmetric or not.

I always thought this was a big liability for such a software. Symmetric
eigensolver would be way faster, more stable, eigenvectors would have
guaranteed orthogonality and eigenvalues would be real.

Below a Matlab script that shows that when there is a cluster of
eigenvalue, eigs (ARPACK) returns the correct invariant subspace but the
eigenvectors are not orthogonalized. (logic for a nonsymmetric eigensolver
...)

	>> clear
	>> n=10; V=orth(randn(n)); D=1:n; D(n-1)=n; A=V*diag(D)*V';
	>> opts.disp=0;
	>> [Veigs,Deigs] = eigs(A,2,'LM',opts);
	>> Deigs

	Deigs =

	    10     0
	     0    10

	>> Veigs'*Veigs

	ans =

	    1.0000    0.9545
	    0.9545    1.0000

Conclusion having a true symmetric eigensolver (e.g. lobpcg) for the SciPy
eigs would be clearly a big plus with respect to Matlab.

* They are four LAPACK symmetric eigensolvers. SYEV (symmetric QR), SYEVD
(divide and conquer), SYEVD (divide and conquer), SYEVX (bisection and
inverse iteration) and SYEVR (MRRR). Matlab uses SYEV. The most robust but
by far the slowest. We have worked on all those algoirthms and they are
now all robust. For MRRR, this is new. The version in 3.1 is really
stable, the one in 3.0, well .... In term of performance, there is no
comparison SYEVR and SYEVX are way faster than SYEV. Performance of SYEVX
are really matrix dependent. Comparing SYEVD and SYEVR is hard.

* SYEVX and SYEVR supports subset computation (same thing: subset
computation in SYEVR is from 3.1). So in some sense if you really want to
compute only 10 of the eigenvalue of a matrix with LAPACK, you can save a
good bunch of FLOPS. This option is not available from Matlab. Note that
the cost is still O(n^3) (due to tridiagonalization).

* LAPACK-3.1 nonsymmetric eigensolver is about 3-5 times faster than the
LAPACK-3.0 one. Matlab has still not yet upgraded to 3.1, as far as I
know.

* For the question whether to use eigs or eig; or equivallently when to
use ARPACK or LAPACK. This depend on the size of the matrix, the sparsity
of the matrix and the number of eigenvalues desired. Just to check
attached is a plot of comparison betwenn eig and eigs for dense matrices
going from n=100:100:1000, for k=[1 10].
Conclusion: eigs behaves as an N^2 algorithm as expected, eig behaves as a
N^3 algorithm as expected. On my machine, if you want one eigenvalue, you
better off with LAPACK for n<150. If you want 10 eigenvalues, you better
off with LAPACK for n<450. Note that the scale is logarithmic so when
n=1000, the n^3 factor completely kills LAPACK and it is about 5 times
slower than eigs.

-j



On Wed, 21 Mar 2007, Andrew Knyazev wrote:

> Robert Cimrman wrote:
>> Andrew Knyazev wrote:
>>> LAPACK is indeed for full matrices only. If you want all 
>>> eigenvectors, the matrix of all eigenvectors is usually full even
>>> if your original matrix is sparse. Since you loose the sparsity
>>> at the end anyway, the standard procedure is to convert your
>>> sparse input matrix into full and plug it into LAPACK. Or may be
>>> you try to do something non orthodox that I do not know.
>> 
>> OK. What SciPy needs is an equivalent of Matlab's eigs. I have
>> briefly looked at it, and it switches to eig( full( A ) ) only if
>> all eigenvalues are requested. I would say that if one wants e.g. k
>> = 999 eigenvalues out of n = 1000, making a full matrix would still
>> be preferable. Is there any recommendation on for which k to switch
>> to the full matrix mode, especially when using BLOPEX?
>> 
>> Thank you for your time and answers!
>> 
>> r.
> 
> I never did a comparison between MATLAB's eig (which is an interface
> to LAPACK) and eigs (which is an interface to ARPACK). I would guess
> that under normal circumstances, e.g., having enough memory, eig
> should be faster than eigs if >50% eigenvectors are requested, may be
> less. With BLOPEX, I would guess that eig is faster for >10% (if no
> preconditioning is used). The main advantage of BLOPEX is that it
> allows using preconditioning to accelerate the convergence. The
> preconditioning must be provided by the user, though.
> 
> I copy Julien Langou, one of the LAPACK developers, who may be able
> to add some more to this.
>