# Bayesian PLS!
Approximate Bayes Optimal Pseudo-Label Selection (PLS)

Julian Rodemann Department of Statistics
Ludwig-Maximilians-Universität (LMU)
Munich, Germany  Jann Goschenhofer  Department of Statistics
Ludwig-Maximilians-Universität (LMU)
Munich, Germany  Munich Center for Machine Learning (MCML)
Munich, Germany  Fraunhofer Institute of Integrated Circuits (IIS)
Erlangen, Germany  Emilio Dorigatti  Department of Statistics
Ludwig-Maximilians-Universität (LMU)
Munich, Germany  Munich Center for Machine Learning (MCML)
Munich, Germany  Institute of Computational Biology
Helmholtz Zentrum
Neuherberg, Germany  Thomas Nagler  Department of Statistics
Ludwig-Maximilians-Universität (LMU)
Munich, Germany  Munich Center for Machine Learning (MCML)
Munich, Germany  Thomas Augustin  Department of Statistics
Ludwig-Maximilians-Universität (LMU)
Munich, Germany

###### Abstract

Semi-supervised learning by self-training heavily relies on pseudo-label
selection (PLS). The selection often depends on the initial model fit on
labeled data. Early overfitting might thus be propagated to the final model by
selecting instances with overconfident but erroneous predictions, often
referred to as confirmation bias. This paper introduces BPLS, a Bayesian
framework for PLS that aims to mitigate this issue. At its core lies a
criterion for selecting instances to label: an analytical approximation of the
posterior predictive of pseudo-samples. We derive this selection criterion by
proving Bayes optimality of the posterior predictive of pseudo-samples. We
further overcome computational hurdles by approximating the criterion
analytically. Its relation to the marginal likelihood allows us to come up
with an approximation based on Laplace’s method and the Gaussian integral. We
empirically assess BPLS for parametric generalized linear and non-parametric
generalized additive models on simulated and real-world data. When faced with
high-dimensional data prone to overfitting, BPLS outperforms traditional PLS
methods.111Code available at: https://anonymous.4open.science/r/Bayesian-pls

##  1 INTRODUCTION

Labeled data are scarce in many learning settings. This can be due to a
variety of reasons such as restrictions on time, knowledge, or financial
resources. Unlabeled data, however, are often much more accessible. This has
given rise to the paradigm of semi-supervised learning (SSL), where
information from unlabeled data is integrated into model training to improve
predictions in a supervised learning framework. Within SSL, an intuitive and
widely used approach is referred to as self-training or pseudo-labeling [Shi
et al., 2018, Lee et al., 2013, McClosky et al., 2006]. The idea is to fit an
initial model to labeled data and iteratively assign pseudo-labels to
unlabeled data according to the model’s predictions. The latter requires a
criterion (sometimes called confidence measure) for pseudo-label selection
(PLS), that is, the selection of instances to be pseudo-labeled and added to
the training data.

By design, self-training strongly relies on the initial model fit and the way
instances are selected to be pseudo-labeled. Everything hinges upon the
interplay between the selection criterion and the initial model’s
generalization performance. If the initial model generalizes poorly, initial
misconceptions can propagate throughout the process, only making things worse.
High-dimensional data prone to overfitting are particularly sensitive to such
confirmation bias [Arazo et al., 2020]. Usually, self-training’s sweet spot
lies somewhere else: When the labeled data allow the model to learn
sufficiently well while still leaving some room for improvement. Generally,
the poorer the initial generalization, the harder it is to select sensible
pseudo-labels to improve generalization, i.e., the more crucial the role of
the selection criterion. Note that SSL is applied to data with high shares
(typically over $80\%$ [Sohn et al., 2020, Arazo et al., 2020]) of unlabeled
data, where initial overfitting is likely for high-dimensional models, while
final overfitting is not.

###  1.1 Motivation

Accordingly, we strive for a selection criterion that is robust with respect
to the initial model fit, i.e., its learned parameters. At the same time, it
should still exploit the information in the labeled data. Such a measure calls
for disentangling the uncertainty contributions of the data and the model’s
parameters. This is in line with recent work in uncertainty quantification
(UQ) that suggests decomposing epistemic uncertainty into approximation
uncertainty driven by (a lack of) data and modeling uncertainty driven by
(primarily parametric) assumptions [Hüllermeier and Waegeman, 2021]. Bayesian
inference offers a sound and consistent framework for this distinction. Its
rationale of technically modeling not only data but also parameters as random
variables has proven to offer much insight into UQ for machine learning
[Hüllermeier and Waegeman, 2021] and deep learning [Abdar et al., 2021,
Malinin and Gales, 2018].

We exploit the Bayesian framework for pinpointing uncertainty with regard to
data and parameters in PLS. Our approach of Bayesian pseudo-label selection
(BPLS) enables us to choose pseudo-labels that are likely given the observed
labeled data but not necessarily likely given the estimated parameters of the
fitted model. What is more, BPLS allows to include prior information not only
for predicting but also for selecting pseudo-labels. Notably, BPLS is flexible
enough to be applied to any kind of predictive model whose likelihood and
Fisher-information are accessible, including non-Bayesian models. BPLS entails
a Bayes optimal selection criterion, the pseudo posterior predictive (PPP).
Its intuition is straightforward yet effective: By averaging over all
parameter values, PPP is more robust towards the initial fit compared to the
predictive distribution based on a single optimal parameter vector. Our
approximate version of the PPP is simple and computationally cheap to
evaluate:
$\ell(\hat{\theta})-\frac{1}{2}\log\lvert\mathcal{I(\hat{\theta})}\rvert$ with
$\ell(\hat{\theta})$ being the likelihood and $\mathcal{I(\hat{\theta})}$ the
Fisher-information matrix at the fitted parameter vector $\hat{\theta}$. As an
approximation of the joint PPP, it does not require an i.i.d. assumption,
rendering it applicable to a wide range of applied learning setups.

###  1.2 Main Contributions

(1) We derive PPP by formalizing PLS as a decision problem and show222Proofs
of all theorems in this paper can be found in the supplementary material. that
PPP corresponds to the Bayes criterion, rendering selecting instances with
regard to it Bayes optimal, see sections 2.1 and 2.2.

(2) Since our selection criterion includes a possibly intractable integral, we
provide analytical approximations, exploiting Laplace’s method and the
Gaussian integral, both for uninformative and informative priors. Using
varying levels of accuracy, we balance the trade-off between computational
feasibility and precision, see Section 3.

(3) We provide empirical evidence333 Implementations of the proposed methods
as well as reproducible scripts for the experiments are provided in the
anonymous repository named Bayesian-pls (Bayesian, please!), see abstract.
for BPLS’ superiority over traditionally predominant PLS methods in case of
semi-supervised generalized additive models (GAMs) and generalized linear
models (GLMs) faced with high-dimensional data prone to overfitting, see
Section 4.

##  2 BAYESIAN PLS

Most semi-supervised methods deal with classification or clustering tasks [Van
Engelen and Hoos, 2020, Chapelle et al., 2006]. Loosely leaning on [Triguero
et al., 2015], we formalize SSL as follows. Consider labeled data

$\mathcal{D}=\left\\{\left(x_{i},y_{i}\right)\right\\}_{i=1}^{n}\in\left(\mathcal{X}\times\mathcal{Y}\right)^{n}$
(1)

and unlabeled data

$\mathcal{U}=\left\\{\left(x_{i},\mathcal{Y}\right)\right\\}_{i=n+1}^{m}\in\left(\mathcal{X}\times
2^{\mathcal{Y}}\right)^{m-n}$ (2)

from the same data generation process, where $\mathcal{X}$ is the feature
space and $\mathcal{Y}$ is the categorical target space. The aim of SSL is to
learn a predictive classification function $f$ such that
$f(x)=\hat{y}\in\mathcal{Y}$ utilizing both $\mathcal{D}$ and $\mathcal{U}$.
As is customary in self-training, we start by fitting a model with unknown
parameter vector $\theta\in\Theta$, $\Theta$ compact with $\dim(\Theta)=q$, on
labeled data
$\mathcal{D}=\left\\{\left(x_{i},y_{i}\right)\right\\}_{i=1}^{n}$. Our goal is
– as usual – to learn the conditional distribution of $p(y\mid x)$ through
$\theta$ from observing features ${x}=(x_{1},\dots,x_{n})\in\mathcal{X}^{n}$,
and responses ${y}=(y_{1},\dots,y_{n})\in\mathcal{Y}^{n}$ in $\mathcal{D}$.
Adopting the Bayesian perspective, we can state a prior function over $\theta$
as $\pi(\theta)$. The prior can represent information on $\theta$ but may also
be uninformative.

Within existing frameworks for self-training (see Section 5) in SSL, one could
deploy such a Bayesian setting for predicting unknown labels of
$\mathcal{U}=\left\\{\left(x_{i},\mathcal{Y}\right)\right\\}_{i=1}^{m}$ as
well as for the final predictions on unseen test data. However, we aim at a
Bayesian framework for selecting pseudo-labels. This is beneficial for two
reasons. First and foremost, considering the Bayesian posterior predictive
distribution in PLS will turn out to be more robust towards the initial fit on
$\mathcal{D}$ than classical selection criteria. Second, the Bayesian engine
brings along the usual benefit of allowing to explicitly account for prior
knowledge when selecting instances to be labeled. Notably, our framework of
Bayesian pseudo-label selection is unrelated to how pseudo-labels are
predicted.

