Sparsity-Aware Distributed Learning for Gaussian Processes with Linear Multiple Kernel

GP regression models an unknown function with a covariance kernel ¹. The grid spectral mixture (GSM) kernel ² approximates any stationary kernel's spectral density with a Gaussian mixture whose means and variances are fixed to preselected grid points, leaving only the mixture weights to be optimized. This turns the GSM kernel into a linear multiple kernel (LMK): a weighted sum of one-dimensional sub-kernels. The GSM kernel was designed for one-dimensional inputs. Extending it naively to P-dimensional data by fixing grid points in ℝ^P makes the number of components grow exponentially with P^a, which the authors' own earlier work on a multi-dimensional GSM kernel ³ already ran into.

This paper's grid spectral mixture product (GSMP) kernel avoids that blow-up by taking the product of one-dimensional sub-kernels across dimensions, rather than fixing grid points directly in the joint input space, while still using only Q weights. The resulting spectral density is itself a Gaussian mixture, so GSMP retains the universal approximation property of the original GSM and spectral mixture (SM) kernels⁴. In a two-dimensional worked example, GSMP's spectral density covers the same modes as the earlier grid-based GSM kernel plus additional ones lying along the main diagonal that GSM structurally cannot reach, and this gap widens with dimension P, so GSMP reaches the same expressiveness with far fewer hyperparameters.

Fitting the GSM/GSMP weights means minimizing the negative log marginal likelihood, a difference-of-convex (DC) program solved by successive convex approximation (SCA)⁶: at each step, a convex surrogate from a first-order Taylor expansion of the concave part is minimized with an off-the-shelf convex solver⁷. This paper proves every local minimum of this problem for the GSMP kernel is sparse, with or without observation noise, so only a handful of sub-kernels end up with non-negligible weight. Centralized SCA still scales as O(Qn³) and needs the whole dataset in one place at every iteration, which is impractical for large or privacy-sensitive data.

SLIM-KL splits the dataset across N agents, each holding a local shard, and asks them to agree on a shared set of GSMP weights without pooling their data. It builds on ADMM⁸, but quantizes every transmitted value with a stochastic quantizer⁹ that rounds each weight up or down with probability proportional to its distance from the nearest grid point, keeping the quantized value unbiased in expectation. Each round has agents average their quantized local weights into a global estimate, broadcast it, then re-optimize their own local weights against it before updating their dual variable. Because quantization rounds many already-small GSMP weights to exactly zero, it reinforces rather than fights the kernel's sparsity, which the paper argues can also help steer agents away from poor local minima^b.

The local subproblem each agent solves at every ADMM round is itself non-convex and non-separable across weights, so solving it directly does not scale to large Q. The distributed SCA (DSCA) algorithm splits an agent's weights into s blocks, builds one convex surrogate per block from the same Taylor-expansion trick as centralized SCA, and optimizes all s blocks in parallel on separate computing units. Each block's surrogate is reformulated as a second-order cone program so it can be solved efficiently with MOSEK ⁷ rather than a generic convex solver. Splitting the work this way reduces the per-iteration computational complexity of hyperparameter optimization from O(Qn³) to O(Qn³/(sN³)), so adding agents and computing units directly cuts the cost of training a large number of hyperparameters.

The local subproblem each agent solves at every round is itself non-convex and non-separable across weights, so the distributed SCA (DSCA) algorithm splits an agent's weights into s blocks, builds one convex surrogate per block, and optimizes all s blocks in parallel on separate computing units, each reformulated as a second-order cone program solved with MOSEK⁷. This cuts the per-iteration cost of hyperparameter optimization from O(Qn³) to O(Qn³/(sN³)), so adding agents and computing units directly cuts the cost of training a large number of hyperparameters.

All experiments use MATLAB with MOSEK ⁷ as the convex solver. Code and data are available in the paper's repository.

