Without any specific optimization for MNIST it can achieve a score of 97.89 percent.

This is a small, dependency-free Odin implementation of a multi-output Stochastic Variational Gaussian Process ( SVGP ) for f64 inputs and f64 outputs.

It implements:

• RBF, Matérn 3/2 and Matérn 5/2 kernels.
• Whitened inducing-variable SVGP: q( v ) = N( mu, diag( exp( log_s ) ) ).
• Stochastic mini-batch training with the Gaussian-likelihood ELBO. Is the core mathematical "loss function" that your model is trying to maximize during training. ELBO stands for Evidence Lower Bound. When you couple it with a Gaussian likelihood, it means you are assuming that your target data (the outputs) contains continuous, normally distributed noise.
• Adam optimizer, gradient clipping and numerical Cholesky jitter.
• CPU parallel feature computation, gradient accumulation and inference over all available processor cores.
• MNIST example with gzip or uncompressed IDX files.

Important limitation: the MNIST example treats classification as multi-output Gaussian regression on smoothed one-hot labels, then takes argmax. This is intentionally simple and fast. True GP classification would require a non-Gaussian likelihood and quadrature/variational bounds.

svgp/svgp.odin             reusable SVGP package
examples/mnist/main.odin   MNIST CLI example

From the project root:

odin run examples/mnist -- ../data

For speed:

odin build examples/mnist -o:speed -microarch:native -out:mnist_svgp
./mnist_svgp ../data --inducing 512 --epochs 10 --batch 128 --lr 0.002

../data should contain either these MNIST DataSet files:

train-images-idx3-ubyte
train-labels-idx1-ubyte
t10k-images-idx3-ubyte
t10k-labels-idx1-ubyte

or their .gz versions.

This version keeps the same public API, f64 arithmetic, kernels, optimizer, and MNIST example behavior, but removes several generic hot-path costs:

• RBF kernel evaluation no longer computes an unused square root.
• Kernel vector construction switches on the kernel once per row instead of once per inducing point.
• Kzz construction is parallelized across worker threads before LAPACK Cholesky.
• Triangular solves can use OpenBLAS dtrsv_ for larger inducing counts.
• Training fuses phi computation and gradient accumulation, so each batch creates one wave of workers instead of separate phi and gradient waves.
• exp( q_log_s ) is cached once per batch instead of recomputed inside every sample/output/inducing loop.
• Mean-only prediction skips variance terms and avoids exp( q_log_s ), which accelerates classification / evaluation without changing predicted means.

For best CPU performance, build with the supplied Makefile and set OpenBLAS threading deliberately. If your SVGP worker count is already using all cores, try OPENBLAS_NUM_THREADS=1 during training to avoid oversubscription. If you mostly benchmark rebuild_kernel/Cholesky, allowing OpenBLAS to use all cores can be faster.

# 1. Create the data directory and enter it
mkdir -p data
cd data

# 2. Define the reliable Google Cloud storage base URL
BASE_URL="https://storage.googleapis.com/cvdf-datasets/mnist"

# 3. Download all 4 parts of the dataset using curl
curl -O "$BASE_URL/train-images-idx3-ubyte.gz"
curl -O "$BASE_URL/train-labels-idx1-ubyte.gz"
curl -O "$BASE_URL/t10k-images-idx3-ubyte.gz"
curl -O "$BASE_URL/t10k-labels-idx1-ubyte.gz"

# 4. Extract the files (Most native binary readers expect raw uncompressed bytes)
gunzip *.gz

./mnist_svgp ../data \
  --inducing 512 \
  --epochs 20 \
  --batch 128 \
  --lr 0.001 \
  --workers 16 \
  --lengthscale 9.8 \
  --noise 0.03 \
  --train-limit 60000

The most important accuracy/speed tradeoff is --inducing. More inducing points usually improves accuracy but costs roughly O( batch * inducing ^ 2 ) per training batch due to the Cholesky-whitened feature solve.

import svgp "path/to/svgp"

D := 784
O := 10
M := 256
cfg := svgp.default_config(D, O, M)
cfg.worker_count   = 0 // use all cores
cfg.lengthscale    = 9.8
cfg.noise_variance = 0.03

model := svgp.init( cfg )
defer svgp.deinit( & model )

svgp.init_inducing_from_data( & model, train_x, train_count )

tc := svgp.default_train_config( )
tc.epochs        = 10
tc.batch_size    = 256
tc.learning_rate = 0.002

svgp.fit( & model, train_x, train_y, train_count, tc )

means := make( [ ]f64, test_count * O )
vars  := make( [ ]f64, test_count * O )
svgp.predict( & model, test_x, test_count, means, vars )

The model uses inducing variables u with a whitened representation v = L^-1 u, where Kzz = L L^T, and trains a diagonal Gaussian variational posterior:

q( v ) = N( mu, diag( exp( log_s ) ) )

For an input x, it computes:

phi( x ) = L^-1 k( Z, x )
mean[ f( x ) ] = phi( x )^T mu
var[ f( x ) ]  = k( x,x ) - || phi( x ) ||^2 + phi( x )^T diag( exp( log_s ) ) phi( x )

The optimizer minimizes the negative stochastic ELBO for a Gaussian likelihood:

N / B * sum_batch E_q[ -log p( y | f ) ] + KL( q( v ) || N( 0, I ) )

For raw MNIST pixels, an RBF GP is not the best possible model; a CNN or feature extractor followed by SVGP would be much stronger. For generic smooth regression with moderate input dimension, this implementation is much closer to the intended use case and should be much more competitive.

MIT Open Source

Best regards,
Joao Carvalho