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. > library(expm) Loading required package: Matrix Attaching package: 'expm' The following object is masked from 'package:Matrix': expm > > (sI <- sessionInfo()) R Under development (unstable) (2024-08-17 r87027 ucrt) Platform: x86_64-w64-mingw32/x64 Running under: Windows Server 2022 x64 (build 20348) Matrix products: default locale: [1] LC_COLLATE=C LC_CTYPE=German_Germany.utf8 [3] LC_MONETARY=C LC_NUMERIC=C [5] LC_TIME=C time zone: Europe/Berlin tzcode source: internal attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] expm_1.0-0 Matrix_1.7-0 loaded via a namespace (and not attached): [1] compiler_4.5.0 grid_4.5.0 lattice_0.22-6 > packageDescription("Matrix") Package: Matrix Version: 1.7-0 VersionNote: do also bump src/version.h, inst/include/Matrix/version.h Date: 2024-03-16 Priority: recommended Title: Sparse and Dense Matrix Classes and Methods Description: A rich hierarchy of sparse and dense matrix classes, including general, symmetric, triangular, and diagonal matrices with numeric, logical, or pattern entries. Efficient methods for operating on such matrices, often wrapping the 'BLAS', 'LAPACK', and 'SuiteSparse' libraries. License: GPL (>= 2) | file LICENCE URL: https://Matrix.R-forge.R-project.org BugReports: https://R-forge.R-project.org/tracker/?atid=294&group_id=61 Contact: Matrix-authors@R-project.org Authors@R: c(person("Douglas", "Bates", role = "aut", comment = c(ORCID = "0000-0001-8316-9503")), person("Martin", "Maechler", role = c("aut", "cre"), email = "mmaechler+Matrix@gmail.com", comment = c(ORCID = "0000-0002-8685-9910")), person("Mikael", "Jagan", role = "aut", comment = c(ORCID = "0000-0002-3542-2938")), person("Timothy A.", "Davis", role = "ctb", comment = c(ORCID = "0000-0001-7614-6899", "SuiteSparse libraries", "collaborators listed in dir(system.file(\"doc\", \"SuiteSparse\", package=\"Matrix\"), pattern=\"License\", full.names=TRUE, recursive=TRUE)")), person("George", "Karypis", role = "ctb", comment = c(ORCID = "0000-0003-2753-1437", "METIS library", "Copyright: Regents of the University of Minnesota")), person("Jason", "Riedy", role = "ctb", comment = c(ORCID = "0000-0002-4345-4200", "GNU Octave's condest() and onenormest()", "Copyright: Regents of the University of California")), person("Jens", "Oehlschlägel", role = "ctb", comment = "initial nearPD()"), person("R Core Team", role = "ctb", comment = "base R's matrix implementation")) Depends: R (>= 4.4.0), methods Imports: grDevices, graphics, grid, lattice, stats, utils Suggests: MASS, datasets, sfsmisc, tools Enhances: SparseM, graph LazyData: no LazyDataNote: not possible, since we use data/*.R and our S4 classes BuildResaveData: no Encoding: UTF-8 NeedsCompilation: yes Packaged: 2024-03-19 17:15:14 UTC; maechler Author: Douglas Bates [aut] (), Martin Maechler [aut, cre] (), Mikael Jagan [aut] (), Timothy A. Davis [ctb] (, SuiteSparse libraries, collaborators listed in dir(system.file("doc", "SuiteSparse", package="Matrix"), pattern="License", full.names=TRUE, recursive=TRUE)), George Karypis [ctb] (, METIS library, Copyright: Regents of the University of Minnesota), Jason Riedy [ctb] (, GNU Octave's condest() and onenormest(), Copyright: Regents of the University of California), Jens Oehlschlägel [ctb] (initial nearPD()), R Core Team [ctb] (base R's matrix implementation) Maintainer: Martin Maechler Repository: CRAN Date/Publication: 2024-04-26 12:03:02 UTC Built: R 4.5.0; x86_64-w64-mingw32; 2024-08-18 22:58:10 UTC; windows Archs: x64 -- File: D:/RCompile/CRANpkg/lib/4.5/Matrix/Meta/package.rds > packageDescription("expm") Package: expm Type: Package Title: Matrix Exponential, Log, 'etc' Version: 1.0-0 Date: 2024-08-19 Authors@R: c(person("Martin", "Maechler", role=c("aut","cre"), email="maechler@stat.math.ethz.ch", comment = c(ORCID = "0000-0002-8685-9910")) , person("Christophe","Dutang", role = "aut", comment = c(ORCID = "0000-0001-6732-1501")) , person("Vincent", "Goulet", role = "aut", comment = c(ORCID = "0000-0002-9315-5719")) , person("Douglas", "Bates", role = "ctb", comment = "cosmetic clean up, in svn r42") , person("David", "Firth", role = "ctb", comment = "expm(method= \"PadeO\" and \"TaylorO\")") , person("Marina", "Shapira", role = "ctb", comment = "expm(method= \"PadeO\" and \"TaylorO\")") , person("Michael", "Stadelmann", role = "ctb", comment = "\"Higham08*\" methods, see ?expm.Higham08...") ) Contact: expm-developers@lists.R-forge.R-project.org Description: Computation of the matrix exponential, logarithm, sqrt, and related quantities, using traditional and modern methods. Depends: Matrix Imports: methods Suggests: RColorBrewer, sfsmisc, Rmpfr BuildResaveData: no License: GPL (>= 2) URL: https://R-Forge.R-project.org/projects/expm/ BugReports: https://R-forge.R-project.org/tracker/?atid=472&group_id=107 Encoding: UTF-8 NeedsCompilation: yes Packaged: 2024-08-19 09:24:41 UTC; maechler Author: Martin Maechler [aut, cre] (), Christophe Dutang [aut] (), Vincent Goulet [aut] (), Douglas Bates [ctb] (cosmetic clean up, in svn r42), David Firth [ctb] (expm(method= "PadeO" and "TaylorO")), Marina Shapira [ctb] (expm(method= "PadeO" and "TaylorO")), Michael Stadelmann [ctb] ("Higham08*" methods, see ?expm.Higham08...) Maintainer: Martin Maechler Built: R 4.5.0; x86_64-w64-mingw32; 2024-08-19 13:34:16 UTC; windows Archs: x64 -- File: D:/RCompile/CRANincoming/R-devel/lib/expm/Meta/package.rds > > 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). > > ##-------- expm() ------------------- > > z3 <- T3 * (1 + 1i) > Lz3 <- expmAll(z3) > str(Lz3) List of 13 $ Higham08.b : cplx [1:3, 1:3] -0.427+1.3i -1.878+4.82i -1.113+4.26i ... $ Higham08 : cplx [1:3, 1:3] -0.427+1.3i -1.878+4.82i -1.113+4.26i ... $ AlMohy-Hi09 : chr "matexp_MH09(.) is _not yet_ implemented for complex matrices" $ Ward77 : chr "invalid argument: not a numeric matrix" $ PadeRBS : chr "expm(, method='PadeRBS') is not (yet) implemented" $ Pade : chr "expm(, method='Pade') is not (yet) implemented" $ Taylor : chr "expm(, method='Taylor') is not (yet) implemented" $ PadeO : chr "expm(, method='PadeO') is not (yet) implemented" $ TaylorO : chr "expm(, method='TaylorO') is not (yet) implemented" $ R_Eigen : cplx [1:3, 1:3] -0.427+1.3i -1.878+4.82i -1.113+4.26i ... $ R_Pade : cplx [1:3, 1:3] -0.427+1.3i -1.878+4.82i -1.113+4.26i ... $ R_Ward77 : cplx [1:3, 1:3] -0.427+1.3i -1.878+4.82i -1.113+4.26i ... $ hybrid_Eigen_Ward: chr "invalid argument: not a numeric matrix" > Lz3. <- Lz3[.methComplex] > str(allEq(Lz3., tol=0)) # -> max seen (Lnx 64): 1.3376e-12 List of 4 $ Higham08: chr "Mean relative Mod difference: 7.236117e-13" $ R_Eigen : chr "Mean relative Mod difference: 5.542393e-13" $ R_Pade : chr "Mean relative Mod difference: 1.337624e-12" $ R_Ward77: chr "Mean relative Mod difference: 2.757525e-13" > stopifnot(unlist(allEq(Lz3.))) > > > ## ---------------------------- > ## 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.39 0.25 1.57