# Does pymc3 have a cholesky decomposition function?

I have a matrix representing the distance between random normal variables. I’d like to estimate an additional variable representing the rate of decline in co-variance between these variables based on their distance. The centered implementation looks like this:

``````distance_matrix = np.array([[0,1,2],
[1,0,1],
[2,1,0]])

# centered implementation
rho_sq = pm.HalfNormal('rho_sq', sd=1) # rate of decline by distance
cov = pm.math.exp(-rho_sq * distance_matrix**2)
pm.MvNormal('correlated_vars', mu=0, cov=cov, shape=len(distance_matrix))
``````

The non centered implementation requires a cholesky decomposition. How can I do that in pymc3?

``````distance_matrix = np.array([[0,1,2],
[1,0,1],
[2,1,0]])

# non centered implementation
rho_sq = pm.HalfNormal('rho_sq', sd=1) # rate of decline by distance
cov = pm.math.exp(-rho_sq * distance_matrix**2)

# this does not work
chol = np.linalg.cholesky(cov)

uncorrelated_vars = pm.Normal('uncorrelated_vars', mu=0, sd=1, shape=len(distance_matrix))
correlated_vars = pm.Deterministic('correlated_vars', tt.dot(chol, uncorrelated_vars.T))``````

Hi there I think this notebook answers your question - it’s the adaptation to PyMC3 of McElreath’s Rethinking chapter 13.
Up until `Code 13.28`, it deals with categorical varying effects, and then it uses a Gaussian process to include the distances between islands as a varying effect

Thanks for your response. The example in that book uses the centered implementation where you pass a covariance matrix to a multi variate normal distribution. I’d like to figure out a non centered or alternative implementation that samples more quickly.

`theano.tensor.slinalg.cholesky` should do what you want.
See http://deeplearning.net/software/theano/library/tensor/slinalg.html

1 Like

Not entirely sure, but I think the cell titled "A reparameterization of m13.7 " in the NB I linked is using a non-centered parametrization:

``````gmu = pm.Normal('gmu', 0., 1., shape=Nsociety)
g = pm.Deterministic('g', pm.math.dot(tt.slinalg.cholesky(Kij), gmu))
``````