https://github.com/GPflow/GPflow
Tip revision: 5ab5aaee8696bcb7bcafe6a8cb9334e5ff9186ad authored by mqbssaby on 11 January 2017, 14:12:37 UTC
contrast full and diag model bounds
contrast full and diag model bounds
Tip revision: 5ab5aae
gplvm.py
import tensorflow as tf
import numpy as np
from .model import GPModel
from .gpr import GPR
from . import kullback_leiblers
from .param import Param
from .mean_functions import Zero
from . import likelihoods
from .tf_wraps import eye
from . import transforms
from . import kernels
from ._settings import settings
float_type = settings.dtypes.float_type
def PCA_reduce(X, Q):
"""
A helpful function for linearly reducing the dimensionality of the data X
to Q.
:param X: data array of size N (number of points) x D (dimensions)
:param Q: Number of latent dimensions, Q < D
:return: PCA projection array of size N x Q.
"""
assert Q <= X.shape[1], 'Cannot have more latent dimensions than observed'
evecs, evals = np.linalg.eigh(np.cov(X.T))
i = np.argsort(evecs)[::-1]
W = evals[:, i]
W = W[:, :Q]
return (X - X.mean(0)).dot(W)
class GPLVM(GPR):
"""
Standard GPLVM where the likelihood can be optimised with respect to the latent X.
"""
def __init__(self, Y, latent_dim, X_mean=None, kern=None, mean_function=Zero()):
"""
Initialise GPLVM object. This method only works with a Gaussian likelihood.
:param Y: data matrix, size N (number of points) x D (dimensions)
:param X_mean: latent positions (N x Q), for the initialisation of the latent space.
:param kern: kernel specification, by default RBF
:param mean_function: mean function, by default None.
"""
if kern is None:
kern = kernels.RBF(latent_dim, ARD=True)
if X_mean is None:
X_mean = PCA_reduce(Y, latent_dim)
assert X_mean.shape[1] == latent_dim, \
'Passed in number of latent ' + str(latent_dim) + ' does not match initial X ' + str(X_mean.shape[1])
self.num_latent = X_mean.shape[1]
assert Y.shape[1] >= self.num_latent, 'More latent dimensions than observed.'
GPR.__init__(self, X_mean, Y, kern, mean_function=mean_function)
del self.X # in GPLVM this is a Param
self.X = Param(X_mean)
class BayesianGPLVM(GPModel):
def __init__(self, X_mean, X_std, Y, kern, M, Z=None, X_prior_mean=None, X_prior_var=None):
"""
Initialise Bayesian GPLVM object. This method only works with a Gaussian likelihood.
:param X_mean: initial latent positions, size N (number of points) x Q (latent dimensions).
:param X_std: Chol of latent variances (N x Q or N x Q x Q), for the initialisation of the latent space.
:param Y: data matrix, size N (number of points) x D (dimensions)
:param kern: kernel specification, by default RBF
:param M: number of inducing points
:param Z: matrix of inducing points, size M (inducing points) x Q (latent dimensions). By default
random permutation of X_mean.
:param X_prior_mean: prior mean used in KL term of bound. By default 0. Same size as X_mean.
:param X_prior_var: pripor variance used in KL term of bound. By default 1.
"""
GPModel.__init__(self, X_mean, Y, kern, likelihood=likelihoods.Gaussian(), mean_function=Zero())
del self.X # in GPLVM this is a Param
self.X_mean = Param(X_mean)
# diag_transform = transforms.DiagMatrix(X_var.shape[1])
# self.X_var = Param(diag_transform.forward(transforms.positive.backward(X_var)) if X_var.ndim == 2 else X_var,
# diag_transform)
if X_std.ndim == 2:
self.X_std = Param(X_std, transforms.positive)
elif X_std.ndim == 3:
self.X_std = Param(X_std)
else:
raise AssertionError("Incorrect number of dimensions for X_std.")
self.num_data = X_mean.shape[0]
self.output_dim = Y.shape[1]
if X_std.ndim == 2:
assert np.all(X_mean.shape == X_std.shape)
elif X_std.ndim == 3:
assert X_mean.shape[1] == X_std.shape[1] == X_std.shape[2]
assert X_mean.shape[0] == Y.shape[0], 'X mean and Y must be same size.'
assert X_std.shape[0] == Y.shape[0], 'X var and Y must be same size.'
# inducing points
if Z is None:
# By default we initialize by subset of initial latent points
Z = np.random.permutation(X_mean.copy())[:M]
else:
assert Z.shape[0] == M
self.Z = Param(Z)
self.num_latent = Z.shape[1] # TODO: Better name for num_latent.
assert X_mean.shape[1] == self.num_latent
# deal with parameters for the prior mean variance of X
if X_prior_mean is None:
X_prior_mean = np.zeros((self.num_data, self.num_latent))
self.X_prior_mean = X_prior_mean
if X_prior_var is None:
X_prior_var = np.eye(self.num_latent) # For now I set the X prior variance to one
self.X_prior_var = X_prior_var
assert X_prior_mean.shape[0] == self.num_data
assert X_prior_mean.shape[1] == self.num_latent
assert X_prior_var.shape[0] == self.num_latent
assert X_prior_var.shape[1] == self.num_latent
@property
def _X_var(self):
if self.X_std.get_shape().ndims == 3:
return tf.matmul(self.X_std, self.X_std, transpose_b=True)
elif self.X_std.get_shape().ndims == 2:
return tf.square(self.X_std)
def build_likelihood(self):
"""
Construct a tensorflow function to compute the bound on the marginal
likelihood.
"""
num_inducing = tf.shape(self.Z)[0]
psi0 = tf.reduce_sum(self.kern.eKdiag(self.X_mean, self._X_var), 0)
psi1 = self.kern.eKxz(self.Z, self.X_mean, self._X_var)
psi2 = tf.reduce_sum(self.kern.eKzxKxz(self.Z, self.X_mean, self._X_var), 0)
Kuu = self.kern.K(self.Z) + eye(num_inducing) * settings.numerics.jitter_level
L = tf.cholesky(Kuu)
sigma2 = self.likelihood.variance
sigma = tf.sqrt(sigma2)
# Compute intermediate matrices
A = tf.matrix_triangular_solve(L, tf.transpose(psi1), lower=True) / sigma
tmp = tf.matrix_triangular_solve(L, psi2, lower=True)
AAT = tf.matrix_triangular_solve(L, tf.transpose(tmp), lower=True) / sigma2
B = AAT + eye(num_inducing)
LB = tf.cholesky(B)
log_det_B = 2. * tf.reduce_sum(tf.log(tf.diag_part(LB)))
c = tf.matrix_triangular_solve(LB, tf.matmul(A, self.Y), lower=True) / sigma
# KL[q(x) || p(x)]
if self.X_std.get_shape().ndims == 3:
KL = kullback_leiblers.gauss_kl(tf.transpose(self.X_prior_mean - self.X_mean),
tf.transpose(self.X_std, [1, 2, 0]), self.X_prior_var)
elif self.X_std.get_shape().ndims == 2:
KL = kullback_leiblers.gauss_kl_diag(tf.transpose(self.X_prior_mean - self.X_mean),
tf.transpose(self.X_std), self.X_prior_var)
# compute log marginal bound
D = tf.cast(tf.shape(self.Y)[1], float_type)
ND = tf.cast(tf.size(self.Y), float_type)
bound = -0.5 * ND * tf.log(2 * np.pi * sigma2)
bound += -0.5 * D * log_det_B
bound += -0.5 * tf.reduce_sum(tf.square(self.Y)) / sigma2
bound += 0.5 * tf.reduce_sum(tf.square(c))
bound += -0.5 * D * (tf.reduce_sum(psi0) / sigma2 -
tf.reduce_sum(tf.diag_part(AAT)))
bound -= KL
return bound
def build_predict(self, Xnew, full_cov=False):
"""
Compute the mean and variance of the latent function at some new points.
Note that this is very similar to the SGPR prediction, for which
there are notes in the SGPR notebook.
:param Xnew: Point to predict at.
"""
num_inducing = tf.shape(self.Z)[0]
psi1 = self.kern.eKxz(self.Z, self.X_mean, self._X_var)
psi2 = tf.reduce_sum(self.kern.eKzxKxz(self.Z, self.X_mean, self._X_var), 0)
Kuu = self.kern.K(self.Z) + eye(num_inducing) * 1e-6
Kus = self.kern.K(self.Z, Xnew)
sigma2 = self.likelihood.variance
sigma = tf.sqrt(sigma2)
L = tf.cholesky(Kuu)
A = tf.matrix_triangular_solve(L, tf.transpose(psi1), lower=True) / sigma
tmp = tf.matrix_triangular_solve(L, psi2, lower=True)
AAT = tf.matrix_triangular_solve(L, tf.transpose(tmp), lower=True) / sigma2
B = AAT + eye(num_inducing)
LB = tf.cholesky(B)
c = tf.matrix_triangular_solve(LB, tf.matmul(A, self.Y), lower=True) / sigma
tmp1 = tf.matrix_triangular_solve(L, Kus, lower=True)
tmp2 = tf.matrix_triangular_solve(LB, tmp1, lower=True)
mean = tf.matmul(tf.transpose(tmp2), c)
if full_cov:
var = self.kern.K(Xnew) + tf.matmul(tf.transpose(tmp2), tmp2) \
- tf.matmul(tf.transpose(tmp1), tmp1)
shape = tf.pack([1, 1, tf.shape(self.Y)[1]])
var = tf.tile(tf.expand_dims(var, 2), shape)
else:
var = self.kern.Kdiag(Xnew) + tf.reduce_sum(tf.square(tmp2), 0) \
- tf.reduce_sum(tf.square(tmp1), 0)
shape = tf.pack([1, tf.shape(self.Y)[1]])
var = tf.tile(tf.expand_dims(var, 1), shape)
return mean + self.mean_function(Xnew), var