Nothing
tests/exact-ex.R
In expm: Matrix Exponential, Log, 'etc'
#### Examples where we know the result "exactly"
library(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
mSource("demo", "exact-fn.R")# -> nilA3(), facMat(), m2ex3(), ....
doExtras
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))
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()
## 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)))
epsC <- .Machine$double.eps
cat("relErr(expm(.,Pade)) - relErr(expm(.,Ward77)) in Machine_eps units:\n")
print(summary(c(re - RE)) / epsC)
## 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()
###--- 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(gsub("[^[:alnum:]]", '', sub("\\(.*", '', osVersion)), 12))
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
## 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()
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
## 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")
t(apply(nDig, 2, summary))
## Now look at that:
eaxis <- if(requireNamespace("sfsmisc")) sfsmisc::eaxis else axis
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")
eaxis(1); grid(lty = 3); mtext(myRversion, adj=1, cex=3/4)
par(op)
showProc.time()
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)))
}
## 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
eps <- 10^-(1:18)
t.m2 <- t(sapply(eps, re.m2ex3, EXPMlist = expmList))
## --> 3 error messages from solve(V), 5 error messages from try(. "R_Eigen" ...)
showProc.time()
cbind(sort(apply(log(t.m2),2, median, na.rm=TRUE)))
## '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))
## 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)
## ==> 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")
}
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}), "]")))
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))
me$expA * exp(1) ## the correct value ; numerically identical to simple matrix:
## 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)
## 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)
## 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()
###--- 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*
E.m <- expm::expm(m, method="Pade")
as(E.m, "sparseMatrix")
(dN <- 9*7*320) # 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)
stopifnot(all.equal(E.m, Em.xct, check.attributes = FALSE, tolerance= 1e-13))
re.x <- sapply(expmList, function(EXPM) relErr(Em.xct, EXPM(m)))
## 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 ?!?)
(re.x <- re.x[!is.na(re.x)])
## 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))
ems <- sapply(meths, function(met)
tryCatch(expm::expm(m., method=met), error=identity))
ok <- !sapply(ems, is, class2="error")
meths[ok] # now most... are
showProc.time()
## 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" "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))
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")
t(apply(nDig, 2, summary))
## times (in millisec):
print(1000*t(apply(fREZ["c.time",,], 1, summary)))
## 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' <complex> eigen-ok 20 x 20 matrix")
eaxis(1); grid(lty = 3); mtext(myRversion, adj=1, cex=3/4)
par(op)
showProc.time()
Try the expm package in your browser
Any scripts or data that you put into this service are public.
expm documentation built on Sept. 11, 2024, 8:58 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
#### Examples where we know the result "exactly"
library(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
mSource("demo", "exact-fn.R")# -> nilA3(), facMat(), m2ex3(), ....
doExtras
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))
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()
## 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)))
epsC <- .Machine$double.eps
cat("relErr(expm(.,Pade)) - relErr(expm(.,Ward77)) in Machine_eps units:\n")
print(summary(c(re - RE)) / epsC)
## 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()
###--- 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(gsub("[^[:alnum:]]", '', sub("\\(.*", '', osVersion)), 12))
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
## 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()
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
## 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")
t(apply(nDig, 2, summary))
## Now look at that:
eaxis <- if(requireNamespace("sfsmisc")) sfsmisc::eaxis else axis
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")
eaxis(1); grid(lty = 3); mtext(myRversion, adj=1, cex=3/4)
par(op)
showProc.time()
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)))
}
## 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
eps <- 10^-(1:18)
t.m2 <- t(sapply(eps, re.m2ex3, EXPMlist = expmList))
## --> 3 error messages from solve(V), 5 error messages from try(. "R_Eigen" ...)
showProc.time()
cbind(sort(apply(log(t.m2),2, median, na.rm=TRUE)))
## '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))
## 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)
## ==> 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")
}
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}), "]")))
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))
me$expA * exp(1) ## the correct value ; numerically identical to simple matrix:
## 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)
## 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)
## 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()
###--- 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*
E.m <- expm::expm(m, method="Pade")
as(E.m, "sparseMatrix")
(dN <- 9*7*320) # 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)
stopifnot(all.equal(E.m, Em.xct, check.attributes = FALSE, tolerance= 1e-13))
re.x <- sapply(expmList, function(EXPM) relErr(Em.xct, EXPM(m)))
## 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 ?!?)
(re.x <- re.x[!is.na(re.x)])
## 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))
ems <- sapply(meths, function(met)
tryCatch(expm::expm(m., method=met), error=identity))
ok <- !sapply(ems, is, class2="error")
meths[ok] # now most... are
showProc.time()
## 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" "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))
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")
t(apply(nDig, 2, summary))
## times (in millisec):
print(1000*t(apply(fREZ["c.time",,], 1, summary)))
## 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' <complex> eigen-ok 20 x 20 matrix")
eaxis(1); grid(lty = 3); mtext(myRversion, adj=1, cex=3/4)
par(op)
showProc.time()
Try the expm package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.