TF Cholesky fail on RBF pos-def matrix

cholesky-issues.py

import gpflow as gp
import numpy as np
import pylab as plt
import tensorflow as tf

X = np.array([[  1.16527441e+09],
 [  1.16527442e+09],
 [  1.16527443e+09],
 [  1.16527443e+09],
 [  1.16527444e+09],
 [  1.16527445e+09],
 [  1.16527446e+09],
 [  1.16527447e+09],
 [  1.16527447e+09],
 [  1.16527448e+09],
 [  1.16527449e+09],
 [  1.16527450e+09],
 [  1.16527451e+09],
 [  1.16527451e+09],
 [  1.16527452e+09],
 [  1.16527453e+09],
 [  1.16527454e+09],
 [  1.16527455e+09],
 [  1.16527455e+09],
 [  1.16527456e+09],
 [  1.16527457e+09],
 [  1.16527458e+09],
 [  1.16527459e+09],
 [  1.16527459e+09],
 [  1.16527460e+09],
 [  1.16527461e+09],
 [  1.16527462e+09],
 [  1.16527463e+09],
 [  1.16527463e+09],
 [  1.16527464e+09],
 [  1.16527465e+09],
 [  1.16527466e+09],
 [  1.16527467e+09],
 [  1.16527467e+09],
 [  1.16527468e+09],
 [  1.16527469e+09],
 [  1.16527470e+09],
 [  1.16527471e+09],
 [  1.16527471e+09],
 [  1.16527472e+09],
 [  1.16527473e+09],
 [  1.16527474e+09],
 [  1.16527475e+09],
 [  1.16527475e+09],
 [  1.16527476e+09],
 [  1.16527477e+09],
 [  1.16527478e+09],
 [  1.16527479e+09],
 [  1.16527479e+09],
 [  1.16527480e+09],
 [  1.16527481e+09],
 [  1.16527482e+09],
 [  1.16527483e+09],
 [  1.16527483e+09],
 [  1.16527484e+09],
 [  1.16527485e+09],
 [  1.16527486e+09],
 [  1.16527487e+09],
 [  1.16527487e+09],
 [  1.16527488e+09],
 [  1.16527489e+09],
 [  1.16527490e+09],
 [  1.16527491e+09],
 [  1.16527491e+09],
 [  1.16527492e+09],
 [  1.16527493e+09],
 [  1.16527494e+09],
 [  1.16527495e+09],
 [  1.16527495e+09],
 [  1.16527496e+09],
 [  1.16527497e+09],
 [  1.16527498e+09],
 [  1.16527499e+09],
 [  1.16527499e+09],
 [  1.16527500e+09],
 [  1.16527501e+09],
 [  1.16527502e+09],
 [  1.16527503e+09],
 [  1.16527503e+09],
 [  1.16527504e+09],
 [  1.16527505e+09],
 [  1.16527506e+09],
 [  1.16527507e+09],
 [  1.16527507e+09],
 [  1.16527508e+09],
 [  1.16527509e+09],
 [  1.16527510e+09],
 [  1.16527511e+09],
 [  1.16527511e+09],
 [  1.16527512e+09],
 [  1.16527513e+09],
 [  1.16527514e+09],
 [  1.16527515e+09],
 [  1.16527516e+09],
 [  1.16527516e+09],
 [  1.16527517e+09],
 [  1.16527518e+09],
 [  1.16527519e+09],
 [  1.16527520e+09],
 [  1.16527520e+09],
 [  1.16527521e+09],
 [  1.16527522e+09],
 [  1.16527523e+09],
 [  1.16527524e+09],
 [  1.16527524e+09],
 [  1.16527525e+09],
 [  1.16527526e+09],
 [  1.16527527e+09],
 [  1.16527528e+09],
 [  1.16527528e+09],
 [  1.16527529e+09],
 [  1.16527530e+09],
 [  1.16527531e+09],
 [  1.16527532e+09],
 [  1.16527532e+09],
 [  1.16527533e+09],
 [  1.16527534e+09],
 [  1.16527535e+09],
 [  1.16527536e+09],
 [  1.16527536e+09]])
