R Under development (unstable) (2024-01-10 r85799 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > #### Examples where we know the result "exactly" > > library(expm) Loading required package: Matrix Attaching package: 'expm' The following object is masked from 'package:Matrix': expm > > options(digits = 4, width = 90, keep.source = FALSE) > > mSource <- function(file, ...) + source(system.file(file, ..., package = "expm", mustWork=TRUE)) > mSource("test-tools.R")## -> assertError(), rMat(), .. doExtras Loading required package: tools > mSource("demo", "exact-fn.R") > doExtras [1] FALSE > > re.nilA3 <- function(xyz, EXPMlist) + { + stopifnot(is.list(EXPMlist)) + r <- do.call(nilA3, as.list(xyz)) + sapply(EXPMlist, function(Efn) relErr(r$expA, Efn(r$A))) + } > > re.facMat <- function(n, EXPMlist, rFUN = rnorm, ...) + { + stopifnot(is.list(EXPMlist)) + r <- facMat(n, rFUN, ...) + vapply(EXPMlist, function(EXPM) { + ct <- system.time(E <- EXPM(r$A), gc = FALSE)[[1]] + c(relErr = relErr(r$expA, E), c.time = ct) + }, double(2)) + } > > re.m2ex3 <- function(eps, EXPMlist) + { + stopifnot(is.list(EXPMlist)) + r <- m2ex3(eps) + sapply(EXPMlist, function(EXPM) relErr(r$expA, EXPM(r$A))) + } > > ## check 1x1 matrices > stopifnot( + ## these had failed before 2017-03-28 (= Liselotte's 58-th birthday): + all.equal(as.matrix(sqrtm(matrix(4))), matrix(2)) , + all.equal(as.matrix(logm (matrix(pi))), matrix(log(pi))) , + ## these had "always" worked : + all.equal(as.matrix(expm (matrix(0))), matrix(1)) , + all.equal(as.matrix(expm (matrix(1))), matrix(exp(1))) + ) > > > set.seed(321) > re <- replicate(1000, + c(re.nilA3(rlnorm(3),list(function(x)expm(x,"Pade"))), + re.nilA3(rnorm(3), list(function(x)expm(x,"Pade"))))) > > summary(t(re)) V1 V2 Min. :0.00e+00 Min. :0.00e+00 1st Qu.:0.00e+00 1st Qu.:0.00e+00 Median :2.15e-17 Median :1.47e-17 Mean :3.15e-17 Mean :2.27e-17 3rd Qu.:4.98e-17 3rd Qu.:3.61e-17 Max. :2.71e-16 Max. :1.29e-16 > stopifnot(rowMeans(re) < 1e-15, + apply(re, 1, quantile, 0.80) < 1e-16, + apply(re, 1, quantile, 0.90) < 2e-15, + apply(re, 1, max) < c(4e-14, 6e-15)) > > showProc.time() Time (user system elapsed): 0.36 0.02 0.38 > > ## Check *many* random nilpotent 3 x 3 matrices: > set.seed(321) > RE <- replicate(1000, + c(re.nilA3(rlnorm(3), list(function(x) expm(x, "Ward77"))), + re.nilA3(rnorm(3), list(function(x) expm(x, "Ward77"))))) > stopifnot(rowMeans(RE) < 1e-15, + apply(RE, 1, quantile, 0.80) < 1e-16, + apply(RE, 1, quantile, 0.90) < 2e-15, + apply(RE, 1, max) < c(4e-14, 6e-15)) > > print(summary(t(RE))) V1 V2 Min. :0.00e+00 Min. :0.00e+00 1st Qu.:0.00e+00 1st Qu.:0.00e+00 Median :0.00e+00 Median :0.00e+00 Mean :1.49e-17 Mean :7.03e-18 3rd Qu.:2.70e-17 3rd Qu.:5.20e-18 Max. :1.87e-16 Max. :7.52e-17 > epsC <- .Machine$double.eps > cat("relErr(expm(.,Pade)) - relErr(expm(.,Ward77)) in Machine_eps units:\n") relErr(expm(.,Pade)) - relErr(expm(.,Ward77)) in Machine_eps units: > print(summary(c(re - RE)) / epsC) Min. 1st Qu. Median Mean 3rd Qu. Max. -0.8429 0.0000 0.0000 0.0727 0.1065 1.2215 > ## Min. 1st Qu. Median Mean 3rd Qu. Max. > ## -0.6183442 0.0000000 0.0000000 1.3650410 0.1399719 94.9809161 > ## nb-mm3; ditto lynne (both x64), 2014-09-11: > ## Min. 1st Qu. Median Mean 3rd Qu. Max. > ## -0.8422 0.0000 0.0000 0.0725 0.1067 1.2205 > ## 32-bit [i686, florence, Linx 3.14.8-100.fc19..]: > ## Min. 1st Qu. Median Mean 3rd Qu. Max. > ## -0.62 0.00 0.00 1.36 0.14 95.93 > > > showProc.time() Time (user system elapsed): 0.24 0 0.23 > > ###--- A second group --- where we know the diagonalization of A --- > > if(!require("Matrix")) + q('no') > ## ------ the rest really uses 'Matrix' > ##---> now use expm::expm() since Matrix has its own may mask the expm one > ## ^^^^^^^^^^^^ > > if(!dev.interactive(orNone=TRUE)) pdf("expm_exact-ex.pdf") > > > ## rMat() relies on Matrix::rcond(): > ## Now with the change default rcondMin, this "works" > R40 <- rMat(40) > R80 <- rMat(80) > showProc.time() Time (user system elapsed): 0.06 0 0.07 > > expm.safe.Eigen <- function(x, silent = FALSE) { + r <- try(expm::expm(x, "R_Eigen"), silent = silent) + if(inherits(r, "try-error")) NA else r + } > > ## the S4 generic > Matrix::expm standardGeneric for "expm" defined from package "Matrix" function (x) standardGeneric("expm") Methods may be defined for arguments: x Use showMethods(expm) for currently available ones. > ## the dgeMatrix method - adapted to Matrix changes, had *more versatile* ..2dge() : > expm.Matr.dge <- function(x) getDataPart(getMethod("expm", "dgeMatrix"))(Matrix::.m2dense(x)) > expmList <- + list(Matr = Matrix::expm, + Matr.d = expm.Matr.dge, + Ward = function(x) expm::expm(x, "Ward77"), + s.P.s = function(x) expm::expm(x, "Pade"), + s.P.sO= function(x) expm::expm(x, "PadeO"), + s.P.sRBS= function(x) expm::expm(x, "PadeRBS"), + sPs.H08.= function(x) expm:: expm.Higham08(x, balancing=FALSE), + sPs.H08b= function(x) expm:: expm.Higham08(x, balancing= TRUE), + AmHi09.06= function(x) expm:::expm.AlMoHi09(x, p = 6), + AmHi09.08= function(x) expm:::expm.AlMoHi09(x, p = 8), + AmHi09.10= function(x) expm:::expm.AlMoHi09(x, p = 10), + AmHi09.12= function(x) expm:::expm.AlMoHi09(x, p = 12), + AmHi09.13= function(x) expm:::expm.AlMoHi09(x, p = 13), + s.T.s = function(x) expm::expm(x, "Taylor"), + s.T.sO= function(x) expm::expm(x, "TaylorO"), + Eigen = expm.safe.Eigen, + hybrid= function(x) expm::expm(x, "hybrid") + ) > > > set.seed(12) > fRE <- replicate(if(doExtras) 100 else 20, + re.facMat(20, expmList)) > cat("Number of correct decimal digits for facMat(20, rnorm):\n") Number of correct decimal digits for facMat(20, rnorm): > summary(-log10(t(fRE["relErr",,]))) Matr Matr.d Ward s.P.s s.P.sO Min. :13.7 Min. :13.7 Min. :13.7 Min. :12.8 Min. :12.8 1st Qu.:14.2 1st Qu.:14.2 1st Qu.:14.2 1st Qu.:13.1 1st Qu.:13.1 Median :14.3 Median :14.3 Median :14.3 Median :13.2 Median :13.2 Mean :14.3 Mean :14.3 Mean :14.3 Mean :13.2 Mean :13.2 3rd Qu.:14.5 3rd Qu.:14.5 3rd Qu.:14.5 3rd Qu.:13.3 3rd Qu.:13.3 Max. :14.7 Max. :14.7 Max. :14.7 Max. :13.6 Max. :13.6 s.P.sRBS sPs.H08. sPs.H08b AmHi09.06 AmHi09.08 Min. :13.9 Min. :14.4 Min. :14.4 Min. : 6.61 Min. : 9.78 1st Qu.:14.3 1st Qu.:14.8 1st Qu.:14.8 1st Qu.: 7.69 1st Qu.:11.27 Median :14.4 Median :14.9 Median :14.9 Median : 8.24 Median :11.72 Mean :14.4 Mean :14.9 Mean :14.9 Mean : 8.47 Mean :12.19 3rd Qu.:14.6 3rd Qu.:15.0 3rd Qu.:15.0 3rd Qu.: 9.41 3rd Qu.:13.56 Max. :14.8 Max. :15.2 Max. :15.2 Max. :10.35 Max. :14.40 AmHi09.10 AmHi09.12 AmHi09.13 s.T.s s.T.sO Min. :13.4 Min. :14.5 Min. :14.5 Min. :13.0 Min. :13.0 1st Qu.:14.6 1st Qu.:14.7 1st Qu.:14.7 1st Qu.:13.4 1st Qu.:13.4 Median :14.7 Median :14.8 Median :14.8 Median :13.5 Median :13.5 Mean :14.6 Mean :14.8 Mean :14.8 Mean :13.5 Mean :13.5 3rd Qu.:14.8 3rd Qu.:14.9 3rd Qu.:14.9 3rd Qu.:13.7 3rd Qu.:13.7 Max. :15.0 Max. :15.2 Max. :15.2 Max. :13.9 Max. :13.9 Eigen hybrid Min. :14.2 Min. :14.2 1st Qu.:14.4 1st Qu.:14.5 Median :14.5 Median :14.5 Mean :14.5 Mean :14.5 3rd Qu.:14.6 3rd Qu.:14.6 Max. :14.8 Max. :14.8 > > > ## Now look at that: > boxplot(t(fRE["relErr",,]), log="y", notch=TRUE, ylim = c(8e-16, 1e-8), + main = "relative errors for 'random' eigen-ok 20 x 20 matrix") Warning message: In (function (z, notch = FALSE, width = NULL, varwidth = FALSE, : some notches went outside hinges ('box'): maybe set notch=FALSE > > showProc.time() Time (user system elapsed): 0.22 0.03 0.25 > > if(doExtras) { + str(rf100 <- replicate(20, re.facMat(100, expmList))) + print(1000*t(apply(rf100["c.time",,], 1, summary))) + ## lynne {Linux 2.6.34.7-56.fc13.x86_64 --- AMD Phenom II X4 925}: + ## Min. 1st Qu. Median Mean 3rd Qu. Max. + ## Ward 23 24 24.5 24.4 25.0 25 + ## s.P.s 107 109 109.0 109.0 109.0 112 + ## s.P.sO 188 190 191.0 192.0 193.0 198 + ## s.P.sRBS 17 18 19.0 18.9 19.2 21 + ## sPs.H08. 15 17 18.0 17.6 18.0 19 + ## sPs.H08b 18 18 19.0 23.4 20.0 107 + ## s.T.s 44 45 45.0 45.6 46.0 48 + ## s.T.sO 96 98 99.0 100.0 100.0 116 + ## Eigen 18 19 20.0 24.4 21.0 109 + ## hybrid 40 42 42.0 47.1 44.0 133 + + ##--> take out the real slow ones for the subsequent tests: + `%w/o%` <- function(x, y) x[!x %in% y] #-- x without y + print(nms.swift <- names(expmList) %w/o% + c("s.P.s", "s.P.sO", "s.T.s", "s.T.sO")) + expmL.swift <- expmList[nms.swift] + + set.seed(18) ## 12 replicates is too small .. but then it's too slow otherwise: + rf400 <- replicate(12, re.facMat(400, expmL.swift)) + print(1000*t(apply(rf400["c.time",,], 1, summary))) + ## lynne: + ## Min. 1st Qu. Median Mean 3rd Qu. Max. + ## Ward 1740 1790 1830 1820 1860 1900 + ## s.P.sRBS 1350 1420 1440 1430 1450 1460 + ## sPs.H08. 1020 1030 1130 1140 1210 1290 + ## sPs.H08b 1120 1130 1220 1220 1300 1390 + ## Eigen 962 977 989 992 1000 1030 + ## hybrid 2740 2800 2840 2840 2890 2910 + + showProc.time() + + }## if(doExtras) only > > ## Now try an example with badly conditioned "random" M matrix... > ## ... > ## ... (not yet) > > > ### m2ex3() --- The 2x2 example with bad condition , see A3 in ./ex2.R > > RE <- re.m2ex3(1e-8, expmList) > sort(RE)# Ward + both sps.H08 are best; s.P.s fair, Eigen (and hybrid): ~1e-9 Matr Matr.d Ward s.P.sRBS sPs.H08. sPs.H08b AmHi09.13 AmHi09.12 s.P.s 5.584e-33 5.584e-33 5.584e-33 5.584e-33 5.584e-33 5.584e-33 1.509e-16 2.515e-16 2.113e-15 s.P.sO s.T.s s.T.sO AmHi09.10 Eigen hybrid AmHi09.08 AmHi09.06 2.113e-15 2.113e-15 2.113e-15 2.325e-12 1.836e-09 1.836e-09 4.819e-09 3.156e-06 > > eps <- 10^-(1:18) > t.m2 <- t(sapply(eps, re.m2ex3, EXPMlist = expmList)) Error in solve.default(V) : system is computationally singular: reciprocal condition number = 1.11022e-16 Error in solve.default(V) : system is computationally singular: reciprocal condition number = 1.11022e-16 Error in solve.default(V) : system is computationally singular: reciprocal condition number = 1.11022e-16 > ## --> 3 error messages from solve(V), 5 error messages from try(. "R_Eigen" ...) > > showProc.time() Time (user system elapsed): 0.09 0 0.09 > > cbind(sort(apply(log(t.m2),2, median, na.rm=TRUE))) [,1] sPs.H08. -86.05 sPs.H08b -86.05 Matr -83.62 Matr.d -83.62 Ward -83.62 s.P.sRBS -81.20 AmHi09.13 -36.43 AmHi09.12 -35.92 s.P.s -33.79 s.P.sO -33.79 s.T.s -33.79 s.T.sO -33.79 AmHi09.10 -26.79 hybrid -23.74 Eigen -20.12 AmHi09.08 -19.15 AmHi09.06 -12.67 > ## 'na.rm=TRUE' needed for Eigen which blows up for the last 3 eps > t.m2.ranks <- sort(rowMeans(apply(t.m2, 1, rank))) > cbind(signif(t.m2.ranks, 3)) [,1] sPs.H08. 3.19 sPs.H08b 3.19 s.P.sRBS 3.92 Matr 4.94 Matr.d 4.94 Ward 4.94 AmHi09.13 5.61 AmHi09.12 7.31 s.T.s 10.10 s.T.sO 10.10 s.P.s 11.10 s.P.sO 11.10 hybrid 12.50 AmHi09.10 13.70 Eigen 14.80 AmHi09.08 15.10 AmHi09.06 16.40 > ## lynne (x86_64, Linux 3.14.4-100; Intel i7-4765T), 2014-09: > ## sPs.H08. 2.67 > ## sPs.H08b 2.67 > ## s.P.sRBS 3.06 > ## Ward 4.03 > ## AmHi09.13 4.33 <<- still not close to H08 ! > ## AmHi09.12 5.86 > ## s.T.s 8.33 > ## s.T.sO 8.33 > ## s.P.s 9.11 > ## s.P.sO 9.11 > ## hybrid 10.80 > ## AmHi09.10 11.70 << astonishingly bad > ## Eigen 12.60 > ## AmHi09.08 13.10 > ## AmHi09.06 14.40 > > print(t.m2[, names(t.m2.ranks)[1:8]], digits = 3) sPs.H08. sPs.H08b s.P.sRBS Matr Matr.d Ward AmHi09.13 AmHi09.12 [1,] 9.99e-17 9.99e-17 2.00e-16 6.02e-16 6.02e-16 6.02e-16 3.50e-16 5.52e-16 [2,] 1.01e-16 1.01e-16 2.01e-16 5.53e-16 5.53e-16 5.53e-16 1.01e-16 2.51e-16 [3,] 5.03e-17 5.03e-17 2.51e-16 3.02e-16 3.02e-16 3.02e-16 2.51e-16 3.02e-16 [4,] 2.51e-16 2.51e-16 3.75e-25 1.51e-16 1.51e-16 1.51e-16 1.51e-16 2.51e-16 [5,] 3.52e-16 3.52e-16 2.01e-16 1.51e-16 1.51e-16 1.51e-16 1.01e-16 1.01e-16 [6,] 2.01e-16 2.01e-16 3.52e-16 3.52e-16 3.52e-16 3.52e-16 2.01e-16 3.52e-16 [7,] 3.52e-16 3.52e-16 3.52e-16 6.04e-16 6.04e-16 6.04e-16 7.15e-31 2.01e-16 [8,] 5.58e-33 5.58e-33 5.58e-33 5.58e-33 5.58e-33 5.58e-33 1.51e-16 2.51e-16 [9,] 4.36e-35 4.36e-35 4.36e-35 4.36e-35 4.36e-35 4.36e-35 1.51e-16 2.51e-16 [10,] 0.00e+00 0.00e+00 6.82e-37 0.00e+00 0.00e+00 0.00e+00 1.51e-16 2.51e-16 [11,] 0.00e+00 0.00e+00 5.33e-39 5.33e-39 5.33e-39 5.33e-39 1.51e-16 2.51e-16 [12,] 4.16e-41 4.16e-41 0.00e+00 4.16e-41 4.16e-41 4.16e-41 1.51e-16 2.51e-16 [13,] 0.00e+00 0.00e+00 6.50e-43 6.50e-43 6.50e-43 6.50e-43 1.51e-16 2.51e-16 [14,] 0.00e+00 0.00e+00 5.08e-45 5.08e-45 5.08e-45 5.08e-45 1.51e-16 2.51e-16 [15,] 3.97e-47 3.97e-47 7.94e-47 7.94e-47 7.94e-47 7.94e-47 1.51e-16 2.51e-16 [16,] 0.00e+00 0.00e+00 6.20e-49 6.20e-49 6.20e-49 6.20e-49 1.51e-16 2.51e-16 [17,] 0.00e+00 0.00e+00 4.84e-51 0.00e+00 0.00e+00 0.00e+00 1.51e-16 2.51e-16 [18,] 3.78e-53 3.78e-53 0.00e+00 3.78e-53 3.78e-53 3.78e-53 1.51e-16 2.51e-16 > ## ==> 1st class: H08 (both) and (but slightly better than) Ward > ## 2nd class s.T.s and s.P.s > ## "bad" : hybrid and Eigen > > ## ??? AmHi09 - methods, up to order = 10 are worse ! > if(require(RColorBrewer)) { + ## Bcol <- brewer.pal(ncol(t.m2),"Dark2") + Bcol <- brewer.pal(min(9, ncol(t.m2)), "Set1") + Bcol <- Bcol[sqrt(colSums(col2rgb(Bcol)^2)) < 340] + ## FIXME: more colors ==> ~/R/MM/GRAPHICS/color-palettes.R + } else { + ## 7 from Dark2 + ## Bcol <- c("#1B9E77", "#D95F02", "#7570B3", "#E7298A", + ## "#66A61E", "#E6AB02", "#A6761D") + ## Rather: those from "Set1" + Bcol <- c("#E41A1C", "#377EB8", "#4DAF4A", + "#984EA3", "#FF7F00", # too bright: "#FFFF33", + "#A65628", "#F781BF", "#999999") + } Loading required package: RColorBrewer > > matplot(eps, t.m2, type = "b", log = "xy", col=Bcol, lty = 1:9, pch=1:15, + axes=FALSE, frame = TRUE, + xlab = expression(epsilon), ylab = "relative error", + main = expression(expm(A, method == "*") *" relative errors for " * + A == bgroup("[", atop({-1} *" "* 1, {epsilon^2} *" "*{-1}), "]"))) Warning messages: 1: In xy.coords(x, y, xlabel, ylabel, log = log, recycle = TRUE) : 21 y values <= 0 omitted from logarithmic plot 2: In xy.coords(x, y, xlabel, ylabel, log) : 2 y values <= 0 omitted from logarithmic plot > legend("bottomright",colnames(t.m2), col=Bcol, lty = 1:9, pch=1:15, + inset = 0.02) > if(require("sfsmisc")) { + sfsmisc::eaxis(1, labels=FALSE) + sfsmisc::eaxis(1, at = eps[c(TRUE,FALSE)]) + sfsmisc::eaxis(2) + ## sfsmisc::eaxis(2, labels=FALSE) + ## op <- par(las=2) + ## sfsmisc::eaxis(2, at = axTicks(2,log=TRUE)[c(TRUE,FALSE,FALSE)]) + ## par(op) + } else { + axis(1) + axis(2) + } Loading required package: sfsmisc Attaching package: 'sfsmisc' The following objects are masked _by_ '.GlobalEnv': relErr, relErrV > > ## typical case: > ep <- 1e-10 > (me <- m2ex3(ep)) $A [,1] [,2] [1,] -1e+00 1 [2,] 1e-20 -1 $expA [,1] [,2] [1,] 3.679e-01 0.3679 [2,] 3.679e-21 0.3679 > me$expA * exp(1) ## the correct value ; numerically identical to simple matrix: [,1] [,2] [1,] 1e+00 1 [2,] 1e-20 1 > ## identical() not fulfilled e.g. on Solaris > stopifnot(all.equal(me$expA * exp(1), + rbind(c( 1, 1), + c(ep^2, 1)), + tolerance = 1e-14)) > ## The relative error (matrices): > lapply(expmList, function(EXPM) 1 - EXPM(me$A)/me$expA) $Matr 2 x 2 Matrix of class "dgeMatrix" [,1] [,2] [1,] 0 0 [2,] 0 0 $Matr.d 2 x 2 Matrix of class "dgeMatrix" [,1] [,2] [1,] 0 0 [2,] 0 0 $Ward [,1] [,2] [1,] 0 0 [2,] 0 0 $s.P.s [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 1.665e-16 $s.P.sO [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 1.665e-16 $s.P.sRBS [,1] [,2] [1,] 0.00e+00 0 [2,] -2.22e-16 0 $sPs.H08. [,1] [,2] [1,] 0 0 [2,] 0 0 $sPs.H08b [,1] [,2] [1,] 0 0 [2,] 0 0 $AmHi09.06 [,1] [,2] [1,] -1.153e-06 7.160e-06 [2,] 7.160e-06 -1.153e-06 $AmHi09.08 [,1] [,2] [1,] -1.414e-09 1.163e-08 [2,] 1.163e-08 -1.414e-09 $AmHi09.10 [,1] [,2] [1,] -5.702e-13 5.836e-12 [2,] 5.836e-12 -5.702e-13 $AmHi09.12 [,1] [,2] [1,] 1.110e-16 4.441e-16 [2,] 7.772e-16 1.110e-16 $AmHi09.13 [,1] [,2] [1,] 1.11e-16 1.11e-16 [2,] 2.22e-16 1.11e-16 $s.T.s [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 0 $s.T.sO [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 0 $Eigen [,1] [,2] [1,] 0.00e+00 -5.109e-07 [2,] -1.12e-06 0.000e+00 $hybrid [,1] [,2] [1,] 0.000e+00 -5.109e-07 [2,] -2.419e-07 -2.220e-16 > > ## Average number of correct digits [less "extreme" than plot above] > nDig <- sapply(expmList, function(EXPM) -log10(mean(abs(1 - EXPM(me$A)/me$expA)))) > round(nDig, 2) Matr Matr.d Ward s.P.s s.P.sO s.P.sRBS sPs.H08. sPs.H08b AmHi09.06 Inf Inf Inf 14.65 14.65 16.26 Inf Inf 5.38 AmHi09.08 AmHi09.10 AmHi09.12 AmHi09.13 s.T.s s.T.sO Eigen hybrid 8.19 11.49 15.44 15.86 14.65 14.65 6.39 6.73 > ## Ward s.P.s s.P.sO s.T.s s.T.sO Eigen hybrid > ## 16.26 14.65 14.65 14.65 14.65 6.20 6.39 [AMD Opteron 64-bit] > ## Inf 14.65 14.65 14.65 14.65 6.74 6.33 [Pentium-M (32-bit)] > > ###--- rnilMat() : random upper triangular (zero-diagonal) nilpotent n x n matrix > > set.seed(17) > m <- rnilMat(10) > (m. <- as(m,"sparseMatrix"))# for nicer printing - and more *below* 10 x 10 sparse Matrix of class "dtCMatrix" [1,] . 3 10 7 3 4 9 5 9 6 [2,] . . 5 4 3 . 5 6 3 6 [3,] . . . 5 7 7 3 7 5 6 [4,] . . . . 3 7 6 8 2 7 [5,] . . . . . 9 5 2 7 6 [6,] . . . . . . 8 5 4 6 [7,] . . . . . . . 5 5 3 [8,] . . . . . . . . 3 5 [9,] . . . . . . . . . 3 [10,] . . . . . . . . . . > E.m <- expm::expm(m, method="Pade") > as(E.m, "sparseMatrix") 10 x 10 sparse Matrix of class "dtCMatrix" [1,] 1 3 17.5 50.5 110.9 464.7 1224.0 2169.2 3394.8 5938.8 [2,] . 1 5.0 16.5 39.0 172.8 476.1 865.6 1376.4 2438.8 [3,] . . 1.0 5.0 14.5 78.5 251.7 498.6 836.0 1540.0 [4,] . . . 1.0 3.0 20.5 77.5 170.2 299.4 578.0 [5,] . . . . 1.0 9.0 41.0 97.0 180.5 357.1 [6,] . . . . . 1.0 8.0 25.0 51.5 112.3 [7,] . . . . . . 1.0 5.0 12.5 30.5 [8,] . . . . . . . 1.0 3.0 9.5 [9,] . . . . . . . . 1.0 3.0 [10,] . . . . . . . . . 1.0 > > (dN <- 9*7*320) # 20160 [1] 20160 > stopifnot(abs(round(E.m * dN) - (E.m * dN)) < 9e-6) > EmN <- matrix(c(dN, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 3*dN, dN, 0, 0, 0, 0, 0, 0, 0, 0, + 352800, 5*dN, dN, 0, 0, 0, 0, 0, 0, 0, + 1018080, 332640, 5*dN, dN, 0, 0, 0, 0, 0, 0, + 2235240, 786240, 292320, 3*dN, dN, 0, 0, 0, 0, 0, + 9368520, 3483480, 1582560, 413280, 181440, dN, 0, 0, 0, 0, + 24676176, 9598680, 5073600, 1562400, 826560, 161280, dN, 0,0,0, + 43730160, 17451000, 10051440, 3430560, 1955520, 504000, + 5*dN, dN, 0, 0, + 68438436, 27747480, 16853760, 6036240, 3638880, 1038240, + 252000, 3*dN, dN, 0, + 119725855, 49165892, 31046760, 11652480, 7198800, 2264640, + 614880, 191520, 3*dN, dN), + 10, 10) > > Em.xct <- EmN / dN > > stopifnot(all.equal(E.m, Em.xct, + check.attributes = FALSE, tolerance= 1e-13)) > re.x <- sapply(expmList, function(EXPM) relErr(Em.xct, EXPM(m))) Error in solve.default(V) : Lapack routine dgesv: system is exactly singular: U[3,3] = 0 > ## with error message from "safe.Eigen" --> Eigen is NA here > > ## result depends quite a bit on platform > which(is.na(re.x)) # gives c(Eigen = 16L) (but not everywhere ?!?) Eigen 16 > (re.x <- re.x[!is.na(re.x)]) Matr Matr.d Ward s.P.s s.P.sO s.P.sRBS sPs.H08. sPs.H08b AmHi09.06 4.660e-17 4.660e-17 4.660e-17 3.989e-17 3.989e-17 4.424e-17 1.103e-16 1.103e-16 1.509e-11 AmHi09.08 AmHi09.10 AmHi09.12 AmHi09.13 s.T.s s.T.sO hybrid 8.671e-17 8.671e-17 8.671e-17 8.671e-17 8.435e-17 8.435e-17 4.660e-17 > > ## Pentium-M 32-bit ubuntu gave > ## Ward s.P.s s.P.sO sPs.H08. sPs.H08b s.T.s s.T.sO hybrid > ## 1.079e-16 4.505e-14 4.503e-14 9.379e-17 9.379e-17 3.716e-17 7.079e-18 1.079e-16 > ## 32-bit Quad-Core AMD Opteron 2380 (Linux 2.6.30.10-105.2.23.fc11.i686.PAE): > ## Ward s.P.s s.P.sO sPs.H08. sPs.H08b s.T.s s.T.sO hybrid > ## 1.079e-16 4.505e-14 4.503e-14 9.379e-17 9.379e-17 3.716e-17 7.079e-18 1.079e-16 > > ## "Ward77": again more accurate than s+Pade+s, but s+Taylor+s is even more accurate > > ## but on 64-bit AMD Opterons > ## Ward s.P.s s.P.sO sPs.H08. sPs.H08b s.T.s s.T.sO hybrid > ## 4.42e-17 3.99e-17 3.99e-17 1.10e-16 1.10e-16 8.44e-17 8.44e-17 4.42e-17 > ## > ## even more astonishing the result on Mac OSX (x86_32_mac; R-forge, R 2.9.0 patch.) > ## Ward s.P.s s.P.sO sPs.H08. sPs.H08b s.T.s s.T.sO hybrid > ## 5.13e-17 3.99e-17 3.99e-17 1.84e-15 1.84e-15 8.44e-17 8.44e-17 5.13e-17 > > ## 2014-09: AmHi09 are very good (64bit: 8e-17) for p >= 8 (p=6 has 1.5e-11) > > not.09.06 <- which(names(re.x) != "AmHi09.06") > stopifnot(re.x[c("Ward", "s.T.s", "s.T.sO")] < 3e-16, + re.x[["AmHi09.06"]] < 9e-11, # x64 & 686(lnx): = 1.509e-11 + re.x[not.09.06] < 4e-13)# max: 686(32b): 4.52e-14, x64(lnx): 1.103e-16 > > ##-- Looking at *sparse* matrices: [C,Fortran "dense" code based methods will fail]: > (meths <- eval(formals(expm)$method)) [1] "Higham08.b" "Higham08" "AlMohy-Hi09" "Ward77" [5] "PadeRBS" "Pade" "Taylor" "PadeO" [9] "TaylorO" "R_Eigen" "R_Pade" "R_Ward77" [13] "hybrid_Eigen_Ward" > ems <- sapply(meths, function(met) + tryCatch(expm::expm(m., method=met), error=identity)) coercing to dense matrix, as required by method "Higham08.b" coercing to dense matrix, as required by method "AlMohy-Hi09" coercing to dense matrix, as required by method "Ward77" coercing to dense matrix, as required by method "PadeRBS" coercing to dense matrix, as required by method "Pade" coercing to dense matrix, as required by method "Taylor" coercing to dense matrix, as required by method "PadeO" coercing to dense matrix, as required by method "TaylorO" coercing to dense matrix, as required by method "R_Ward77" coercing to dense matrix, as required by method "hybrid_Eigen_Ward" > ok <- !sapply(ems, is, class="error") > meths[ok] # now most... are [1] "Higham08.b" "Higham08" "AlMohy-Hi09" "Ward77" [5] "PadeRBS" "Pade" "Taylor" "PadeO" [9] "TaylorO" "R_Pade" "hybrid_Eigen_Ward" > > > showProc.time() Time (user system elapsed): 0.13 0 0.12 > > proc.time() user system elapsed 2.14 0.09 2.21