R Under development (unstable) (2024-08-17 r87027 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 = 99, 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")# -> nilA3(), facMat(), m2ex3(), .... > 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), gcFirst = 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.63 0.01 0.64 > > ## 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.41 0 0.4 > > ###--- 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 > ## ^^^^^^^^^^^^ > > osV <- abbreviate(sub("\\(.*", "", osVersion), 10) > if(!dev.interactive(TRUE)) pdf(paste0("expm_exact-ex_", osV, ".pdf"), width = 9, height=5) > ## > myRversion <- paste(R.version.string, "--", osVersion) > if((mach <- Sys.info()[["machine"]]) != "x86_64") + myRversion <- paste0(myRversion, "_", mach) > if(!capabilities("long.double")) + myRversion <- paste0(myRversion, "_no_LDbl") > myRversion [1] "R Under development (unstable) (2024-08-17 r87027 ucrt) -- Windows Server 2022 x64 (build 20348)_x86-64" > ## in plots use mtext(myRversion, adj=1, cex=3/4) > > > ## 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.1 0 0.11 > > 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"), + Ward1b = function(x) expm::expm(x, "Ward77", preconditioning = "1bal"), + RWard6 = function(x) expm::expm(x, "R_Ward77", order = 6), + RWard7 = function(x) expm::expm(x, "R_Ward77", order = 7), + RWard8 = function(x) expm::expm(x, "R_Ward77", order = 8), # default + RWard9 = function(x) expm::expm(x, "R_Ward77", order = 9), + s.P.s7 = function(x) expm::expm(x, "Pade", order = 7), + s.POs7 = function(x) expm::expm(x, "PadeO",order = 7), + s.P.s8 = function(x) expm::expm(x, "Pade" ,order = 8), # default + s.POs8 = function(x) expm::expm(x, "PadeO",order = 8), # default + s.P.s9 = function(x) expm::expm(x, "Pade", order = 9), + s.POs9 = function(x) expm::expm(x, "PadeO",order = 9), + s.P.sRBS = function(x) expm::expm(x, "PadeRBS"), + Rs.P.s7 = function(x) expm::expm(x, "R_Pade", order = 7), + Rs.P.s8 = function(x) expm::expm(x, "R_Pade", order = 8), # default + Rs.P.s9 = function(x) expm::expm(x, "R_Pade", order = 9), + 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.07= function(x) expm:::expm.AlMoHi09(x, p = 7), + AmHi09.08= function(x) expm:::expm.AlMoHi09(x, p = 8), # default -- really sub optimal + AmHi09.09= function(x) expm:::expm.AlMoHi09(x, p = 9), + AmHi09.10= function(x) expm:::expm.AlMoHi09(x, p = 10), + AmHi09.11= function(x) expm:::expm.AlMoHi09(x, p = 11), + AmHi09.12= function(x) expm:::expm.AlMoHi09(x, p = 12), + AmHi09.13= function(x) expm:::expm.AlMoHi09(x, p = 13), + s.T.s7 = function(x) expm::expm(x, "Taylor", order = 7), + s.TOs7 = function(x) expm::expm(x, "TaylorO",order = 7), + s.T.s8 = function(x) expm::expm(x, "Taylor", order = 8), # default + s.TOs8 = function(x) expm::expm(x, "TaylorO",order = 8), # default + s.T.s9 = function(x) expm::expm(x, "Taylor", order = 9), + s.TOs9 = function(x) expm::expm(x, "TaylorO",order = 9), + Eigen = expm.safe.Eigen, # "R_Eigen" + hybrid = function(x) expm::expm(x, "hybrid") + ) > > > ## set.seed(12) > ## facMchk <- replicate(if(doExtras) 100 else 20, facMat(20, rnorm)) > set.seed(12) > fRE <- replicate(if(doExtras) 100 else 20, + re.facMat(20, expmList)) # if(doExtras) gives one "No Matrix found ..." warning > nDig <- -log10(t(fRE["relErr",,])) > cat("Number of correct decimal digits for facMat(20, rnorm):\n") Number of correct decimal digits for facMat(20, rnorm): > t(apply(nDig, 2, summary)) Min. 1st Qu. Median Mean 3rd Qu. Max. Matr 13.681 14.19 14.29 14.32 14.53 14.73 Matr.d 13.681 14.19 14.29 14.32 14.53 14.73 Ward 13.681 14.19 14.29 14.32 14.53 14.73 Ward1b 13.681 14.19 14.29 14.32 14.53 14.73 RWard6 13.243 13.82 13.98 13.95 14.09 14.32 RWard7 13.629 13.82 13.96 13.96 14.11 14.32 RWard8 13.333 13.80 13.91 13.91 14.11 14.22 RWard9 13.521 13.74 13.82 13.90 14.12 14.23 s.P.s7 12.582 13.05 13.24 13.22 13.42 13.50 s.POs7 12.582 13.05 13.24 13.22 13.42 13.50 s.P.s8 12.825 13.15 13.24 13.25 13.32 13.63 s.POs8 12.825 13.15 13.24 13.25 13.32 13.63 s.P.s9 12.764 13.15 13.30 13.23 13.37 13.50 s.POs9 12.764 13.15 13.30 13.23 13.37 13.50 s.P.sRBS 13.918 14.29 14.45 14.42 14.56 14.77 Rs.P.s7 13.452 13.79 13.88 13.91 14.06 14.28 Rs.P.s8 13.359 13.80 13.91 13.92 14.05 14.31 Rs.P.s9 13.364 13.73 13.90 13.90 14.04 14.35 sPs.H08. 14.379 14.81 14.92 14.89 14.99 15.23 sPs.H08b 14.379 14.83 14.94 14.92 15.04 15.23 AmHi09.07 8.618 10.09 10.54 11.03 12.36 13.39 AmHi09.08 9.782 11.27 11.72 12.19 13.56 14.40 AmHi09.09 12.000 13.89 14.24 14.07 14.63 15.01 AmHi09.10 13.434 14.57 14.71 14.64 14.84 15.01 AmHi09.11 14.512 14.67 14.77 14.81 14.91 15.15 AmHi09.12 14.512 14.67 14.77 14.81 14.91 15.17 AmHi09.13 14.512 14.67 14.77 14.81 14.91 15.17 s.T.s7 13.020 13.43 13.51 13.53 13.71 13.90 s.TOs7 13.020 13.43 13.51 13.53 13.71 13.90 s.T.s8 13.020 13.43 13.51 13.53 13.71 13.90 s.TOs8 13.020 13.43 13.51 13.53 13.71 13.90 s.T.s9 13.020 13.43 13.51 13.53 13.71 13.90 s.TOs9 13.020 13.43 13.51 13.53 13.71 13.90 Eigen 14.206 14.44 14.49 14.50 14.57 14.76 hybrid 14.197 14.45 14.52 14.50 14.58 14.77 > > ## Now look at that: > eaxis <- if(requireNamespace("sfsmisc")) sfsmisc::eaxis else axis Loading required namespace: sfsmisc > op <- par(mar=.1 + c(5,4 + 1.5, 4,2)) > boxplot(t(fRE["relErr",,]), log="x", xaxt="n", + notch=TRUE, ylim = c(8e-16, 4e-9), + horizontal=TRUE, las = 1, + 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 > eaxis(1); grid(lty = 3); mtext(myRversion, adj=1, cex=3/4) > par(op) > > showProc.time() Time (user system elapsed): 0.52 0.05 0.57 > > if(doExtras) withAutoprint({ # also "large" n = 100 ------------------------------------------ + str(rf100 <- replicate(20, re.facMat(100, expmList))) + 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.POs 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.TOs 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 + + nDig <- -log10(t(rf100["relErr",,])) + cat("Number of correct decimal digits for facMat(100, rnorm):\n") + t(apply(nDig, 2, summary)) + + ##--> take out the real slow ones for the subsequent tests: + (not.slow <- grep("^s\\.[PT]", names(expmList), invert = TRUE, value = TRUE)) + + set.seed(18) ## 12 replicates is too small .. but then it's too slow otherwise: + rf400 <- replicate(12, re.facMat(400, expmList[not.slow])) + + showProc.time() + 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 + + nDig <- -log10(t(rf400["relErr",,])) + cat("Number of correct decimal digits for (12 rep. of) facMat(400, rnorm):\n") + t(apply(nDig, 2, summary)) + + }) else { # *not* (doExtras) ----------------------------------------------------------------- + ## times (in millisec): + print(1000*t(apply(fRE["c.time",,], 1, summary))) + } Min. 1st Qu. Median Mean 3rd Qu. Max. Matr 0 0 0 0.0 0 0 Matr.d 0 0 0 0.0 0 0 Ward 0 0 0 0.0 0 0 Ward1b 0 0 0 0.0 0 0 RWard6 0 0 0 1.5 0 20 RWard7 0 0 0 0.5 0 10 RWard8 0 0 0 0.5 0 10 RWard9 0 0 0 1.0 0 10 s.P.s7 0 0 0 1.0 0 20 s.POs7 0 0 0 0.5 0 10 s.P.s8 0 0 0 0.0 0 0 s.POs8 0 0 0 2.0 0 20 s.P.s9 0 0 0 0.0 0 0 s.POs9 0 0 0 1.0 0 20 s.P.sRBS 0 0 0 0.0 0 0 Rs.P.s7 0 0 0 0.0 0 0 Rs.P.s8 0 0 0 1.0 0 20 Rs.P.s9 0 0 0 0.0 0 0 sPs.H08. 0 0 0 0.5 0 10 sPs.H08b 0 0 0 0.0 0 0 AmHi09.07 0 0 0 0.5 0 10 AmHi09.08 0 0 0 0.0 0 0 AmHi09.09 0 0 0 0.0 0 0 AmHi09.10 0 0 0 0.0 0 0 AmHi09.11 0 0 0 0.0 0 0 AmHi09.12 0 0 0 0.0 0 0 AmHi09.13 0 0 0 0.0 0 0 s.T.s7 0 0 0 0.0 0 0 s.TOs7 0 0 0 1.0 0 20 s.T.s8 0 0 0 0.0 0 0 s.TOs8 0 0 0 0.5 0 10 s.T.s9 0 0 0 0.0 0 0 s.TOs9 0 0 0 1.0 0 20 Eigen 0 0 0 1.0 0 20 hybrid 0 0 0 1.0 0 20 > > ## Now try an example with badly conditioned "random" M matrix... > ## ... > ## ... (not yet -- TODO?) > > ### 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 RWard6 RWard7 RWard8 RWard9 Matr Matr.d Ward Ward1b s.P.sRBS 0.000e+00 0.000e+00 0.000e+00 0.000e+00 5.584e-33 5.584e-33 5.584e-33 5.584e-33 5.584e-33 Rs.P.s9 sPs.H08. sPs.H08b AmHi09.13 AmHi09.12 Rs.P.s8 Rs.P.s7 s.P.s8 s.POs8 5.584e-33 5.584e-33 5.584e-33 1.509e-16 2.515e-16 7.545e-16 8.048e-16 2.113e-15 2.113e-15 s.P.s9 s.POs9 s.T.s8 s.TOs8 s.T.s9 s.TOs9 s.P.s7 s.POs7 s.T.s7 2.113e-15 2.113e-15 2.113e-15 2.113e-15 2.113e-15 2.113e-15 2.213e-15 2.213e-15 2.616e-15 s.TOs7 AmHi09.11 AmHi09.10 AmHi09.09 Eigen hybrid AmHi09.08 AmHi09.07 2.616e-15 1.373e-14 2.325e-12 6.032e-11 1.836e-09 1.836e-09 4.819e-09 6.717e-08 > > 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.22 0.02 0.23 > > cbind(sort(apply(log(t.m2),2, median, na.rm=TRUE))) [,1] RWard6 -Inf RWard7 -Inf RWard8 -Inf RWard9 -Inf sPs.H08. -86.05 sPs.H08b -86.05 Rs.P.s9 -85.70 Matr -83.62 Matr.d -83.62 Ward -83.62 Ward1b -83.62 s.P.sRBS -81.20 AmHi09.13 -36.43 AmHi09.12 -35.92 Rs.P.s8 -34.82 Rs.P.s7 -34.76 s.P.s8 -33.79 s.POs8 -33.79 s.P.s9 -33.79 s.POs9 -33.79 s.T.s8 -33.79 s.TOs8 -33.79 s.T.s9 -33.79 s.TOs9 -33.79 s.P.s7 -33.74 s.POs7 -33.74 s.T.s7 -33.58 s.TOs7 -33.58 AmHi09.11 -31.92 AmHi09.10 -26.79 hybrid -23.74 AmHi09.09 -23.53 Eigen -20.12 AmHi09.08 -19.15 AmHi09.07 -16.52 > ## '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] RWard8 3.31 RWard6 3.53 RWard7 3.61 RWard9 3.78 sPs.H08. 6.89 sPs.H08b 6.89 s.P.sRBS 8.36 Matr 10.20 Matr.d 10.20 Ward 10.20 Ward1b 10.20 AmHi09.13 10.90 Rs.P.s9 11.00 AmHi09.12 13.30 Rs.P.s8 14.70 Rs.P.s7 15.90 s.T.s9 19.30 s.TOs9 19.30 s.T.s8 19.70 s.TOs8 19.70 s.P.s9 21.70 s.POs9 21.70 s.P.s8 22.60 s.POs8 22.60 s.P.s7 25.40 s.POs7 25.40 s.T.s7 25.90 s.TOs7 25.90 hybrid 27.00 AmHi09.11 29.40 AmHi09.10 30.70 Eigen 31.80 AmHi09.09 31.80 AmHi09.08 33.10 AmHi09.07 34.20 > ## 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.TOs 8.33 > ## s.P.s 9.11 > ## s.POs 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) RWard8 RWard6 RWard7 RWard9 sPs.H08. sPs.H08b s.P.sRBS Matr [1,] 1.00e-16 1.01e-16 1.51e-16 2.01e-16 9.99e-17 9.99e-17 2.00e-16 6.02e-16 [2,] 5.03e-17 2.01e-16 5.03e-17 5.03e-17 1.01e-16 1.01e-16 2.01e-16 5.53e-16 [3,] 5.03e-17 0.00e+00 1.01e-16 5.03e-17 5.03e-17 5.03e-17 2.51e-16 3.02e-16 [4,] 7.50e-25 0.00e+00 0.00e+00 7.50e-25 2.51e-16 2.51e-16 3.75e-25 1.51e-16 [5,] 5.86e-27 0.00e+00 1.17e-26 1.17e-26 3.52e-16 3.52e-16 2.01e-16 1.51e-16 [6,] 4.57e-29 5.03e-17 5.03e-17 0.00e+00 2.01e-16 2.01e-16 3.52e-16 3.52e-16 [7,] 5.03e-17 1.01e-16 1.01e-16 5.03e-17 3.52e-16 3.52e-16 3.52e-16 6.04e-16 [8,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 5.58e-33 5.58e-33 5.58e-33 5.58e-33 [9,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 4.36e-35 4.36e-35 4.36e-35 4.36e-35 [10,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 6.82e-37 0.00e+00 [11,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 5.33e-39 5.33e-39 [12,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 4.16e-41 4.16e-41 0.00e+00 4.16e-41 [13,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 6.50e-43 6.50e-43 [14,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 5.08e-45 5.08e-45 [15,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 3.97e-47 3.97e-47 7.94e-47 7.94e-47 [16,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 6.20e-49 6.20e-49 [17,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 4.84e-51 0.00e+00 [18,] 0.00e+00 0.00e+00 0.00e+00 0.00e+00 3.78e-53 3.78e-53 0.00e+00 3.78e-53 > ## ==> 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) : 74 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(requireNamespace("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) + } > > ## 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 $Ward1b [,1] [,2] [1,] 0 0 [2,] 0 0 $RWard6 [,1] [,2] [1,] 0 0 [2,] 0 0 $RWard7 [,1] [,2] [1,] 0 0 [2,] 0 0 $RWard8 [,1] [,2] [1,] 0 0 [2,] 0 0 $RWard9 [,1] [,2] [1,] 0 0 [2,] 0 0 $s.P.s7 [,1] [,2] [1,] -2.22e-15 -2.442e-15 [2,] -2.22e-15 -2.220e-15 attr(,"accuracy") [1] 2.998e-15 $s.POs7 [,1] [,2] [1,] -2.22e-15 -2.442e-15 [2,] -2.22e-15 -2.220e-15 attr(,"accuracy") [1] 2.998e-15 $s.P.s8 [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 1.665e-16 $s.POs8 [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 1.665e-16 $s.P.s9 [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 4.19e-15 $s.POs9 [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 4.19e-15 $s.P.sRBS [,1] [,2] [1,] 0.00e+00 0 [2,] -2.22e-16 0 $Rs.P.s7 [,1] [,2] [1,] -6.661e-16 -8.882e-16 [2,] -8.882e-16 -6.661e-16 $Rs.P.s8 [,1] [,2] [1,] -6.661e-16 -6.661e-16 [2,] -8.882e-16 -6.661e-16 $Rs.P.s9 [,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.07 [,1] [,2] [1,] 1.999e-08 -1.615e-07 [2,] -1.615e-07 1.999e-08 $AmHi09.08 [,1] [,2] [1,] -1.414e-09 1.163e-08 [2,] 1.163e-08 -1.414e-09 $AmHi09.09 [,1] [,2] [1,] 1.497e-11 -1.510e-10 [2,] -1.510e-10 1.497e-11 $AmHi09.10 [,1] [,2] [1,] -5.702e-13 5.836e-12 [2,] 5.836e-12 -5.702e-13 $AmHi09.11 [,1] [,2] [1,] 2.998e-15 -3.508e-14 [2,] -3.486e-14 2.998e-15 $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.s7 [,1] [,2] [1,] 1.776e-15 -4.219e-15 [2,] -4.219e-15 1.776e-15 attr(,"accuracy") [1] 2.22e-15 $s.TOs7 [,1] [,2] [1,] 1.776e-15 -4.219e-15 [2,] -4.219e-15 1.776e-15 attr(,"accuracy") [1] 2.22e-15 $s.T.s8 [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 0 $s.TOs8 [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 0 $s.T.s9 [,1] [,2] [1,] -2.22e-15 -2.22e-15 [2,] -2.22e-15 -2.22e-15 attr(,"accuracy") [1] 0 $s.TOs9 [,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 Ward1b RWard6 RWard7 RWard8 RWard9 s.P.s7 Inf Inf Inf Inf Inf Inf Inf Inf 14.64 s.POs7 s.P.s8 s.POs8 s.P.s9 s.POs9 s.P.sRBS Rs.P.s7 Rs.P.s8 Rs.P.s9 14.64 14.65 14.65 14.65 14.65 16.26 15.11 15.14 16.26 sPs.H08. sPs.H08b AmHi09.07 AmHi09.08 AmHi09.09 AmHi09.10 AmHi09.11 AmHi09.12 AmHi09.13 Inf Inf 7.04 8.19 10.08 11.49 13.72 15.44 15.86 s.T.s7 s.TOs7 s.T.s8 s.TOs8 s.T.s9 s.TOs9 Eigen hybrid 14.52 14.52 14.65 14.65 14.65 14.65 6.39 6.73 > ## Ward s.P.s s.POs s.T.s s.TOs 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)] > > showProc.time() Time (user system elapsed): 0.01 0.04 0.06 > > ###--- 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 > all.equal(E.m, Em.xct, check.attributes = FALSE, tolerance = 0) [1] "Mean relative difference: 1.549e-16" > 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 34 > (re.x <- re.x[!is.na(re.x)]) Matr Matr.d Ward Ward1b RWard6 RWard7 RWard8 RWard9 s.P.s7 4.660e-17 4.660e-17 4.660e-17 4.660e-17 8.435e-17 7.020e-17 5.368e-17 8.435e-17 7.145e-17 s.POs7 s.P.s8 s.POs8 s.P.s9 s.POs9 s.P.sRBS Rs.P.s7 Rs.P.s8 Rs.P.s9 7.145e-17 3.989e-17 3.989e-17 1.638e-16 1.638e-16 4.424e-17 7.020e-17 5.368e-17 8.435e-17 sPs.H08. sPs.H08b AmHi09.07 AmHi09.08 AmHi09.09 AmHi09.10 AmHi09.11 AmHi09.12 AmHi09.13 1.103e-16 1.103e-16 2.377e-16 8.671e-17 8.671e-17 8.671e-17 8.671e-17 8.671e-17 8.671e-17 s.T.s7 s.TOs7 s.T.s8 s.TOs8 s.T.s9 s.TOs9 hybrid 8.435e-17 8.435e-17 8.435e-17 8.435e-17 8.435e-17 8.435e-17 4.660e-17 > > ## Pentium-M 32-bit ubuntu gave > ## Ward s.P.s s.POs sPs.H08. sPs.H08b s.T.s s.TOs 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.POs sPs.H08. sPs.H08b s.T.s s.TOs 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.POs sPs.H08. sPs.H08b s.T.s s.TOs 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.POs sPs.H08. sPs.H08b s.T.s s.TOs 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.s8", "s.TOs8")] < 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, class2="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.08 0 0.08 > > ## Complex Matrices > re.facMat.Z <- function(n, EXPMlist, rFUN = function(n) rnorm(n) + 1i*rnorm(n), ...) + { + stopifnot(is.list(EXPMlist)) + r <- facMat(n, rFUN, ...) + vapply(EXPMlist, function(EXPM) { + ct <- system.time(E <- EXPM(r$A), gcFirst = FALSE)[[1]] + c(relErr = relErr(r$expA, E), c.time = ct) + }, double(2)) + } > > (nmL <- names(expmList)) [1] "Matr" "Matr.d" "Ward" "Ward1b" "RWard6" "RWard7" "RWard8" [8] "RWard9" "s.P.s7" "s.POs7" "s.P.s8" "s.POs8" "s.P.s9" "s.POs9" [15] "s.P.sRBS" "Rs.P.s7" "Rs.P.s8" "Rs.P.s9" "sPs.H08." "sPs.H08b" "AmHi09.07" [22] "AmHi09.08" "AmHi09.09" "AmHi09.10" "AmHi09.11" "AmHi09.12" "AmHi09.13" "s.T.s7" [29] "s.TOs7" "s.T.s8" "s.TOs8" "s.T.s9" "s.TOs9" "Eigen" "hybrid" > ## [1] "Matr" "Matr.d" "Ward" "Ward1b" "s.P.s" "s.POs" "s.P.s7" > ## [8] "s.POs7" "s.P.s9" "s.POs9" "s.P.sRBS" "sPs.H08." "sPs.H08b" "AmHi09.06" > ## [15] "AmHi09.07" "AmHi09.08" "AmHi09.09" "AmHi09.10" "AmHi09.12" "AmHi09.13" "s.T.s" > ## [22] "s.TOs" "s.T.s7" "s.TOs7" "s.T.s9" "s.TOs9" "Eigen" "hybrid" > > ## dropping "Matr", "Matr.d" (which gives "dgeMatrix" currently --> mean(.) fails ... > ## "Ward" "Ward1b" and "hybrid" error "not a numeric Matrix" > ## "AmHi09": C code currently only for double precision ((FIXME)) > (cplxN <- grep("^(Matr|Ward|hybr|AmHi09|s\\.[PT])", nmL, invert = TRUE, value = TRUE)) [1] "RWard6" "RWard7" "RWard8" "RWard9" "Rs.P.s7" "Rs.P.s8" "Rs.P.s9" "sPs.H08." [9] "sPs.H08b" "Eigen" > > rr <- re.facMat.Z(4, expmList[cplxN]) > > set.seed(47) > fREZ <- replicate(if(doExtras) 64 else 12, + re.facMat.Z(15, expmList[cplxN])) > nDig <- -log10(t(fREZ["relErr",,])) > cat("Number of correct decimal digits for facMat(20, rnorm + i*rnorm):\n") Number of correct decimal digits for facMat(20, rnorm + i*rnorm): > t(apply(nDig, 2, summary)) Min. 1st Qu. Median Mean 3rd Qu. Max. RWard6 13.45 13.75 13.86 13.85 13.96 14.17 RWard7 13.45 13.74 13.79 13.80 13.89 14.18 RWard8 13.56 13.71 13.80 13.81 13.91 14.07 RWard9 13.57 13.75 13.87 13.85 13.93 14.21 Rs.P.s7 13.57 13.75 13.86 13.86 13.93 14.16 Rs.P.s8 13.59 13.78 13.85 13.84 13.92 14.06 Rs.P.s9 13.56 13.72 13.76 13.80 13.95 14.00 sPs.H08. 14.67 14.82 14.85 14.85 14.89 15.08 sPs.H08b 14.67 14.85 14.86 14.87 14.89 15.08 Eigen 14.18 14.36 14.37 14.38 14.39 14.55 > > ## times (in millisec): > print(1000*t(apply(fREZ["c.time",,], 1, summary))) Min. 1st Qu. Median Mean 3rd Qu. Max. RWard6 0 0 0 0.0000 0 0 RWard7 0 0 0 0.0000 0 0 RWard8 0 0 0 0.0000 0 0 RWard9 0 0 0 0.8333 0 10 Rs.P.s7 0 0 0 0.0000 0 0 Rs.P.s8 0 0 0 0.0000 0 0 Rs.P.s9 0 0 0 0.0000 0 0 sPs.H08. 0 0 0 0.0000 0 0 sPs.H08b 0 0 0 1.6667 0 20 Eigen 0 0 0 1.6667 0 20 > > ## Now look at that: > op <- par(mar=.1 + c(5,4 + 1.5, 4,2)) > boxplot(t(fREZ["relErr",,]), log="x", xaxt="n", + notch=TRUE, # ylim = c(8e-16, 4e-9), + horizontal=TRUE, las = 1, + 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 > eaxis(1); grid(lty = 3); mtext(myRversion, adj=1, cex=3/4) > par(op) > > > > showProc.time() Time (user system elapsed): 0.13 0 0.13 > > proc.time() user system elapsed 3.31 0.32 3.61