22-519

# Uncertainty Propagation and Filtering via the Koopman Operator in
Astrodynamics

Simone Servadio Postdoctoral Associate, Department of Aeronautics and
Astronautics, Massachusetts Institute of Technology, MA 02139, USA. email:
<EMAIL_ADDRESS>William Parker SM Candidate, Department of Aeronautics and
Astronautics, Massachusetts Institute of Technology, MA 02139, USA. email:
<EMAIL_ADDRESS>and Richard Linares Boeing Assistant Professor, Department of
Aeronautics and Astronautics, Massachusetts Institute of Technology, MA 02139,
USA. email<EMAIL_ADDRESS>
###### Abstract

The Koopman Operator (KO) provides an analytical solution of dynamical systems
in terms of orthogonal polynomials. This work exploits this representation to
include the propagation of uncertainties, where the polynomials are modified
to work with stochastic variables. Thus, a new uncertainty quantification
technique is proposed, where the KO solution is expanded to include the
prediction of central moments, up to an arbitrary order. The propagation of
uncertainties is then expanded to develop a new filtering algorithm, where
measurements are considered as additional observables in the KO mathematics.
Numerical simulations in astrodynamics assess the accuracy and performance of
the new methodologies.

##  1 Introduction

The development of capable general methods for propagating the orbital motion
of a satellite subject to strongly nonlinear dynamics is an important
challenge in modern astrodynamics. Along with predicted orbital motion itself,
an awareness of the accumulated uncertainty in the prediction is useful for a
wide variety of operational applications relating to analysis, estimation, and
control tasks. Uncertainty propagation is usually followed by the measurement
update, where the system receives noisy observations from the sensors, and
corrects its predicted state into its estimate [1]. The final result
constitutes a complete filtering algorithm. For the linear and Gaussian case,
the Kalman formulation [2, 3] represents the optimal solution, since the
conditional probability density function (PDF) remains Gaussian at all times
and it is evaluated through Bayes’s rule [4].

However, most of the systems in astrodynamics are affected by nonlinearities,
which leads to a sub-optimal solution for the evaluation of the conditional
PDF [5]. The first solution to address nonlinearities is the Extended Kalman
Filter (EKF) [6], where the estimation error is linearized around the most
current estimate and the Kalman formulation is applied. However, it has been
shown that the linearization assumption is insufficient to correctly propagate
uncertainties for orbit determination applications [7]. A more robust
technique is offered by the Unscented Transformation (UT) [8], where carefully
selected sigma points undergo the full nonlinear dynamics to propagate central
moments as weighted means. The Unscented Kalman Filter (UKF) [9] applies UT
for state and measurement prediction, obtaining superior performance than the
EKF.

The propagation of uncertainties in orbital mechanics has been widely
addressed. Park and Scheeres [10] propagate mean and high order central
moments using state transition tensors (STT) through arbitrary nonlinear
dynamics. They, later [11], complete a filtering algorithm able to include
measurement update. Majji, Turner and Junkins [12] expand this work to update
all the higher central moments. Afterwards, Valli et al. are able to recreate
these results using differential algebra (DA) techniques [13], with
applications to orbit determination and rendezvous maneuvers [14, 15].

This work proposes a novel uncertainty propagation technique, based on the
Koopman Operator (KO). KO enables the re-formulation of strongly nonlinear
dynamics into a linear polynomial framework. The KO theory was initially
developed by Bernard Koopman and further developed by John Von Neumann in the
1930s [16, 17]. Koopman and Von Neumann applied the operator theoretic
perspective used in the development of quantum mechanics [18, 19] to classical
mechanics problems. They introduced the Koopman Operator: a linear operator
that can reformulate a nonlinear dynamical system of finite dimension into a
linear system with an infinite number of dimensions. To leverage the benefit
of working with linear systems, a finite-dimensional approximation of the
dynamics can be attained by observing the projection on a finite subspace
spanned by a set of basis functions. Recent interest in studying the evolution
of nonlinear dynamical systems in fluid mechanics [20], control [21, 22], and
attitude dynamics [23] has led to a fruitful re-discovery of the KO theory for
spectral analysis and linear approximations of these dynamical systems.

More recently, KO has been employed to address computational challenges in
astrodynamics with promising results. The operator theoretic framework was
first introduced to model satellite motion in the proximity of libration
points [24]. Koopman-based perturbation theory has been used to model the
motion of a satellite about an oblate planet [25, 26] and an approximate
analytical solution to the zonal harmonics problem was identified [27]. In
this work, KO approximates PDFs via their central moments, which are
propagated using the KO eigendecompostion of the system with orthogonal
polynomials. The KO analytical solution of the system provides a global
analysis of the dynamics, such that it can overcome the drawbacks shown by
local approximations, such as DA and the Taylor polynomial solutions, where
large uncertainties need to be split, using Automatic Domain Splitting (ADS)
[28] or Gaussian Multiple Models (GMM) [29], to maintain accurate performance.

The paper is divided as follows. The KO theory is presented in the next
section. Afterwards, the new uncertainty propagation technique and the
filtering algorithm are presented. Then, uncertainty quantification
methodology is applied to the complex dynamics of satellite motion under the
circular restricted three body problem (CRTBP) [30], with applications to
orbit determination with limited measurements. Lastly, conclusions are drawn.

##  2 Koopman Operator Theory

A nonlinear dynamical system can be expressed by the initial value problem,
which introduces a set of coupled autonomous ordinary differential equations
(ODEs) of the form

