MatrixExponential
MatrixExponential.Rmd
Matrix exponential in c++tsem
c++tsem relies heavily on the RcppArmadillo library. This can sometimes lead to results which deviate from ctsem. Especially the matrix exponential can be prone to numerical precision problems.
Example:
library(regCtsem)
#> Loading required package: ctsemOMX
#> Loading required package: ctsem
#> Loading required package: Rcpp
#> ctsem also changes in time, for manual run ctDocs(), for blog see https://cdriver.netlify.app/, for citation info run citation('ctsem'), for original OpenMx functionality install.packages('ctsemOMX'), and for discussion https://github.com/cdriveraus/ctsem/discussions
#> Warning in doTryCatch(return(expr), name, parentenv, handler): The following
#> important packages for ctsem are out of date: rstan, ctsem
#> Loading required package: OpenMx
#>
#> Attaching package: 'ctsemOMX'
#> The following objects are masked from 'package:ctsem':
#>
#> ctFit, ctIndplot
set.seed(53455)
randMat <- matrix(stats::rnorm(n = 64,
mean = 0,
sd = 10),
nrow = 8, ncol = 8)
regCtsem:::armaExpmat(randMat) - OpenMx::expm(randMat)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] -768.5128 232.49915 -961.5552 -342.9163 125.22721 1258.8476
#> [2,] -1209.0286 365.76894 -1512.7241 -539.4778 197.00813 1980.4258
#> [3,] 408.5591 -123.60188 511.1849 182.3022 -66.57366 -669.2322
#> [4,] 378.3700 -114.46876 473.4127 168.8316 -61.65444 -619.7817
#> [5,] 1343.3712 -406.41174 1680.8121 599.4224 -218.89891 -2200.4831
#> [6,] 2867.5182 -867.51385 3587.8090 1279.5085 -467.25479 -4697.0826
#> [7,] 287.3833 -86.94242 359.5710 128.2326 -46.82837 -470.7426
#> [8,] 1806.3732 -546.48431 2260.1154 806.0175 -294.34392 -2958.8947
#> [,7] [,8]
#> [1,] -252.35180 1060.8234
#> [2,] -397.00122 1668.8932
#> [3,] 134.15601 -563.9581
#> [4,] 124.24302 -522.2864
#> [5,] 441.11446 -1854.3340
#> [6,] 941.58920 -3958.2036
#> [7,] 94.36626 -396.6920
#> [8,] 593.14757 -2493.4430
Note that there are huge differences in the matrix exponential using RcppArmadillo and OpenMx.