R help - Apr 2012 - Solve an ordinary or generalized eigenvalue problem in R

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
>