$\left\\{\begin{tabular}[]{l}$\displaystyle\frac{d}{dt}{\bf x}(t)={\bf f}({\bf
x})$\\\ ${\bf x}(t_{0})={\bf x_{0}}$\end{tabular}\right.,$ (1)

where ${\bf x}\in\mathbb{R}^{d}$ is the state which depends on the time
evolution $t$, ${\bf f}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is the
nonlinear dynamics model, $d$ is the number of dimensions in which the problem
is defined, and ${\bf x_{0}}$ is the initial condition of the system at time
$t_{0}$. The Koopman Operator $(\mathcal{K})$, propagates $\mathcal{G}({\bf
x})$, the functions describing observables of the state.

The KO is an infinite-dimensional operator. A result of this property is that
the vector space of the observable functions, defined here as $\mathcal{F}$
(where $\mathcal{G}({\bf x})\in\mathcal{F}$) is also infinite-dimensional. It
follows that if $g\subseteq\mathcal{G}({\bf x})$ is a given observable in this
space, the evolution of $g$ in the dynamical system can be found by

$\mathcal{K}\left(g({\bf x})\right)=\frac{d}{dt}g({\bf x})=\left(\nabla_{{\bf
x}}g({\bf x})\right)\frac{d}{dt}{\bf x}(t)=\left(\nabla_{{\bf x}}g({\bf
x})\right){\bf f}({\bf x}),$ (2)

where $\nabla_{{\bf x}}g=(\partial g/\partial x_{1},\partial g/\partial
x_{2},\dots,\partial g/\partial x_{d})$. More generally, any observable
governed by the dynamical system can be evolved using the KO in this way;

$\mathcal{K}\left(\cdot\right)=\left(\nabla_{{\bf x}}\cdot\right){\bf f}({\bf
x}).$ (3)

Koopman’s initial work [16] outlined the linearity property of the KO, which
implies

$\mathcal{K}\left(\beta_{1}g_{1}({\bf x})+\beta_{2}g_{2}({\bf
x})\right)=\beta_{1}\mathcal{K}\left(g_{1}({\bf
x})\right)+\beta_{2}\mathcal{K}\left(g_{2}({\bf x})\right),$ (4)

for any pair of observables $g_{1}\subseteq\mathcal{G}({\bf x})$ and
$g_{2}\subseteq\mathcal{G}({\bf x})$ and any arbitrary constants $\beta_{1}$
and $\beta_{2}$. The linearity of the KO allows for efficient computation, but
the infinite-dimensionality of the operator makes it impossible to capture the
evolution of all measurement functions in a Hilbert space [21]. This issue can
be mitigated, however, by instead capturing the evolution of the system on a
finite subspace spanned by a finite set of basis functions. Performing this
step truncates the KO to a finite subspace $\mathcal{F}_{D}$ of dimension $m$,
where $\mathcal{F}_{D}\in\mathcal{F}$. This subspace $\mathcal{F}_{D}$ can be
spanned by any set of eigenfunctions $\phi_{i}\in\mathcal{F}_{D}$, with
$i\in\\{1,2,\dots,m\\}$, defined as

$\mathcal{K}\left(\phi_{i}({\bf x})\right)=\frac{d}{dt}\phi_{i}({\bf
x})=\lambda_{i}\phi_{i}({\bf x}),$ (5)

where $\lambda_{i}$ are the eigenvalues associated with the eigenfunctions
$\phi_{i}$, and $m$ is the number of eigenfunctions chosen to represent the
space. To linearize the system, a transformation of variables is performed
using the Koopman eigenfunctions, ${\bf\Phi}({\bf x})=\left(\phi_{1}({\bf
x}),\dots,\phi_{m}({\bf x})\right)^{T}\in\mathcal{F}_{D}$. Then, using the
relation in Eq. (5), the evolution of ${\bf\Phi}$ can be expressed as

$\mathcal{K}\left(\bf\Phi\right)=\frac{d}{dt}{\bf\Phi}=\Lambda{\bf\Phi},$ (6)

where $\Lambda=\text{diag}([\lambda_{1},\dots,\lambda_{m}])$ is the diagonal
matrix containing the eigenvalues of the system in $\mathcal{F}_{D}$. The
solution of Eq. (6) is

${\bf\Phi}(t)=\exp(\Lambda t){\bf\Phi}(t_{0}),$ (7)

where ${\bf\Phi}(t_{0})$ is the value of the eigenfunctions at the initial
time $t_{0}$. Once the KO eigenfunctions are computed, Eq. (7) evolves them in
time. For the general case, there is interest in the identity observable, that
is, ${\bf g}({\bf x})={\bf x}$.

###  2.1 Computing the Koopman Eigenfunctions via Galerkin Method

Following the theory developed in [25], we noticed that the KO gives rise to a
set of linear Partial Differential Equations (PDE) for the eigenfunctions:

$\mathcal{K}(\phi_{i})=\left(\nabla_{x}\phi_{i}({\bf x})\right){\bf f}({\bf
x})=\lambda_{i}.$ (8)

The Galerkin method offers an approximation of the solution of this PDE. Let
us select the Legendre polynomials as the set of basis functions, since they
create an orthogonal basis and have computational advantages in the
computation of the Koopman matrix [31]. The main idea of the methodology is to
represent any function of the space as a linear combination of the basis
functions. This is achieved by defining an inner product:

$\langle f,g\rangle=\int_{\Omega}f({\bf x})g({\bf x})w({\bf x})d{\bf x},$ (9)

where $f$ and $g$ be two arbitrary functions, and $w({\bf x})$ is a positive
weighting function defined on the space domain $\Omega$. For the case of
Legendre polynomials, the weighting function is a constant $w({\bf x})=1$, and
the domain for each variable ranges between $[-1,1]$. In addition, the
normalized Legendre polynomials are defined such that:

$\langle L_{i},L_{j}\rangle=\int_{\Omega}L_{i}({\bf x})L_{j}({\bf x})w({\bf
x})d{\bf x}=\delta_{ij},$ (10)

where $L_{i}$ and $L_{j}$ with $\\{i,j\\}\in\\{1,\dots,m\\}$ are two given
normalized Legendre polynomials from the set of basis functions selected, and
$\delta_{ij}$ is Kronecker’s delta. By using this set of orthogonal
multivariate polynomials, the KO eigenfunctions can be represented as a series
expansion in terms of this set of basis functions, particularly:

$\phi_{i}({\bf x})=\sum_{\ell=1}^{\infty}c_{i\ell}L_{\ell}({\bf
x})\approx\sum_{\ell=1}^{m}c_{i\ell}L_{\ell}({\bf x})={\bf c}_{i}{\bf
L}^{T}({\bf x})={\bf L}({\bf x}){\bf c}_{i}^{T},$ (11)

where $c_{i\ell}$ is a constant coefficient associated with the eigenfunction
$\phi_{i}$ and the basis $L_{\ell}$, and ${\bf c}_{i}$ and ${\bf L}$ are two
row vectors containing the set of coefficients $c_{i\ell}$ and the whole set
of basis functions respectively. Note that although the series is infinite, a
truncation was performed using $m$ different basis functions, and thus, this
represents an approximation of the eigenfunctions.

The PDE for the computation of the Koopman eigenfunction, Eq. (8), can be
approximated using a Galerkin method and the series expansion from Eq. (11),
where each function that undergoes the KO is projected onto the
$m$-dimensional subspace $\mathcal{F}_{D}$. The time evolution of each $i$-th
eigenfunction can be grouped into a diagonal system, with dimensions $m$.
After some mathematical manipulations described in [25], [31], and [29], the
diagonalization of the system is obtained as:

$CK=\Lambda C,$ (12)

where $C$ is the matrix containing all the coefficients $c_{ij}$, $\Lambda$ is
the diagonal matrix of eigenvalues $\lambda_{i}$, and $K$ is the Koopman
matrix, a $m\times m$ sized matrix whose components are:

$K_{ij}=\langle\nabla_{\bf x}L_{i}({\bf x}){\bf f}({\bf x}),L_{j}({\bf
x})\rangle.$ (13)

Note that Eq. (12) is in fact defining a eigendecomposition of matrix $K$,
where $C$ is the matrix containing the left eigenvectors of $K$, and $\Lambda$
the associated eigenvalues. Therefore, it is possible to first obtain $K$ in
closed form solution from Eq (13), and then perform the eigendecomposition to
obtain $C$, and thus, the approximated eigenfunctions $\phi_{i}$ of the system
through Eq. (11).

###  2.2 Koopman modes

The Koopman modes are the representation of the observables into the set of
eigenfunctions ${\bf\Phi}$ of the system. Thus, let ${\bf g}({\bf x})$ be the
observables of interest, function of the state. Using the Galerkin method,
these observables can be approximated as an expansion in the set of
eigenfunctions

$g_{i}({\bf x})\approx\sum_{\ell=1}^{m}b_{i\ell}\phi_{\ell}({\bf x}),$ (14)

where $b_{i\ell}$ are the coefficients, computed by projecting the $i$-th
observable onto the $\ell$-th basis functions:

$b_{i\ell}=\left\langle g_{i},\phi_{\ell}\right\rangle.$ (15)

Therefore, the observables can be represented in matrix notation as ${\bf
g}({\bf x})=B{\bf\Phi}({\bf x})$, where $B$ is a matrix of size $q\times m$
containing the coefficients $b_{i\ell}$, and $q$ represents the number of
observables. The time evolution of ${\bf g}({\bf x})$ is expressed by the time
evolution of the eigenfunctions, or, equivalently, in terms of the basis
functions:

${\bf g}({\bf x}(t))=B{\bf\Phi}({\bf x}(t))=B\exp(\Lambda t){\bf\Phi}({\bf
x}(t_{0}))=B\exp(\Lambda t)C{\bf L}({\bf x}(t_{0})),$ (16)

where ${\bf x}(t_{0})$ are the initial condition values of the states.

Note that the observables can be also projected onto the set of basis
functions ${\bf L}$:

$g_{i}({\bf x})\approx\sum_{\ell=1}^{m}a_{i\ell}L_{\ell}({\bf
x})\rightarrow{\bf g}({\bf x})=A{\bf L}({\bf x})\quad\text{where}\quad
a_{i\ell}=\left\langle g_{i},L_{\ell}\right\rangle,$ (17)

where $A_{ij}=a_{ij}$, whose solution under the system of differential
equations is

${\bf g}({\bf x}(t))=A{\bf L}({\bf x}(t))=AC^{-1}\exp(\Lambda t)C{\bf L}({\bf
x}(t_{0})).$ (18)

For a more detailed analysis of the analytical KO solution with orthogonal
polynomials, with a step-by-step explanation of the methodology, the reader
can refer to [32, 25, 31, 30].

##  3 Uncertainty Propagation with the Koopman Operator

The solution of the system can be formulated using the identity observable
(the observable g(x)=x) and the KO eigendecomposition as in the following
expression:

${\bf g}({\bf x}(t_{0}),t)=AC^{-1}\exp(\Lambda t)C{\bf L}({\bf x}(t_{0}))$
(19)

Note that the solution is a function of the initial state vector at time,
$t_{0}$, and the current time, $t$. When propagating the state up to the final
time $t_{f}$, the Koopman solution does not require any numerical integration,
and the final state ${\bf x}_{f}={\bf x}(t_{f})$, is directly connected to the
value of the state at the initial time ${\bf x}_{0}={\bf x}(t_{0})$,

${\bf x}_{f}({\bf x}_{0})=AC^{-1}\exp(\Lambda t_{f})C{\bf L}({\bf x}_{0}).$
(20)

which gives the polynomial dependence of the propagated states on the initial
condition. Equation (20) is a polynomial in the variable ${\bf x_{0}}$ and it
constitutes the forward map in time of the dynamics of the system.

When propagating PDFs, it is more convenient to work with the deviations from
the mean than with the state itself as the variable [33]. This way, the
central moments are easier to calculate and provide a more intuitive
understanding of the uncertainties of the system when compared to raw moments
[34]. Let us assume $\hat{\bf x}_{0}$ is the known mean of the initial PDF
that will be propagated, in terms of its central moments, to time $t_{f}$. It
is possible to express the Koopman solution of the state as a function of the
deviation vector $\delta{\bf{x}}_{0}$ from the initial condition through the
variable shift ${\bf{x}}_{0}=\hat{\bf x}_{0}+\delta{\bf{x}}_{0}$:

$\displaystyle\tilde{\bf x}_{f}(\delta{\bf{x}}_{0})$
$\displaystyle=AC^{-1}\exp(\Lambda t_{f})C{\bf L}(\hat{\bf x}_{0}+\delta{\bf
x}_{0})$ (21) $\displaystyle=AC^{-1}\exp(\Lambda t_{f})C\tilde{\bf
L}(\delta{\bf x}_{0})$ (22)

This state polynomial represents the same map as Equation (20), but with the
variable centered at the current mean. Thus, the Legendre polynomials ${\bf
L}$ and $\tilde{\bf L}$ describe the solution, but their coefficients are
different such that their variable dependency has been shifted by a value
$\hat{\bf x}_{0}$. It can be observed that the equality ${\bf L}(\hat{\bf
x}_{0})=\tilde{\bf L}({\bf 0})$ holds.

The deviations are determining all the possible outcomes of the propagated
state around the center $\bar{\bf x}_{f}=\tilde{\bf x}_{f}({\bf 0})={\bf
x}_{f}(\hat{\bf x}_{0})$, which is the time propagation of the initial mean
$\hat{\bf x}_{0}$. In the particular case of linear dynamics and Gaussian
distributions, the mean of the propagated PDF, $\hat{\bf x}_{f}$, coincides
with the propagation of the initial mean and center of the map: $\hat{\bf
x}_{f}=\bar{\bf x}_{f}$. However, due to non-linearities, the polynomial map
that connects deviations from $t_{0}$ to $t_{f}$ loses its Gaussian property,
and the propagated central moments need to be evaluated using expectations.

Let us assume an initial Gaussian PDF of the state, perfectly described by its
mean $\hat{\bf x}_{0}$ and covariance ${\bf P}_{0}$:

${\bf x}_{0}\sim\mathcal{N}(\hat{\bf x}_{0},{\bf P}_{0})$ (23)

Thanks to the Isserlis formulation [35], it is possible to evaluate all
central moments of the Gaussian distribution up to order $\psi$. The mean of
the state PDF at time $t_{f}$ is calculated by applying the expected value
operator $\mathbb{E}[\cdot]$ on the polynomial $\tilde{\bf
x}_{f}(\delta{\bf{x}}_{0})$

$\hat{\bf x}_{f}=\mathbb{E}[\tilde{\bf{x}}_{f}(\delta{\bf{x}}_{0})]=\bar{\bf
x}_{f}+\sum\chi_{i}\mathbb{E}[\delta{\bf x}_{0,i}]$ (24)

where $\chi_{i}$ indicates the coefficient of the $i$-th monomial in the state
polynomial, and $\delta{\bf x}_{0,i}$ the relative deviation. Through linear
operator property, the evaluation of the mean works directly on each single
monomial of the polynomial. The expected value is evaluated by substituting,
for each expectation of deviations, the relative central moment previously
obtained through the Isserlis formulation [35]. Thus, having the state
polynomial as a function of the deviations from the mean leads to a trivial
evaluation of the expected values, where different powers of the deviations
are easily connected to the relative central moments of the distribution.

The same procedure can be adopted to evaluate central moments of any order.
Starting from the propagated covariance, the expectation is given by:

${\bf P}_{f}=\mathbb{E}[(\tilde{\bf{x}}_{f}(\delta{\bf{x}}_{0})-\hat{\bf
x}_{f})(\tilde{\bf{x}}_{f}(\delta{\bf{x}}_{0})-\hat{\bf x}_{f})^{T}]$ (25)

The argument of the expectation is a polynomial, where deviations get, once
again, substituted by the initial central moments of the original state PDF,
approximated as Gaussian. Therefore, it is possible to approximate the
propagated central moments up to any arbitrary order. Let us indicate with
$\tilde{\bf{x}}_{f,i}$ the $i$-th component of the state polynomial vector;
then, each component of the propagated state skewness ${\bf S}_{f}$ (the
third-order central moment) and kurtosis ${\bf\Sigma}_{f}$ (fourth-order
central moment) can be approximated as:

$\displaystyle{\bf S}_{f,ijl}$
$\displaystyle=\mathbb{E}[(\tilde{\bf{x}}_{f,i}(\delta{\bf{x}}_{0})-\hat{\bf
x}_{f,i})(\tilde{\bf{x}}_{f,j}(\delta{\bf{x}}_{0})-\hat{\bf
x}_{f,j})(\tilde{\bf{x}}_{f,l}(\delta{\bf{x}}_{0})-\hat{\bf x}_{f,l})]$ (26)
$\displaystyle{\bf\Sigma}_{f,ijlm}$
$\displaystyle=\mathbb{E}[(\tilde{\bf{x}}_{f,i}(\delta{\bf{x}}_{0})-\hat{\bf
x}_{f,i})(\tilde{\bf{x}}_{f,j}(\delta{\bf{x}}_{0})-\hat{\bf
x}_{f,j})(\tilde{\bf{x}}_{f,l}(\delta{\bf{x}}_{0})-\hat{\bf
x}_{f,l})(\tilde{\bf{x}}_{f,m}(\delta{\bf{x}}_{0})-\hat{\bf x}_{f,m})]$ (27)

Therefore, the state PDF is propagated in time in terms of its central
moments, evaluated using polynomial substitution for each state component.

It is worth noticing that performing the propagation of moments does not
require any numerical integration. The dynamics have been solved and
represented as a linear combination of well-defined eigenfunctions through the
KO solution, which are connected to the basis polynomials. The propagated
states, and relative moments, are evaluated directly from the Legendre
polynomials once the final time is selected. This eliminates the computational
expensive propagation required by other methods.

It is possible now to complete a filtering algorithm by adding a linear update
(or polynomial) to the propagated state probability density function according
to the Kalman formulation [2, 3].

##  4 The Koopman Operator Filter

The propagation of uncertainties can be completed into a full filtering
algorithm by adding the measurement update. As such, the propagation of the
state mean and covariance can be discretized for the prediction of the state
distribution from time step $\\{t_{k},{\bf{x}}_{k},{\bf P}_{k}\\}$ to time
step $\\{t_{k+1},{\bf{x}}_{k+1},{\bf P}_{k+1}\\}$, according to the
propagation equations described beforehand:

$\displaystyle\tilde{\bf x}_{k+1}(\delta{\bf{x}}_{k})$
$\displaystyle=AC^{-1}\exp(\Lambda t_{k+1})C{\bf L}(\hat{\bf x}_{k}+\delta{\bf
x}_{k})=AC^{-1}\exp(\Lambda t_{k+1})C\tilde{\bf L}(\delta{\bf x}_{k})$ (28)
$\displaystyle\hat{\bf x}_{k+1}^{-}$
$\displaystyle=\mathbb{E}[\tilde{\bf{x}}_{k+1}(\delta{\bf{x}}_{k})]$ (29)
$\displaystyle{\bf P}_{k+1}^{-}$
$\displaystyle=\mathbb{E}[(\tilde{\bf{x}}_{k+1}(\delta{\bf{x}}_{k})-\hat{\bf
x}_{k+1}^{-})(\tilde{\bf{x}}_{k+1}(\delta{\bf{x}}_{k})-\hat{\bf
x}_{k+1}^{-})^{T}]$ (30)

where $\hat{\bf x}_{k+1}^{-}$ and ${\bf P}_{k+1}^{-}$ represent the predicted
mean and covariance, respectively. The state polynomial at time step $t_{k+1}$
is evaluated efficiently using the eigendecomposition of the dynamics of the
system performed previously, (offline), during the Koopman spectral analysis.
Therefore, the state polynomials of deviations of each state component at any
time step are obtained directly from the basis functions of the system.

The correction step (or measurement update) takes its name from modifying the
predicted state mean and covariance, $\hat{\bf x}_{k+1}^{-}$ and ${\bf
P}_{k+1}^{-}$, into their updated values, $\hat{\bf x}_{k+1}^{+}$ and ${\bf
P}_{k+1}^{+}$. The filter receives information about the state of the system
through a series of observations ${\bf y}_{k+1}$, which relationship with the
state itself is known;

${\bf y}_{k+1}={\bf h}({\bf x}_{k+1})+{\boldsymbol{\eta}}_{k+1}$ (31)

where ${\bf h}(\cdot)$ is the measurement function and
${\boldsymbol{\eta}}_{k+1}$ is the measurement noise, with known covariance
matrix ${\bf R}_{k+1}$. The aim of the algorithm is to best merge information
from the prediction step and the observations to obtain the most accurate
estimate of the system by “filtering out” the noise.

The Koopman Operator Filter (KOF) develops the measurement update with the
same polynomial algebra applied for the propagation of the state
uncertainties. After all, the measurements ${\bf y}_{k+1}$ are just but a more
complicated set of observables that can be represented by the set of basis
functions ${\bf L}$, or $\tilde{\bf L}$. In an analogous way matrix $A$ was
calculated for the identity observables to get the state of the system, a new
matrix $A_{y}$ is evaluated to express in a polynomial way the measurement
function, as a linear combination of the basis functions: ${\bf y}({\bf
x})=A_{y}{\bf L}({\bf x})$. Each entry of $A_{y}$ is evaluated using the
Galerkin method with the inner product between functions, likewise any other
Koopman matrix presented in the paper. After shifting the center of the
polynomials using $\tilde{\bf L}$ over ${\bf L}$, the final result is the
measurement polynomials as a function of the state deviation vector

$\displaystyle\tilde{\bf y}_{k+1}(\delta\bf{x}_{k})$
$\displaystyle=A_{y}C^{-1}\exp(\Lambda t_{k+1})C{\bf L}(\hat{\bf
x}_{k}^{+}+\delta{\bf x}_{k})=A_{y}C^{-1}\exp(\Lambda t_{k+1})C\tilde{\bf
L}(\delta{\bf x}_{k})$ (32)

These polynomials express the distribution of the measurements as a function
of the deviations of the state from the original time step $t_{k}$, therefore,
they represent the Koopman polynomial map of deviations from the state space
to the measurement space. This map directly connects to any possible state
deviation, the relative outcome in the measurement space: the map tells how
the measurement PDF is connected, and shaped, with respect to the state PDF.

In order to apply the Galerkin method, and evaluate the inner product between
the basis functions and the measurement observables, the measurement equations
need to be expressed as a polynomial approximation of ${\bf h}(\cdot)$. For
small deviations, a Taylor approximation furnish a very accurate
representation of the measurement map around the current center. However, any
other polynomial representation is a valid solution and its selection depends
on the application.

Following the Kalman update formulation of an estimator with a linear
dependence on the measurements, the Koopman Operator Filter evaluates the
predicted measurement mean by applying the expectation operator directly on
the measurement polynomials

$\displaystyle\hat{\bf y}_{k+1}$
$\displaystyle=\mathbb{E}[\tilde{\bf{y}}_{k+1}(\delta{\bf{x}}_{k})]$ (33)

where, once again, the mean of each component is calculated by substituting
the state central moments, evaluated with Isserlis, on each monomial of the
series. Continuing the analogy with the prediction step, the measurement
covariance is evaluated as

$\displaystyle{\bf P}_{yy}$
$\displaystyle=\mathbb{E}[(\tilde{\bf{y}}_{k+1}(\delta{\bf{x}}_{k})-\hat{\bf
y}_{k+1})(\tilde{\bf{y}}_{k+1}(\delta{\bf{x}}_{k})-\hat{\bf
y}_{k+1})^{T}]+{\bf R}_{k+1}$ (34)

having considered the influence of the measurement noise. The state-
measurement cross-covariance is evaluated as

$\displaystyle{\bf P}_{xy}$
$\displaystyle=\mathbb{E}[(\tilde{\bf{x}}_{k+1}(\delta{\bf{x}}_{k})-\hat{\bf
x}_{k+1}^{-})(\tilde{\bf{y}}_{k+1}(\delta{\bf{x}}_{k})-\hat{\bf
y}_{k+1})^{T}]$ (35)

so that the Kalman gain is

$\displaystyle G$ $\displaystyle={\bf P}_{xy}{\bf P}_{yy}^{-1}$ (36)

as it has been proven by the orthogonality principle [33]. The update of the
state estimate and the relative covariance is performed using the Kalman
update equations

$\displaystyle\hat{\bf x}_{k+1}^{+}$ $\displaystyle=\hat{\bf
x}_{k+1}^{-}+G({\bf y}_{k+1}-\hat{\bf y}_{k+1})$ (37) $\displaystyle{\bf
P}_{k+1}^{+}$ $\displaystyle={\bf P}_{k+1}^{-}-G{\bf P}_{yy}G^{T}$ (38)

The filtering algorithm has ended and it is ready to start the next iteration,
having approximated the updated state distribution as a Gaussian PDF with mean
the estimate $\hat{\bf x}_{k+1}^{+}$ and covariance ${\bf P}_{k+1}^{+}$. The
following step will start by evaluating the state polynomial at the new time
step and shifting the basis functions at the current estimate. The
expectations will be evaluated by applying the Isserlis formulation with the
updated covariance matrix.

The polynomial nature of the KOF makes it suitable for the implementation to
nonlinear estimator, where the updated estimate is not a linear function of
the measurement, but it has a high-order polynomial dependence on the
observation that increases accuracy [34, 36].

##  5 The Circular Restricted Three Body Problem

The circular restricted three body problem (CRTBP) equations can be expressed
as:

$\displaystyle\ddot{x}$ $\displaystyle=\frac{\partial\Omega}{\partial
x}+2\dot{y},$ (39a) $\displaystyle\ddot{y}$
$\displaystyle=\frac{\partial\Omega}{\partial y}-2\dot{x},$ (39b)
$\displaystyle\ddot{z}$ $\displaystyle=\frac{\partial\Omega}{\partial z},$
(39c)

where $\Omega=1/2\left(x^{2}+y^{2}\right)+(1-\mu)/r_{1}+\mu/r_{2}$, $\mu$ is
the ratio between the mass of the secondary over the total mass of the system,
and $r_{1}$ and $r_{2}$ are the distances to the primary and secondary masses.
The distance $\gamma$ from the libration points, $\mathcal{L}_{j}$, to the
primary body is found by solving the Euler quintic equation:

$\displaystyle\gamma_{j}^{5}\mp\left(3-\mu\right)\gamma_{j}^{4}+\left(3-2\mu\right)\gamma_{j}^{3}-\mu\gamma_{j}^{2}\pm
2\mu\gamma_{j}-\mu=0,\,j=1,2$ (40)

Following [37], the equations of motion for a satellite moving near
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ can be expressed by translating the
origin to the equilibrium point:

$\bar{x}=\frac{x-1+\mu\pm\gamma}{\gamma},\quad\quad\bar{y}=\frac{y}{\gamma},\quad\quad\bar{z}=\frac{z}{\gamma},$
(41)

where, in the $\pm$ symbol, the upper sign relates to the dynamic about
$\mathcal{L}_{1}$, while the lower sign is for $\mathcal{L}_{2}$. The CRTBP
position variables $\\{x,y,z\\}$ have been shifted and scaled to a new set of
coordinates $\\{\bar{x},\bar{y},\bar{z}\\}$, where the new origin is either
$\mathcal{L}_{1}$ or $\mathcal{L}_{2}$, and they are normalized such that the
distance between $\mathcal{L}_{1}$ or $\mathcal{L}_{2}$ and the primary is 1.
For convenience of notation, $x$, $y$, $z$ will refer to $\bar{x}$, $\bar{y}$,
and $\bar{z}$ variables henceforth.

Richardson [37] showed that the full nonlinear CRTBP equations can be expanded
in terms of polynomials, especially Legendre polynomials, $P_{n}$, for their
computational advantages. Following the developments of [37], the CRTBP around
$\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ can be represented as:

$\displaystyle\ddot{x}-2\dot{y}-(1+2c_{2})x$
$\displaystyle=\frac{\partial}{\partial x}\sum_{n\geqslant
3}c_{n}\rho^{n}P_{n}\left(\frac{x}{\rho}\right),$ (42a)
$\displaystyle\ddot{y}+2\dot{x}+(c_{2}-1)y$
$\displaystyle=\frac{\partial}{\partial y}\sum_{n\geqslant
3}c_{n}\rho^{n}P_{n}\left(\frac{x}{\rho}\right),$ (42b)
$\displaystyle\ddot{x}+c_{2}z$ $\displaystyle=\frac{\partial}{\partial
z}\sum_{n\geqslant 3}c_{n}\rho^{n}P_{n}\left(\frac{x}{\rho}\right),$ (42c)

where $\rho^{2}=x^{2}+y^{2}+z^{2}$ and the coefficients $c_{n}$ are given by

$c_{n}=\frac{1}{\gamma^{3}}\left((\pm)^{n}\mu+(-1)^{n}\frac{(1-\mu)\gamma^{n+1}}{(1\mp\gamma)^{n+1}}\right),$
(43)

with upper sign for $\mathcal{L}_{1}$ and the lower one for $\mathcal{L}_{2}$.
This expression can be expressed using an elegant formulation [38], that uses
the recursive property of creating Legendre polynomials. Considering the
function

$T_{n}(x,y,z)=\rho^{n}P_{n}\left(\frac{x}{\rho}\right),$ (44)

where $T_{n}$ is a homogeneous polynomial of degree $n$ in variable $x$, $y$,
and $z$, and using the Legendre polynomial three-term recursive formula, it
can be shown that:

$T_{n}=\frac{2n-1}{n}xT_{n-1}-\frac{n-1}{n}\left(x^{2}+y^{2}+z^{2}\right)T_{n-2},$
(45)

where the recursion starts with $T_{0}=1$ and $T_{1}=x$. This work uses this
formula to define the polynomial representation of the CRTBP model.

###  5.1 Hamiltonian in Normal Form

Jorba and Masdemont [39] provide a summary for the evaluation of the normal
form of the CRTBP Hamiltonian. The equations of motion (42) come from applying
Hamilton’s equation to the Hamiltonian

$H=\frac{1}{2}\left(p_{x}^{2}+p_{y}^{2}+p_{z}^{2}\right)+yp_{x}-xp_{y}-\sum_{n\geq
2}c_{n}\rho^{n}P_{n}\left(\frac{x}{\rho}\right),$ (46)

where the variables $p_{x}=\dot{x}-y$, $p_{y}=\dot{y}+x$, and $p_{z}=\dot{z}$
are defined as pseudo momenta.

The normal form is obtained by applying a symplectic linear change of
variables, where parameters are found by studying the linearized terms of the
CRTBP dynamics, in terms of its eigenvalues and eigenvectors. Following the
procedure described in [39], after applying linear transformations, the
Hamiltonian reaches the form

$H_{2}=\lambda_{1}xp_{x}+\frac{\omega_{1}}{2}\left(y^{2}+p_{y}^{2}\right)+\frac{\omega_{2}}{2}\left(y^{2}+p_{y}^{2}\right)-\sum_{n\geq
2}c_{n}\rho^{n}P_{n}\left(\frac{x}{\rho}\right),$ (47)

where the constants are defined as

$\lambda_{1}^{2}=\frac{c_{2}-2-\sqrt{9c_{2}^{2}-8c_{2}}}{2},\quad\omega_{1}^{2}=\frac{c_{2}-2+\sqrt{9c_{2}^{2}-8c_{2}}}{2},\quad\omega_{2}^{2}=c_{2}.$
(48)

The evaluation of these constants comes from the study of the characteristic
equation of the linearized system, where the $x$-$y$ motion characteristic
equation has four roots, two real and two imaginary roots given by
$(\pm\lambda_{1},\pm\omega_{1})$, while the $z$ direction motion (decoupled)
is characterized by two imaginary eigenvalues $\pm\omega_{2}$.

Finally, the Hamiltonian in normal form is obtained by applying a complex
normal transformation:

$\displaystyle x=q_{1},\quad y=\frac{q_{2}+\sqrt{-1}p_{2}}{\sqrt{2}},\quad
z=\frac{q_{3}+\sqrt{-1}p_{3}}{\sqrt{2}},$ (49a) $\displaystyle
p_{x}=p_{1},\quad p_{y}=p_{2},\quad
p_{z}=\frac{\sqrt{-1}q_{3}+p_{3}}{\sqrt{2}}.$ (49b)

Using this new set of variables, the complex normal form Hamiltonian is given
by

$H_{3}=\lambda_{2}q_{1}p_{1}+\sqrt{-1}\omega_{1}q_{2}p_{2}+\sqrt{-1}\omega_{2}q_{3}p_{3}-\sum_{n\geq
2}c_{n}\rho^{n}P_{n}\left(\frac{x}{\rho}\right),$ (50)

where $x$ and $\rho$ can be computed in terms of the new variables using the
previous transformations. Hamiltonian $H_{3}$ decouples the variables in the
linear leading term and it decouples of effects of the stable and unstable
manifolds. The KOF uses the equations of motion evaluated from $H_{3}$
(obtained using Hamilton’s equations) to perform the propagation of
uncertainties, and subsequently it transforms variables back in the original
set. However, every truth model used for performance and comparison purposes
has been embedded with the original set of ODE described in system (39).

###  5.2 Prediction of Central Moments

The correct propagation of uncertainties in time is offered on an
$\mathcal{L}_{1}$ halo orbit application for the Earth-Moon system,
characterized by $\mu=0.012153281419431$ and $\gamma=0.150944856782713$. The
novel technique based on the KO adapts the spectral decomposition of the
dynamics in terms of its eigenfunctions to accurately predict the state
central moments at any given time (without performing any numerical
propagation). Consider the halo orbit connected to the initial condition

$\displaystyle{\bf{x}}_{0}=\left[0.823376807050253\quad 0\quad
0.001386166961157\quad 0\quad 0.126366690232230\quad 0\right]^{T}$ (51)

and with given initial covariance matrix ${\bf P}_{0}=\sigma_{0}^{2}{\bf
I}_{6\times 6}$, with $\sigma_{0}=0.0001$. This initial distribution shows the
exponential divergence nature of the CRTBP. Figure 1 portraits the evolution
of $10^{4}$ orbits randomly chosen form the Gaussian PDF
$\mathcal{N}({\bf{x}}_{0},{\bf P}_{0})$, up to the final time $t_{f}=2$.

Figure 1: Monte Carlo propagation of $10^{4}$ orbits for $t_{f}=2$ given the
initial PDF $\mathcal{N}({\bf{x}}_{0},{\bf P}_{0})$. Projection on the $(x,y)$
plane.

The projection of these orbits into the $(x,y)$ plane illustrates the
divergence of the states. The initial small covariance matrix rapidly
increases and the figure shows the family of orbits form the conical shave
connected to a chaotic system. Therefore, even before reaching half of a
revolution, most of the orbits already strongly separate from the mean
(periodic) orbit, in black. Thus, in just a short time span, small
perturbations from the mean at the initial time lead to extremely large
deviations, that make the prediction of uncertainties a challenging problem.

The newly developed technique has been applied to the CRTBP in the Earth-Moon
system and results have been validated through a Monte Carlo analysis, where
the predicted central moments of the state have been compared to the effective
moments evaluated from all of the Monte Carlo runs. Figure 2 shows the time
evolution of the uncertainties represented with the state raw moments. Indeed,
the continuous lines represents the true evolution of the central moments of
the system: the mean in black, the covariance, as $3{\boldsymbol{\sigma}}$
standard deviation, in blue, the skewness in cyan and the kurtosis in magenta.
The values of skewness and kurtosis are evaluated similarly to the covariance
of the system, with an analogy to the standard deviation:

$\displaystyle\sigma_{j}=\left(\dfrac{\sum_{i=1}^{N}\left(x_{i,j}-\hat{x}_{j}\right)^{2}}{N-1}\right)^{1/2}\
\ \forall j=1,\dots,n$ (52)
$\displaystyle\sigma_{skew,j}=\left(\dfrac{\sum_{i=1}^{N}\left(x_{i,j}-\hat{x}_{j}\right)^{3}}{N-1}\right)^{1/3}\
\ \forall j=1,\dots,n$ (53)
$\displaystyle\sigma_{kurt,j}=\left(\dfrac{\sum_{i=1}^{N}\left(x_{i,j}-\hat{x}_{j}\right)^{4}}{N-1}\right)^{1/4}\
\ \forall j=1,\dots,n$ (54)

where ${N}$ is the number of Monte Carlo runs. Therefore, for each component
of the state, the high order central moment is evaluated by analysing all of
the $10^{4}$ simulations. On the contrary, the KO uncertainty propagation
predicts the central moments of the state PDF every 0.2 years. The values of
the predicted moments are reported, with the same colors used for the Monte
Carlo, as points. The correct prediction of the mean and uncertainties of the
system is asserted by the overlapping of the points over the continuous lines.
The accuracy of the estimation of the moments decreases as the integrating
time increases, meaning that, for extremely large propagation times, the
prediction of the state PDF is unfeasible. The CRTBP is strongly diverging:
thus, high order central moments after long time steps require an elevated
amount of runs in order to have correct indicating values. Moreover, the
figure shows the high level of accuracy in the first half of the simulation,
for times shorter than 1 year, zoomed in in the figure. Station keeping
applications, in filtering, usually receive measurements in the section where
KO shows an extremely accurate prediction of moments. The KO technique can
predict central moments up to any arbitrary high order. However, due to the
fact that odd moments can be negative, the accuracy of the prediction of even
moments is higher than their odd counterpart. Indeed, some values of the
skewness, in the figure, are slightly off by the end of the simulation, where
the state PDF largely expands without control. In the remaining sectors, the
prediction of odd moments outperforms any other uncertainty prediction
technique based on the Gaussian approximation of the state distribution (such
as UT), since they provide a null value for any odd moment.

Figure 2: Prediction of raw moments of the state PDF for a Halo orbit orbiting
around $\mathcal{L}_{1}$ in the earth-Moon system.

The central moments estimation can be appreciated by reporting the
uncertainties as around their mean. Indeed, while Figure 2 shows the correct
prediction of the family of orbits, in terms of its overall possible
trajectories, Figure 3 focuses on the spread around the mean. Therefore,
Figure 3 displays similar results as Figure 2, but with a new prospective
centered at the current mean. Figure 3 shows the chaotic nature of the
problem, as the uncertainties, in terms of their $\sigma$ value, increase
exponentially and the spread of all the possible resolutions of the state
creates a conical shape. The KO prediction is able to keep track of the
central moments of the system, especially during the first half of the
simulation. Both the state covariance and kurtosis expand without any
boundaries. The star points show the correct estimation performed by the KO
technique. The new methodology is keeping an accurate prediction of the state
distribution even for a long period of time, even if slightly reducing the
accuracy. The figure has a logarithmic scale for the values of moments,
highlighting the diverging nature of the CRTBP dynamics. Thus, it is worth
noticing the exponential increase of the spread of the state PDF, with
standard deviations that are two order of magnitude bigger after merely three
quarters of a full revolution.

Figure 3: Prediction of central moments of the state PDF for a Halo orbit
orbiting around $\mathcal{L}_{1}$ in the earth-Moon system.

This application assesses the correct prediction of uncertainties and opens
the road to the KOF, a complete filtering algorithm, that completes the
propagation of mean and covariance with the knowledge obtained externally
though noisy observations.

###  5.3 Orbit Determination

While the previous problem focused on the corrected prediction of the moments
of the state PDF under an highly nonlinear dynamics, a station keeping
application is now proposed, where the system receives external information
from a set of sensors, and the whole KOF algorithm is implemented. Therefore,
let consider the Sun-Earth system, characterized by
$\mu=3.003410642560030e-06$ and $\gamma=0.009970325504020$, and a Lyapunov
orbit around the $\mathcal{L}_{1}$ libration point. The initial true state is
selected randomly from a Gaussian distribution with initial state mean

$\displaystyle{\bf{x}}_{0}=\left[0.989826595322\quad 0\quad 0\quad 0\quad
0.00137295958289\quad 0\right]^{T}$ (55)

and covariance matrix ${\bf P}_{0}=\sigma^{2}{\bf I}_{6\times 6}$ with
$\sigma=10^{-4}$. The KOF scales and translates uncertainties and central
moments to be consistent with the normal form representation of the
Hamiltonian $H_{3}$. The propagation of the truth, used to evaluate the
performance of the filter in terms of accuracy error, is propagated
numerically using and high precision Runge-Kutta 7/8 integrator. Accuracy
errors and estimates are then evaluated by scaling back the results obtained
from the KOF, in order to agree with the classic representation of the
dynamics.

Nonlinear angle measurements, azimuth and elevation from the center of the
Earth, are received by the filter every 0.4 year, which translates to about a
total of 15 observations each 2 revolutions;

$\displaystyle y_{1}$
$\displaystyle=\arctan\left(\dfrac{x_{2}}{x_{1}-1+\mu}\right)+\eta_{1}$ (56)
$\displaystyle y_{2}$
$\displaystyle=\arcsin\left(\dfrac{x_{3}}{\sqrt{(x_{1}-1+\mu)^{2}+x_{2}^{2}+x_{3}^{2}}}\right)+\eta_{2}$
(57)

where each $\eta$ is assumed Gaussian noise with zero mean and a standard
deviation of 10 arcsec.

A Monte Carlo analysis is performed to assess the performance of the filter.
The initial level of uncertainties is so large that the accuracy level at
steady state is two orders of magnitude smaller than the initial standard
deviation. Therefore, in order to show convergence in a more pleasant manner,
a logarithmic scale for the errors has been applied. Thus, Figure 4 shows, for
each component of the state, the 500 runs Monte Carlo analysis performed with
the KOF on the selected Lyapunov orbit. Indeed, each gray line represents the
errors connected to the $j$-th run of the Monte Carlo, for each $i$-th state
component, evaluated as

$\displaystyle\epsilon_{i,j}$ $\displaystyle=\hat{x}_{i,j}-x_{i,j}$ (58)

The mean of the errors is portrayed with black lines, calculated as

$\displaystyle\bar{\epsilon}_{i}$
$\displaystyle=\dfrac{1}{N}\sum_{j=1}^{N}\epsilon_{i,j}$ (59)

where $N$ is the number of total Monte Carlo runs. The black lines settle on
the zero value, proving that the KOF is an unbiased filter, as expected from
the Minimum Mean Square Error (MMSE) principle on which the update of the
filter is based on. The Monte Carlo analysis provides also information about
the spread of the state error and the level of uncertainties. Figure 4 has
multiple pairs of dashed and continuous blue lines. Dashed blue lines
represents the level of standard deviation, as $\pm 3\sigma_{eff}$, obtained
from the Monte Carlo runs and calculated using Equation (52): they represent
the effective level of accuracy of the filter. These dashed lines are
calculated, for each time step, by considering the values of the state errors
among all the 500 simulations and they represent the actual behaviour of the
filter. On the contrary, the continuous blue lines report the estimated level
of uncertainties, again as $\pm 3\sigma_{pred}$, predicted by the filter.
Thus, for each time step, these lines are calculated extracting the error std
directly from the update covariance of the filter, as
$\sigma_{pred,i}=\sqrt{{\bf P}^{+}_{ii}}$. By looking at Figure 4, we can
assess the consistency of the KOF, since the dashed and continuous blue lines
overlap, meaning that the filter is able to correctly predict its own
uncertainty. The KOF shows convergence and consistency, where steady state
accuracy levels are reached rapidly after a few updates. The large initial
uncertainties are rapidly reduced around the most current estimate, and the
updated covariance of the filter correctly represents the spread of the state
error PDF.

Figure 4: 500 runs Monte Carlo analysis of the KOF for a Lyapunov orbit around
$\mathcal{L}_{1}$.

The performance of the KOF can now be compared to other common filters: the
EKF, the iterated extended Kalman filter (IKF), and the UKF have been selected
as benchmark. The comparison in accuracy among different estimators is
performed by analyzing the updated level of uncertainties, with a consistency
check obtained from the Monte Carlo. Therefore, the 500 runs Monte Carlo
analysis has been repeated for the EKF, IKF and UKF. Standard deviations for
position and velocity are calculated and compared separately. The predicted
stds are obtained directly from the trace of the updated covariance matrix,
for each time step:

$\displaystyle\sigma_{pos,pred}=\sqrt{{\bf P}^{+}_{11}+{\bf P}^{+}_{22}+{\bf
P}^{+}_{33}}$ (60) $\displaystyle\sigma_{vel,pred}=\sqrt{{\bf P}^{+}_{44}+{\bf
P}^{+}_{55}+{\bf P}^{+}_{66}}$ (61)

while the effective stds levels, $\sigma_{pos,eff}$ and $\sigma_{vel,eff}$,
are evaluated directly from the Monte Carlo runs. Thus, a consistent filtering
algorithm has matching between $\sigma_{eff}$ and $\sigma_{pred}$. Figure 5
reports the accuracy levels, in terms of error stds, of the four filters,
displaying estimated values with continuous lines and the actual behaviour
with dashed lines. The EKF is shown in red. The dashed red lines, both for
position and velocity, are out of scale when compared to the continuous
counterparts, which are overlapping with the green continuous lines. Indeed,
the green lines are connected to the IKF. The EKF and the IKF share the same
updated covariance matrix, while their estimate changes. Thus, the overlapping
between the green and the red continuous lines matches with what is expected
from theory: the two estimators approximate the posterior distribution of the
state as Gaussian, with the same level of uncertainty (same stds), but
different means. Indeed, while the EKF is a MMSE filter, the IKF is a filter
based on the Maximum A Posteriori (MAP) principle, which outputs, as its
estimate, the most likely state of the posterior distribution. However, the
IFK and EKF share the same prediction step, based on the linearization of the
dynamics, which is not sufficient to achieve an accurate approximation of the
state prior distribution. Therefore, both filters are inconsistent and their
dashed lines settle orders of magnitude above the continuous ones. The EKF and
IKF believe that they are achieving higher accuracy than their actual results,
and they are overconfident in their performance The difference between the
green dashed lines and the red dashed lines is connected to the different
update step performed by the two filtering techniques: a MMMSE one for the
EKF, and a MAP one for the EKF. On the contrary, in blue, the UKF applies the
unscented transformation in its prediction step to obtain a more accurate
state prior distribution. Thus, the dashed blue line, which shows the
effective accuracy level of the UKF, settles below the EKF and the IKF ones.
However, the matching with the predicted uncertainty is missing and the filter
is inconsistent as well. The UKF estimates a performance similar to the KOF,
in black, shown by the overlapping between the black and the blue continuous
lines. Therefore, the KOF is the only filter that converges with consistency.
The KOF is able to correctly predict its own uncertainties and its estimate is
the most accurate among the other filtering techniques. The dashed black lines
are the lowest dashed lines for both position and velocity, and they overlap
their continuous counterparts along the whole time-length of the simulation.

Figure 5: Position and Velocity standard deviation analysis among EKF, IKF,
UKF, and KOF.

##  6 Conclusions

The propagation of uncertainties under highly nonlinear dynamics has been
performed in the KO framework, where the solution of the system is
approximated using a set of well defined orthogonal polynomials. Thanks to
this representation, central moments of the state PDF can be predicted up to
an arbitrary order. The new methodology has been paired with the Kalman
measurement update to develop the KOF: a nonlinear filtering algorithm where
both the state and measurement distributions are represented using the KO
eigenfunctions. Numerical applications have assessed the robustness and
accuracy of the new theory. Uncertainty prediction in the CRTBP has shown that
the new technique is able to track the exponential spread of the spacecraft
position and velocity PDF, providing estimates about the system’s raw and
central moments. Whenever information is provided in terms of external
observations, the KOF is able to shrink the state uncertainties around the
most current estimate. Monte Carlo analyses have assessed the accuracy,
consistency, and robustness of the KOF, which has been proven to outperform
other filtering techniques reported for comparison purposes.

## Acknowledgments

The authors want to acknowledge the support of this work by the Air Force’s
Office of Scientific Research under Contract Number FA9550-22-1-0092.

## References

  * [1] Y. Bar-Shalom, X. R. Li, and T. Kirubarajan, Estimation with applications to tracking and navigation: theory algorithms and software.  John Wiley & Sons, 2004.
  * [2] R. E. Kalman, “A new approach to linear filtering and prediction problems,” 1960\.
  * [3] R. E. Kalman and R. S. Bucy, “New results in linear filtering and prediction theory,” 1961.
  * [4] D. Simon, Optimal state estimation: Kalman, H infinity, and nonlinear approaches.  John Wiley & Sons, 2006.
  * [5] B. Schutz, B. Tapley, and G. H. Born, Statistical orbit determination.  Elsevier, 2004.
  * [6] A. Gelb et al., Applied optimal estimation.  MIT press, 1974.
  * [7] J. L. Junkins and P. Singla, “How nonlinear is it? A tutorial on nonlinearity of orbit and attitude dynamics,” 2004.
  * [8] S. Julier, J. Uhlmann, and H. F. Durrant-Whyte, “A new method for the nonlinear transformation of means and covariances in filters and estimators,” IEEE Transactions on automatic control, Vol. 45, No. 3, 2000, pp. 477–482.
  * [9] S. J. Julier and J. K. Uhlmann, “Unscented filtering and nonlinear estimation,” Proceedings of the IEEE, Vol. 92, No. 3, 2004, pp. 401–422.
  * [10] R. S. Park and D. J. Scheeres, “Nonlinear mapping of Gaussian statistics: theory and applications to spacecraft trajectory design,” Journal of guidance, Control, and Dynamics, Vol. 29, No. 6, 2006, pp. 1367–1375.
  * [11] R. S. Park and D. J. Scheeres, “Nonlinear semi-analytic methods for trajectory estimation,” Journal of guidance, control, and dynamics, Vol. 30, No. 6, 2007, pp. 1668–1676.
  * [12] M. Majji, J. L. Junkins, and J. D. Turner, “A high order method for estimation of dynamic systems,” The Journal of the Astronautical Sciences, Vol. 56, No. 3, 2008, pp. 401–440.
  * [13] M. Valli, R. Armellin, P. Di Lizia, and M. Lavagna, “Nonlinear mapping of uncertainties in celestial mechanics,” Journal of Guidance, Control, and Dynamics, Vol. 36, No. 1, 2013, pp. 48–63.
  * [14] S. Servadio, F. Cavenago, P. Di Lizia, and M. Massari, “Nonlinear Prediction in Marker-Based Spacecraft Pose Estimation with Polynomial Transition Maps,” Journal of Spacecraft and Rockets, 2021, pp. 1–13.
  * [15] F. Cavenago, P. Di Lizia, M. Massari, S. Servadio, and A. Wittig, “DA-based nonlinear filters for spacecraft relative state estimation,” 2018 Space Flight Mechanics Meeting, 2018, p. 1964.
  * [16] B. O. Koopman, “Hamiltonian systems and transformation in Hilbert space,” Proceedings of the National Academy of Sciences, Vol. 17, No. 5, 1931, pp. 315–318, 10.1073/pnas.17.5.315.
  * [17] J. v. Neumann, “Zur Operatorenmethode in der klassischen Mechanik,” Annals of Mathematics, 1932, pp. 587–642, 10.2307/1968537.
  * [18] A. Naylor and G. Sell, Linear Operator Theory in Engineering and Science.  Applied Mathematical Sciences, Springer, 2000.
  * [19] I. Prigogine, Non–equilibrium statistical mechanics.  Dover Publications, 1962.
  * [20] I. Mezić, “Analysis of fluid flows via spectral properties of the Koopman operator,” Annual Review of Fluid Mechanics, Vol. 45, 2013, pp. 357–378, 10.1146/annurev-fluid-011212-140652.
  * [21] S. L. Brunton, B. W. Brunton, J. L. Proctor, and J. N. Kutz, “Koopman invariant subspaces and finite linear representations of nonlinear dynamical systems for control,” PloS one, Vol. 11, No. 2, 2016, p. e0150171, 10.1371/journal.pone.0150171.
  * [22] A. Surana, “Koopman operator based observer synthesis for control-affine nonlinear systems,” Decision and Control (CDC), 2016 IEEE 55th Conference on, IEEE, 2016, pp. 6492–6499.
  * [23] T. Chen and J. Shan, “Koopman-operator-based attitude dynamics and control on SO (3),” Journal of Guidance, Control, and Dynamics, Vol. 43, No. 11, 2020, pp. 2112–2126.
  * [24] R. Linares, “Koopman Operator Theory Applied to the Motion of Satellites,” Advances in the Astronautical Sciences, Vol. 171, AAS/AIAA, 2019, pp. AAS 19–821.
  * [25] D. Arnas, R. Linares, and K. Alfriend, “An Analysis of Koopman-based Perturbation Theory Applied to the Motion About an Oblate Planet,” 31st AAS/AIAA Space Flight Mechanics Meeting, 2021.
  * [26] D. Arnas and R. Linares, “Koopman operator theory applied to the J2 perturbed motion around an oblate planet,” 31st AAS/AIAA Space Flight Mechanics Meeting, AAS/AIAA, 2021.
  * [27] D. Arnas and R. Linares, “Approximate analytical solution to the zonal harmonics problem using Koopman operator theory,” Journal of Guidance, Control, and Dynamics, Vol. 44, No. 11, 2021, pp. 1909–1923.
  * [28] A. Wittig, P. Di Lizia, R. Armellin, K. Makino, F. Bernelli-Zazzera, and M. Berz, “Propagation of large uncertainty sets in orbital dynamics by automatic domain splitting,” Celestial Mechanics and Dynamical Astronomy, Vol. 122, No. 3, 2015, pp. 239–261.
  * [29] S. Servadio and R. Zanetti, “Differential Algebra-Based Multiple Gaussian Particle Filter for Orbit Determination,” Journal of Optimization Theory and Applications, Vol. 191, No. 2, 2021, pp. 459–485.
  * [30] S. Servadio, D. Arnas, and R. Linares, “Dynamics Near the Three-Body Libration Points via Koopman Operator Theory,” arXiv preprint arXiv:2110.12119, 2021\.
  * [31] D. Arnas and R. Linares, “Approximate Analytical Solution to the Zonal Harmonics Problem Using Koopman Operator Theory,” Journal of Guidance, Control, and Dynamics, 2021, accepted for publication, 10.2514/1.G005864.
  * [32] S. Servadio, D. Arnas, and R. Linares, “A Koopman Operator Tutorial with Othogonal Polynomials,” arXiv preprint arXiv:2111.07485, 2021.
  * [33] S. Servadio and R. Zanetti, “Recursive polynomial minimum mean-square error estimation with applications to orbit determination,” Journal of Guidance, Control, and Dynamics, Vol. 43, No. 5, 2020, pp. 939–954.
  * [34] S. Servadio, R. Zanetti, and B. A. Jones, “Nonlinear filtering with a polynomial series of Gaussian random variables,” IEEE Transactions on Aerospace and Electronic Systems, Vol. 57, No. 1, 2020, pp. 647–658.
  * [35] L. Isserlis, “On a formula for the product-moment coefficient of any order of a normal frequency distribution in any number of variables,” Biometrika, Vol. 12, No. 1/2, 1918, pp. 134–139.
  * [36] S. Servadio and R. Zanetti, “Estimation Of The Conditional State And Covariance with Taylor Polynomials,” Journal of Advances in Information Fusion, Vol. 6, No. 2, 2021, pp. 126–142.
  * [37] D. L. Richardson, “Analytic construction of periodic orbits about the collinear points,” Celestial mechanics, Vol. 22, No. 3, 1980, pp. 241–253.
  * [38] W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, “Dynamical systems, the three-body problem and space mission design,” Equadiff 99: (In 2 Volumes), pp. 1167–1181, World Scientific, 2000.
  * [39] A. Jorba and J. Masdemont, “Dynamics in the center manifold of the collinear points of the restricted three body problem,” Physica D: Nonlinear Phenomena, Vol. 132, No. 1-2, 1999, pp. 189–213.


# An analysis of over-sampling labeled data in semi-supervised learning with
FixMatch

Miquel Martí i Rabadán Corresponding Author<EMAIL_ADDRESS>KTH Royal
Institute of Technology, Sweden  Univrses AB, Sweden  Sebastian Bujwid  KTH
Royal Institute of Technology, Sweden  Alessandro Pieropan  Univrses AB,
Sweden  Hossein Azizpour  KTH Royal Institute of Technology, Sweden  Atsuto
Maki  KTH Royal Institute of Technology, Sweden

###### Abstract

Most semi-supervised learning methods over-sample labeled data when
constructing training mini-batches. This paper studies whether this common
practice improves learning and how. We compare it to an alternative setting
where each mini-batch is uniformly sampled from all the training data, labeled
or not, which greatly reduces direct supervision from true labels in typical
low-label regimes. However, this simpler setting can also be seen as more
general and even necessary in multi-task problems where over-sampling labeled
data would become intractable. Our experiments on semi-supervised CIFAR-10
image classification using FixMatch show a performance drop when using the
uniform sampling approach which diminishes when the amount of labeled data or
the training time increases. Further, we analyse the training dynamics to
understand how over-sampling of labeled data compares to uniform sampling. Our
main finding is that over-sampling is especially beneficial early in training
but gets less important in the later stages when more pseudo-labels become
correct. Nevertheless, we also find that keeping some true labels remains
important to avoid the accumulation of confirmation errors from incorrect
pseudo-labels.

##  1 Introduction

Semi-supervised learning has recently established itself as an effective
approach to greatly reduce the annotation costs in developing deep learning
models. With labels only for a small portion of the entire dataset, the
performance approaches that of fully-supervised learning, e.g. giving
competitive results in the context of image classification using less than 1%
of labeled samples [3, 22, 18, 23].

Semi-supervised methods have explored different flavours of consistency or
perturbation-based regularization [21] combined with entropy minimization [6]
or pseudo-labeling [9] by adding one or more extra terms to the training loss
computed on the unlabeled samples.

The consistency-based regularization losses aim at exploiting the smoothness
assumption [21], i.e. the predictions of a model should not depend on small
perturbations in the input such as injecting noise [16, 2, 8], adversarial
perturbations [11] or random transformations [17, 8, 3, 22, 18]. Pseudo-
labeling and entropy minimization methods enforce instead high-confidence
(low-entropy) predictions for unlabeled samples, directly exploiting the low-
density assumption [21], i.e. decision boundaries should lie on low-density
regions. These approaches are closely related to self-training [13, 25, 1,
23], where pseudo-labels are extracted from a model pre-trained on labeled
data and then used for training in a second step or in multiple iterative
steps.

|
---|---
(a) Explicit setting | (b) Implicit setting
Figure 1: (a) and (b) A toy dataset with 20 samples and a task with 2 classes
- dog (green) and cat (red) - with 2 labeled samples per class, and a ratio of
labeled samples of $20\%$. Sampled into 5 training mini-batches of 4 samples
each, with over-sampling of labeled data in (a) and uniform sampling in (b).
In (a) we choose to have an equal number of labeled and unlabeled samples in
each mini-batch. Due to this, each labeled sample (red/green, #0-3) is seen 4x
more often than each unlabeled sample (gray, #4-19). In (b), each mini-batch
can contain different amounts of labeled samples as it is uniformly sampled
regardless of labels and each sample is seen exactly once. Labeled samples are
seen 2.5x more often in the explicit setting.

Yet, all semi-supervised learning methods face a common issue. Since they
bootstrap their own predictions in order to use them as targets, many of these
will be incorrect and can lead to so-called confirmation bias [20, 1], i.e.
fitting incorrect labels that do not represent the true task and thus slowing
down or stopping its learning. This is more likely to happen at the start of
training since the predictions of the model are not much better than random
and the amount of misclassified samples is large.

Most works try to alleviate this issue by using a ramp-up from zero of the
weight for the semi-supervised loss term over training. This greatly reduces
the influence of possibly wrong predictions on unlabeled data in the early
phases of training. Thus, the only way to learn at the start is from true
labels, but since these are scarce there is little supervision available and
learning is slow. This in turn makes it more difficult to get useful
supervision through the semi-supervised loss.

Repeating labeled samples more often than unlabeled ones during training is a
common mini-batch sampling strategy in semi-supervised learning and directly
and artificially increases the influence of labeled data over learning [16,
14, 3, 20, 11, 1, 22, 18, 23]. To the best of our knowledge, this mechanism
was never well studied with the goal of enhancing the training signal from
true labels and it likely just happened as a side effect of common
implementation choices when constructing batches from separate datasets for
labeled and unlabeled samples.

In this paper, we aim at studying if it is essential to repeat labeled samples
more frequently than unlabeled samples or over-sample labeled data for state-
of-the-art semi-supervised methods to work well. We call this way of sampling
mini-batches the explicit semi-supervised learning setting due to its
distinction between labeled and unlabeled data in how the data is presented.
We compare it to an alternative setting in which all data is sampled
uniformly, regardless of being labeled or not. We call the latter the implicit
semi-supervised learning setting.111In the rest of the paper we will use both
nomenclatures for the sampling approaches interchangeably. Figure 1
illustrates both approaches.

Exploring the simpler implicit setting will highlight the benefits which the
common over-sampling of labeled data in the explicit setting brings and which
mechanisms allow to achieve these. Understanding the differences between both
settings could then shed light on how sampling affects the learning process in
semi-supervised methods and could ultimately guide future work into developing
more efficient and better-performing methods.

Moreover, there are concrete reasons to avoid the explicit setting. Most
importantly, trying to use it for multi-task problems easily leads into a
roadblock. Samples cannot be said to be labeled or unlabeled in general, since
this property can only be assigned to a sample in relation to each task. Thus,
samples can be labeled for one or more tasks but unlabeled for others, labeled
for no tasks or for all. This leads to different sets of samples labeled for
each task, which might or might not overlap.

|
---|---
(a) Multi-task extension of explicit setting | (b) Implicit setting
Figure 2: (a) and (b) The labeling of the dataset in Figure 1 is extended with
a new task with 2 classes - black and white - with 2 labeled samples per class
as well. In (a) we extend the over-sampling of labeled data to the multi-task
case so that in each batch there is the same number of labeled and unlabeled
samples for each task. In (b) the uniform sampling means the sampling is
independent of the labels, just as for the single task case.

In Figure 2, we compare an extension of the explicit setting for multi-task
scenarios to the implicit setting. One possible direct way to recreate the
explicit setting for multi-task learning is to divide each mini-batch into
multiple sub-parts, one for each possible configuration of available task
labels. Such a mini-batch can include parts where all samples are labeled or
unlabeled, parts with samples labeled only for one or more tasks, etc.
Instead, the implicit setting can be used in multi-task problems without any
modifications since the sampling is independent of the labels.

The number of sub-parts in a mini-batch constructed using this extension of
the explicit setting grows exponentially with the number of tasks $T$ as
$2^{T}$ if all label configurations are present in the dataset. Therefore, if
the number of tasks is large, the number of mini-batch parts can quickly
become larger than the total batch size that can be used in practice. Even
with fewer tasks, the possible size of each mini-batch part becomes a new
hyper-parameter which can be difficult to control since it is limited by the
total batch size and the number of parts. Setting the number of tasks to one
recovers the original explicit setting, with one part for labeled data for the
only task in the problem and one for unlabeled data. Other sampling approaches
for semi-supervised multi-task learning problems are likely possible, but our
aim is to highlight the issues with the explicit setting that make it not
suitable for this class of problems.

To study the impact of different approaches to sampling training mini-batches
on semi-supervised performance, we compare the training characteristics of
both the implicit and explicit settings. We run semi-supervised image
classification experiments on CIFAR-10 and CIFAR-100 [7], varying the number
of available labels. As a semi-supervised method we adopt FixMatch [18], a
standard, state-of-the-art method. FixMatch enforces consistency between
pseudo-labels obtained from weakly transformed versions of one sample and the
predictions obtained from strongly transformed versions of the same sample,
incorporating in the same method pseudo-labeling and consistency
regularization. We find that both the implicit and explicit settings work for
reasonably long training runs, but the latter performs better. This difference
is smaller in higher-label regimes or with longer training budgets and we
partially attribute it to the low efficiency of learning from true labels at
the start of training when using the implicit setting, which becomes more
important for shorter training runs.

In sum, this paper (i) emphasizes the importance of the choice of mini-batch
sampling strategy in semi-supervised learning, (ii) analyses the differences
between the common over-sampling of labeled data and standard uniform sampling
throughout the training process using FixMatch, and (iii) demonstrates the
viability of uniform sampling with long training budgets, enabling its use in
cases where over-sampling labeled data is not possible, such as in multi-task
learning.

##  2 Related work

Most semi-supervised learning methods use some mechanism that minimizes the
effect of confirmation bias, either at the earlier stages of training or
throughout. One such mechanism directly ramps up the weight for the
unsupervised loss term during training, starting from zero. It was originally
proposed in Grandvalet and Bengio [6] and is reported to be used in most works
[9, 2, 8, 20, 3]. FixMatch does not use the weight ramp-up as it instead sets
a threshold to use only confident predictions as pseudo-labels and averages
the loss on these over the total number of unlabeled samples instead of over
the number of pseudo-labels used. This achieves a similar effect to the ramp-
up since it decreases the contribution of the unsupervised loss when the model
produces low confidence predictions early in the training.

A second mechanism is to over-sample labeled data, i.e. repeating labeled
samples more often than unlabeled ones when constructing mini-batches. This is
usually implemented by splitting the labeled and unlabeled samples as two
separate datasets and sampling the labeled and unlabeled parts of the mini-
batch from each. Then, the first set that gets depleted is repeated as many
times as needed until the other one is used up. Since the labeled dataset is
typically much smaller, this results in the repetition of labeled samples at a
higher rate. This approach was already used in Rasmus et al. [16]222See:
https://github.com/CuriousAI/ladder/blob/
5a8daa1760535ec4aa25c20c531e1cc31c76d911/run.py#L74 although without any
explicit mention in the paper. Most recent works in semi-supervised learning
use this setting, both using consistency-based, pseudo-labeling or self-
training methods, but did not study the impact of sampling mini-batches in
this way. Some methods set equal sizes for the labeled and unlabeled parts of
each training mini-batch [16, 14, 3], but some use a larger unlabeled batch
size [20, 11, 1, 22, 18, 23].

Both FixMatch [18] and Noisy Student [23] include ablation experiments on this
ratio and conclude that a larger ratio of unlabeled samples in each training
mini-batch leads to better performance. Arazo et al. [1] is the only work we
are aware of that studied the effect of keeping a minimum amount of true
labels in each mini-batch in pseudo-labeling, showing that it is indeed an
effective way of reducing confirmation bias. Only Temporal Ensembling [8] and
one of the experiments in Mean Teacher [20] use the implicit setting. However,
the latter mentions that repeating labeled samples can be beneficial since
”the supervised training signal is strong enough early on to train quickly and
prevent getting stuck into uncertainty” [5]. All methods relying on the
explicit setting can technically be adapted to use the implicit setting
instead.

##  3 Explicit vs. implicit setting

We aim at understanding how essential the explicit setting is for good
performance in semi-supervised learning, and whether the implicit one can
achieve competitive results. To do so, we will compare the final
classification performance of both settings.

###  3.1 Experimental setting

We compare the explicit and implicit settings head to head when using FixMatch
[18], which uses asymmetric perturbations in the input space in a siamese
network configuration with two branches. The teacher branch gets weakly
augmented samples (only flip-and-shift operations) and generates pseudo-labels
from them, filtering out those below a certain confidence threshold. The
student branch gets strongly augmented samples via RandAugment [4] and is
trained using the standard cross-entropy loss against both true labels (if
available) and the pseudo-labels generated by the teacher.

We use the semi-supervised image classification problem on the CIFAR-10 [7]
dataset for the main experiments, as it is a standard benchmark in the deep
semi-supervised learning literature [14]. We use three different amounts of
labeled samples - 40, 250 and 4000 - simulating respectively a problem with
extremely limited labeled data at only four labels per class, a scarcely-
labeled problem with approximately 0.5% labeled samples, and a semi-supervised
learning problem with a larger amount of labeled samples with approximately 9%
labeled samples. We reserve 10% of the total samples in the training set for
validation and use the rest as training data, which we split into labeled and
unlabeled subsets. All splits are done preserving the class balance. The
validation set is used for early stopping and for hyper-parameter tuning based
on random search with a budget of 100 trials per experiment. We pick the
hyper-parameters leading to the highest validation accuracy on a single seed
and data split. We tune the following hyper-parameters for each case: learning
rate, batch size, weight of the semi-supervised loss, and the ratio of labeled
data in each batch for the explicit setting,. Finally, we report the final
classification accuracy values on the test set as the mean and standard
deviation over five different random seeds and five different splits of the
labeled samples in the training set.

In addition, we perform experiments on CIFAR-100 with 2500 and 10000 labels,
following the same setting as for CIFAR-10 but not tuning the hyper-
parameters. Instead, we use the ones reported in FixMatch [18] for the
explicit setting and our best single guess of a decent set of hyper-parameters
for the implicit setting.

For CIFAR-10 experiments, we use a Wide ResNet-28-2 [24] (WRN-28-2) as the
base network with a linear classifier after the pooling layer and train it via
stochastic gradient descent with Nesterov momentum of $0.9$ [15, 12, 19], a
cosine learning rate decay schedule [10], and weight decay of $0.0005$. For
CIFAR-100, we use WRN-28-5 and a smaller weight decay of $0.0001$. Evaluation
is done on an exponential moving average of the model parameters with decay
$0.999$ [18], which we found gave more stable results and significantly better
performance early in the training.

Each training run has a maximum budget for the total number of samples that
each model sees during training. We set this budget to be equivalent to $1000$
epochs on the full training set, i.e. $45$ million samples, after removing
5000 samples from the training set for validation. Therefore we refer to an
epoch as one pass through as many samples as the size of the training set,
regardless of whether each sample has an associated label or not, or how many
times it has been repeated.

Other works [14, 18] have used a fixed number of labeled samples or a fixed
number of training steps instead. We believe this approach leads to unfair
comparisons between supervised baselines or when using different label ratios
in each batch since for higher ratios the models will effectively see many
more samples over training. Moreover, comparing the total amount of samples
seen aligns better with the amount of compute used for each training run.
Using this definition of training length, MixMatch [3] training runs go
through $\sim 134$ million samples since they use a total batch size of 128
(64 for labeled samples and 64 for unlabeled samples) but they only count the
labeled ones towards their definition of total length. FixMatch, using an
unlabeled batch size seven times larger, uses a total batch size of $512$ and
the total number of samples seen during training is $\sim 536$ million.

We choose to use a smaller budget on purpose to better reflect what would
happen in practical applications and to keep contained the computational
budget required for our experiments. We noticed that for smaller amounts of
labeled data longer training budgets are important, so we include longer
experiments in such cases as well.

###  3.2 Experimental results

CIFAR-10 40 labels 250 labels 4000 labels All labels Supervised $36.77\pm
4.48$ $59.88\pm.73$ $87.39\pm.20$ $96.54\pm.11$ FixMatch(E) $83.44\pm 6.68$
$93.01\pm.58$ $94.89\pm.16$ - FixMatch(I) $73.90\pm 8.04$ $91.42\pm.98$
$94.84\pm.05$ - FixMatch(E) 6x $85.40\pm 2.83$ $94.40\pm.78$ - - FixMatch(I)
6x $86.59\pm 4.14$ $94.00\pm.54$ - - CIFAR-100 2500 labels 10000 labels All
labels Supervised $47.26$ $67.58$ $82.38$ FixMatch(E) $62.86$ $74.61$ -
FixMatch(I) $55.53$ $71.83$ -

Table 1: Classification accuracy on test data for supervised baselines and
FixMatch with both implicit (I) and explicit (E) settings using different
amounts of labels. Shown as mean and standard deviation on 5 runs with
different seeds and data splits for each CIFAR-10 experiment and a single seed
and data split for CIFAR-100. 6x indicates six times longer training budget.

Table 1 shows the classification accuracy of our models trained on CIFAR-10
and CIFAR-100 with different amounts of labels when using FixMatch in either
the explicit or implicit settings.

We observe that both settings outperform supervised baselines by a significant
margin in CIFAR-10. For 4000 labels both settings perform similarly, however
we see a larger gap favouring the explicit setting by $1.6$ percentage points
when only using 250 labels. A significantly larger gap of $9.5$ percentage
points appears when using 40 labels, although both settings show a much larger
standard deviation. Comparing individual runs on each data split for both
settings, we can see that the large standard deviations are caused mainly by
the choice of labeled samples in each data split but also that the implicit
setting is more sensitive to it.

We hypothesize that this gap is due to the much smaller amount of labeled
samples present (on average) in each training mini-batch for the lower-label
regimes in the implicit setting, where for example around only one labeled
sample appears for every 200 unlabeled samples for the 250-label case. Due to
this lack of direct supervision at the start, the model learns at a very slow
pace. This is because the model is initially relying only on labeled data as
the predictions on unlabeled data are below the confidence threshold. A longer
training budget of 6000 epochs improves the performance of both settings. The
gap in test accuracy between them reduces over training, and both settings end
up performing on par.

We see similar results for CIFAR-100. The gap between settings is larger at
lower label regimes, reaching $7.3$ percentage points better accuracy for the
explicit setting with 2500 labels and $2.8$ percentage points with 10000
labels, and both perform better than the supervised baselines.

|  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 3: Training curves for CIFAR-10 with 250 labels (red) or 4000 labels
(blue) using FixMatch in the explicit (solid line) or implicit (dashed)
settings. (a) Top: training error on labeled samples. Bottom: test error on
the unseen test set. Values at the end of each epoch. (b) Top: evolution over
training of cross-entropy loss on labeled training data. Bottom: semi-
supervised cross-entropy loss between pseudo-labels and predictions on
unlabeled data. (c) Top: average confidence of predictions on unlabeled
samples in each batch. Bottom: ratio of predictions used as pseudo-labels. (d)
Top: Ratio of correct predictions on unlabeled data. Bottom: Ratio of correct
pseudo-labels on each batch.

##  4 Analysis

Training dynamics We analyze in Figure 3 the training curves of models in both
settings for 250 and 4000 labels to illustrate the impact of their different
training dynamics. We can see that the explicit setting does indeed fit
quicker the labeled data; training error, shown in (a, top), and supervised
loss, in (b, top), both decrease quickly thanks to more direct supervision
from repeating labeled samples. The models start giving confident predictions
earlier, as shown in (c, top), and thus the number of pseudo-labels, in (c,
bottom), and the semi-supervised loss, in (b, bottom), grow quicker. On the
other hand, for the implicit setting, the model takes much longer to fit the
training data and to produce high-confidence predictions that will be used as
pseudo-labels. Therefore, the semi-supervised loss grows much slower, having
no effect for the first 100 epochs while learning is very slow and test
accuracy is barely better than random. This clearly points to a lack of direct
supervision from true labels early in training as one cause for the implicit
setting lagging behind in terms of final test error with our relatively short
training budget with fewer labeled samples.

Moreover, we make use of privileged information (the labels dropped to
simulate unlabeled samples) to compute and show in Figure 3 (d) the ratio of
correct predictions on unlabeled samples and the ratio of the correct pseudo-
labels, i.e. the subset of predictions on unlabeled data whose confidence was
above the threshold of $0.95$. We see that not only are pseudo-labels used
earlier in the explicit setting but are also cleaner, i.e. more of them are
correct, from early on and can thus help guide learning of the true
classification task better. Also, we can see that the usage of confident
predictions on unlabeled data as pseudo-labels is effective since the subset
of chosen pseudo-labels has a ratio of correctness consistently higher than
that of all predictions on unlabeled data.

This analysis holds for both amounts of labeled samples shown, although the
differences are much smaller for the 4000-label experiments. In this last
case, the implicit setting has a lower test error earlier during training
which we hypothesize might be due to the larger over-fitting to labeled
samples in the explicit case, before the loss on pseudo-labels starts having a
regularizing effect and the model can generalize to unseen data.

Training without direct supervision  We hypothesize that supervision from true
labels loses importance when the model can already generate cleaner pseudo-
labels, which could explain why either using more labels or longer training
budgets reduces the performance gap between settings.

To verify this hypothesis, we run FixMatch with the semi-supervised loss only
but starting training from different supervised checkpoints. All models
improve their test accuracy at first - the best goes from $87.39\%$ to
$92.9\%$ \- but their accuracy ends up dropping suddenly to levels close to
random while fitting the pseudo-labels. This can be attributed to the
accumulation of confirmation errors and the lack of supervision from true
labels to correct for these. We see this happen as well for long experiments
with the implicit setting and few labels, but with learning recovering after
the dip. Thus, we conclude that seeing some labeled samples is still required
to guide learning to the true task even when few pseudo-labels are incorrect.

##  5 Conclusion

Our experiments in semi-supervised image classification show that using the
implicit setting incurs a penalty in final accuracy compared to using the
explicit setting. How large this gap becomes depends, however, on the amount
of labels and training budget. Both settings still outperform supervised
baselines by a large margin. The explicit setting is significantly more
efficient in the beginning but the difference between both settings is reduced
over training and both perform on par when more labeled data is available and
for longer training runs. However, seeing labeled samples remains important
even in the later stages of training, as it avoids falling into confirmation
bias. The more general aspects like the required amount of labels and the most
effective way to use them are yet to be studied. Overall, these are important
results since the implicit setting is much simpler to use, potentially
allowing semi-supervised learning on settings like partially-labeled multi-
task problems which would not be possible in the conventional explicit
setting.

This analysis was done in the context of FixMatch, which we believe is a good
representative of modern semi-supervised deep learning methods combining
consistency regularization and pseudo-labeling approaches. Future research
will be directed to see if the results also apply to semi-supervised methods
that are substantially different.

Finally, we suggest two possible alternative approaches to mini-batch sampling
in the multi-task scenario. First, if there is a large enough amount of
samples labeled for all tasks, a version of the explicit setting could be
constructed using fully labeled samples versus partially labeled samples.
Second, we could construct a variation of the implicit setting with a
guaranteed minimum number of samples labeled for each task so that there is
always some direct supervision from true labels for each task, in each mini-
batch. Studying these and other possible sampling methods will be an
interesting direction for future work.

## Acknowledgments

This work was partially supported by the Wallenberg AI, Autonomous Systems and
Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation.

The support by the Swedish Research Council through grant agreement no.
2016-04022 is also gratefully acknowledged.

The computations were partially enabled by resources provided by the Swedish
National Infrastructure for Computing (SNIC) at Chalmers Centre for
Computational Science and Engineering (C3SE) partially funded by the Swedish
Research Council through grant agreement no. 2018-05973.

## References

  * Arazo et al. [2020] E. Arazo, D. Ortego, P. Albert, N. E. O’Connor, and K. McGuinness.  Pseudo-labeling and confirmation bias in deep semi-supervised learning.  In _2020 International Joint Conference on Neural Networks (IJCNN)_ , pages 1–8, 2020.  10.1109/IJCNN48605.2020.9207304.
  * Bachman et al. [2014] P. Bachman, O. Alsharif, and D. Precup.  Learning with pseudo-ensembles.  In Z. Ghahramani, M. Welling, C. Cortes, N. D. Lawrence, and K. Q. Weinberger, editors, _Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada_ , pages 3365–3373, 2014.  URL https://proceedings.neurips.cc/paper/2014/hash/66be31e4c40d676991f2405aaecc6934-Abstract.html.
  * Berthelot et al. [2019] D. Berthelot, N. Carlini, I. J. Goodfellow, N. Papernot, A. Oliver, and C. Raffel.  Mixmatch: A holistic approach to semi-supervised learning.  In H. M. Wallach, H. Larochelle, A. Beygelzimer, F. d’Alché-Buc, E. B. Fox, and R. Garnett, editors, _Advances in Neural Information Processing Systems 32: Annual Conference on Neural Information Processing Systems 2019, NeurIPS 2019, December 8-14, 2019, Vancouver, BC, Canada_ , pages 5050–5060, 2019.  URL https://proceedings.neurips.cc/paper/2019/hash/1cd138d0499a68f4bb72bee04bbec2d7-Abstract.html.
  * Cubuk et al. [2020] E. D. Cubuk, B. Zoph, J. Shlens, and Q. Le.  Randaugment: Practical automated data augmentation with a reduced search space.  In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.  URL https://proceedings.neurips.cc/paper/2020/hash/d85b63ef0ccb114d0a3bb7b7d808028f-Abstract.html.
  * [5] CuriousAI.  Curiousai/mean-teacher: A state-of-the-art semi-supervised method for image recognition.  URL https://github.com/CuriousAI/mean-teacher/blob/51d20acd65/README.md#tips-for-choosing-hyperparameters-and-other-tuning.
  * Grandvalet and Bengio [2004] Y. Grandvalet and Y. Bengio.  Semi-supervised learning by entropy minimization.  In _Advances in Neural Information Processing Systems 17 [Neural Information Processing Systems, NIPS 2004, December 13-18, 2004, Vancouver, British Columbia, Canada]_ , pages 529–536, 2004.  URL https://proceedings.neurips.cc/paper/2004/hash/96f2b50b5d3613adf9c27049b2a888c7-Abstract.html.
  * Krizhevsky [2009] A. Krizhevsky.  Learning multiple layers of features from tiny images.  Technical report, 2009.  URL https://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf.
  * Laine and Aila [2017] S. Laine and T. Aila.  Temporal ensembling for semi-supervised learning.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net, 2017.  URL https://openreview.net/forum?id=BJ6oOfqge.
  * Lee et al. [2013] D.-H. Lee et al.  Pseudo-label: The simple and efficient semi-supervised learning method for deep neural networks.  In _Workshop on challenges in representation learning, ICML_ , volume 3-2, page 896, 2013.  URL https://storage.googleapis.com/kaggle-forum-message-attachments/8924/pseudo_label_final.pdf.
  * Loshchilov and Hutter [2017] I. Loshchilov and F. Hutter.  SGDR: stochastic gradient descent with warm restarts.  In _5th International Conference on Learning Representations, ICLR 2017, Toulon, France, April 24-26, 2017, Conference Track Proceedings_. OpenReview.net, 2017.  URL https://openreview.net/forum?id=Skq89Scxx.
  * Miyato et al. [2019] T. Miyato, S. I. Maeda, M. Koyama, and S. Ishii.  Virtual adversarial training: A regularization method for supervised and semi-supervised learning.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 41:1979–1993, 8 2019.  ISSN 19393539.  10.1109/TPAMI.2018.2858821.
  * Nesterov [1983] Y. E. Nesterov.  A method for solving the convex programming problem with convergence rate $o(1/k^{2})$.  _Dokl. Akad. Nauk SSSR_ , 269:543–547, 1983.
  * Nigam and Ghani [2000] K. Nigam and R. Ghani.  Analyzing the effectiveness and applicability of co-training.  In _Proceedings of the Ninth International Conference on Information and Knowledge Management_ , CIKM ’00, page 86–93, New York, NY, USA, 2000. Association for Computing Machinery.  ISBN 1581133200.  10.1145/354756.354805.
  * Oliver et al. [2018] A. Oliver, A. Odena, C. Raffel, E. D. Cubuk, and I. J. Goodfellow.  Realistic evaluation of deep semi-supervised learning algorithms.  In S. Bengio, H. M. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, _Advances in Neural Information Processing Systems 31: Annual Conference on Neural Information Processing Systems 2018, NeurIPS 2018, December 3-8, 2018, Montréal, Canada_ , pages 3239–3250, 2018.  URL https://proceedings.neurips.cc/paper/2018/hash/c1fea270c48e8079d8ddf7d06d26ab52-Abstract.html.
  * Polyak [1964] B. Polyak.  Some methods of speeding up the convergence of iteration methods.  _USSR Computational Mathematics and Mathematical Physics_ , 4(5):1–17, 1964.  ISSN 0041-5553.  10.1016/0041-5553(64)90137-5.
  * Rasmus et al. [2015] A. Rasmus, M. Berglund, M. Honkala, H. Valpola, and T. Raiko.  Semi-supervised learning with ladder networks.  In C. Cortes, N. D. Lawrence, D. D. Lee, M. Sugiyama, and R. Garnett, editors, _Advances in Neural Information Processing Systems 28: Annual Conference on Neural Information Processing Systems 2015, December 7-12, 2015, Montreal, Quebec, Canada_ , pages 3546–3554, 2015.  URL https://proceedings.neurips.cc/paper/2015/hash/378a063b8fdb1db941e34f4bde584c7d-Abstract.html.
  * Sajjadi et al. [2016] M. Sajjadi, M. Javanmardi, and T. Tasdizen.  Regularization with stochastic transformations and perturbations for deep semi-supervised learning.  In D. D. Lee, M. Sugiyama, U. von Luxburg, I. Guyon, and R. Garnett, editors, _Advances in Neural Information Processing Systems 29: Annual Conference on Neural Information Processing Systems 2016, December 5-10, 2016, Barcelona, Spain_ , pages 1163–1171, 2016.  URL https://proceedings.neurips.cc/paper/2016/hash/30ef30b64204a3088a26bc2e6ecf7602-Abstract.html.
  * Sohn et al. [2020] K. Sohn, D. Berthelot, N. Carlini, Z. Zhang, H. Zhang, C. Raffel, E. D. Cubuk, A. Kurakin, and C. Li.  Fixmatch: Simplifying semi-supervised learning with consistency and confidence.  In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.  URL https://proceedings.neurips.cc/paper/2020/hash/06964dce9addb1c5cb5d6e3d9838f733-Abstract.html.
  * Sutskever et al. [2013] I. Sutskever, J. Martens, G. E. Dahl, and G. E. Hinton.  On the importance of initialization and momentum in deep learning.  In _Proceedings of the 30th International Conference on Machine Learning, ICML 2013, Atlanta, GA, USA, 16-21 June 2013_ , volume 28 of _JMLR Workshop and Conference Proceedings_ , pages 1139–1147. JMLR.org, 2013.  URL http://proceedings.mlr.press/v28/sutskever13.html.
  * Tarvainen and Valpola [2017] A. Tarvainen and H. Valpola.  Mean teachers are better role models: Weight-averaged consistency targets improve semi-supervised deep learning results.  In I. Guyon, U. von Luxburg, S. Bengio, H. M. Wallach, R. Fergus, S. V. N. Vishwanathan, and R. Garnett, editors, _Advances in Neural Information Processing Systems 30: Annual Conference on Neural Information Processing Systems 2017, December 4-9, 2017, Long Beach, CA, USA_ , pages 1195–1204, 2017.  URL https://proceedings.neurips.cc/paper/2017/hash/68053af2923e00204c3ca7c6a3150cf7-Abstract.html.
  * van Engelen and Hoos [2020] J. E. van Engelen and H. H. Hoos.  A survey on semi-supervised learning.  _Machine Learning_ , 109:373–440, 2 2020.  ISSN 15730565.  10.1007/s10994-019-05855-6.
  * Xie et al. [2020a] Q. Xie, Z. Dai, E. H. Hovy, T. Luong, and Q. Le.  Unsupervised data augmentation for consistency training.  In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020a.  URL https://proceedings.neurips.cc/paper/2020/hash/44feb0096faa8326192570788b38c1d1-Abstract.html.
  * Xie et al. [2020b] Q. Xie, M. Luong, E. H. Hovy, and Q. V. Le.  Self-training with noisy student improves imagenet classification.  In _2020 IEEE/CVF Conference on Computer Vision and Pattern Recognition, CVPR 2020, Seattle, WA, USA, June 13-19, 2020_ , pages 10684–10695. Computer Vision Foundation / IEEE, 2020b.  10.1109/CVPR42600.2020.01070.
  * Zagoruyko and Komodakis [2016] S. Zagoruyko and N. Komodakis.  Wide residual networks.  In E. R. H. Richard C. Wilson and W. A. P. Smith, editors, _Proceedings of the British Machine Vision Conference (BMVC)_ , pages 87.1–87.12. BMVA Press, 9 2016.  ISBN 1-901725-59-6.  10.5244/C.30.87.
  * Zoph et al. [2020] B. Zoph, G. Ghiasi, T. Lin, Y. Cui, H. Liu, E. D. Cubuk, and Q. Le.  Rethinking pre-training and self-training.  In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, _Advances in Neural Information Processing Systems 33: Annual Conference on Neural Information Processing Systems 2020, NeurIPS 2020, December 6-12, 2020, virtual_ , 2020.  URL https://proceedings.neurips.cc/paper/2020/hash/27e9661e033a73a6ad8cefcde965c54d-Abstract.html.