###  2.1 The Case for the Posterior Predictive in PLS

For any model with parameters $\theta\in\Theta$, the likelihood function for
observed features ${x}$ and labels ${y}$ is commonly defined as
$\mathcal{L}_{y\mid x}(\theta)=f_{\theta}(y\mid x),$ where $f_{\theta}(\cdot)$
is from a parameterized family of probability density functions. In the
Bayesian universe, parameters $\theta$ are more than just functional arguments
[Murphy, 2012]. They are random quantities themselves, allowing us to
condition on them: $\mathcal{L}_{y\mid x}(\theta)=p({y}\mid{x},\theta)$.
Recall that we have specified a prior $\pi(\theta)$ on the parameters
beforehand. After observing data, it can be updated to a posterior following
Bayes’ Theorem
$p(\theta\mid{y},{x})=p({y}\mid{x},\theta)\,\pi(\theta)/p({y}\mid{x}),$ where
the denominator is the marginal likelihood

$p({y}\mid{x})=\int_{\Theta}p({y}\mid{x},\theta)\,\pi(\theta)\,d\theta,$ (3)

or, more colloquially, Bayesian evidence [Lotfi et al., 2022, Barber, 2012].
For previously unseen data $(\tilde{y},\tilde{x})$, the posterior predictive
distribution is defined as

$p(\tilde{y}\mid\tilde{x},{y},{x})=\int_{\Theta}p(\tilde{y}\mid\tilde{x},\theta)\,p(\theta\mid{y},{x})\,d\theta.$
(4)

The posterior predictive closely resembles the marginal likelihood in case we
include $(\tilde{y},\tilde{x})$ in the data – a fact that we will exploit for
our approximations in Section 3. Both marginalize the likelihood over
$\theta$. The difference is the weight: The marginal likelihood integrates out
$\theta$ with regard to the prior, while the posterior predictive integrates
out $\theta$ with regard to the posterior. Accordingly, both can be considered
PLS criteria that are robust towards the initial fit: They average over all
possible $\theta$-values instead of relying on one estimated $\hat{\theta}$
from the trained model. 444The probabilistic interpretation of the marginal
likelihood – in the words of [Lotfi et al., 2022] – is: The probability that
we would generate a dataset with a model if we randomly sample from a prior
over its parameters. The posterior predictive, analogously, is the probability
that we would generate data with a model if we randomly sample from a
posterior over its parameters. Computational issues aside, the posterior
predictive of pseudo-labeled data thus encapsulates a perfectly natural
selection criterion for self-training: It selects pseudo-labels that are most
likely conditioned on the true observed $\mathcal{D}$, the assumed model and
all plausible parameters from the prior or posterior, respectively.

Both the data and the estimated parameters (as functions of the data) will
change throughout the process of self-training. We argue that conditioning the
choice of unlabeled instances solely on the estimated parameters in early
iterations over-emphasizes the influence of the initial model. This optimistic
reliance can be harmful in case of small $n$ and high $q$, where overfitting
is likely. Selecting instances by the posterior predictive mitigates this.

###  2.2 Bayes Optimality of Pseudo Posterior Predictive

In the following, we show that selecting pseudo-labels with regard to their
posterior predictive is Bayes optimal. We further show the same holds for
selection with regard to the marginal likelihood in case of a non-updated
prior. To this end, we formalize the selection of data points to be pseudo-
labeled as a canonical decision problem, where an action corresponds to the
selection of an instance from the set of unlabeled data $\mathcal{U}$.

######  Definition 1 (PLS as Decision Problem)

Consider the decision-theoretic triple $(\mathcal{U},\Theta,u(\cdot))$ with an
action space of unlabeled data555We assume absence of tied observations for
simplicity such that we can understand $\mathcal{U}$ as set. to be selected,
i.e., instances $(x_{i},\mathcal{Y})$ as actions, a space of unknown states of
nature (parameters) $\Theta$ and a utility function
$u:\mathcal{U}\times\Theta\to\mathbb{R}$.

Loosely inspired by [Cattaneo, 2007], we now define the utility of a selected
data point $(x_{i},\mathcal{Y})$ as the plausibility of being generated
jointly with $\mathcal{D}$ by a model with parameters $\theta\in\Theta$ if we
include it with pseudo-label $\hat{y}_{i}\in\mathcal{Y}$ (obtained through any
predictive model) in $\mathcal{D}\cup(x_{i},\hat{y}_{i})$. This is
incorporated by the likelihood of $\mathcal{D}\cup(x_{i},\hat{y}_{i})$, which
shall be called pseudo-label likelihood and written as
$p(\mathcal{D}\cup(x_{i},\hat{y}_{i})\mid\theta)$. We thus condition the
selection problem on a model class as well as on already predicted pseudo-
labels. The former conditioning is not required (see the extension in Section
6) for the well-definedness of the pseudo-likelihood while the latter is.

######  Definition 2 (Pseudo-Label Likelihood as Utility)

Let $(x_{i},\mathcal{Y})$ be any decision (selection) from $\mathcal{U}$. We
assign utility to each $(x_{i},\mathcal{Y})$ given $\mathcal{D}$ and pseudo-
labels $\hat{y}\in\mathcal{Y}$ by the following measurable utility function

$\displaystyle u\colon\mathcal{U}\times\Theta$ $\displaystyle\to\mathbb{R}$
$\displaystyle((x_{i},\mathcal{Y}),\theta)$ $\displaystyle\mapsto
u((x_{i},\mathcal{Y}),\theta)=p(\mathcal{D}\cup(x_{i},\hat{y}_{i})\mid\theta),$

which is said to be the pseudo-label likelihood.

This utility function is a natural probabilistic choice to assign utilities to
selected pseudo-labels given the predicted pseudo-labels. With a prior
$\pi(\theta)$, we get the following result.

###### Theorem 1

In the decision problem $(\mathcal{U},\Theta,u(\cdot))$ (Definition 1) with
the pseudo-label likelihood as utility function (Definition 2) and a prior
$\pi(\theta)$ on $\Theta$, the standard Bayes criterion

$\displaystyle\Phi(\cdot,\pi)\colon\mathcal{U}$ $\displaystyle\to\mathbb{R}$
$\displaystyle a$
$\displaystyle\mapsto\Phi(a,\pi)=\mathbb{E}_{\pi}(u(a,\theta))$

corresponds to the pseudo marginal likelihood
$p(\mathcal{D}\leavevmode\nobreak\ \cup\leavevmode\nobreak\
(x_{i},\hat{y}_{i}))$.

###### Corollary 1

For any prior $\pi(\theta)$ on $\Theta$, the action
$a_{m}^{*}=\operatorname*{arg\,max}_{i}p(\mathcal{D}\cup(x_{i},\hat{y}_{i}))$
is Bayes optimal.

Taking the observed labeled data $\mathcal{D}$ into account by updating the
prior $\pi(\theta)$ to a posterior $p(\theta\mid\mathcal{D})$, we end up with
an analogous result for the pseudo posterior predictive. The Theorem requires
only the Proposition by [Berger, 1985, section 4.4.1] stating that posterior
loss equals prior risk. That is, conditional Bayes optimality equals
unconditional Bayes optimality.

###### Theorem 2

In the decision problem $(\mathcal{U},\Theta,u(\cdot))$ and the pseudo-label
likelihood as utility function as in Theorem 1 but with the prior updated by
the posterior $\pi(\theta)=p(\theta\mid\mathcal{D})$ on $\Theta$, the standard
Bayes criterion
$\Phi(\cdot,\pi)\colon\mathcal{U}\to\mathbb{R};\,a\mapsto\Phi(a,\pi)=\mathbb{E}_{\pi}(u(a,\theta))$
corresponds to the pseudo posterior predictive
$p(\mathcal{D}\cup(x_{i},\hat{y}_{i})\mid\mathcal{D})$.

###### Corollary 2

Action
$a_{p}^{*}=\operatorname*{arg\,max}_{i}p(\mathcal{D}\cup(x_{i},\hat{y}_{i})\mid\mathcal{D})$
is Bayes optimal for any updated prior $\pi(\theta)=p(\theta\mid\mathcal{D})$.

Further note that directly maximizing the likelihood with regard to $a$
corresponds to the optimistic max-max-criterion, see Theorem 3.

###### Theorem 3

In the decision problem $(\mathcal{U},\Theta,u(\cdot))$ with the pseudo-label
likelihood as utility function as in Theorem 1, the max-max criterion

$\displaystyle\Phi(\cdot)\colon\mathcal{U}$ $\displaystyle\to\mathbb{R};$
$\displaystyle a$ $\displaystyle\mapsto\Phi(a)=\max_{\theta}(u(a,\theta))$

corresponds to the (full) likelihood.

The max-max-criterion advocates deciding for an action (here: selection of
pseudo-labeled data) with the highest utility (here: likelihood) according to
the most favorable state of nature $\theta$, e.g. see [Rapoport, 1998]. It can
hardly be seen as a rational criterion, as it reflects wishful thinking
[Rapoport, 1998, page 57]. We thus abstain from it in what follows. Our
roughest approximation of the PPP in Section 3, however, will correspond to
this case as well as the more general concept of optimistic superset learning
(OSL) [Hüllermeier, 2014, Rodemann et al., 2022].

