Created
March 1, 2016 20:43
-
-
Save juancamilog/f7ba342435c767f2dacd to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import theano | |
import theano.tensor as T | |
from theano.sandbox.linalg import psd,matrix_inverse,det,cholesky | |
import numpy as np | |
from time import time | |
def SEard(loghyp,X): | |
''' Squared exponential kernel with diagonal scaling matrix (one lengthscale per dimension)''' | |
N,idims = X.shape | |
M = T.diag(T.exp(-2*loghyp[:idims])) | |
XM = X.dot(M); | |
sumXsq = T.sum(XM*X,1) | |
D = sumXsq[:,None] + sumXsq[None,:] - 2*XM.dot(X.T) | |
K = T.exp(2*loghyp[idims] - 0.5*D) + T.eye(N)*T.exp(2*loghyp[idims+1]) | |
return K | |
def nlml_(X,Y, loghyps): | |
N = X.shape[0] | |
K = SEard(loghyps,X) | |
iK = matrix_inverse(psd(K)) | |
beta = iK.dot(Y) | |
return 0.5*(Y.T.dot(beta) + T.log(det(psd(K))) + N*np.log(2*np.pi)) | |
def nlml_chol_(X,Y, loghyps): | |
N = X.shape[0] | |
K = SEard(loghyps,X) | |
cholK = cholesky(K) | |
Yc = matrix_inverse(cholK).dot(Y) # this should get optimized into a call to solve. Didn't call solve directly since it does not provide a gradient | |
return 0.5*(Yc.T.dot(Yc) + 2*T.sum(T.log(T.diag(cholK))) + N*np.log(2*np.pi)) | |
# generate some test data | |
np.set_printoptions(linewidth=200) | |
n_samples = 1000 | |
idims = 4 | |
def f(X): | |
return np.sin(X.sum(1))*np.cos(X.sum(1))*np.exp(-0.5*np.sum(X**2,1)) | |
Xd = np.random.multivariate_normal(np.zeros((idims,)),10*np.eye(idims),n_samples) | |
Yd = f(Xd) | |
# set some initial hyper parameters | |
lh = np.zeros((idims+2,)) | |
lh[:idims] = np.var(Xd,0,ddof=1) | |
lh[idims] = np.std(Yd,ddof=1) | |
lh[idims+1] = 0.1*np.std(Yd,ddof=1) | |
lh = np.log(lh) | |
X = theano.shared(Xd) | |
Y = theano.shared(Yd) | |
loghyp1 = theano.shared(lh, borrow=False) | |
loghyp2 = theano.shared(lh, borrow=False) | |
nlml = nlml_(X,Y,loghyp1) | |
dnlml = T.grad(nlml,loghyp1) | |
nlml_chol = nlml_chol_(X,Y,loghyp2) | |
dnlml_chol = T.grad(nlml_chol,loghyp2) | |
f1 = theano.function([],nlml, name='nlml') | |
f2 = theano.function([],nlml_chol, name='nlml_chol') | |
df1 = theano.function([],dnlml, name='dnlml') | |
df2 = theano.function([],dnlml_chol, name='dnlml_chol') | |
t_nlml_ = 0 | |
t_nlml_chol_ = 0 | |
n_tests = 100 | |
for i in xrange(n_tests): | |
print '--- %d ---'%(i) | |
# test nlml | |
start = time() | |
ret = df1() | |
t_nlml = time()-start | |
loghyp1.set_value(loghyp1.get_value()-1e-5*ret) | |
print f1(),ret,'Elapsed: %f'%(t_nlml) | |
t_nlml_ += t_nlml | |
# test nlml_chol | |
start = time() | |
ret = df2() | |
t_nlml_chol = time()-start | |
loghyp2.set_value(loghyp2.get_value()-1e-5*ret) | |
print f2(),ret,'Elapsed: %f'%(t_nlml_chol) | |
t_nlml_chol_ += t_nlml_chol | |
print 'Average for nlml using matrix inverse: %f'%(t_nlml_/n_tests) | |
print 'Average for nlml using cholesky: %f'%(t_nlml_chol_/n_tests) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment