# [R] A function for raising a matrix to a power?

Atte Tenkanen attenka at utu.fi
Wed May 9 16:07:15 CEST 2007

```Hello,

Thanks for many replys. I tested all the functions presented.
I'm a beginner in linear algebra, but now I have a serious question ;-)
Here is a matrix A, which determinant is 3, so it is nonsingular.
Then there are similar computer runs done with each function proposed.
I have calculated powers of A between A^274 - A^277.
In the first version (see the outputs lower), when n=277,
there comes "zero in the upper right corner".
Can you say, if the first function is more reliable than others?
What can you say about the accuracy of calculations?

-Atte

#-------------------------------------------------------------#

# Matrix A:
A=rbind(c(-3,-4,-2),c(3,5,1),c(2,-1,4))

A
det(A)

# 1st version:
"%^%"<-function(A,n){
if(n==1) A else {B<-A; for(i in (2:n)){A<-A%*%B}}; A
}

for(i in 274:277){print(i);print(A%^%i)}

# 2nd version:
"%^%" <- function(A, n) if(n == 1) A else A %*% (A %^% (n-1))

for(i in 274:277){print(i);print(A%^%i)}

# 3rd version:

mp <- function(mat,pow){
ans <- mat
for ( i in 1:(pow-1)){
ans <- mat%*%ans
}
return(ans)
}

for(i in 274:277){print(i);print(mp(A,i))}

# 4th version:

library(Malmig)

mtx.exp
for(i in 274:277){print(i);print(mtx.exp(A,i))}

# 5th version:

matrix.power <- function(mat, n)
{
# test if mat is a square matrix
# treat n < 0 and n = 0 -- this is left as an exercise
# trap non-integer n and return an error
if (n == 1) return(mat)
result <- diag(1, ncol(mat))
while (n > 0) {
if (n %% 2 != 0) {
result <- result %*% mat
n <- n - 1
}
mat <- mat %*% mat
n <- n / 2
}
return(result)
}

for(i in 274:277){print(i);print(matrix.power(A,i))}

# 6th version:

expM.sd <- function(X,e){Xsd <- eigen(X); Xsd\$vec %*% diag(Xsd\$val^e) %*%t(Xsd\$vec)}

for(i in 274:277){print(i);print(expM.sd(A,i))}

#-----------------OUTPUTS---------------------------#

> A=rbind(c(-3,-4,-2),c(3,5,1),c(2,-1,4))
>
> A
[,1] [,2] [,3]
[1,]   -3   -4   -2
[2,]    3    5    1
[3,]    2   -1    4
> det(A)
[1] 3
>
> # 1st version:
> "%^%"<-function(A,n){
+   if(n==1) A else {B<-A; for(i in (2:n)){A<-A%*%B}}; A
+   }
>
> for(i in 274:277){print(i);print(A%^%i)}
[1] 274
[,1]           [,2]           [,3]
[1,]  1.615642e+131  3.231283e+131 -3.940201e+115
[2,] -5.385472e+130 -1.077094e+131  4.925251e+114
[3,] -3.769831e+131 -7.539661e+131  7.880401e+115
[1] 275
[,1]           [,2]           [,3]
[1,]  4.846925e+131  9.693850e+131 -1.182060e+116
[2,] -1.615642e+131 -3.231283e+131   0.000000e+00
[3,] -1.130949e+132 -2.261898e+132  4.728241e+116
[1] 276
[,1]           [,2]           [,3]
[1,]  1.454078e+132  2.908155e+132 -1.576080e+116
[2,] -4.846925e+131 -9.693850e+131   0.000000e+00
[3,] -3.392848e+132 -6.785695e+132  1.576080e+117
[1] 277
[,1]           [,2]           [,3]
[1,]  4.362233e+132  8.724465e+132   0.000000e+00
[2,] -1.454078e+132 -2.908155e+132 -1.576080e+116
[3,] -1.017854e+133 -2.035709e+133  3.782593e+117
>
>
> # 2nd version:
> "%^%" <- function(A, n) if(n == 1) A else A %*% (A %^% (n-1))
>
> for(i in 274:277){print(i);print(A%^%i)}
[1] 274
[,1]           [,2]           [,3]
[1,]  1.615642e+131  3.231283e+131 -3.101125e+114
[2,] -5.385472e+130 -1.077094e+131  1.533306e+114
[3,] -3.769831e+131 -7.539661e+131  5.664274e+114
[1] 275
[,1]           [,2]           [,3]
[1,]  4.846925e+131  9.693850e+131 -8.158398e+114
[2,] -1.615642e+131 -3.231283e+131  4.027430e+114
[3,] -1.130949e+132 -2.261898e+132  1.492154e+115
[1] 276
[,1]           [,2]           [,3]
[1,]  1.454078e+132  2.908155e+132 -2.147761e+115
[2,] -4.846925e+131 -9.693850e+131  1.058350e+115
[3,] -3.392848e+132 -6.785695e+132  3.934193e+115
[1] 277
[,1]           [,2]           [,3]
[1,]  4.362233e+132  8.724465e+132 -5.658504e+115
[2,] -1.454078e+132 -2.908155e+132  2.782660e+115
[3,] -1.017854e+133 -2.035709e+133  1.038290e+116
>
>
> # 3rd version:
>
> mp <- function(mat,pow){
+ ans <- mat
+ for ( i in 1:(pow-1)){
+ ans <- mat%*%ans
+ }
+ return(ans)
+ }
>
> for(i in 274:277){print(i);print(mp(A,i))}
[1] 274
[,1]           [,2]           [,3]
[1,]  1.615642e+131  3.231283e+131 -3.101125e+114
[2,] -5.385472e+130 -1.077094e+131  1.533306e+114
[3,] -3.769831e+131 -7.539661e+131  5.664274e+114
[1] 275
[,1]           [,2]           [,3]
[1,]  4.846925e+131  9.693850e+131 -8.158398e+114
[2,] -1.615642e+131 -3.231283e+131  4.027430e+114
[3,] -1.130949e+132 -2.261898e+132  1.492154e+115
[1] 276
[,1]           [,2]           [,3]
[1,]  1.454078e+132  2.908155e+132 -2.147761e+115
[2,] -4.846925e+131 -9.693850e+131  1.058350e+115
[3,] -3.392848e+132 -6.785695e+132  3.934193e+115
[1] 277
[,1]           [,2]           [,3]
[1,]  4.362233e+132  8.724465e+132 -5.658504e+115
[2,] -1.454078e+132 -2.908155e+132  2.782660e+115
[3,] -1.017854e+133 -2.035709e+133  1.038290e+116
>
>
> # 4th version
>
> library(Malmig)
>
> mtx.exp
function (X, n)
{
if (n != round(n)) {
n <- round(n)
warning("rounding exponent `n' to", n)
}
phi <- diag(nrow = nrow(X))
pot <- X
while (n > 0) {
if (n%%2)
phi <- phi %*% pot
n <- n%/%2
pot <- pot %*% pot
}
return(phi)
}
> for(i in 274:277){print(i);print(mtx.exp(A,i))}
[1] 274
[,1]           [,2]           [,3]
[1,]  1.615642e+131  3.231283e+131 -2.156709e+114
[2,] -5.385472e+130 -1.077094e+131  1.218501e+114
[3,] -3.769831e+131 -7.539661e+131  3.460636e+114
[1] 275
[,1]           [,2]           [,3]
[1,]  4.846925e+131  9.693850e+131 -5.325150e+114
[2,] -1.615642e+131 -3.231283e+131  3.083014e+114
[3,] -1.130949e+132 -2.261898e+132  8.310626e+114
[1] 276
[,1]           [,2]           [,3]
[1,]  1.454078e+132  2.908155e+132 -1.297786e+115
[2,] -4.846925e+131 -9.693850e+131  7.750249e+114
[3,] -3.392848e+132 -6.785695e+132  1.950919e+115
[1] 277
[,1]           [,2]           [,3]
[1,]  4.362233e+132  8.724465e+132 -3.108580e+115
[2,] -1.454078e+132 -2.908155e+132  1.932685e+115
[3,] -1.017854e+133 -2.035709e+133  4.433079e+115
>
> # 5th version
>
> matrix.power <- function(mat, n)
+ {
+   # test if mat is a square matrix
+   # treat n < 0 and n = 0 -- this is left as an exercise
+   # trap non-integer n and return an error
+   if (n == 1) return(mat)
+   result <- diag(1, ncol(mat))
+   while (n > 0) {
+     if (n %% 2 != 0) {
+       result <- result %*% mat
+       n <- n - 1
+     }
+     mat <- mat %*% mat
+     n <- n / 2
+   }
+   return(result)
+ }
>
> for(i in 274:277){print(i);print(matrix.power(A,i))}
[1] 274
[,1]           [,2]           [,3]
[1,]  1.615642e+131  3.231283e+131 -2.156709e+114
[2,] -5.385472e+130 -1.077094e+131  1.218501e+114
[3,] -3.769831e+131 -7.539661e+131  3.460636e+114
[1] 275
[,1]           [,2]           [,3]
[1,]  4.846925e+131  9.693850e+131 -5.325150e+114
[2,] -1.615642e+131 -3.231283e+131  3.083014e+114
[3,] -1.130949e+132 -2.261898e+132  8.310626e+114
[1] 276
[,1]           [,2]           [,3]
[1,]  1.454078e+132  2.908155e+132 -1.297786e+115
[2,] -4.846925e+131 -9.693850e+131  7.750249e+114
[3,] -3.392848e+132 -6.785695e+132  1.950919e+115
[1] 277
[,1]           [,2]           [,3]
[1,]  4.362233e+132  8.724465e+132 -3.108580e+115
[2,] -1.454078e+132 -2.908155e+132  1.932685e+115
[3,] -1.017854e+133 -2.035709e+133  4.433079e+115
>
>
> # 6th versio
>
> expM.sd <- function(X,e){Xsd <- eigen(X); Xsd\$vec %*% diag(Xsd\$val^e) %*%t(Xsd\$vec)}
>
> for(i in 274:277){print(i);print(expM.sd(A,i))}
[1] 274
[,1]           [,2]           [,3]
[1,]  8.215127e+129 -2.738376e+129 -1.916863e+130
[2,] -2.738376e+129  9.127919e+128  6.389543e+129
[3,] -1.916863e+130  6.389543e+129  4.472680e+130
[1] 275
[,1]           [,2]           [,3]
[1,]  2.464538e+130 -8.215127e+129 -5.750589e+130
[2,] -8.215127e+129  2.738376e+129  1.916863e+130
[3,] -5.750589e+130  1.916863e+130  1.341804e+131
[1] 276
[,1]           [,2]           [,3]
[1,]  7.393614e+130 -2.464538e+130 -1.725177e+131
[2,] -2.464538e+130  8.215127e+129  5.750589e+130
[3,] -1.725177e+131  5.750589e+130  4.025412e+131
[1] 277
[,1]           [,2]           [,3]
[1,]  2.218084e+131 -7.393614e+130 -5.175530e+131
[2,] -7.393614e+130  2.464538e+130  1.725177e+131
[3,] -5.175530e+131  1.725177e+131  1.207624e+132
>
>

