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

From: Atte Tenkanen <attenka_at_utu.fi>
Date: Wed, 09 May 2007 17:07:15 +0300

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}}

______________________________________________
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. Received on Wed 09 May 2007 - 14:10:57 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 09 May 2007 - 15:31:14 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.