1##------OVERVIEW---------------------------------------------------------------- 2 3## Input: A; nxn Matrix, no eigenvalues <=0, not singular 4## Output: log(A); Matrixlogarithm; nxn Matrix 5 6 7## Function for Calculation of log(A) with the Inverse Scaling&Squaring Method 8 9## Step 0: Schur Decompostion Tr 10## Step 1: Scaling (root of Tr) 11## Step 2: Padé-Approximation 12## Step 3: Squaring 13## Step 4: Reverse Schur Decomposition 14 15## R-Implementation of Higham's Algorithm from the Book 16## "Functions of Matrices - Theory and Computation", Chapter 11, Algorithm 11.9 17 18##-------CODE------------------------------------------------------------------- 19 20## The coefficients for the Padé-approximation can be computed at install time: 21 22## r: exponents are in (-51):(-56) 23## p: exponents are in c((-47):(-53), -56) 24 25 26logm.H08.r <- 27 rbind(c(5003999585967230*2^(-54), 8006399337547537*2^(-54), 5/18, 0,0,0,0), 28 c(5640779706068081*2^(-51), 8899746432686114*2^(-53), 29 8767290225458872*2^(-54), 6733946100265013*2^(-55), 0,0,0), 30 c(5686538473148996*2^(-51), 4670441098084653*2^(-52), 31 5124095576030447*2^(-53), 5604406634440294*2^(-54), 32 8956332917077493*2^(-56), 0,0), 33 c(5712804453675980*2^(-51), 4795663223967718*2^(-52), 34 5535461316768070*2^(-53), 6805310445892841*2^(-54), 35 7824302940658783*2^(-55), 6388318485698934*2^(-56), 0), 36 c(5729264333934497*2^(-51), 4873628951352824*2^(-52), 37 5788422587681293*2^(-53), 7529283295392226*2^(-54), 38 4892742764696865*2^(-54), 5786545115272933*2^(-55), 39 4786997716777457*2^(-56))) 40 41logm.H08.p <- 42 - rbind(c(7992072898328873*2^(-53), 1/2, 8121010851296995*2^(-56), 0,0,0,0), 43 c(8107950463991866*2^(-49), 6823439817291852*2^(-51), 44 6721885580294475*2^(-52), 4839623620596807*2^(-52), 0,0,0), 45 c(6000309411699298*2^(-48), 4878981751356277*2^(-50), 2, 46 5854649940415304*2^(-52), 4725262033344781*2^(-52),0,0), 47 c(8336234321115872*2^(-48), 6646582649377394*2^(-50), 48 5915042177386279*2^(-51), 7271968136730531*2^(-52), 49 5422073417188307*2^(-52), 4660978705505908*2^(-52), 0), 50 c(5530820008925390*2^(-47), 8712075454469181*2^(-50), 51 7579841581383744*2^(-51), 4503599627370617*2^(-51), 52 6406963985981958*2^(-52), 5171999978649488*2^(-52), 53 4621190647118544*2^(-52))) 54 55 56logm.Higham08 <- function(x) { 57 ## work with "Matrix" too: x<-as.matrix(x) 58 59 ##MM: No need to really check here; we get correct error msg later anyway 60 ## and don't need to compute det() here, in the good cases ! 61 ## if (det(x) == 0) stop("'x' is singular") 62 63 ##-------Step 0: Schur Decomposition----------------------------------------- 64 65 ## Schur() checks for square matrix also: 66 Sch.x <- Schur(Matrix(x, sparse=FALSE)) 67 ## FIXME 'sparse=FALSE' is workaround - good as long Matrix has no sparse Schur() 68 ev <- Sch.x@EValues 69 if(getOption("verbose") && any(abs(Arg(ev) - pi) < 1e-7)) 70## Let's see what works: temporarily *NOT* stop()ping : 71 message(sprintf("'x' has negative real eigenvalues; maybe ok for %s", "logm()")) 72 n <- Sch.x@Dim[1] 73 Tr <- as.matrix(Sch.x@T) 74 Q <- as.matrix(Sch.x@Q) 75 76 ##----- Step 1: [Inverse] Scaling ------------------------------------------- 77 I <- diag(n) 78 thMax <- 0.264 79 theta <- c(0.0162, 0.0539, 0.114, 0.187, thMax) 80 p <- k <- 0 ; t.o <- -1 81 ## NB: The following could loop forever, e.g., for logm(Diagonal(x=1:0)) 82 repeat{ 83 t <- norm(Tr - I, "1") # norm(x, .) : currently x is coerced to dgeMatrix 84 if(is.na(t)) { 85 warning(sprintf(ngettext(k, 86 "NA/NaN from %s after %d step.\n", 87 "NA/NaN from %s after %d steps.\n"), 88 " || Tr - I || ", k), 89 "The matrix logarithm may not exist for this matrix.") 90 return(array(t, dim=dim(Tr))) 91 } 92 if (t < thMax) { 93 ## FIXME: use findInterval() 94 j2 <- which.max( t <= theta) 95 j1 <- which.max( (t/2) <= theta) 96 if ((j2-j1 <= 1) || ((p <- p+1) == 2)) { 97 m <- j2 ## m := order of the Padé-approximation 98 break 99 } 100 } else if(k > 20 && abs(t.o - t) < 1e-7*t) { 101 ## 102 warning(sprintf("Inverse scaling did not work (t = %g).\n", t), 103 "The matrix logarithm may not exist for this matrix.", 104 "Setting m = 3 arbitrarily.") 105 m <- 3 106 break 107 } 108 Tr <- rootS(Tr)##--> Matrix Square root of Jordan T 109 ## ----- [see below; compare with ./sqrtm.R 110 t.o <- t 111 k <- k+1 112 } 113 if(getOption("verbose")) 114 message(sprintf("logm.Higham08() -> (k, m) = (%d, %d)", k,m)) 115 116 ##------ Step 2: Padé-Approximation ----------------------------------------- 117 118 ## of order m : 119 r.m <- logm.H08.r[m,] 120 p.m <- logm.H08.p[m,] 121 X <- 0 122 Tr <- Tr-I 123 for (s in 1:(m+2)) { 124 X <- X + r.m[s]*solve(Tr - p.m[s]*I, Tr) 125 } 126 127 ##--- Step 3 & 4: Squaring & reverse Schur Decomposition ----------------- 128 129 2^k* Q %*% X %*% solve(Q) 130} 131 132### --- was rootS.r ----------------------------------------------------------- 133### ~~~~~~~ 134 135##------OVERVIEW---------------------------------------------------------------- 136 137## Input: UT; nxn upper triangular block matrix (real Schur decomposition) 138## Output: root of matrix UT, nxn upper triangular Matrix 139 140## Function for calculation of UT^(1/2), which is used for the logarithm function 141 142## Step 0: Analyse block structure 143## Step 1: Calculate diagonal elements/blocks 144## Step 2: Calculate superdiagonal elements/blocks 145 146## R-Implementation of Higham's Algorithm from the Book 147## "Functions of Matrices - Theory and Computation", Chapter 6, Algorithm 6.7 148 149## NB: Much in parallel with sqrtm() in ./sqrtm.R <<< keep in sync 150## ~~~~~ ~~~~~~~ 151rootS <- function(x) { 152 ## Generate Basic informations of Matrix x 153 stopifnot(length(d <- dim(x)) == 2, is.numeric(d), 154 (n <- d[1]) == d[2], n >= 1) 155 ## FIXME : should work for "Matrix" too: not S <- as.matrix(x) 156 S <- x 157 158 ##------- STEP 0: Analyse block structure ---------------------------------- 159 if(n > 1L) { 160 ## Count 2x2 blocks (as Schur(x) is the real Schur Decompostion) 161 J.has.2 <- S[cbind(2:n, 1:(n-1))] != 0 162 k <- sum(J.has.2) ## := number of non-zero SUB-diagonals 163 } else k <- 0L 164 165 ## Generate Blockstructure and save it as R.index 166 R.index <- vector("list",n-k) 167 l <- 1L 168 i <- 1L 169 while(i < n) { ## i advances by 1 or 2, depending on 1- or 2- Jordan Block 170 if (S[i+1L,i] == 0) { 171 R.index[[l]] <- i 172 } 173 else { 174 i1 <- i+1L 175 R.index[[l]] <- c(i,i1) # = i:(i+1) 176 i <- i1 177 } 178 i <- i+1L 179 l <- l+1L 180 } 181 if (is.null(R.index[[n-k]])) { # needed; FIXME: should be able to "know" 182 ##message(sprintf("R.index[n-k = %d]] is NULL, set to n=%d", n-k,n)) 183 R.index[[n-k]] <- n 184 } 185 186 ##---------STEP 1: Calculate diagonal elements/blocks------------------------ 187 ## Calculate the root of the diagonal blocks of the Schur Decompostion S 188 I <- diag(2) 189 X <- matrix(0,n,n) 190 for (j in seq_len(n-k)) { 191 ij <- R.index[[j]] 192 if (length(ij) == 1L) { 193 ## Sij <- S[ij,ij] 194 ## if(Sij < 0) 195 ## ## FIXME(?) : in sqrtm(), we take *complex* sqrt() if needed : 196 ## ## ----- but afterwards norm(Tr - I, "1") fails with complex 197 ## ## Sij <- complex(real = Sij, imaginary = 0) 198 ## stop("negative diagonal entry -- matrix square does not exist") 199 ## X[ij,ij] <- sqrt(Sij) 200 X[ij,ij] <- sqrt(S[ij,ij]) 201 } 202 else { 203 ## "FIXME"(better algorithm): only need largest eigen value 204 ev1 <- eigen(S[ij,ij], only.values=TRUE)$values[1] 205 r1 <- Re(sqrt(ev1)) ## sqrt(<complex>) ... 206 X[ij,ij] <- 207 r1*I + 1/(2*r1)*(S[ij,ij] - Re(ev1)*I) 208 } 209 } 210 211### ___ FIXME __ code re-use: All the following is identical to 'STEP 3' in sqrtm() 212### ----- and almost all of STEP 1 above is == 'STEP 2' of sqrtm() 213 214 ##---------STEP 2: Calculate superdiagonal elements/blocks------------------- 215 216 ## Calculate the remaining, not-diagonal blocks 217 if (n-k > 1L) for (j in 2L:(n-k)) { 218 ij <- R.index[[j]] 219 for (i in (j-1L):1L) { 220 ii <- R.index[[i]] 221 sumU <- 0 222 223 ## Calculation for 1x1 Blocks 224 if (length(ij) == 1L & length(ii) == 1L ) { 225 if (j-i > 1L) for (l in (i+1L):(j-1L)) { 226 il <- R.index[[l]] 227 sumU <- sumU + { 228 if (length(il) == 2 ) X[ii,il]%*%X[il,ij] 229 else X[ii,il] * X[il,ij] 230 } 231 } 232 X[ii,ij] <- solve(X[ii,ii]+X[ij,ij],S[ii,ij]-sumU) 233 } 234 235 ## Calculation for 1x2 Blocks 236 else if (length(ij) == 2 & length(ii) == 1L ) { 237 if (j-i > 1L) for (l in(i+1L):(j-1L)) { 238 il <- R.index[[l]] 239 sumU <- sumU + { 240 if (length(il) == 2) X[ii,il]%*%X[il,ij] 241 else X[ii,il] * X[il,ij] 242 } 243 } 244 X[ii,ij] <- solve(t(X[ii,ii]*I + X[ij,ij]), 245 as.vector(S[ii,ij] - sumU)) 246 } 247 248 ## Calculation for 2x1 Blocks 249 else if (length(ij) == 1L & length(ii) == 2 ) { 250 if (j-i > 1L) for (l in(i+1L):(j-1L)) { 251 il <- R.index[[l]] 252 sumU <- sumU + { 253 if (length(il) == 2 ) X[ii,il]%*%X[il,ij] 254 else X[ii,il] * X[il,ij] 255 } 256 } 257 X[ii,ij] <- solve(X[ii,ii]+X[ij,ij]*I,S[ii,ij]-sumU) 258 } 259 ## Calculation for 2x2 Blocks with special equation for solver 260 else if (length(ij) == 2 & length(ii) == 2 ) { 261 if (j-i > 1L) for (l in(i+1L):(j-1L)) { 262 il <- R.index[[l]] 263 sumU <- sumU + { 264 if (length(il) == 2 ) X[ii,il] %*% X[il,ij] 265 else X[ii,il] %*% t(X[il,ij]) 266 } 267 } 268 tUii <- matrix(0,4,4) 269 tUii[1:2,1:2] <- X[ii,ii] 270 tUii[3:4,3:4] <- X[ii,ii] 271 tUjj <- matrix(0,4,4) 272 tUjj[1:2,1:2] <- t(X[ij,ij])[1L,1L]*I 273 tUjj[3:4,3:4] <- t(X[ij,ij])[2L,2L]*I 274 tUjj[1:2,3:4] <- t(X[ij,ij])[1L,2L]*I 275 tUjj[3:4,1:2] <- t(X[ij,ij])[2L,1L]*I 276 X[ii,ij] <- solve(tUii+tUjj,as.vector(S[ii,ij]-sumU)) 277 } 278 } 279 } 280 X 281} 282 283