From: Simon Wood <s.wood_at_bath.ac.uk>

Date: Thu, 19 Jun 2008 10:25:23 +0100

Date: Thu, 19 Jun 2008 10:25:23 +0100

*>
*

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

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

--Received on Thu 19 Jun 2008 - 09:24:29 GMT

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

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Thu 19 Jun 2008 - 10:31:16 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*