From: Martin Maechler <maechler_at_stat.math.ethz.ch>

Date: Thu, 19 Jun 2008 19:05:57 +0200

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 - 17:15:01 GMT

Date: Thu, 19 Jun 2008 19:05:57 +0200

>>>>> "KateM" == Katharine Mullen <kate_at_few.vu.nl> >>>>> 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.

In the mean time (I've suffered from a computer crash and other interruptions) I see Gabor has already mentioned the following

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

>> > 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_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 - 17:15:01 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 - 19:30:46 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.
*