I’ll try to give examples of efficient gaussian process computation here, like the vec trick (Kronecker product trick), efficient toeliptz and circulant matrix computations, RTS smoothing and Kalman filtering using state space representations, and so on.

Fast Toeplitz Matrix-Vector Products and Solving

Circulant matrix-vector products are super-fast due to the way circulant matrices can be decomposed using their fourier transforms. Toeplitz matrices can be ‘embedded’ into a circulant matrix and their matrix-vector products can be computed efficiently too. I’ll add details later. Sample code is shown below.

This can then be used to compute multivariate log-densities quickly. Furthermore, this multiplication can then be used in conjugate gradient solvers to efficiently compute inverse-matrix-vector products in $O(n \log n)$ time and $O(n)$ space.

Python Code
import numpy as np
from scipy.linalg import toeplitz

def toeplitz_matmul(c, r = None, v = None):
	if v is None: v = np.zeros(c.shape)
	if r is None: r = c.conj()
        if len(v.shape) == 1: v = v.reshape(-1, 1)

	n = c.shape[0]; assert c.shape == r.shape
	embeded_col = np.hstack((c, np.flip(r[1:])))
	padded_v = np.vstack((v, np.zeros((n - 1, v.shape[1]))))
	fft_T = np.fft.fft(embeded_col, axis = 0).reshape(-1, 1)
	fft_v = np.fft.fft(padded_v, axis = 0)
	return np.fft.ifft(fft_T*fft_v, axis = 0).real[:n, :]

Toeplitz Matrix Cholesky Decomposition

… and also circulant matrix solving in the comments (using scipy and ctypes).

I got the toeplitz_cholesky library from here and compiled it. I’m going to check out toeblitz in the future.

C/Python Code
# in bash: clang -shared -fpic toeplitz_cholesky.c -o tc.dylib -O3

import time
import numpy as np
import matplotlib.pyplot as plt

plt.style.use("ggplot"); plt.ion()
# sc.linalg.solve_circulant, append t[(n - 2):0:-1]

dll = np.ctypeslib.load_library("tc", ".")
np_poin = np.ctypeslib.ndpointer

def kernel(n = 100):
	k = np.linspace(0, 5, n)
	k = np.exp(-(k - k[0])**2)
	k[0] += 1e-10
	return k

def toep_chol_prepare(n):
	type_input_1 = np.ctypeslib.ctypes.c_int64
	type_input_2 = np_poin(dtype = np.double, ndim = 1, shape = n)
	type_output = np_poin(dtype = np.double, ndim = 2, shape = (n, n))

	dll.t_cholesky_lower.argtypes = [type_input_1, type_input_2]
	dll.t_cholesky_lower.restype = type_output
	return dll.t_cholesky_lower

def timer(n = 100):
	func_ptr = toep_chol_prepare(n)
	tic = time.time()
	L = func_ptr(n, kernel(n))
	toc = time.time()
	return toc - tic, L

if __name__ == "__main__":
	runtime, L = timer(10000) # half a second!
	x = np.matmul(L.T, np.random.normal(size = 10000))
	input("Press the enter key to quit.")

Efficient Inference of GP Covariances

Sometimes, it is easier to do parameter inference in the frequency domain. Here, I use SymPy to get the theoretical spectrum of a GP (using a Fourier transform of the covariance - note that to get from samples to the PSD, the PSD is defined as the expected value of the series squared due to Wiener-Khinchin) and we use the fact that the empirical spectrum divided by the theoretical spectrum has an \(Exp(1)\) distribution (\(\chi^2_2 \stackrel{d}{=} 0.5Exp(0.5)\)) to get to the likelihood.

I was writing a Gaussian Process Vocoder that synthesizes speech from mel spectrograms (using an LSTM to get from the mel spectrograms to the spectral kernel’s parameters), but the whole thing looks too similar to Tokuda & Zen (Directly Modelling Speech Waveforms …) - which I discovered after writing a good chuck of the code. I might complete it at some point, it uses the spectral kernel to get a zero mean GP of the right frequencies, another GP for amplitude modulation and block-stationary treatments of the non-stationary GP (so synthesis also happens blockwise, each block is conditioned on the previous one).

Sample python Code
import sympy as sy
import numpy as np
from tqdm import tqdm
import tensorflow as tf
import matplotlib.pyplot as plt
from scipy.io import wavfile as wav
from scipy.linalg import toeplitz

md = sy.functions.Abs
sess = tf.InteractiveSession()
plt.style.use("ggplot"); plt.ion()

def spectrum(x, u = 1, return_freq = True):
	n = len(x)
	intm = np.fft.fft(x)
	intm = (intm * intm.conjugate()).real/n
	intm = intm * u/n 
	psd = intm[1:int(np.ceil(0.5 * (n + 1)))]
	if return_freq:
		freq = range(1, len(psd) + 1) 
		freq = np.array(freq, dtype = float)
		freq *= 0.5*n / (u * len(psd))
		return freq, psd
		return psd

