apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-15 06:52:39 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

fix Linear kernel with SSGPLVM

Browse source
This commit is contained in:
Zhenwen Dai 2014-06-09 14:21:56 +01:00
parent 6119c33078
commit 695bad63c5
6 changed files with 33 additions and 117 deletions
Download patch file Download diff file Expand all files Collapse all files

2
GPy/core/gp.py
View file

@ -40,7 +40,7 @@ class GP(Model):
assert Y.ndim == 2
self.Y = ObsAr(Y)
assert Y.shape[0] == self.num_data
# assert Y.shape[0] == self.num_data
_, self.output_dim = self.Y.shape
#TODO: check the type of this is okay?

11
GPy/inference/latent_function_inference/posterior.py
View file

@ -153,9 +153,14 @@ class Posterior(object):
$$
"""
if self._woodbury_inv is None:
self._woodbury_inv, _ = dpotri(self.woodbury_chol, lower=1)
#self._woodbury_inv, _ = dpotrs(self.woodbury_chol, np.eye(self.woodbury_chol.shape[0]), lower=1)
symmetrify(self._woodbury_inv)
if self._woodbury_chol is not None:
self._woodbury_inv, _ = dpotri(self._woodbury_chol, lower=1)
#self._woodbury_inv, _ = dpotrs(self.woodbury_chol, np.eye(self.woodbury_chol.shape[0]), lower=1)
symmetrify(self._woodbury_inv)
elif self._covariance is not None:
B = self._K - self._covariance
tmp, _ = dpotrs(self.K_chol, B)
self._woodbury_inv, _ = dpotrs(self.K_chol, tmp.T)
return self._woodbury_inv
@property

6
GPy/kern/_src/linear.py
View file

@ -11,8 +11,8 @@ from ...core.parameterization import Param
from ...core.parameterization.transformations import Logexp
from ...util.caching import Cache_this
from ...core.parameterization import variational
from psi_comp import linear_psi_comp
from ...util.config import *
from .psi_comp import PSICOMP_Linear
class Linear(Kern):
"""
@ -53,10 +53,8 @@ class Linear(Kern):
self.variances = Param('variances', variances, Logexp())
self.add_parameter(self.variances)
self.psicomp = PSICOMP_Linear()
def set_for_SpikeAndSlab(self):
self.psicomp = linear_psi_comp.PSICOMP_SSLinear()
@Cache_this(limit=2)
def K(self, X, X2=None):
if self.ARD:

35
GPy/kern/_src/psi_comp/__init__.py
View file

@ -6,20 +6,7 @@ from GPy.util.caching import Cache_this
from ....core.parameterization import variational
import rbf_psi_comp
import ssrbf_psi_comp
class PSICOMP(Pickleable):
def psicomputations(self, variance, Z, variational_posterior):
"""
Compute psi-statistics
"""
pass
def psiDerivativecomputations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, Z, variational_posterior):
"""
Compute the derivatives of parameters by combing dL_dpsi and dpsi_dparam
"""
pass
import sslinear_psi_comp
class PSICOMP_RBF(Pickleable):
@ -38,5 +25,25 @@ class PSICOMP_RBF(Pickleable):
return rbf_psi_comp.psiDerivativecomputations(dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, lengthscale, Z, variational_posterior)
elif isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
return ssrbf_psi_comp.psiDerivativecomputations(dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, lengthscale, Z, variational_posterior)
else:
raise ValueError, "unknown distriubtion received for psi-statistics"
class PSICOMP_Linear(Pickleable):
@Cache_this(limit=1, ignore_args=(0,))
def psicomputations(self, variance, Z, variational_posterior):
if isinstance(variational_posterior, variational.NormalPosterior):
raise NotImplementedError
elif isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
return sslinear_psi_comp.psicomputations(variance, Z, variational_posterior)
else:
raise ValueError, "unknown distriubtion received for psi-statistics"
@Cache_this(limit=1, ignore_args=(0,1,2,3))
def psiDerivativecomputations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, Z, variational_posterior):
if isinstance(variational_posterior, variational.NormalPosterior):
raise NotImplementedError
elif isinstance(variational_posterior, variational.SpikeAndSlabPosterior):
return sslinear_psi_comp.psiDerivativecomputations(dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, Z, variational_posterior)
else:
raise ValueError, "unknown distriubtion received for psi-statistics"

93
GPy/kern/_src/psi_comp/linear_psi_comp.py
View file

