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] (<https://orcid.org/0000-0001-8316-9503>),
Martin Maechler [aut, cre]
(<https://orcid.org/0000-0002-8685-9910>), Mikael Jagan [aut]
(<https://orcid.org/0000-0002-3542-2938>), Timothy A. Davis
[ctb] (<https://orcid.org/0000-0001-7614-6899>, SuiteSparse
libraries, collaborators listed in dir(system.file("doc",
"SuiteSparse", package="Matrix"), pattern="License",
full.names=TRUE, recursive=TRUE)), George Karypis [ctb]
(<https://orcid.org/0000-0003-2753-1437>, METIS library,
Copyright: Regents of the University of Minnesota), Jason Riedy
[ctb] (<https://orcid.org/0000-0002-4345-4200>, 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 <mmaechler+Matrix@gmail.com>
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-09
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-09 16:43:33 UTC; maechler
Author: Martin Maechler [aut, cre]
(<https://orcid.org/0000-0002-8685-9910>), Christophe Dutang
[aut] (<https://orcid.org/0000-0001-6732-1501>), Vincent Goulet
[aut] (<https://orcid.org/0000-0002-9315-5719>), 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 <maechler@stat.math.ethz.ch>
Built: R 4.5.0; x86_64-w64-mingw32; 2024-08-19 08:55:09 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(<complex>) -------------------
>
> 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(<complex>, method='PadeRBS') is not (yet) implemented"
$ Pade : chr "expm(<complex>, method='Pade') is not (yet) implemented"
$ Taylor : chr "expm(<complex>, method='Taylor') is not (yet) implemented"
$ PadeO : chr "expm(<complex>, method='PadeO') is not (yet) implemented"
$ TaylorO : chr "expm(<complex>, 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.32 0.25 1.56