Re: [R] highest eigenvalues of a matrix

From: Martin Maechler <maechler_at_stat.math.ethz.ch>
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

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

  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 ?

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

list of date sections of archive