power of a matrix

Thank you very much! Thanks also to the authors of this function,
Vincente Canto Cassola and Martin Maechler!

This is exactly what I hoped for.

Best,
Roland
-----Original Message-----
From: Dimitris Rizopoulos 
[mailto:dimitris.rizopoulos at med.kuleuven.be] 
Sent: Wednesday, August 17, 2005 1:38 PM
To: Rau, Roland
Cc: r-help at stat.math.ethz.ch
Subject: Re: [R] power of a matrix

look at function ?mtx.exp() in the Malmig package, e.g.,

A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)
n.zero <- c(1,0,0)
##################
library(Malmig)

all.equal(A %*% (A %*% n.zero), mtx.exp(A, 2) %*% n.zero)
all.equal(A %*% (A %*% (A %*% n.zero)), mtx.exp(A, 3) %*% n.zero)

mtx.exp(A, 15) %*% n.zero

I hope it helps.

Best,
Dimitris

----
Dimitris Rizopoulos
Ph.D. Student
Biostatistical Centre
School of Public Health
Catholic University of Leuven

Address: Kapucijnenvoer 35, Leuven, Belgium
Tel: +32/16/336899
Fax: +32/16/337015
Web: http://www.med.kuleuven.be/biostat/
     http://www.student.kuleuven.be/~m0390867/dimitris.htm

----- Original Message ----- 
From: "Rau, Roland" <Rau at demogr.mpg.de>
To: "R-Help" <r-help at stat.math.ethz.ch>
Sent: Wednesday, August 17, 2005 12:15 PM
Subject: [R] power of a matrix

Dear all,

I have a population with three age-classes, at time t=0 the 
population
is:
n.zero <- c(1,0,0)
I have a transition matrix A which denotes "fertility" and 
"survival":

A <- matrix(c(0,1,5, 0.3,0,0, 0,0.5,0), ncol=3, byrow=TRUE)

To obtain the population at t=1, I calculate:
A %*% n.zero

To obtain the population t=2, I calculate:
A %*% (A %*% n.zero)

... and so on ...

I thought now to obtain the population at time x, I should simply 
do:
A^x %*% n.zero

But this, of course, does not work since the following two 
statements
are not equivalent for matrices:
A %*% A
A * A

Is there something like a "powermatrix"-function?

Thanks,
Roland

P.S. The example is taken from:
Example 3.1 ("A linear, time-invariant model") from Keyfitz/Caswell:
"Applied Mathematical Demography", 3rd Edition, 2005, p. 52f

+++++
This mail has been sent through the MPI for Demographic 
Rese...{{dropped}}

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! 
http://www.R-project.org/posting-guide.html

+++++
This mail has been sent through the MPI for Demographic Rese...{{dropped}}
"Rau, Roland" <Rau at demogr.mpg.de> writes:
Thank you very much! Thanks also to the authors of this function,
Vincente Canto Cassola and Martin Maechler!

This is exactly what I hoped for.
....
look at function ?mtx.exp() in the Malmig package, e.g.,
Also, there is mexp() in the Matrix package. I'm not sure about the
relative merits. mexp() is one of the less dubious implementations of
matrix exponentials, but it does require to and from class "Matrix".
mtx.exp is a bit unfortunately named as it appears to calculate matrix
*powers* (which in this case is what you need).
O__  ---- Peter Dalgaard             ??ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907
Hi, Peter:

	  I couldn't find "mexp" in the Matrix package, but I did find it in 
fMultivar and in Lindsey's rmutil.  These are different functions, but 
produced essentially the same answer for mexp(array(1:4, dim=c(2,2))). 
While hunting for that, I also also found reference by Doug Bates in a 
previous interchange on r-help to "a classic paper ... I would recommend 
reading":

  Moler C., van Loan C., (2003); _Nineteen dubious ways to compute
      the exponential of a matrix,  twenty-five years later_, SIAM
      Review 45, 3-49.

	  This paper was cited in the help page for mexp in fMultivar but not 
in rmutil.

	  spencer graves

"Rau, Roland" <Rau at demogr.mpg.de> writes:

Thank you very much! Thanks also to the authors of this function,
Vincente Canto Cassola and Martin Maechler!

This is exactly what I hoped for.
....

look at function ?mtx.exp() in the Malmig package, e.g.,

Also, there is mexp() in the Matrix package. I'm not sure about the
relative merits. mexp() is one of the less dubious implementations of
matrix exponentials, but it does require to and from class "Matrix".
mtx.exp is a bit unfortunately named as it appears to calculate matrix
*powers* (which in this case is what you need).

Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915
Its expm.
Hi, Peter:

         I couldn't find "mexp" in the Matrix package, but I did find it in