##  3 Approximate Bayes Optimal PLS

Since the pseudo posterior predictive (PPP)
$p(\mathcal{D}\cup(x_{i},\hat{y}_{i})\mid\mathcal{D})=p(\hat{y}\mid x,y,x)$
(Theorem 2) is computationally costly to evaluate via Markov Chain Monte Carlo
(MCMC), we aim at approximating it analytically. In light of the general
computational complexity of BPLS involving model refitting, see Section 4,
this appears particularly crucial. We will approximate the joint PPP
directly.666 For i.i.d. data we could focus on the single PPP contributions
$p(y_{i}\mid x_{i},\mathcal{D})$ instead of the joint. Still, we would have to
deal with a possibly intractable integral and end up with similar
computational hustle. We thus opt for approximating the joint directly.
Moreover, considering the joint quantities instead of the distributions
implies no loss of generality, with possible extensions for dependent data in
mind. Our method hence does not need an i.i.d. assumption, which makes it very
versatile.

Due to the aforementioned similarity of the PPP and the marginal likelihood,
we are in the fortunate position of borrowing from some classical marginal
likelihood approximations, see [Llorente et al., 2023]. Especially popular are
approximations based on Laplace’s method as in [Schwarz, 1978]. Our main
motivation, however, is to obtain a Gaussian integral [Gauß, 1877], which we
can then compute explicitly.

###  3.1 Approximation of the PPP

We will start by transferring Laplace’s method to the PPP. Recall that the
predictive posterior of a pseudo-sample $(x_{i},\hat{y}_{i})$ (the PPP) given
data $\mathcal{D}$ is defined as

$\displaystyle p(\mathcal{D}\cup(x_{i},\hat{y}_{i})|\mathcal{D})$
$\displaystyle=\int_{\Theta}p(\mathcal{D}\cup(x_{i},\hat{y}_{i})\mid\theta)p(\theta\mid\mathcal{D})d\theta,$

where Bayes’ theorem gives

$\displaystyle
p(\theta\mid\mathcal{D})=p(\mathcal{D}\mid\theta)\pi(\theta)/p(\mathcal{D}).$

Denoting $\ell_{\mathcal{D}}(\theta)=\log p(\mathcal{D}\mid\theta)$ and
$\tilde{\ell}(\theta)=\ell_{\mathcal{D}\cup(x_{i},\hat{y}_{i})}(\theta)+\ell_{\mathcal{D}}(\theta)$,
we can write the integrand as

$\displaystyle
p(\mathcal{D}\cup(x_{i},\hat{y}_{i})\mid\theta)p(\theta\mid\mathcal{D})=\exp[\tilde{\ell}(\theta)\bigr{)}]\pi(\theta)/p(\mathcal{D}).$

Let $\mathcal{I}(\theta)=-\tilde{\ell}^{\prime\prime}(\theta)/n$ denote the
observed Fisher information matrix. Further denote by
$\tilde{\theta}=\arg\max_{\theta}\tilde{\ell}(\theta)$ the maximizer of
$\tilde{\ell}(\theta)$. It holds $\tilde{\ell}^{\prime}(\theta)=0$ by
definition of $\tilde{\theta}$. A Taylor expansion around $\tilde{\theta}$
thus gives

$\displaystyle\tilde{\ell}(\theta)\approx\tilde{\ell}(\tilde{\theta})-\frac{n}{2}(\theta-\tilde{\theta})^{\prime}\mathcal{I}(\tilde{\theta})(\theta-\tilde{\theta}).$

The integrand decays exponentially in $n\|\theta-\tilde{\theta}\|$, so we can
approximate it locally around $\tilde{\theta}$ by also taking
$\pi(\theta)\leavevmode\nobreak\ \approx\leavevmode\nobreak\
\pi(\tilde{\theta})$ inside the integral with an analogous Taylor series. We
refer to [Miller, 2006, Section 3.7] and [Łapiński, 2019, Theorem 2] for a
rigorous treatment of the remainder terms and regularity conditions.

We can eventually approximate
$p(\mathcal{D}\cup(x_{i},\hat{y}_{i})|\mathcal{D})$ by

$\displaystyle\frac{\exp[\tilde{\ell}(\tilde{\theta})]\pi(\tilde{\theta})}{p(\mathcal{D})}\int_{\Theta}\exp\biggl{[}-\frac{n}{2}(\theta-\tilde{\theta})^{\prime}\mathcal{I}(\tilde{\theta})(\theta-\tilde{\theta})\biggr{]}d\theta,$

The integral on the right is a Gaussian integral. Defining
$\Sigma=[n\mathcal{I}(\tilde{\theta})]^{-1}$ and $\phi_{\Sigma}$ as the
density of the $\mathcal{N}(0,\Sigma)$ distribution, it equals

$\displaystyle(2\pi)^{q/2}|\Sigma|^{1/2}\int_{\Theta}\phi_{\Sigma}(\theta)d\theta=\biggl{(}\frac{2\pi}{n}\biggr{)}^{q/2}|\mathcal{I}(\tilde{\theta})|^{-1/2}.$

Altogether, we have shown that

$\displaystyle p(\mathcal{D}\cup(x_{i},\hat{y}_{i})|\mathcal{D})$
$\displaystyle\approx\biggl{(}\frac{2\pi}{n}\biggr{)}^{q/2}\frac{\exp[\tilde{\ell}(\tilde{\theta})]\pi(\tilde{\theta})}{|\mathcal{I}(\tilde{\theta})|^{1/2}p(\mathcal{D})}.$
(5)

###  3.2 Approximate selection criteria

To find the pseudo-sample $(x_{i},\hat{y}_{i})$ maximizing the PPP, we can
equivalently maximize its logarithm, i.e. maximize

$\begin{split}\frac{q}{2}\,\log\biggl{(}\frac{2\pi}{n}\biggr{)}+\tilde{\ell}(\tilde{\theta})+\log\pi(\tilde{\theta})-\frac{1}{2}\log|\mathcal{I}(\tilde{\theta})|-\log
p(\mathcal{D}).\end{split}$

Dropping all terms that do not depend on $(x_{i},\hat{y}_{i})$ leads to the
selection criterion

$\displaystyle\tilde{\ell}(\tilde{\theta})-\frac{1}{2}\log|\mathcal{I}(\tilde{\theta})|+\log\pi(\tilde{\theta}).$
(6)

The term

$\tilde{\ell}(\theta)=\ell_{\mathcal{D}\cup(x_{i},\hat{y}_{i})}(\theta)+\ell_{\mathcal{D}}(\theta),$

quantifies how well the pseudo-sample $(x_{i},\hat{y}_{i})$ conforms with the
data set $\mathcal{D}$ given a parameter $\theta$, e.g. the optimal (argmax)
parameter $\tilde{\theta}$ in Equation (6). It is curious that samples in
$\mathcal{D}$ contribute twice to $\tilde{\ell}$, but $(x_{i},\hat{y}_{i})$
only once. However, this is irrelevant when comparing two pseudo-samples
$(x_{i},\hat{y}_{i})$ and $(x_{j},\hat{y}_{j})$. To see this, we expand
$\ell_{\mathcal{D}}$ around its maximizer $\hat{\theta}$, so that
$\ell_{\mathcal{D}}(\tilde{\theta})=\ell_{\mathcal{D}}(\hat{\theta})+O(\|\hat{\theta}-\tilde{\theta}\|^{2})$.
Since $\mathcal{D}\cup(x_{i},\hat{y}_{i})$ and $\mathcal{D}$ differ in only
one sample, the difference $\hat{\theta}-\tilde{\theta}$ is of order
$O(n^{-1})$. Thus,

$\tilde{\ell}(\theta)=\ell_{\mathcal{D}\cup(x_{i},\hat{y}_{i})}(\theta)+\ell_{\mathcal{D}}(\hat{\theta})+O(n^{-2}).$

The remainder is negligible compared to the other terms in (6) and
$\ell_{\mathcal{D}}(\hat{\theta})$ does not depend on the pseudo-sample
$(x_{i},\hat{y}_{i})$. This suggests the simplified _informative BPLS
criterion_

$\displaystyle\operatorname{iBPLS}=\ell_{\mathcal{D}\cup(x_{i},\hat{y}_{i})}(\tilde{\theta})-\frac{1}{2}\log|\mathcal{I}(\tilde{\theta})|+\log\pi(\tilde{\theta}).$
(7)

Equivalence of (6) and (7) is verified numerically for small $n$ by
experiments on real-world and simulated data in Supplement F.

The ability to incorporate prior information into the selection is generally a
strength of our criterion. By default, however, we cannot assume that such
information is available. We can instead choose an uninformative prior where
$\pi(\theta)$ is constant with respect to $\theta$. Recall that we assume
$\Theta$ to be compact, which allows us to specify a uniform prior as
uninformative prior. Then (7) simplifies to the _uninformative BPLS criterion_

$\displaystyle\operatorname{uBPLS}=\ell_{\mathcal{D}\cup(x_{i},\hat{y}_{i})}(\tilde{\theta})-\frac{1}{2}\log|\mathcal{I}(\tilde{\theta})|.$
(8)

