RcppCore / RcppArmadillo

Rcpp integration for the Armadillo templated linear algebra library
192 stars 56 forks source link

warning: chol(): given matrix is not symmetric #344

Closed UseRSteph closed 3 years ago

UseRSteph commented 3 years ago

Hello @conradsnicta @eddelbuettel I recently double checked some code I wrote some time ago to apply the bvarsv package to estimate Bayesian time varying parameter vector autoregressions, which calls on the RcppArmadillo package. With one pair of series, the code previously had some convergence issues, such that only about half of the specifications based on the assumed parameterizations reported results; I was probably using RcppArmadillo 9.2 and R 3.5. But when I now run the code using the parameters that previously worked with RcppArmadillo 10.6 and R 4.1, or RcppArmadillo 9.2 and R 3.5, I get "warning: chol(): given matrix is not symmetric". I think it may be because of the RcppArmadillo updates, although bvarsv also calls on the coda, R2jags, R2OpenBUGS, Rcpp and vars packages. It will still work if I increase the training sample size for the estimator; in the code "tau=42" produces the error message, while "tau=48" still works. But I wanted to know if there's something I can do to get it to work for "tau=42". Many thanks. GitHub Data.csv GitHub Example.txt

eddelbuettel commented 3 years ago

Thanks for reporting it. You can actually edit code quite nicely here (see the help on github editing, eg here which makes your code example be

library(coda)
library(bvarsv)
library(R2jags)
library(R2OpenBUGS)
library(Rcpp)
library(RcppArmadillo)
library(vars)

data <- read.table("H:/GitHub Data.csv",header=T,sep=",")   # non-portable path, use just "github_data.csv" 

y1 <- ts(data$y1, deltat=1/12)
y2 <- ts(data$y2, deltat=1/12)

datamat <- ts(matrix(cbind(y2,y1),102,2,byrow=F), deltat=1/12)
colnames(datamat) <- c("y2","y1")

set.seed(1900)
btvpsvar <- bvar.sv.tvp(datamat,p=6,tau=42,k_Q=0.01, k_S=0.01, k_W=0.001)

and that is just about the text book example of what one should not do. I really cherish minimally complete verifiable examples so could your possibly create something that does not depend on five other packages we have nothing to do with?

UseRSteph commented 3 years ago

Thank you. I believe your correction is the minimally complete verifiable example that generates the error. Maybe it's not RcppArmadillo.

eddelbuettel commented 3 years ago

That has happened to many of us before: In the process of whittling an issue down to the (initially suspected) core, we find that something along the way may not have been as invariant, or innocent, or ... as we thought. So I will close this; should you find something reproducible please feel free to reopen.

UseRSteph commented 3 years ago

@conradsnicta many thanks for following up and apologies. I discovered that the repeated warnings were not actually errors and that the code for that pair of series actually still runs. The bvarsv package by default does 50000 MCMC draws, and reports progress after every 10000. I had not scrolled down far enough to see that the updates after every 10000 were still appearing; so you'd see 9999 consecutive warnings, and then at each multiple of 10000 draw you get an update. [Also, to run the bvarsv example above, only Rcpp and RcppArmadillo are dependencies; the other packages were necessary for diagnostic purposes.]

eddelbuettel commented 3 years ago

Monte Carlo simulations are a great way to find boundary conditions. I used a lot of try/catch back when I needed it :)

UseRSteph commented 3 years ago

Thank you @conradsnicta. I'll look into this. I checked for that problem with the raw series, but I think the warning message concerns Monte Carlo draws. Seems to me that the warnings, which I didn't get when I first tried the bvarsv package, actually help with the model selection, as I found that a larger training sample size helps eliminate the warnings, or possibly using fewer lags (p < 6) or different priors (k_Q, k_S and k_W in the code).