[R] highest eigenvalues of a matrix

[baptiste Auguié]

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