Incremental Structure Discovery of Classification via Sequential Monte Carlo

Let $\mathcal{D}:=(\mathbf{X},\mathbf{y})$ be a dataset for binary classification, where $\mathbf{X}:=[\mathbf{x}_{1},\ldots,\mathbf{x}_{n}]$ is the $n$ data points and $\mathbf{y}:=[y_{1},\ldots,y_{n}]$ is the corresponding labels, i.e., $\mathbf{x}_{i}\in\mathcal{X}$ and $y_{i}\in\{0,1\}$ for each $i\in[n]$, where $\mathcal{X}$ stands for the feature space of data points. GPC typically samples a function $f:\mathcal{X}\to\mathbb{R}$ from a Gaussian Process prior $f\sim\mathrm{GP}(m,k)$, where $m:\mathcal{X}\to\mathbb{R}$ is a mean function and $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ is a covariance function, i.e., a kernel [28]. The probabilistic model for classification $P_{\mathbf{X}}(\mathbf{y};m,k)$ is then be formulated by

where $\sigma:\mathbb{R}\to[0,1]$ is a sigmoid function, e.g., the logistic function or the probit function. To infer the label $y_{*}$ on a new data point $\mathbf{x}_{*}$, i.e., to reason about the posterior distribution $P_{\mathbf{X},\mathbf{x}_{*}}(y_{*}=1\mid\mathbf{y})$, we treat $f(\mathbf{X})$ as a vector $\mathbf{f}:=[f_{1},\ldots,f_{n}]$ of latent variables and carry out the inference in two steps: (i) first derive the distribution of the latent variable $f_{*}$ as

and (ii) then derive the posterior of $y_{*}$ by integrating out $f_{*}$ as

Note that the key part is to account for the posterior distribution of latent variables $\mathbf{P}_{\mathbf{X}}(\mathbf{f}\mid k,\mathbf{y})$.

Usually, the mean $m$ is pre-determined to be the constant-zero function, whereas the kernel $k$ is parameterized by a vector $\theta=[\theta_{1},\ldots,\theta_{d(k)}]\in\Theta_{k}\subseteq\mathbb{R}^{d(k)}$, where $d(k)$ is the number of real-valued parameters in $k$. Let us denote $k_{\theta}:=(k,\theta)$ and reuse the same notation for the actual covariance function derived from $k$ and $\theta$. While many GP-based methods pre-determine $k$ and treat $\theta$ as hyper-parameters, in this paper, we aim to characterize both $k$ and $\theta$ as latent information of the classification model, as well as develop a method that is adaptive in both the structure and parameters of the kernel in the following sections.

To allow a rich prior distribution over kernel structures $k$, we define a sample space $\mathcal{K}$ using a probabilistic context-free grammar (PCFG), following the practice of a line of prior work [9, 26, 25]:

The non-terminal $B$ stands for an extensible collection of basic kernels. Two kernels can be combined with a binary operator $\oplus$ that computes the pointwise addition or multiplication of the two kernels. The PCFG also assigns a probability to each production rule in (5)-(7), thus formulating a prior distribution on $\mathcal{K}$. The meaning of a kernel expression $k\in\mathcal{K}$ is defined inductively as follows.

With the PCFG, we extend the standard GPC model (c.f. (1)-(2)) to treat both the structure $k$ and the parameters $\theta$ of the kernel as latent variables, resulting in a probabilistic model $P_{\mathbf{X}}(k,\theta,\epsilon,\mathbf{y})$:

The latent variable $\epsilon$ stands for the noise in the GP. Note that the model samples kernel parameters from the standard Normal distribution, but basic kernels may impose constraints on its parameters. We apply standard transformations to obtain constrained parameters and omit the details here; for example, $z\mapsto\exp(z)$ for positive parameters and $z\mapsto 2/(1+\exp(z))$ for parameters in $(0,2)$.

We further apply a standard reparameterization trick to the model above via Cholesky decomposition. That is, instead of sampling from a multivariate Normal distribution, we sample $n+1$ i.i.d. auxiliary variables from the standard Normal distribution and compute the latent vector $\mathbf{f}$ as follows.

