A Factor Graph Approach to Scalable Multi-Output Gaussian Process Regression

This paper presents a factor graph approach to scalable multi-output Gaussian process (MOGP) regression, addressing the computational challenges of traditional kernel-matrix methods.

The paper addresses multi-output Gaussian process regression where the goal is to infer a vector-valued function f: X → R D with D correlated outputs given N training inputs and a noisy observation matrix Y with Gaussian noise. A binary mask O records which entries are observed, with arbitrary missingness allowed. The authors target exact Gaussian inference under the chain-induced factor-graph model with two requirements: (1) inference cost should grow linearly in the number of training points N for fixed output count D and fixed latent rank, and (2) the missingness mask O should be absorbed into the model itself so that arbitrary missingness patterns produce no additional restructuring cost.

The proposed method, called SS-LMC (State-Space Linear Model of Coregionalization), combines three key components:

1. Nearest-neighbor chain ordering: For M-dimensional inputs, a greedy nearest-neighbor heuristic constructs a one-dimensional chain ordering of a fixed candidate set of size C. The chain distances serve as time steps in a state-space model, creating a chain-induced kernel approximation where the kernel is evaluated at chain distances rather than Euclidean distances.

2. State-space Gaussian processes: Matérn-class covariance functions are expressed as linear stochastic differential equations, converted to discrete-time state-space models. The full state z i ∈ R(2L) concatenates all L latent states with block-diagonal system matrices.

3. Linear Model of Coregionalization (LMC): The LMC mixes L latent GPs through a mixing matrix W ∈ R(D×L), expressed as a deterministic mixing factor with per-output scalar observation factors.

The generative model is expressed as a Forney-style factor graph where each output at each chain position is an independent scalar observation factor. Crucially, entries with O i,d = 0 contribute no factor so that predictions flow through the graph unaffected.

The paper provides several theoretical results:

• Lemma 1 bounds the difference between the Euclidean kernel matrix and the chain-induced kernel matrix: ∥K l − K π,l∥ 2 ≤ C α l ε π, where ε π is the maximum chain distortion and α l is the Lipschitz constant of the Matérn kernel.

• Theorem 1 extends this to the multi-output LMC case, bounding the difference in posterior parameters (mean and covariance) between the exact and chain-induced formulations. The bound scales as η ≤ C ε π Σ l ∥w l∥ 22 α l.

• Proposition 1 bounds the variation across different chain starting points: T(π1) ≤ (½⌈log2C⌉ + 1)T(π2), showing the chain length varies by at most O(log C) across starting points.

The per-step inference cost is O(C(DL2 + L3)) after chain construction. The paper notes that "the block-diagonal structure of A i and Q i keeps the latent dynamics sparse, but the LMC observation matrix H couples the L latent processes, so exact inference maintains dense Gaussian beliefs over the joint 2L-dimensional state." The D per-output scalar observation factors contribute O(DL2) work per chain position, and the dense Gaussian state updates contribute the O(L3) term.

For partial observations, the effective D at each chain position may be smaller, reducing the cost further. The factor graph handles arbitrary masks natively: "Each potential observation y i,d corresponds to an independent scalar factor, and the local update at chain position i only fires the factors for which O i,d = 1. No covariance-matrix restructuring is needed and no code path changes."

The paper sweeps input dimension M ∈ 2, 4, 8, 16, 32 on a synthetic sensor network benchmark with C=2000 candidates, D=3 outputs, L=2 latents, and 50/50 train/test split. Results show:

• At low input dimension (M=2), SS-LMC closely tracks the exact kernel-matrix posterior: held-out RMSE is 0.061 versus 0.058 at M=2 and 0.055 versus 0.046 at M=4

• As M increases, the chain stretches (mean ∆ rises from 0.10 to 5.6) and all methods degrade gracefully

• At M=32, SS-LMC reaches RMSE 0.531 vs. 0.417 for exact KM-LMC, with the gap staying ≤ 0.12 throughout

• Cost separates methods sharply: SS-LMC takes 0.015–0.017 s (flat in M), vs. 1.27–1.36 s for KM-LMC, 0.07–0.17 s for SVGP-LMC, and 0.41–0.51 s for NNGP-LMC

The paper forecasts the ETTh1 dataset with M=3 inputs and D=4 correlated outputs, sweeping window length N ∈ 500, 1000, 2000, 4000, 8000 at dropout p=0.3 and dropout rate p ∈ 0.1, 0.3, 0.5, 0.7 at N=2000. Results show:

• At N=2000: held-out RMSE is 8.17 for SS-LMC against 7.98 for exact KM-LMC, 8.10 for SVGP-LMC and 7.98 for NNGP-LMC

