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 :slight_smile: 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))

which corroborates Adrian’s answer

Thank you both! That is what I was looking for.

1 Like