Our novel PLS criteria provide great intuition.

  * •

The first term is the joint likelihood of the pseudo-sample
$(x_{i},\hat{y}_{i})$ and $\mathcal{D}$ under the optimal parameter
$\tilde{\theta}$. It measures how well the pseudo-sample complies with the
previous model and previously seen data $\mathcal{D}$. It tells the value of
this joint likelihood at its maximum. Loosely speaking, this maximum height of
the likelihood can be seen as a very rough approximation of the area under it,
i.e., the integral with uniform weights.777Technically, we also need that
$\lambda(\Theta)=1$ with $\lambda$ a Lebesgue-measure for this interpretation.

  * •

The second term penalizes high curvature of the pseudo-label likelihood
function $\ell_{\mathcal{D}\cup(x_{i},\hat{y}_{i})}(\theta)$ at its peak
$\tilde{\theta}$, since the Fisher-information is its second derivative. Due
to the negative sign, the criterion prefers pseudo-samples that lead to
flatter maxima of the likelihood. In line with recent insights into sharp and
flat minima of loss surfaces [Dinh et al., 2017, Li et al., 2018,
Andriushchenko and Flammarion, 2022], such a penalty can be expected to
improve generalization. The lower the curvature, the more probability mass
(area under the likelihood) is expected on $\Theta\setminus
B_{\epsilon}(\tilde{\theta})$ with
$B_{\epsilon}=\\{\theta\in\Theta\mid\|\theta-\tilde{\theta}\|<\epsilon\\}$ an
$\epsilon$-ball for fixed $\epsilon>0$ around $\tilde{\theta}$ in the
uninformative case. Intuitively, this corrects the very rough approximation of
the area under the likelihood by the likelihood’s maximal height, see above.

  * •

The third term in the informative BPLS criterion adjusts the selection for our
prior beliefs $\pi$ about $\theta$. Here, the effect of $(x_{i},\hat{y}_{i})$
is only implicit, because it affects the maximizer $\tilde{\theta}$. The more
likely the updated parameter $\tilde{\theta}$ is under $\pi$, the higher the
PPP.

In summary, our approximation of the PPP grows in the absolute value of the
likelihood’s peak, decreases in its curvature at this point, and increases in
the prior likelihood of the updated parameter.

When $n\to\infty$, the criteria iBPLS and uBPLS are dominated by the
likelihood, thus

$\log
p(\mathcal{D}\cup(x_{i},\hat{y}_{i})|\mathcal{D})\overset{n\to\infty}{\propto}\,\ell_{\mathcal{D}\cup(x_{i},\hat{y}_{i})}(\tilde{\theta}).$

This approximation is computationally cheaper to evaluate, as it does not
involve the Fisher-information. However, this comes at the cost of poor
accuracy in case of small $n$. Selection with regard to this rough
approximation of the PPP corresponds to selection with regard to the
likelihood. As pointed out in Section 2.2, this corresponds to the overly
optimistic max-max-criterion.

Figure 1: Results from 8 classification tasks based on real-world data [Dua
and Graff, 2017] in descending difficulty (measured by supervised test
accuracy), where $p$ denotes the number of features here and the share of
unlabeled data is 0.8. Accuracy averaged over 100 repetitions.

##  4 EXPERIMENTS

Algorithmic Procedure: For all predicted pseudo-labels, we refit the model on
$\mathcal{D}\leavevmode\nobreak\ \cup\leavevmode\nobreak\ (x_{i},\hat{y}_{i})$
and evaluate its PPP by means of the derived approximations iBPLS and uBPLS to
select one instance to be added to the training data. Detailed pseudo code for
BPLS can be found in Supplement A. The computational complexity depends on the
evaluation of the PPP. With $\lvert\mathcal{U}\rvert=m$ unlabeled data points
and no stopping criterion, $m+(m-1)+\dots+1=\frac{m^{2}+m}{2}$ PPPs have to be
evaluated (that is, approximated). Hence, BPLS’ complexity depends on the
model’s complexity and the amount of unlabeled data.

###### Hypothesis 1

(a) PPP with uninformative prior outperforms traditional PLS on data prone to
initial overfitting (i.e., with high ratio of features to data $\frac{q}{n}$
and poor initial generalization). (b) For low $\frac{q}{n}$ and high initial
generalization, BPLS is outperformed by traditional PLS.

###### Hypothesis 2

(a) Among all PLS methods, the pseudo-label likelihood (max-max-action)
reinforces the initial model fit the most and (b) hardly improves
generalization.

###### Hypothesis 3

PPP with informative prior outperforms traditional PLS methods universally.

Experimental Setup: We formulate three hypotheses beforehand. Hypothesis 1
corresponds to the main motivation behind BPLS; its second part is a logical
consequence thereof: If we are sceptical towards the initial model in case it
generalizes well, we expect to select pseudo-labels in a worse way than when
trusting the initial model. Hypothesis 2 is based on the decision-theoretic
insights regarding PLS by the likelihood, see section 2.2: It embodies an
optimistic reliance on the initial model and is thus expected to pick data
that fits best into that model. We further expect (Hypothesis 3) BPLS to
unambiguously outperform non-Bayesian selection methods in case the prior
provides actual information about the data generating process – the latter is
simply not available for non-Bayesian PLS. We benchmark semi-supervised
(parametric) generalized linear models (GLMs) and (non-parametric) generalized
additive models (GAMs) [Hastie and Tibshirani, 1987, Hastie, 2017] with PPP
and pseudo-label likelihood against two common selection criteria (probability
score and predictive variance) [Triguero et al., 2015] as well as a supervised
baseline. For the latter, we abstain from self-training and only use the
labeled data for training. Experiments are run on simulated binomially
distributed data as well as on eight data sets for binary classification from
the UCI repository [Dua and Graff, 2017]. The binomially distributed data was
simulated through a linear predictor consisting of normally distributed
features. Details on the simulations as well as on the data sets can be found
in Supplement C and Supplement H. The share of unlabeled data was set to $0.8$
and $0.9$. PLS methods were compared w.r.t. to (inductive) accuracy of
prediction on unseen test. All data sets were found to be fairly balanced
except for the EEG data (minority share: $0.29$).

Results: Figures 1 and 2 as well as Table 1 summarize the results in the
uninformative case (grey figures) for real-world and simulated data,
respectively. Oracle stopping in Table 1 refers to comparing PLS methods with
regard to their overall best accuracy as opposed to final comparisons after
the whole data set was labeled. Figure 2 sheds further light on results for
simulated data, while Figure 3 displays results from benchmarking BPLS to
classical PLS methods in the informative case (black figures). Detailed
figures displaying results from all experiments can be found in the
supplementary material.

Figure 2: Results from simulated data. Accuracy averaged over 100 repetitions.
Legend: see Figure 1.

Interpretation: At first sight, comparing the accuracy gains in Figure 1 on
different data sets (in order of ascending baseline performance) clearly
supports Hypothesis 1: For harder tasks like EEG or sonar with relatively high
ratio of features to data $\frac{q}{n}$, Bayesian PPP outperforms traditional
PLS, whilst being dominated by the probability score in case of easier tasks
like banknote or breast cancer. For data sets with intermediate difficulty
(mushrooms and ionosphere), PPP and other PLS methods compete head-to-head.
The results on abalone data underpin a general fact in SSL (see section 1):
Successful self-training requires at least some baseline supervised
performance. Results on simulated data (Table 1) further support the role of
$\frac{q}{n}$ in Hypothesis 1. Their visualization (Figure 2) nicely
illustrates the inner working of selection by PPP: By not trusting the initial
model, PPP affects the model’s test accuracy the most. While $n=400$ leaves
some room for improvement through mitigating the overfitting by pseudo-labeled
data, PPP leads to a noisy performance in case of $n=100$ close to $p$. Here,
even the final model still overfits. These promising results should not hide
an inconsistency: The fact that PPP is superior on the cars task but not on
the ionosphere task contradicts Hypothesis 1 (a), since cars is harder than
ionosphere, while having almost identical $\frac{q}{n}$. We find Hypothesis 2
to be partially supported by the results. While 2 (a) holds for both the
majority of simulated (see supplementary material) and real-world data
(likelihood generally the closest to supervised performance), 2 (b) is
challenged by considerable generalization performance gain on ionosphere and
breast cancer data.

Figure 3 clearly supports Hypothesis 3: When using informative priors based on
the true data-generating process, BPLS clearly outperforms traditional PLS
methods. Results in Supplement D.3 further back this finding. This comes at no
big surprise, since non-Bayesian PLS simply lack ways to incorporate such
prior knowledge. From this perspective, the uninformative case (Hypothesis 1)
corresponds to raising the bar and clearly is the theoretically more
interesting benchmarking setup. However, many practical applications of SSL
entail a myriad of pre-existing knowledge, e.g., radio spectrum identification
[Camelo et al., 2019]. For practical purposes, thus, the informative situation
might even be more relevant.