• At N=8000 (beyond KM-LMC's reach): SS-LMC reaches RMSE 5.15 vs. 4.95 for SVGP-LMC and 5.19 for NNGP-LMC

• SS-LMC's cost grows linearly, reaching 0.16 s at N=8000, 7.6× below the inducing-point SVGP-LMC (1.25 s) and 38× below NNGP-LMC (6.16 s)

• SS-LMC's cost is essentially invariant to the missingness rate (∼0.036 s across p) while KM-LMC's cost falls with sparser data but stays tied to dense factorization

The paper notes that chain construction is "preprocessing that the kernel-matrix and inducing-point baselines avoid, but it depends only on the candidate inputs, and not on Y, the mask O, or the hyperparameters. So it is built once and reused by every later inference call."

Hyperparameter learning is discussed: "The mixing matrix enters only through H and so touches the C observation blocks alone, whereas l l and γ l2 enter through A i and Q i and re-parameterize every transition factor: unlike the missingness mask, a hyperparameter update is not a local graph edit."

Limitations acknowledged include: (1) the greedy nearest-neighbor chain is path-dependent (though empirically negligible, with held-out RMSE varying by <0.1% across 20 starting points), (2) systematic evaluation at higher input dimensions on real data is left to future work, (3) the theoretical bounds are illustrative from a theoretical perspective but do not furnish error estimates one would use in practice, and (4) the method is currently limited to Matérn-class covariances that admit state-space representations.

The paper concludes that low-input-dimensional multi-output GP inference as message passing on a Forney-style factor graph, combining state-space GPs, LMC, and reactive message passing tracks the exact kernel-matrix posterior at low input dimension and degrades gracefully as the chain stretches with input dimension. On real ETTh1 forecasting, SS-LMC matches both baselines in forecast accuracy to within a small margin, while its inference cost grows only linearly in the window length and is invariant to the dropout rate. The factor graph handles partial observations without covariance-matrix restructuring.

Improvements to AI Systems Based on This Paper:

1. Scalable Multi-Output Regression with Arbitrary Missing Data

• Implement SS-LMC as a drop-in module for AI systems that must predict multiple correlated targets (e.g., sensor networks, financial indicators, climate variables) from high-dimensional inputs.

• The system can now handle incomplete observation matrices natively—no data imputation or covariance restructuring needed—while maintaining linear inference cost in the number of training points.

2. Real-Time Forecasting with Long Time Series

• Use the state-space factor graph to build forecasting systems (e.g., energy demand, traffic, stock prices) that process windows of 8,000+ time steps in under 0.2 seconds, unlike kernel-matrix methods that become computationally intractable.

• The system can update predictions incrementally as new data arrives, with cost invariant to the fraction of missing values (e.g., sensor dropouts).

3. Uncertainty-Aware Decision Making

• Leverage the exact Gaussian posterior (mean and covariance) from the factor graph to provide calibrated uncertainty estimates for each output, enabling risk-sensitive AI applications (e.g., autonomous navigation, medical diagnostics) where knowing confidence is critical.

• The chain-induced approximation preserves posterior quality at low input dimensions (RMSE within 5% of exact methods), so the system can trust its uncertainty bounds.

4. Modular Hyperparameter Learning

• The factor graph’s separation of mixing matrix (observation blocks) from kernel hyperparameters (transition blocks) allows efficient gradient-based optimization: only local blocks need re-computation per hyperparameter update, reducing retraining cost for AI systems that must adapt to new data distributions.

5. Handling High-Dimensional Inputs Gracefully

• The nearest-neighbor chain ordering provides a principled way to compress high-dimensional input spaces into a one-dimensional sequence, with theoretical bounds on posterior error.

• An improved system can automatically choose the chain candidate set size C to trade off accuracy vs. speed, degrading gracefully as input dimension increases (RMSE gap ≤ 0.12 even at M=32).

6. Preprocessing Reuse Across Tasks

• Since the chain construction depends only on input locations (not on outputs, missingness, or hyperparameters), an AI system can build the chain once and reuse it for multiple tasks (e.g., different target variables, different missingness patterns, or repeated model retraining) with zero additional cost.

7. Memory-Efficient Batch Processing

• The block-diagonal state-space structure reduces memory footprint from O(N2) kernel matrices to O(N·L2) state representations, enabling deployment on edge devices or in-memory processing of datasets that would otherwise exceed RAM limits.

What the Improved AI System Can Do:

• Process multi-output regression tasks with millions of training points and arbitrary missing data in linear time, on commodity hardware.

• Provide real-time, uncertainty-calibrated forecasts for long-horizon time series (e.g., 8,000-step windows) with sub-second latency, even under 70% data dropout.

• Automatically adapt to new input dimensions or missingness patterns without code changes or covariance matrix recomputation.

• Support both batch and streaming inference, with the same factor graph code path, for applications ranging from sensor fusion to financial risk modeling.