Skip to contents

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.