class SimulationData:
	def __init__(self):
		self.num_data = None
		self.time = None
		self.upper_lim = None
		self.kernel_func = None
		self.cholesky_fac = None
	def gen_grid(self, num_data = 500, upper_lim = 5):
		self.num_data = num_data
		self.upper_lim = upper_lim
		self.time = np.linspace(0, upper_lim, num_data)
	def gen_chol(self, kernel_func = None):
		if kernel_func is None:
			kernel_func = self.kernel_func
		C = toeplitz(kernel_func(self.time))
		C += np.identity(len(self.time))*1e-7
		C = np.linalg.cholesky(C)
		self.cholesky_fac = C
	def simulate(self):
		C = self.cholesky_fac
		z = np.random.normal(size = self.num_data)
		return C @ z
	def gen_expcos_kern(self):
		sm = dict(r = sy.symbols("r", real = True),
			w = sy.symbols("w", real = True, positive = True),
			l = sy.symbols("l", real = True, positive = True),
			p = sy.symbols("p", real = True, positive = True),
			s = sy.symbols("s", real = True, positive = True))
		self.symbols = sm

def spectrum_check(l = 0.5, s = 2.0, p = 0.5, n = 10000):
	r, w = sy.symbols("r, w", real = True)
	K = s**2 * sy.exp(-(r/l)**2) * sy.cos(2*p*sy.pi*r)
	S = sy.fourier_transform(K, r, w)
	S = sy.lambdify(w, S)
	K = sy.lambdify(r, K)

	data = SimulationData()

	freq, psd = spectrum(data.simulate(), data.upper_lim)
	for i in range(n - 1):
		psd += spectrum(data.simulate(), data.upper_lim)[1]

	s_t = S(data.time)
	plt.plot(data.time, s_t, label = "Theoretical")
	plt.xlim(-0.5, data.upper_lim)
	plt.plot(freq, psd/n, label = "Empirical")

if __name__ == "__main__":

	spectrum_check() # To see if my spectrum function is correct

	l, s, p, u = 0.5, 2.0, 0.5, 5.0
	r, w = sy.symbols("r, w", real = True)
	s_p, p_p, l_p = sy.symbols("s, p, l", real = True, positive = True)

	K = s**2 * sy.exp(-(r/l)**2) * sy.cos(2*p*sy.pi*r)
	K_t = s_p**2 * sy.exp(-(r/l_p)**2) * sy.cos(2*p_p*sy.pi*r)
	S = sy.fourier_transform(K_t, r, w)
	S = sy.lambdify([w, l_p, p_p, s_p], S, "tensorflow")
	K = sy.lambdify(r, K)

	data = SimulationData()
	data.gen_grid(upper_lim = u)

	l_param = tf.clip_by_value(tf.Variable(1.1), 1e-5, 1e5)
	p_param = tf.clip_by_value(tf.Variable(2.2), 1e-5, 1e5)
	s_param = tf.clip_by_value(tf.Variable(0.9), 1e-5, 1e5)

	freq_obs, spec_obs = spectrum(data.simulate(), u)
	freq_obs = np.array(freq_obs, dtype = "float32")
	spec_obs = np.array(spec_obs, dtype = "float32")

	spec_const = tf.constant(spec_obs)

	add_term = S(freq_obs, l_param, p_param, s_param) + 1e-10
	loss = tf.reduce_sum((spec_const - add_term)**2) # likelihood is too unstable
	step = tf.train.AdamOptimizer().minimize(loss)


	for i in tqdm(range(10000)):
		if i % 50 == 0:
			freq_obs, spec_obs = spectrum(data.simulate(), u)
			freq_obs = np.array(freq_obs, dtype = "float32")
			spec_obs = np.array(spec_obs, dtype = "float32")
		sess.run(step, feed_dict = {spec_const: spec_obs})
	print(np.round(sess.run([l_param, p_param, s_param]), 2))
	print(np.round([l, p, s], 2))


Efficient Gaussian Process Computation

I’ll try to give examples of efficient gaussian process computation here, like the vec trick (Kronecker product trick), efficient toeliptz and circulant matrix computations, RTS smoothing and Kalman filtering using state space representations, and so on.

3 min read

Gaussian Process Speech Synthesis (Draft)

Very untidy first working draft of the idea mentioned on the efficient computation page. Here, I fit a spectral mixture to some audio data to build a “generative model” for audio. I’ll implement efficient sampling later, and I’ll replace the arbitrary way this is trained with an LSTM-RNN to go straight from text/spectrograms to waveforms.

6 min read
Back to Top ↑


Gaussian Process Middle C

First of my experiments on audio modelling using gaussian processes. Here, I construct a GP that, when sampled, plays middle c the way a grand piano would.

1 min read
Back to Top ↑


Back to Top ↑