# Solve an ordinary or generalized eigenvalue problem in R?

21 messages
12
Open this post in threaded view
|
Report Content as Inappropriate

## Solve an ordinary or generalized eigenvalue problem in R?

 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.htmlAny 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! --j -- 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: [hidden email], Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html        [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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.htmlR 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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 <[hidden email]> 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: [hidden email], Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html        [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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 <[hidden email]>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? > > --j > > > > On Fri, Apr 20, 2012 at 1:37 AM, Berend Hasselman <[hidden email]> 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: [hidden email], 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: [hidden email], Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html        [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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 <[hidden email]>  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 <[hidden email]>wrote: > 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 <[hidden email]>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? >> >> --j >> >> >> >> On Fri, Apr 20, 2012 at 1:37 AM, Berend Hasselman <[hidden email]> 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: [hidden email], 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: [hidden email], 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: [hidden email], Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html        [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

 In reply to this post by Jonathan Greenberg-4 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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: [hidden email], Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html        [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

 In reply to this post by Jonathan Greenberg-4 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 <[hidden email]> > 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: # 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) Berend ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

 In reply to this post by Jonathan Greenberg-4 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]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

 In reply to this post by Berend Hasselman On 21-04-2012, at 11:40, Berend Hasselman wrote: > > ..... > See this: > > > # 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

 In reply to this post by Jonathan Greenberg-4 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: [hidden email]  Priv: [hidden email] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

Open this post in threaded view
|
Report Content as Inappropriate

## Re: Solve an ordinary or generalized eigenvalue problem in R?

 In reply to this post by Peter Dalgaard-2 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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: [hidden email]  Priv: [hidden email] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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.htmlSome of these are already implemented, but clearly not all of them. --j On Sat, Apr 21, 2012 at 1:31 PM, Berend Hasselman <[hidden email]> 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 > > ______________________________________________ > [hidden email] 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: [hidden email], Gchat: jgrn307, Skype: jgrn3007 http://www.geog.illinois.edu/people/JonathanGreenberg.html        [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|
Report Content as Inappropriate

## Re: 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.