>
>
> > Paul,
> >
> > Your solution based on SVD does not work
>
> Ooops. I really am getting rusty. The idea is based on eigenvalues
> which, of course, are not always the same as singular values.
>
> Paul
> even for the matrix in your example
> > (the reason it worked for e=3, was that it is an odd power and
> since P is a
> > permutation matrix.  It will also work for all odd powers, but
> not for even
> > powers).
> >
> > However, a spectral decomposition will work for all symmetric
> matrices (SVD
> > based exponentiation doesn't work for any matrix).
> >
> > Here is the function to do exponentiation based on spectral
> decomposition:>
> > expM.sd <- function(X,e){Xsd <- eigen(X); Xsd\$vec %*%
> diag(Xsd\$val^e) %*%
> > t(Xsd\$vec)}
> >
> >> P <- matrix(c(.3,.7, .7, .3), ncol=2)
> >> P
> >      [,1] [,2]
> > [1,]  0.3  0.7
> > [2,]  0.7  0.3
> >> P %*% P
> >      [,1] [,2]
> > [1,] 0.58 0.42
> > [2,] 0.42 0.58
> >> expM(P,2)
> >      [,1] [,2]
> > [1,] 0.42 0.58
> > [2,] 0.58 0.42
> >> expM.sd(P,2)
> >      [,1] [,2]
> > [1,] 0.58 0.42
> > [2,] 0.42 0.58
> >
> > Exponentiation based on spectral decomposition should be
> generally more
> > accurate (not sure about the speed).  A caveat is that it will
> only work for
> > real, symmetric or complex, hermitian matrices.  It will not work
> for> asymmetric matrices.  It would work for asymmetric matrix if the
> > eigenvectors are made to be orthonormal.
> >
> > Ravi.
> >
> > ------------------------------------------------------------------
> ----------
> > -------
> >
> >
> > Assistant Professor, The Center on Aging and Health
> >
> > Division of Geriatric Medicine and Gerontology
> >
> > Johns Hopkins University
> >
> > Ph: (410) 502-2619
> >
> > Fax: (410) 614-9625
> >
> > Email: rvaradhan at jhmi.edu
> >
> > Webpage:
> >
> >
> > ------------------------------------------------------------------
> ----------
> > --------
> >
> >
> > -----Original Message-----
> > From: r-help-bounces at stat.math.ethz.ch
> > [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Paul Gilbert
> > Sent: Monday, May 07, 2007 5:16 PM
> > To: Atte Tenkanen
> > Cc: r-help at stat.math.ethz.ch
> > Subject: Re: [R] A function for raising a matrix to a power?
> >
> > You might try this, from 9/22/2006 with subject line Exponentiate
> a matrix:
> >
> > I am getting a bit rusty on some of these things, but I seem to
> recall
> > that there is a numerical advantage (speed and/or accuracy?) to
> > diagonalizing:
> >
> >  > expM <- function(X,e) { v <- La.svd(X); v\$u %*% diag(v\$d^e)
> %*% v\$vt }
> >
> >  > P <- matrix(c(.3,.7, .7, .3), ncol=2)
> >  > P %*% P %*% P
> >       [,1]  [,2]
> > [1,] 0.468 0.532
> > [2,] 0.532 0.468
> >  > expM(P,3)
> >       [,1]  [,2]
> > [1,] 0.468 0.532
> > [2,] 0.532 0.468
> >
> > I think this also works for non-integer, negative, large, and
> complex
> > exponents:
> >
> >  > expM(P, 1.5)
> >           [,1]      [,2]
> > [1,] 0.3735089 0.6264911
> > [2,] 0.6264911 0.3735089
> >  > expM(P, 1000)
> >      [,1] [,2]
> > [1,]  0.5  0.5
> > [2,]  0.5  0.5
> >  > expM(P, -3)
> >         [,1]    [,2]
> > [1,] -7.3125  8.3125
> > [2,]  8.3125 -7.3125
> >  > expM(P, 3+.5i)
> >                   [,1]              [,2]
> > [1,] 0.4713+0.0141531i 0.5287-0.0141531i
> > [2,] 0.5287-0.0141531i 0.4713+0.0141531i
> >  >
> >
> > Paul Gilbert
> >
> > Doran, Harold wrote:
> >
> >  > Suppose I have a square matrix P
> >  >
> >  > P <- matrix(c(.3,.7, .7, .3), ncol=2)
> >  >
> >  > I know that
> >  >
> >  >
> >  >> P * P
> >  >
> >  > Returns the element by element product, whereas
> >  >
> >  >
> >  >
> >  >> P%*%P
> >  >>
> >  >
> >  > Returns the matrix product.
> >  >
> >  > Now, P2 also returns the element by element product. But, is
> there a
> >  > slick way to write
> >  >
> >  > P %*% P %*% P
> >  >
> >  > Obviously, P3 does not return the result I expect.
> >  >
> >  > Thanks,
> >  > Harold
> >  >
> >
> >
> > Atte Tenkanen wrote:
> >> Hi,
> >>
> >> Is there a function for raising a matrix to a power? For example
> if you
> > like to compute A%*%A%*%A, is there an abbreviation similar to A^3?
> >> Atte Tenkanen
> >>
> >>> A=rbind(c(1,1),c(-1,-2))
> >>> A
> >>      [,1] [,2]
> >> [1,]    1    1
> >> [2,]   -1   -2
> >>> A^3
> >>      [,1] [,2]
> >> [1,]    1    1
> >> [2,]   -1   -8
> >>
> >> But:
> >>
> >>> A%*%A%*%A
> >>      [,1] [,2]
> >> [1,]    1    2
> >> [2,]   -2   -5
> >>
> >> ______________________________________________
> >> R-help at stat.math.ethz.ch mailing list
> >> https://stat.ethz.ch/mailman/listinfo/r-help
> > http://www.R-project.org/posting-guide.html
> >> and provide commented, minimal, self-contained, reproducible code.
> >
> ============================================================================> ========
> >
> > La version française suit le texte anglais.
> >
> > ------------------------------------------------------------------
> ----------
> > --------
> >
> > This email may contain privileged and/or confidential
> inform...{{dropped}}>
> > ______________________________________________
> > R-help at stat.math.ethz.ch mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> guide.html> and provide commented, minimal, self-contained,
> reproducible code.
> ====================================================================================
>
> La version française suit le texte anglais.
>
> --------------------------------------------------------------------
> ----------------
>
> This email may contain privileged and/or confidential info...{{dropped}}

```