On synthetic two-dimensional data generated from an SM-kernel GP with a known spectral density, the GSMP kernel is compared against the earlier multi-dimensional GSM kernel ³ across five spectral-density configurations, using Q = 50 grid points for both. When the ground-truth density's modes lie on the main diagonal, the grid-based GSM kernel cannot represent them and its GP prediction MSE is far worse than the GSMP kernel's, for example 0.9431 versus 0.3769 in one configuration. When the modes instead lie on the anti-diagonal, a direction the GSM kernel's fixed grid already supports, the two kernels perform comparably, with GSM sometimes slightly ahead. The GSMP kernel matches or beats the GSM kernel's prediction accuracy in most configurations without requiring any additional grid points.

The GSMP-kernel GP trained with SLIM-KL (s = 4 computing units, N = 2 agents, quantization resolution Δ = 0.01) is compared on nine real datasets, spanning input dimension 1 to 11, against the original GSM-kernel GP ³, an SM-kernel GP ⁴, a sparse spectrum GP (SSGP) ⁵, a squared-exponential GP ¹, an LSTM ¹⁰, and the Informer Transformer forecaster ¹¹. GSMPGP achieves the lowest prediction MSE on every one of the nine datasets, and its solutions stay sparse (typically well under 5% of weights nonzero) at accuracy comparable to or better than the also-sparse GSMGP baseline, while the dense SMGP baseline needs 100% of its weights and still trails both. The improvement is largest on the multi-dimensional datasets, consistent with the approximation-capability results above.

Varying the number of DSCA computing units s ∈ {1, 4, 10, 50, 100} shows prediction MSE does not consistently worsen as s grows, and for several datasets more computing units even improves it, while computation time drops sharply as s increases, matching the O(Qn³/(sN³)) complexity bound. Varying the number of agents N ∈ {2, 4, 8, 10} against a centralized baseline, SLIM-KL with a moderate number of agents (N = 2) matches or beats the centralized case on nearly every dataset, but performance degrades once N grows large relative to a dataset's size, since each agent then has too little local data to train on. On a version of the CCCP dataset enlarged to 9,500 training points, SLIM-KL (N = 10, s = 4) keeps the same prediction MSE as the smaller version, showing it scales to substantially larger data.

Stochastic quantization is compared against deterministic quantization at Δ = 0.01, repeating each run 5 times. Stochastic quantization achieves a consistently lower mean prediction MSE than deterministic quantization, with its error bars (±2 standard deviations) showing the variability stays bounded.

Sweeping the quantization resolution Δ ∈ {0.001, 0.01, 0.1, 1, 5} against no quantization at all shows quantization does not necessarily hurt accuracy: on the CO2 dataset, Δ = 1 gives the lowest MSE of any setting while cutting transmitted bits by 4.8× on average, and on the Water dataset quantization saves communication by more than 200× while keeping prediction MSE essentially unchanged. Coarser resolutions eventually do degrade accuracy on some datasets (for example ECG at Δ = 5), so the paper reports Δ = 0.1 as a practical balance between communication savings and prediction performance.

On five larger-scale forecasting benchmarks (Weather and four ETT datasets), SLIM-KL is scaled up to s = 100 computing units and N = 100 agents with Δ = 0.01, and compared against recent Transformer- and RNN-based forecasters: SegRNN ¹², PatchTST ¹³, MICN ¹⁴, and TiDE ¹⁵. GSMPGP finishes in the top two on four of the five datasets and first on two of them (ETTh1 and ETTm1), while remaining competitive on Weather where it trails the best model only slightly. Against distributed and random-feature multiple-kernel-learning baselines, DOMKL ¹⁶ and CoKle/BoKle ¹⁷, SLIM-KL gets the lowest MSE on all five datasets, showing that scaling to 100 agents does not come at the cost of accuracy.

Hyperparameter optimization for this kernel yields sparse weights, which quantized ADMM then exploits, shrinking communication and, in principle, the model class the predictor uses. Deriving whether those zeros improve prediction, or mainly cut communication while leaving the same function class, remains open.