1### ===== File part of R package expm ===== 2### 3### Function to compute the matrix exponential 4### 5### exp(M) = sum(n = 0:Inf; M^n / n!), 6### 7### where M is an (n x n) matrix. 8### 9 10expm.s.Pade.s <- function(x, order, n=nrow(x)) { 11 ## no checking here; this is not to be called by the user 12 ## try have this work with "mpfrMatrix" <==> solve(<mpfrArray>) 13 14 ## s := the number of [s]quarings 15 e <- ceiling(log2(max(rowSums(abs(x))))) 16 s <- max(e+1, 0) 17 18 ## preconditions x : 19 x <- x / (2^s) 20 21 ## Pade approximation for exp(x) 22 c <- .5 23 D <- diag(1, n) 24 E <- D + c*x 25 D <- D - c*x 26 X <- x 27 p <- TRUE 28 for(k in 2:order) { 29 c <- c * (order-k+1) / (k*(2*order-k+1)) 30 X <- x %*% X # now X = x ^ k 31 cX <- c*X 32 E <- E + cX 33 D <- if(p) D + cX else D - cX 34 p <- !p 35 } 36 E <- solve(D, E) 37 38 ## Undo the scaling by repeated squaring : 39 for(k in seq_len(s)) 40 E <- E %*% E 41 E 42} 43 44 45expm.methSparse <- c("Higham08", "R_Eigen", "R_Pade") 46## keep this list up-to-date - test by setting R_EXPM_NO_DENSE_COERCION 47## but NOTE: It may make sense to *keep* the message() about coercion to dense (memory blow up!) 48 49expm <- function(x, method = c("Higham08.b", "Higham08", 50 "AlMohy-Hi09", 51 "Ward77", "PadeRBS", "Pade", "Taylor", "PadeO", "TaylorO", 52 "R_Eigen", "R_Pade", "R_Ward77", "hybrid_Eigen_Ward"), 53 order = 8, 54 trySym = TRUE, tol = .Machine$double.eps, 55 do.sparseMsg = TRUE, 56 preconditioning = c("2bal", "1bal", "buggy")) 57{ 58 ## some methods work for "matrix", "Matrix", or "mpfrMatrix", iff solve(.,.) worked: 59 stopifnot(is.numeric(x) || (isM <- inherits(x, "dMatrix")) || inherits(x, "mpfrMatrix")) 60 if(length(d <- dim(x)) != 2) stop("argument is not a matrix") 61 if (d[1] != d[2]) stop("matrix not square") 62 method <- match.arg(method) 63 checkSparse <- !nzchar(Sys.getenv("R_EXPM_NO_DENSE_COERCION")) 64 isM <- !is.numeric(x) && isM 65 if(isM && checkSparse) { # i.e., a "dMatrix" 66 if(!(method %in% expm.methSparse) && is(x, "sparseMatrix")) { 67 if(do.sparseMsg) 68 message("coercing to dense matrix, as required by method ", 69 dQuote(method)) 70 x <- as(x, "denseMatrix") 71 } 72 } 73 switch(method, 74 "AlMohy-Hi09" = expm.AlMoHi09(x, p = order) 75 , 76 "Higham08.b" = expm.Higham08(x, balancing = TRUE) 77 , 78 "Higham08" = expm.Higham08(x, balancing = FALSE) 79 , 80 "Ward77" = { 81 ## AUTHORS: Christophe Dutang, Vincent Goulet at act ulaval ca 82 ## built on "Matrix" package, built on 'octave' code 83 ## Martin Maechler, for the preconditioning etc 84 if(!is.numeric(x)) x <- as(x, "matrix") 85 switch(match.arg(preconditioning), 86 "2bal" = .Call(do_expm, x, "Ward77"), 87 "1bal" = .Call(do_expm, x, "Ward77_1"), 88 "buggy"= .Call(do_expm, x, "buggy_Ward77"), 89 stop("invalid 'preconditioning'")) 90 }, 91 "R_Eigen" = { 92 ## matrix exponential using eigenvalues / spectral decomposition : 93 ## == Dubious Way 'Method 14' : is 94 ## good for 'symmetric' or 'orthogonal' (or other 'normal' : A'A = AA' ): 95 96 ## MM: improved from mexp2() with 'trySym' and isSymmetric() 97 isSym <- if(trySym) isSymmetric.matrix(x) else FALSE 98 z <- eigen(x, symmetric = isSym) 99 V <- z$vectors 100 Vi <- if(isSym) t(V) else solve(V) 101 Re(V %*% ( exp(z$values) * Vi)) ## == 102 ##(V %*% diag(exp(z$values)) %*% Vi) 103 }, 104 "hybrid_Eigen_Ward" = { 105 ## AUTHOR: Christophe Dutang 106 ## matrix exponential using eigenvalues / spectral decomposition and 107 ## Ward(1977) algorithm if x is numerically non diagonalisable 108 if(!is.numeric(x)) x <- as(x, "matrix") 109 .Call(do_expm_eigen, x, tol) 110 }, 111 "R_Pade"= { ## use scaling + Pade + squaring with R code: 112 113 ## matrix exponential using a scaling and squaring algorithm 114 ## with a Pade approximation 115 ## source code translated from expmdemo1.m in Matlab 116 ## by Stig Mortensen <sbm@imm.dtu.dk>, 117 ## prettified by MM -- works for "matrix" or "Matrix" matrices ! 118 119 stopifnot(order >= 2) 120 expm.s.Pade.s(x, order, n=d[1]) 121 }, 122 "R_Ward77" = { ## R implementation of "Ward(1977)" 123 ## also works for "Matrix" matrices 124 stopifnot(order >= 2) 125 n <- d[1] 126 ## Preconditioning Step 1: shift diagonal by average diagonal 127 trShift <- sum(d.x <- diag(x)) 128 if(trShift) { 129 trShift <- trShift/n 130 diag(x) <- d.x - trShift 131 } 132 133 ## Preconditioning Step 2: balancing with balance. 134 ## ------ 135 ## For now, do like the octave implementation 136 ## TODO later: use "B" (faster; better condition of result) 137 baP <- balance(x, "P") 138 baS <- balance(baP$z, "S") 139 140 x <- expm.s.Pade.s(baS$z, order) 141 ## ------------- scaling + Pade + squaring ------ 142 ## i.e., entails Preconditioning Step 3 (and its reverse) 143 144 ## Reverse step 2: ------------------ 145 ## 146 ## Step 2 b: apply inverse scaling 147 d <- baS$scale 148 x <- x * (d * rep(1/d, each = n)) 149 ## 150 ## Step 2 a: apply inverse permutation (of rows and columns): 151 pp <- as.integer(baP$scale) 152 if(baP$i1 > 1) { ## The lower part 153 for(i in (baP$i1-1):1) { # 'p1' in *reverse* order 154 tt <- x[,i]; x[,i] <- x[,pp[i]]; x[,pp[i]] <- tt 155 tt <- x[i,]; x[i,] <- x[pp[i],]; x[pp[i],] <- tt 156 } 157 } 158 if(baP$i2 < n) { ## The upper part 159 for(i in (baP$i2+1):n) { # 'p2' in *forward* order 160 ## swap i <-> pp[i] both rows and columns 161 tt <- x[,i]; x[,i] <- x[,pp[i]]; x[,pp[i]] <- tt 162 tt <- x[i,]; x[i,] <- x[pp[i],]; x[pp[i],] <- tt 163 } 164 } 165 166 ## reverse step 1 (diagonal shift) 167 if(trShift) { 168 exp(trShift) * x 169 } 170 else x 171 }, 172 "PadeRBS" = { ## the "expofit" method by Roger B. Sidje (U.Queensland, AU) 173 if(!is.numeric(x)) x <- as(x, "matrix") 174 stopifnot((order <- as.integer(order)) >= 1) 175 if(!is.double(x)) storage.mode(x) <- "double" 176 Fobj <- .Fortran(matexpRBS, 177 order, # IDEG 1 178 as.integer(d[1]), # M 2 179 T = 1., # T 3 180 H = x, # H 4 181 iflag = integer(1) # IFLAG 5 182 )[c("H","iflag")] 183 if(Fobj[["iflag"]] < 0) 184 stop("Unable to determine matrix exponential") 185 Fobj[["H"]] 186 } 187 , { ## the "mexp" methods by 188 ## AUTHORS: Marina Shapira and David Firth -------------- 189 if(!is.numeric(x)) x <- as(x, "matrix") 190 if(!is.double(x)) storage.mode(x) <- "double" 191 order <- as.integer(order) 192 ## MM: a "silly" way to code the method / order 193 ntaylor <- npade <- 0L 194 if (substr(method,1,4) == "Pade") 195 npade <- order else ntaylor <- order 196 res <- if(identical(grep("O$", method), 1L)) 197 .Fortran(matrexpO, 198 X = x, 199 size = d[1], 200 ntaylor, 201 npade, 202 accuracy = double(1))[c("X", "accuracy")] 203 else 204 .Fortran(matrexp, 205 X = x, 206 size = d[1], 207 ntaylor, 208 npade, 209 accuracy = double(1))[c("X", "accuracy")] 210 structure(res$X, accuracy = res$accuracy) 211 })## end{switch} 212} 213 214