Error Doing Cholesky Decomposition

[Arya Pourzanjani]

I have the following function getting called in a model block and failing at the Cholesky decompose line with the error "cholesky_decompose: Cholesky decomposition of m failed"

 

  matrix chol_solve(matrix K, int M, matrix A, int N) {

    matrix[M,M] L;

    matrix[M,N] K_inv_A;

    

    print("K=", K);

    L = cholesky_decompose(K);

    K_inv_A = mdivide_left_tri_low(K, A);

    

    return K_inv_A;

  }

with the matrix

> K
          [,1]      [,2]     [,3]
[1,]  5.006080  23.39400 0.108915
[2,] 23.394000 146.81100 4.484610
[3,]  0.108915   4.48461 0.423998
> chol(K)
         [,1]      [,2]        [,3]
[1,] 2.237427 10.455760 0.048678681
[2,] 0.000000  6.122752 0.649321995
[3,] 0.000000  0.000000 0.003054872

As illustrated here, the Cholesky decomposition of the matrix in R works just fine. What's the best way to debug this sort of issue?

[Daniel Lee]

Hi,

It's really hard to help without the rest of the program. It's impossible to guess as to how you're using you're function and depending on what you're doing will determine how to help.

Some examples of what you could be doing:

- that function is used to set transformed data

- it's used within transformed parameters

- it's used within a loop

- there's some other call to cholesky_decompose() that's actually causing the problem

If you want help, please provide the whole Stan program (even if the problem is right within that little section). And since you're printing the matrix, can you show what's being printed?

What version of RStan are you using? Have you tried exposing the Stan function and passing the same K into the function from within R?

Daniel

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.
To post to this group, send email to stan-...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Arya Pourzanjani]

Ah, sorry about that. My model file is sort of hard to read which is why I hesitated, but I'll include it now along with my data initial MCMC point (theta_0), you just have to load the Rdata file. I used the following call

fit <- stan(file = 'gplvm_spgp.stan', data = synth_dat, 

            iter = 10, chains = 1, init = theta_0)

I'm using Rstan 2.11.1 and here is what's being printed

"K=[[5.00608,23.394,0.108915],[23.394,146.811,4.48461],[0.108915,4.48461,0.423998]]"  

How do I expose the Stan function and call it from R? That sounds like a useful debugging technique.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.

[Daniel Lee]

Ok... this is a little tricky. I'm running it out of CmdStan to see if I can see what's causing it.

(It is indeed where you pointed first, but it really does help to see the rest of the program.)

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Arya Pourzanjani]

How do you guys (Stan development team) usually debug? Is there an easy way to use GDB? I'd like to learn the control flow of Stan so I can understand where the bottlenecks are in my models, and maybe even contribute. What's the best way to do that?

[Daniel Lee]

I'm usually always using CmdStan. It's such a lightweight wrapper that it's easy to change a bit of the math library and see the changes percolate all the way through without needing to do much to rebuild.

If you're wanting to get that set up, follow these directions: https://github.com/stan-dev/cmdstan/wiki/Getting-Started-with-CmdStan

You can do the same with the math library (that was probably written before we split out the math submodule).

I haven't been successful getting gdb running, but that might be cause I'm on a Mac. I think Ben might have gotten it up and running.

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Arya Pourzanjani]

So do you debug with print statements and just rebuild every time you add in or take away a print statement? Is rebuilding pretty quick? Also is there a document or diagram that says which files get called when you run a command and which files get called from those files or should I just look through what the makefile is doing?

[Daniel Lee]

Yes, exactly. Rebuilding is quick. In this case, I've tracked it down to cholesky_decompose() within the math library, so I'm building a unit test for it. At that point, to build that particular test is 13 seconds. (Not instantaneous, but much faster than building everything.) And there isn't a diagram. The makefile for CmdStan isn't quite as robust as it needs to be, so the dependencies isn't worked all the way through.

