Online Gradient Computation for Warping Gaussian Process Transformations

# Online Gradient Computation for Warping Gaussian Process Transformations
Source: [https://arxiv.org/html/2609.16472](https://arxiv.org/html/2609.16472)
Emilio Ruiz\-MorenoAffiliation:SIGIPRO Department, Simula Metropolitan Center for Digital Engineering, Oslo, Norway;Affiliation:SURE\-AI Center, Simula Research Laboratory, Oslo, Norway;Konstantinos SlavakisAffiliation:Institute of Science Tokyo, Yokohama, Japan\. This work was supported by the DISCO grant 338740, the DRIVE grant 360486, and the SURE\-AI Center grant 357482 from the Research Council of Norway\. code:[https://github\.com/SIGIPRO/onlineWGP](https://github.com/SIGIPRO/onlineWGP)Baltasar Beferull\-LozanoAffiliation:SIGIPRO Department, Simula Metropolitan Center for Digital Engineering, Oslo, Norway;Affiliation:SURE\-AI Center, Simula Research Laboratory, Oslo, Norway;

###### Abstract

Warped Gaussian processes \(GPs\) handle non\-Gaussian observations by mapping them into a latent standard GP via a parametric transformation called warping\. Existing streaming variants, however, either optimize the warping parameters periodically or sacrifice analytical tractability for a higher model capacity\. To bridge this gap, we show that the gradient of the instantaneous negative log\-likelihood of a warped GP admits an exact recursive computation\. Based on this result, we propose a novel online method for warped GPs that jointly updates the latent GP moments and optimizes the warping parameters\.

###### Index Terms:

Gaussian process, warping transformation, online learning

## 1Introduction

Gaussian processes \(GPs\) provide a non\-parametric Bayesian framework for probabilistic function approximation\[[16](https://arxiv.org/html/2609.16472#bib.bib1)\], proving especially useful when the true underlying function is analytically unknown or expensive to query\[[22](https://arxiv.org/html/2609.16472#bib.bib8),[19](https://arxiv.org/html/2609.16472#bib.bib9)\]\. By placing a prior directly over the space of functions of interest, GPs offer a principled approach to uncertainty quantification, making them particularly attractive for applications where predictive confidence is as critical as point accuracy; for instance, in robotics\[[5](https://arxiv.org/html/2609.16472#bib.bib7)\]or geostatistics\[[7](https://arxiv.org/html/2609.16472#bib.bib5),[6](https://arxiv.org/html/2609.16472#bib.bib6)\]\. Crucially, the mathematical tractability of GPs allows for exact Bayesian inference, yielding closed\-form predictive distributions when paired with Gaussian function observation likelihoods\.

Despite their widespread success, standard GPs are fundamentally limited by their core assumption that the function observations are adequately modeled by a joint Gaussian distribution\. Indeed, this assumption often breaks down in practice, as real\-world observations may exhibit heavy tails, distinct skewness, or physical constraints \(such as strict positivity\)\. Consequently, applying standard GPs to such non\-Gaussian observations can produce highly inaccurate predictions and poorly calibrated uncertainty bounds\.

Motivated by these limitations, warped GPs\[[21](https://arxiv.org/html/2609.16472#bib.bib3)\]extend standard GPs by applying a parametric transformation—termed warping—to non\-Gaussian observations, mapping them into latent targets where standard GP assumptions hold\. By jointly estimating latent GP moments and warping parameters, warped GPs accommodate complex observation likelihoods while retaining the analytical tractability of standard GPs\.

However, the implementation of warped GPs in memory\-constrained or real\-time tasks remains a significant challenge\. This is because, to the best of our knowledge, there is no online method to jointly update the latent GP moments and optimize the warping parameters of warped GPs\. The methods arguably closest to this work in the literature either optimize the warping parameters periodically\[[12](https://arxiv.org/html/2609.16472#bib.bib23)\]or trade analytical tractability for higher model capacity\[[2](https://arxiv.org/html/2609.16472#bib.bib24)\]\.

In this paper, we show that the gradient of the negative log\-likelihood \(NLL\) of a warped GP can be computed recursively in a exact manner\. Based on this, we propose a novel online method for warped GPs that jointly updates the latent GP moments and optimizes the warping parameters\.

## 2Background

Figure 1:GP prior \(top\) and GP posterior \(bottom\), i\.e\., the GP prior conditioned on observations \(indicated by \+ markers\)\. The solid lines represent sample functions drawn from the GP\. The shaded area covers±1\.96\\pm 1\.96standard deviations from the mean\. Example adapted from\[[16](https://arxiv.org/html/2609.16472#bib.bib1), Ch\. 2\]\.Consider the problem of inferring an unknown, real\-valued functionffover an arbitrarydd\-dimensional input space𝒳⊆ℝd\\mathcal\{X\}\\subseteq\\mathbb\{R\}^\{d\}, based on prior beliefs about its function space and a set ofnn\(possibly corrupted\) function observationsy1,…,yny\_\{1\},\\dots,y\_\{n\}at corresponding input locations𝒙1,…,𝒙n\\bm\{x\}\_\{1\},\\dots,\\bm\{x\}\_\{n\}\.

Such a problem can be addressed analytically, provided certain assumptions hold, as we detail next\.

### 2\.1Gaussian processes

A GP model typically assumes the functionffis a realization of a zero\-mean GP prior and the function observations are corrupted by white Gaussian noise\. That is,

f⁡\(𝒙\)\\displaystyle f\(\\bm\{x\}\)∼𝒢​𝒫​\(0,κ⁡\(𝒙,𝒙′\)\),\\displaystyle\\sim\\mathcal\{GP\}\\left\(0,\\kappa\(\\bm\{x\},\\bm\{x\}^\{\\prime\}\)\\right\),\(1a\)yi\\displaystyle y\_\{i\}=f⁡\(𝒙i\)\+ϵi,\\displaystyle=f\(\\bm\{x\}\_\{i\}\)\+\\epsilon\_\{i\},\(1b\)where the kernelκ:𝒳×𝒳→ℝ\\kappa:\\mathcal\{X\}\\times\\mathcal\{X\}\\to\\mathbb\{R\}is the covariance function of the GP prior, and eachϵi∼𝒩⁡\(0,σ2\)\\epsilon\_\{i\}\\sim\\mathcal\{N\}\(0,\\sigma^\{2\}\)is an observation noise term with standard deviationσ\\sigma\.

Thanks to Gaussianity and the closure properties of GPs\[[16](https://arxiv.org/html/2609.16472#bib.bib1)\], the posterior distribution over functions conditioned on a set ofnnobservations collected in the observation vector𝒚n=\[y1,…,yn\]⊤∈ℝn\\bm\{y\}\_\{n\}=\[y\_\{1\},\\dots,y\_\{n\}\]^\{\\top\}\\in\\mathbb\{R\}^\{n\}is also Gaussian\. Specifically,

f⁡\(𝒙\)\|𝒚n∼𝒩⁡\(mn​\(𝒙\),vn​\(𝒙\)\),f\(\\bm\{x\}\)\|\\bm\{y\}\_\{n\}\\sim\\mathcal\{N\}\\left\(m\_\{n\}\(\\bm\{x\}\),v\_\{n\}\(\\bm\{x\}\)\\right\),\(2\)with closed\-form posterior mean and variance

mn​\(𝒙\)\\displaystyle m\_\{n\}\(\\bm\{x\}\)=𝒌n​\(𝒙\)⊤​𝜶n,\\displaystyle=\\bm\{k\}\_\{n\}\(\\bm\{x\}\)^\{\\top\}\\bm\{\\alpha\}\_\{n\},\(3a\)vn​\(𝒙\)\\displaystyle v\_\{n\}\(\\bm\{x\}\)=κ⁡\(𝒙,𝒙\)−𝒌n​\(𝒙\)⊤​𝛀n​𝒌n​\(𝒙\),\\displaystyle=\\kappa\(\\bm\{x\},\\bm\{x\}\)\-\\bm\{k\}\_\{n\}\(\\bm\{x\}\)^\{\\top\}\\bm\{\\Omega\}\_\{n\}\\bm\{k\}\_\{n\}\(\\bm\{x\}\),\(3b\)respectively, where

𝜶n\\displaystyle\\bm\{\\alpha\}\_\{n\}=𝛀n​𝒚n∈ℝn,\\displaystyle=\\bm\{\\Omega\}\_\{n\}\\bm\{y\}\_\{n\}\\in\\mathbb\{R\}^\{n\},\(4a\)𝛀n\\displaystyle\\bm\{\\Omega\}\_\{n\}=\(𝑲n\+σ2​𝑰n\)−1∈𝕊≥0n,,\\displaystyle=\\left\(\\bm\{K\}\_\{n\}\+\\sigma^\{2\}\\bm\{I\}\_\{n\}\\right\)^\{\-1\}\\in\\mathbb\{S\}^\{n\}\_\{\\geq 0\},,\(4b\)with𝒌n​\(𝒙\)∈ℝn\\bm\{k\}\_\{n\}\(\\bm\{x\}\)\\in\\mathbb\{R\}^\{n\}and𝑲n∈𝕊≥0n\\bm\{K\}\_\{n\}\\in\\mathbb\{S\}^\{n\}\_\{\\geq 0\}constructed as\[𝒌n​\(𝒙\)\]i=κ⁡\(𝒙i,𝒙\)\[\\bm\{k\}\_\{n\}\(\\bm\{x\}\)\]\_\{i\}=\\kappa\(\\bm\{x\}\_\{i\},\\bm\{x\}\)and\[𝑲n\]i,j=κ⁡\(𝒙i,𝒙j\)\[\\bm\{K\}\_\{n\}\]\_\{i,j\}=\\kappa\(\\bm\{x\}\_\{i\},\\bm\{x\}\_\{j\}\), for alli,j∈\{1,…,n\}i,j\\in\\\{1,\\dots,n\\\}\.

To aid conceptual understanding, an illustration of the prior \([1a](https://arxiv.org/html/2609.16472#S2.E1.1)\) and posterior \([2](https://arxiv.org/html/2609.16472#S2.E2)\) of a GP is provided in Fig\.[1](https://arxiv.org/html/2609.16472#S2.F1)\.

#### 2\.1\.1Recursive Gaussian process updates

Direct computation of the terms𝜶n\\bm\{\\alpha\}\_\{n\}and𝛀n\\bm\{\\Omega\}\_\{n\}for the posterior moments in \([3](https://arxiv.org/html/2609.16472#S2.E3)\) scales poorly with the number of observationsnn, mostly due to the matrix inversion on the right\-hand side of \([4b](https://arxiv.org/html/2609.16472#S2.E4.2)\)\. Fortunately, one can still compute these terms exactly, but at a reduced computational load, by exploiting the following recursive structure

𝜶n\\displaystyle\\bm\{\\alpha\}\_\{n\}=\[𝜶n−10\]−yn−mn−1​\(𝒙n\)sn​\[𝝎n−1\],\\displaystyle=\\begin\{bmatrix\}\\bm\{\\alpha\}\_\{n\-1\}\\\\ 0\\end\{bmatrix\}\-\\frac\{y\_\{n\}\-m\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\}\{s\_\{n\}\}\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}\\\\ \-1\\end\{bmatrix\},\(5a\)𝛀n\\displaystyle\\bm\{\\Omega\}\_\{n\}=\[𝛀n−1𝟎n−1𝟎n−1⊤0\]\+1sn​\[𝝎n−1\]​\[𝝎n⊤−1\],\\displaystyle=\\begin\{bmatrix\}\\bm\{\\Omega\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{s\_\{n\}\}\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\},\(5b\)wheresn=vn−1​\(𝒙n\)\+σ2∈ℝs\_\{n\}=v\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\+\\sigma^\{2\}\\in\\mathbb\{R\}, and𝝎n=𝛀n−1​𝒌n−1​\(𝒙n\)∈ℝn−1\\bm\{\\omega\}\_\{n\}=\\bm\{\\Omega\}\_\{n\-1\}\\bm\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\in\\mathbb\{R\}^\{n\-1\}\. Thus, the GP posterior \([2](https://arxiv.org/html/2609.16472#S2.E2)\) can be updated recursively as each new observationyny\_\{n\}and associated input location𝒙n\\bm\{x\}\_\{n\}become available, simply by maintaining𝜶n\\bm\{\\alpha\}\_\{n\}and𝛀n\\bm\{\\Omega\}\_\{n\}as state variables\.

#### 2\.1\.2Sparse Gaussian process updates

Indeed, the recursive updates in Sec\.[2\.1\.1](https://arxiv.org/html/2609.16472#S2.SS1.SSS1)effectively avoid the computational complexity bottleneck of ann×nn\\times nmatrix inversion\. However, computing the posterior moments in \([3](https://arxiv.org/html/2609.16472#S2.E3)\) still requires evaluating ann×1n\\times 1vector dot product for the mean and ann×nn\\times nquadratic form for the variance\. That is, the computational and memory requirements scale as𝒪⁡\(n\)\\mathcal\{O\}\(n\)for the mean, and𝒪⁡\(n2\)\\mathcal\{O\}\(n^\{2\}\)for the variance\.

This unbounded growth with the number of observationsnn, is known in the literature as the “curse of kernelization”\[[23](https://arxiv.org/html/2609.16472#bib.bib14),[3](https://arxiv.org/html/2609.16472#bib.bib15)\]and is typically addressed via kernel approximation methods such as Nyström\[[24](https://arxiv.org/html/2609.16472#bib.bib17)\]and random Fourier features\[[15](https://arxiv.org/html/2609.16472#bib.bib16),[13](https://arxiv.org/html/2609.16472#bib.bib18)\], or via dictionary\-based sparsification methods including forgetting mechanisms\[[20](https://arxiv.org/html/2609.16472#bib.bib19)\]and the approximate linear dependency \(ALD\) technique\[[8](https://arxiv.org/html/2609.16472#bib.bib10),[9](https://arxiv.org/html/2609.16472#bib.bib11)\]\.

In this work, we rely on the ALD technique, as it naturally integrates with the recursive updates in Sec\.[2\.1\.1](https://arxiv.org/html/2609.16472#S2.SS1.SSS1)\(see Sec\.[4](https://arxiv.org/html/2609.16472#S4)for more details\)\.

### 2\.2Warped Gaussian processes

Warped GP models assume that observations follow a GP model only after applying a monotonic, parametric, and differentiable transformation referred to as warping\[[21](https://arxiv.org/html/2609.16472#bib.bib3)\]\. Specifically, given a warping transformationggofrrparameters𝜽∈Θ⊆ℝr\\bm\{\\theta\}\\in\\Theta\\subseteq\\mathbb\{R\}^\{r\}, the observation model in \([1b](https://arxiv.org/html/2609.16472#S2.E1.2)\) becomes

zi:=g⁡\(yi,𝜽\)=f⁡\(𝒙i\)\+ϵi,z\_\{i\}:=g\(y\_\{i\};\\bm\{\\theta\}\)=f\(\\bm\{x\}\_\{i\}\)\+\\epsilon\_\{i\},\(6\)where eachzi∈ℝz\_\{i\}\\in\\mathbb\{R\}serves as a latent target\. Accordingly, the term𝜶n\\bm\{\\alpha\}\_\{n\}, defined in \([4a](https://arxiv.org/html/2609.16472#S2.E4.1)\), becomes𝜶n=𝛀n​𝒛n∈ℝn\\bm\{\\alpha\}\_\{n\}=\\bm\{\\Omega\}\_\{n\}\\bm\{z\}\_\{n\}\\in\\mathbb\{R\}^\{n\}, with𝒛n=\[z1,…,zn\]⊤∈ℝn\\bm\{z\}\_\{n\}=\[z\_\{1\},\\dots,z\_\{n\}\]^\{\\top\}\\in\\mathbb\{R\}^\{n\}\.

Warped GPs can be seen as a generalization of GPs\. In fact, they are typically non\-Gaussian and even asymmetric in the observation space\.

By invoking the change of variables formula\[[4](https://arxiv.org/html/2609.16472#bib.bib2), Ch\. 2\], the joint probability densityppof anynnobservations𝒚n\\bm\{y\}\_\{n\}can be expressed in terms of the corresponding densityqqof latent targets𝒛n\\bm\{z\}\_\{n\}, and the warping transformationggand warping parameters𝜽\\bm\{\\theta\}\. Specifically,

p⁡\(𝒚n\|𝜽\)=q⁡\(𝒈⁡\(𝒚n,𝜽\)\)​∏i=1n∂g⁡\(yi,𝜽\)∂yi,p\(\\bm\{y\}\_\{n\}\|\\bm\{\\theta\}\)=q\(\\bm\{g\}\(\\bm\{y\}\_\{n\};\\bm\{\\theta\}\)\)\\prod\_\{i=1\}^\{n\}\\frac\{\\partial g\(y\_\{i\};\\bm\{\\theta\}\)\}\{\\partial y\_\{i\}\},\(7\)where𝒈\\bm\{g\}applies the warping transformationggelementwise, and the joint probability of latent targets follows

q⁡\(𝒛\)=\(\(2​π\)n​\|𝛀n−1\|\)−12​exp⁡\(−12​𝒛n⊤​𝛀n​𝒛n\)q\(\\bm\{z\}\)=\\left\(\(2\\pi\)^\{n\}\\left\|\\bm\{\\Omega\}\_\{n\}^\{\-1\}\\right\|\\right\)^\{\-\\frac\{1\}\{2\}\}\\exp\\left\(\-\\frac\{1\}\{2\}\\bm\{z\}\_\{n\}^\{\\top\}\\bm\{\\Omega\}\_\{n\}\\bm\{z\}\_\{n\}\\right\)\(8\)by construction from \([6](https://arxiv.org/html/2609.16472#S2.E6)\)\.

One of the greatest strengths of warped GPs is that the warping parameters𝜽\\bm\{\\theta\}can be estimated directly from the observations𝒚n\\bm\{y\}\_\{n\}without requiring explicit knowledge of their joint probability densitypp\. As long as the latent densityqqand the warping transformationggare defined, \([7](https://arxiv.org/html/2609.16472#S2.E7)\) provides a tractable likelihood that can be optimized\.

Algorithm 1Proposed online method for warped GPs\.1:% tilde superscripts distinguish sparsified variables from their exact counterparts\.

2:Choose

gg,

Θ\\Theta,

κ\\kappa,

\{σn\}\\\{\\sigma\_\{n\}\\\}, the ALD sparsification threshold

ν\\nu, and the projected optimizer

ΠΘ\\Pi\_\{\\Theta\}\.

3:Initialize

𝜽0∈Θ\\bm\{\\theta\}\_\{0\}\\in\\Theta,

𝑲~0−1=1/κ⁡\(𝒙1,𝒙1\)\\tilde\{\\bm\{K\}\}^\{\-1\}\_\{0\}=1/\\kappa\(\\bm\{x\}\_\{1\},\\bm\{x\}\_\{1\}\),

𝜶~0=𝑪~0=0\\tilde\{\\bm\{\\alpha\}\}\_\{0\}=\\tilde\{\\bm\{C\}\}\_\{0\}=0,

𝑩~0=𝟎r\\tilde\{\\bm\{B\}\}\_\{0\}=\\bm\{0\}\_\{r\}, and the dictionary

𝒟0=\{𝒙1\}\\mathcal\{D\}\_\{0\}=\\\{\\bm\{x\}\_\{1\}\\\}\.

4:for

n=1,2,…n=1,2,\\dotsdo

5:Observe

𝒙n\\bm\{x\}\_\{n\}and

yny\_\{n\}
6:Compute

𝒂^n=𝑲~n−1−1​𝒌~n−1​\(𝒙t\)\\hat\{\\bm\{a\}\}\_\{n\}=\\tilde\{\\bm\{K\}\}^\{\-1\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{t\}\), and

7:

δn=κ⁡\(𝒙n,𝒙n\)−𝒌~n−1​\(𝒙n\)⊤​𝒂^n\\delta\_\{n\}=\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\hat\{\\bm\{a\}\}\_\{n\}
8:Warp

zn=g⁡\(yn,𝜽n−1\)z\_\{n\}=g\(y\_\{n\};\\bm\{\\theta\}\_\{n\-1\}\)and

𝒛˙n=∇𝜽g​\(yn,𝜽n−1\)\\dot\{\\bm\{z\}\}\_\{n\}=\\nabla\_\{\\bm\{\\theta\}\}\\,g\(y\_\{n\};\\bm\{\\theta\}\_\{n\-1\}\)
9:Compute

e~n=zn−𝒌~n−1​\(𝒙n\)⊤​𝜶~n−1\\tilde\{e\}\_\{n\}=z\_\{n\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\}, and

10:

𝒃~n=𝒛˙n−𝑩~n−1​𝒌~n−1​\(𝒙n\)\\tilde\{\\bm\{b\}\}\_\{n\}=\\dot\{\\bm\{z\}\}\_\{n\}\-\\tilde\{\\bm\{B\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)
11:if

δn\>ν\\delta\_\{n\}\>\\nuthen

12:Update

𝑲~n−1=1δn​\[δn​𝑲~n−1−1\+𝒂^n​𝒂^n⊤−𝒂^n−𝒂^n⊤1\]\\tilde\{\\bm\{K\}\}^\{\-1\}\_\{n\}=\\frac\{1\}\{\\delta\_\{n\}\}\\begin\{bmatrix\}\\delta\_\{n\}\\tilde\{\\bm\{K\}\}^\{\-1\}\_\{n\-1\}\+\\hat\{\\bm\{a\}\}\_\{n\}\\hat\{\\bm\{a\}\}\_\{n\}^\{\\top\}&\-\\hat\{\\bm\{a\}\}\_\{n\}\\\\ \-\\hat\{\\bm\{a\}\}\_\{n\}^\{\\top\}&1\\end\{bmatrix\}
13:Compute

𝒄~n=\[−𝑪~n−1​𝒌~n−1​\(𝒙n\)1\]\\tilde\{\\bm\{c\}\}\_\{n\}=\\begin\{bmatrix\}\-\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ 1\\end\{bmatrix\}, and

14:

s~n=κ⁡\(𝒙n,𝒙n\)\+σn2−𝒌~n−1​\(𝒙n\)⊤​𝑪~n−1​𝒌~n−1​\(𝒙n\)\\tilde\{s\}\_\{n\}\\\!=\\\!\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\\\!\+\\\!\\sigma^\{2\}\_\{n\}\\\!\-\\\!\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)
15:Reshape

𝑪~n−1=\[𝑪~n−1𝟎Dn−1𝟎Dn−1⊤0\]\\tilde\{\\bm\{C\}\}\_\{n\-1\}=\\begin\{bmatrix\}\\tilde\{\\bm\{C\}\}\_\{n\-1\}&\\bm\{0\}\_\{D\_\{n\}\-1\}\\\\ \\bm\{0\}\_\{D\_\{n\}\-1\}^\{\\top\}&0\\end\{bmatrix\},

16:

𝑩~n−1=\[𝑩~n−1,𝟎r\]\\tilde\{\\bm\{B\}\}\_\{n\-1\}=\[\\tilde\{\\bm\{B\}\}\_\{n\-1\},\\bm\{0\}\_\{r\}\], and

𝜶~n−1=\[𝜶~n−10\]\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\}=\\begin\{bmatrix\}\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\}\\\\ 0\\end\{bmatrix\}
17:Update

𝒟n=𝒟n−1∪\{𝒙n\}\\mathcal\{D\}\_\{n\}=\\mathcal\{D\}\_\{n\-1\}\\cup\\\{\\bm\{x\}\_\{n\}\\\}
18:else

19:Update

𝑲~n=𝑲~n−1\\tilde\{\\bm\{K\}\}\_\{n\}=\\tilde\{\\bm\{K\}\}\_\{n\-1\}
20:Compute

𝒄~n=𝒂^n−𝑪~n−1​𝒌~n−1​\(𝒙n\)\\tilde\{\\bm\{c\}\}\_\{n\}=\\hat\{\\bm\{a\}\}\_\{n\}\-\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\), and

s~n=\\tilde\{s\}\_\{n\}=
21:

𝒌~n−1​\(𝒙n\)⊤​𝒂^n\+σn2−𝒌~n−1​\(𝒙n\)⊤​𝑪~n−1​𝒌~n−1​\(𝒙n\)\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\hat\{\\bm\{a\}\}\_\{n\}\\\!\+\\\!\\sigma\_\{n\}^\{2\}\\\!\-\\\!\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)
22:Update

𝒟n=𝒟n−1\\mathcal\{D\}\_\{n\}=\\mathcal\{D\}\_\{n\-1\}
23:endif

24:Update

𝜶~n=𝜶~n−1\+e~ns~n​𝒄~n\\tilde\{\\bm\{\\alpha\}\}\_\{n\}\\\!=\\\!\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\}\\\!\+\\\!\\frac\{\\tilde\{e\}\_\{n\}\}\{\\tilde\{s\}\_\{n\}\}\\tilde\{\\bm\{c\}\}\_\{n\},

𝑩~n=𝑩~n−1\+1s~n​𝒃~n​𝒄~n⊤\\tilde\{\\bm\{B\}\}\_\{n\}\\\!=\\\!\\tilde\{\\bm\{B\}\}\_\{n\-1\}\\\!\+\\\!\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\tilde\{\\bm\{b\}\}\_\{n\}\\tilde\{\\bm\{c\}\}\_\{n\}^\{\\top\},

25:and

𝑪~n=𝑪~n−1\+1s~n​𝒄~n​𝒄~n⊤\\tilde\{\\bm\{C\}\}\_\{n\}=\\tilde\{\\bm\{C\}\}\_\{n\-1\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\tilde\{\\bm\{c\}\}\_\{n\}\\tilde\{\\bm\{c\}\}\_\{n\}^\{\\top\}
26:Compute

∇𝜽ℓ~n​\(𝜽n−1\)=e~ns~n​𝒃~n−∇𝜽​log​\(∂zn∂yn\)\\nabla\_\{\\bm\{\\theta\}\}\\tilde\{\\ell\}\_\{n\}\(\\bm\{\\theta\}\_\{n\-1\}\)=\\frac\{\\tilde\{e\}\_\{n\}\}\{\\tilde\{s\}\_\{n\}\}\\tilde\{\\bm\{b\}\}\_\{n\}\-\\nabla\_\{\\bm\{\\theta\}\}\\log\(\\frac\{\\partial z\_\{n\}\}\{\\partial y\_\{n\}\}\)
27:Update

𝜽n=ΠΘ​\(𝜽n−1,∇𝜽ℓ~n​\(𝜽n−1\)\)\\bm\{\\theta\}\_\{n\}=\\Pi\_\{\\Theta\}\(\\bm\{\\theta\}\_\{n\-1\},\\nabla\_\{\\bm\{\\theta\}\}\\tilde\{\\ell\}\_\{n\}\(\\bm\{\\theta\}\_\{n\-1\}\)\)
28:endfor

29:Return

𝒟n\\mathcal\{D\}\_\{n\},

𝜶~n\\tilde\{\\bm\{\\alpha\}\}\_\{n\},

𝑪~n\\tilde\{\\bm\{C\}\}\_\{n\}, and

𝜽n\\bm\{\\theta\}\_\{n\}

![Refer to caption](https://arxiv.org/html/2609.16472v1/experiment.png)\(a\)True process, warped GP, and distribution slices \(gray dashed line\)
\(b\)Forward and inverse warping transformation
Figure 2:In Fig\.[2a](https://arxiv.org/html/2609.16472#S2.F2.sf1), the observations are marked with gray dots\. The dotted lines show the true generating distributionp⁡\(𝒚\)p\(\\bm\{y\}\)and the solid lines show the warped GP predictionp⁡\(𝒚\|𝜽\)p\(\\bm\{y\}\|\\bm\{\\theta\}\)\. The triplets of lines represent the median, along with the2\.52\.5th and97\.597\.5th percentiles in each case\. The cross\-section shows the probability densities atx=−π12x=\-\\frac\{\\pi\}\{12\}; that is,p⁡\(y\|x=−π12CLOSEp\(y\|x=\-\\frac\{\\pi\}\{12\}\) andp⁡\(y\|𝜽,x=−π12\)p\(y\|\\bm\{\\theta\},x=\-\\frac\{\\pi\}\{12\}\)\. Fig\.[2b](https://arxiv.org/html/2609.16472#S2.F2.sf2)shows the forward\-inverse pass of the estimated warping transformation as well as its inversion error\.

## 3Recursive gradient computation

Thenn\-observation joint NLL

Ln​\(𝜽\)=−log⁡p⁡\(𝒚n\|𝜽\)L\_\{n\}\(\\bm\{\\theta\}\)=\-\\log p\(\\bm\{y\}\_\{n\}\|\\bm\{\\theta\}\)\(9\)measures the discrepancy between the observations𝒚n\\bm\{y\}\_\{n\}and the joint probability densitypp\(implicitly\) described by the warping parameters𝜽\\bm\{\\theta\}\. Using the multiplication rule of probability\[[17](https://arxiv.org/html/2609.16472#bib.bib12)\], one can readily get that

Ln​\(𝜽\)=Ln−1​\(𝜽\)\+ℓn​\(𝜽\),L\_\{n\}\(\\bm\{\\theta\}\)=L\_\{n\-1\}\(\\bm\{\\theta\}\)\+\\ell\_\{n\}\(\\bm\{\\theta\}\),\(10\)where the instantaneous NLL

ℓn​\(𝜽\)=−log⁡p⁡\(yn\|𝒚n−1,𝜽\),\\ell\_\{n\}\(\\bm\{\\theta\}\)=\-\\log p\(y\_\{n\}\|\\bm\{y\}\_\{n\-1\},\\bm\{\\theta\}\),\(11\)quantifies how unexpected the last observationyny\_\{n\}was given the previous observations𝒚n−1\\bm\{y\}\_\{n\-1\}\.

In many applications, the set of observationsy1,…,yny\_\{1\},\\dots,y\_\{n\}is too large to fit in memory, or is streamed, rendering standard batch processing intractable\. As a result, estimating the warping parameters by directly minimizing the NLL in \([9](https://arxiv.org/html/2609.16472#S3.E9)\) with respect to𝜽\\bm\{\\theta\}may not be possible in practice\. On the other hand, by exploiting the additive decomposition of the NLL in \([10](https://arxiv.org/html/2609.16472#S3.E10)\), one can naturally transition from a batch to an online estimation framework\. Specifically, the warping parameter estimate can be iteratively refined utilizing the gradient of the instantaneous NLL, e\.g\., by a first\-order method\[[1](https://arxiv.org/html/2609.16472#bib.bib13)\]\.

Although evaluating every instantaneous NLLℓn​\(𝜽\)\\ell\_\{n\}\(\\bm\{\\theta\}\)depends on the previous observations𝒚n−1\\bm\{y\}\_\{n\-1\}, Theorem[3\.1](https://arxiv.org/html/2609.16472#S3.Thmtheorem1)demonstrates that its gradient can be computed recursively\.

###### Theorem 3\.1\(Recursive gradient computation\)\.

The gradient of the instantaneous NLL in \([11](https://arxiv.org/html/2609.16472#S3.E11)\) with respect to the warping parameters𝛉\\bm\{\\theta\}can be computed as

∇𝜽ℓn​\(𝜽\)=ensn​𝒃n−∇𝜽​log​\(∂zn∂yn\),\\nabla\_\{\\bm\{\\theta\}\}\\ell\_\{n\}\(\\bm\{\\theta\}\)=\\frac\{e\_\{n\}\}\{s\_\{n\}\}\\bm\{b\}\_\{n\}\-\\nabla\_\{\\bm\{\\theta\}\}\\log\\left\(\\frac\{\\partial z\_\{n\}\}\{\\partial y\_\{n\}\}\\right\),\(12\)where the termsen∈ℝe\_\{n\}\\in\\mathbb\{R\}and𝐛n∈ℝr\\bm\{b\}\_\{n\}\\in\\mathbb\{R\}^\{r\}correspond to

en\\displaystyle e\_\{n\}=𝒛n−mn−1​\(𝒙n\),and\\displaystyle=\\bm\{z\}\_\{n\}\-m\_\{n\-1\}\(\\bm\{x\}\_\{n\}\),\\text\{ and \}\(13a\)𝒃n\\displaystyle\\bm\{b\}\_\{n\}=𝒛˙n−𝑩n−1​𝒌n−1​\(𝒙n\),\\displaystyle=\\dot\{\\bm\{z\}\}\_\{n\}\-\\bm\{B\}\_\{n\-1\}\\bm\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\),\(13b\)with𝐳˙n=∇𝛉g​\(yn,𝛉\)∈ℝr\\dot\{\\bm\{z\}\}\_\{n\}=\\nabla\_\{\\bm\{\\theta\}\}g\(y\_\{n\};\\bm\{\\theta\}\)\\in\\mathbb\{R\}^\{r\}, and the term𝐁n∈ℝr×n\\bm\{B\}\_\{n\}\\in\\mathbb\{R\}^\{r\\times n\}acts as a state variable updated according to

𝑩n=\[𝑩n−1𝟎r\]−𝒛˙n−𝑩n−1​𝒌n−1​\(𝒙n\)sn​\[𝝎n⊤−1\]\.\\bm\{B\}\_\{n\}=\\begin\{bmatrix\}\\bm\{B\}\_\{n\-1\}&\\bm\{0\}\_\{r\}\\end\{bmatrix\}\-\\frac\{\\dot\{\\bm\{z\}\}\_\{n\}\-\\bm\{B\}\_\{n\-1\}\\bm\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\}\{s\_\{n\}\}\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\}\.\(14\)

###### Proof\.

See the supplementary material[S1](https://arxiv.org/html/2609.16472#S1a)\. ∎

As a result, the gradient of the instantaneous NLL∇𝜽ℓn​\(𝜽\)\\nabla\_\{\\bm\{\\theta\}\}\\ell\_\{n\}\(\\bm\{\\theta\}\)can be updated recursively as each new observationyny\_\{n\}and associated input location𝒙n\\bm\{x\}\_\{n\}become available by maintaining𝜶n\\bm\{\\alpha\}\_\{n\},𝛀n\\bm\{\\Omega\}\_\{n\}\(as in Sec\.[2\.1\.1](https://arxiv.org/html/2609.16472#S2.SS1.SSS1)\), and𝑩n\\bm\{B\}\_\{n\}as state variables\.

## 4Proposed online method

Algorithm[1](https://arxiv.org/html/2609.16472#alg1)outlines the proposed online method for warped GPs, which jointly updates the latent GP moments and optimizes the warping parameters\. It integrates the recursive GP updates from Sec\.[2\.1\.1](https://arxiv.org/html/2609.16472#S2.SS1.SSS1), the recursive gradient computation of the instantaneous NLL from Sec\.[3](https://arxiv.org/html/2609.16472#S3), and the ALD sparsification technique introduced in Sec\.[2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2)\. As a result, the update rule \([12](https://arxiv.org/html/2609.16472#S3.E12)\) yields an approximate gradient\.

The per\-step computational and memory complexity is𝒪⁡\(Dn2\+r​Dn\)\\mathcal\{O\}\(D\_\{n\}^\{2\}\+rD\_\{n\}\)\. We refer the reader to the supplementary material[S2](https://arxiv.org/html/2609.16472#S2a)and[S3](https://arxiv.org/html/2609.16472#S3a)for further details\.

## 5Revisiting the 1D regression task

To evaluate our proposed method, we adapt the simple 1D regression task introduced in\[[21](https://arxiv.org/html/2609.16472#bib.bib3), Sec\. 4\]to an online setting\. This task is specifically designed to generate non\-Gaussian observations with sharp transitions, causing standard GPs to fail; therefore, relying entirely on the modeling capabilities of the warped GP\.

Experimental setup\. The input domain is the real subset𝒳=\[−π,π\]⊂ℝ\\mathcal\{X\}=\[\-\\pi,\\pi\]\\subset\\mathbb\{R\}\. The observational data are generated from a sinusoidal signal corrupted by white Gaussian noise followed by a cubic root transformation\. That is, everynnth observation is generated asyn=\(sin⁡\(xn\)\+εn\)13y\_\{n\}=\(\\sin\(x\_\{n\}\)\+\\varepsilon\_\{n\}\)^\{\\frac\{1\}\{3\}\}withεn∼𝒩⁡\(0,132\)\\varepsilon\_\{n\}\\sim\\mathcal\{N\}\(0,\\frac\{1\}\{3^\{2\}\}\)\. Then, we generate10001000realizations of101101uniformly spaced observations each\. That is, a total of101000101000observations, input location pairs that are fed sequentially to the warped GP model\.

Model configuration\. The choice of warping transformation follows that of the revisited task\. That is,g⁡\(y,𝜽\)=y\+∑i=1tθ1,i​tanh⁡\(θ2,i​\(y\+θ3,i\)\)g\(y;\\bm\{\\theta\}\)=y\+\\sum^\{t\}\_\{i=1\}\\theta\_\{1,i\}\\tanh\(\\theta\_\{2,i\}\(y\+\\theta\_\{3,i\}\)\), where𝜽=\[𝜽1⊤,𝜽2⊤,𝜽3⊤\]⊤∈Θ\\bm\{\\theta\}=\[\\bm\{\\theta\}\_\{1\}^\{\\top\},\\bm\{\\theta\}\_\{2\}^\{\\top\},\\bm\{\\theta\}\_\{3\}^\{\\top\}\]^\{\\top\}\\in\\Theta, andΘ=\{𝜽∈ℝ3​t:𝜽1,𝜽2⪰𝟎t,𝟏t⊤𝜽1\>0,and𝟏t⊤𝜽2\>0\}\\Theta=\\\{\\bm\{\\theta\}\\in\\mathbb\{R\}^\{3t\}:\\bm\{\\theta\}\_\{1\},\\bm\{\\theta\}\_\{2\}\\succeq\\bm\{0\}\_\{t\},\\bm\{1\}\_\{t\}^\{\\top\}\\bm\{\\theta\}\_\{1\}\>0,\\text\{ and \}\\bm\{1\}\_\{t\}^\{\\top\}\\bm\{\\theta\}\_\{2\}\>0\\\}, to ensure monotonicity\. We uset=10t=10terms, which leads to a total ofr=30r=30warping parameters\. Similarly, the choice of the covariance function remains the same, i\.e\.,κ⁡\(x,x′\)=ka​exp⁡\(−12​kw2​\(x−x′\)2\)\\kappa\(x,x^\{\\prime\}\)=k\_\{a\}\\exp\(\-\\frac\{1\}\{2k\_\{w\}^\{2\}\}\(x\-x^\{\\prime\}\)^\{2\}\), withka=2k\_\{a\}=2, andkw=2⋅2​π101≈0\.124k\_\{w\}=2\\cdot\\frac\{2\\pi\}\{101\}\\approx 0\.124\. Lastly, the noise standard deviation of the latent observation model is set toσ=3\\sigma=3\.

Implementation details\. We use the Adam111Pytorch\[[14](https://arxiv.org/html/2609.16472#bib.bib21)\]default configuration\.optimizer\[[11](https://arxiv.org/html/2609.16472#bib.bib20)\], followed by a projection onto the feasible setΘ\\Theta\. We set the ALD threshold toν=0\.1\\nu=0\.1\.

Results\. The results of the experiment are summarized in Fig\.[2](https://arxiv.org/html/2609.16472#S2.F2)\. Fig\.[2a](https://arxiv.org/html/2609.16472#S2.F2.sf1)shows that the predictions of the warped GP closely align with the true generating distribution\. On the other hand, Fig\.[2b](https://arxiv.org/html/2609.16472#S2.F2.sf2)shows that the learned warped transformation closely resembles the true generating cubic transformations in the region of interest where most observationsyylie, i\.e\.,−1≲y≲1\-1\\lesssim y\\lesssim 1\. It also shows that its inversion is numerically stable\. The \(ALD sparsification\) dictionary ends up with6868atoms\. Finally, the warped GP achieves an empirical NLL of−0\.28\-0\.28nats per observation\. This is a0\.760\.76nat improvement over a standard GP under an identical model configuration except for a well\-calibrated noise variance\.

## 6Conclusion and future work

We presented a recursive gradient evaluation for the instantaneous NLL of warped GPs, requiring only one extra state variable over standard recursive implementations\. Building on this result, we introduced an online algorithm for jointly updating latent GP moments and optimizing warping parameters, demonstrating its performance on a warped GP benchmark\. Future work includes online learning of the kernel parameters \(or directly a suitable kernel\[[18](https://arxiv.org/html/2609.16472#bib.bib22)\]\) and the observation model noise variance; evaluations on broader benchmarks are underway\.

## References

- \[1\]\(2017\)First\-order methods in optimization\.SIAM\.Cited by:[§3](https://arxiv.org/html/2609.16472#S3.p2.1)\.
- \[2\]T\. Bui, D\. Hernández\-Lobato, J\. Hernandez\-Lobato, Y\. Li, and R\. Turner\(2016\)Deep gaussian processes for regression using approximate expectation propagation\.InInternational conference on machine learning,pp\.1472–1481\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p4.1)\.
- \[3\]D\. Calandriello, A\. Lazaric, and M\. Valko\(2017\)Efficient second\-order online kernel learning with adaptive embedding\.Advances in Neural Information Processing Systems30\.Cited by:[§2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2.p2.1)\.
- \[4\]G\. Casella and R\. Berger\(2024\)Statistical inference\.CRC press\.Cited by:[§2\.2](https://arxiv.org/html/2609.16472#S2.SS2.p3.1)\.
- \[5\]M\. Deisenroth and C\. E\. Rasmussen\(2011\)PILCO: a model\-based and data\-efficient approach to policy search\.InProceedings of the 28th International Conference on machine learning \(ICML\-11\),pp\.465–472\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p1.1)\.
- \[6\]M\. Dewey, L\. Wilcox, B\. Samset, and A\. Ekman\(2026\)Deep\-aerogp: deep kernel learning for projecting the regional climate response to anthropogenic aerosol emission changes\.Technical reportCopernicus Meetings\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p1.1)\.
- \[7\]P\. J\. Diggle, J\. A\. Tawn, and R\. A\. Moyeed\(1998\)Model\-based geostatistics\.Journal of the Royal Statistical Society Series C: Applied Statistics47\(3\),pp\.299–350\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p1.1)\.
- \[8\]Y\. Engel, S\. Mannor, and R\. Meir\(2004\)The kernel recursive least\-squares algorithm\.IEEE Transactions on signal processing52\(8\),pp\.2275–2285\.Cited by:[§2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2.p2.1)\.
- \[9\]Y\. Engel\(2005\)Algorithms and representations for reinforcement learning\.Hebrew University of Jerusalem Jerusalem\.Cited by:[§2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2.p2.1)\.
- \[10\]R\. A\. Horn and C\. R\. Johnson\(2012\)Matrix analysis\.Cambridge university press\.Cited by:[§S3\.1](https://arxiv.org/html/2609.16472#S3.SS1.p2.4),[§S3\.2](https://arxiv.org/html/2609.16472#S3.SS2.p2.2),[§S3\.2](https://arxiv.org/html/2609.16472#S3.SS2.p6.1)\.
- \[11\]D\. P\. Kingma and J\. Ba\(2014\)Adam: a method for stochastic optimization\.arXiv preprint arXiv:1412\.6980\.Cited by:[§5](https://arxiv.org/html/2609.16472#S5.p4.1)\.
- \[12\]P\. Kou, F\. Gao, and X\. Guan\(2013\)Sparse online warped gaussian process for wind power probabilistic forecasting\.Applied energy108,pp\.410–428\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p4.1)\.
- \[13\]R\. T\. Money, J\. P\. Krishnan, and B\. Beferull\-Lozano\(2023\)Sparse online learning with kernels using random features for estimating nonlinear dynamic graphs\.IEEE Transactions on Signal Processing71,pp\.2027–2042\.Cited by:[§2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2.p2.1)\.
- \[14\]A\. Paszke, S\. Gross, F\. Massa, A\. Lerer, J\. Bradbury, G\. Chanan, T\. Killeen, Z\. Lin, N\. Gimelshein, L\. Antiga,et al\.\(2019\)Pytorch: an imperative style, high\-performance deep learning library\.Advances in neural information processing systems32\.Cited by:[footnote 1](https://arxiv.org/html/2609.16472#footnote1)\.
- \[15\]A\. Rahimi and B\. Recht\(2007\)Random features for large\-scale kernel machines\.Advances in neural information processing systems20\.Cited by:[§2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2.p2.1)\.
- \[16\]C\. E\. Rasmussen\(2003\)Gaussian processes in machine learning\.InSummer school on machine learning,pp\.63–71\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p1.1),[Figure 1](https://arxiv.org/html/2609.16472#S2.F1),[§2\.1](https://arxiv.org/html/2609.16472#S2.SS1.p2.1)\.
- \[17\]S\. M\. Ross, S\. M\. Ross, S\. M\. Ross, S\. M\. Ross, and E\. Mathématicien\(2014\)A first course in probability\.Vol\.8,Pearson London\.Cited by:[§3](https://arxiv.org/html/2609.16472#S3.p1.2)\.
- \[18\]E\. Ruiz\-Moreno and B\. Beferull\-Lozano\(2023\)An online multiple kernel parallelizable learning scheme\.IEEE Signal Processing Letters31,pp\.121–125\.Cited by:[§6](https://arxiv.org/html/2609.16472#S6.p1.1)\.
- \[19\]E\. Ruiz\-Moreno and B\. Beferull\-Lozano\(2025\)Doubly truncated mode kriging\.In2025 IEEE Statistical Signal Processing Workshop \(SSP\),pp\.306–310\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p1.1)\.
- \[20\]K\. Slavakis, P\. Bouboulis, and S\. Theodoridis\(2014\)Online learning in reproducing kernel hilbert spaces\.InAcademic Press Library in Signal Processing,Vol\.1,pp\.883–987\.Cited by:[§2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2.p2.1)\.
- \[21\]E\. Snelson, Z\. Ghahramani, and C\. Rasmussen\(2003\)Warped gaussian processes\.Advances in neural information processing systems16\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p3.1),[§2\.2](https://arxiv.org/html/2609.16472#S2.SS2.p1.1),[§5](https://arxiv.org/html/2609.16472#S5.p1.1)\.
- \[22\]L\. P\. Swiler, M\. Gulian, A\. L\. Frankel, C\. Safta, and J\. D\. Jakeman\(2020\)A survey of constrained gaussian process regression: approaches and implementation challenges\.Journal of Machine Learning for Modeling and Computing1\(2\)\.Cited by:[§1](https://arxiv.org/html/2609.16472#S1.p1.1)\.
- \[23\]Z\. Wang, K\. Crammer, and S\. Vucetic\(2012\)Breaking the curse of kernelization: budgeted stochastic gradient descent for large\-scale svm training\.The Journal of Machine Learning Research13\(1\),pp\.3103–3131\.Cited by:[§2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2.p2.1)\.
- \[24\]C\. Williams and M\. Seeger\(2000\)Using the nyström method to speed up kernel machines\.Advances in neural information processing systems13\.Cited by:[§2\.1\.2](https://arxiv.org/html/2609.16472#S2.SS1.SSS2.p2.1)\.

Supplementary material for “Online Gradient Computation for Warping Gaussian Process Transformations”

## Proof of Theorem[3\.1](https://arxiv.org/html/2609.16472#S3.Thmtheorem1)

The NLL in \([9](https://arxiv.org/html/2609.16472#S3.E9)\) corresponds to

Ln​\(𝜽\)\\displaystyle L\_\{n\}\(\\bm\{\\theta\}\)=−log⁡p⁡\(𝒚n\|𝜽\)\\displaystyle=\-\\log p\(\\bm\{y\}\_\{n\}\|\\bm\{\\theta\}\)\(s1a\)=−logq\(𝒈\(𝒚n;𝜽\)\)−log∏i=1n∂g⁡\(yi,𝜽\)∂yi\\displaystyle=\-\\log q\(\\bm\{g\}\(\\bm\{y\}\_\{n\};\\bm\{\\theta\}\)\)\-\\log\\prod^\{n\}\_\{i=1\}\\frac\{\\partial g\(y\_\{i\};\\bm\{\\theta\}\)\}\{\\partial y\_\{i\}\}\(s1b\)=12​\(2​π\)n​\|𝛀n−1\|\+12​𝒈​\(𝒚n,𝜽\)⊤​𝛀n​𝒈​\(𝒚n,𝜽\)−∑i=1nlog⁡\(∂g⁡\(yi,𝜽\)∂yi\)\.\\displaystyle=\\frac\{1\}\{2\}\(2\\pi\)^\{n\}\\left\|\\bm\{\\Omega\}\_\{n\}^\{\-1\}\\right\|\+\\frac\{1\}\{2\}\\bm\{g\}\(\\bm\{y\}\_\{n\};\\bm\{\\theta\}\)^\{\\top\}\\,\\bm\{\\Omega\}\_\{n\}\\,\\bm\{g\}\(\\bm\{y\}\_\{n\};\\bm\{\\theta\}\)\-\\sum^\{n\}\_\{i=1\}\\log\\left\(\\frac\{\\partial g\(y\_\{i\};\\bm\{\\theta\}\)\}\{\\partial y\_\{i\}\}\\right\)\.\(s1c\)Its gradient with respect to the warping parameters𝜽\\bm\{\\theta\}is thus,

∇𝜽Ln​\(𝜽\)\\displaystyle\\nabla\_\{\\bm\{\\theta\}\}L\_\{n\}\(\\bm\{\\theta\}\)=∇𝜽12​\(2​π\)n​\|𝛀n−1\|\+12​∇𝜽⊤𝒈​\(𝒚n,𝜽\)​𝛀n​𝒈​\(𝒚n,𝜽\)−∑i=1n∇𝜽​log​\(∂g⁡\(yi,𝜽\)∂yi\)\\displaystyle=\\cancel\{\\nabla\_\{\\bm\{\\theta\}\}\\frac\{1\}\{2\}\(2\\pi\)^\{n\}\\left\|\\bm\{\\Omega\}\_\{n\}^\{\-1\}\\right\|\}\+\\frac\{1\}\{2\}\\nabla\_\{\\bm\{\\theta\}\}^\{\\top\}\\bm\{g\}\(\\bm\{y\}\_\{n\};\\bm\{\\theta\}\)\\,\\bm\{\\Omega\}\_\{n\}\\,\\bm\{g\}\(\\bm\{y\}\_\{n\};\\bm\{\\theta\}\)\-\\sum^\{n\}\_\{i=1\}\\nabla\_\{\\bm\{\\theta\}\}\\log\\left\(\\frac\{\\partial g\(y\_\{i\};\\bm\{\\theta\}\)\}\{\\partial y\_\{i\}\}\\right\)\(s2a\)=𝒁˙n​𝛀n​𝒛n−∑i=1n∇𝜽​log​\(∂zi∂yi\),\\displaystyle=\\dot\{\\bm\{Z\}\}\_\{n\}\\bm\{\\Omega\}\_\{n\}\\bm\{z\}\_\{n\}\-\\sum^\{n\}\_\{i=1\}\\nabla\_\{\\bm\{\\theta\}\}\\log\\left\(\\frac\{\\partial z\_\{i\}\}\{\\partial y\_\{i\}\}\\right\),\(s2b\)where

𝒛n\\displaystyle\\bm\{z\}\_\{n\}=𝒈⁡\(𝒚n,𝜽\)=\[z1zn\]∈ℝn,\\displaystyle=\\bm\{g\}\(\\bm\{y\}\_\{n\};\\bm\{\\theta\}\)=\\begin\{bmatrix\}z\_\{1\}\\\\ \\vdots\\\\ z\_\{n\}\\end\{bmatrix\}\\in\\mathbb\{R\}^\{n\},\(s3a\)𝒛˙i\\displaystyle\\dot\{\\bm\{z\}\}\_\{i\}=∇𝜽g\(yi;𝜽\)∈ℝr,for alli∈\{1,2,…,n\},and\\displaystyle=\\nabla\_\{\\bm\{\\theta\}\}g\(y\_\{i\};\\bm\{\\theta\}\)\\in\\mathbb\{R\}^\{r\},\\text\{ for all \}i\\in\\\{1,2,\\dots,n\\\},\\text\{ and \}\(s3b\)𝒁˙n\\displaystyle\\dot\{\\bm\{Z\}\}\_\{n\}=∇𝜽⊤​𝒈​\(𝒚n,𝜽\)=\[𝒛˙1𝒛˙2…𝒛˙n\]=\[𝒁˙n−1𝒛˙n\]∈ℝr×n\.\\displaystyle=\\nabla\_\{\\bm\{\\theta\}\}^\{\\top\}\\bm\{g\}\(\\bm\{y\}\_\{n\};\\bm\{\\theta\}\)=\\begin\{bmatrix\}\\dot\{\\bm\{z\}\}\_\{1\}&\\dot\{\\bm\{z\}\}\_\{2\}&\\dots&\\dot\{\\bm\{z\}\}\_\{n\}\\end\{bmatrix\}=\\begin\{bmatrix\}\\dot\{\\bm\{Z\}\}\_\{n\-1\}&\\dot\{\\bm\{z\}\}\_\{n\}\\end\{bmatrix\}\\in\\mathbb\{R\}^\{r\\times n\}\.\(s3c\)
The first term in \([s2b](https://arxiv.org/html/2609.16472#S1.E2.2)\) can be expanded by using the recursive representation of𝛀n\\bm\{\\Omega\}\_\{n\}described in \([5b](https://arxiv.org/html/2609.16472#S2.E5.2)\) as

𝒁˙n​𝛀n​𝒛n\\displaystyle\\dot\{\\bm\{Z\}\}\_\{n\}\\bm\{\\Omega\}\_\{n\}\\bm\{z\}\_\{n\}=\[𝒁˙n−1𝒛˙n\]⁡\(\[𝛀n−1𝟎n−1𝟎n−1⊤0\]\+1sn​\[𝝎n−1\]​\[𝝎n⊤−1\]\)​\[𝒛n−1zn\]\\displaystyle=\\begin\{bmatrix\}\\dot\{\\bm\{Z\}\}\_\{n\-1\}&\\dot\{\\bm\{z\}\}\_\{n\}\\end\{bmatrix\}\\left\(\\begin\{bmatrix\}\\bm\{\\Omega\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{s\_\{n\}\}\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\}\\right\)\\begin\{bmatrix\}\\bm\{z\}\_\{n\-1\}\\\\ z\_\{n\}\\end\{bmatrix\}\(s4a\)=𝒁˙n−1​𝛀n−1​𝒛n−1\+1sn​\(𝒁˙n−1​𝝎n−𝒛˙n\)​\(𝝎n⊤​𝒛n−1−zn\)\\displaystyle=\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\Omega\}\_\{n\-1\}\\bm\{z\}\_\{n\-1\}\+\\frac\{1\}\{s\_\{n\}\}\\left\(\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\omega\}\_\{n\}\-\\dot\{\\bm\{z\}\}\_\{n\}\\right\)\\left\(\\bm\{\\omega\}\_\{n\}^\{\\top\}\\bm\{z\}\_\{n\-1\}\-z\_\{n\}\\right\)\(s4b\)=𝒁˙n−1​𝛀n−1​𝒛n−1\+ensn​𝒃n,\\displaystyle=\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\Omega\}\_\{n\-1\}\\bm\{z\}\_\{n\-1\}\+\\frac\{e\_\{n\}\}\{s\_\{n\}\}\\bm\{b\}\_\{n\},\(s4c\)where

en\\displaystyle e\_\{n\}=zn−𝝎n⊤​𝒛n−1∈ℝ,and\\displaystyle=z\_\{n\}\-\\bm\{\\omega\}\_\{n\}^\{\\top\}\\bm\{z\}\_\{n\-1\}\\in\\mathbb\{R\},\\text\{ and \}\(s5a\)𝒃n\\displaystyle\\bm\{b\}\_\{n\}=𝒛˙n−𝒁˙n−1​𝝎n∈ℝr\.\\displaystyle=\\dot\{\\bm\{z\}\}\_\{n\}\-\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\omega\}\_\{n\}\\in\\mathbb\{R\}^\{r\}\.\(s5b\)Note thatene\_\{n\}already admits a recursive update rule since it is expressed in terms of the current latent targetznz\_\{n\}and its corresponding location𝒙n\\bm\{x\}\_\{n\}, and the previous state terms𝜶n−1\\bm\{\\alpha\}\_\{n\-1\}and𝛀n−1\\bm\{\\Omega\}\_\{n\-1\}\. Explicitly,

en\\displaystyle e\_\{n\}=zn−𝒌n−1​\(𝒙n\)⊤​𝛀n−1​𝒛n−1\\displaystyle=z\_\{n\}\-\\bm\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{\\Omega\}\_\{n\-1\}\\bm\{z\}\_\{n\-1\}\(s6a\)=zn−𝒌n−1​\(𝒙n\)⊤​𝜶n−1\\displaystyle=z\_\{n\}\-\\bm\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{\\alpha\}\_\{n\-1\}\(s6b\)=zn−mn−1​\(𝒙n\)\.\\displaystyle=z\_\{n\}\-m\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\.\(s6c\)Regarding𝒃n\\bm\{b\}\_\{n\}, we maintain an additional state term𝑩n=𝒁˙n​𝛀n∈ℝr×n\\bm\{B\}\_\{n\}=\\dot\{\\bm\{Z\}\}\_\{n\}\\bm\{\\Omega\}\_\{n\}\\in\\mathbb\{R\}^\{r\\times n\}to enable recursive updates\. Specifically,

𝑩n\\displaystyle\\bm\{B\}\_\{n\}=\[𝒁˙n−1𝒛˙n\]⁡\(\[𝛀n−1𝟎n−1𝟎n−1⊤0\]\+1sn​\[𝝎n−1\]​\[𝝎n⊤−1\]\)\\displaystyle=\\begin\{bmatrix\}\\dot\{\\bm\{Z\}\}\_\{n\-1\}&\\dot\{\\bm\{z\}\}\_\{n\}\\end\{bmatrix\}\\left\(\\begin\{bmatrix\}\\bm\{\\Omega\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{s\_\{n\}\}\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\}\\right\)\(s7a\)=\[𝒁˙n−1​𝛀n−1𝟎r\]\+1sn​\(𝒁˙n−1​𝝎n−𝒛˙n\)​\[𝝎n⊤−1\]\\displaystyle=\\begin\{bmatrix\}\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\Omega\}\_\{n\-1\}&\\bm\{0\}\_\{r\}\\end\{bmatrix\}\+\\frac\{1\}\{s\_\{n\}\}\\left\(\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\omega\}\_\{n\}\-\\dot\{\\bm\{z\}\}\_\{n\}\\right\)\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\}\(s7b\)=\[𝑩n−1𝟎r\]−1sn​\(𝒛˙n−𝑩n−1​𝒌n−1​\(𝒙n\)\)​\[𝝎n⊤−1\],\\displaystyle=\\begin\{bmatrix\}\\bm\{B\}\_\{n\-1\}&\\bm\{0\}\_\{r\}\\end\{bmatrix\}\-\\frac\{1\}\{s\_\{n\}\}\\left\(\\dot\{\\bm\{z\}\}\_\{n\}\-\\bm\{B\}\_\{n\-1\}\\bm\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\right\)\\begin\{bmatrix\}\\bm\{\\omega\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\},\(s7c\)and accordingly,

𝒃n\\displaystyle\\bm\{b\}\_\{n\}=𝒛˙n−𝒁˙n−1​𝛀n−1​𝒌n−1​\(𝒙n\)\\displaystyle=\\dot\{\\bm\{z\}\}\_\{n\}\-\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\Omega\}\_\{n\-1\}\\bm\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\(s8a\)=𝒛˙n−𝑩n−1​𝒌n−1​\(𝒙n\)\.\\displaystyle=\\dot\{\\bm\{z\}\}\_\{n\}\-\\bm\{B\}\_\{n\-1\}\\bm\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\.\(s8b\)
Substituting \([s4c](https://arxiv.org/html/2609.16472#S1.E4.3)\) in \([s2b](https://arxiv.org/html/2609.16472#S1.E2.2)\) we get

∇𝜽Ln​\(𝜽\)\\displaystyle\\nabla\_\{\\bm\{\\theta\}\}L\_\{n\}\(\\bm\{\\theta\}\)=𝒁˙n−1​𝛀n−1​𝒛n−1\+ensn​𝒃n−∑i=1n∇𝜽​log​\(∂zi∂yi\)\\displaystyle=\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\Omega\}\_\{n\-1\}\\bm\{z\}\_\{n\-1\}\+\\frac\{e\_\{n\}\}\{s\_\{n\}\}\\bm\{b\}\_\{n\}\-\\sum^\{n\}\_\{i=1\}\\nabla\_\{\\bm\{\\theta\}\}\\log\\left\(\\frac\{\\partial z\_\{i\}\}\{\\partial y\_\{i\}\}\\right\)\(s9a\)=𝒁˙n−1​𝛀n−1​𝒛n−1−∑i=1n−1∇𝜽​log​\(∂zi∂yi\)\+ensn​𝒃n−∇𝜽​log​\(∂zn∂yn\)\\displaystyle=\\dot\{\\bm\{Z\}\}\_\{n\-1\}\\bm\{\\Omega\}\_\{n\-1\}\\bm\{z\}\_\{n\-1\}\-\\sum^\{n\-1\}\_\{i=1\}\\nabla\_\{\\bm\{\\theta\}\}\\log\\left\(\\frac\{\\partial z\_\{i\}\}\{\\partial y\_\{i\}\}\\right\)\+\\frac\{e\_\{n\}\}\{s\_\{n\}\}\\bm\{b\}\_\{n\}\-\\nabla\_\{\\bm\{\\theta\}\}\\log\\left\(\\frac\{\\partial z\_\{n\}\}\{\\partial y\_\{n\}\}\\right\)\(s9b\)=∇𝜽Ln−1​\(𝜽\)\+ensn​𝒃n−∇𝜽​log​\(∂zn∂yn\)\.\\displaystyle=\\nabla\_\{\\bm\{\\theta\}\}L\_\{n\-1\}\(\\bm\{\\theta\}\)\+\\frac\{e\_\{n\}\}\{s\_\{n\}\}\\bm\{b\}\_\{n\}\-\\nabla\_\{\\bm\{\\theta\}\}\\log\\left\(\\frac\{\\partial z\_\{n\}\}\{\\partial y\_\{n\}\}\\right\)\.\(s9c\)Finally, by taking the gradient with respect to the warping parameters𝜽\\bm\{\\theta\}in both sides of \([10](https://arxiv.org/html/2609.16472#S3.E10)\) and identifying terms with \([s9](https://arxiv.org/html/2609.16472#S1.E9)\) we obtain

∇𝜽ℓn​\(𝜽\)=ensn​𝒃n−∇𝜽​log​\(∂zn∂yn\)\.\\nabla\_\{\\bm\{\\theta\}\}\\ell\_\{n\}\(\\bm\{\\theta\}\)=\\frac\{e\_\{n\}\}\{s\_\{n\}\}\\bm\{b\}\_\{n\}\-\\nabla\_\{\\bm\{\\theta\}\}\\log\\left\(\\frac\{\\partial z\_\{n\}\}\{\\partial y\_\{n\}\}\\right\)\.\(s10\)

## ALD sparsification

Consider a feature mappingϕ:𝒳→ℱ\\phi:\\mathcal\{X\}\\to\\mathcal\{F\}such thatκ⁡\(𝒙,𝒙′\)=⟨ϕ⁡\(𝒙\),ϕ⁡\(𝒙′\)⟩ℱ\\kappa\(\\bm\{x\},\\bm\{x\}^\{\\prime\}\)=\\langle\\phi\(\\bm\{x\}\),\\phi\(\\bm\{x\}^\{\\prime\}\)\\rangle\_\{\\mathcal\{F\}\}for all𝒙,𝒙′∈𝒳\\bm\{x\},\\bm\{x\}^\{\\prime\}\\in\\mathcal\{X\}, an ongoing trajectory of input locations𝒙1,𝒙2,…,𝒙n−1\\bm\{x\}\_\{1\},\\bm\{x\}\_\{2\},\\dots,\\bm\{x\}\_\{n\-1\}, and a dictionary \(sparse set\)𝒟n−1⊆𝒳\\mathcal\{D\}\_\{n\-1\}\\subseteq\\mathcal\{X\}of representative input locations collected across that trajectory\.

Every time a new input location𝒙n\\bm\{x\}\_\{n\}is presented, we check whether its feature representationϕ⁡\(𝒙n\)\\phi\(\\bm\{x\}\_\{n\}\)is approximately linearly dependent on the feature representation of the previously collected input locations𝒙~1,𝒙~2,…,𝒙~Dn−1\\tilde\{\\bm\{x\}\}\_\{1\},\\tilde\{\\bm\{x\}\}\_\{2\},\\dots,\\tilde\{\\bm\{x\}\}\_\{D\_\{n\-1\}\}in𝒟n−1\\mathcal\{D\}\_\{n\-1\}whereDn=\|𝒟n\|D\_\{n\}=\|\\mathcal\{D\}\_\{n\}\|\. That is, we check whether the squared distanceδn\\delta\_\{n\}betweenϕ⁡\(𝒙n\)\\phi\(\\bm\{x\}\_\{n\}\)andspan​\{ϕ⁡\(𝒙~1\),ϕ⁡\(𝒙~2\),…,ϕ⁡\(𝒙~Dn−1\)\}\\text\{span\}\\\{\\phi\(\\tilde\{\\bm\{x\}\}\_\{1\}\),\\phi\(\\tilde\{\\bm\{x\}\}\_\{2\}\),\\dots,\\phi\(\\tilde\{\\bm\{x\}\}\_\{D\_\{n\-1\}\}\)\\\}is less than or equal to a user\-defined thresholdν∈ℝ≥0\\nu\\in\\mathbb\{R\}\_\{\\geq 0\}\. If it is, the dictionary remains the same, i\.e\.,𝒟n=𝒟n−1\\mathcal\{D\}\_\{n\}=\\mathcal\{D\}\_\{n\-1\}\. If not, we update the dictionary to include the new input location, i\.e\.,𝒟n=𝒟n−1∪\{𝒙n\}\\mathcal\{D\}\_\{n\}=\\mathcal\{D\}\_\{n\-1\}\\cup\\\{\\bm\{x\}\_\{n\}\\\}\.

Explicitly,

δn\\displaystyle\\delta\_\{n\}=min𝒂∈ℝDn−1​‖ϕ⁡\(𝒙n\)−∑i=1Dn−1ai​ϕ​\(𝒙~i\)‖ℱ2\\displaystyle=\\underset\{\\bm\{a\}\\in\\mathbb\{R\}^\{D\_\{n\-1\}\}\}\{\\text\{ min \}\}\\left\\\|\\phi\(\\bm\{x\}\_\{n\}\)\-\\sum^\{D\_\{n\-1\}\}\_\{i=1\}a\_\{i\}\\phi\(\\tilde\{\\bm\{x\}\}\_\{i\}\)\\right\\\|^\{2\}\_\{\\mathcal\{F\}\}\(s11a\)=min𝒂∈ℝDn−1​⟨ϕ⁡\(𝒙n\)−∑i=1Dn−1ai​ϕ​\(𝒙~i\),ϕ⁡\(𝒙n\)−∑j=1Dn−1aj​ϕ​\(𝒙~j\)⟩ℱ\\displaystyle=\\underset\{\\bm\{a\}\\in\\mathbb\{R\}^\{D\_\{n\-1\}\}\}\{\\text\{ min \}\}\\left\\langle\\phi\(\\bm\{x\}\_\{n\}\)\-\\sum^\{D\_\{n\-1\}\}\_\{i=1\}a\_\{i\}\\phi\(\\tilde\{\\bm\{x\}\}\_\{i\}\),\\phi\(\\bm\{x\}\_\{n\}\)\-\\sum^\{D\_\{n\-1\}\}\_\{j=1\}a\_\{j\}\\phi\(\\tilde\{\\bm\{x\}\}\_\{j\}\)\\right\\rangle\_\{\\mathcal\{F\}\}\(s11b\)=min𝒂∈ℝDn−1​⟨ϕ⁡\(𝒙n\),ϕ⁡\(𝒙n\)⟩ℱ−2​∑i=1Dn−1ai​⟨ϕ⁡\(𝒙n\),ϕ⁡\(𝒙~i\)⟩ℱ\+∑i,j=1Dn−1ai​aj​⟨ϕ⁡\(𝒙~i\),ϕ⁡\(𝒙~j\)⟩ℱ\\displaystyle=\\underset\{\\bm\{a\}\\in\\mathbb\{R\}^\{D\_\{n\-1\}\}\}\{\\text\{ min \}\}\\left\\langle\\phi\(\\bm\{x\}\_\{n\}\),\\phi\(\\bm\{x\}\_\{n\}\)\\right\\rangle\_\{\\mathcal\{F\}\}\-2\\sum^\{D\_\{n\-1\}\}\_\{i=1\}a\_\{i\}\\left\\langle\\phi\(\\bm\{x\}\_\{n\}\),\\phi\(\\tilde\{\\bm\{x\}\}\_\{i\}\)\\right\\rangle\_\{\\mathcal\{F\}\}\+\\sum^\{D\_\{n\-1\}\}\_\{i,j=1\}a\_\{i\}a\_\{j\}\\left\\langle\\phi\(\\tilde\{\\bm\{x\}\}\_\{i\}\),\\phi\(\\tilde\{\\bm\{x\}\}\_\{j\}\)\\right\\rangle\_\{\\mathcal\{F\}\}\(s11c\)=min𝒂∈ℝDn−1​κ​\(𝒙n,𝒙n\)−2​∑i=1Dn−1ai​κ​\(𝒙n,𝒙~i\)\+∑i,j=1Dn−1ai​aj​κ​\(𝒙~i,𝒙~j\)\\displaystyle=\\underset\{\\bm\{a\}\\in\\mathbb\{R\}^\{D\_\{n\-1\}\}\}\{\\text\{ min \}\}\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\-2\\sum^\{D\_\{n\-1\}\}\_\{i=1\}a\_\{i\}\\kappa\(\\bm\{x\}\_\{n\},\\tilde\{\\bm\{x\}\}\_\{i\}\)\+\\sum^\{D\_\{n\-1\}\}\_\{i,j=1\}a\_\{i\}a\_\{j\}\\kappa\(\\tilde\{\\bm\{x\}\}\_\{i\},\\tilde\{\\bm\{x\}\}\_\{j\}\)\(s11d\)=min𝒂∈ℝDn−1​κ​\(𝒙n,𝒙n\)−2​𝒌~n−1​\(𝒙n\)⊤​𝒂\+𝒂⊤​𝑲~n−1​𝒂,\\displaystyle=\\underset\{\\bm\{a\}\\in\\mathbb\{R\}^\{D\_\{n\-1\}\}\}\{\\text\{ min \}\}\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\-2\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{a\}\+\\bm\{a\}^\{\\top\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\},\(s11e\)where𝑲~n∈ℝDn×Dn\\tilde\{\\bm\{K\}\}\_\{n\}\\in\\mathbb\{R\}^\{D\_\{n\}\\times D\_\{n\}\}and𝒌~n​\(𝒙\)∈ℝDn\\tilde\{\\bm\{k\}\}\_\{n\}\(\\bm\{x\}\)\\in\\mathbb\{R\}^\{D\_\{n\}\}are constructed as\[𝑲~n\]i,j=κ⁡\(𝒙~i,𝒙~j\)\[\\tilde\{\\bm\{K\}\}\_\{n\}\]\_\{i,j\}=\\kappa\(\\tilde\{\\bm\{x\}\}\_\{i\},\\tilde\{\\bm\{x\}\}\_\{j\}\)and\[𝒌~n​\(𝒙\)\]i=κ⁡\(𝒙~i,𝒙\)\[\\tilde\{\\bm\{k\}\}\_\{n\}\(\\bm\{x\}\)\]\_\{i\}=\\kappa\(\\tilde\{\\bm\{x\}\}\_\{i\},\\bm\{x\}\)for all𝒙~i,𝒙~j∈𝒟n\\tilde\{\\bm\{x\}\}\_\{i\},\\tilde\{\\bm\{x\}\}\_\{j\}\\in\\mathcal\{D\}\_\{n\}and𝒙∈𝒳\\bm\{x\}\\in\\mathcal\{X\}\. Thus, findingδn\\delta\_\{n\}consists of solving a convex quadratic optimization problem with a closed\-form solution

δn\\displaystyle\\delta\_\{n\}=κ⁡\(𝒙n,𝒙n\)−𝒌~n−1​\(𝒙n\)⊤​𝒂^n,where\\displaystyle=\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\hat\{\\bm\{a\}\}\_\{n\},\\text\{ where \}\(s12a\)𝒂^n\\displaystyle\\hat\{\\bm\{a\}\}\_\{n\}=𝑲~n−1−1​𝒌~n−1​\(𝒙𝒏\)\.\\displaystyle=\\tilde\{\\bm\{K\}\}^\{\-1\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\_\{n\}\}\)\.\(s12b\)
From here, if𝒙n\\bm\{x\}\_\{n\}is included in𝒟n\\mathcal\{D\}\_\{n\}\(meaning thatδn\>ν\\delta\_\{n\}\>\\nu\), we set the approximation coefficients as𝒂n=\[0,…,0,1\]⊤∈\{0\}Dn−1×\{1\}⊂ℝDn\\bm\{a\}\_\{n\}=\[0,\\dots,0,1\]^\{\\top\}\\in\\\{0\\\}^\{D\_\{n\}\-1\}\\times\\\{1\\\}\\subset\\mathbb\{R\}^\{D\_\{n\}\}sinceϕ⁡\(𝒙n\)\\phi\(\\bm\{x\}\_\{n\}\)can be exactly represented by itself\. If not \(meaning thatδn≤ν\\delta\_\{n\}\\leq\\nu\), we set𝒂n=𝒂^n\\bm\{a\}\_\{n\}=\\hat\{\\bm\{a\}\}\_\{n\}with𝒂^n\\hat\{\\bm\{a\}\}\_\{n\}as in \([s12b](https://arxiv.org/html/2609.16472#S2.E12.2)\)\.

Let us construct the matrix𝚽n=\[ϕ⁡\(𝒙1\),ϕ⁡\(𝒙2\),…,ϕ⁡\(𝒙n\)\]\\bm\{\\Phi\}\_\{n\}=\[\\phi\(\\bm\{x\}\_\{1\}\),\\phi\(\\bm\{x\}\_\{2\}\),\\dots,\\phi\(\\bm\{x\}\_\{n\}\)\]of sizedim​\(ℱ\)×n\\text\{dim\}\(\\mathcal\{F\}\)\\times n\. By the sparsification procedure, we know that everyiith featureϕ⁡\(𝒙i\)\\phi\(\\bm\{x\}\_\{i\}\)is approximately linearly dependent on the features of the input locations in𝒟i\\mathcal\{D\}\_\{i\}, i\.e\.,ϕ⁡\(𝒙i\)≃∑j=1Di\[𝒂i\]j​ϕ​\(𝒙~j\)\\phi\(\\bm\{x\}\_\{i\}\)\\simeq\\sum^\{D\_\{i\}\}\_\{j=1\}\[\\bm\{a\}\_\{i\}\]\_\{j\}\\phi\(\\tilde\{\\bm\{x\}\}\_\{j\}\)\. Based on this, we can approximate𝚽n≃𝚽~n​𝑨n⊤\\bm\{\\Phi\}\_\{n\}\\simeq\\tilde\{\\bm\{\\Phi\}\}\_\{n\}\\bm\{A\}\_\{n\}^\{\\top\}where𝚽~n=\[ϕ⁡\(𝒙~1\),ϕ⁡\(𝒙~2\),…,ϕ⁡\(𝒙~Dn\)\]\\tilde\{\\bm\{\\Phi\}\}\_\{n\}=\[\\phi\(\\tilde\{\\bm\{x\}\}\_\{1\}\),\\phi\(\\tilde\{\\bm\{x\}\}\_\{2\}\),\\dots,\\phi\(\\tilde\{\\bm\{x\}\}\_\{D\_\{n\}\}\)\]is a matrix of sizedim​\(ℱ\)×Dn\\text\{dim\}\(\\mathcal\{F\}\)\\times D\_\{n\}, and𝑨n\\bm\{A\}\_\{n\}is a matrix of sizen×Dnn\\times D\_\{n\}for which eachiith row is constructed by padding zeros to𝒂i\\bm\{a\}\_\{i\}until dimensionDnD\_\{n\}\. On the other hand, we know that𝑲n=𝚽n⊤​𝚽n\\bm\{K\}\_\{n\}=\\bm\{\\Phi\}\_\{n\}^\{\\top\}\\bm\{\\Phi\}\_\{n\}by construction\. As a result, we can approximate

𝑲n\\displaystyle\\bm\{K\}\_\{n\}=𝚽n⊤​𝚽n\\displaystyle=\\bm\{\\Phi\}\_\{n\}^\{\\top\}\\bm\{\\Phi\}\_\{n\}\(s13a\)≃𝑨n​𝚽~n⊤​𝚽~n​𝑨n⊤\\displaystyle\\simeq\\bm\{A\}\_\{n\}\\tilde\{\\bm\{\\Phi\}\}\_\{n\}^\{\\top\}\\tilde\{\\bm\{\\Phi\}\}\_\{n\}\\bm\{A\}\_\{n\}^\{\\top\}\(s13b\)=𝑨n​𝑲~n​𝑨n⊤\.\\displaystyle=\\bm\{A\}\_\{n\}\\tilde\{\\bm\{K\}\}\_\{n\}\\bm\{A\}\_\{n\}^\{\\top\}\.\(s13c\)Similarly,

𝒌n​\(𝒙\)\\displaystyle\\bm\{k\}\_\{n\}\(\\bm\{x\}\)=\[ϕ​\(𝒙1\)⊤​ϕ​\(𝒙\)ϕ​\(𝒙2\)⊤​ϕ​\(𝒙\)⋯ϕ​\(𝒙n\)⊤​ϕ​\(𝒙\)\]=𝚽n⊤​ϕ​\(𝒙\)\\displaystyle=\\begin\{bmatrix\}\\phi\(\\bm\{x\}\_\{1\}\)^\{\\top\}\\phi\(\\bm\{x\}\)&\\phi\(\\bm\{x\}\_\{2\}\)^\{\\top\}\\phi\(\\bm\{x\}\)&\\cdots&\\phi\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\phi\(\\bm\{x\}\)\\end\{bmatrix\}=\\bm\{\\Phi\}^\{\\top\}\_\{n\}\\phi\(\\bm\{x\}\)\(s14a\)≃\(𝚽~n​𝑨n⊤\)⊤​ϕ​\(𝒙\)=𝑨n​𝚽~n⊤​ϕ​\(𝒙\)\\displaystyle\\simeq\\left\(\\tilde\{\\bm\{\\Phi\}\}\_\{n\}\\bm\{A\}\_\{n\}^\{\\top\}\\right\)^\{\\top\}\\phi\(\\bm\{x\}\)=\\bm\{A\}\_\{n\}\\tilde\{\\bm\{\\Phi\}\}\_\{n\}^\{\\top\}\\phi\(\\bm\{x\}\)\(s14b\)=𝑨n​𝒌~n​\(𝒙\)\.\\displaystyle=\\bm\{A\}\_\{n\}\\tilde\{\\bm\{k\}\}\_\{n\}\(\\bm\{x\}\)\.\(s14c\)
By extension, plugging \([s13](https://arxiv.org/html/2609.16472#S2.E13)\) and \([s14](https://arxiv.org/html/2609.16472#S2.E14)\) into \([3](https://arxiv.org/html/2609.16472#S2.E3)\) yields the following approximation of the posterior moments

mn​\(𝒙\)\\displaystyle m\_\{n\}\(\\bm\{x\}\)≃m~n​\(𝒙\)=𝒌~n​\(𝒙\)⊤​𝜶~n,and\\displaystyle\\simeq\\tilde\{m\}\_\{n\}\(\{\\bm\{x\}\}\)=\\tilde\{\\bm\{k\}\}\_\{n\}\(\\bm\{x\}\)^\{\\top\}\\tilde\{\\bm\{\\alpha\}\}\_\{n\},\\text\{ and \}\(s15a\)vn​\(𝒙\)\\displaystyle v\_\{n\}\(\\bm\{x\}\)≃v~n​\(𝒙\)=κ⁡\(𝒙,𝒙\)−𝒌~n​\(𝒙\)⊤​𝑪~n​𝒌~n​\(𝒙\),\\displaystyle\\simeq\\tilde\{v\}\_\{n\}\(\\bm\{x\}\)=\\kappa\(\\bm\{x\},\\bm\{x\}\)\-\\tilde\{\\bm\{k\}\}\_\{n\}\(\\bm\{x\}\)^\{\\top\}\\tilde\{\\bm\{C\}\}\_\{n\}\\tilde\{\\bm\{k\}\}\_\{n\}\(\\bm\{x\}\),\(s15b\)where

𝜶~n\\displaystyle\\tilde\{\\bm\{\\alpha\}\}\_\{n\}=𝑨n⊤​𝛀~n​𝒚n∈ℝDn,\\displaystyle=\\bm\{A\}\_\{n\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\}\\bm\{y\}\_\{n\}\\in\\mathbb\{R\}^\{D\_\{n\}\},\(s16a\)𝛀~n\\displaystyle\\tilde\{\\bm\{\\Omega\}\}\_\{n\}=\(𝑨n​𝑲~n​𝑨n⊤\+𝚺n\)−1∈ℝn×n,and\\displaystyle=\\left\(\\bm\{A\}\_\{n\}\\tilde\{\\bm\{K\}\}\_\{n\}\\bm\{A\}\_\{n\}^\{\\top\}\+\\bm\{\\Sigma\}\_\{n\}\\right\)^\{\-1\}\\in\\mathbb\{R\}^\{n\\times n\},\\text\{ and \}\(s16b\)𝑪~n\\displaystyle\\tilde\{\\bm\{C\}\}\_\{n\}=𝑨n⊤​𝛀~n​𝑨n∈ℝDn×Dn\.\\displaystyle=\\bm\{A\}^\{\\top\}\_\{n\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\}\\bm\{A\}\_\{n\}\\in\\mathbb\{R\}^\{D\_\{n\}\\times D\_\{n\}\}\.\(s16c\)

## ALD sparse and recursive Gaussian process updates

The goal is updating the approximated posterior moments in \([s15](https://arxiv.org/html/2609.16472#S2.E15)\) from the current observationyny\_\{n\}, input location𝒙n\\bm\{x\}\_\{n\}, and the previous state variables𝜶~n−1\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\}and𝑪~n−1\\tilde\{\\bm\{C\}\}\_\{n\-1\}\. Recall from Sec\.[S2](https://arxiv.org/html/2609.16472#S2a)that everynnth input location𝒙n\\bm\{x\}\_\{n\}may either be left out of the dictionary, in which case𝒟n=𝒟n−1\\mathcal\{D\}\_\{n\}=\\mathcal\{D\}\_\{n\-1\}, or added to it, in which case𝒟n=𝒟n−1∪\{𝒙n\}\\mathcal\{D\}\_\{n\}=\\mathcal\{D\}\_\{n\-1\}\\cup\\\{\\bm\{x\}\_\{n\}\\\}\. Depending on the case, the recursive update will vary accordingly\.

### S3\.1Case𝒟n=𝒟n−1\\mathcal\{D\}\_\{n\}=\\mathcal\{D\}\_\{n\-1\}\(meaning thatδn≤ν\\delta\_\{n\}\\leq\\nu\)

Since the dictionary remains unchanged

𝑲~n\\displaystyle\\tilde\{\\bm\{K\}\}\_\{n\}=𝑲~n−1∈ℝDn×Dn,\\displaystyle=\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\in\\mathbb\{R\}^\{D\_\{n\}\\times D\_\{n\}\},\(s17a\)𝑨n\\displaystyle\\bm\{A\}\_\{n\}=\[𝑨n−1𝒂n⊤\]∈ℝn×Dn,and\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}\\\\ \\bm\{a\}\_\{n\}^\{\\top\}\\end\{bmatrix\}\\in\\mathbb\{R\}^\{n\\times D\_\{n\}\},\\text\{ and \}\(s17b\)𝒂n\\displaystyle\\bm\{a\}\_\{n\}=𝑲~n−1−1​𝒌~n−1​\(𝒙n\)\.\\displaystyle=\\tilde\{\\bm\{K\}\}^\{\-1\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\.\(s17c\)
Then,

𝛀~n−1\\displaystyle\\tilde\{\\bm\{\\Omega\}\}\_\{n\}^\{\-1\}=\[𝑨n−1𝒂n⊤\]​𝑲~n−1​\[𝑨n−1⊤𝒂n\]\+\[σ2​𝑰n−1𝟎n−1𝟎n−1σ2\]\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}\\\\ \\bm\{a\}\_\{n\}^\{\\top\}\\end\{bmatrix\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}&\\bm\{a\}\_\{n\}\\end\{bmatrix\}\+\\begin\{bmatrix\}\\sigma^\{2\}\\bm\{I\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}&\\sigma^\{2\}\\end\{bmatrix\}\(s18a\)=\[𝑨n−1​𝑲~n−1​𝑨n−1⊤𝑨n−1​𝑲~n−1​𝒂n𝒂n⊤​𝑲~n−1​𝑨n−1⊤𝒂n⊤​𝑲~n−1​𝒂n\]\+\[σ2​𝑰n−1𝟎n−1𝟎n−1σ2\]\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}^\{\\top\}&\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\}\_\{n\}\\\\ \\bm\{a\}\_\{n\}^\{\\top\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{A\}^\{\\top\}\_\{n\-1\}&\\bm\{a\}\_\{n\}^\{\\top\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\}\_\{n\}\\end\{bmatrix\}\+\\begin\{bmatrix\}\\sigma^\{2\}\\bm\{I\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}&\\sigma^\{2\}\\end\{bmatrix\}\(s18b\)=\[𝛀~n−1𝑨n−1​𝑲~n−1​𝒂n\(𝑨n−1​𝑲~n−1​𝒂n\)⊤𝒂n⊤​𝑲~n−1​𝒂n\+σ2\]\.\\displaystyle=\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\}^\{\-1\}&\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\}\_\{n\}\\\\ \\left\(\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\}\_\{n\}\\right\)^\{\\top\}&\\bm\{a\}\_\{n\}^\{\\top\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\}\_\{n\}\+\\sigma^\{2\}\\end\{bmatrix\}\.\(s18c\)Notice that,

𝑨n−1​𝑲~n−1​𝒂n\\displaystyle\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\}\_\{n\}=𝑨n−1​𝑲~n−1​𝑲~n−1−1​𝒌~n−1​\(𝒙n\)\\displaystyle=\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}^\{\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\(s19a\)=𝑨n−1​𝒌~n−1​\(𝒙n\),and\\displaystyle=\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\),\\text\{ and \}\(s19b\)𝒂n⊤​𝑲~n−1​𝒂n\\displaystyle\\bm\{a\}^\{\\top\}\_\{n\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\}\_\{n\}=𝒌~n−1​\(𝒙n\)⊤​𝑲~n−1−1​𝑲~n−1​𝒂n\\displaystyle=\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}^\{\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{a\}\_\{n\}\(s19c\)=𝒌~n−1​\(𝒙n\)⊤​𝒂n\.\\displaystyle=\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{a\}\_\{n\}\.\(s19d\)Therefore,

𝛀~n−1=\[𝛀~n−1−1𝑨n−1​𝒌~n−1​\(𝒙n\)\(𝑨n−1​𝒌~n−1​\(𝒙n\)\)⊤𝒌~n−1​\(𝒙n\)⊤​𝒂n\+σ2\],\\tilde\{\\bm\{\\Omega\}\}\_\{n\}^\{\-1\}=\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}^\{\-1\}&\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ \\left\(\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\right\)^\{\\top\}&\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{a\}\_\{n\}\+\\sigma^\{2\}\\end\{bmatrix\},\(s20\)and applying the partitioned \(symmetric and positive semi\-definite\) matrix inverse formula\[[10](https://arxiv.org/html/2609.16472#bib.bib4), Ch\. 0\.7\]we get

𝛀~n=\[𝛀~n−1𝟎n−1𝟎n−1⊤0\]\+1s~n​\[𝝎~n−1\]​\[𝝎~n⊤−1\],\\tilde\{\\bm\{\\Omega\}\}\_\{n\}=\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\},\(s21\)where

𝝎~n\\displaystyle\\tilde\{\\bm\{\\omega\}\}\_\{n\}=𝛀~n−1​𝑨n−1​𝒌~n−1​\(𝒙n\)∈ℝn−1,and\\displaystyle=\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\in\\mathbb\{R\}^\{n\-1\},\\text\{ and \}\(s22a\)s~n\\displaystyle\\tilde\{s\}\_\{n\}=k~n−1​\(𝒙n\)⊤​𝒂n\+σ2−𝝎~n⊤​𝑨n−1​𝒌~n−1​\(𝒙n\)∈ℝ\.\\displaystyle=\\tilde\{k\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{a\}\_\{n\}\+\\sigma^\{2\}\-\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\in\\mathbb\{R\}\.\(s22b\)
Now,

𝑪~n\\displaystyle\\tilde\{\\bm\{C\}\}\_\{n\}=\[𝑨n−1⊤𝒂n\]⁡\(\[𝛀~n−1𝟎n−1𝟎n−1⊤0\]\+1s~n​\[𝝎~n−1\]​\[𝝎~n⊤−1\]\)​\[𝑨n−1𝒂n⊤\]\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}&\\bm\{a\}\_\{n\}\\end\{bmatrix\}\\left\(\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\}\\right\)\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}\\\\ \\bm\{a\}\_\{n\}^\{\\top\}\\end\{bmatrix\}\(s23a\)=𝑨n−1⊤​𝛀~n−1​𝑨n−1\+1s~n​\(𝑨n−1⊤​𝝎~n−𝒂n\)​\(𝝎~n⊤​𝑨n−1−𝒂n⊤\)\\displaystyle=\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\left\(\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\-\\bm\{a\}\_\{n\}\\right\)\\left\(\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}\\bm\{A\}\_\{n\-1\}\-\\bm\{a\}\_\{n\}^\{\\top\}\\right\)\(s23b\)=𝑪~n−1\+1s~n​𝒄~n​𝒄~n⊤,\\displaystyle=\\tilde\{\\bm\{C\}\}\_\{n\-1\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\tilde\{\\bm\{c\}\}\_\{n\}\\tilde\{\\bm\{c\}\}\_\{n\}^\{\\top\},\(s23c\)where

𝒄~n=𝒂n−𝑨n−1⊤​𝒘~n∈ℝDn\.\\tilde\{\\bm\{c\}\}\_\{n\}=\\bm\{a\}\_\{n\}\-\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{w\}\}\_\{n\}\\in\\mathbb\{R\}^\{D\_\{n\}\}\.\(s24\)
In the same way,

𝜶~n\\displaystyle\\tilde\{\\bm\{\\alpha\}\}\_\{n\}=\[𝑨n−1⊤𝒂n\]⁡\(\[𝛀~n−1𝟎n−1𝟎n−1⊤0\]\+1s~n​\[𝝎~n−1\]​\[𝝎~n⊤−1\]\)​\[𝒚n−1yn\]\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}&\\bm\{a\}\_\{n\}\\end\{bmatrix\}\\left\(\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\}\\right\)\\begin\{bmatrix\}\\bm\{y\}\_\{n\-1\}\\\\ y\_\{n\}\\end\{bmatrix\}\(s25a\)=𝑨n−1⊤​𝛀~n−1​𝒚n−1\+1s~n​\(𝑨n−1⊤​𝝎~n−𝒂n\)​\(𝝎~n⊤​𝒚n−1−yn\)\\displaystyle=\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{y\}\_\{n\-1\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\left\(\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\-\\bm\{a\}\_\{n\}\\right\)\\left\(\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}\\bm\{y\}\_\{n\-1\}\-y\_\{n\}\\right\)\(s25b\)=𝜶~n−1\+e~ns~n​𝒄~n,\\displaystyle=\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\}\+\\frac\{\\tilde\{e\}\_\{n\}\}\{\\tilde\{s\}\_\{n\}\}\\tilde\{\\bm\{c\}\}\_\{n\},\(s25c\)where

e~n=yn−𝝎~n⊤​𝒚n−1∈ℝ\.\\tilde\{e\}\_\{n\}=y\_\{n\}\-\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}\\bm\{y\}\_\{n\-1\}\\in\\mathbb\{R\}\.\(s26\)
Finally,

𝒄~n\\displaystyle\\tilde\{\\bm\{c\}\}\_\{n\}=𝒂n−𝑨n−1⊤​𝛀~n−1​𝑨n−1​𝒌~n−1​\(𝒙n\)\\displaystyle=\\bm\{a\}\_\{n\}\-\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\(s27a\)=𝒂n−𝑪~n−1​𝒌~n−1​\(𝒙n\),\\displaystyle=\\bm\{a\}\_\{n\}\-\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\),\(s27b\)e~n\\displaystyle\\tilde\{e\}\_\{n\}=yn−𝒌~n−1​\(𝒙n\)⊤​𝑨n−1⊤​𝛀~n−1​𝒚n−1\\displaystyle=y\_\{n\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{y\}\_\{n\-1\}\(s27c\)=yn−𝒌~n−1​\(𝒙n\)⊤​𝜶~n−1,and\\displaystyle=y\_\{n\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\},\\text\{ and \}\(s27d\)s~n\\displaystyle\\tilde\{s\}\_\{n\}=𝒌~n−1​\(𝒙n\)⊤​𝒂n\+σ2−𝒌~n−1​\(𝒙n\)​𝑨n−1⊤​𝛀~n−1​𝑨n−1​𝒌~n−1​\(𝒙n\)\\displaystyle=\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{a\}\_\{n\}\+\\sigma^\{2\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\(s27e\)=𝒌~n−1​\(𝒙n\)⊤​𝒂n\+σ2−𝒌~n−1​\(𝒙n\)⊤​𝑪~n−1​𝒌~n−1​\(𝒙n\)\.\\displaystyle=\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{a\}\_\{n\}\+\\sigma^\{2\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\.\(s27f\)

### S3\.2Case𝒟n=𝒟n−1∪\{𝒙n\}\\mathcal\{D\}\_\{n\}=\\mathcal\{D\}\_\{n\-1\}\\cup\\\{\\bm\{x\}\_\{n\}\\\}\(meaning thatδn\>ν\\delta\_\{n\}\>\\nu\)

In this case,

𝑲~n\\displaystyle\\tilde\{\\bm\{K\}\}\_\{n\}=\[𝑲~n−1𝒌~n−1​\(𝒙n\)𝒌~n−1​\(𝒙n\)⊤κ⁡\(𝒙n,𝒙n\)\]∈ℝDn×Dn,\\displaystyle=\\begin\{bmatrix\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}&\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ \\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}&\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\\end\{bmatrix\}\\in\\mathbb\{R\}^\{D\_\{n\}\\times D\_\{n\}\},\(s28a\)𝑨n\\displaystyle\\bm\{A\}\_\{n\}=\[\[𝑨n−1𝟎n−1\]𝒂n\]∈ℝn×Dn,and\\displaystyle=\\begin\{bmatrix\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\end\{bmatrix\}\\\\ \\bm\{a\}\_\{n\}\\end\{bmatrix\}\\in\\mathbb\{R\}^\{n\\times D\_\{n\}\},\\text\{ and \}\(s28b\)𝒂n\\displaystyle\\bm\{a\}\_\{n\}=\[0⋯01\]∈\{0\}Dn−1×\{1\}⊂ℝDn\.\\displaystyle=\\begin\{bmatrix\}0&\\cdots&0&1\\end\{bmatrix\}\\in\\\{0\\\}^\{D\_\{n\}\-1\}\\times\\\{1\\\}\\subset\\mathbb\{R\}^\{D\_\{n\}\}\.\(s28c\)
Then,

𝛀~n−1\\displaystyle\\tilde\{\\bm\{\\Omega\}\}\_\{n\}^\{\-1\}=\[\[𝑨n−1𝟎n−1\]𝒂n\]​\[𝑲~n−1𝒌~n−1​\(𝒙n\)𝒌~n−1​\(𝒙n\)⊤κ⁡\(𝒙n,𝒙n\)\]​\[\[𝑨n−1⊤𝟎n−1⊤\]𝒂n\]\+\[σ2​𝑰n−1𝟎n−1𝟎n−1⊤σ2\]\\displaystyle=\\begin\{bmatrix\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\end\{bmatrix\}\\\\ \\bm\{a\}\_\{n\}\\end\{bmatrix\}\\begin\{bmatrix\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}&\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ \\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}&\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\\end\{bmatrix\}\\begin\{bmatrix\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}&\\bm\{a\}\_\{n\}\\end\{bmatrix\}\+\\begin\{bmatrix\}\\sigma^\{2\}\\bm\{I\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&\\sigma^\{2\}\\end\{bmatrix\}\(s29a\)=\[𝑨n−1​𝑲~n−1𝑨n−1​𝒌~n−1​\(𝒙n\)𝒌~n−1​\(𝒙n\)⊤κ⁡\(𝒙n,𝒙n\)\]​\[\[𝑨n−1⊤𝟎n−1⊤\]𝒂n\]\+\[σ2​𝑰n−1𝟎n−1𝟎n−1⊤σ2\]\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}&\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ \\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}&\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\\end\{bmatrix\}\\begin\{bmatrix\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}&\\bm\{a\}\_\{n\}\\end\{bmatrix\}\+\\begin\{bmatrix\}\\sigma^\{2\}\\bm\{I\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&\\sigma^\{2\}\\end\{bmatrix\}\(s29b\)=\[𝑨n−1​𝑲~n−1​𝑨n−1⊤\+σ2​𝑰n−1𝑨n−1​𝒌~n−1​\(𝒙n\)𝒌~n−1​\(𝒙n\)⊤​𝑨n−1⊤κ⁡\(𝒙n,𝒙n\)\+σ2\]\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\+\\sigma^\{2\}\\bm\{I\}\_\{n\-1\}&\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ \\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{A\}\_\{n\-1\}^\{\\top\}&\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\+\\sigma^\{2\}\\end\{bmatrix\}\(s29c\)=\[𝛀~n−1−1𝑨n−1​𝒌~n−1​\(𝒙n\)\(𝑨n−1​𝒌~n−1​\(𝒙n\)\)⊤κ⁡\(𝒙n,𝒙n\)\+σ2\],\\displaystyle=\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}^\{\-1\}&\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ \\left\(\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\right\)^\{\\top\}&\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\+\\sigma^\{2\}\\end\{bmatrix\},\(s29d\)and applying the partitioned \(symmetric and positive semi\-definite\) matrix inverse formula\[[10](https://arxiv.org/html/2609.16472#bib.bib4), Ch\. 0\.7\]we get

𝛀~n=\[𝛀~n−1𝟎n−1𝟎n−1⊤0\]\+1s~n​\[𝝎~n−1\]​\[𝝎~n⊤−1\],\\tilde\{\\bm\{\\Omega\}\}\_\{n\}=\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\},\(s30\)where

𝝎~n\\displaystyle\\tilde\{\\bm\{\\omega\}\}\_\{n\}=𝛀~n−1​𝑨n−1​𝒌~n−1​\(𝒙n\)∈ℝn−1,and\\displaystyle=\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\in\\mathbb\{R\}^\{n\-1\},\\text\{ and \}\(s31a\)s~n\\displaystyle\\tilde\{s\}\_\{n\}=κ⁡\(𝒙n,𝒙n\)\+σ2−𝝎~n⊤​𝑨n−1​𝒌~n−1​\(𝒙n\)∈ℝ\.\\displaystyle=\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\+\\sigma^\{2\}\-\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\in\\mathbb\{R\}\.\(s31b\)
From here, we are able to expand

𝑪~n\\displaystyle\\tilde\{\\bm\{C\}\}\_\{n\}=\[\[𝑨n−1⊤𝟎n−1⊤\]𝒂n\]⁡\(\[𝛀~n−1𝟎n−1𝟎n−1⊤0\]\+1s~n​\[𝝎~n−1\]​\[𝝎~n⊤−1\]\)​\[\[𝑨n−1𝟎n−1\]𝒂n\]\\displaystyle=\\begin\{bmatrix\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}&\\bm\{a\}\_\{n\}\\end\{bmatrix\}\\left\(\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\}\\right\)\\begin\{bmatrix\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\end\{bmatrix\}\\\\ \\bm\{a\}\_\{n\}\\end\{bmatrix\}\(s32a\)=\[𝑨n−1⊤𝟎n−1⊤\]​𝛀~n−1​\[𝑨n−1𝟎n−1\]\+1s~n​\(\[𝑨n−1⊤𝟎n−1⊤\]​𝝎~n−𝒂n\)​\(𝝎~n⊤​\[𝑨n−1𝟎n−1\]−𝒂n⊤\)\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\end\{bmatrix\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\left\(\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\-\\bm\{a\}\_\{n\}\\right\)\\left\(\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\end\{bmatrix\}\-\\bm\{a\}\_\{n\}^\{\\top\}\\right\)\(s32b\)=\[𝑪~n−1𝟎Dn−1𝟎Dn−1⊤0\]\+1s~n​𝒄~n​𝒄~n⊤\\displaystyle=\\begin\{bmatrix\}\\tilde\{\\bm\{C\}\}\_\{n\-1\}&\\bm\{0\}\_\{D\_\{n\}\-1\}\\\\ \\bm\{0\}\_\{D\_\{n\}\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\tilde\{\\bm\{c\}\}\_\{n\}\\tilde\{\\bm\{c\}\}\_\{n\}^\{\\top\}\(s32c\)where

𝒄~n=𝒂n−\[𝑨n−1⊤𝟎n−1⊤\]​𝒘~n∈ℝDn\.\\tilde\{\\bm\{c\}\}\_\{n\}=\\bm\{a\}\_\{n\}\-\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}\\tilde\{\\bm\{w\}\}\_\{n\}\\in\\mathbb\{R\}^\{D\_\{n\}\}\.\(s33\)
Similarly,

𝜶~n\\displaystyle\\tilde\{\\bm\{\\alpha\}\}\_\{n\}=\[\[𝑨n−1⊤𝟎n−1⊤\]𝒂n\]⁡\(\[𝛀~n−1𝟎n−1𝟎n−1⊤0\]\+1s~n​\[𝝎~n−1\]​\[𝝎~n⊤−1\]\)​\[𝒚n−1yn\]\\displaystyle=\\begin\{bmatrix\}\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}&\\bm\{a\}\_\{n\}\\end\{bmatrix\}\\left\(\\begin\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}&\\bm\{0\}\_\{n\-1\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}&0\\end\{bmatrix\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\\\\ \-1\\end\{bmatrix\}\\begin\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}&\-1\\end\{bmatrix\}\\right\)\\begin\{bmatrix\}\\bm\{y\}\_\{n\-1\}\\\\ y\_\{n\}\\end\{bmatrix\}\(s34a\)=\[𝑨n−1⊤𝟎n−1⊤\]​𝛀~n−1​𝒚n−1\+1s~n​\(\[𝑨n−1⊤𝟎n−1⊤\]​𝝎~n−𝒂n\)​\(𝝎~n⊤​𝒚n−1−yn\)\\displaystyle=\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{y\}\_\{n\-1\}\+\\frac\{1\}\{\\tilde\{s\}\_\{n\}\}\\left\(\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\\\ \\bm\{0\}\_\{n\-1\}^\{\\top\}\\end\{bmatrix\}\\tilde\{\\bm\{\\omega\}\}\_\{n\}\-\\bm\{a\}\_\{n\}\\right\)\\left\(\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}\\bm\{y\}\_\{n\-1\}\-y\_\{n\}\\right\)\(s34b\)=\[𝜶~n−10\]\+e~ns~n​𝒄~n,\\displaystyle=\\begin\{bmatrix\}\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\}\\\\ 0\\end\{bmatrix\}\+\\frac\{\\tilde\{e\}\_\{n\}\}\{\\tilde\{s\}\_\{n\}\}\\tilde\{\\bm\{c\}\}\_\{n\},\(s34c\)where

e~n=yn−𝝎~n⊤​𝒚n−1∈ℝ\.\\tilde\{e\}\_\{n\}=y\_\{n\}\-\\tilde\{\\bm\{\\omega\}\}\_\{n\}^\{\\top\}\\bm\{y\}\_\{n\-1\}\\in\\mathbb\{R\}\.\(s35\)
Finally,

𝒄~n\\displaystyle\\tilde\{\\bm\{c\}\}\_\{n\}=𝒂n−\[𝑨n−1⊤​𝛀~n−1​𝑨n−1​𝒌~n−1​\(𝒙n\)0\]\\displaystyle=\\bm\{a\}\_\{n\}\-\\begin\{bmatrix\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ 0\\end\{bmatrix\}\(s36a\)=\[−𝑪~n−1​𝒌~n−1​\(𝒙n\)1\],\\displaystyle=\\begin\{bmatrix\}\-\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\\\\ 1\\end\{bmatrix\},\(s36b\)e~n\\displaystyle\\tilde\{e\}\_\{n\}=yn−𝒌~n−1​\(𝒙n\)⊤​𝑨n−1⊤​𝛀~n−1​𝒚n−1\\displaystyle=y\_\{n\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{y\}\_\{n\-1\}\(s36c\)=yn−𝒌~n−1​\(𝒙n\)⊤​𝜶~n−1,and\\displaystyle=y\_\{n\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\tilde\{\\bm\{\\alpha\}\}\_\{n\-1\},\\text\{ and \}\(s36d\)s~n\\displaystyle\\tilde\{s\}\_\{n\}=κ⁡\(𝒙n,𝒙n\)\+σ2−𝒌~n−1​\(𝒙n\)⊤​𝑨n−1⊤​𝛀~n−1​𝑨n−1​𝒌~n−1​\(𝒙n\)\\displaystyle=\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\+\\sigma^\{2\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\bm\{A\}\_\{n\-1\}^\{\\top\}\\tilde\{\\bm\{\\Omega\}\}\_\{n\-1\}\\bm\{A\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\(s36e\)=κ⁡\(𝒙n,𝒙n\)\+σ2−𝒌~n−1​\(𝒙n\)⊤​𝑪~n−1​𝒌~n−1​\(𝒙n\)\.\\displaystyle=\\kappa\(\\bm\{x\}\_\{n\},\\bm\{x\}\_\{n\}\)\+\\sigma^\{2\}\-\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)^\{\\top\}\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)\.\(s36f\)
It is worth noting that by applying the partitioned matrix inverse formula\[[10](https://arxiv.org/html/2609.16472#bib.bib4), Ch\. 0\.7\]on \([s28a](https://arxiv.org/html/2609.16472#S3.E28.1)\) one can get

𝑲~n−1−1=1δn​\[δn​𝑲~n−1−1\+𝒂^n​𝒂^n⊤−𝒂^n−𝒂^n⊤1\],\\tilde\{\\bm\{K\}\}\_\{n\-1\}^\{\-1\}=\\frac\{1\}\{\\delta\_\{n\}\}\\begin\{bmatrix\}\\delta\_\{n\}\\tilde\{\\bm\{K\}\}\_\{n\-1\}^\{\-1\}\+\\hat\{\\bm\{a\}\}\_\{n\}\\hat\{\\bm\{a\}\}\_\{n\}^\{\\top\}&\-\\hat\{\\bm\{a\}\}\_\{n\}\\\\ \-\\hat\{\\bm\{a\}\}\_\{n\}^\{\\top\}&1\\end\{bmatrix\},\(s37\)withδn\\delta\_\{n\}and𝒂^n\\hat\{\\bm\{a\}\}\_\{n\}as in \([s12](https://arxiv.org/html/2609.16472#S2.E12)\)\.

### S3\.3Comments on the computational and memory complexity

The per\-step computational and memory complexity of Algorithm[1](https://arxiv.org/html/2609.16472#alg1)is𝒪⁡\(Dn2\+r​Dn\)\\mathcal\{O\}\(D\_\{n\}^\{2\}\+rD\_\{n\}\), whereDnD\_\{n\}is the ALD sparsification dictionary size, andrris the number of warping parameters\.

The following operations dominate this complexity:

- •Evaluating the kernel vector𝒌~n−1​\(𝒙n\)\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)is𝒪⁡\(Dn\)\\mathcal\{O\}\(D\_\{n\}\), since it evaluates the kernel between𝒙n\\bm\{x\}\_\{n\}and each of theDnD\_\{n\}dictionary elements\.
- •Computing the ALD projection coefficients𝒂^n=𝑲~n−1−1​𝒌~n−1​\(𝒙n\)\\hat\{\\bm\{a\}\}\_\{n\}=\\tilde\{\\bm\{K\}\}^\{\-1\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)is𝒪⁡\(Dn2\)\\mathcal\{O\}\(D\_\{n\}^\{2\}\)since𝑲~n−1\\tilde\{\\bm\{K\}\}^\{\-1\}\_\{n\}is updated via a rank\-1 block inversion\.
- •Evaluating the warping transformation used in Sec\.[5](https://arxiv.org/html/2609.16472#S5)g⁡\(yn,𝜽\)g\(y\_\{n\};\\bm\{\\theta\}\)\(and its gradient\) is𝒪⁡\(r\)\\mathcal\{O\}\(r\)due to basic arithmetic andtanh\\tanhoperations\.
- •Computing the matrix\-vector products𝑩~n−1​𝒌~n−1​\(𝒙n\)\\tilde\{\\bm\{B\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)and𝑪~n−1​𝒌~n−1​\(𝒙n\)\\tilde\{\\bm\{C\}\}\_\{n\-1\}\\tilde\{\\bm\{k\}\}\_\{n\-1\}\(\\bm\{x\}\_\{n\}\)is𝒪⁡\(r​Dn\)\\mathcal\{O\}\(rD\_\{n\}\)and𝒪⁡\(Dn2\)\\mathcal\{O\}\(D^\{2\}\_\{n\}\), respectively, due to the matrix and vector dimensions\.
- •Finally, updating𝑪~n\\tilde\{\\bm\{C\}\}\_\{n\}and𝑩~n\\tilde\{\\bm\{B\}\}\_\{n\}is𝒪⁡\(Dn2\)\\mathcal\{O\}\(D\_\{n\}^\{2\}\)and𝒪⁡\(r​Dn\)\\mathcal\{O\}\(rD\_\{n\}\), respectively, due to the rank\-1 outer products𝒄~n​𝒄~n⊤\\tilde\{\\bm\{c\}\}\_\{n\}\\tilde\{\\bm\{c\}\}\_\{n\}^\{\\top\}and𝒃~n​𝒄~n⊤\\tilde\{\\bm\{b\}\}\_\{n\}\\tilde\{\\bm\{c\}\}^\{\\top\}\_\{n\}in their update rules\.