Table 1: Best performing PLS method (uninformative) on simulated data n | q | ORACLE STOPPING | FINAL
---|---|---|---
60 | 60 | PPP | PPP
100 | 60 | PPP | Supervised Learning
400 | 60 | PPP | PPP
1000 | 60 | Probability Score | Probability Score
Figure 3: Results of PPP with informative priors and non-parametric GAMs on
simulated data with different shares of unlabeled data. Accuracy averaged over
100 repetitions.

##  5 RELATED WORK

Robust PLS: Robustness of PLS is a widely discussed issue in the self-training
literature. [Aminian et al., 2022] propose information-theoretic PLS robust
towards covariate shift. [Lienen and Hüllermeier, 2021] label instances in the
form of sets of probability distributions (credal sets), weakening the
reliance on a single distribution. [Vandewalle et al., 2013] aim at robustness
to modeling assumptions by allowing model selection through the deviance
information criterion during semi-supervised learning. [Rizve et al., 2020]
propose uncertainty-aware pseudo-label selection which proves to compete with
state-of-the-art SSL based on consistency regulation. The idea is to select
pseudo-labeled instances whose probability score and predictive uncertainty
are above (tunable) thresholds. The latter is operationalized by the
prediction’s variance, and thus, unlike BPLS, fails to decompose approximation
and modeling uncertainty, see Section 1. Both predictive variance and
probability score serve as benchmarks in Section 4.

Bayesian Self-Training: There is a broad body of research on deploying
Bayesian predictions in SSL and particularly in self-training [Gordon and
Hernández-Lobato, 2020, Ng et al., 2018, Cai et al., 2011, Adams and
Ghahramani, 2009]. The same holds for explicit likelihood-based inference,
such as weighted likelihood [Sokolovska et al., 2008], conditional likelihood
[Grandvalet and Bengio, 2004], and joint mixture likelihood [Amini and
Gallinari, 2002]. Most of them use Bayesian models for predicting pseudo-
labels. In contrast, we prove that the argmax of the PPP is the Bayes optimal
selection of pseudo-labels given any predictive model. Regarding Bayesian or
likelihood-based selection of pseudo-labels, there exists only little
(Bayesian) or hardly any (likelihood-based) work. [Li et al., 2020] quantify
the uncertainties of pseudo-labels by mixtures of predictive distributions of
a neural net, applying MC dropout. This could be seen as an expensive MC-based
approximation of the PPP. Very recently, [Patel et al., 2023] proposed PLS
with regard to (a sampling-based approximation of) the entropy of the pseudo-
labels’ posterior predictive distribution. The entropy is considered a measure
of total uncertainty (aleatoric and epistemic) and often considered as
regularization for PLS, see [Saporta et al., 2020, Liu et al., 2021] for
instance. Abstaining from the entropy – as we do – effectively means not
considering the aleatoric uncertainty. While including aleatoric uncertainty
(e.g. measurement noise) generally makes sense, we consider it of minor
importance in the concrete problem of initial overfitting, where we aim at
disentangling epistemic uncertainty with regard to data and parameters: We
want to choose pseudo-labels that are likely given the observed labeled data
but not necessarily likely given the estimated parameters of the (over-)fitted
model.

##  6 DISCUSSION

Extensions: We briefly discuss two venues for future work. The first extension
loosens the restriction on one particular model class by performing model
selection and PLS simultaneously. The idea would be to select these instances
that can be best explained by the simplest learner (i.e., the one with least
parameters). Further recall that both the framework of BPLS and our
approximation of the PPP do not require data to be i.i.d distributed. Applying
BPLS on dependent observations, such as in auto-correlated data like time
series, is thus another promising line of further research.

Limitations: BPLS’ strength of being applicable to any learner can imply high
computational costs in case of expensive-to-train models such as neural nets,
because PPP approximations require refitting the model. Additionally, it might
be difficult for practitioners to assess the risk of overfitting to the
initial data set beforehand and opt for BPLS in response. Given the fact that
BPLS is outperformed by traditional PLS in cases with no overfitting, this
might be considered a drawback for practical application. However, Section 4
demonstrated that $\frac{q}{n}$ and the baseline supervised performance (both
easily accessible) provide sound proxies for initial overfitting scenarios
that can induce a confirmation bias in PLS. These proxies can (alongside
cross-validation) help practitioners to identify such scenarios.

Conclusion: BPLS renders self-training more robust with respect to the initial
model. This improves final performance if the latter overfits and harms it if
not. Identifying overfitting scenarios is thus crucial for BPLS’ usage. What
is more, BPLS allows incorporating prior knowledge, with the help of which
substantial performance gains can be achieved. Besides, our insights from
formalizing PLS as a decision problem clear the way for promising future work
exploiting rich literature on Bayesian decision theory. Ultimately, we
conclude that a Bayesian view can add great value not only to predicting but
also to selecting data for self-training.

#### Acknowledgements

Thomas Augustin gratefully acknowledges support by the Federal Statistical
Office of Germany within the co-operation project "Machine Learning in
Official Statistics".

