Solve an ordinary or generalized eigenvalue problem in R?

[Jonathan Greenberg-4]

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!

--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: [phone]
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Berend Hasselman]

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

______________________________________________
[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 Greenberg-4]

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: [phone]
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Jonathan Greenberg-4]

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: [phone]
> 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: [phone]
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Jonathan Greenberg-4]

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: [phone]
>> 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: [phone]
> 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: [phone]
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Berend Hasselman]

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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Jonathan Greenberg-4]

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: [phone]
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[PaulJohnson32gmail]

On Fri, Apr 20, 2012 at 2:18 PM, Jonathan Greenberg <[hidden email]> 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 <[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: [phone]
>>> 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: [phone]
>> 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: [phone]
> 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-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.

--
Paul E. Johnson
Professor, Political Science    Assoc. Director
1541 Lilac Lane, Room 504     Center for Research Methods
University of Kansas               University of Kansas
http://pj.freefaculty.org            http://quant.ku.edu

______________________________________________
[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.

[Berend Hasselman]

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:

<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

______________________________________________
[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.

[Luke Hartigan]

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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Berend Hasselman]

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

______________________________________________
[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.

[Peter Dalgaard-2]

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: [phone]
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[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.

[Rui Barradas]

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

______________________________________________
[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.

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

[Berend Hasselman]

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

______________________________________________
[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.

[Berend Hasselman]

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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Peter Dalgaard-2]

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: [phone]
Email: [hidden email]  Priv: [hidden email]

______________________________________________
[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.

[Berend Hasselman]

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 Greenberg-4]

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 <[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: [phone]
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Berend Hasselman]

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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.

[Jonathan Greenberg-4]

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 <[hidden email]>
# 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 <[hidden email]> 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: [phone]
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-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.