[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
>
>
>
>
> Ravi Varadhan wrote:
> > 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.
> >
> > ------------------------------------------------------------------
> ----------
> > -------
> >
> > Ravi Varadhan, Ph.D.
> >
> > 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:
> http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html>
> >
> >
> > ------------------------------------------------------------------
> ----------
> > --------
> >
> >
> > -----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
> >> PLEASE do read the posting guide
> > 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
> > PLEASE do read the posting guide 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 info...{{dropped}}
More information about the R-help
mailing list