amdlibm_trig_op_seglfault.py

import theano
import theano.tensor as T
from theano.sandbox.linalg import psd,matrix_inverse,det,cholesky
from theano.tensor.slinalg import solve_lower_triangular
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].astype(theano.config.floatX)
    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*T.log(2*np.pi)).astype(theano.config.floatX)

def nlml_ss_(X,Y, loghyps,n_basis=100):
    # sample initial unscaled spectral points
    w_ = np.random.randn(n_basis,idims).astype(theano.config.floatX)
    w = theano.shared(w_,borrow=True)
    sr = w*T.exp(-loghyps[:idims])

    N = X.shape[0].astype(theano.config.floatX)
    M = sr.shape[0].astype(theano.config.floatX)
    sf2 = T.exp(2*loghyps[idims])
    sn2 = T.exp(2*loghyps[idims+1])
    # sr.T.dot(x) for all sr ( n_basis x *D) and X (N x D). size n_basis x N
    srdotX = sr.dot(X.T)
    # convert to sin cos
    phi_f = T.vertical_stack(T.sin(srdotX), T.cos(srdotX))
     
    A = phi_f.dot(phi_f.T) + (M*sn2/sf2)*T.eye(2*sr.shape[0])
    iA = matrix_inverse(psd(A))
    beta_ss = iA.dot(phi_f.dot(Y))
    nlml_ss = 0.5*( Y.dot(Y) - phi_f.dot(Y).dot(beta_ss) )/sn2 + 0.5*T.log(det(psd(A))) - M*T.log((M*sn2/sf2)) + 0.5*N*T.log(2*np.pi*sn2)
    return nlml_ss
    
# generate some test data
np.set_printoptions(linewidth=200)
n_samples = 500
idims = 4
def f(X):
    return 10*np.sin(X.sum(1))*np.cos(X.sum(1))*np.exp(-0.005*np.sum(X**2,1))
Xd = np.random.multivariate_normal(np.zeros((idims,)),10*np.eye(idims),n_samples).astype(theano.config.floatX)
Yd = f(Xd).astype(theano.config.floatX)

# 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).astype(theano.config.floatX)

X = theano.shared(Xd)
Y = theano.shared(Yd)
loghyp1 = theano.shared(lh, borrow=False)

#nlml =  nlml_(X,Y,loghyp1)   <-- this one works with amdlibm
nlml =  nlml_ss_(X,Y,loghyp1) # <-- this crashes with amdlibm

f1 = theano.function([],nlml, name='nlml')
print f1()