There's a bug somewhere, I think. I can see that choleky_decompose() fails from the Stan program, but the math library doesn't. I'm headed to sleep now, but will hunt it down in the morning.

Daniel

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

It varies by member. But all the usual advice about
debugging any computer program holds. I usually debug
everything (including models) by printing. The first
step is always to get a minimal reproducible example that
isolates the bug. We almost always then turn that into
a unit test. The first suspect's always the math library
(stan-dev/math), then if I can verify the problem's not there
(because the unit test passes), I work up to the Stan library
(stan-dev/stan), and then up to CmdStan (stan-dev/cmdstan).

I work in emacs and a unix shell and almost always debug
by printing. Even when I can get something like gdb
working, I've never been good at getting the checkpointing
and stepping working, so I almost always use a divide-and-conquer
approach with print statements and tests within the code that
verify things do what I think they're going to do.

It's a pain that our templated C++ code takes so long to
rebuild, but it's better at the math library level than
trying to compile and run a Stan program. The unit tests
are set up so they can be run individually.

- Bob

[Daniel Lee]

So... the issue is that there's an issue with the Eigen C++ library for that particular K. In C++ (on a Mac, using clang++ 703.0.29), this matrix will get through the cholesky decomposition just fine:

  Eigen::MatrixXd K(3, 3);

   K <<

      5.006077533027154,   23.39404564220951,   0.1089147679987995,

     23.39404564220951,   146.8106822579209,   4.484612487682218,

     0.1089147679987995,  4.484612487682218,   0.4239981299549255;

and this won't:

  Eigen::MatrixXd K(3, 3);

  K <<

    5.0060775330271543,  23.394045642209505, 0.10891476799879951,

    23.394045642209505,  146.81068225792086,  4.4846124876822184,

    0.10891476799879951,  4.4846124876822184, 0.42399812995492547;

I don't know why it fails, but now that I have a reproducible example, I'm going to add an issue to the math library, add a test case, and hopefully we'll have a fix soon at the math library.

For the time being, we should try to find you a workaround. Is there a way to slightly perturb your starting point for now and see if that moves things forward?

Daniel

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

> To post to this group, send email to stan-...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Daniel Lee]

I take it back... R will also say that the particular matrix is not positive definite:

K = matrix(c(5.0060775330271543,  23.394045642209505, 0.10891476799879951,

             23.394045642209505,  146.81068225792086,  4.4846124876822184,

             0.10891476799879951,  4.4846124876822184, 0.42399812995492547),

            3, 3);

chol(K)

output:

Error in chol.default(K) : 

  the leading minor of order 3 is not positive definite

[Rob Trangucci]

Arya, 

The matrix you're printing out is rounded, and not precisely the matrix that's getting decomposed while the model is running. The matrix that is printed above is poorly conditioned (i.e. abs(max(eigen(K)$values))/abs(min(eigen(K)$values)) is quite large, on the order 1e7). It's not a problem with cholesky_decompose, but rather the fact that the matrices you're decomposing are close to singular.

Rob

[Arya Pourzanjani]

Rob, that's an astute point about the rounding. That's why the Cholesky solve works in R, but not Stan. You're absolutely right, the non-rounded version should NOT work because it's a 3x3 matrix with rank only equal to 2. 

The reason the rank is only 2 is because it's a linear kernel. So if X is a data matrix with N rows, and D columns, the linear kernel would be X*X^T, which is NxN and has rank at most D. So if D < N the matrix would not be invertible, as was the case with me because I had D=2, N=3. If on the other hand I had a different kernel such as the RBF then D would be uncountably infinite and hence much larger than N, making my kernel invertible. 

This sort of illustrates a hand-wavy way to see why GP's have identifiability issues, as Michael Betancourt has mentioned on here in the past. The "feature-space" is of much higher dimension than the amount of data you have. 

Thank you Daniel and Bob for the debugging tips! I think using the RBF kernel instead of the linear kernel should fix things. I'm going to try that next.

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/ckDe3MJycT4/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+unsubscribe@googlegroups.com.