The term $\mathrm{Cholesky}(\mathbf{M})$ performs Cholesky decomposition of the covariance matrix $\mathbf{M}$, i.e., computes a lower-triangular matrix $\mathbf{L}$ such that $\mathbf{L}\mathbf{L}^{\top}=\mathbf{M}$. The reparmeterized model specifies a joint distribution $P_{\mathbf{X}}(k,\theta,\epsilon,\beta,\boldsymbol{\eta})$. To simplify the notation, we define $\phi:=(\theta,\epsilon)$ and $\psi:=(\beta,\boldsymbol{\eta})$ to denote the kernel parameters and auxiliary variables, respectively.

By the Bayes’ rule, the posterior distribution on $k,\phi,\psi$ given a dataset $(\mathbf{X},\mathbf{y})$ is given by

The goal of our method is to generate and maintain a finite set of $M\geq 1$ weighted particles

each of which consists a weight $w^{i}>0$ and a tuple of latent variables $(k^{i},\phi^{i},\psi^{i})$. These particles are intended to approximate the posterior distribution given in (21), so as to compute the expectation of some test function $\tau$ with respect to the posterior distribution as

Recall that in (3)-(4) we reviewed how to condition a standard GPC model on a dataset to make predictions. We generalize the idea to define a test function $\tau^{\mathsf{prob}}_{\mathbf{x}_{*}}$ to compute the classification probability of a new data point $\mathbf{x}_{*}$, given all the latent information $k,\phi,\psi$:

Different from (3), the posterior distribution on $f_{*}$ becomes analytically tractable because the latent vector $\mathbf{f}$ is determined by the given latent information $k,\phi,\psi$. To see that, we derive and simplify the posterior distribution on $f_{*}$ as

which is a Normal distribution because of the multivariate Normal joint distribution (c.f. (15)). However, (24) might not be analytically tractable due to the choice of the sigmoid function $\sigma$. Fortunately, to predict the label for a new data point $\mathbf{x}_{*}$, (24) is essentially a uni-variate integral, so we resort to use Monte Carlo estimation. In some case, e.g., when $\sigma$ is the probit function, the integral has a closed-form solution so that $\tau^{\mathsf{prob}}_{\mathbf{x}_{*}}$ can be evaluated easily.

Algorithm 1 shows our method of applying SMC to adaptive and incremental learning on the GPC models present in Section 2. It follows a standard reweight-resample-rejuvenate pipeline for implementing SMC samplers [4]. Our method reweights and resamples the particle set $\{(k^{i}_{j},\phi^{i},\psi^{i}_{j})\}_{i=1}^{M}$ as shown in line 1 to line 1. Line 1 and line 1 are the process of reweighting; it implicitly chooses the forward Markov kernel $K$ to be the identity kernel and the backward Markov kernel $L$ to be its time reversal. As a result, it calculates the new weight for each particle based on the joint probability of model’s latent variables $(k,\phi,\psi)$ and the first $j$ data points $(\mathbf{X}_{1:j},\mathbf{y}_{1:j})$. The process of resampling in lines 1 to 1 adopts a standard machinery of adaptive resampling, i.e., only initiate the resample process when the effective sample size (ESS) drops below a threshold. Finally, lines 1 to 1 implement the rejuvenation loop. In each rejuvenation iteration, we follow the design of AutoGP [25] to first evolve the kernel structure of each particle via IMCMC (reviewed in the previous section) and then apply HMC to evolve the real-valued latent kernel parameters $\phi$ and auxiliary variables $\psi$.

Algorithm 1 can already to applied to the online setting, in the sense that the step $j$ stands for the order when a data point arrives. In the offline setting, although algorithm 1 can also work, we can make it more effective by allowing batch tempering, i.e., we can divide a dataset $\mathcal{D}=(\mathbf{X},\mathbf{y})$ into $T$ batches, each of which contains $n/T$ data points $\mathcal{D}_{j}=(\mathbf{X}_{j\cdot n/T+1:(j+1)\cdot n/T},\mathbf{y}_{j\cdot n/T+1:(j+1)\cdot n/T})$, for $j=0,\ldots,T-1$. Then in the iteration for step $j$, instead of incorporating one single data point, we use the whole batch $\mathcal{D}_{j}$ in the reweighting and rejuvenation processes.

