helper class for Gaussians

[karlnapf]

Hi there, many samplers are based on Gaussian distributions, that need to be sampled from and whose log-pdf needs to be evaluated. If the Gaussian's covariance is not diagonal, this requires some computation. I think in the current pymc3 branch, things are sometimes re-computed (see #592), so what about adding a helper class for Gaussians that just decompose the covariance once, and then cheap to sample/log-pdf from?

I use:

class Gaussian():
    """
    Helper class to sample from and evaluate log-pdf of a multivariate Gaussian,
    using efficient Cholesky based representation (Cholesky only computed once)
    """
    def __init__(self, mu, Sigma, is_cholesky=False):
        self.mu = mu
        self.is_cholesky = is_cholesky

        if self.is_cholesky:
            self.L = Sigma
        else:
            self.L = np.linalg.cholesky(Sigma)

        self.dimension = len(mu)

    def log_pdf(self, X):
        # duck typing for shape
        if len(np.shape(X)) == 1:
            X = X.reshape(1, len(X))

        log_determinant_part = -sum(np.log(np.diag(self.L)))

        quadratic_parts = np.zeros(len(X))
        for i in range(len(X)):
            x = X[i] - self.mu

            # solve y=K^(-1)x = L^(-T)L^(-1)x
            y = sp.linalg.solve_triangular(self.L, x.T, lower=True)
            y = sp.linalg.solve_triangular(self.L.T, y, lower=False)
            quadratic_parts[i] = -0.5 * x.dot(y)

        const_part = -0.5 * len(self.L) * np.log(2 * np.pi)

        return const_part + log_determinant_part + quadratic_parts

    def sample(self, n=1):
        V = np.random.randn(self.dimension, n)

        # map to our desired Gaussian and transpose to have row-wise vectors
        return self.L.dot(V).T + self.mu

[karlnapf]

BTW theano does not support pre-computing linear solves, see https://github.com/Theano/Theano/issues/2032 Thats why I use scipy for triangular solves based on Cholesky factors

[twiecki]

@karlnapf Just getting back to this. The new Wishart can just provide the Cholesky decomposed cov matrix. Using that together with this transform Gaussian logp should speed things up if understand correctly (so there's a benefit besides the faster sampling). Correct?

[karlnapf]

Yes, great. The Kameleon MCMC produces a proposal that is a Gaussian which is given by a DxN matrix L, such that the covariance is C = L.dot(L.T) For the adaptive Metropolis Hastings, we need to allow for low-rank updates of the Cholesky factor. I am currently busy (I know I always say this), but really hope to get back on this soon as I am working on another kernel based sampler which needs to be compared against Kameleon.

Gaussians with precomputed Cholesky factors is extremely useful for any ML code

[twiecki]

Ah, I see. If I understand you correctly, your logpdf can't be used with pymc3 yet because theano does not support triangular solve yet.

[karlnapf]

yep, exactly

[jsalvatier]

@karlnapf I think you could implement a triangular solve pretty easily in theano by wrapping scipy.linalg.triangular_solve with as_op while simultaneously pushing the theano devs to add a C version. Now that we're getting automatic transformations it would be great to have a transform based on a cholesky factorization.

[twiecki]

Good point, there are lots of examples of doing that here: https://github.com/Theano/Theano/blob/master/theano/tensor/slinalg.py

[karlnapf]

Ok I will have a look

[springcoil]

How does this look like based on new improvements in Theano?

[springcoil]

Is this resolved now by some of the other changes @twiecki ?

[springcoil]

Since there's been no commentary on this I vote we close it.