## References

  * [Abdar et al., 2021] Abdar, M., Pourpanah, F., Hussain, S., Rezazadegan, D., Liu, L., Ghavamzadeh, M., Fieguth, P., Cao, X., Khosravi, A., Acharya, U. R., et al. (2021).  A review of uncertainty quantification in deep learning: Techniques, applications and challenges.  Information Fusion, 76:243–297.
  * [Adams and Ghahramani, 2009] Adams, R. P. and Ghahramani, Z. (2009).  Archipelago: nonparametric Bayesian semi-supervised learning.  In 26th International Conference on Machine Learning, pages 1–8.
  * [Amini and Gallinari, 2002] Amini, M.-R. and Gallinari, P. (2002).  Semi-supervised logistic regression.  In 15th European Conference on Artificial Intelligence (ECAI 2002), pages 390–394.
  * [Aminian et al., 2022] Aminian, G., Abroshan, M., Khalili, M. M., Toni, L., and Rodrigues, M. (2022).  An information-theoretical approach to semi-supervised learning under covariate-shift.  In International Conference on Artificial Intelligence and Statistics, pages 7433–7449. PMLR.
  * [Andriushchenko and Flammarion, 2022] Andriushchenko, M. and Flammarion, N. (2022).  Towards understanding sharpness-aware minimization.  In International Conference on Machine Learning, pages 639–668.
  * [Arazo et al., 2020] Arazo, E., Ortego, D., Albert, P., O’Connor, N. E., and McGuinness, K. (2020).  Pseudo-labeling and confirmation bias in deep semi-supervised learning.  In 2020 International Joint Conference on Neural Networks, pages 1–8. IEEE.
  * [Barber, 2012] Barber, D. (2012).  Bayesian reasoning and machine learning.  Cambridge University Press.
  * [Berger, 1985] Berger, J. O. (1985).  Statistical decision theory and Bayesian analysis. Springer, Berlin., 2nd edition.
  * [Cai et al., 2011] Cai, R., Zhang, Z., and Hao, Z. (2011).  BASSUM: A Bayesian semi-supervised method for classification feature selection.  Pattern Recognition, 44(4):811–820.
  * [Camelo et al., 2019] Camelo, M., Shahid, A., Fontaine, J., de Figueiredo, F. A. P., De Poorter, E., Moerman, I., and Latre, S. (2019).  A semi-supervised learning approach towards automatic wireless technology recognition.  In 2019 IEEE International Symposium on Dynamic Spectrum Access Networks (DySPAN), pages 1–10.
  * [Cattaneo, 2007] Cattaneo, M. E. (2007).  Statistical decisions based directly on the likelihood function.  PhD thesis, ETH Zurich.
  * [Chapelle et al., 2006] Chapelle, O., Schölkopf, B., and Zien, A. (2006).  Semi-supervised learning. Adaptive computation and machine learning series.  MIT Press.
  * [Dinh et al., 2017] Dinh, L., Pascanu, R., Bengio, S., and Bengio, Y. (2017).  Sharp minima can generalize for deep nets.  In International Conference on Machine Learning, pages 1019–1028. PMLR.
  * [Dua and Graff, 2017] Dua, D. and Graff, C. (2017).  UCI machine learning repository.  http://archive.ics.uci.edu/ml.
  * [Gauß, 1877] Gauß, C. F. (1877).  Theoria motus corporum coelestium in sectionibus conicis solem ambientium, volume 7.  FA Perthes.
  * [Gordon and Hernández-Lobato, 2020] Gordon, J. and Hernández-Lobato, J. M. (2020).  Combining deep generative and discriminative models for Bayesian semi-supervised learning.  Pattern Recognition, 100:107156.
  * [Grandvalet and Bengio, 2004] Grandvalet, Y. and Bengio, Y. (2004).  Semi-supervised learning by entropy minimization.  In Advances in Neural Information Processing Systems, volume 17.
  * [Hastie, 2017] Hastie, T. (2017).  Generalized additive models.  In Chambers, J. M. and Hastie, T., editors, Statistical models in S, pages 249–307. Routledge.
  * [Hastie and Tibshirani, 1987] Hastie, T. and Tibshirani, R. (1987).  Generalized additive models: some applications.  Journal of the American Statistical Association, 82(398):371–386.
  * [Hüllermeier, 2014] Hüllermeier, E. (2014).  Learning from imprecise and fuzzy observations: Data disambiguation through generalized loss minimization.  International Journal of Approximate Reasoning, 55:1519–1534.
  * [Hüllermeier and Waegeman, 2021] Hüllermeier, E. and Waegeman, W. (2021).  Aleatoric and epistemic uncertainty in machine learning: An introduction to concepts and methods.  Machine Learning, 110(3):457–506.
  * [Łapiński, 2019] Łapiński, T. M. (2019).  Multivariate Laplace’s approximation with estimated error and application to limit theorems.  Journal of Approximation Theory, 248:105305.
  * [Lee et al., 2013] Lee, D.-H. et al. (2013).  Pseudo-label: The simple and efficient semi-supervised learning method for deep neural networks.  In Workshop on challenges in representation learning, International Conference on Machine Learning, volume 3, page 896.
  * [Li et al., 2018] Li, H., Xu, Z., Taylor, G., Studer, C., and Goldstein, T. (2018).  Visualizing the loss landscape of neural nets.  In Advances in Neural Information Processing Systems, volume 31.
  * [Li et al., 2020] Li, S., Wei, Z., Zhang, J., and Xiao, L. (2020).  Pseudo-label selection for deep semi-supervised learning.  In 2020 IEEE International Conference on Progress in Informatics and Computing (PIC), pages 1–5. IEEE.
  * [Lienen and Hüllermeier, 2021] Lienen, J. and Hüllermeier, E. (2021).  Credal self-supervised learning.  In Advances in Neural Information Processing Systems, volume 34, pages 14370–14382.
  * [Liu et al., 2021] Liu, H., Wang, J., and Long, M. (2021).  Cycle self-training for domain adaptation.  In Advances in Neural Information Processing Systems, volume 34, pages 22968–22981.
  * [Llorente et al., 2023] Llorente, F., Martino, L., Delgado, D., and Lopez-Santiago, J. (2023).  Marginal likelihood computation for model selection and hypothesis testing: an extensive review.  SIAM Review, 65(1):3–58.
  * [Lotfi et al., 2022] Lotfi, S., Izmailov, P., Benton, G., Goldblum, M., and Wilson, A. G. (2022).  Bayesian model selection, the marginal likelihood, and generalization.  In International Conference on Machine Learning, pages 14223–14247.
  * [Malinin and Gales, 2018] Malinin, A. and Gales, M. (2018).  Predictive uncertainty estimation via prior networks.  In Advances in Neural Information Processing Systems, volume 31.
  * [McClosky et al., 2006] McClosky, D., Charniak, E., and Johnson, M. (2006).  Effective self-training for parsing.  In Proceedings of the Human Language Technology Conference of the NAACL, Main Conference, pages 152–159.
  * [Miller, 2006] Miller, P. D. (2006).  Applied asymptotic analysis, volume 75.  American Mathematical Soc.
  * [Murphy, 2012] Murphy, K. P. (2012).  Machine learning: a probabilistic perspective.  MIT press.
  * [Ng et al., 2018] Ng, Y. C., Colombo, N., and Silva, R. (2018).  Bayesian semi-supervised learning with graph gaussian processes.  In Advances in Neural Information Processing Systems, volume 31.
  * [Patel et al., 2023] Patel, G., Allebach, J. P., and Qiu, Q. (2023).  Seq-ups: Sequential uncertainty-aware pseudo-label selection for semi-supervised text recognition.  In Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 6180–6190.
  * [Rapoport, 1998] Rapoport, A. (1998).  Decision theory and decision behaviour.  Springer.
  * [Rizve et al., 2020] Rizve, M. N., Duarte, K., Rawat, Y. S., and Shah, M. (2020).  In defense of pseudo-labeling: An uncertainty-aware pseudo-label selection framework for semi-supervised learning.  In International Conference on Learning Representations, 2020.
  * [Rodemann et al., 2022] Rodemann, J., Kreiss, D., Hüllermeier, E., and Augustin, T. (2022).  Levelwise data disambiguation by cautious superset classification.  In International Conference on Scalable Uncertainty Management, pages 263–276. Springer.
  * [Saporta et al., 2020] Saporta, A., Vu, T.-H., Cord, M., and Pérez, P. (2020).  ESL: Entropy-guided self-supervised learning for domain adaptation in semantic segmentation.  arXiv preprint arXiv:2006.08658.
  * [Schwarz, 1978] Schwarz, G. (1978).  Estimating the dimension of a model.  The Annals of Statistics, 6(2):461–464.
  * [Shi et al., 2018] Shi, W., Gong, Y., Ding, C., Tao, Z. M., and Zheng, N. (2018).  Transductive semi-supervised deep learning using min-max features.  In European Conference on Computer Vision, pages 299–315.
  * [Sohn et al., 2020] Sohn, K., Berthelot, D., Carlini, N., Zhang, Z., Zhang, H., Raffel, C. A., Cubuk, E. D., Kurakin, A., and Li, C.-L. (2020).  Fixmatch: Simplifying semi-supervised learning with consistency and confidence.  In Advances in Neural Information Processing Systems, volume 33.
  * [Sokolovska et al., 2008] Sokolovska, N., Cappé, O., and Yvon, F. (2008).  The asymptotics of semi-supervised learning in discriminative probabilistic models.  In 25th International Conference on Machine learning, pages 984–991.
  * [Triguero et al., 2015] Triguero, I., García, S., and Herrera, F. (2015).  Self-labeled techniques for semi-supervised learning: taxonomy, software and empirical study.  Knowledge and Information systems, 42(2):245–284.
  * [Van Engelen and Hoos, 2020] Van Engelen, J. E. and Hoos, H. H. (2020).  A survey on semi-supervised learning.  Machine Learning, 109(2):373–440.
  * [Vandewalle et al., 2013] Vandewalle, V., Biernacki, C., Celeux, G., and Govaert, G. (2013).  A predictive deviance criterion for selecting a generative model in semi-supervised classification.  Computational Statistics & Data Analysis, 64:220–236.


# FPGA Deployment of LFADS for Real-time Neuroscience Experiments

Xiaohan Liu<EMAIL_ADDRESS>University of WashingtonSeattleWA 98195USA ,
ChiJui Chen<EMAIL_ADDRESS>National Yang Ming Chiao Tung
UniversityHsinchuTaiwan ,  YanLun Huang<EMAIL_ADDRESS>National Yang
Ming Chiao Tung UniversityHsinchuTaiwan ,  LingChi Yang
<EMAIL_ADDRESS>National Yang Ming Chiao Tung
UniversityHsinchuTaiwan ,  Elham E Khoda<EMAIL_ADDRESS>University of
WashingtonSeattleWA 98195USA ,  Yihui Chen<EMAIL_ADDRESS>University of
WashingtonSeattleUSA ,  Scott Hauck<EMAIL_ADDRESS>University of
WashingtonSeattleWA 98195USA ,  Shih-Chieh Hsu<EMAIL_ADDRESS>University of
WashingtonSeattleWA 98195USA and  Bo-Cheng Lai<EMAIL_ADDRESS>National
Yang Ming Chiao Tung UniversityHsinchuTaiwan

(2023)

###### Abstract.

Large-scale recordings of neural activity are providing new opportunities to
study neural population dynamics. A powerful method for analyzing such high-
dimensional measurements is to deploy an algorithm to learn the low-
dimensional latent dynamics. LFADS (Latent Factor Analysis via Dynamical
Systems) is a deep learning method for inferring latent dynamics from high-
dimensional neural spiking data recorded simultaneously in single trials. This
method has shown a remarkable performance in modeling complex brain signals
with an average inference latency in milliseconds. As our capacity of
simultaneously recording many neurons is increasing exponentially, it is
becoming crucial to build capacity for deploying low-latency inference of the
computing algorithms. To improve the real-time processing ability of LFADS, we
introduce an efficient implementation of the LFADS models onto Field
Programmable Gate Arrays (FPGA). Our implementation shows an inference latency
of 41.97 $\mu$s for processing the data in a single trial on a Xilinx U55C.

neuroscience, latent variable model, recurrent neural network, behavioural
neuroscience, hls4ml, FPGA

††copyright: acmcopyright††journalyear: 2023††conference: ; November 2, 2023;
San Francisco, CA

##  1\. Introduction

Figure 1. LFADS architecture used for this study.

Over the past decade, the ability to record large-scale neural activity has
improved dramatically, providing us with new opportunities to study neural
population dynamics. A powerful strategy for analyzing such high-dimensional
measurements is to learn the low-dimensional latent dynamics that explain much
of the variance in the measurements. Analysis of neural signals also widely
relies on neural networks (Yang and Wang, 2020). LFADS (Sussillo et al., 2016;
Pandarinath et al., 2018) model is a cutting-edge neural network architecture
in the field of neural modeling and computational neuroscience. LFADS
leverages the capacity of recurrent neural networks (RNNs) for uncovering
hidden patterns and combines them with variational inference techniques to
infer the latent factors that drive the observed data. Although LFADS
demonstrates satisfactory performance in modeling brain signals, processing
the extensive neural recordings in real time still poses significant
computational challenges. Large-scale neural recordings are commonly involved
in novel neuroscience experiments (Cunningham and Yu, 2014). To process the
large-scale neural recordings in real-time, Field Programmable Gate Arrays
(FPGAs) are used to accelerate the inference.

