Re: [R] 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-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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
Building in Berend's suggestions I think this function should work for most
people (I'm going to wrap it into a package but figured people may want to
grab this directly):
# Please see http://www.netlib.org/lapack/double/dggev.f for a description
of inputs and outputs.
Rdggev - function(A,B,JOBVL=FALSE,JOBVR=TRUE)
{
# R implementation of the DGGEV LAPACK function (with generalized
eigenvalue computation)
# See http://www.netlib.org/lapack/double/dggev.f
# coded by Jonathan A. Greenberg i...@estarcion.net
# Contributions from Berend Hasselman.
if( .Platform$OS.type == windows ) {
Lapack.so -
file.path(R.home(bin),paste(Rlapack,.Platform$dynlib.ext,sep=))
} else {
Lapack.so -
file.path(R.home(modules),paste(lapack,.Platform$dynlib.ext,sep=))
}
dyn.load(Lapack.so)
if(JOBVL)
{
JOBVL=V
} else
{
JOBVL=N
}
if(JOBVR)
{
JOBVR=V
} else
{
JOBVR=N
}
if(!is.matrix(A)) stop(Argument A should be a matrix)
if(!is.matrix(B)) stop(Argument B should be a matrix)
dimA - dim(A)
if(dimA[1]!=dimA[2]) stop(A must be a square matrix)
dimB - dim(B)
if(dimB[1]!=dimB[2]) stop(B must be a square matrix)
if(dimA[1]!=dimB[1]) stop(A and B must have the same dimensions)
if( is.complex(A) ) stop(A may not be complex)
if( is.complex(B) ) stop(B may not be complex)
# Input parameters
N=dim(A)[[1]]
LDA=N
LDB=N
LDVL=N
LDVR=N
LWORK=as.integer(max(1,8*N))
Rdggev_out - .Fortran(dggev, JOBVL, JOBVR, N, A, LDA, B, LDB,
double(N), double(N), double(N),
array(data=0,dim=c(LDVL,N)), LDVL, array(data=0,dim=c(LDVR,N)), LDVR,
double(max(1,LWORK)), LWORK, integer(1))
names(Rdggev_out)=c(JOBVL,JOBVR,N,A,LDA,B,LDB,ALPHAR,ALPHAI,
BETA,VL,LDVL,VR,LDVR,WORK,LWORK,INFO)
# simplistic calculation of eigenvalues (see caveat in
http://www.netlib.org/lapack/double/dggev.f)
if( all(Rdggev_out$ALPHAI==0) )
Rdggev_out$GENEIGENVALUES - Rdggev_out$ALPHAR/Rdggev_out$BETA
else
Rdggev_out$GENEIGENVALUES - complex(real=Rdggev_out$ALPHAR,
imaginary=Rdggev_out$ALPHAI)/Rdggev_out$BETA
return(Rdggev_out)
}
Thanks!
--j
On Sun, Apr 22, 2012 at 2:27 PM, Berend Hasselman b...@xs4all.nl wrote:
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
--
Jonathan A. Greenberg, PhD
Assistant Professor
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 415-763-5476
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
http://www.geog.illinois.edu/people/JonathanGreenberg.html
[[alternative HTML version deleted]]
__
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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
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. --j On Sat, Apr 21, 2012 at 1:31 PM, Berend Hasselman b...@xs4all.nl wrote: On 21-04-2012, at 20:20, peter dalgaard wrote: The eigenvalues are identical upto the printed 9 digits but the eigenvectors appear to be quite different. Maybe this is what Luke meant. Berend They look quite similar to me: ev - eigen(solve(B,A) )$vectors ge - geigen(A, B, TRUE , TRUE) ev / ge$vl [,1] [,2] [,3] [1,] 0.9324603 0.813422 -0.7423694 [2,] 0.9324603 0.813422 -0.7423694 [3,] 0.9324603 0.813422 -0.7423694 ev / ge$vr [,1] [,2] [,3] [1,] 0.9324603 0.813422 -0.7423694 [2,] 0.9324603 0.813422 -0.7423694 [3,] 0.9324603 0.813422 -0.7423694 (and of course, eigenvectors of any sort are only defined up to a constant multiplier) Correct. I should have checked your way and not optically. 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. -- Jonathan A. Greenberg, PhD Assistant Professor Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 607 South Mathews Avenue, MC 150 Urbana, IL 61801 Phone: 415-763-5476 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html [[alternative HTML version deleted]] __ 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On 20-04-2012, at 21:18, Jonathan Greenberg wrote:
Ok, I figured out a solution and I'd like to get some feedback on this from
the R-helpers as to how I could modify the following to be package
friendly -- the main thing I'm worried about is how to dynamically set the
dyn.load statement below correctly (obviously, its hard coded to my
particular install, and would only work with windows since I'm using a
.dll):
Rdggev - function(JOBVL=F,JOBVR=T,A,B)
{
# R implementation of the DGGEV LAPACK function (with generalized
eigenvalue computation)
# See http://www.netlib.org/lapack/double/dggev.f
# coded by Jonathan A. Greenberg j...@illinois.edu
dyn.load(C:\\Program Files\\R\\R-2.14.2\\bin\\x64\\Rlapack.dll)
if(JOBVL)
{
JOBVL=V
} else
{
JOBVL=N
}
if(JOBVR)
{
JOBVR=V
} else
{
JOBVR=N
}
# Note, no error checking is performed.
# Input parameters
N=dim(A)[[1]]
LDA=N
LDB=N
LDVL=N
LDVR=N
LWORK=as.integer(max(1,8*N))
Rdggev_out - .Fortran(dggev, JOBVL, JOBVR, N, A, LDA, B, LDB,
double(N), double(N), double(N),
array(data=0,dim=c(LDVL,N)), LDVL, array(data=0,dim=c(LDVR,N)), LDVR,
double(max(1,LWORK)), LWORK, integer(1))
names(Rdggev_out)=c(JOBVL,JOBVR,N,A,LDA,B,LDB,ALPHAR,ALPHAI,
BETA,VL,LDVL,VR,LDVR,WORK,LWORK,INFO)
Rdggev_out$GENEIGENVALUES=(Rdggev_out$ALPHAR+Rdggev_out$ALPHAI)/Rdggev_out$BETA
return(Rdggev_out)
}
See this:
start R code
# This works on Mac OS X
# Change as needed for other systems
# or compile geigen into a standalone shared object.
dyn.load(file.path(R.home(lib),libRlapack.dylib))
geigen - function(A,B,jobvl=FALSE,jobvr=FALSE) {
# simplistic interface to Lapack dggev
# for generalized eigenvalue problem
# general matrices
# for symmetric matrices use dsygv
if(!is.matrix(A)) stop(Argument A should be a matrix)
if(!is.matrix(B)) stop(Argument B should be a matrix)
dimA - dim(A)
if(dimA[1]!=dimA[2]) stop(A must be a square matrix)
dimB - dim(B)
if(dimB[1]!=dimB[2]) stop(B must be a square matrix)
if(dimA[1]!=dimB[1]) stop(A and B must have the same dimensions)
if( is.complex(A) ) stop(A may not be complex)
if( is.complex(B) ) stop(B may not be complex)
jobvl.char - if(jobvl) V else N
jobvr.char - if(jobvr) V else N
n - dimA[1]
# minimum amount of work memory
# for performance this needs to be set properly
# (see source dggev)
lwork - as.integer(8*n)
work - numeric(lwork)
if( jobvl jobvr )
z - .Fortran(dggev, jobvl.char, jobvr.char, n, A, n, B, n,
alphar=numeric(n), alphai=numeric(n),
beta=numeric(n), vl=matrix(0,nrow=n,ncol=n), n,
vr=matrix(0,nrow=n,ncol=n), n,
work, lwork,info=integer(1L))
else if(jobvl !jobvr )
z - .Fortran(dggev, jobvl.char, jobvr.char, n, A, n, B, n,
alphar=numeric(n), alphai=numeric(n),
beta=numeric(n), vl=matrix(0,nrow=n,ncol=n), n,
numeric(1), 1L,
work, lwork,info=integer(1L))
else if(!jobvl jobvr )
z - .Fortran(dggev, jobvl.char, jobvr.char, n, A, n, B, n,
alphar=numeric(n), alphai=numeric(n),
beta=numeric(n), numeric(1),1L,
vr=matrix(0,nrow=n,ncol=n), n,
work, lwork,info=integer(1L))
else
z - .Fortran(dggev, jobvl.char, jobvr.char, n, A, n, B, n,
alphar=numeric(n), alphai=numeric(n),
beta=numeric(n), numeric(1), 1L, numeric(1), 1L,
work, lwork,info=integer(1L))
if( z$info 0 ) stop(paste(Lapack dggev fails with info=,z$info))
# simplistic calculation of eigenvalues (see caveat in source dggev)
if( all(z$alphai==0) )
values - z$alphar/z$beta
else
values - complex(real=z$alphar, imaginary=z$alphai)/z$beta
if( jobvl jobvr )
return(list(values=values, vl=z$vl, vr=z$vr))
else if( jobvl )
return(list(values=values, vl=z$vl))
else if( jobvr )
return(list(values=values, vr=z$vr))
else
return(list(values=values))
}
A - matrix(c(1457.738, 1053.181, 1256.953,
1053.181, 1213.728, 1302.838,
1256.953, 1302.838, 1428.269), nrow=3, byrow=TRUE)
B - matrix(c(4806.033, 1767.480, 2622.744,
1767.480, 3353.603, 3259.680,
2622.744, 3259.680, 3476.790), nrow=3, byrow=TRUE)
A
B
geigen(A, B)
geigen(A, B, TRUE , TRUE)
geigen(A, B, TRUE , FALSE)
geigen(A, B, FALSE, TRUE)
end R code
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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
Hi all, In my experience, using eigen to solve generalized eigenvalue / eigenvector problems only gives correct looking eigenvalues while the eigenvectors seem to be wrong (in comparison to results from MATLAB's 'eig' function for example). However, I think it is possible to solve generalized eigenvalue / eigenvectors problems in R. The way I found that works in my case is to use the simultaneous diagonalization method which is effectively two applications of the 'svd' function. I have used this myself to get results effectively the same (down to 6th decimal place for example) as those from MATLAB. Kind regards, Luke [[alternative HTML version deleted]] __ 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On 21-04-2012, at 11:40, Berend Hasselman wrote: . See this: start R code # This works on Mac OS X # Change as needed for other systems # or compile geigen into a standalone shared object. dyn.load(file.path(R.home(lib),libRlapack.dylib)) Replacing the dyn.load line with dyn.load(file.path(R.home(modules),lapack.so)) lets run geigen1.R run on Mac OS X and Ubuntu unchanged. Hopefully this is the proper way to load Lapack as provided by R. How to do it on Windows, I can't test. 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On Apr 21, 2012, at 15:22 , Luke Hartigan wrote: Hi all, In my experience, using eigen to solve generalized eigenvalue / eigenvector problems only gives correct looking eigenvalues while the eigenvectors seem to be wrong (in comparison to results from MATLAB's 'eig' function for example). Could you please document that? There are many misconceptions about when eigenvectors are correct and platform dependencies too. As far as I can tell, both R and Matlab use the same LAPACK routines. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
Hello,
Berend Hasselman wrote
On 21-04-2012, at 11:40, Berend Hasselman wrote:
.
See this:
start R code
# This works on Mac OS X
# Change as needed for other systems
# or compile geigen into a standalone shared object.
dyn.load(file.path(R.home(lib),libRlapack.dylib))
Replacing the dyn.load line with
dyn.load(file.path(R.home(modules),lapack.so))
lets run geigen1.R run on Mac OS X and Ubuntu unchanged.
Hopefully this is the proper way to load Lapack as provided by R.
How to do it on Windows, I can't test.
Berend
__
R-help@ 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.
My system is a Windows 7 and the following function solved the dyn.load.
dynlib.load - function(x, dir){
dynname - paste(x, .Platform$dynlib.ext, sep=)
dynname - file.path(R.home(dir), dynname)
dyn.load(dynname)
}
# In windows it's these names
dynlib.load(Rlapack, bin)
The rest worked at the first try.
Note that this function is independent of the sub-architecture, like the
help page for R.home() says:
The return value for modules and on Windows bin is to a
sub-architecture-specific location.
(R.home(bin) returns .../bin/i386 or .../bin/x64)
Hope this helps
Rui Barradas
--
View this message in context:
http://r.789695.n4.nabble.com/Solve-an-ordinary-or-generalized-eigenvalue-problem-in-R-tp4571823p4576615.html
Sent from the R help mailing list archive at Nabble.com.
__
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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On 21-04-2012, at 17:37, Rui Barradas wrote:
Hello,
Berend Hasselman wrote
On 21-04-2012, at 11:40, Berend Hasselman wrote:
.
See this:
start R code
# This works on Mac OS X
# Change as needed for other systems
# or compile geigen into a standalone shared object.
dyn.load(file.path(R.home(lib),libRlapack.dylib))
Replacing the dyn.load line with
dyn.load(file.path(R.home(modules),lapack.so))
lets run geigen1.R run on Mac OS X and Ubuntu unchanged.
Hopefully this is the proper way to load Lapack as provided by R.
How to do it on Windows, I can't test.
Berend
__
R-help@ 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.
My system is a Windows 7 and the following function solved the dyn.load.
dynlib.load - function(x, dir){
dynname - paste(x, .Platform$dynlib.ext, sep=)
dynname - file.path(R.home(dir), dynname)
dyn.load(dynname)
}
# In windows it's these names
dynlib.load(Rlapack, bin)
The rest worked at the first try.
Note that this function is independent of the sub-architecture, like the
help page for R.home() says:
The return value for modules and on Windows bin is to a
sub-architecture-specific location.
(R.home(bin) returns .../bin/i386 or .../bin/x64)
Thank you.
I've now changed the dyn.load line to this
if( .Platform$OS.type == windows ) {
Lapack.so - file.path(R.home(bin),paste0(Rlapack,.Platform$dynlib.ext))
} else {
Lapack.so -
file.path(R.home(modules),paste0(lapack,.Platform$dynlib.ext))
}
dyn.load(Lapack.so)
Hopefully that will do.
BTW: lapack.so on Mac OS X is located in a sub-architecture-specific location.
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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On 21-04-2012, at 16:40, peter dalgaard wrote: On Apr 21, 2012, at 15:22 , Luke Hartigan wrote: Hi all, In my experience, using eigen to solve generalized eigenvalue / eigenvector problems only gives correct looking eigenvalues while the eigenvectors seem to be wrong (in comparison to results from MATLAB's 'eig' function for example). Could you please document that? There are many misconceptions about when eigenvectors are correct and platform dependencies too. As far as I can tell, both R and Matlab use the same LAPACK routines. The OP posted two matrices: A - matrix(c(1457.738, 1053.181, 1256.953, 1053.181, 1213.728, 1302.838, 1256.953, 1302.838, 1428.269), nrow=3, byrow=TRUE) B - matrix(c(4806.033, 1767.480, 2622.744, 1767.480, 3353.603, 3259.680, 2622.744, 3259.680, 3476.790), nrow=3, byrow=TRUE) I've tried eigen(solve(B)%*%A) (which is probably not the best way to handle this problem numerically) to approximate the generalized ev problem. And the geigen function geigen(A,B,TRUE,TRUE) The eigenvalues are identical upto the printed 9 digits but the eigenvectors appear to be quite different. Maybe this is what Luke meant. 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On Apr 21, 2012, at 19:13 , Berend Hasselman wrote: On 21-04-2012, at 16:40, peter dalgaard wrote: On Apr 21, 2012, at 15:22 , Luke Hartigan wrote: Hi all, In my experience, using eigen to solve generalized eigenvalue / eigenvector problems only gives correct looking eigenvalues while the eigenvectors seem to be wrong (in comparison to results from MATLAB's 'eig' function for example). Could you please document that? There are many misconceptions about when eigenvectors are correct and platform dependencies too. As far as I can tell, both R and Matlab use the same LAPACK routines. The OP posted two matrices: A - matrix(c(1457.738, 1053.181, 1256.953, 1053.181, 1213.728, 1302.838, 1256.953, 1302.838, 1428.269), nrow=3, byrow=TRUE) B - matrix(c(4806.033, 1767.480, 2622.744, 1767.480, 3353.603, 3259.680, 2622.744, 3259.680, 3476.790), nrow=3, byrow=TRUE) I've tried eigen(solve(B)%*%A) (which is probably not the best way to handle this problem numerically) to approximate the generalized ev problem. And the geigen function geigen(A,B,TRUE,TRUE) The eigenvalues are identical upto the printed 9 digits but the eigenvectors appear to be quite different. Maybe this is what Luke meant. Berend They look quite similar to me: ev - eigen(solve(B,A) )$vectors ge - geigen(A, B, TRUE , TRUE) ev / ge$vl [,1] [,2] [,3] [1,] 0.9324603 0.813422 -0.7423694 [2,] 0.9324603 0.813422 -0.7423694 [3,] 0.9324603 0.813422 -0.7423694 ev / ge$vr [,1] [,2] [,3] [1,] 0.9324603 0.813422 -0.7423694 [2,] 0.9324603 0.813422 -0.7423694 [3,] 0.9324603 0.813422 -0.7423694 (and of course, eigenvectors of any sort are only defined up to a constant multiplier) -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd@cbs.dk Priv: pda...@gmail.com __ 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On 21-04-2012, at 20:20, peter dalgaard wrote: The eigenvalues are identical upto the printed 9 digits but the eigenvectors appear to be quite different. Maybe this is what Luke meant. Berend They look quite similar to me: ev - eigen(solve(B,A) )$vectors ge - geigen(A, B, TRUE , TRUE) ev / ge$vl [,1] [,2] [,3] [1,] 0.9324603 0.813422 -0.7423694 [2,] 0.9324603 0.813422 -0.7423694 [3,] 0.9324603 0.813422 -0.7423694 ev / ge$vr [,1] [,2] [,3] [1,] 0.9324603 0.813422 -0.7423694 [2,] 0.9324603 0.813422 -0.7423694 [3,] 0.9324603 0.813422 -0.7423694 (and of course, eigenvectors of any sort are only defined up to a constant multiplier) Correct. I should have checked your way and not optically. 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On 19-04-2012, at 20:50, Jonathan Greenberg wrote: Folks: I'm trying to port some code from python over to R, and I'm running into a wall finding R code that can solve a generalized eigenvalue problem following this function model: http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html Any ideas? I don't want to call python from within R for various reasons, I'd prefer a native R solution if one exists. Thanks! An old thread is here: http://tolstoy.newcastle.edu.au/R/help/05/06/6802.html R does provide eigen(). R has the Lapack routines dggev and dgges in its library. You'd have to write your own .Fortran interface to those routines. 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
So I am a complete noob at doing this type of linking. When I write this statement (all the input values are assigned, but I have to confess I don't know what to do with the output variables -- in this call they are all assigned NA): .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA, BETA=NA, + VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE=base) I get: Error in .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR = NA, : C symbol name La_dgeev not in DLL for package base I'm running this on Windows 7 x64 version of R 2.14.2. The R_ext/Lapack.h file states: /* DGEEV - compute for an N-by-N real nonsymmetric matrix A, the */ /* eigenvalues and, optionally, the left and/or right eigenvectors */ La_extern void F77_NAME(dgeev)(const char* jobvl, const char* jobvr, const int* n, double* a, const int* lda, double* wr, double* wi, double* vl, const int* ldvl, double* vr, const int* ldvr, double* work, const int* lwork, int* info); Any ideas? --j On Fri, Apr 20, 2012 at 1:37 AM, Berend Hasselman b...@xs4all.nl wrote: On 19-04-2012, at 20:50, Jonathan Greenberg wrote: Folks: I'm trying to port some code from python over to R, and I'm running into a wall finding R code that can solve a generalized eigenvalue problem following this function model: http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html Any ideas? I don't want to call python from within R for various reasons, I'd prefer a native R solution if one exists. Thanks! An old thread is here: http://tolstoy.newcastle.edu.au/R/help/05/06/6802.html R does provide eigen(). R has the Lapack routines dggev and dgges in its library. You'd have to write your own .Fortran interface to those routines. Berend -- Jonathan A. Greenberg, PhD Assistant Professor Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 607 South Mathews Avenue, MC 150 Urbana, IL 61801 Phone: 415-763-5476 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html [[alternative HTML version deleted]] __ 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
Incidentally, if you want to test this out: A [,1] [,2] [,3] [1,] 1457.738 1053.181 1256.953 [2,] 1053.181 1213.728 1302.838 [3,] 1256.953 1302.838 1428.269 B [,1] [,2] [,3] [1,] 4806.033 1767.480 2622.744 [2,] 1767.480 3353.603 3259.680 [3,] 2622.744 3259.680 3476.790 I THINK this is correct for the other parameters: JOBVL=N JOBVR=N N=dim(A)[[1]] LDA=N LDB=N LDVL=N LDVR=N LWORK=max(1,8*N) .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA, BETA=NA, VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE=base) LAPACK's documentation is: http://www.netlib.org/lapack/double/dggev.f --j On Fri, Apr 20, 2012 at 11:58 AM, Jonathan Greenberg j...@illinois.eduwrote: So I am a complete noob at doing this type of linking. When I write this statement (all the input values are assigned, but I have to confess I don't know what to do with the output variables -- in this call they are all assigned NA): .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA, BETA=NA, + VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE=base) I get: Error in .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR = NA, : C symbol name La_dgeev not in DLL for package base I'm running this on Windows 7 x64 version of R 2.14.2. The R_ext/Lapack.h file states: /* DGEEV - compute for an N-by-N real nonsymmetric matrix A, the */ /* eigenvalues and, optionally, the left and/or right eigenvectors */ La_extern void F77_NAME(dgeev)(const char* jobvl, const char* jobvr, const int* n, double* a, const int* lda, double* wr, double* wi, double* vl, const int* ldvl, double* vr, const int* ldvr, double* work, const int* lwork, int* info); Any ideas? --j On Fri, Apr 20, 2012 at 1:37 AM, Berend Hasselman b...@xs4all.nl wrote: On 19-04-2012, at 20:50, Jonathan Greenberg wrote: Folks: I'm trying to port some code from python over to R, and I'm running into a wall finding R code that can solve a generalized eigenvalue problem following this function model: http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html Any ideas? I don't want to call python from within R for various reasons, I'd prefer a native R solution if one exists. Thanks! An old thread is here: http://tolstoy.newcastle.edu.au/R/help/05/06/6802.html R does provide eigen(). R has the Lapack routines dggev and dgges in its library. You'd have to write your own .Fortran interface to those routines. Berend -- Jonathan A. Greenberg, PhD Assistant Professor Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 607 South Mathews Avenue, MC 150 Urbana, IL 61801 Phone: 415-763-5476 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html -- Jonathan A. Greenberg, PhD Assistant Professor Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 607 South Mathews Avenue, MC 150 Urbana, IL 61801 Phone: 415-763-5476 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html [[alternative HTML version deleted]] __ 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
Ok, I figured out a solution and I'd like to get some feedback on this from
the R-helpers as to how I could modify the following to be package
friendly -- the main thing I'm worried about is how to dynamically set the
dyn.load statement below correctly (obviously, its hard coded to my
particular install, and would only work with windows since I'm using a
.dll):
Rdggev - function(JOBVL=F,JOBVR=T,A,B)
{
# R implementation of the DGGEV LAPACK function (with generalized
eigenvalue computation)
# See http://www.netlib.org/lapack/double/dggev.f
# coded by Jonathan A. Greenberg j...@illinois.edu
dyn.load(C:\\Program Files\\R\\R-2.14.2\\bin\\x64\\Rlapack.dll)
if(JOBVL)
{
JOBVL=V
} else
{
JOBVL=N
}
if(JOBVR)
{
JOBVR=V
} else
{
JOBVR=N
}
# Note, no error checking is performed.
# Input parameters
N=dim(A)[[1]]
LDA=N
LDB=N
LDVL=N
LDVR=N
LWORK=as.integer(max(1,8*N))
Rdggev_out - .Fortran(dggev, JOBVL, JOBVR, N, A, LDA, B, LDB,
double(N), double(N), double(N),
array(data=0,dim=c(LDVL,N)), LDVL, array(data=0,dim=c(LDVR,N)), LDVR,
double(max(1,LWORK)), LWORK, integer(1))
names(Rdggev_out)=c(JOBVL,JOBVR,N,A,LDA,B,LDB,ALPHAR,ALPHAI,
BETA,VL,LDVL,VR,LDVR,WORK,LWORK,INFO)
Rdggev_out$GENEIGENVALUES=(Rdggev_out$ALPHAR+Rdggev_out$ALPHAI)/Rdggev_out$BETA
return(Rdggev_out)
}
On Fri, Apr 20, 2012 at 12:01 PM, Jonathan Greenberg j...@illinois.eduwrote:
Incidentally, if you want to test this out:
A
[,1] [,2] [,3]
[1,] 1457.738 1053.181 1256.953
[2,] 1053.181 1213.728 1302.838
[3,] 1256.953 1302.838 1428.269
B
[,1] [,2] [,3]
[1,] 4806.033 1767.480 2622.744
[2,] 1767.480 3353.603 3259.680
[3,] 2622.744 3259.680 3476.790
I THINK this is correct for the other parameters:
JOBVL=N
JOBVR=N
N=dim(A)[[1]]
LDA=N
LDB=N
LDVL=N
LDVR=N
LWORK=max(1,8*N)
.Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA,
BETA=NA,
VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE=base)
LAPACK's documentation is:
http://www.netlib.org/lapack/double/dggev.f
--j
On Fri, Apr 20, 2012 at 11:58 AM, Jonathan Greenberg j...@illinois.eduwrote:
So I am a complete noob at doing this type of linking. When I write this
statement (all the input values are assigned, but I have to confess I don't
know what to do with the output variables -- in this call they are all
assigned NA):
.Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA,
BETA=NA,
+ VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE=base)
I get:
Error in .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR = NA,
:
C symbol name La_dgeev not in DLL for package base
I'm running this on Windows 7 x64 version of R 2.14.2. The
R_ext/Lapack.h file states:
/* DGEEV - compute for an N-by-N real nonsymmetric matrix A, the */
/* eigenvalues and, optionally, the left and/or right eigenvectors */
La_extern void
F77_NAME(dgeev)(const char* jobvl, const char* jobvr,
const int* n, double* a, const int* lda,
double* wr, double* wi, double* vl, const int* ldvl,
double* vr, const int* ldvr,
double* work, const int* lwork, int* info);
Any ideas?
--j
On Fri, Apr 20, 2012 at 1:37 AM, Berend Hasselman b...@xs4all.nl wrote:
On 19-04-2012, at 20:50, Jonathan Greenberg wrote:
Folks:
I'm trying to port some code from python over to R, and I'm running
into a
wall finding R code that can solve a generalized eigenvalue problem
following this function model:
http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html
Any ideas? I don't want to call python from within R for various
reasons,
I'd prefer a native R solution if one exists. Thanks!
An old thread is here:
http://tolstoy.newcastle.edu.au/R/help/05/06/6802.html
R does provide eigen().
R has the Lapack routines dggev and dgges in its library.
You'd have to write your own .Fortran interface to those routines.
Berend
--
Jonathan A. Greenberg, PhD
Assistant Professor
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 415-763-5476
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
http://www.geog.illinois.edu/people/JonathanGreenberg.html
--
Jonathan A. Greenberg, PhD
Assistant Professor
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 415-763-5476
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
http://www.geog.illinois.edu/people/JonathanGreenberg.html
--
Jonathan A. Greenberg, PhD
Assistant Professor
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South Mathews Avenue, MC 150
Urbana, IL 61801
Phone: 415-763-5476
AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On 20-04-2012, at 18:58, Jonathan Greenberg wrote: So I am a complete noob at doing this type of linking. When I write this statement (all the input values are assigned, but I have to confess I don't know what to do with the output variables -- in this call they are all assigned NA): .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA, BETA=NA, + VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE=base) I get: Error in .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR = NA, : C symbol name La_dgeev not in DLL for package base I'm running this on Windows 7 x64 version of R 2.14.2. The R_ext/Lapack.h file states: /* DGEEV - compute for an N-by-N real nonsymmetric matrix A, the */ /* eigenvalues and, optionally, the left and/or right eigenvectors */ La_extern void F77_NAME(dgeev)(const char* jobvl, const char* jobvr, const int* n, double* a, const int* lda, double* wr, double* wi, double* vl, const int* ldvl, double* vr, const int* ldvr, double* work, const int* lwork, int* info); Any ideas? Yes. Stop this. You are completely on the wrong track and getting yourself into a real muddle. 1. you are mixing up dgeev and dggev. You are using the dggev arguments to call dgeev. Even if you got it all correct this would result in nasty errors. 2. La_dgeev doesn't exist in the R sources (nor does it exist in Lapack) and you don't need to use dgeev. See below. Since you are a noob I would advise you to not try and interface with Fortran or C at this moment. It is possible to interface with dggev for generalized eigenvalues but it is not a straightforward and easy job to do. I've had a look at the R and C code used in the R sources to interface with dgeev. It would need a lot of TLC. As I told you before: R provides a function eigen for calculating normal eigenvalues end eigenvectors. It's probably a good idea to explain why you need generalized eigenvalues. Then someone may be able to come up with suggestions and/or alternatives. 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
My apologies: When reviewing my initial email I made a typo -- what I needed was dggev, not dgeev. I have confirmed that my function reproduces the scipy outputs I mentioned in my first email. The function is part of lapack. I'm not an R noob, I just haven't dealt with external Fortran or C calls before. The reason I need it is that I am trying to port an existing python/IDL algorithm to R (iMad/RADCAL, located at http://mcanty.homepage.t-online.de/software.html). The algorithm is fairly complex so I'm trying to reproduce it line-by-line as closely as possible before going back and improving its efficiency using various R functions (checking the inputs and the outputs between R and the python code each step). This algorithm does a weighted canonical correlation analysis which does not apparently exist in other CCA R packages. So, my main question at this point is given I can call dggev via .Fortran() and get the outputs I want, what is the best way to package this up with a dyn.load() statement that would work on any R install that has the lapack libraries available to it? --j Yes. Stop this. You are completely on the wrong track and getting yourself into a real muddle. 1. you are mixing up dgeev and dggev. You are using the dggev arguments to call dgeev. Even if you got it all correct this would result in nasty errors. 2. La_dgeev doesn't exist in the R sources (nor does it exist in Lapack) and you don't need to use dgeev. See below. Since you are a noob I would advise you to not try and interface with Fortran or C at this moment. It is possible to interface with dggev for generalized eigenvalues but it is not a straightforward and easy job to do. I've had a look at the R and C code used in the R sources to interface with dgeev. It would need a lot of TLC. As I told you before: R provides a function eigen for calculating normal eigenvalues end eigenvectors. It's probably a good idea to explain why you need generalized eigenvalues. Then someone may be able to come up with suggestions and/or alternatives. Berend -- Jonathan A. Greenberg, PhD Assistant Professor Department of Geography and Geographic Information Science University of Illinois at Urbana-Champaign 607 South Mathews Avenue, MC 150 Urbana, IL 61801 Phone: 415-763-5476 AIM: jgrn307, MSN: jgrn...@hotmail.com, Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html [[alternative HTML version deleted]] __ 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.
Re: [R] Solve an ordinary or generalized eigenvalue problem in R?
On Fri, Apr 20, 2012 at 2:18 PM, Jonathan Greenberg j...@illinois.edu wrote:
Ok, I figured out a solution and I'd like to get some feedback on this from
the R-helpers as to how I could modify the following to be package
friendly -- the main thing I'm worried about is how to dynamically set the
There's documentation on this in the Writing R Extensions manual,
See 1.2.1 Using Makevars and Section 6.7 Numerical analysis
subrouties. They recommend looking at the package fastICA. It appears
to me you will use some platform neutral Makefile variables.
If you give a full working example from windows--your function, data
to use it--I'll see what I can do on the Linux side. I've never had
need to worry about this question before, but it is about time I
learned.
Before you think about packaging something like this, I'd say step one
is to clean up your code so it is easier to read. R CMD check won't
let you get bye with T or F in place of TRUE and FALSE, and you seem
to mix - and = in your code, which makes it difficult to read.
pj
dyn.load statement below correctly (obviously, its hard coded to my
particular install, and would only work with windows since I'm using a
.dll):
Rdggev - function(JOBVL=F,JOBVR=T,A,B)
{
# R implementation of the DGGEV LAPACK function (with generalized
eigenvalue computation)
# See http://www.netlib.org/lapack/double/dggev.f
# coded by Jonathan A. Greenberg j...@illinois.edu
dyn.load(C:\\Program Files\\R\\R-2.14.2\\bin\\x64\\Rlapack.dll)
if(JOBVL)
{
JOBVL=V
} else
{
JOBVL=N
}
if(JOBVR)
{
JOBVR=V
} else
{
JOBVR=N
}
# Note, no error checking is performed.
# Input parameters
N=dim(A)[[1]]
LDA=N
LDB=N
LDVL=N
LDVR=N
LWORK=as.integer(max(1,8*N))
Rdggev_out - .Fortran(dggev, JOBVL, JOBVR, N, A, LDA, B, LDB,
double(N), double(N), double(N),
array(data=0,dim=c(LDVL,N)), LDVL, array(data=0,dim=c(LDVR,N)), LDVR,
double(max(1,LWORK)), LWORK, integer(1))
names(Rdggev_out)=c(JOBVL,JOBVR,N,A,LDA,B,LDB,ALPHAR,ALPHAI,
BETA,VL,LDVL,VR,LDVR,WORK,LWORK,INFO)
Rdggev_out$GENEIGENVALUES=(Rdggev_out$ALPHAR+Rdggev_out$ALPHAI)/Rdggev_out$BETA
return(Rdggev_out)
}
On Fri, Apr 20, 2012 at 12:01 PM, Jonathan Greenberg j...@illinois.eduwrote:
Incidentally, if you want to test this out:
A
[,1] [,2] [,3]
[1,] 1457.738 1053.181 1256.953
[2,] 1053.181 1213.728 1302.838
[3,] 1256.953 1302.838 1428.269
B
[,1] [,2] [,3]
[1,] 4806.033 1767.480 2622.744
[2,] 1767.480 3353.603 3259.680
[3,] 2622.744 3259.680 3476.790
I THINK this is correct for the other parameters:
JOBVL=N
JOBVR=N
N=dim(A)[[1]]
LDA=N
LDB=N
LDVL=N
LDVR=N
LWORK=max(1,8*N)
.Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA,
BETA=NA,
VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE=base)
LAPACK's documentation is:
http://www.netlib.org/lapack/double/dggev.f
--j
On Fri, Apr 20, 2012 at 11:58 AM, Jonathan Greenberg
j...@illinois.eduwrote:
So I am a complete noob at doing this type of linking. When I write this
statement (all the input values are assigned, but I have to confess I don't
know what to do with the output variables -- in this call they are all
assigned NA):
.Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR=NA, ALPHAI=NA,
BETA=NA,
+ VL=NA, LDVL, VR=NA, LDVR, WORK=NA, LWORK, INFO=NA, PACKAGE=base)
I get:
Error in .Call(La_dgeev, JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR = NA,
:
C symbol name La_dgeev not in DLL for package base
I'm running this on Windows 7 x64 version of R 2.14.2. The
R_ext/Lapack.h file states:
/* DGEEV - compute for an N-by-N real nonsymmetric matrix A, the */
/* eigenvalues and, optionally, the left and/or right eigenvectors */
La_extern void
F77_NAME(dgeev)(const char* jobvl, const char* jobvr,
const int* n, double* a, const int* lda,
double* wr, double* wi, double* vl, const int* ldvl,
double* vr, const int* ldvr,
double* work, const int* lwork, int* info);
Any ideas?
--j
On Fri, Apr 20, 2012 at 1:37 AM, Berend Hasselman b...@xs4all.nl wrote:
On 19-04-2012, at 20:50, Jonathan Greenberg wrote:
Folks:
I'm trying to port some code from python over to R, and I'm running
into a
wall finding R code that can solve a generalized eigenvalue problem
following this function model:
http://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.eig.html
Any ideas? I don't want to call python from within R for various
reasons,
I'd prefer a native R solution if one exists. Thanks!
An old thread is here:
http://tolstoy.newcastle.edu.au/R/help/05/06/6802.html
R does provide eigen().
R has the Lapack routines dggev and dgges in its library.
You'd have to write your own .Fortran interface to those routines.
Berend
--
Jonathan A. Greenberg, PhD
Assistant Professor
Department of Geography and Geographic Information Science
University of Illinois at Urbana-Champaign
607 South
