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 can be 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. Sample code is shown below.

Main reference:

J. Dongarra, P. Koev, and X. Li, “Matrix-Vector and Matrix-Matrix Multiplications”. In Z. Bai, J. Demmel, J. Dongarra, A. Ruhe, and H. van der Vorst, editors. “Templates for the Solution of Algebraic Eigenvalue Problems: A Practical Guide”. SIAM, Philadelphia, 2000. Available online.

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

I’m trying to contribute some code to SciPy, pull request can be found here.

Toeplitz matrices can be inverted in quadratic time using the Levinson algorithm, that also has the capability of returning a determinant as a side product using no more computation (as I understand).

Some code that calculates this log determinant can be found here:

Marano S, Edwards B, Ferrari G and Fah D (2017), “Fitting Earthquake Spectra: Colored Noise and Incomplete Data”, Bulletin of the Seismological Society of America., January, 2017. Vol. 107(1), pp. 276-291. Implementation available online..

The fast matrix-vector multiplication can be used with a preconditioned conjugate gradient algorithm to implement fast inverse-vector products.

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. I believe that the following algorithm calculates the cholesky decomposition in $$O(n^2)$$ time.

C/Python Code

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 obtain a zero mean GP with the right frequencies and block-stationary treatments of these non-stationary GPs to synthesize audio. Tokuda & Zen used an LSTM-RNN to obtain the cepstral coefficients, while I used it to obtain spectral parameters.

Have a look at Gaussian Process Speech Synthesis (Draft) for code.

Sampling Periodic GPs using the State Space Representation

The representation is due to Solin & Sarkka (2014).

Python Code

Efficient inverse-matrix multiplication

2020

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.

Minimal Working Examples

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.

Random Projects

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.

Consider: