Log-determinant in RTMB?

[phillipbvetter]

Hi,

Thanks for the fantastic work with TMB / RTMB.

I imagine the current answer is no, but is it possible to calculate a log-determinant using RTMB? I'm currently using atomic::logdet in my TMB C++ function, for evaluating a multivariate normal. I do not use dmvnorm because I like to have the ability to apply a regulariser (e.g. pseudo huber) to the (covariance normalised) sum of squared residuals.

These two approaches will fail.

fn = function(p){
    nll = determinant(p$A,logarithm=TRUE)
    nll = log(prod(eigen(p$A,only.values=TRUE)))
  return(nll)
}
nll = RTMB::MakeADFun(fn, parameters=list(A=diag(2)))

[kaskr]

Yes, log determinant is not yet available in RTMB. I've added a new feature to implement such missing atomics from R. It's in the master branch under ?ADjoint with logdet as an example.

[phillipbvetter]

That is awesome. Thanks a lot!

Update: For anyone else interested in this, the logdet-code copied from the new feature branch mentioned above reads

logdet <- RTMB::ADjoint(
   function(x) {
       dim(x) <- rep(sqrt(length(x)), 2)
       determinant(x, log=TRUE)$modulus
   },
   function(x, y, dy) {
       dim(x) <- rep(sqrt(length(x)), 2)
       t(RTMB::solve(x)) * dy
   },
   name = "logdet")

so that can be run after grabbing the package from source via

remotes::install_github("https://github.com/kaskr/RTMB", subdir="RTMB")

and log-determinants can then be calculated in the likelihood function, and parsed to MakeADFun without any errors

neg.log.like = function(p){
    nll = logdet(p$A)
  return(nll)
}
nll = RTMB::MakeADFun(neg.log.like, parameters=list(A=diag(2)))

[phillipbvetter]

Just a heads up - I may have discovered an issue. Trying to calculate the fixed effects hessian causes my session to crash i.e.

neg.log.like = function(p){
    nll = logdet(p$A)
  return(nll)
}
nll = RTMB::MakeADFun(neg.log.like, parameters=list(A=diag(2)))
nll$he(nll$par)

This is my session info

R version 4.3.2 (2023-10-31)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.5

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Copenhagen
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] gtable_0.3.4                 TMB_1.9.11                   xfun_0.41                    ggplot2_3.5.0               
 [5] sdeTMB_0.5                   remotes_2.4.2.1              htmlwidgets_1.6.4            devtools_2.4.5              
 [9] lattice_0.21-9               vctrs_0.6.5                  tools_4.3.2                  generics_0.1.3              
[13] tibble_3.2.1                 fansi_1.0.6                  pkgconfig_2.0.3              Matrix_1.6-4                
[17] data.table_1.14.10           checkmate_2.3.1              webshot_0.5.5                lifecycle_1.0.4             
[21] compiler_4.3.2               metR_0.15.0                  stringr_1.5.1                munsell_0.5.0               
[25] codetools_0.2-19             httpuv_1.6.13                usethis_2.2.2                htmltools_0.5.7             
[29] yaml_2.3.8                   pracma_2.4.4                 urlchecker_1.0.1             later_1.3.2                 
[33] pillar_1.9.0                 RTMB_1.4                     MASS_7.3-60                  ellipsis_0.3.2              
[37] cachem_1.0.8                 sessioninfo_1.2.2            mime_0.12                    tidyselect_1.2.1            
[41] rvest_1.0.3                  digest_0.6.35                stringi_1.8.3                purrr_1.0.2                 
[45] dplyr_1.1.4                  reshape2_1.4.4               latex2exp_0.9.6              fastmap_1.1.1               
[49] grid_4.3.2                   colorspace_2.1-0             cli_3.6.2                    PredSecondaryClarifier_1.1.1
[53] magrittr_2.0.3               patchwork_1.2.0              pkgbuild_1.4.4               utf8_1.2.4                  
[57] clipr_0.8.0                  promises_1.2.1               scales_1.3.0                 backports_1.4.1             
[61] rmarkdown_2.25               httr_1.4.7                   kableExtra_1.3.4             memoise_2.0.1               
[65] shiny_1.8.0                  evaluate_0.23                knitr_1.45                   miniUI_0.1.1.1              
[69] viridisLite_0.4.2            profvis_0.3.8                rlang_1.1.3                  Rcpp_1.0.12                 
[73] isoband_0.2.7                xtable_1.8-4                 glue_1.7.0                   xml2_1.3.6                  
[77] pkgload_1.3.4                svglite_2.1.3                rstudioapi_0.15.0            R6_2.5.1                    
[81] plyr_1.8.9                   systemfonts_1.0.5            fs_1.6.3         

[kaskr]

@phillipbvetter Could it be that you 'forgot' to run library(RTMB)? That would be needed for the correct solve function method to be selected... If I do not load RTMB I can reproduce your crash and get:

Process R aborted (core dumped) at Mon Apr 22 14:43:21 2024
Error in solve.default(x) : 
  Lapack routine zgesv: system is exactly singular: U[2,2] = 0
terminate called after throwing an instance of 'Rcpp::LongjumpException'

Edit: I can alternatively avoid the crash by using RTMB::solve rather than solve.

[phillipbvetter]

Ah yeah that was it! I did not think about that at all...Thanks again, @kaskr.

[kaskr]

The crash was definitely a bug in RTMB though - should now be fixed