Y = np.array([[ 2.46963311],
 [ 2.81275076],
 [ 2.69000406],
 [ 2.58296666],
 [ 2.31684845],
 [ 2.18883278],
 [ 2.25738763],
 [ 2.39459771],
 [ 2.12522259],
 [ 2.49173322],
 [ 2.78297639],
 [ 2.65098281],
 [ 2.57148568],
 [ 2.77789879],
 [ 2.56483398],
 [ 2.53930404],
 [ 2.45846044],
 [ 2.08432404],
 [ 2.17751587],
 [ 2.06752078],
 [ 2.0319239 ],
 [ 1.89906843],
 [ 1.96754216],
 [ 1.69002888],
 [ 2.14023237],
 [ 2.14404539],
 [ 2.13234382],
 [ 1.82146425],
 [ 1.79754275],
 [ 1.87372486],
 [ 2.21795348],
 [ 2.2174718 ],
 [ 2.24379265],
 [ 2.3705007 ],
 [ 2.21843136],
 [ 2.12549238],
 [ 2.02982958],
 [ 2.11328605],
 [ 2.21238481],
 [ 2.21103339],
 [ 2.97721817],
 [ 2.67457913],
 [ 2.72743356],
 [ 2.68134139],
 [ 3.13029503],
 [ 3.43839799],
 [ 3.29713706],
 [ 3.46492849],
 [ 3.67350701],
 [ 3.88400577],
 [ 3.82812032],
 [ 4.05832468],
 [ 4.43553132],
 [ 4.41363113],
 [ 4.41936664],
 [ 4.70351457],
 [ 4.30679838],
 [ 4.00992616],
 [ 4.06255067],
 [ 4.06596992],
 [ 3.75273823],
 [ 3.97826996],
 [ 3.82667734],
 [ 4.0148188 ],
 [ 3.8492277 ],
 [ 4.20674274],
 [ 4.37250744],
 [ 4.2806856 ],
 [ 4.19154178],
 [ 4.11226701],
 [ 3.93415102],
 [ 3.88161289],
 [ 3.6367137 ],
 [ 3.51130442],
 [ 3.58040664],
 [ 3.38126723],
 [ 3.48296338],
 [ 3.37396197],
 [ 3.1863685 ],
 [ 3.49161607],
 [ 2.98974131],
 [ 3.2761566 ],
 [ 3.37750468],
 [ 3.08078927],
 [ 3.02741803],
 [ 3.18978172],
 [ 3.39560613],
 [ 3.6608262 ],
 [ 3.39282486],
 [ 3.26884529],
 [ 3.31234319],
 [ 3.43635524],
 [ 3.38701946],
 [ 3.2540391 ],
 [ 3.15468288],
 [ 3.15481631],
 [ 3.52582369],
 [ 2.88244772],
 [ 3.22489644],
 [ 3.43033441],
 [ 3.47408322],
 [ 3.70102469],
 [ 3.83587249],
 [ 3.90001012],
 [ 4.01125156],
 [ 3.76123542],
 [ 3.82831775],
 [ 4.42422725],
 [ 4.46235979],
 [ 4.69925034],
 [ 5.01875834],
 [ 5.18805289],
 [ 5.54333723],
 [ 5.4755843 ],
 [ 5.42707204],
 [ 5.75181478],
 [ 5.72793423],
 [ 6.01085997],
 [ 6.11502704],
 [ 6.41143633]])

plt.plot(X.flatten(),Y.flatten())
plt.show()

K_np = np.exp(-(X - X.T)**2/200.**2/2.) + np.eye(X.shape[0])*0.01
L_np = np.linalg.cholesky(K_np)

K_tf = tf.constant(K_np[None,:,:])
L_tf = tf.cholesky(K_tf)
sess = tf.Session()
L_tf = sess.run(L_tf)
assert np.allclose(L_tf, L_np)
#The difference of ~1e-14 may be enough to cause divergent behaviour 
#between numpy and tf

with gp.defer_build():
    k = gp.kernels.RBF(1,lengthscales=[50]) 
    k.lengthscales.transform = gp.transforms.Logistic(20,500)
    white = gp.kernels.White(1,variance=0.01)
    white.variance.set_trainable(False)
    kern = k + white
    mean = gp.mean_functions.Linear()
    m = gp.models.GPR(X, Y, kern=kern, mean_function=mean)
    m.compile()
    
o = gp.train.ScipyOptimizer(method='BFGS')

print(o.minimize(m,maxiter=1000))
print(m)
ystar,varstar = m.predict_f(X)
plt.plot(X.flatten(),Y.flatten())
plt.plot(X.flatten(),ystar.flatten())
plt.show()