Thank you all for your interest in this. The ARPACK library looks
like a very interesting option, I wish I had some experience in
interfacing R with such libraries. I'll try to learn on a simple
example first, and if successful I will try to adapt the igraphlib
interface to provide the desired options ( for the record, the
original problem was to find the first 3 eigenvalues (decreasing
magnitude) and corresponding eigenvectors of a 200x200 Hermitian
matrix by an iterative process). Maybe a more general approach would
be better, but I'm already definitely out of my league here (not to
mention the incompatibility issue with LAPACK already mentioned).
Many thanks,
baptiste
On 19 Jun 2008, at 19:08, Katharine Mullen wrote:
On Thu, 19 Jun 2008, Martin Maechler wrote:
"KateM" == Katharine Mullen <[EMAIL PROTECTED]>
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
______________________________________________
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.
