[R] Solve an ordinary or generalized eigenvalue problem in R

[John C Nash]

This thread reveals that R has some holes in the solution of some of the linear algebra 
problems that may arise. It looks like Jim Ramsay used a quick and dirty approach to the 
generalized eigenproblem by using B^(-1) %*% A, which is usually not too successful due to 
issues with condition of B and making a symmetric/Hermitian problem unsymmetric.

In short, the problem is stated as follows:

  Find the eigenvalues e and vectors v that satisfy

     A v  =  e B v

  There is a matrix form (I think   A V = B V E is usual, though I tend to work on them 
one at a time in my mind.)

The real trouble is that A and B can have different forms e.g., real, complex, Hermitian
and B may or may not be positive definite. I published a code for complex Hermitian case 
with posdef (but complex) metric B in 1974 in Computer Physics Communications. Around the 
same time there was Linda Kaufman's LZ code (TOMS / ACM 496) and various QZ codes (van 
Loan and Smith), see TOMS / ACM 535. There are descendants of these in LAPACK. And I 
suspect there are some sparse variants too.

The "nice" place for these would be in the Matrix package, but as a shorter-term solution, 
it could be useful to put together a suite for dense matrices. Using existing Fortran 
codes would allow this to be done fairly quickly and it could be a nice summer or term 
project for a student (Google Summer of Code 2013?). I'll be happy to kibbitz and mentor, 
but I'd rather stay on the sidelines of actual package building, as I no longer have a 
direct need. My application was the quantum electronic structure of polymers in 1972/3. 
However, I think the code is still reasonable.

Best, JN


On 04/23/2012 06:00 AM, r-help-request at r-project.org wrote:
> Message: 36 Date: Sun, 22 Apr 2012 21:27:23 +0200 From: Berend Hasselman <bhh at xs4all.nl>
> To: Jonathan Greenberg <jgrn at illinois.edu> Cc: R Project Help <r-help at r-project.org>
> Subject: Re: [R] Solve an ordinary or generalized eigenvalue problem in R? Message-ID:
> <378C084B-5CF5-4087-9ABA-116C5155D12E at xs4all.nl> Content-Type: text/plain;
> charset=us-ascii On 22-04-2012, at 21:08, Jonathan Greenberg wrote:
>> >  Thanks all (particularly to you, Berend) -- I'll push forward with these solutions and integrate them into my code.  I did come across geigen while rooting around in the CCA code but its not formally documented (it just says "for internal use" or something along those lines) and as you found out above, it does not produce the same solution as the dggev.  It would be nice to have a more complete set of formal packages for doing LA in R (rather than having to hand-write .Fortran calls) but I'll leave that to someone with more expertise in linear algebra than me.  Something that perhaps matches the SciPy set of functions (both in terms of input and output):
>> >
>> >  http://docs.scipy.org/doc/scipy/reference/linalg.html
>> >
>> >  Some of these are already implemented, but clearly not all of them.
> Package CCA has package fda as dependency.
> And package fda defines a function geigen.
> The first 14 lines of this function are
>
> geigen<- function(Amat, Bmat, Cmat)
> {
>    #  solve the generalized eigenanalysis problem
>    #
>    #    max {tr L'AM / sqrt[tr L'BL tr M'CM] w.r.t. L and M
>    #
>    #  Arguments:
>    #  AMAT ... p by q matrix
>    #  BMAT ... order p symmetric positive definite matrix
>    #  CMAT ... order q symmetric positive definite matrix
>    #  Returns:
>    #  VALUES ... vector of length s = min(p,q) of eigenvalues
>    #  LMAT   ... p by s matrix L
>    #  MMAT   ... q by s matrix M
>
> It's not clear to me how it is used and exactly what it is doing and how that compares with Lapack.
>
> Berend
>