FPGAs enable customized data processing logic and have gained widespread
adoption for achieving highly parallel dataflow processing with minimal
latency. Thus, low-latency, low-power, and high-throughput model inference can
be achieved on FPGAs, which makes large-scale real-time neural experiments
possible. For instance, Low-latency LFADS can be employed to create real-time
closed-loop experiments that decode neural activity to manipulate external
devices. Furthermore, high-throughput machine learning (ML) in neural-related
experiments improves the ability to examine extensive neural recordings. This
capability is increasingly crucial in modern neuroscience experiments, as the
integration of large-scale neural recordings has become a commonplace and
essential component.

In this paper, we present an efficient implementation of the LFADS model in
High-Level Synthesis (HLS) for the hls4ml package (FastML Team, 2021). This
HLS implementation will opens the door for wider low-latency and high-
throughput applications of the LFADS models for neural sequential data.
Several optimizations for input/output (IO) ports and data access mechanisms
are done while implementing the HLS representation of the LFADS model.

##  2\. Core Concepts

As LFADS is one of the promising models for modeling complex brain activities
by inferring smooth dynamics based on the collected neural spiking data, it is
considered for this study. Typically, data collection by the NeuroPort System
(Blackrock Microsystems) and neural control and cueing tasks on the
Simulink/xPC real-time platform (MathWorks) can be performed within 10 ms
(Pandarinath et al., 2018). The data collection and experimental control can
be done much faster by using hardware like FPGAs. The focus of the work is to
demonstrate the ability to run low-latency inference of LFADS-like models on
an FPGA. A modified LFADS model is used for this study as described in Sec.
2.1.

###  2.1. Model Description

LFADS is a sequential model based on a variational autoencoder (VAE). The
encoder has a bidirectional GRU layer, which takes the single-trial neural
spikes as input and converts them into a low-dimensional latent space
representation. This representation acts as the initial state for the decoder
GRU layer, which tries to regenerate the input data. LFADS assumes the
observed spikes are samples from a Poisson process with firing rates
($r_{t}$). So, instead of estimating the neural spikes, the decoder learns the
firing rates ($r_{t}$) as a function of time. The decoder (or generator) GRU
is expected to learn the underlying dynamics. So, the decoder is trained to
infer a reduced set of latent dynamic factors, $f_{t}$. The firing rates can
be constructed from the latent factors.

In this work, we have studied an autoencoder (AE)-based model to reduce the
complexity of FPGA deployment. Furthermore, our studies show that removing the
Gaussian-sampling has almost no impact on the model performance. The original
LFADS model uses a controller network to predict the dynamics of the system in
the presence of an external input. To simplify the model, we have studied an
LFADS, which can only predict the autonomous dynamics. Figure 1 presents the
overview of the LFADS architecture used in this work. The model structure is
adopted from (Hurwitz et al., 2021). In our model, a bidirectional GRU encoder
with 64 units for both forward and backward layers, compresses the input
spikes into a latent vector of dimension 64. Then the latent variables are
used for initializing the unidirectional GRU decoder, which has 64 GRU units.
The decoder produces a set of low-dimensional temporal factors of dimension 4.
These four latent dynamic factors are then passed through a fully connected
layer with 70 units to produce the log firing rates, $\log(r_{t})$,
corresponding to the input spike. The model is trained with Keras (et al.,
2015) and TensorFlow (Developers, 2023).

###  2.2. Dataset and Evaluation Metrics

In this paper, we have used the data collected from a study of monkey reaching
tasks (Perich et al., 2018). The input dataset contains a total of 170 trials.
Out of which, 136 trials (80%) are used for training, while 17 trials are used
for each, validation and test. Each trial consists of 70 recording channels,
with 73 discrete time steps per channel, resulting a shape of (136,73,70) for
the training dataset. A detailed training and evaluation methodologies are
presented in Ref. (Hurwitz et al., 2021).

Two metrics are used to evaluate the performance of LFADS model: negative
Poisson log-likelihood (NPLL) and the coefficient of determination (R2). As
LFADS estimates the firing rates, assuming spiking variability follows a
Poisson distribution, it uses NPLL as a loss function in training. It should
be noted that the main goal of this work is to deploy LFADS onto an FPGA, not
to optimize its performance. Figure 2 shows the NPLL comparison of a VAE-based
LFADS model with the AE-based model used in this work. The R2 score is
calculated by fitting the reconstructed temporal factors $f_{t}$ to the
measured behavioral data (hand position). The closer the score is to 1, the
better the factors align with the behavioral data, indicating a stronger
correlation between the model’s predictions and the observed behavior.

Figure 2. The figure shows the NPLL values of the AE-based LFADS (red) and a
VAE-based LFADS (blue) with a similar architecture.

##  3\. Implementation

One of the main focuses of this work is to implement the LFADS model in HLS
and integrate it into hls4ml for rapid development on FPGA. Currently, almost
all components of the LFADS model are available in hls4ml, except for the
Bidirectional wrapper for the GRU and the quantized GRU (QGRU).

###  3.1. HLS implementation: Keras Model

The HLS bidirectional wrapper is used for the Bidirectional layer of the Keras
LFADS model. To implement the Bidirectional wrapper, a reverse function is
applied to invert the order of the input sequence. During inference, the
original and reversed inputs pass through corresponding GRU layers and their
outputs are concatenated.

###  3.2. QKeras Model

We have also studied quantization-aware training (QAT) using the QKeras
package (et al., 2019) . We have built our QKeras model using the QGRU,
QBidirectional, QDense layers and different quantized form of sigmoid and tanh
activation functions. Since these two activation functions are non-linear,
outputs of the original and the quantized function could be different. So,
this has been properly optimized in our implementation. For QBidirectional, we
use two QGRUs, one for the forward and one for the backward input sequences,
and concatenate the two outputs to form the QBidirectional output. It is
important to note the precision of the inputs also needs to be adjusted by
using a quantizer to make different operations consistent.

Most of the hyperparameters used for QAT are similar to the floating-point
models. For the QKeras model, a step learning rate schedule is used, with an
initial learning rate of 0.01, patience of 10 epochs, decay factor of 0.5, and
a minimum learning rate of $10^{-5}$. The SGD optimizer is used instead of
Adam when the total bits are less than 6 bits for more stable training. After
through scanning, we have chosen to use 3 integer bits for all activations and
1 integer bit (only sign bit) for weights and biases, since they are
initialized with the Lecun uniform initializer (Klambauer et al., 2017).

Figure 3. The structure of 64-units quantized GRU cell

###  3.3. HLS implementation: QKeras Model

Figure 3 displays the architecture of a 64-unit quantized GRU cell. There are
four weight quantizers, two activation quantizers (quantized sigmoid, and
quantized tanh) and one state quantizer. During the hls4ml compilation, the
precisions for all the following calculations are automatically determined
based on those quantizers. For example, the colored arrows in Figure 3 show
the required bit-width for the addition before the sigmoid.

For the choice of quantized activation, we use hard quantized activation.
Figure 4 shows the curves of the real activations and the quantized (hard)
activations. In HLS, it requires a look-up table to achieve the non-linearity
in quantized activations. However, as the bit-width increases, the size of the
look-up table exponentially grows to maintain accuracy. In contrast, the hard
sigmoid calculation can be supported with simple wiring in hardware, without
requiring a multiplier or look-up table. This statement is also true for hard
tanh as it can be derived from hard sigmoid.

(a)

(b)

Figure 4. The plots of 4-bits quantized activations, quantized hard
activations, and real activations (a) quantized sigmoid and quantized hard
sigmoid ranged from 0 to 0.9375. (b) quantized tanh and quantized hard tanh
ranged from -1 to 0.875.

The current version of hls4ml has two options for the datatype used in the
dataflow: IO-parallel and IO-stream. To maximize the bandwidth, fully-
partitioned arrays were utilized in IO-parallel, which can be resource-
intensive. As a result, we use streams since they are synthesized as more
resource-efficient FIFOs. Furthermore, we have modified the hls4ml compiler to
use array-of streams datatype, offering more flexibility in terms of bandwidth
compared to the default packed struct datatype. The GRU layers only support
IO-parallel, so we have applied a simple wrapper layer to the GRU to convert
the stream type to an array type for IO. Similarly, for the bidirectional
layer, we have employed the same approach to handle IO.

##  4\. Results

The LFADS model described in Sec. 2.1 is translated into an HLS model by
utilizing the hls4ml tool. The implementation is tested using Vivado HLS
2020.1 with a Xilinx Alveo U55C FPGA as the target. In this work, we have
carefully optimized quantization and resource utilization on the FPGA.

###  4.1. Quantization Results

The quantization process is to reduce the precision of the model parameters as
well as the inputs. Typically the model parameters, such as weights and
biases, are stored as 32-bit floating-point numbers in the Keras models. As
the floating-point number takes up a lot of resources on the FPGA, it is
preferred to use fixed-point representations of the model parameters. The
fixed-point format is referred to Vivado HLS ap_fixed type (Xilinx, 2023a),
where the sign bit is included in the integer bits. In this is paper, the
ap_fixed numbers are written in <total bits, integer bits> format. hls4ml
converts the input model parameters into ap_fixed while implementing the HLS
model representation. This process is called post-training quantization (PTQ).
In our study, we varied the fractional bits between 2 and 16 while keeping the
precision of the integer part fixed to 4, 6, and 8. The NPLL and $R^{2}$
scores are calculated for different quantization settings and the values are
compared with that of the floating-point model as shown in Figure 5. Different
colors in the figure represent integer bits of 4 (orange), 6 (green) or 8
(red). We observe that at least 6 integer bits and 10 fractional bits, <16,6>,
are needed to achieve a similar performance as the floating-point model.

