[R] highest eigenvalues of a matrix

[Katharine Mullen]

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.
>

ARPACK (http://www.caam.rice.edu/software/ARPACK/) uses a Lanczos method
for symmetric matrics; otherwise it uses an Arnoldi iteration.
Development of an R interface to ARPACK would be a nice project (but
unfortunately one I don't have time for for a while).  Maybe one of the
maintainers of a package for sparse matrices would be interested.

> 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?
> >

dsyevr.f will work with symmetric real matrices only.  When the range
argument of dysevr is set to 'I', arguments il and iu seem to specify the
range of eigenvalue indices you want in ascending order (lowest to
highest, not highest to lowest).  If you look at
https://svn.r-project.org/R/trunk/src/modules/lapack/Lapack.c you see that
range is always set to 'A'.

> > 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
> >
> > ______________________________________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> > http://www.R-project.org/posting-guide.html and provide commented, minimal,
> > self-contained, reproducible code.
>
> --
> > Simon Wood, Mathematical Sciences, University of Bath, Bath, BA2 7AY UK
> > +44 1225 386603  www.maths.bath.ac.uk/~sw283
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>