@ -1,93 +0,0 @@
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
"""
The package for the Psi statistics computation of the linear kernel for SSGPLVM
"""
import numpy as np
from . import PSICOMP
from GPy.util.caching import Cache_this
class PSICOMP_SSLinear(PSICOMP):
@Cache_this(limit=1, ignore_args=(0,))
def psicomputations(self, variance, Z, variational_posterior):
"""
Compute psi-statistics for ss-linear kernel
"""
# here are the "statistics" for psi0, psi1 and psi2
# Produced intermediate results:
# psi0 N
# psi1 NxM
# psi2 MxM
mu = variational_posterior.mean
S = variational_posterior.variance
gamma = variational_posterior.binary_prob
psi0 = np.einsum('q,nq,nq->n',variance,gamma,np.square(mu)+S)
psi1 = np.einsum('nq,q,mq,nq->nm',gamma,variance,Z,mu)
mu2 = np.square(mu)
variances2 = np.square(variance)
tmp = np.einsum('nq,q,mq,nq->nm',gamma,variance,Z,mu)
psi2 = np.einsum('nq,q,mq,oq,nq->mo',gamma,variances2,Z,Z,mu2+S)+\
np.einsum('nm,no->mo',tmp,tmp) - np.einsum('nq,q,mq,oq,nq->mo',np.square(gamma),variances2,Z,Z,mu2)
return psi0, psi1, psi2
@Cache_this(limit=1, ignore_args=(0,1,2,3))
def psiDerivativecomputations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, variance, Z, variational_posterior):
mu = variational_posterior.mean
S = variational_posterior.variance
gamma = variational_posterior.binary_prob
dL_dvar, dL_dgamma, dL_dmu, dL_dS, dL_dZ = self._psi2computations(dL_dpsi2, variance, Z, mu, S, gamma)
# Compute for psi0 and psi1
mu2S = np.square(mu)+S
dL_dvar += np.einsum('n,nq,nq->q',dL_dpsi0,gamma,mu2S) + np.einsum('nm,nq,mq,nq->q',dL_dpsi1,gamma,Z,mu)
dL_dgamma += np.einsum('n,q,nq->nq',dL_dpsi0,variance,mu2S) + np.einsum('nm,q,mq,nq->nq',dL_dpsi1,variance,Z,mu)
dL_dmu += np.einsum('n,nq,q,nq->nq',dL_dpsi0,gamma,2.*variance,mu) + np.einsum('nm,nq,q,mq->nq',dL_dpsi1,gamma,variance,Z)
dL_dS += np.einsum('n,nq,q->nq',dL_dpsi0,gamma,variance)
dL_dZ += np.einsum('nm,nq,q,nq->mq',dL_dpsi1,gamma, variance,mu)
return dL_dvar, dL_dZ, dL_dmu, dL_dS, dL_dgamma
def _psi2computations(self, dL_dpsi2, variance, Z, mu, S, gamma):
"""
Z - MxQ
mu - NxQ
S - NxQ
gamma - NxQ
"""
# here are the "statistics" for psi1 and psi2
# Produced intermediate results:
# _psi2_dvariance Q
# _psi2_dZ MxQ
# _psi2_dgamma NxQ
# _psi2_dmu NxQ
# _psi2_dS NxQ
mu2 = np.square(mu)
gamma2 = np.square(gamma)
variance2 = np.square(variance)
mu2S = mu2+S # NxQ
common_sum = np.einsum('nq,q,mq,nq->nm',gamma,variance,Z,mu) # NxM
dL_dvar = np.einsum('mo,nq,q,mq,oq->q',dL_dpsi2,2.*(gamma*mu2S-gamma2*mu2),variance,Z,Z)+\
np.einsum('mo,nq,mq,nq,no->q',dL_dpsi2,gamma,Z,mu,common_sum)+\
np.einsum('mo,nq,oq,nq,nm->q',dL_dpsi2,gamma,Z,mu,common_sum)
dL_dgamma = np.einsum('mo,q,mq,oq,nq->nq',dL_dpsi2,variance2,Z,Z,(mu2S-2.*gamma*mu2))+\
np.einsum('mo,q,mq,nq,no->nq',dL_dpsi2,variance,Z,mu,common_sum)+\
np.einsum('mo,q,oq,nq,nm->nq',dL_dpsi2,variance,Z,mu,common_sum)
dL_dmu = np.einsum('mo,q,mq,oq,nq,nq->nq',dL_dpsi2,variance2,Z,Z,mu,2.*(gamma-gamma2))+\
np.einsum('mo,nq,q,mq,no->nq',dL_dpsi2,gamma,variance,Z,common_sum)+\
np.einsum('mo,nq,q,oq,nm->nq',dL_dpsi2,gamma,variance,Z,common_sum)
dL_dS = np.einsum('mo,nq,q,mq,oq->nq',dL_dpsi2,gamma,variance2,Z,Z)
dL_dZ = 2.*(np.einsum('om,nq,q,mq,nq->oq',dL_dpsi2,gamma,variance2,Z,mu2S)+np.einsum('om,nq,q,nq,nm->oq',dL_dpsi2,gamma,variance,mu,common_sum)
-np.einsum('om,nq,q,mq,nq->oq',dL_dpsi2,gamma2,variance2,Z,mu2))
return dL_dvar, dL_dgamma, dL_dmu, dL_dS, dL_dZ

3
GPy/util/linalg.py
View file

@ -16,8 +16,7 @@ import warnings
import os
from config import *
if np.all(np.float64((scipy.__version__).split('.')[:2]) >= np.array([0, 12])):
#import scipy.linalg.lapack.clapack as lapack
if float('.'.join((scipy.__version__).split('.')[:2])) >= 0.12:
from scipy.linalg import lapack
else:
from scipy.linalg.lapack import flapack as lapack