Gaussian Process regression of a 3D scalar function

Questions
1

Hi all,

I am new to Pymc3 and Gaussian Processes so I apologize if this question seems trivial or I am missing something fundamental. I am trying to model a scalar function with multiple inputs. I have successfully done this using Pymc3 for a 2d toy function ( f(x,y) where x and y are coordinates ). However, I am unable to extend my model for a function with three inputs. I am trying to model now a toy function of the form f(x,y,z)=cos(x)+cos(y)+cos(z) + 2*(x2 + y2 + z**2). The code I am using is:

import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt

import theano
import theano.tensor.nlinalg
import theano.tensor as tt
import pymc3 as pm

#Define Test Function

def true_funct(x,y,z):
val = np.cos(x) + np.cos(y) + np.cos(z) + 2*(x2 + y2 + z**2)

return val

Create Observations

n_points = 40
noise_true = 0.001
data_x = np.random.uniform(-1,1,n_points)
data_y = np.random.uniform(-1,1,n_points)
data_z = np.random.uniform(-1,1,n_points)

X_train = np.concatenate((data_x[:,None],data_y[:,None],data_z[:,None]),axis=1)
Y_true = true_funct(data_x,data_y,data_z)[:,None]
Y_train = Y_true + np.random.normal(0,noise_true, Y_true.shape[0])[:,None]

sigma_noise = noise_true

with pm.Model() as model:

#List of Prior parameters
#-------------------------------

l_scale = pm.HalfCauchy('l_scale', 1, shape=3) #length scale parameter
sigma_f = pm.HalfCauchy('sigma_f', 1) #amplitude parameter
sigma_noise = pm.HalfCauchy('sigma_noise', .1)
  
#build the mean and cov matrix
cov_funct=pm.gp.cov.ExpQuad(3,l_scale)*sigma_f**2.
cov_funct=cov_funct(X_train)
K = tt.fill_diagonal(cov_funct, sigma_f**2. + sigma_noise**2.)
mean_func = pm.gp.mean.Zero()

obs = pm.MvNormal('obs', mean_func(X_train), K, observed=Y_train)
trace = pm.sample(1000,chains=1,cores=1, step=pm.NUTS())

pm.traceplot(trace, var_names=[‘sigma_f’,‘l_scale’]);

The issue I am having is that the length scales I get are extremely large and I pretty much get a constant response when sampling the GP at new locations. I have used successfully the above code for a function with: f(x,y) = sin(norm(x,y)) but I am having difficulty now with the three input case and I do not see why.

Thank you very much in advance for your time and apologies if this has been asked before.

2

I wanted to do something similar to you, so I used your example (reduced to 2D), but I can’t get either your example or my reduced version to work. I’m using theano 1.0.5 & pymc 3.9.3 on Linux (conda-forge install).

import numpy as np
from numpy.linalg import inv
import matplotlib.pyplot as plt

import theano
import theano.tensor.nlinalg
import theano.tensor as tt
import pymc3 as pm

#Define Test Function

def true_funct(x,y):
    val = np.cos(x) + np.cos(y) + 2*(x*x + y*y)
    return val

# Create Observations

n_points = 40
noise_true = 0.001
data_x = np.random.uniform(-1,1,n_points)
data_y = np.random.uniform(-1,1,n_points)

X_train = np.concatenate((data_x[:,None],data_y[:,None]),axis=1)
Y_true = true_funct(data_x,data_y)[:,None]
Y_train = Y_true + np.random.normal(0,noise_true, Y_true.shape[0])[:,None]
sigma_noise = noise_true

with pm.Model() as model:
    #List of Prior parameters
    #-------------------------------

    l_scale = pm.HalfCauchy('l_scale', 1, shape=2) #length scale parameter
    sigma_f = pm.HalfCauchy('sigma_f', 1) #amplitude parameter
    sigma_noise = pm.HalfCauchy('sigma_noise', .1)
    
    # build the mean and cov matrix
    cov_funct=pm.gp.cov.ExpQuad(2,l_scale)*sigma_f**2.
    cov_funct=cov_funct(X_train)
    K = tt.fill_diagonal(cov_funct, sigma_f**2. + sigma_noise**2.)
    mean_func = pm.gp.mean.Zero()

    obs = pm.MvNormal('obs', mean_func(X_train), K, observed=Y_train)
    trace = pm.sample(1000,chains=2,cores=2, step=pm.NUTS())

pm.traceplot(trace, var_names=["sigma_f","l_scale"]);

This code basically stops with an error:

~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/step_methods/hmc/quadpotential.py in raise_ok()
    262                               ' is zero.'.format(*name_slc[ii]))
--> 263             raise ValueError('\n'.join(errmsg))
    264 

ValueError: Mass matrix contains zeros on the diagonal. 
The derivative of RV `sigma_f_log__`.ravel()[0] is zero.

The above exception was the direct cause of the following exception:

RuntimeError                              Traceback (most recent call last)
<ipython-input-2-44d3d3d36dd9> in <module>
     41 
     42     obs = pm.MvNormal('obs', mean_func(X_train), K, observed=Y_train)
---> 43     trace = pm.sample(1000,chains=2,cores=2, step=pm.NUTS())
     44 
     45 pm.traceplot(trace, var_names=["sigma_f","l_scale"]);

~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/sampling.py in sample(draws, step, init, n_init, start, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, return_inferencedata, idata_kwargs, mp_ctx, pickle_backend, **kwargs)                                                       
    543         _print_step_hierarchy(step)
    544         try:
--> 545             trace = _mp_sample(**sample_args, **parallel_args)
    546         except pickle.PickleError:
    547             _log.warning("Could not pickle model, sampling singlethreaded.")

~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/sampling.py in _mp_sample(draws, tune, step, chains, cores, chain, random_seed, start, progressbar, trace, model, callback, discard_tuned_samples, mp_ctx, pickle_backend, **kwargs)
   1484         try:
   1485             with sampler:
-> 1486                 for draw in sampler:
   1487                     trace = traces[draw.chain - chain]
   1488                     if trace.supports_sampler_stats and draw.stats is not None:

~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/parallel_sampling.py in __iter__(self)
    490 
    491         while self._active:
--> 492             draw = ProcessAdapter.recv_draw(self._active)
    493             proc, is_last, draw, tuning, stats, warns = draw
    494             self._total_draws += 1

~/miniconda/envs/python3/lib/python3.7/site-packages/pymc3/parallel_sampling.py in recv_draw(processes, timeout)
    363             else:
    364                 error = RuntimeError("Chain %s failed." % proc.chain)
--> 365             raise error from old_error
    366         elif msg[0] == "writing_done":
    367             proc._readable = True

RuntimeError: Chain 0 failed.