[R] highest eigenvalues of a matrix

[Katharine Mullen]

On Thu, 19 Jun 2008, Martin Maechler wrote:

> >>>>> "KateM" == Katharine Mullen <kate at few.vu.nl>
> >>>>>     on Thu, 19 Jun 2008 11:50:22 +0200 (CEST) writes:
>
>     KateM> On Thu, 19 Jun 2008, Simon Wood wrote:
>     >> >
>     >> > I happily use eigen() to compute the eigenvalues and
>     >> eigenvectors of > a fairly large matrix (200x200, say),
>     >> but it seems over-killed as its > rank is limited to
>     >> typically 2 or 3. I sort of remember being taught > that
>     >> numerical techniques can find iteratively decreasing
>     >> eigenvalues > and corresponding orthogonal eigenvectors,
>     >> which would provide a nice > alternative (once I have the
>     >> first 3, say, I stop the search).
>     >>
>     >> Lanczos iteration will do this efficiently (see
>     >> e.g. Golub & van Loan "Matrix Computations"), but I don't
>     >> think that there are any such routines built into R or
>     >> LAPACK (although I haven't checked the latest LAPACK
>     >> release). When I looked it seemed that the LAPACK options
>     >> that allow you to select eigen-values/vectors still
>     >> depend on an initial O(n^3) decomposition of the matrix,
>     >> rather than the O(n^2) that a Lanczos based method would
>     >> require.
>     >>
>
>      KateM> ARPACK (http://www.caam.rice.edu/software/ARPACK/) uses a
>      KateM> Lanczos method for symmetric matrics; otherwise it uses an
>      KateM> Arnoldi iteration.  Development of an R interface to ARPACK
>      KateM> would be a nice project (but unfortunately one I don't have
>      KateM> time for for a while).  Maybe one of the maintainers of a
>      KateM> package for sparse matrices would be interested.
>
> Yes, thank you, Kate, we have been.
> The nice thing in ARPACK is its concept of "reverse
> communication", such that  we should be able to use it both
> for sparse and for dense matrices.
> We would only have to provide a "function" for computing
> "A %*% x" and hence could use different ones, depending on the
> class of A.... all in theory at least.
>
> One problem I see with the  with that code is that its README file contains
>
>   *** NOTE ***  Unless the LAPACK library on your system is version 2.0,
>    we strongly recommend that you install the LAPACK routines provided with
>    ARPACK. Note that the current LAPACK release is version 3.0; if you are
>    not sure which version of LAPACK is installed, pleaase compile and link
>    to the subset of LAPACK included with ARPACK.
>
> i.e. you need to use old/outdated Lapack with ARPACK instead of
> the current Lapack .  Something that sheds a bad light to ARPACK
> in my view.
>

Dear Martin,

I agree this is not ideal.  I wonder why it has not been kept up to date
with LAPACK and how broken it is with version 3.1.1. (but can't experiment
now).

> In the mean time (I've suffered from a computer crash and other
> interruptions) I see Gabor has already mentioned the following
>
>   Then, we have seen that the 'igraph' package als uses a port of
>   (part of) ARPACK as part of the C++ 'igraphlib' library.
>
>
>
>     >> My `mgcv' package (see cran) uses Lanczos iteration for
>     >> setting up low rank bases for smoothing. The source code
>     >> is in mgcv/src/matrix.c:lanczos_spd, but I'm afraid that
>     >> there is no direct R interface, although it would not be
>     >> too hard to write a suitable wrapper. It requires the
>     >> matrix to be symmetric.
>     >>
>     >> > Looking at the R source code for eigen and some posts
>     >> on this list, > it seems that the function uses a LAPACK
>     >> routine, but obviously all > the options are not
>     >> available through the R wrapper. I'm not > experienced
>     >> enough to try and make my own interface with Fortran >
>     >> code, so here are two questions:
>     >> >
>     >> > - is this option (choosing a desired number of
>     >> eigenvectors) already > implemented in some function /
>     >> package that I missed?  --- In the symmetric case you can
>     >> use `svd' which lets you select (although you'd need to
>     >> fix up the signs of the singular values to get
>     >> eigen-values if the matrix is not +ve definite). But this
>     >> answer is pretty useless as it will be slower than using
>     >> `eigen' and getting the full decomposition.
>     >>
>     >> Of course if you know that your matrix is low rank
>     >> because it's a product of non-square matrices then
>     >> there's usually some way of getting at the
>     >> eigen-decomposition efficiently. E.g. if A=B'B where B is
>     >> 3 by 1000, then the cost can easily be kept down to
>     >> O(1000^2) in R...
>     >>
>     >> best, Simon
>     >>
>     >> > - is the "range of indices" option in DSYEVR.f <
>     >> http:// > www.netlib.org/lapack/double/dsyevr.f > what I
>     >> think, the indices of > the desired eigenvalues ordered
>     >> from the highest to lowest?
>     >> >
>
>     KateM> dsyevr.f will work with symmetric real matrices only.  When the range
>     KateM> argument of dysevr is set to 'I', arguments il and iu
>     KateM> seem to specify the range of eigenvalue indices you
>     KateM> want in ascending order (lowest to highest, not
>     KateM> highest to lowest).  If you look at
>     KateM> https://svn.r-project.org/R/trunk/src/modules/lapack/Lapack.c
>     KateM> you see that range is always set to 'A'.
>
> yes, but do you agree that using 'I' (and a range) does not
> really help much, since the factorization used is the "full" one
> anyway ?
>

Now I see -- I had read

Eigenvalues and eigenvectors can be
*  selected by specifying either a range of values or a range of
*  indices for the desired eigenvalues.

with _selected_ magically replaced by _computed_.  Thanks.

>
>     >> > Many thanks in advance for any piece of advice,
>
>     >> > Sincerely,
>     >> >
>     >> > Baptiste
>     >> >
>     >> > dummy example if needed:
>     >> >
>     >> > test <- matrix(c(1, 2, 0, 4, 5, 6, 1.00001, 2, 0),
>     >> ncol=3) > eigen(test)
>     >> >
>     >> >
>     >> >
>     >> >
>     >> > _____________________________
>     >> >
>     >> > Baptiste Auguié
>     >> >
>     >> > Physics Department > University of Exeter > Stocker
>     >> Road, > Exeter, Devon, > EX4 4QL, UK
>     >> >
>     >> > Phone: +44 1392 264187
>     >> >
>     >> > http://newton.ex.ac.uk/research/emag >
>     >> http://projects.ex.ac.uk/atto
>     [............]
>
>
>     >> --
>     >> > Simon Wood, Mathematical Sciences, University of Bath,
>     >> Bath, BA2 7AY UK > +44 1225 386603
>     >> www.maths.bath.ac.uk/~sw283
>