Righto - but do they need different priors? (which is what LKJCholskyCov controls)
For example this snippet would give you a 2D covariance matrix of 2 features for each of n observations, but they have a shared prior
sd_dist = pm.InverseGamma.dist(alpha=101., beta=100.)
chol, corr_, stds_ = pm.LKJCholeskyCov('lkjcc', n=2, eta=2.,
sd_dist=sd_dist, compute_corr=True)
mvn_dist = pm.MvNormal.dist(mu=tt.zeros(2), chol=chol, shape=(n_obs, 2))
mvn_like = pm.Potential('mvn_like', (mvn_dist.logp(y)))