R Under development (unstable) (2023-11-28 r85645 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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. > library(expm) Loading required package: Matrix Attaching package: 'expm' The following object is masked from 'package:Matrix': expm > > source(system.file("test-tools.R", package= "expm"), keep.source=FALSE) Loading required package: tools > > ## Note that these results are achieved with the default > ## settings order=8, method="Pade" -- accuracy could > ## presumably be improved still further by some tuning > ## of these settings. > > ### Latest ATLAS (for BDR on F 36; R-devel Oct.2023) has much worse precision: > ##==> use much larger tolerance in such cases: > > ## Simplified (needs R 3.4.0 and newer, from robustbase/inst/xtraR/platform-sessionInfo.R ) : > BLAS <- extSoftVersion()[["BLAS"]] > Lapack <- La_library() > is.BLAS.Lapack <- identical(BLAS, Lapack) > ## A cheap check (that works on KH's debian-gcc setup, 2019-05): > if(!length(BLAS.is.openBLAS <- grepl("openblas", BLAS, ignore.case=TRUE))) + BLAS.is.openBLAS <- NA > if(!length(Lapack.is.openBLAS <- grepl("openblas", Lapack, ignore.case=TRUE))) + Lapack.is.openBLAS <- NA > > (maybeATLAS <- is.BLAS.Lapack && !BLAS.is.openBLAS) [1] TRUE > > > ## ---------------------------- > ## Test case 1 from Ward (1977) > ## ---------------------------- > T1 <- rbind(c(4, 2, 0), + c(1, 4, 1), + c(1, 1, 4)) > (m1 <- expm(T1, method="Pade")) [,1] [,2] [,3] [1,] 147.8666 183.7651 71.79703 [2,] 127.7811 183.7651 91.88257 [3,] 127.7811 163.6796 111.96811 attr(,"accuracy") [1] 9.942847e-12 > (m1O <- expm(T1, method="PadeO"))# very slightly different [,1] [,2] [,3] [1,] 147.8666 183.7651 71.79703 [2,] 127.7811 183.7651 91.88257 [3,] 127.7811 163.6796 111.96811 attr(,"accuracy") [1] 9.942847e-12 > (m1T <- expm(T1, method="Taylor")) [,1] [,2] [,3] [1,] 147.8666 183.7651 71.79703 [2,] 127.7811 183.7651 91.88257 [3,] 127.7811 163.6796 111.96811 attr(,"accuracy") [1] 3.907339e-13 > (m1TO <- expm(T1, method="TaylorO")) [,1] [,2] [,3] [1,] 147.8666 183.7651 71.79703 [2,] 127.7811 183.7651 91.88257 [3,] 127.7811 163.6796 111.96811 attr(,"accuracy") [1] 3.907339e-13 > ## True Eigenvalue Decomposition of T1 > s2 <- sqrt(2) > eV1 <- matrix(c(s2,s2,s2, -2,1,1, 2,-1,-1) / sqrt(6), + 3,3) > L1 <- diag(lm1 <- c(6, 3, 3)) > stopifnot( + all.equal(eV1 %*% L1, T1 %*% eV1, tolerance=1e-15) + ) > ## However, eV1 is not orthogonal, but of rank 2 > > if(FALSE) { ## require("Rmpfr")) { ## 200 bit precision version of that + S2 <- sqrt(mpfr(2,200)) + E1 <- c(S2,S2,S2, -2,1,1, 2,-1,-1) / sqrt(mpfr(6,200)) + dim(E1) <- c(3,3) + print(E1 %*% L1) + print(E1) + } > > ## "true" result > m1.t <- matrix(c(147.866622446369, 127.781085523181, 127.781085523182, + 183.765138646367, 183.765138646366, 163.679601723179, + 71.797032399996, 91.8825693231832, 111.968106246371), 3,3) > stopifnot(all.equal(m1.t, m1, check.attributes=FALSE, tolerance = 1e-13), + all.equal(m1.t, m1O, check.attributes=FALSE, tolerance = 1e-13), + all.equal(m1.t,m1T, check.attributes=FALSE, tolerance = 1e-13), + all.equal(m1.t,m1TO, check.attributes=FALSE, tolerance = 1e-13), + all.equal(m1.t, expm(T1,"Ward77"), tolerance = 1e-13), + all.equal(m1.t, expm(T1,"R_Pade"), tolerance = 1e-13), + all.equal(m1.t, expm(T1,"R_Ward77"), tolerance = 1e-13)) > ## -- these agree with ward (1977, p608) > ## > m1.2 <- try( expm(T1, "R_Eigen") ) ## 32-bit: gives an error from solve; 64-bit "ok" > if(!inherits(m1.2, "try-error")) { + if(FALSE)## with libatlas R_Eigen is "sehr eigen" + stopifnot(all.equal(m1.t, m1.2, check.attributes=FALSE)) + ## but it's less accurate: + print( all.equal(m1.t, m1.2, check.attributes=FALSE, tolerance= 1e-12)) + ##-> rel.diff = 6.44e-10 / 6.2023e-10 + ##__ BUT 0.1228099 + ##__ with libatlas (ubuntu 12.04 libatlas-base-dev Version: 3.8.4-3build1) + } [1] "Mean relative difference: 6.202299e-10" > > ## > ## ---------------------------- > ## Test case 2 from Ward (1977) > ## ---------------------------- > T2 <- t(matrix(c( + 29.87942128909879, .7815750847907159, -2.289519314033932, + .7815750847907159, 25.72656945571064, 8.680737820540137, + -2.289519314033932, 8.680737820540137, 34.39400925519054), + 3, 3)) > (m2 <- expm(T2, method="Pade")) [,1] [,2] [,3] [1,] 5.496314e+15 -1.823188e+16 -3.047577e+16 [2,] -1.823188e+16 6.060523e+16 1.012918e+17 [3,] -3.047577e+16 1.012918e+17 1.692944e+17 attr(,"accuracy") [1] 53010 > ## [,1] [,2] [,3] > ##[1,] 5496313853692357 -18231880972009844 -30475770808580828 > ##[2,] -18231880972009852 60605228702227024 101291842930256144 > ##[3,] -30475770808580840 101291842930256144 169294411240859072 > ## -- which agrees with Ward (1977) to 13 significant figures > (m2O <- expm(T2, method="PadeO")) [,1] [,2] [,3] [1,] 5.496314e+15 -1.823188e+16 -3.047577e+16 [2,] -1.823188e+16 6.060523e+16 1.012918e+17 [3,] -3.047577e+16 1.012918e+17 1.692944e+17 attr(,"accuracy") [1] 53010 > (m2T <- expm(T2,method="Taylor")) [,1] [,2] [,3] [1,] 5.496314e+15 -1.823188e+16 -3.047577e+16 [2,] -1.823188e+16 6.060523e+16 1.012918e+17 [3,] -3.047577e+16 1.012918e+17 1.692944e+17 attr(,"accuracy") [1] 55.1362 > (m2TO <- expm(T2,method="TaylorO")) [,1] [,2] [,3] [1,] 5.496314e+15 -1.823188e+16 -3.047577e+16 [2,] -1.823188e+16 6.060523e+16 1.012918e+17 [3,] -3.047577e+16 1.012918e+17 1.692944e+17 attr(,"accuracy") [1] 55.1362 > > m2.t <- matrix(c(5496313853692216, -18231880972008932, -30475770808579672, + -18231880972008928, 60605228702222480, 101291842930249776, + -30475770808579672, 101291842930249808, 169294411240850528), + 3, 3) > > ## -- in this case a very similar degree of accuracy -- even Taylor is ok > stopifnot(all.equal(m2.t, m2, check.attributes=FALSE, tolerance = 1e-12), + all.equal(m2.t, m2O,check.attributes=FALSE, tolerance = 1e-12), + all.equal(m2.t,m2T, check.attributes=FALSE, tolerance = 1e-12), + all.equal(m2.t,m2TO,check.attributes=FALSE, tolerance = 1e-12), + all.equal(m2.t, expm(T2,"Ward77"), tolerance = 1e-12), + all.equal(m2.t, expm(T2,"R_Ward77"), tolerance = 1e-12), + all.equal(m2.t, expm(T2,"R_Pade"), tolerance = 1e-12), + TRUE) > > ## ---------------------------- > ## Test case 3 from Ward (1977) > ## ---------------------------- > T3 <- t(matrix(c( + -131, 19, 18, + -390, 56, 54, + -387, 57, 52), 3, 3)) > (m3 <- expm(T3, method="Pade")) [,1] [,2] [,3] [1,] -1.509644 0.3678794 0.1353353 [2,] -5.632571 1.4715178 0.4060058 [3,] -4.934938 1.1036383 0.5413411 attr(,"accuracy") [1] 4.375779e-11 > ## [,1] [,2] [,3] > ##[1,] -1.5096441587713636 0.36787943910439874 0.13533528117301735 > ##[2,] -5.6325707997970271 1.47151775847745725 0.40600584351567010 > ##[3,] -4.9349383260294299 1.10363831731417195 0.54134112675653534 > ## -- agrees to 10dp with Ward (1977), p608. > (m3O <- expm(T3, method="PadeO")) [,1] [,2] [,3] [1,] -1.509644 0.3678794 0.1353353 [2,] -5.632571 1.4715178 0.4060058 [3,] -4.934938 1.1036383 0.5413411 attr(,"accuracy") [1] 4.375779e-11 > (m3T <- expm(T3,method="Taylor")) [,1] [,2] [,3] [1,] -1.509644 0.3678794 0.1353353 [2,] -5.632571 1.4715178 0.4060058 [3,] -4.934938 1.1036383 0.5413411 attr(,"accuracy") [1] 0 > (m3TO <- expm(T3,method="TaylorO")) [,1] [,2] [,3] [1,] -1.509644 0.3678794 0.1353353 [2,] -5.632571 1.4715178 0.4060058 [3,] -4.934938 1.1036383 0.5413411 attr(,"accuracy") [1] 0 > > m3.t <- matrix(c(-1.50964415879218, -5.6325707998812, -4.934938326092, + 0.367879439109187, 1.47151775849686, 1.10363831732856, + 0.135335281175235, 0.406005843524598, 0.541341126763207), + 3,3) > > stopifnot(all.equal(m3.t, m3, check.attributes=FALSE, tolerance = 3e-11), + # ^^^^^ + # 1.2455e-11 for libatlas (above) + all.equal(m3.t, m3T, check.attributes=FALSE, tolerance = 1e-11), + all.equal(m3.t, m3O, check.attributes=FALSE, tolerance = 8e-11),# M1: 1.39e-11 + all.equal(m3.t, m3TO, check.attributes=FALSE, tolerance = 1e-11), + all.equal(m3.t, expm(T3,"R_Eigen"), tolerance = 1e-11), + all.equal(m3.t, expm(T3,"Ward77"), tolerance = 1e-11), + all.equal(m3.t, expm(T3,"R_Ward"), tolerance = 1e-11), + all.equal(m3.t, expm(T3,"R_Pade"), tolerance = 1e-11), + TRUE) > ## -- in this case, a similar level of agreement with Ward (1977). > > ## ---------------------------- > ## Test case 4 from Ward (1977) > ## ---------------------------- > T4 <- + array(c(0, 0, 0, 0, 0, 0, 0, 0, 0, 1e-10, + 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 1, 0), + dim = c(10, 10)) > (m4 <- expm(T4, method="Pade")) [,1] [,2] [,3] [,4] [,5] [1,] 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 [2,] 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 [3,] 2.480159e-15 2.755567e-16 1.000000e+00 1.000000e+00 5.000000e-01 [4,] 1.984127e-14 2.480059e-15 2.755582e-16 1.000000e+00 1.000000e+00 [5,] 1.388889e-13 1.984077e-14 2.480067e-15 2.755581e-16 1.000000e+00 [6,] 8.333333e-13 1.388869e-13 1.984080e-14 2.480067e-15 2.755581e-16 [7,] 4.166667e-12 8.333273e-13 1.388870e-13 1.984080e-14 2.480067e-15 [8,] 1.666667e-11 4.166654e-12 8.333274e-13 1.388869e-13 1.984079e-14 [9,] 5.000000e-11 1.666667e-11 4.166667e-12 8.333333e-13 1.388889e-13 [10,] 1.000000e-10 5.000000e-11 1.666667e-11 4.166667e-12 8.333333e-13 [,6] [,7] [,8] [,9] [,10] [1,] 8.333333e-03 1.388889e-03 1.984127e-04 2.480159e-05 2.755732e-06 [2,] 4.166667e-02 8.333333e-03 1.388889e-03 1.984127e-04 2.480159e-05 [3,] 1.666667e-01 4.166667e-02 8.333333e-03 1.388889e-03 1.984127e-04 [4,] 5.000000e-01 1.666667e-01 4.166667e-02 8.333333e-03 1.388889e-03 [5,] 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 8.333333e-03 [6,] 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 [7,] 2.755583e-16 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 [8,] 2.480064e-15 2.755577e-16 1.000000e+00 1.000000e+00 5.000000e-01 [9,] 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 [10,] 1.388889e-13 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 attr(,"accuracy") [1] 1.440257e-16 > (m4O <- expm(T4, method="PadeO")) [,1] [,2] [,3] [,4] [,5] [1,] 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 [2,] 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 [3,] 2.480159e-15 2.755567e-16 1.000000e+00 1.000000e+00 5.000000e-01 [4,] 1.984127e-14 2.480059e-15 2.755582e-16 1.000000e+00 1.000000e+00 [5,] 1.388889e-13 1.984077e-14 2.480067e-15 2.755581e-16 1.000000e+00 [6,] 8.333333e-13 1.388869e-13 1.984080e-14 2.480067e-15 2.755581e-16 [7,] 4.166667e-12 8.333273e-13 1.388870e-13 1.984080e-14 2.480067e-15 [8,] 1.666667e-11 4.166654e-12 8.333274e-13 1.388869e-13 1.984079e-14 [9,] 5.000000e-11 1.666667e-11 4.166667e-12 8.333333e-13 1.388889e-13 [10,] 1.000000e-10 5.000000e-11 1.666667e-11 4.166667e-12 8.333333e-13 [,6] [,7] [,8] [,9] [,10] [1,] 8.333333e-03 1.388889e-03 1.984127e-04 2.480159e-05 2.755732e-06 [2,] 4.166667e-02 8.333333e-03 1.388889e-03 1.984127e-04 2.480159e-05 [3,] 1.666667e-01 4.166667e-02 8.333333e-03 1.388889e-03 1.984127e-04 [4,] 5.000000e-01 1.666667e-01 4.166667e-02 8.333333e-03 1.388889e-03 [5,] 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 8.333333e-03 [6,] 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 [7,] 2.755583e-16 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 [8,] 2.480064e-15 2.755577e-16 1.000000e+00 1.000000e+00 5.000000e-01 [9,] 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 [10,] 1.388889e-13 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 attr(,"accuracy") [1] 1.440257e-16 > (m4T <- expm(T4,method="Taylor")) [,1] [,2] [,3] [,4] [,5] [1,] 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 [2,] 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 [3,] 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 [4,] 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 [5,] 1.388889e-13 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 [6,] 8.333333e-13 1.388889e-13 1.984127e-14 2.480159e-15 2.755732e-16 [7,] 4.166667e-12 8.333333e-13 1.388889e-13 1.984127e-14 2.480159e-15 [8,] 1.666667e-11 4.166667e-12 8.333333e-13 1.388889e-13 1.984127e-14 [9,] 5.000000e-11 1.666667e-11 4.166667e-12 8.333333e-13 1.388889e-13 [10,] 1.000000e-10 5.000000e-11 1.666667e-11 4.166667e-12 8.333333e-13 [,6] [,7] [,8] [,9] [,10] [1,] 8.333333e-03 1.388889e-03 1.984127e-04 2.480159e-05 2.755732e-06 [2,] 4.166667e-02 8.333333e-03 1.388889e-03 1.984127e-04 2.480159e-05 [3,] 1.666667e-01 4.166667e-02 8.333333e-03 1.388889e-03 1.984127e-04 [4,] 5.000000e-01 1.666667e-01 4.166667e-02 8.333333e-03 1.388889e-03 [5,] 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 8.333333e-03 [6,] 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 [7,] 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 [8,] 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 [9,] 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 [10,] 1.388889e-13 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 attr(,"accuracy") [1] 2.50637e-18 > (m4TO <- expm(T4,method="TaylorO")) [,1] [,2] [,3] [,4] [,5] [1,] 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 [2,] 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 [3,] 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 [4,] 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 [5,] 1.388889e-13 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 [6,] 8.333333e-13 1.388889e-13 1.984127e-14 2.480159e-15 2.755732e-16 [7,] 4.166667e-12 8.333333e-13 1.388889e-13 1.984127e-14 2.480159e-15 [8,] 1.666667e-11 4.166667e-12 8.333333e-13 1.388889e-13 1.984127e-14 [9,] 5.000000e-11 1.666667e-11 4.166667e-12 8.333333e-13 1.388889e-13 [10,] 1.000000e-10 5.000000e-11 1.666667e-11 4.166667e-12 8.333333e-13 [,6] [,7] [,8] [,9] [,10] [1,] 8.333333e-03 1.388889e-03 1.984127e-04 2.480159e-05 2.755732e-06 [2,] 4.166667e-02 8.333333e-03 1.388889e-03 1.984127e-04 2.480159e-05 [3,] 1.666667e-01 4.166667e-02 8.333333e-03 1.388889e-03 1.984127e-04 [4,] 5.000000e-01 1.666667e-01 4.166667e-02 8.333333e-03 1.388889e-03 [5,] 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 8.333333e-03 [6,] 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 4.166667e-02 [7,] 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 1.666667e-01 [8,] 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 5.000000e-01 [9,] 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 1.000000e+00 [10,] 1.388889e-13 1.984127e-14 2.480159e-15 2.755732e-16 1.000000e+00 attr(,"accuracy") [1] 2.50637e-18 > > ## ATLAS on BDR's gannet (Fedora 26; gcc-13; R-devel 2023-10-24) > tol1 <- if(maybeATLAS) 4e-7 else 5e-15 # (m4, m4O) gave "Mean relative difference: 1.242879e-07" > > stopifnot(all.equal(m4 [,10], 1/gamma(10:1), tolerance=1e-14), + all.equal(m4O [,10], 1/gamma(10:1), tolerance=1e-14), + all.equal(m4T [,10], 1/gamma(10:1), tolerance=1e-14), + all.equal(m4TO[,10], 1/gamma(10:1), tolerance=1e-14), + all.equal(m4, m4O, check.attributes=FALSE, tolerance=tol1), + all.equal(m4, m4T, check.attributes=FALSE, tolerance=tol1), + all.equal(m4, m4TO,check.attributes=FALSE, tolerance=tol1), + all.equal(m4, expm(T4,"Ward77"), check.attributes=FALSE, tolerance = 1e-14), + all.equal(m4, expm(T4,"R_Ward"), check.attributes=FALSE, tolerance = 1e-14), + all.equal(m4, expm(T4,"R_Pade"), check.attributes=FALSE, tolerance = 1e-14), + max(abs(m4 - expm(T4,"R_Eigen"))) < 1e-7) > ## here expm(., EV ) is accurate only to 7 d.p., whereas > ## expm(.,Pade) is correct to at least 14 d.p. > > ### Test case with diagonalizable matrix with multiple Eigen values: > A4 <- rbind(c(-1, 3, -1), + c(-3, 5, -1), + c(-3, 3, 1)) > Ea4 <- eigen(A4) > stopifnot(all.equal(Ea4$values, (lam4 <- c(2,2,1)))) > ## However, the eigen values don't show the nice property: > V4 <- Ea4$vectors > crossprod(V4) [,1] [,2] [,3] [1,] 1.0000000 0.7742943 0.9258201 [2,] 0.7742943 1.0000000 0.9408751 [3,] 0.9258201 0.9408751 1.0000000 > ## i.e., they are *not* orthogonal > ## but still diagonalize: > stopifnot(all.equal(A4, + V4 %*% diag(x=lam4) %*% solve(V4))) > ## whereas this diagonalizes *and* looks nice > W4 <- rbind(c(1, 1, -1), + c(1, 1, 0), + c(1, 0, 3)) > (sW4 <- solve(W4)) [,1] [,2] [,3] [1,] 3 -3 1 [2,] -3 4 -1 [3,] -1 1 0 > assert.EQ(diag(x = c(1,2,2)), solve(W4) %*% A4 %*% W4, giveRE=TRUE) > assert.EQ(A4, logm(expm(A4)), tol = 3e-13, giveRE=TRUE) Mean relative difference: 3.848773e-15 > ## seen 5.5e-14 with R's own matprod > > proc.time() user system elapsed 1.48 0.14 1.62