matStig.RdStig Mortensen wrote on Oct 22, 2007 to the authors of the Matrix
package with subject “Strange result from expm”.
There, he presented the following \(8 \times 8\) matrix for
which the Matrix expm() gave a “strange” result.
As we later researched, the result indeed was wrong: the correct
entries were wrongly permuted. The reason has been in the underlying
source code in Octave from which it had been ported to Matrix.
data(matStig)data(matStig)
as(matStig, "sparseMatrix") # since that prints more nicely.
#> 8 x 8 sparse Matrix of class "dgCMatrix"
#>
#> [1,] 1.725 -0.765 -15.0000 . . . . .
#> [2,] -0.795 0.765 . . . . . .
#> [3,] . . 0.3949 -0.3949 . . 343.4851 .
#> [4,] . . . 0.1456 . . . .
#> [5,] . . . . -1.725 0.795 . .
#> [6,] . . . . 0.765 -0.765 . .
#> [7,] . . . . 15.000 . -0.3949 .
#> [8,] . . . . . . 0.3949 -0.1456
## For more compact printing:
op <- options(digits = 4)
E1 <- expm(matStig, "Ward77", preconditioning="buggy") # the wrong result
as(E1, "sparseMatrix")
#> 8 x 8 sparse Matrix of class "dgCMatrix"
#>
#> [1,] 1.1567 . . . . . . .
#> [2,] -1.8012 3.124 17.170 . 3.459e+03 458.0991 1399.4256 -3.163
#> [3,] -0.5188 . 1.484 . 1.649e+03 397.4250 352.4825 .
#> [4,] . . . 0.8645 1.534e+00 0.3907 0.3022 .
#> [5,] . . . . 2.590e-01 0.2623 . .
#> [6,] . . . . 2.524e-01 0.5757 . .
#> [7,] . . . . 6.101e+00 2.4747 0.6737 .
#> [8,] 7.1255 -3.044 -51.693 . -1.601e+04 -2577.2615 -5376.6619 6.944
str(E2 <- expm(matStig, "Pade"))# the correct one (has "accuracy" attribute)
#> num [1:8, 1:8] 6.94 -3.16 0 0 0 ...
#> - attr(*, "accuracy")= num 9.88e-09
as(E2, "sparseMatrix")
#> 8 x 8 sparse Matrix of class "dgCMatrix"
#>
#> [1,] 6.944 -3.044 -51.693 7.1255 -1.601e+04 -2577.2615 -5376.6619 .
#> [2,] -3.163 3.124 17.170 -1.8012 3.459e+03 458.0991 1399.4256 .
#> [3,] . . 1.484 -0.5188 1.649e+03 397.4250 352.4825 .
#> [4,] . . . 1.1567 . . . .
#> [5,] . . . . 2.590e-01 0.2623 . .
#> [6,] . . . . 2.524e-01 0.5757 . .
#> [7,] . . . . 6.101e+00 2.4747 0.6737 .
#> [8,] . . . . 1.534e+00 0.3907 0.3022 0.8645
attr(E2,"accuracy") <- NULL # don't want it below
E3 <- expm(matStig, "R_Eigen") # even that is fine here
all.equal(E1,E2) # not at all equal (rel.difference >~= 1.)
#> [1] "Mean relative difference: 1.512"
stopifnot(all.equal(E3,E2)) # ==
##________ The "proof" that "Ward77" is wrong _________
M <- matStig
Et1 <- expm(t(M), "Ward77", precond= "buggy")
Et2 <- expm(t(M), "Pade"); attr(Et2,"accuracy") <- NULL
all.equal(Et1, t(E1)) # completely different (rel.diff ~ 1.7 (platform dep.))
#> [1] "Mean relative difference: 1.786"
stopifnot(all.equal(Et2, t(E2))) # the same (up to tolerance)
options(op)