Skip to content

matrix exponential: M0

2 messages · Vicente Canto Casasola, Martin Maechler

Original
Vicente Canto Casasola
# 
H i, all!

First of all, I'd like to apologize for my poor English. It's for years 
I don't use it.

This is a R-version of a function I wrote a long ago for my HP48 calculator.
It works with the binary expression of the power and just need to 
duplicate the mem used by X.


Hope this helps.

mtx.exp<-function(X,n)
#Function to calculate the n-th power of a matrix X;
{
phi <- diag(rep(1,length(X[1,])))
pot <- X #This is the first power of the matrix.

while (n > 0)
  {
  if (n%%2)
    {
    phi <- phi%*%pot;
    }
    n <- n%/%2;
    pot <- pot %*% pot;
  }
return(phi);
}



#Here is some output:

 > xx <- matrix(c(1,0,1,1),2,2)
 > xx
     [,1] [,2]
[1,]    1    1
[2,]    0    1
 > mtx.exp(xx,3)
     [,1] [,2]
[1,]    1    3
[2,]    0    1
 > mtx.exp(xx,4)
     [,1] [,2]
[1,]    1    4
[2,]    0    1
 > mtx.exp(xx,6)
     [,1] [,2]
[1,]    1    6
[2,]    0    1
 > mtx.exp(xx,10)
     [,1] [,2]
[1,]    1   10
[2,]    0    1
 > mtx.exp(xx,1000)
     [,1] [,2]
[1,]    1 1000
[2,]    0    1


Vicente D. Canto Casasola
Martin Maechler
# 
Vicente> H i, all!  First of all, I'd like to apologize for
    Vicente> my poor English. It's for years I don't use it.

no problem at all.

    Vicente> This is a R-version of a function I wrote a long
    Vicente> ago for my HP48 calculator.  It works with the
    Vicente> binary expression of the power and just need to
    Vicente> duplicate the mem used by X.

excellent. This is really the way to solve the problem 
I think. 

As I've mentioned earlier in this thread,
computing a matrix "power" is really much easier than the
matrix exponential.

Hence I wouldn't use exponential in the function name.
Also note that trailing ";" are considered as `dirty' (they are
completely superfluous).

These slight modifications (+ initial "test") 
give

matPower <- function(X,n)
## Function to calculate the n-th power of a matrix X
{
    if(n != round(n)) {
        n <- round(n)
        warning("rounding exponent `n' to", n)
    }
    phi <- diag(nrow = nrow(X))
    pot <- X # the first power of the matrix.

    while (n > 0)
    {
        if (n %% 2)
            phi <- phi %*% pot

        n <- n %/% 2
        pot <- pot %*% pot
    }
    return(phi)
}


Regards,

Martin Maechler <maechler at stat.math.ethz.ch>	http://stat.ethz.ch/~maechler/
Seminar fuer Statistik, ETH-Zentrum  LEO C16	Leonhardstr. 27
ETH (Federal Inst. Technology)	8092 Zurich	SWITZERLAND
phone: x-41-1-632-3408		fax: ...-1228			<><