From: baptiste Auguié <ba208_at_exeter.ac.uk>

Date: Thu, 19 Jun 2008 15:12:50 +0100

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 - 14:27:23 GMT

Date: Thu, 19 Jun 2008 15:12:50 +0100

On 19 Jun 2008, at 10:50, Katharine Mullen wrote:

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

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 - 14:27:23 GMT

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 - 15:00:56 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.
*