segfault when some columns of Y are all 0

[davharris]

code to reproduce on my mac, with R 3.2.0:

set.seed(1)
Y = matrix(rbinom(150, size = 1, prob = .1), ncol = 10)
colSums(Y)
library(BayesComm)
bc = BC(Y = Y, model = "community")

The easiest fix is probably to throw an error at the start of BC or BC.fit. What do you think?

[davharris]

Actually, I'm getting a different error now for some reason, not a segfault:

error: chol(): failed to converge

libc++abi.dylib: terminating with uncaught exception of type std::runtime_error: chol(): failed to converge
Abort trap: 6

[stevencarlislewalker]

Same ... on the off chance it helps. This could still ultimately be a memory bug.

> set.seed(1)
> Y = matrix(rbinom(150, size = 1, prob = .1), ncol = 10)
> colSums(Y)
 [1] 2 2 0 0 1 1 1 2 2 1
> library(BayesComm)
Loading required package: Rcpp
Loading required package: RcppArmadillo
Loading required package: coda
Loading required package: lattice
> sessionInfo()
R version 3.2.0 (2015-04-16)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.9.5 (Mavericks)

locale:
[1] C

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

other attached packages:
[1] BayesComm_0.1-0         coda_0.16-1             lattice_0.20-31        
[4] RcppArmadillo_0.4.400.0 Rcpp_0.11.6             setup_0.0-1            

loaded via a namespace (and not attached):
[1] grid_3.2.0
> bc = BC(Y = Y, model = "community")

error: chol(): failed to converge

libc++abi.dylib: terminating with uncaught exception of type std::runtime_error: chol(): failed to converge

Process R abort trap: 6 at Wed Jul 15 17:41:28 2015

[goldingn]

I suspect this is just a non-positive definite matrix trying to be inverted - there should be a work around (e.g. by adding a small amount of jitter to the diagonal).

I'll be honest, this package isn't a major priority for me any more. Partly as I'm planning to supercede it with a more general Gaussian process package which I'll be focussing on from next year. Fixing this, adhering to the new Rcpp interface and making it all pass CRAN checks will take a bit of time.

I'd be happy to do it if it would be useful - @stevencarlislewalker do you have something you need the package for?

[stevencarlislewalker]

No I didn't have anything I need BayesComm for. I just like keeping up with what you guys are doing, so when github told me about the issue I thought I'd reproduce on my system in case it would be useful. I'm looking forward to seeing the GP stuff though!.

About jittering the diagonal ... under most circumstances and in an ideal world I would rather either (1) test for PSDness and throw an error at an appropriate place or (2) build in explicit regularization/penalization.

[goldingn]

Cool, thanks for checking it out.

Re. failed Choleskies - that's helpful to know.

If I were to go down the jittering approach, I'd issue a warning notifying the user how much jitter was added, and then throw an error if it got too* much. I'd personally rather the thing ran than throw a PSD error to users who might give up when there are solutions to their problems...

Re. option 2 - any references for those approaches?

*to be determined

[goldingn]

Also, I thought for some reason that this was an old thread - I only just realised it was newly opened by @davharris. Dave, presumably you're applying the package to something then?

(sorry, jetlagged brain)

[davharris]

Yeah, but I can deal with malformed inputs manually. No need to update the package on my behalf.

Thanks for double-checking.

On Jul 16, 2015, at 3:41 PM, Nick Golding notifications@github.com wrote:

Also, I thought for some reason that this was an old thread - I only just realised it was newly opened by @davharris https://github.com/davharris. Dave, presumably you're applying the package to something then?

— Reply to this email directly or view it on GitHub https://github.com/goldingn/BayesComm/issues/13#issuecomment-122123613.