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-requ...@r-project.org wrote:
Message: 36 Date: Sun, 22 Apr 2012 21:27:23 +0200 From: Berend Hasselman
<b...@xs4all.nl>
To: Jonathan Greenberg <j...@illinois.edu> Cc: R Project Help
<r-help@r-project.org>
Subject: Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
Message-ID:
<378c084b-5cf5-4087-9aba-116c5155d...@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
______________________________________________
R-help@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.