With the particle set $\{(w^{i},(k^{i},\phi^{i},\psi^{i}))\}_{i=1}^{M}$ that approximates the posterior distribution $P_{\mathbf{X}}(k,\phi,\psi\mid\mathbf{y})$, we can use it to make predictions for new data points. Following (23)-(25), for a new data point $\mathbf{x}_{*}$, we approximately compute the predictive probability of $y_{*}$ being $1$ as

where $K$ is the number of Monte Carlo estimations to approximate $\tau^{\mathsf{prob}}_{\mathbf{x}_{*}}$ and $f^{i}_{1},\ldots,f^{i}_{K}$ are $K$ i.i.d. samples from the posterior $P_{\mathbf{X},\mathbf{x}_{*}}(\cdot\mid\mathbf{f}^{i})$ for each particle $i$.

The toy datasets used in our experiments are generated by Scikit-Learn [23]. Those toy datasets with 2-dimensional inputs are easy to be visualized with weights; thus, they provide an intuitive way to illustrate the characterization of our method.

For RQ1, the first experiment is designed to demonstrate the automatic adaptability of our method in selecting appropriate kernels for varying datasets. Figure 1 illustrates the comparison between our method and GPs using pre-selected kernels—Linear and SquaredExponential—across three datasets. GPs used in this experiment are implemented using the same machinery as Algorithm 1, except that the kernel structure $k$ is fixed. This comparison demonstrates that our method can exhibit characteristics of different kernels.

To take a closer look at on how our method combines properties of different kernels, Figure 2 plots the kernel’s behavior of different particles after running Algorithm 1 on a dataset that can be linearly separated. It shows that our method could discover different kernel structures with similar performance; thus, the method would be robust to incorporate future unseen data points.

Furthermore, to evaluate our method’s performance on pattern shifts and simulate an online setting described in RQ2, we initially run our method to learn on a portion of the dataset, followed by subsequently learning on the remaining data in another batch. Figure 3 shows that our method is able to adjust kernel structures to adapt to the shift of pattern in an online setting.

In summary, experiments on the toy datasets show that our method can discover various kernel structures, as well as learn and adapt to pattern shifts in the online setting.

To test our method in practice, we select three real-world datasets in different domain. Some datasets [8, 31] were used in previous studies [17]. The details of datasets can be found at Appendix.

We use two GPs with pre-determined kernels, Random Forest, and three other classification methods for evaluating RQ1. The two GPs are set up in the same manner as Section 4.1. We use a Julia implementation of Random Forest [27]. To further test the effectiveness of our method, we select three methods provided by Scikit-Learn [23]. We choose Passive Aggressive Classifier (PAC) [5], Stochastic Gradient Descent (SGD) [32], and Naive Bayes in this experiment.

In order to test the ability of online learning for RQ2, we adapt a similar setting to Section 4.1 but make it more extreme and biased: in the first batch of learning, we only include one class of the data and then incorporate batches of remaining data with the other class.

Results are shown in Table 1. Accuracy of each offline task is calculated in a standard way and reported. For online task, average accuracy of each task is calculated by $\frac{1}{B}\sum^{B}_{i=1}A_{i}$, where $B$ is the total number of batches and $A_{i}$ is accuracy after each batch of learning. For the offline setting, our method outperforms all other methods. It shows that our method achieves the goal of reducing requirement of prior knowledge by automatic selection of kernel. For the online setting, our method achieves highest accuracy on two out of three datasets. The results of the online setting demonstrate our method is able to adapt kernel structures in an incremental manner.

In summary, experiments on the real-world datasets show that our method has ability to seamlessly integrate features of different kernels and accurately adjust itself against pattern shifts within the dataset.