(a) PTQ NPLL

(b) PTQ R2

Figure 5. Shows (a) NPLL and (b) R2 score as a function of fractional bits.
The blue line in each figure represents the floating-point, whereas the lines
correspond to inter bits of 4 (orange), 6 (green), or 8 (red).

(a) QAT NPLL

(b) QAT R2

Figure 6.  Shows (a) NPLL and (b) R2 score as a function of total bits for the
quantized-aware trained (Q) models using. Results from the baseline floating-
point (FP) model are shown in the right subplots.

The performance of the QAT model is estimated for different total bit widths
between 4 and 16, as shown in Figure 6. We see a noticeable degradation in
performance below total width of 10 bits. A similar observation can be seen in
the behavior reconstruction shown in Figure 7. The figure shows the hand
movement trajectories in the 2D $x-y$ plan. The movements performed by the
subject in the same direction are grouped together and denoted by the same
color. The dotted lines represent the target movement direction, where the
reconstructed movement trajectories are shown in solid lines. The QAT models
trained with less than 10 total widths show significant degradation in
performance. The performance of the 12-bit QAT model is much closer to the
floating-point model. Based on the results, we have selected the QAT model
with 10 bits for subsequent steps.

(a) 4 bits

(b) 8 bits

(c) 12 bits

(d) floating point

Figure 7. The visualizations compare the reconstructed and target behavior.
The reconstructed hand movement trajectories are predicted by the QAT models
trained with (a) 4, (b) 8, and (c) 12 total bits and (d) floating-point model
with the target behavior. The movements in different directions are grouped
and shown in different colors. The solid lines denote the reconstructed
movements, while the dotted lines represent the target movements.

###  4.2. Resource Utilization

Figure 8. Resource utilization for different total bits with post-training
quantization. The target FPGA is Xilinx Alveo U250.

We have used Alveo U250, which is currently the largest board in the Xilinx
Alveo FPGA series, as the target device to ensure the the available resources
are sufficient at a high bit width. The LFADS model is synthesized using
Vivado HLS with different quantization precision to obtain a resource
estimation. For each case, we estimate the utilization of different FPGA
resources like memory (BRAM), digital signal processing units (DSPs), flip-
flops (FFs), and lookup tables (LUTs). Figure 8 shows the resource consumption
of the LFADs model with post-training quantization in relation to the total
available resources. The four curves, representing different resources,
exhibit a downward trend as the total bit width decreases. The limitation of
FPGA inference in high bitwidth is DSPs. The DSP consumption decreases
significantly from 32 bits to 22 bits, followed by a stable trend until 11
bits. This is due to the maximum input size for multiplication in DSP48E2 is
27 × 18 (Xilinx, 2023b). If an input exceeds this limit, two DSPs will be
applied to perform the multiplication. For bit widths below 12 bits, the DSP
utilization decreases to nearly zero percent since the multiplication is
carried out by LUTs, resulting in an increase in LUT consumption from the
utilization of 12 bits to 10 bits.

###  4.3. FPGA Latency

To quantify latency, we selected the optimal PTQ model with precision
ap_fixed¡16,6¿, comprising six integer bits and ten fractional bits. Given the
constraints of our experimental setup, we opted for a more compact FPGA board,
the Xilinx Alveo U55C, as it suits our requirements. Despite its smaller size
compared to the Alveo U250, the U55C board comfortably accommodates the best
PTQ model, with minimal anticipated alterations in latency results. Table 1
shows the configuration used for this test. The actual running frequency on
the U55C is 200 MHz, and the estimated latency is 41.97 $\mu$s.

Table 1. Deployment Configuration Target FPGA | Xilinx Alveo U55C
---|---
Precision | ap_fixed¡16,6¿
Frequency | 200 MHz
Latency | 41.97$\mu$s

##  5\. Summary and Outlook

We develop an automated design workflow based on hls4ml to deploy LFADS onto
an FPGA to accelerate the model for potential large-scale real-time
neuroscience experiments. For efficient implementation of the model, we
studied both PTQ and QAT. Using QAT it is possible to reduce total number of
bits required to 10 bits, compared to the 16 bits required for PTQ, while
incurring negligible loss compared to the floating-point model. For the on-
board evaluation, the 16-bit model demonstrates a latency of 41.97 $\mu$s for
processing the data in a single trial, which is relatively small compared to
the $\mathcal{O}(ms)$ sampling interval commonly used in neuroscience
experiments. The acceleration enables large-scale real-time experiments and
enhances the capacity to process extensive neural recordings. An AE-based
LFADS model is studied for this work. The VAE-based LFADS should be supported
by hls4ml in the near future.

###### Acknowledgements.

We acknowledge the Fast Machine Learning collective and NSF A3D3 team as an
open community of multi-domain experts and collaborators. This communities
were important for the development of this project. We thank Eli Shlizerman,
Amy Orsborn, Cole Hurwitz, and Michael Nolan for valuable discussions at
different stages of the project. This research was funded in part by National
Science Foundation (NSF) grants No. 1934360 and 2117997, and 1934360.

## CODE AVAILABILITY STATEMENT

We re-implemented the LFADS in Keras functional API. The model can be found in
the ”Tensorflow 2 Keras LFADS” repository available at
https://github.com/XLSMALL/TF2Keras-LFADs, originally developed by Hurwitz et
al. The original work can be found at https://github.com/HennigLab/tndm. The
PTQ model conversion was done using the hls4ml package found at
https://github.com/XLSMALL/hls4ml/tree/LFADs, and we intend to contribute the
QAT code to hls4ml in the near future.

## References

  * (1)
  * Cunningham and Yu (2014) John P Cunningham and Byron M Yu. 2014.  Dimensionality reduction for large-scale neural recordings.  _Nature neuroscience_ 17, 11 (2014), 1500–1509.
  * Developers (2023) TensorFlow Developers. 2023\.  _TensorFlow_.  https://doi.org/10.5281/zenodo.8381573 Specific TensorFlow versions can be found in the ”Versions” list on the right side of this page.¡br¿See the full list of authors ¡a href=”htt ps://github.com/tensorflow/tensorflow/graphs/contr ibutors”¿on GitHub¡/a¿..
  * et al. (2019) Claudionor N. Coelho Jr et al. 2019\.  _QKeras_.  https://github.com/google/qkeras
  * et al. (2015) F. Chollet et al. 2015\.  _Keras_.  https://https://keras.io/
  * FastML Team (2021) FastML Team. 2021\.  _fastmachinelearning/hls4ml_.  https://doi.org/10.5281/zenodo.1201549
  * Hurwitz et al. (2021) Cole Hurwitz, Akash Srivastava, Kai Xu, Justin Jude, Matthew G. Perich, Lee E. Miller, and Matthias H. Hennig. 2021.  Targeted Neural Dynamical Modeling.  arXiv:2110.14853 [q-bio.NC]
  * Klambauer et al. (2017) Günter Klambauer, Thomas Unterthiner, Andreas Mayr, and Sepp Hochreiter. 2017.  Self-Normalizing Neural Networks.  _CoRR_ abs/1706.02515 (2017).  arXiv:1706.02515 http://arxiv.org/abs/1706.02515
  * Pandarinath et al. (2018) Chethan Pandarinath et al. 2018\.  Inferring single-trial neural population dynamics using sequential auto-encoders.  _Nat Methods_ 15, 10 (Oct 2018), 805–815.  https://doi.org/10.1038/s41592-018-0109-9
  * Perich et al. (2018) Matthew G. Perich, Juan A. Gallego, and Lee E. Miller. 2018\.  A Neural Population Mechanism for Rapid Learning.  _Neuron_ 100, 4 (Nov 2018), 964–976.e7.  https://doi.org/10.1016/j.neuron.2018.09.030
  * Sussillo et al. (2016) David Sussillo, Rafal Jozefowicz, L. F. Abbott, and Chethan Pandarinath. 2016.  LFADS - Latent Factor Analysis via Dynamical Systems.  _arXiv_ (22 Aug 2016).  http://arxiv.org/abs/1608.06315 Accessed: Aug. 03, 2023.
  * Xilinx (2023a) Xilinx. 2023a.  _Overview of Arbitrary Precision Fixed-Point Data Types_.  https://docs.xilinx.com/r/en-US/ug1399-vitis-hls/Overview-of-Arbitrary-Precision-Fixed-Point-Data-Types
  * Xilinx (2023b) Xilinx. 2023b.  _Xilinx DSP48E2 block_.  https://docs.xilinx.com/r/en-US/ug958-vivado-sysgen-ref/DSP48E2
  * Yang and Wang (2020) Guangyu Robert Yang and Xiao-Jing Wang. 2020.  Artificial Neural Networks for Neuroscientists: A Primer.  _Neuron_ 107, 6 (2020), 1048–1070.  https://doi.org/10.1016/j.neuron.2020.09.005