Dear all,
Thank you for the suggestions and pointers. It looks like I'll need
to do some interface with Fortran/C code. The igraph package seems to
provide an interface with ARPACK, albeit not with the options I need,
so it could be a good starting point.
Best regards,
baptiste
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@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@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@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@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.