fMultivar and in Lindsey's rmutil.  These are different functions, but
produced essentially the same answer for mexp(array(1:4, dim=c(2,2))).
While hunting for that, I also also found reference by Doug Bates in a
previous interchange on r-help to "a classic paper ... I would recommend
reading":

 Moler C., van Loan C., (2003); _Nineteen dubious ways to compute
     the exponential of a matrix,  twenty-five years later_, SIAM
     Review 45, 3-49.

         This paper was cited in the help page for mexp in fMultivar but not
in rmutil.

         spencer graves

Peter Dalgaard wrote:

"Rau, Roland" <Rau at demogr.mpg.de> writes:

Thank you very much! Thanks also to the authors of this function,
Vincente Canto Cassola and Martin Maechler!

This is exactly what I hoped for.
....

look at function ?mtx.exp() in the Malmig package, e.g.,

Also, there is mexp() in the Matrix package. I'm not sure about the
relative merits. mexp() is one of the less dubious implementations of
matrix exponentials, but it does require to and from class "Matrix".
mtx.exp is a bit unfortunately named as it appears to calculate matrix
*powers* (which in this case is what you need).

--
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Gabor Grothendieck <ggrothendieck at gmail.com> writes:
Its expm.
Doh. I had even looked up its help page...
On 8/17/05, Spencer Graves <spencer.graves at pdf.com> wrote:
Hi, Peter:

         I couldn't find "mexp" in the Matrix package, but I did find it in
fMultivar and in Lindsey's rmutil.  These are different functions, but
produced essentially the same answer for mexp(array(1:4, dim=c(2,2))).
While hunting for that, I also also found reference by Doug Bates in a
previous interchange on r-help to "a classic paper ... I would recommend
reading":

 Moler C., van Loan C., (2003); _Nineteen dubious ways to compute
     the exponential of a matrix,  twenty-five years later_, SIAM
     Review 45, 3-49.

         This paper was cited in the help page for mexp in fMultivar but not
in rmutil.
As far as I remember, that's because rmutil is using one of the more
dubious methods, namely eigendecomposition. It breaks down if there
are eigenvalues with multiplicity greater than 1, which does happen
(systematically) in some statistical models.
         spencer graves

Peter Dalgaard wrote:

"Rau, Roland" <Rau at demogr.mpg.de> writes:

Thank you very much! Thanks also to the authors of this function,
Vincente Canto Cassola and Martin Maechler!

This is exactly what I hoped for.
....

look at function ?mtx.exp() in the Malmig package, e.g.,

Also, there is mexp() in the Matrix package. I'm not sure about the
relative merits. mexp() is one of the less dubious implementations of
matrix exponentials, but it does require to and from class "Matrix".
mtx.exp is a bit unfortunately named as it appears to calculate matrix
*powers* (which in this case is what you need).

--
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

O__  ---- Peter Dalgaard             ??ster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark          Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)                  FAX: (+45) 35327907
Hi, Gabor:  Thanks.  spencer graves

Its expm.

On 8/17/05, Spencer Graves <spencer.graves at pdf.com> wrote:

Hi, Peter:

        I couldn't find "mexp" in the Matrix package, but I did find it in
fMultivar and in Lindsey's rmutil.  These are different functions, but
produced essentially the same answer for mexp(array(1:4, dim=c(2,2))).
While hunting for that, I also also found reference by Doug Bates in a
previous interchange on r-help to "a classic paper ... I would recommend
reading":

Moler C., van Loan C., (2003); _Nineteen dubious ways to compute
    the exponential of a matrix,  twenty-five years later_, SIAM
    Review 45, 3-49.

        This paper was cited in the help page for mexp in fMultivar but not
in rmutil.

        spencer graves

Peter Dalgaard wrote:

"Rau, Roland" <Rau at demogr.mpg.de> writes:

Thank you very much! Thanks also to the authors of this function,
Vincente Canto Cassola and Martin Maechler!

This is exactly what I hoped for.
....

look at function ?mtx.exp() in the Malmig package, e.g.,

Also, there is mexp() in the Matrix package. I'm not sure about the
relative merits. mexp() is one of the less dubious implementations of
matrix exponentials, but it does require to and from class "Matrix".
mtx.exp is a bit unfortunately named as it appears to calculate matrix
*powers* (which in this case is what you need).

--
Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915

______________________________________________
R-help at stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html

Spencer Graves, PhD
Senior Development Engineer
PDF Solutions, Inc.
333 West San Carlos Street Suite 700
San Jose, CA 95110, USA

spencer.graves at pdf.com
www.pdf.com <http://www.pdf.com>
Tel:  408-938-4420
Fax: 408-280-7915