# Accelerating Multilevel Markov Chain Monte Carlo Using Machine Learning
Models

Sohail Reddy
Lawrence Livermore National Laboratory
Liveremore, CA 94550
<EMAIL_ADDRESS>
& Hillary Fairbanks
Lawrence Livermore National Laboratory
Liveremore, CA 94550
<EMAIL_ADDRESS>

###### Abstract

This work presents an efficient approach for accelerating multilevel Markov
Chain Monte Carlo (MCMC) sampling for large-scale problems using low-fidelity
machine learning models. While conventional techniques for large-scale
Bayesian inference often substitute computationally expensive high-fidelity
models with machine learning models, thereby introducing approximation errors,
our approach offers a computationally efficient alternative by augmenting
high-fidelity models with low-fidelity ones within a hierarchical framework.
The multilevel approach utilizes the low-fidelity machine learning model (MLM)
for inexpensive evaluation of proposed samples thereby improving the
acceptance of samples by the high-fidelity model. The hierarchy in our
multilevel algorithm is derived from geometric multigrid hierarchy. We utilize
an MLM to acclerate the coarse level sampling. Training machine learning model
for the coarsest level significantly reduces the computational cost associated
with generating training data and training the model. We present an MCMC
algorithm to accelerate the coarsest level sampling using MLM and account for
the approximation error introduced. We provide theoretical proofs of detailed
balance and demonstrate that our multilevel approach constitutes a consistent
MCMC algorithm. Additionally, we derive conditions on the accuracy of the
machine learning model to facilitate more efficient hierarchical sampling. Our
technique is demonstrated on a standard benchmark inference problem in
groundwater flow, where we estimate the probability density of a quantity of
interest using a four-level MCMC algorithm. Our proposed algorithm accelerates
multilevel sampling by a factor of two while achieving similar accuracy
compared to sampling using the standard multilevel algorithm.

_Keywords_ Bayesian Inference, Machine Learning, Multilevel Markov Chain Monte
Carlo, Stochastic PDE

##  1 Introduction

In the realm of statistical inference, Bayesian methods offer a powerful
framework for decision making under uncertainty by integrating prior knowledge
with observational data. However, the computational complexity of exact
Bayesian inference often becomes prohibitive, particularly for models with
high-dimensional parameter spaces or complex likelihood functions. To address
this challenge, Markov Chain Monte Carlo (MCMC) methods have emerged as
indispensable tools [1]. MCMC algorithms provide a means to approximate
posterior distributions by generating samples from them, circumventing the
need for analytical solutions which may be intractable or computationally
expensive. This has led to widespread use across many disciplines. Although
MCMC sampling has alleviated the challenges of analytically computing the
posterior distribution, the method’s accuracy and sampling efficiency is still
directly governed by the likelihood and the forward model. Employing numerical
models for likelihood computations requires high-fidelity solutions, often on
a finely resolved space-time discretization, which greatly increases the
computational cost of the MCMC sampling.

To enable MCMC sampling with computationally expensive forward maps, reduced
order or surrogate models have been extensively employed [2]. These surrogate
models, in the form of Gaussian processes [3, 4], deep neural networks [5, 6]
and radial basis functions [7], have been employed to reduce the computational
cost of the forward model at the expense of accuracy. Physics constrained
reduced order models, such as those based on proper orthogonal decomposition
[8, 9], have been investigated to yield more robust predictions at a fraction
of the cost. While surrogate and reduced model greatly reduce the cost of the
forward map, the approximation error they introduce can lead to inefficient
sampling and errors in the statistics. To address this, delayed acceptance
MCMC , as developed in [10, 11, 12, 13], uses a combination of low and high
fidelity models whereby the low fidelity model is used to inexpensively
evaluate candidate proposals before evaluating with a high fidelity model.
Although, this two level approach improves sampling efficiency and accuracy by
reducing the impact of approximation error, it does not treat the inefficiency
of sampling in high dimensional space. To address this issue, multilevel
delayed acceptance (MLDA) sampling approaches based on hierarchical
decomposition of the sampling space have become an area of active research
[14, 15, 16]. This approach uses sampling in a low dimensional space to inform
high dimensional sampling and has been demonstrated on problems in subsurface
flows. The extensive work done on surrogate model accelerated MCMC and
multilevel MCMC has naturally led to a coupled approach of machine learning
model augmented multilevel MCMC, where surrogate models are used to further
accelerate levels in the hierarchy. Approximation errors introduced by
inaccurate machine learning models can lead to a divergence of the chain to
regions of low posterior densities, thereby ill-informing sampling on
subsequent finer levels which leads to significant increase in computational
cost and decrease in overall efficiency. Lykkegaard et al.[17] used a deep
neural network (DNN) to accelerate the lowest-fidelity (coarsest) of a two
level hierarchy. They reported difficulty in training the DNN models to
predict the forward map with sufficient accuracy, as the coefficient (input
vector) dimension was smaller than that of the observational data (output
vector) on which they trained their model. They remedied this issue by
decomposing the observational domain into subdomains and training a separate
DNN model to approximate the forward map on each subdomain.

In this work, we consider a geometric multigrid hierarchy of varying level of
refinement and present a multilevel MCMC algorithm that increases the
efficiency of the coarsest level by augmenting the hierarchy with a DNN model.
Our approach utilizes a two stage approach on the coarsest level where the
first stage performs an inexpensive Metropolis-Hastings step using the DNN
model while the second state performs a filtering step using the finite
element forward model. This two stage algorithm prevents divergence of the
coarse level chain due to the approximation error of the DNN model and is
embedded within a multilevel delayed acceptance algorithm to extend the
sampling over all the levels of the multigrid hierarchy. The technique is
demonstrated on a problem in subsurface flow, namely Darcy flow with uncertain
permeability coefficient, whereby both the coefficient and Darcy pressure
solution are solved within the mixed finite element framework. Dissimilar to
the work of [17], we use a PDE-based approach to generate coefficient data
used in the MLM training, as opposed to the commonly used Karhunen-Loève (KL)
expansion [18]. Consequently, our coefficient (input vector) dimension is
large enough to mitigate the need for multiple DNN models to predict the
observational output.

The remainder of the paper is organized as follows. The Bayesian inverse
problem is introduced in Section 2. This includes a description hierarchical
Gaussian random field (GRF) sampling, and our filtered, multilevel MCMC
algorithm. We prove that our hierarchical approach preserves detailed balance
and is a consistent MCMC algorithm. In Section 3 we present the problem
formulation for Bayesian inference in subsurface flows including the structure
of our machine learning model. Finally, in Section 4 we present the results of
our machine learning accelerated MLMCMC algorithm and compare the filtered
algorithm against the unfiltered algorithms. Our results demonstrate a factor
of two speed up over the standard PDE-only hierarchy with improved statistics
while maintaining high accuracy.

##  2 Bayesian Inverse Problem: Markov Chain Monte Carlo

Bayesian inverse problems aim to infer the unknown parameters of a system or
model by incorporating prior knowledge and observed data in a probabilistic
framework. By leveraging Bayes’ theorem, which updates prior information
encoded in the prior distribution with observational data,
$\mathrm{\boldsymbol{y}}_{obs}$, to obtain posterior distributions. Bayesian
inverse problems offer a principled approach to estimating parameters and
assessing their uncertainty, thereby facilitating informed decision-making and
improving the reliability of predictions. While there exist several
approaches, such as maximum-a-posteriori (MAP) estimate [19] and variational
inference [20], for approximating the posterior distribution and quantifying
uncertainty in model parameters, we employ MCMC sampling.

This work considers multilevel MCMC sampling from a posterior distribution. To
that end, we define a hierarchy of maps
$\left\\{\mathcal{F}_{\ell}\right\\}_{\ell=0}^{L}$ with increasing fidelity.
Each
$\mathcal{F}_{\ell}:=\mathcal{D}_{\ell}\circ\mathcal{S}_{\ell}(\zeta_{\ell})$
maps the parameters being inferred $\zeta$ to model output
$\mathrm{\boldsymbol{y}}$ and is a composition of
$\mathcal{S}_{\ell}:\zeta_{\ell}\mapsto\theta_{\ell}$ corresponding to the
sampler for generating GRFs and
$\mathcal{D}_{\ell}:\theta_{\ell}\mapsto\mathrm{\boldsymbol{y}}_{\ell}$
corresponding to forward model (e.g. Darcy solver). Although a hierarchy of
models can be constructed by considering models with increasing complex
physics, we consider it a hierarchy of increasingly refined discretizations of
a single PDE, that naturally stem from multigrid methods. Hence, $\ell$
denotes the order of refinement. For the remainder of the manuscript, we
denote variables pertaining to the surrogate model by $\widehat{\cdot}$.

###  2.1 Hierarchical Sampling of Gaussian Random Fields

Gaussian random fields have found application in several fields including
statistical physics, biomedical imaging, machine learning, and parameterizing
material distribution [21]. Consider a domain
$\Omega\subset\mathbb{R}^{d}~{}(d=2,3)$ with a boundary
$\Gamma:=\partial\Omega$, and a probability space
$(\Omega_{\mathcal{S}},\mathcal{F},\mathbb{P})$ where $\Omega_{\mathcal{S}}$
is the sample space (the set of all possible events), $\mathcal{F}$ is the
$\sigma$-algebra of events, and $\mathbb{P}:\mathcal{F}\rightarrow[0,1]$ is a
probability measure. We then seek functions
$\left\\{\theta(\mathrm{\boldsymbol{x}},\omega)\in
L^{2}(\Omega)|\mathrm{\boldsymbol{x}}\in\Omega,\omega\in\Omega_{\mathcal{S}}\right\\}$
that represent realizations of Gaussian random fields that follow
$\theta\sim\mathcal{N}(0,\mathcal{C})$ with covariance $\mathcal{C}$. A
popular approach for sampling realizations of GRF involves a KL expansion of
the covariance kernel [18]. Although applicable for lower dimensional problem
with a spectrally-fast-decaying kernel, the poor scaling of KL decomposition
limits its use for large scale applications. Furthermore, KL expansion is
often truncated to a finite number of modes, hence only span a subspace.
Efforts to include higher frequency modes in the expansion often leads to
numerical artifacts and pollution of the GRFs leading to level-dependent
truncated expansion [22]. To address these drawbacks, we sample random fields
by solving a stochastic partial differential equation (SPDE) using finite
element discretization [23], and scalable multigrid solvers [24, 25, 26]. The
SPDE approach for sampling GRFs considers the following fractional PDE

$\left(\kappa^{2}-\Delta\right)^{\alpha/2}\theta(\mathrm{\boldsymbol{x}},\omega)=g\mathcal{W}(\mathrm{\boldsymbol{x}},\omega),\quad~{}\mathrm{\boldsymbol{x}}\in\mathbb{R}^{d},~{}\alpha=\nu+\dfrac{d}{2},~{}\kappa>0,~{}\nu>0$
(1)

where the scaling factor

$g=(4\pi)^{d/4}\kappa^{\nu}\sqrt{\dfrac{\Gamma(\nu+d/2)}{\Gamma(\nu)}}$

imposes unit marginal variance of $\theta(\mathrm{\boldsymbol{x}},\omega)$ and
the white noise $\mathcal{W}(\mathrm{\boldsymbol{x}},\omega)$ is an
$L^{2}(\Omega)$ generalized function satisfying

$(\mathcal{W}(\mathrm{\boldsymbol{x}},\omega),\phi)\sim\mathcal{N}(0,||\phi||^{2}_{L^{2}(\Omega)}),\quad\forall\phi\in
L^{2}(\Omega)$

Setting $\nu=\frac{1}{2}$ in three dimension, we obtain the standard, integer-
order, reaction-diffusion equation

$\nabla\cdot{\mathrm{\boldsymbol{\rho}}}-\kappa^{2}\theta=-g\mathcal{W}(\mathrm{\boldsymbol{x}},\omega)$
(2a) $\mathrm{\boldsymbol{\rho}}=\nabla{\theta}$ (2b)

the solution to which yields realizations of GRFs with exponential covariance
function

$\mathrm{cov}(\mathrm{\boldsymbol{x}},\mathrm{\boldsymbol{y}})=\sigma^{2}e^{-\kappa\left\|\mathrm{\boldsymbol{x}}-\mathrm{\boldsymbol{y}}\right\|}$
(3)

We solve the integer-order PDE using the mixed finite element method and
consider the function spaces

$\mathrm{\boldsymbol{R}}=H(\mathrm{div},\Omega):=\left\\{\mathrm{\boldsymbol{u}}\in\mathrm{\boldsymbol{L}}^{2}(\Omega)~{}|~{}\mathrm{div}~{}\mathrm{\boldsymbol{u}}\in
L^{2}(\Omega),\mathrm{\boldsymbol{u}}\cdot\mathrm{\boldsymbol{n}}=0~{}\mathrm{on}~{}\Gamma\right\\}$
$\Theta=L^{2}(\Omega)$

where $\mathrm{\boldsymbol{L}}^{2}(\Omega)=[L^{2}(\Omega)]^{d}$ is the
$d$-dimensional vector function space.

Denote by $\Omega_{h}:=\bigcup_{i=1}^{N}K_{i}$ a partition of $\Omega$ into a
finite collection of non-overlapping elements $K_{i}$ with
$h:=\underset{j=1,\ldots,N}{\max}diam(K_{j})$. We take
$\mathrm{\boldsymbol{R}}_{h}\subset
H(\mathrm{div},\Omega_{h})\subset\mathrm{\boldsymbol{R}}$ and
$\Theta_{h}\subset L^{2}(\Omega_{h})\subset\Theta$ to be the lowest order
Raviart-Thomas and piecewise constant basis functions, respectively.
Furthermore, we consider
$\mathrm{\boldsymbol{\rho}}\in\mathrm{\boldsymbol{R}}$ and $\theta\in\Theta$
and denote their finite element representation by
$\mathrm{\boldsymbol{\rho}}_{h}\in\mathrm{\boldsymbol{R}}_{h}$ and
$\theta_{h}\in\Theta_{h}$. Introducing test functions
$\mathrm{\boldsymbol{v}}_{h}\in\mathrm{\boldsymbol{R}}_{h}$ and
$q_{h}\in\Theta_{h}$, the weak, mixed form of Eq. (2) reads

######  Problem 1.

Find
$\left(\mathrm{\boldsymbol{\rho}}_{h},\theta_{h}\right)\in\mathrm{\boldsymbol{R}}_{h}\times\Theta_{h}$
such that

$\displaystyle\left(\nabla\cdot{\mathrm{\boldsymbol{\rho}}_{h}},q_{h}\right)-\kappa^{2}\left(\theta_{h},q_{h}\right)$
$\displaystyle=-g\left(\mathcal{W}(\omega),q_{h}\right),\quad$
$\displaystyle\forall\mathrm{\boldsymbol{v}}_{h}\in\mathrm{\boldsymbol{R}}_{h}$
(4)
$\displaystyle\left(\mathrm{\boldsymbol{\rho}}_{h},\mathrm{\boldsymbol{v}}_{h}\right)+\left(\nabla\cdot{\mathrm{\boldsymbol{v}}_{h}},\theta_{h}\right)$
$\displaystyle=0,\quad$ $\displaystyle\forall q_{h}\in\Theta_{h}$

with boundary conditions
$\mathrm{\boldsymbol{\rho}}_{h}\cdot\mathrm{\boldsymbol{n}}=0$.

with inner products defined as

$\left(\theta,q\right)=\int_{\Omega}{\theta~{}q~{}d\Omega},~{}\mathrm{and}~{}\left(\mathrm{\boldsymbol{\rho}},\mathrm{\boldsymbol{v}}\right)=\int_{\Omega}{\mathrm{\boldsymbol{\rho}}\cdot\mathrm{\boldsymbol{v}}~{}d\Omega}$

For ease of notation, letting
$\zeta_{h}:=\left(\mathcal{W}_{h}(\omega),q_{h}\right)=W_{h}\mathcal{W}_{h}=W_{h}^{1/2}\xi_{h}$,
the resulting system is then given by

$\begin{bmatrix}M_{h}&B_{h}^{T}\\\
B_{h}&-\kappa^{2}W_{h}\end{bmatrix}\begin{bmatrix}\mathrm{\boldsymbol{\rho}}_{h}\\\
\theta_{h}\end{bmatrix}=\begin{bmatrix}\mathrm{\boldsymbol{0}}\\\
-g\zeta_{h}\end{bmatrix}$ (5)

where $\Theta_{h}\ni\xi_{h}\sim\mathcal{N}(0,I)$, $M_{h}$ is the mass matrix
corresponding to Raviart-Thomas space $\mathrm{\boldsymbol{R}}_{h}$, $W_{h}$
is the mass matrix corresponding to the space $\Theta_{h}$, and $B_{h}$ is the
discrete divergence operator. Solution to Eq. (5) with independently sampled
$\mathcal{W}_{h}$ yields a realization of a GRF independent of the other
levels in the hierarchy. For multilevel sampling, we require the fine-level
white noise $\mathcal{W}_{h}$ to be conditioned on the coarse-level white
noise $\mathcal{W}_{H}$ (i.e. $\mathcal{W}_{h}|\mathcal{W}_{h<H}$).

Consider two levels $\ell$ and $L$ in a hierarchy and their corresponding
independent white noise $\mathcal{W}_{\ell}\sim\mathcal{N}(0,W^{-1}_{\ell})$
and $\mathcal{W}_{L}\sim\mathcal{N}(0,W^{-1}_{L})$, respectively, with
$L>\ell$ implying $h<H$ with nested spaces
$\mathrm{\boldsymbol{R}}_{H}\subset\mathrm{\boldsymbol{R}}_{h}$ and
$\Theta_{H}\subset\Theta_{h}$. Here onwards, we replace the subscripts $h$ and
$H$ with $L$ and $\ell$, respectively, with
$\mathcal{G}_{H}\subset\mathcal{G}_{h}\rightarrow\mathcal{G}_{\ell}\subset\mathcal{G}_{L},~{}\mathcal{G}\in\left\\{\mathrm{\boldsymbol{R}},\Theta\right\\}$.
Furthermore, denote by
$\Theta_{L}\ni\mathcal{W}_{\ell}^{L}\sim\mathcal{N}(0,W^{-1}_{\ell})$, the
realization of white noise $\mathcal{W}_{\ell}$ on level $L$. Using standard
multigrid prolongation operator $P$ and restriction operator $\Pi$, we can
write $P\Pi\mathcal{W}_{\ell}^{L}=P\mathcal{W}_{\ell}$. Following the work of
[16], we then construct a hierarchical decomposition of the white noise on
level $L$, conditioned on $\mathcal{W}_{\ell}$, as

$\widetilde{\mathcal{W}}_{L}=\underbrace{P\mathcal{W}_{\ell}}_{\mathcal{Q}_{L}\mathcal{W}_{\ell}^{L}}+(I-\mathcal{Q}_{L})\mathcal{W}_{L}$
(6)

where $\mathcal{Q}_{L}:\Theta_{\ell}\mapsto\Theta_{L}$ is an
$L^{2}$-projection with $\mathcal{Q}_{L}=P\Pi$, and
$\Pi=W_{\ell}^{-1}P^{T}W_{L}$ with $\Pi P=I$. The first term on the right
($P\mathcal{W}_{\ell}\sim\mathcal{N}(0,W_{\ell})$) represents the prolongation
of the coarse level white noise to the fine level, while the second term
represents the white noise complement in $\Theta_{L}\setminus\Theta_{\ell}$.
For ease of implementation, we employ the following form of the decomposition
of the source term

$\displaystyle\widetilde{\zeta}_{L}:=W_{L}\widetilde{\mathcal{W}}_{L}$
$\displaystyle=W_{L}P\mathcal{W}_{\ell}+W_{L}(I-\mathcal{Q}_{L})\mathcal{W}_{L}$
(7)
$\displaystyle=W_{L}PW^{-1}_{\ell}\zeta_{\ell}+W_{L}(I-PW_{\ell}^{-1}P^{T}W_{L})W^{-1}_{L}\zeta_{L}$
$\displaystyle=\Pi^{T}\zeta_{\ell}+(I-\Pi^{T}P^{T})\zeta_{L}$
$\displaystyle=\Pi^{T}W^{1/2}_{\ell}\xi_{\ell}+(I-\Pi^{T}P^{T})W^{1/2}_{L}\xi_{L}\sim\mathcal{N}(0,W_{L})$

The decomposition Eq. (7) can be applied recursively to finer levels to obtain
conditioned white noise, and therefore, realizations of conditioned GRFs on
finer levels. Furthermore, note that in Eq. (7) the matrix square-root is
needed. Since, $\Theta_{\ell}$ is taken to be piecewise constant, $W_{\ell}$
is diagonal. We employ mesh-embedding [25] to avoid boundary artifacts arising
from no-flux boundary conditions on finite domains.

Following this approach, we are now able to sample fine and coarse
realizations of white noise from a joint distribution, which is an essential
component of multilevel sampling.

###  2.2 Multi-Level Markov Chain Monte Carlo

Our goal is to estimate the posterior mean of a scalar quantity of interest
$Q:=Q(\zeta)$

$\mathbb{E}_{\pi}[Q]=\int_{\Omega}{Q(\zeta)\pi(\zeta)~{}d\zeta}$

where the posterior distribution is
$\pi(\zeta)\equiv\pi(\zeta|\mathrm{\boldsymbol{y}}_{obs})\propto\mathcal{L}(\mathrm{\boldsymbol{y}}_{obs}|\zeta)\nu(\zeta)$,
with prior $\nu(\zeta)$ and likelihood
$\mathcal{L}(\mathrm{\boldsymbol{y}}_{obs}|\zeta)$. We estimate the integral
by means of Monte Carlo (MC) sampling by drawing samples from the posterior
using MCMC

$\mathbb{E}_{\pi}[Q]\approx\overline{Q}=\dfrac{1}{N}\sum_{i=1}^{N}Q(\zeta^{(i)})$

where $\zeta^{(i)}\sim\pi(\zeta|\mathrm{\boldsymbol{y}}_{obs})$. Using
$\zeta_{\ell}$ as the discrete approximation of white noise on level $\ell$,
sampled according the posterior distribution on level $\ell$,
$\pi_{\ell}(\zeta)\equiv\pi(\zeta_{\ell}|\mathrm{\boldsymbol{y}}_{obs})$, we
define the discrete QoI approximation on level $\ell$ as
$Q_{\ell}:=Q_{\ell}(\zeta_{\ell})$. Then the multilevel decomposition of the
expectation is given as

$\mathbb{E}_{\pi_{L}}[Q_{L}]=\mathbb{E}_{\pi_{0}}[Q_{0}]+\sum_{\ell=1}^{L}\left(\mathbb{E}_{\pi_{\ell}}[Q_{\ell}]-\mathbb{E}_{\pi_{\ell-1}}[Q_{\ell-1}]\right).$
(8)

We use the estimator $\widehat{Y}_{\ell}^{N_{\ell}}$ for the difference
$\mathbb{E}_{\pi_{\ell}}[Q_{\ell}]-\mathbb{E}_{\pi_{\ell-1}}[Q_{\ell-1}]$,
defined as

$\widehat{Y}_{\ell}^{N_{\ell}}=\dfrac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}Y^{(i)}_{\ell}=\dfrac{1}{N_{\ell}}\sum_{i=1}^{N_{\ell}}\left(Q^{(i)}_{\ell}-Q^{(i)}_{\ell-1}\right),$
(9)

where $Q^{(i)}_{\ell}:=Q_{\ell}(\zeta_{\ell}^{(i)})$ is the $i^{th}$
independent sample on level $\ell$, and $N_{\ell}$ is the number of samples
performed on level $\ell$.

Multilevel Monte Carlo (MLMCM) and Multilevel MCMC (MLMCMC) obtain variance
reduction by sampling each $Y_{\ell}^{(i)}=Q^{(i)}_{\ell}-Q^{(i)}_{\ell-1}$
from a joint distribution, with expectation and variance denoted as
$\mathbb{E}_{\pi_{\ell},\pi_{\ell-1}}[Y_{\ell}]$ and
$\mathbb{V}_{\pi_{\ell},\pi_{\ell-1}}\left[Y_{\ell}\right]$, respectively. As
described in Section 2.1, we obtain these samples via the hierarchical PDE
sampler approach. Let the computational cost of generating an independent
sample of $Y_{\ell}$ be denoted as $\overline{C}_{\ell}$ for $\ell=0,1,\ldots
L$. For a given MSE tolerance $\epsilon^{2}$, the optimal number of samples on
each level, $N_{\ell}$, is calculated as

$N_{\ell}=\dfrac{2}{\epsilon^{2}}\left(\sum_{k=0}^{L}\sqrt{\mathbb{V}_{\pi_{k},\pi_{k-1}}\left[Y_{k}\right]\overline{C}_{k}}\right)\sqrt{\dfrac{\mathbb{V}_{\pi_{\ell},\pi_{\ell-1}}\left[Y_{\ell}\right]}{\overline{C_{\ell}}}}.$
(10)

If the estimate of $\mathbb{V}_{\pi_{\ell},\pi_{\ell-1}}\left[Y_{\ell}\right]$
decays with level refinement, fewer simulations will be needed on the finer
levels. We refer the reader to Appendix A.2 for details on this derivation.

We sample from the posterior using a Multi-Level Delayed Acceptance (MLDA)
algorithm and generate proposals $\mathrm{\boldsymbol{\zeta}}_{\ell}^{P}$ on
level $\ell$, conditioned on the current state
$\mathrm{\boldsymbol{\zeta}}_{\ell}^{C}$ and the coarse levels using
preconditioned Crank-Nicolson

$\mathrm{\boldsymbol{\zeta}}_{\ell}^{P}:=\sqrt{1-\beta^{2}}\mathrm{\boldsymbol{\zeta}}_{\ell}^{C}+\beta\widetilde{\zeta}_{\ell}$

where $\widetilde{\zeta}_{\ell}$ is the white noise on level $\ell$
conditioned on the coarse level white noise (Eq. (7)), and
$\mathrm{\boldsymbol{\zeta}}_{\ell}^{P}\sim
q(\mathrm{\boldsymbol{\zeta}}|\mathrm{\boldsymbol{\zeta}}_{\ell}^{C},\widetilde{\zeta}_{\ell})=\mathcal{N}(\sqrt{1-\beta^{2}}\mathrm{\boldsymbol{\zeta}}_{\ell}^{C},\beta^{2}W_{\ell})$.

######  Remark 1.

The preconditioned Crank-Nicolson proposal $q(\zeta^{P}|\zeta^{C})$ with a
prior $\nu\sim\mathcal{N}(0,W)$ satisfies
$\nu(\zeta^{C})q(\zeta^{P}|\zeta^{C})=\nu(\zeta^{P})q(\zeta^{C}|\zeta^{P})$.

Since the multilevel sampling redistributes the majority of the computational
costs over the coarser levels, we further accelerate the coarsest level using
a machine learning model. This is acheived by first using the machine learning
model to evaluate the generated proposals followed by a filtering stage using
the PDE model on the same level to filter out poor proposals before proposing
them to the subsequent finer levels. This filtering approach, often referred
to as _Surrogate-Transistion_ [27], is presented in Alg 1, where the first
stage is the standard Metropolis-Hastings algorithm, hence, is ergodic and
satisfies detailed balance. Here, we prove the complete, two-stage algorithm
also satisfies detailed balance.

Algorithm 1 Generate a single sample on a level using the filtered Metropolis-
Hastings algorithm.

Input: Current state: $\mathrm{\boldsymbol{\zeta}}^{C}$
Output: Next state: $\mathrm{\boldsymbol{\zeta}}^{\star}$

$\mathrm{\boldsymbol{\zeta}}^{P}\sim
q(\mathrm{\boldsymbol{\zeta}}|\mathrm{\boldsymbol{\zeta}}^{C})$
$\triangleright$ Generate proposal

$\widehat{\alpha}=\min\left\\{1,\dfrac{\widehat{\mathcal{L}}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}^{P})}{\widehat{\mathcal{L}}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}^{C})}\right\\}$
$\triangleright$ Compute acceptance ratio

if $\widehat{\alpha}\geq u\sim\mathcal{U}(0,1)$ then $\triangleright$ Filter
proposed sample if accepted

$\mathrm{\boldsymbol{\zeta}}^{\star}\leftarrow$
Filter($\mathrm{\boldsymbol{\zeta}}^{C}$,$\mathrm{\boldsymbol{\zeta}}^{P}$)

else$\triangleright$ Reject the proposed sample

$\mathrm{\boldsymbol{\zeta}}^{\star}=\mathrm{\boldsymbol{\zeta}}^{C}$

end if

procedure
Filter($\mathrm{\boldsymbol{\zeta}}^{C}$,$\mathrm{\boldsymbol{\zeta}}^{P}$)

$\alpha=\min\left\\{1,\dfrac{\mathcal{L}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}^{P})}{\mathcal{L}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}^{C})}\dfrac{\widehat{\mathcal{L}}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}^{C})}{\widehat{\mathcal{L}}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}^{P})}\right\\}$
$\triangleright$ Compute multi-level acceptance ratio

if $\alpha\geq u\sim\mathcal{U}(0,1)$ then $\triangleright$ Accept the
proposed sample

return $\mathrm{\boldsymbol{\zeta}}^{P}$

else$\triangleright$ Reject the proposed sample

return $\mathrm{\boldsymbol{\zeta}}^{C}$

end if

end procedure

######  Proposition 1.

Algorithm 1 simulates a Markov chain that is in detailed balance with
$\pi(\cdot)$

###### Proof.

Let
$g(\zeta^{P}|\zeta^{C})=q(\zeta^{P}|\zeta^{C})\widehat{\alpha}(\zeta^{P}|\zeta^{C})$
be the transistion kernel of the first stage that is in detailed balance with
$\widehat{\pi}(\cdot)$

$~{}\widehat{\pi}(\zeta^{C})g(\zeta^{P}|\zeta^{C})=~{}\widehat{\pi}(\zeta^{P})g(\zeta^{C}|\zeta^{P})$

Then the second stage has the transistion kernel

$G(\zeta^{P}|\zeta^{C})=g(\zeta^{P}|\zeta^{C})~{}\min\left\\{1,\dfrac{\pi(\zeta^{P})}{\pi(\zeta^{C})}\dfrac{\widehat{\pi}(\zeta^{C})}{\widehat{\pi}(\zeta^{P})}\right\\}=g(\zeta^{P}|\zeta^{C})~{}\min\left\\{1,\dfrac{\mathcal{L}(\mathrm{\boldsymbol{y}}_{obs}|\zeta^{P})}{\mathcal{L}(\mathrm{\boldsymbol{y}}_{obs}|\zeta^{C})}\dfrac{\widehat{\mathcal{L}}(\mathrm{\boldsymbol{y}}_{obs}|\zeta^{C})}{\widehat{\mathcal{L}}(\mathrm{\boldsymbol{y}}_{obs}|\zeta^{P})}\right\\}$

and satisfies the detailed balance condition

$\displaystyle\pi(\zeta^{C})G(\zeta^{P}|\zeta^{C})$
$\displaystyle=\pi(\zeta^{C})g(\zeta^{P}|\zeta^{C})~{}\min\left\\{1,\dfrac{\pi(\zeta^{P})}{\pi(\zeta^{C})}\dfrac{\widehat{\pi}(\zeta^{C})}{\widehat{\pi}(\zeta^{P})}\right\\}$
(11)
$\displaystyle=\dfrac{\widehat{\pi}(\zeta^{P})g(\zeta^{C}|\zeta^{P})\pi(\zeta^{C})}{\widehat{\pi}(\zeta^{C})}~{}\min\left\\{1,\dfrac{\pi(\zeta^{P})}{\pi(\zeta^{C})}\dfrac{\widehat{\pi}(\zeta^{C})}{\widehat{\pi}(\zeta^{P})}\right\\}$
$\displaystyle=\widehat{\pi}(\zeta^{P})g(\zeta^{C}|\zeta^{P})~{}\min\left\\{\dfrac{\pi(\zeta^{C})}{\widehat{\pi}(\zeta^{C})},\dfrac{\pi(\zeta^{P})}{\widehat{\pi}(\zeta^{P})}\right\\}$
$\displaystyle=\pi(\zeta^{P})G(\zeta^{C}|\zeta^{P})$

∎

######  Remark 2.

As the accuracy of the surrogate model tends to that of the true model (i.e.
$\widehat{\mathcal{F}}(\mathrm{\boldsymbol{\zeta}})\rightarrow\mathcal{F}(\mathrm{\boldsymbol{\zeta}})$)
then
$\widehat{\mathcal{L}}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}})\rightarrow\mathcal{L}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}})$,
and the acceptance ratio of the filter $\alpha\rightarrow 1$. Hence, for
$\widehat{\mathcal{F}}(\mathrm{\boldsymbol{\zeta}})=\mathcal{F}(\mathrm{\boldsymbol{\zeta}})$,
the filter is idempotent and the two-stage algorithm is equivalent to a
single-stage Metropolis-Hastings algorithm.

Since the first stage performs an inexpensive evalution of the proposed
sampled with $\widehat{\mathcal{F}}$ followed by a filtering stage with
$\mathcal{F}$, it is desirable for the acceptance probability of the second
stage to be greater than (or equal to) that of the first stage. The filtered-
MCMC sampling of the coarsest level is then coupled to the remaining, finer
levels of the hierarchy using the standard multi-level delayed acceptance
(MLDA) algorithm. This coupled algorithm, given in Alg. 2, further accelerates
the complete MLDA hierarchical sampling. Furthermore, this couple sampling
satisfies the properties of the MCMC samplers, as demonstrated by the
following lemma.

Algorithm 2 Generate a single sample on the fine level $L$ conditioned on the
coarse level $\ell$.

Input: Current state: $\mathrm{\boldsymbol{\zeta}}_{L}^{C}$,
$\mathrm{\boldsymbol{\zeta}}_{\ell}^{C}$
Output: Next state: $\mathrm{\boldsymbol{\zeta}}_{L}^{\star}$

Generate $\mathrm{\boldsymbol{\zeta}}_{\ell}^{\star}$ using (filtered)
Metropolis-Hasting using $\mathrm{\boldsymbol{\zeta}}_{\ell}^{C}$
$\triangleright$ Generate sample on the coarse level

Compute $\widetilde{\zeta}_{L}$ using Eq. (7) and
$\mathrm{\boldsymbol{\zeta}}_{\ell}^{\star}$ $\triangleright$ Generate coarse-
level conditioned white noise

$\mathrm{\boldsymbol{\zeta}}_{L}^{P}\sim
q(\mathrm{\boldsymbol{\zeta}}_{L}|\mathrm{\boldsymbol{\zeta}}_{L}^{C},\widetilde{\zeta}_{L})$
$\triangleright$ Generate fine-level proposal

$\alpha_{ML}=\min\left\\{1,\dfrac{\mathcal{L}_{L}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}_{L}^{P})\mathcal{L}_{\ell}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}_{\ell}^{C})}{\mathcal{L}_{L}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}_{L}^{C})\mathcal{L}_{\ell}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}}_{\ell}^{\star})}\right\\}$
$\triangleright$ Compute multilevel acceptance ratio

if $\alpha_{ML}\geq u\sim\mathcal{U}(0,1)$ then $\triangleright$ Accept
proposed sample

$\mathrm{\boldsymbol{\zeta}}_{L}^{\star}=\mathrm{\boldsymbol{\zeta}}_{L}^{P}$

else$\triangleright$ Reject the proposed sample

$\mathrm{\boldsymbol{\zeta}}_{L}^{\star}=\mathrm{\boldsymbol{\zeta}}_{L}^{C}$

end if

######  Lemma 1.

Algorithm 2 with coarse sampling of $\pi_{\ell}(\cdot)$ by Alg. 1 simulates a
Markov chain that is in detailed balance with $\pi_{L}(\cdot)$

###### Proof.

Since the coarse level sampling is in detailed balance with
$\pi_{\ell}(\cdot)$ by Proposition 1, the remainder of the proof follows from
[15, Lemma 2.5]. ∎

######  Remark 3.

Since the Metropolis-Hasting (first stage of Alg. 1), the filtered Metropolis-
Hasting (Alg. 1) and the MLDA (Alg. 2) algorithms all satisfy detailed balance
and, therefore, detailed balance on each level, and are ergodic on each level,
the posterior distributions on levels will converge asymptotically to the
stationary distribution. Furthermore, all levels sampled using filtered
Metropolis-Hastings and PDE based hierarchy will converge to the same
stationary distribution; this is not guarenteed for unfiltered sampling
approach.

##  3 Subsurface Flow: Darcy’s Equations

We demonstrate our surrogate-augmented MLMCMC approach on a typical problem in
groundwater flow described by Darcy’s Law, which has been employed extensively
in verifying the efficiency of delayed acceptance MCMC algorithms. The
governing equations with the appropriate boundary conditions are written as

$\displaystyle\begin{split}\mathrm{\boldsymbol{u}}(\mathrm{\boldsymbol{x}})+k(\mathrm{\boldsymbol{x}})\nabla{p(\mathrm{\boldsymbol{x}})}=\mathrm{\boldsymbol{f}}&\quad\mathrm{in~{}}\Omega\end{split}$
(12a)
$\displaystyle\begin{split}\nabla\cdot{\mathrm{\boldsymbol{u}}(\mathrm{\boldsymbol{x}})}=0&\quad\mathrm{in~{}}\Omega\end{split}$
(12b)
$\displaystyle\begin{split}p=p_{s}&\quad\mathrm{on~{}}\Gamma_{D}\end{split}$
(12c)
$\displaystyle\begin{split}\mathrm{\boldsymbol{u}}\cdot\mathrm{\boldsymbol{n}}=0&\quad\mathrm{on~{}}\Gamma_{N}\end{split}$
(12d)

with $\Gamma:=\Gamma_{D}\bigcup\Gamma_{N}=\partial\Omega$ where Eq. (12b) is
the incompressibility condition. The permeability coefficient is assumed to
follow log-normal distribution
$k(\mathrm{\boldsymbol{x}})=\exp(\theta(\mathrm{\boldsymbol{x}},\omega))$
where $\theta(\mathrm{\boldsymbol{x}},\omega)$ is a realization of a Gaussian
random field with covariance kernel defined by Eq. (3). We consider
$\mathrm{\boldsymbol{u}}\in\mathrm{\boldsymbol{R}}$ and $p\in\Theta$ and
denote their finite element representation by
$\mathrm{\boldsymbol{u}}_{h}\in\mathrm{\boldsymbol{R}}_{h}$ and
$p_{h}\in\Theta_{h}$. Then the weak, mixed form of Eq. (12) reads

######  Problem 2.

Find
$\left(\mathrm{\boldsymbol{u}}_{h},p_{h}\right)\in\mathrm{\boldsymbol{R}}_{h}\times\Theta_{h}$
such that

$\displaystyle\left(k^{-1}\mathrm{\boldsymbol{u}}_{h},\mathrm{\boldsymbol{v}}_{h}\right)-\left(\nabla\cdot{\mathrm{\boldsymbol{v}}_{h}},p_{h}\right)$
$\displaystyle=\left(\mathrm{\boldsymbol{f}}_{h},\mathrm{\boldsymbol{v}}_{h}\right),\quad$
$\displaystyle\forall\mathrm{\boldsymbol{v}}_{h}\in\mathrm{\boldsymbol{R}}_{h}$
(13)
$\displaystyle\left(\nabla\cdot{\mathrm{\boldsymbol{u}}_{h}},q_{h}\right)$
$\displaystyle=0,\quad$ $\displaystyle\forall q_{h}\in\Theta_{h}$

The resulting saddle point system is then given by

$\begin{bmatrix}M_{h}(k)&B_{h}^{T}\\\
B_{h}&\mathrm{\boldsymbol{0}}\end{bmatrix}\begin{bmatrix}\mathrm{\boldsymbol{u}}_{h}\\\
p_{h}\end{bmatrix}=\begin{bmatrix}\left(\mathrm{\boldsymbol{f}}_{h},\mathrm{\boldsymbol{v}}_{h}\right)\\\
\mathrm{\boldsymbol{0}}\end{bmatrix}$ (14)

We consider the quantity of interest to be the integrated flux across
$\Gamma_{out}\subset\Gamma$

$Q=\dfrac{1}{\left|{\Gamma_{out}}\right|}\int_{\Gamma_{out}}{\mathrm{\boldsymbol{u}}(\cdot,\omega)\cdot\mathrm{\boldsymbol{n}}~{}d\Gamma}$
(15)

##  4 Results

###  4.1 Problem Formulation

We benchmark our approach on a common inference problem in groundwater flow by
estimating the posterior distribution of the mass flux through a domain
boundary. In particular we consider a unit cube discretized using a structured
hexahedral mesh with different spatial resolution at each level in the
hierarchy. In keeping with previous multilevel MCMC studies [16], we consider
a four-level hierarchy, a spatial resolution 16 elements in each direction on
the coarsest level and a uniform refinement factor of two. Hence the hierarchy
consists of 4k, 32.7k, 262k, and 2.1M elements on levels $\ell_{0}$,
$\ell_{1}$, $\ell_{2}$ and $\ell_{3}$, respectively111The exact numbers are
4096, 32768, 262144, and 2097152 elements on levels $\ell_{0}$, $\ell_{1}$,
$\ell_{2}$ and $\ell_{3}$, respectively.. The benchmarking is performed
against reference distributions and statistics obtained using the PDE-only
hierarchy. The hierarchical PDE sampling and the forward Darcy problem is
solved using the scalable, open-source library ParELAGMC, which uses MFEM [28]
to generate the fine grid finite element space, and hypre [29] to implement
the parallel linear algebra. Each simulation is performed using hybridization
algebraic multigrid [26, 30], which is solved with conjugate gradient
preconditioned by hypre’s BoomerAMG. The machine learning capabilities
(training and deployment) are adopted from TensorFlow [31].

We analyze two MLM-augmented hierarchical approach: without filter and with
filtering. A total of five independent chains are generated for each reference
and MLM-augmented approaches and run each chain for 10k samples on the finest
($\ell_{3}$) level. We consider a prior of Gaussian random fields with
correlation length $\lambda=0.3$ and a marginal variance $\sigma^{2}=0.5$. We
generate synthetic observational data
$\mathrm{\boldsymbol{y}}_{obs}\in\mathbb{R}^{25}$ on a reference mesh of 16.7M
elements and assume that
$\mathcal{F}(\mathrm{\boldsymbol{\zeta}})\sim\mathcal{N}(\mathrm{\boldsymbol{y}}_{obs},\sigma_{\eta}I)$
where $I$ is an identity matrix of appropriate dimension and take the noise to
be $\sigma_{\eta}=0.005$; consequently, the likelihood function has the
relation
$\mathcal{L}(\mathrm{\boldsymbol{y}}_{obs}|\mathrm{\boldsymbol{\zeta}})\propto\exp(-\|\mathcal{F}(\mathrm{\boldsymbol{\zeta}})-\mathrm{\boldsymbol{y}}_{obs}\|^{2}/2\sigma_{\eta}^{2})$.

Since the number of parameters and the training cost for fully-connected
neural networks increases rapidly with the dimensionality of the problem, we
employ convolutional neural networks (CNN). For most problem, it is understood
that the size of the network typically grows with the dimensionality of the
mapping, and the size of the training data required typically grow with the
size of the model. Hence, approximating higher levels of the hierarchy with
neural networks becomes infeasible due to the additional cost of training and
memory requirements for generating and storing high-resolution training data.
Therefore, we only apply the surrogate to approximate the maps on the coarsest
level, $\widehat{\mathcal{F}}_{0}$, where $\widehat{\cdot}$ denotes variables
pertaining to the surrogate model. It should be emphasized that although we
employ a surrogate model based on (convolutional) neural networks, alternative
surrogate or reduced order modeling techniques such as Gaussian process models
[32], and linear-subspace [33] and nonlinear manifold [34] reduced order
models can also be used.

The convolutional neural network consists of two convolutional layers, each
with ReLU activation function and 16 and 32 filters, respectively. These
convolutional layers are followed by two dense layers of with 80 and 25 nodes,
respectively. The CNN is trained using 10,000 independently, identically
distributed (i.i.d) samples generated using the PDE samples on the coarsest
level while a further 5,000 samples are used for validation. The training is
performed using ADAM optimizer for 800 epochs with a decaying learning rate
and an additional 200 epochs, with a smaller learning rate are used to fine
tune the CNN model. Figure 1 displays the mean-squared error (MSE) for the
training and validation sets, and suggests a relatively accurate trained model
with MSE on the order of $\mathcal{O}(10^{-3})$.

Figure 1: The convergence in mean-squared error (MSE) for training and
validation set for the training epochs.

###  4.2 Bayesian Inference in Subsurface Flow

The efficiency and accuracy of the MLM-augmented algorithms can be analyze
using metrics such as acceptance rate, number of effective samples, and the
speed-up relative to the reference configuration. Figure 2 shows these metrics
on each level as well as the mean and variance of $Y_{\ell}$ on each level.
Although the mean and variance in $Y_{\ell}$ on the coarsest level (note that
$Y_{0}=Q_{0}$) are comparable for all methods, the variance on the subsequent
coarse level greatly increases when using the unfiltered MCMC. Hence, the
number of effective samples required on the coarsest level for a certain
tolerance ($\epsilon^{2}=0.01$), given by Eq. (10), significantly increases
(Fig. 2(c)). The mean and variance decay of the filtered approach are
comparable to the reference configuration and, as assumed in the MLMC theory,
follow a log-linear trend. Although the acceptance rate (Fig. 2(d)) on the
coarsest level with the unfiltered approach is comparable to that of the
reference configuration, the acceptance rate on the adjacent level decreases
significantly, and indicates that the coarse level sampling with ML forward
map (without filtering) generates poor candidate proposals. The acceptance
rate on the remaining levels recover to those similar to the reference
configuration. The acceptance rate of the filtered algorithm is lower compared
to the other algorithms but the filtering stage filters out these poor
candidates and drastically improves the acceptance rate of the more
computationally expensive second-stage. This improved acceptance rate and the
lower cost-per-sample reduces the computational cost of the complete
hierarchical sampling. Figure 2(e) shows the speed up relative to the
reference configuration for a hierarchy with different number of levels. As
expected, in both cases, as the number of levels increases the overall
efficiency gain decreases, with the unfiltered configuration demonstrating a
faster decay. The four level MCMC sampling using our machine learning-
augmented filtered MCMC approach results in a factor of two speed up over the
standard multilevel PDE hierarchy. Previous work by the authors [16] have
shown a factor of 3.5 speed up over the single level MCMC, hence, our MLM-
augmented approach can deliver a speed of a factor of 7 over the single level
approach.

(a)

(b)

(c)

(d)

(e)

Figure 2: Analysis of chains on each level of a four-level hierarchy showing:
(a) the mean in $Y_{\ell}$, (b) the variance in $Y_{\ell}$, (c) the effective
sample size, and (d) the acceptance ratio. (e) The speed-up relative to the
reference for a hierarchy with different number of levels. The green $\times$
represents the acceptance of the complete two-stage algorithm, whereas green
$\operatorname*{\scalerel*{\cdot}{\bigodot}}$ represents the acceptance rate
of only the second stage.

We define a measure of distance between probability distributions using the
Wasserstein distance. The Wasserstein distance for $Q$ and $Y$ shown in Figs.
3(a) and 3(d) show that the filtering approach yeild PDF _nearer_ to the
reference. Furthermore, we see the distance between PDFs of $Q$ increase with
increasing levels in the hierarchy, while the distance in $Y$ is relatively
constant. It is expected that the distance between $Y$ will decrease with
increasing refinement as the discrete solution tends towards a continuous,
true solution. The large deviation of the posterior PDF of $Y$ on the coarsest
level (Fig. 3(e)) using the unfiltered sampling is clearly evident; the
deviation in PDF of $Q$ around the high density region is also visible. The
inaccuracies, due to the inaccurate coarse-level forward map, leads the
sampling of regions of low posterior densities, which then propogate to the
finer levels as seen in Figs. 3(c) and 3(f).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: The Wasserstein distance (a,d) between the reference posterior
distribution and those obtained using ML-augmented MCMC on each level, the
posterior distributions on the coarsest level (b,e) and the finest level (c,f)
for $Q$ (top) and $Y$ (bottom).

The autocorrelation is indicative of the _mixing_ of the chains, and the
integrated autocorrelation time (IACT) given by Eq. (16) represents the number
of samples required to generate an independent sample. Figure 4 presents the
autocorrelation estimates (according to Eq. (17)) for $Q$ and $Y$ on different
levels. The two-stage algorithm displays slower mixing than the reference, and
hence, suggests that a significant number of samples are being rejected by the
first stage. The second stage displays similar mixing on the coarsest $\ell=0$
level to the other configurations. However, the coarse level sampling using an
inaccurate forward map without filtering leads the chain to areas of low
posterior density and leads to poor mixing on the subsequent finer levels.
This poor mixing propogates to the more computationally-expensive finer level,
and its influence is evident even on the finest level. The filtered approach
displays similar mixing as the reference configuraton on all levels. Table 1
shows the integrated autocorrelation time (IACT) and cost (in seconds) for
generating an independent sample. Although the cost per-independent sample of
the unfiltered approach is significantly lower, the additional number of
samples needed lead to a decrease in efficiency. The cost for the filtered
approach is approximately 1/3 the cost of the reference, and combined with the
accuracy of the filtered approach, leads to an increase in efficiency. It is
important to note that the speed up presented in Fig. 2(e) is for a chain of
fixed length (10k samples on the finest level) and demonstrate a factor of 3.5
and factor 2 speed up due to the unfiltered and filtered approaches,
respectively, over the reference approach. Table 1 presents the theoretical
speed up when considering the number of effective samples needed to tolerance
$\epsilon^{2}=0.01$ (see Appendix A.2) given the cost per independent samples
and variance ($\mathbb{V}[Y_{\ell}]$). Here, we see the divergence of the
coarse level chain in the unfiltered approach leads to an overall decrease in
efficiency, while the filtered approach still yields a theoretical speed of a
factor of 1.2.

(a)

(b)

(c)

(d)

(e)

Figure 4: Autocorrelation estimates of quantity of interest $Q$ (top) and corresponding $Y$ (bottom) on different levels for increasing lag times. In Fig. 4(a), the solid line represent the estimate for only the second stage, while the dashed line represents the estimate of the complete two-stage algorithm. Table 1: Average (over five chains) IACT of $Y$, the cost per independent sample, in seconds, on each level of a four-level hierarchy, total cost and the speed up. Note that $Y_{0}=Q_{0}$; the value in the bracket represents the IACT of the second-stage. | IACT | Cost per Independent Sample | Total Cost | Speed Up
---|---|---|---|---
Models | $\ell_{0}$ | $\ell_{1}$ | $\ell_{2}$ | $\ell_{3}$ | $\ell_{0}$ | $\ell_{1}$ | $\ell_{2}$ | $\ell_{3}$ |  |
Reference | 38 | 5.5 | 3.6 | 3.1 | 1.46 | 0.63 | 1.55 | 14.3 | 49803 | 1.0
Unfiltered | 33 | 49.4 | 4.0 | 3.1 | 0.07 | 4.53 | 3.25 | 14.66 | 75655 | 0.66
Filtered | 77 (39) | 5.4 | 3.4 | 3.1 | 0.91 | 0.58 | 1.49 | 14.08 | 41386 | 1.20

######  Remark 4.

Similar to [16], each $\widehat{Y}_{\ell}$ estimate uses only two levels
$\ell$ and $\ell-1$, whereas [14] uses a hierarchical sampling on all coarser
levels for the estimation.

##  5 Conclusion

This paper introduces an efficient approach for multilevel Markov Chain Monte
Carlo (MCMC), by leveraging a geometric multigrid hierarchy and augmenting the
coarse level with a machine learning model (MLM). By integrating machine
learning model on the coarsest level, we significantly reduce the
computational cost associated with generating training data and training the
model. The MLM allows for a computationally efficient evaluation of coarse
level proposed samples and limits propagation of poor samples to the
computationally more expensive finer levels. Rather than simply replacing a
numerical model with a machine learning model, which introduces approximation
errors, we introduce an additional filtering stage on the coarsest level to
circumvent this approximation error and limit proposal of these poor samples
to the finer levels. We prove that our filtered MCMC algorithm satisfies
detailed balance and is a consistent MCMC sampling algorithm. Furthermore, we
derive conditions on the accuracy of the MLM, compared to the reference model,
to yield more efficient sampling. We demonstrate the practical applicability
of our technique on a large-scale problem in groundwater flow and report a
factor of two speedup over the reference method, while maintaining high levels
of accuracy. The posterior densities and the statistical moments obtained
using our filtered MCMC algorithm were in good agreement with those obtained
using the reference sampler. Furthermore, both our algorithm and the standard
multilevel method exhibit similar decay in mean and variance that align well
with theoretical expectation. These findings underscore the robustness,
efficiency, and theoretical consistency of our proposed multilevel MCMC
approach, offering significant reducation in computational cost and increase
in efficiency.

## Acknowledgements

This work was supported by U.S. Department of Energy’s Advanced Scientific
Computing Research (ASCR) under award 33825/SCW0897 and was performed under
the auspices of the U.S. Department of Energy by Lawrence Livermore National
Laboratory under Contract DE-AC52-07NA27344.

## References

  * [1] Christian P. Robert and Wu Changye.  Markov Chain Monte Carlo Methods, Survey with Some Frequent Misunderstandings.  Wiley, New York, 2021.
  * [2] C.P. Robert, V. Elvira, N. Tawn, and C. Wu.  Accelerating MCMC algorithms.  WIREs Computational Statistics, 10, 2017.
  * [3] L.M. Paun, M. Colebank, M.U. Qureshi, M. Olufsen, N. Hill, and D. Husmeier.  Mcmc with delayed acceptance using a surrogate model with an application to cardiovascular fluid dynamics.  In Proceedings of the International Conference on Statistics: Theory and Applications, Lisbon, Portugal, August 2019.
  * [4] N.T. Hunt-Smith, W. Melnitchouk, F. Ringer, N. Sato, A.W. Thomas, and M.J. White.  Accelerating markov chain monte carlo sampling with diffusion models.  Computer Physics Communications, 296, 2024.
  * [5] Liang Yan and Tao Zhou.  An adaptive surrogate modeling based on deep neural networks for large-scale bayesian inverse problems.  Communications on Computational Physics, 28(5), 2020.
  * [6] H.J. Hortua, R. Volpi, D. Marinelli, and L. Malago.  Accelerating MCMC algorithms through bayesian deep networks.  In Proceedings of the 34th Conference on Neural Information Processing Systems, Canada, December 2020.
  * [7] H.R.B. Orlande, M.J. Colaco, and G.S. Dulikravich.  Approximation of the likelihood function in bayesian technique for the solution of inverse problems.  Inverse Problems in Science and Engineering, 16(6), 2008.
  * [8] X. Meixin, C. Liuhong, M. Ju, and S. Jaemin.  Accelerating the bayesian inference of inverse problems by using data-driven compressive sensing method based on proper orthogonal decomposition.  Electronic Research Archive, 29(5), 2021.
  * [9] V. Charumathi, M. Ramakrishna, and V. Vasudevan.  Fast and accurate proper orthogonal decomposition using efficient sampling and iterative techniques for singular value decomposition.  ACM Tranactions on Mathematical Software, 48(2), 2022.
  * [10] J Andrés Christen and Colin Fox.  Markov chain monte carlo using an approximation.  Journal of Computational and Graphical statistics, 14(4):795–810, 2005.
  * [11] Y. Efendiev, T. Hou, and W. Luo.  Preconditioning Markov chain Monte Carlo simulations using coarse-scale models.  SIAM Journal on Scientific Computing, 28(2):776–803, 2006.
  * [12] A. Golightly, D.A. Henderson, and C. Sherlock.  Delayed acceptance particle MCMC for exact inference in stochastic kinetic models.  Statistics and Computing, 25:1039–2055, 2015.
  * [13] M. Quiroz, M. Tran, M. Villani, and R. Kohn.  Speeding up MCMC by delayed acceptance and data subsampling.  Journal of Computational and Graphical Statistics, 27(1):12–22, 2018.
  * [14] T.J. Dowell, C. Ketelsen, R. Scheichl, and A.L. Teckentrup.  A hierarchical multilevel markov chain monte carlo algorithm with applications to uncertainty quantification in subsurface flows.  SIAM Journal of Uncertainty Quantification, 3:1075–1108, 2015.
  * [15] M.B. Lykkegaard, T.J. Dodwell, C. Fox, G. Mingas, and R. Scheichl.  Multilevel delayed acceptance mcmc.  SIAM Journal of Uncertainty Quantification, 11(1):1–30, 2023.
  * [16] Hillary Fairbanks, Umberto Villa, and Panayot S. Vassilevski.  A multilevel, hierarchical decomposition of finite element white noise with application to multilevel markov chain monte carlo.  SIAM Journal of Scientific Computing, 43(5):S293–S316, 2021.
  * [17] Mikkel B Lykkegaard, Tim J Dodwell, and David Moxey.  Accelerating uncertainty quantification of groundwater flow modelling using a deep neural network proxy.  Computer Methods in Applied Mechanics and Engineering, 383:113895, 2021.
  * [18] Roger G Ghanem and Pol D Spanos.  Stochastic finite elements: a spectral approach.  Courier Corporation, 2003.
  * [19] R. Bassett and J. Derida.  Maximum a posteriori estimators as a limit of bayes estimators.  arXiv:1611.05917v2, 2018.
  * [20] D.M. Blei, A. Kucukelbir, and J.D. McAuliffe.  Variational inference: A review for statisticians.  arXiv:1601.00670v9, 2018.
  * [21] E. Hernandez-Lemus.  Random fields in physics, biology and data science.  Frontiers in Physics, 9:641859, 2021.
  * [22] A.L. Teckentrup, R. Scheichl, M.B. Giles, and E. Ullmann.  Further analysis of multilevel monte carlo methods for elliptic pdes with random coefficients.  Numerische Mathematik, 125(3):569–600, 2013.
  * [23] F. Lindgren, H. Rue, and J. Lindström.  An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial differential equation approach.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(4):423–498, 2011.
  * [24] Sarah Osborn, Panayot S. Vassilevski, and Umberto Villa.  A multilevel, hierarchical sampling technique for spatially correlated random fields.  SIAM Journal of Scientific Computing, 39(5):S5431–S562, 2017.
  * [25] S. Osborn, P. Zulian, T. Benson, U. Villa, R. Krause, and P.S. Vassilevski.  Scalable hierarchical PDE sampler for generating spatially correlated random fields using nonmatching meshes.  Numerical Linear Algebra with Applications, 25(3):e2146, 2018.
  * [26] C.S. Lee and P.S. Vassilevski.  Parallel solver for H(div) problems using hybridization and AMG.  In Domain Decomposition Methods in Science and Engineering XXIII, pages 69–80. Springer, 2017.
  * [27] Jun S. Liu.  Monte Carlo Strategies in Scientific Computing.  Springer, New York, 2008.
  * [28] MFEM: Modular finite element methods library.  mfem.org.
  * [29] HYPRE: High performance preconditioners.  http://www.llnl.gov/CASC/hypre/.
  * [30] V. Dobrev, T. Kolev, C.S. Lee, V. Tomov, and P.S. Vassilevski.  Algebraic hybridization and static condensation with application to scalable H(div) preconditioning.  SIAM Journal on Scientific Computing, 41(3):B425–B447, 2019.
  * [31] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  Software available from tensorflow.org.
  * [32] A. Marrel and B. Iooss.  Probabilistic surrogate modeling by gaussian process: A review on recent insights in estimation and validation.  Reliability Engineering & System Safety, 247:110094, 2024.
  * [33] B.M. de Gooijer, J. Havinga, H.J.M. Geijselaers, and A.H. van den Boogaard.  Evaluation of POD based surrogate models of fields resulting from nonlinear fem simulations.  Advanced Modeling and Simulation in Engineering Sciences, 8:25, 2021\.
  * [34] K. Lee and K.T. Carlberg.  Model reduction of dynamical systems on nonlinear manifolds using deep convolutional autoencoders.  Journal of Computational Physics, 404:108973, 2020.

##  Appendix A Appendix

###  A.1 Integrated Autocorrelation Time

The number of samples between independent samples is determined by the
integrated autocorrelation time (IACT)

$\widehat{\tau}=1+2\sum_{\tau=1}^{M}\widehat{\rho}_{Q}(\tau)$ (16)

where the normalized autocorrelation function is estimated as

$\widehat{\rho}_{Q}(\tau)=\dfrac{1}{N-\tau}\sum_{i=1}^{M}\dfrac{\left(Q^{(i)}-\widehat{\mu}_{Q}\right)\left(Q^{(i+\tau)}-\widehat{\mu}_{Q}\right)}{\widehat{\sigma}^{2}_{Q}}$
(17)

with mean $\widehat{\mu}_{Q}$ and variance $\widehat{\sigma}_{Q}$ of
$\left\\{Q^{(i)}\right\\}_{i=1}^{N}$ and $M\ll N$.

###  A.2 Number of Effective Samples

Consider the effective cost of generating an independent sample on level
$\ell$

$\overline{C_{\ell}}:=\left\lceil\tau_{\ell}\right\rceil\left(C_{\ell}+\sum_{i=1}^{\ell-1}T_{i}C_{i}\right)$
(18)

where $T_{i}=\prod_{j=0}^{i}\left\lceil\tau_{j}\right\rceil$ is the effective
subsampling rate of the coarser levels and the following bound on error of the
multilevel estimator

$\sum_{\ell=0}^{L}\dfrac{\sigma^{2}_{\ell}}{N_{\ell}}\leq\dfrac{\epsilon^{2}}{2}$

where
$\sigma^{2}_{\ell}=\mathbb{V}_{\pi_{\ell},\pi_{\ell-1}}\left[Y_{\ell}\right]$
is the sample variance with respect to the joint distributions $\pi_{\ell}$
and $\pi_{\ell-1}$. The optimization problem for determining the optimal
number of independent samples reads

$\displaystyle\min$
$\displaystyle~{}\sum_{\ell=0}^{L}\overline{C_{\ell}}N_{\ell}$ (19)
$\displaystyle\mathrm{s.t}$
$\displaystyle~{}\sum_{\ell=0}^{L}\dfrac{\sigma^{2}_{\ell}}{N_{\ell}}\leq\dfrac{\epsilon^{2}}{2}$

with the Lagrangian

$\mathbb{L}=\sum_{\ell=0}^{L}\overline{C_{\ell}}N_{\ell}+\lambda\left(\sum_{\ell=0}^{L}\dfrac{\sigma^{2}_{\ell}}{N_{\ell}}-\dfrac{\epsilon^{2}}{2}\right)$

and the stationarity condition

$\dfrac{\partial\mathbb{L}}{\partial
N_{\ell}}=\overline{C_{\ell}}-\lambda\dfrac{\sigma^{2}_{\ell}}{N^{2}_{\ell}}=0\rightarrow\sqrt{\dfrac{\lambda\sigma^{2}_{\ell}}{\overline{C_{\ell}}}}=N_{\ell}$
(20)

Subsituting Eq. (20) into the primal feasibility condition, where the
inequality constraint is treated as a binding, equality constraint

$\sum_{\ell=0}^{L}\dfrac{\sigma^{2}_{\ell}}{N_{\ell}}-\dfrac{\epsilon^{2}}{2}=0$
(21)

yields

$\sum_{\ell=0}^{L}\sqrt{\dfrac{\sigma^{2}_{\ell}\overline{C_{\ell}}}{\lambda}}-\dfrac{\epsilon^{2}}{2}=0\rightarrow\sqrt{\lambda}=\dfrac{2}{\epsilon^{2}}\sum_{\ell=0}^{L}\sqrt{\sigma^{2}_{\ell}\overline{C_{\ell}}}$

Subsituting back into Eq. (20) yields the desired expression for the number of
effective samples as provided in Eq. (10).


# Semi-automatic 3D Object Keypoint Annotation and Detection for the Masses

Kenneth Blomqvist1, Jen Jen Chung1, Lionel Ott1 and Roland Siegwart1  1
Autonomous Systems Lab, ETH Zürich

###### Abstract

Creating computer vision datasets requires careful planning and lots of time
and effort. In robotics research, we often have to use standardized objects,
such as the YCB object set, for tasks such as object tracking, pose
estimation, grasping and manipulation, as there are datasets and pre-learned
methods available for these objects. This limits the impact of our research
since learning-based computer vision methods can only be used in scenarios
that are supported by existing datasets. In this work, we present a full
object keypoint tracking toolkit, encompassing the entire process from data
collection, labeling, model learning and evaluation. We present a semi-
automatic way of collecting and labeling datasets using a wrist mounted camera
on a standard robotic arm. Using our toolkit and method, we are able to obtain
a working 3D object keypoint detector and go through the whole process of data
collection, annotation and learning in just a couple hours of active time.

##  I INTRODUCTION

Most modern computer vision methods use large datasets to learn to predict
features at run time. These have been demonstrated to enable many new
capabilities in robotic object manipulation. While the methods are impressive,
they are data hungry and require sizeable datasets of ground truth annotations
to train. If we could quickly and cheaply create datasets, we could expand to
more environments and enable many downstream tasks.

The data requirements force researchers of downstream robotics tasks to either
use standard objects, for which trained models and computer vision pipelines
have been made available, or a large investment has to be made upfront to
collect and label a dataset. Creating a dataset requires either hand-labeling
thousands of frames one-by-one, having a data collection setup with
environment markers, as done in [1], or a tool such as LabelFusion [2] can be
used to partially automate the annotation process. However, LabelFusion
requires mesh models of the objects. Creating a known model for objects in
turn requires a high-fidelity object scanning setup, which is often
unavailable. It also requires the objects to be rigid, or additional
parameters need to be estimated to model deformation. Additionally, the
objects and environment have to be such that depth sensors are able to
accurately measure depth, excluding reflective or transparent objects.

In this paper, our goal is to track category-level semantic points in an
object’s coordinate frame relative to the camera frame for downstream robotic
manipulation tasks. “Category-level” meaning that objects vary, but the intra-
category semantic meaning of keypoints are the same. Specifically, we want a
system with the following properties:

  1. 1.

Can estimate 3D object keypoints on arbitrary objects

  2. 2.

Requires little effort to handle novel objects

  3. 3.

Can be used in the wild without having to use markers, motion tracking systems
or otherwise modify the environment

  4. 4.

Does not rely on accurate depth sensing

  5. 5.

Can track multiple objects simultaneously in the image frame

Figure 1: StereoLabel, our keypoint labeling tool. The user is presented with
two images of the scene to label. The images are selected to maximize the
orthogonality of the views.

Existing methods, such as PVN3D [3] and kPAM [4] require semantic segmentation
maps to train or they rely on external object instance segmentation. Semantic
segmentation maps are time consuming to annotate, making the systems more
expensive to deploy in new scenarios and for new objects.

In contrast, we present a complete 3D object keypoint tracking system,
including both a learning-based object keypoint algorithm and a method to very
quickly obtain the training labels needed by the algorithm. Our method builds
on the insight that we can forgo using semantic segmentation maps to
distinguish between objects, if we instead introduce a center keypoint and
predict a center map that associates each keypoint with a center keypoint. The
amount of objects in the scene is inferred from the amount of detected center
keypoints. This makes the labeling task a lot faster, as we can simply label
2D keypoints instead of having to also create dense instance segmentation
masks.

We present a way to speed up data collection by capturing many views of the
scene and propagating labels from two labeled viewpoints to all the others. We
show that by calibrating our robot and making use of calibration and the
kinematics of a robot arm, we can forgo using a motion tracking system or
environment markers, as done by previous works. Using our system, data can be
collected in the wild wherever our robot goes. This means that our tools can
be deployed directly on the hardware intended for the downstream robotic task,
streamlining the full problem definition and solution by avoiding additional
steps. Calibrated robots are now commonly available and by using one, the data
collection can be further automated and enables collecting data autonomously.

We show two different versions of our learning-based algorithm that leverages
our data collection pipeline to track keypoints of multiple objects in a
scene. The first one uses both views of a stereo camera. The other one is a
variation of the algorithm that can work with a monocular RGB camera.

We validate our method and tracking pipeline in experiments on two different
object keypoint tracking scenarios. The first one is a single object valve
tracking scenario. The second is a multiple object cup tracking task,
showcasing that we can handle multiple objects simultaneously in a frame. We
show that using only 22.5 minutes of recorded data across 45 sequences, and
using less than 15 minutes of labeling time, we can learn a model that can
track keypoints on objects of interest. We demonstrate that the resulting
tracker is accurate enough to enable manipulation tasks, such as rotating a
valve.

Code for our project is made available at github.com/ethz-
asl/object_keypoints.

##  II RELATED WORK

###  II-A Datasets and Labeling Tools

Several object pose datasets have emerged which use ground truth meshes. The
most commonly used meshes are of the YCB object set [5]. The YCB-video dataset
from [6] provides labeled 6D poses for objects in RGB images. The authors
demonstrated that the dataset was capable of training their PoseCNN 6D object
pose estimator. An initial estimate of the object poses from PoseCNN were used
to generate the YCB-M dataset [7]. This dataset was collected with seven depth
cameras and they used fiducial markers and depth refinement to obtain the
common frame of reference between cameras. While a robot arm was used to
facilitate data collection, the authors did not use the kinematics of the
robot nor did they use hand-eye calibration in the labeling process. Moreover,
both of these datasets are limited to the YCB object set. [8] uses object
models, simulation and rendering to obtain a dataset of ground truth object
poses.

While other labeling tools exist to create datasets of a priori unknown
objects, these often have other limitations. [9] provides a semi-automated
tool for creating 2D and 3D bounding box labels for multi-object scenes in
RGB-D video. Their algorithm uses a GrabCut-based approach [10] to interpolate
annotations over timesteps. However, the user still needs to adjust the
propagated bounding boxes in each of the following frames. LabelFusion [2] can
handle cluttered scenes, however, object meshes are required and must either
be given or created manually using a scanning routine (e.g. with a handheld
scanner or turntable). Our proposed method avoids this requirement altogether.

Finally, several methods train keypoint detectors using only a small set of
labeled data. Simon et al. [11] bootstrap a keypoint detection dataset for
hand pose estimation using multiple views of the scene. The authors ensure
that each iteration introduces new information via multiview geometry.
However, because of this, performance is tied to the number of cameras in the
setup. Multiview geometry is used by [12] for human and animal pose estimation
by deriving a differentiable semi-supervised loss function which is equivalent
to minimizing epipolar divergence. They show that they can train a keypoint
detection network using a large set of unlabeled images and comparatively few
labeled images.

###  II-B Object Keypoints and 6D Pose Estimation

Methods exist which predict keypoints, in 2D or 3D, to calculate the 6D object
pose. Some estimate 2D points on an RGB image and solve for the pose using a
PnP algorithm [13, 14, 15, 16, 17]. Others predict keypoints directly in 3D
space [18, 3]. 6-PACK [19] presents a way to track single objects in real-time
using keypoints which emerge in an unsupervised way. As the keypoints are
learned end-to-end, additional components such as an attention mechanism are
required in their keypoint tracking pipeline. S3K [20] is a self-supervised
approach to learn semantic 3D keypoints. Similar to our approach, the authors
use multiple camera views to propagate labels across images. However, in their
case, they require a four-camera setup while our method is designed to work
with a single camera. NOCS [21] uses a representation shared within an object
category. The authors learn a model to regress to this representation from RGB
and depth maps. PVN3D [3] learns a model which produces semantic segmentation
maps as well as per pixel keypoint and center votes from RGB-D frames.
Similarly to PVN3D, we use a center prediction map to track multiple objects.
However, PVN3D uses ground truth semantic instance segmentation maps to
distinguish objects from each other, which are hard and expensive to label. We
avoid this by associating keypoints directly with their corresponding object’s
center. This also circumvents the need for the expensive clustering step to
aggregate pixel-wise predictions.

KPAM [4, 22] presented a way to track category-level object keypoints.
However, their system can only track single objects due to the integral pose
regression step it relies on [23]. Furthermore, for the same reason, it can’t
deal with many keypoints of the same type. KeyPose [1] is an object keypoint
detection method and dataset, which also uses stereo views of a scene. This
method is only applicable in single object scenes; detected keypoints are not
associated to objects, which makes tracking multiple objects infeasible.
Further, KeyPose only works with objects that have unique keypoints. Modeling
objects such as the valve in our experiments is not possible, as it has
several ambiguous keypoints. This limitation is due to the spatial softmax
operation that is used in the output heatmaps. The dataset collection method
proposed by KeyPose relies on fiduciary tags that are placed in the scan
environment. We propose and demonstrate the feasibility of a method that does
not require modifying the environment and that relies solely on a calibrated
robot with a camera.

##  III METHOD

Figure 2: The components for both of the proposed keypoint tracking pipelines.

Here we describe the components of our framework: the hardware setup, the
procedure to collect video sequences with camera poses, our algorithm to
compute ground truth labels, a stereo multiple object keypoint detection
pipeline and a variation that uses only a monocular RGB camera.

###  III-A Hardware Setup

Our method requires a calibrated image sensor along with a way to control it
into different known viewpoints. For this, we use a StereoLabs ZED Mini stereo
RGB camera, mounted on the wrist of a Franka Emika Panda robot arm. The robot
arm has accurate encoders at each joint which give readings on the position of
each joint. Using a model of the robot and the position readings, we can
accurately compute the position of the wrist, relative to the base frame of
our robot.

We compute the intrinsic parameters and left-to-right-camera transformation of
the stereo camera using Kalibr [24]. The wrist-to-camera frame transformation
we calibrate using the ethz-asl/hand_eye_calibration package [25].

###  III-B Data Collection

For training a neural network to detect keypoints, we need a dataset of image
frames and a set of keypoint locations for each image. To obtain these, we
collect 30s long sequences while our robot scans the target objects, observing
the objects from multiple viewpoints. We save the pose of both camera frames,
relative to the base frame of the robot and the RGB frames from the left and
right camera sensors. In the next section, we describe how we obtain the
keypoints in image coordinates for each image.

###  III-C StereoLabel: Labeling and Generating a Dataset

For each object category of interest, we define a set of keypoints that is
most convenient for manipulating objects from the category. For example, in
the case of coffee cups, we can define the keypoints to be the bottom center,
center top and the outermost point on the handle of the cup (see Fig. 1). This
allows us, if desired, to solve for the orientation of the cup. Or, we can
grasp the cup by approaching the cup from the top center and grasping the side
wall of the cup with a parallel jaw gripper. If there are several ambiguous
keypoints, as is the case for the valve in our experiments (see Fig. 3a), then
all occurrences of the keypoints are labeled and considered to be of the same
type.

We developed a tool, StereoLabel, to label 3D keypoints from a sequence of
images taken from different viewpoints. Fig. 1, shows the tool in use. The
user is shown image frames from two viewpoints. The viewpoints are picked such
that the z-axes of the image frames are as close to perpendicular as possible.
The image frame is defined to be z-axis forward, y down and x to the right of
the image. The user labels 2D keypoints on both frames by clicking on the
keypoint location in the image. In case a specific keypoint is occluded or
otherwise hard to pinpoint, the user can swap out either frame with a new one.

Once corresponding keypoints are labeled, we triangulate their 3D positions in
the base frame of our robot using the homogeneous direct linear transformation
method [26]. We backproject the triangulated points to both frames using each
frame’s projection matrix, so that the user can validate that the point was
appropriately placed. The user can further validate correct placement by
cycling through images in the sequence by pressing a button and checking the
backprojected points. In addition to the labeled keypoints, we augment the set
with one additional 3D keypoint; this is the average of all the other 3D
keypoints which we call the center keypoint.

Once we have all the image sequences labeled and triangulated, we can generate
a dataset for training a computer vision model. Fig. 2 shows the proposed
triangulation-based (stereo) and depth-based (monocular) keypoint tracking
pipelines. For each frame in a sequence, we create a set of ground truth
heatmaps, one heatmap for each type of keypoint. Should there be several
keypoints of the same type, we pack them onto the same heatmap. We compute the
2D image coordinate for each 3D keypoint by backprojection. We place a
Gaussian distribution over the 2D keypoint location computed using an RBF
kernel (output of the prediction step in Fig. 2). Finally, we normalize each
heatmap to have values in the range $[0,1]$.

As there might be multiple objects in a frame, we compute 2D vector fields
with vectors pointing from non-center 2D keypoints to the center keypoint. We
compute one vector for output pixel having a non-zero heatmap value. With the
center maps, we can associate keypoints to objects and detect multiple objects
in a frame.

For the monocular version of our pipeline, we additionally compute a keypoint
depth map containing the z-value of each 3D keypoint for each pixel within a
fixed radius from each keypoint.

###  III-D Learning the Keypoint Network

We use a convolutional neural network (CNN) to predict the heatmaps, along
with center maps. We use CornerNet-Lite [27] as a backbone network. CornerNet-
Lite is a stacked hourglass-style CNN architecture. The input is first
downsampled through a series of convolutional layers and then upsampled
through transposed convolutional layers in an hourglass module. Two hourglass
modules are composed together.

We predict target maps with two prediction modules which take as input the
output of each respective hourglass module. The prediction modules consist of
three convolutional layers with batch normalization, 1$\times$1 kernels with
stride 1 and relu activation functions, except for the last layer. We use
sigmoid activation functions at the heatmap heads and no activation function
for the center map and relu for the depth map.

The input to our network has size 511$\times$511 pixels and the output map
resolution is 64$\times$64\. We initialize the backbone network weights by
pretraining on COCO [28].

We use three types of losses to train our network: a heatmap loss, a center
loss and a depth loss. For the heatmap loss we use binary cross entropy:

$\displaystyle\mathbf{L}_{h}$
$\displaystyle=-\sum_{c=1}^{C}\sum_{i=1}^{H}\sum_{j=1}^{W}y_{cij}\log{p_{cij}}+(1-y_{cij})\log{(1-p_{cij})}.$
(1)

$p_{cij}$ is the predicted heatmap value for a keypoint of type $c$ at output
index $i,j$. $y_{cij}$ is the ground truth heatmap value for keypoint map $c$
at index $i,j$. $C$, $H$ and $W$ denote the amount of keypoint types, and the
height and width of the output maps.

For the center loss, we simply use a smooth L1 loss:

$\displaystyle\mathbf{L}_{c}$
$\displaystyle=\sum_{c=1}^{C}\sum_{i=1}^{H}\sum_{j=1}^{W}\text{smooth\\_L1}(\hat{\mathbf{c}}_{cij}-\mathbf{c}_{cij})\check{y}_{cij},$
(2)

where $\hat{\mathbf{c}}_{cij}$ is the center vector prediction for keypoint
type $c$ at index $i,j$. $\mathbf{c}_{cij}$ is the corresponding ground truth
center vector. $\check{y}_{cij}$ is a binary value denoting whether the
heatmap value for keypoint type $c$ at index $i,j$ is nonzero. The smooth L1
loss is squared below a value of 1 and linear otherwise and is applied
elementwise.

To enable using a monocular camera, we additionally need an estimate of how
far along the z-axis each keypoint is. To do this, we predict a pixelwise
depth estimate for each keypoint type. We learn this using an L1 loss
function:

$\displaystyle\mathbf{L}_{d}$
$\displaystyle=\sum_{c=1}^{C}\sum_{i=1}^{H}\sum_{j=1}^{W}\left\lVert
z_{cij}-\hat{z}_{cij}\right\rVert_{1}\check{y}_{cij},$ (3)

where $z_{cij}$ is the ground truth depth value for keypoint $c$ at location
$i,j$, while $\hat{z}_{cij}$ is the corresponding estimate.

All losses are applied at both stages of the hourglass network and are
combined by weighting parameters $\lambda_{\cdot}$:

$\displaystyle
L=\lambda_{h}(L_{h1}+L_{h2})+\lambda_{c}(L_{c1}+L_{c2})+\lambda_{d}(L_{d1}+L_{d2}).$
(4)

$L_{h1}$ denotes the heatmap loss at the first hourglass, $L_{h2}$ for the
second hourglass, $L_{c1}$ the center loss for the first hourglass and so
forth. When training the triangulation-based pipeline, we set the depth loss
weight $\lambda_{d}$ to 0 to disregard it entirely. We train our network using
the dataset generated in Section III-C.

###  III-E Keypoint Extraction

At runtime, we extract keypoint locations from the heatmaps by first applying
a version of non-maxima supression, where we zero the non-maximum values in
5$\times$5 regions surrounding each location. We then zero out all values
below a threshold of $0.25$. From each of the remaining heatmap values, we
compute keypoint locations by weighing image indices by the predicted heatmap
density in a 5$\times$5 region on the unprocessed heatmap predictions centered
at the maxima location.

For each non-center keypoint, we compute the object center estimate by summing
the center vector with the corresponding image index. We associate each
keypoint with the center keypoint closest to the keypoint’s predicted center
position in pixel coordinates.

###  III-F Keypoint Association and Triangulation

After predicting and extracting keypoints in left and right image frames, we
need to associate each keypoint in the left frame, to its counterpart in the
right frame. To do this, we select the keypoint in the right image where
$\boldsymbol{x^{\prime}}\boldsymbol{F}\boldsymbol{x}$ is below a cutoff value
of $32.0$. $\boldsymbol{F}$ is the fundamental matrix derived from the camera
calibration, $\boldsymbol{x}$ is the homogeneous pixel coordinates of the
keypoint in the left image and $\boldsymbol{x}^{\prime}$ is the homogeneous
pixel coordinates of the keypoint in the right image.

If several keypoints match, which happens when two keypoints are on the
epipolar line, we shift the point by a fixed amount and pick the closest
match. The fixed shift is equivalent to the difference in pixel coordinates
between a point, projected onto both the left and right image frames, that is
60cm in front of the center of the left camera frame.

Finally, we triangulate the 3D location of the keypoints using the same direct
linear transformation method used when creating the dataset.

###  III-G 2D-to-3D

In the monocular version of our pipeline, we use the depth prediction,
combined with the camera matrix $\mathbf{K}$ to compute the 3D point
$\mathbf{X}$ corresponding to the 2D detection $\mathbf{x}$:

$\displaystyle\mathbf{X}=\mathbf{K}^{-1}\mathbf{x}\hat{z},$ (5)

where $\hat{z}$ is the depth estimate for keypoint $\mathbf{x}$.

##  IV EXPERIMENTS

We are interested in the following questions:

  * •

Can we use our method to quickly build up object keypoint tracking datasets?

  * •

Can we train our keypoint tracking method on an amount of data that can be
easily collected by one user?

  * •

Is the object tracking performance good enough to enable robotic manipulation?

###  IV-A Valve: Single Object Tracking

In this experiment, we track a valve with three spokes. We define four
keypoints: one at the center hub of the valve, and three at the front center
points where the spokes meet the rim of the valve, shown in Fig. 3a. The three
keypoints at the rim are indistinguishable from each other, and are thus
considered to be of the same type and packed onto the same heatmap.

We collect 50 sequences of 30s using our data collection method, which we
label using StereoLabel. The sequences differ in object arrangement, clutter,
occlusion, background and lighting conditions. We split the resulting dataset
into 45 sequences for training and 5 sequences for testing.

###  IV-B Label Accuracy

Our semi-automatic labeling approach has a few sources of error:
synchronization between camera frames and joint encoder readings, intrinsics
calibration error and hand-eye calibration error. These all result in some
error in the triangulation and reprojection steps of our pipeline. Without the
ground truth 3D keypoint locations, we instead manually label 2D keypoints
frame-by-frame on 100 randomly sampled frames in our valve dataset and compare
the human labels to ones produced by our system. This allows us to quantify
how much the 2D keypoint labels drift as we observe the target object from
different viewpoints. As the user does not always place the keypoints
perfectly, even the human labels will have some error. We therefore establish
a baseline by doing two manual frame-by-frame labeling passes to measure the
variance.

(a)

(b)

Figure 3: (a) Valve setup showing keypoints for the value. (b) An image from
the cup tracking scene.

###  IV-C KeyPose

We compare our method against KeyPose [1]. KeyPose doesn’t support multiple
objects nor is it possible to detect keypoints on objects with multiple
keypoints of the same type. We therefore can’t run KeyPose on our datasets,
but instead opt to evaluate our method on the KeyPose mugs dataset.

###  IV-D Cups: Multi-object Tracking

In this experiment, we seek to track up to four cups simultaneously in a
scene. We collect a dataset with 100 sequences observing the cups from various
viewpoints, varying the number of cups between 1 and 4 in the scene, changing
the clutter, lighting conditions and background of the scene across sequences.
For each scene, we randomly select between 1 and 4 cups from a set of 25
different cups. We split this dataset into 87 sequences for training and the
rest for testing. We split the sequences such that 2 cups only ever occur in
the test set.

##  V RESULTS

###  V-A Valve: Single Object Tracking

We timed how long it takes to label a sequence of images. Labeling a pair of
valve images took us just under 15 seconds. With 50 sequences, this makes for
a total of $\sim$12 minutes and 30 seconds to label all 19’507 frames in our
dataset.

TABLE I: Results from the valve and cup tracking experiments. Valve
---
Method |  Mean (cm) |  xy (cm) |  $<3$ cm |  25th (cm) |  75th (cm)
GT | 0.39 | 0.18 | 99.3% | 0.17 | 0.30
Stereo | 3.63 | 1.43 | 59.5% | 1.29 | 4.26
Mono | 2.99 | 1.06 | 65.0% | 1.61 | 3.55
Cups
Method |  Mean (cm) |  xy (cm) |  $<3$ cm |  25th (cm) |  75th (cm)
GT | 1.14 | 0.39 | 97.7% | 0.09 | 0.22
Stereo | 6.71 | 2.24 | 68.5% | 1.29 | 3.53
Mono | 3.1 | 1.56 | 62.2% | 1.43 | 4.02

Table I shows the keypoint tracking performance on a held out test set. Mean
refers to the mean error of the 3D keypoints in centimeters. xy is the mean
error, disregarding the depth axis in the left camera frame of reference. $<3$
cm is the percentage of measurements that were within 3 centimeters of the
labeled ground truth location. 25th and 75th respectively denote the 25th and
75th percentiles of the combined keypoint errors. GT refers to the tracking
performance using ground truth heatmaps and center maps as input with the
stereo pipeline, Stereo is the stereo pipeline with a learned model, Mono is
using only the left view of our stereo camera and our monocular pipeline.

Both the stereo- and depth-based pipelines perform reasonably well. For both
pipelines, errors are within the range of the width of a parallel jaw gripper,
and much smaller in scale than the size of the object.

####  V-A1 Valve Manipulation

We deployed our keypoint tracking system on a mobile manipulation system. The
goal of the experiment was to use the system to rotate the valve in Fig. 3a
using the manipulator. In this case, we know the type of the valve and have a
CAD model of it. We first detect the valve using keypoint tracking, and when
we have detected all four keypoints (center and three spokes), we further
refine the pose using ICP to match the depth readings from our camera to the
object model. After refining the object pose, we command the arm to track a
trajectory that rotates the valve. See the supplementary video for a
successful completion of the task.

###  V-B Label Accuracy

Comparing our generated keypoints to a manually labeled dataset, we found that
the mean label difference is 6.3 pixels on average with a standard deviation
of 3.4 on images with a size of 1280x720. Comparing two manually and
separately labeled instances of the same datasets yield a mean difference of
2.9 pixels with a standard deviation of 1.7 pixels. While the manually labeled
examples have slightly less variance, they are of the same order of magnitude
for both methods.

###  V-C KeyPose

TABLE II: Our monocular pipeline on the KeyPose mugs dataset. Method | MAE (cm) |  $<2$ cm
---|---|---
KeyPose | 1.6 | 78.6
Ours monocular | 2.0 | 66.4
Ours Stereo | 1.9 | 69.7

Table II shows results on the KeyPose mugs dataset evaluating on the unseen
$mugs\\_0$ instance. Keypose performs slightly better. We attribute this to
its more restricted problem formulation and the learned stereo image fusion
employed in its network architecture.

###  V-D Cups: Multi-object Tracking

It took us an average of 19 seconds to label a scene with 2 cups. Which makes
for a total of roughly 32 minutes to label the 66’419 frames in our dataset.

Table I shows the accuracy when tracking multiple cups on a held out test set.
Similarly to the valve tracking experiment, the error is larger in the depth
direction. Performance of both the stereo and depth pipelines is acceptable,
i.e. errors are within the width of a parallel jaw gripper. Performance is
slightly worse than on the valve tracking experiment. However, we note that
this is a harder task with several different objects and significantly more
keypoint occlusion.

TABLE III: Time spent on each step of the stereo pipeline. Stage | Mean time (ms)
---|---
Prediction | 32.9
Keypoint extraction | 6.8
Object association | 0.62
Left-to-right association | 0.1
Triangulation | 0.2

Failure modes for both pipelines include misdetecting a keypoint or
associating a keypoint with the wrong object. Failure modes of the stereo
pipeline also include misassociating keypoints from left-to-right and a bad
triangulation due to slightly misdetected 2D keypoints. Additionally, both
approaches fail when two keypoints of the same type align, either from the
same or different objects, and occlude each other. In such cases, the
keypoints will get detected as one and the center prediction might point
toward either object in the case of multiple objects.

Table III shows how long each step of our stereo pipeline takes on average for
our implementation. The measurements were made on a computer with an Nvidia
RTX 2080 GPU and AMD EPYC 7742 CPU.

##  VI DISCUSSION AND CONCLUSIONS

In this paper, we presented a method to quickly collect and label object
keypoint tracking datasets and a system that learns to recover the labels at
runtime on unseen examples. We showed that we fully rely on calibration to
avoid having to place markers in the environment. In experiments, we showed
that we can generate accurate object keypoint labels much quicker than using a
2D labeling approach, while also annotating the z-dimension and without having
to create segmentation maps. We showed that our presented system is able to
detect keypoints on multiple objects simultaneously at real-time rates, using
the produced datasets. We showed that it can be successfully used as part of a
system to solve real world manipulation tasks.

While we presented two different keypoint tracking algorithms, the actual
keypoint and object detection algorithm can be replaced with any other
pipeline, as long as the labels can be derived from our data collection
method. Additional information about the objects could also be used to further
improve the estimated keypoints. The Perspective-n-Point algorithm could be
used for known objects or a category-level object model could be fitted to the
keypoint detections to further improve them, similar to what is used in [29].
Additional computation could be traded for accuracy by predicting keypoint
heatmaps at a higher resolution. The 64x64 pixel output resolution we used is
quite limiting.

One weakness of our data collection method, is that it relies on measurement
timestamps from different sensors. On our platform, these are not hardware
synchronized. Some cameras can be tightly synchronized, while triggering
boards such as the VersaVIS board exist [30] which are able to synchronize
several cameras and IMUs. Extending these to cover other types of sensors,
such as joint encoders, would improve the usability and accuracy of our
proposed method.

Finally, when collecting our data, we manually guide the robot into different
viewpoints. This could be automated. Furthermore, the robot could semi-
autonomously improve upon an initial learned model using an active learning
type approach. Different models and viewpoints could be used to bootstrap a
dataset, similar to what is done in [11].

## References

  * [1] X. Liu, R. Jonschkowski, A. Angelova, and K. Konolige, “KeyPose: Multi-view 3D labeling and keypoint estimation for transparent objects,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 11 602–11 610.
  * [2] P. Marion, P. R. Florence, L. Manuelli, and R. Tedrake, “LabelFusion: A pipeline for generating ground truth labels for real RGBD data of cluttered scenes,” in _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2018, pp. 3235–3242.
  * [3] Y. He, W. Sun, H. Huang, J. Liu, H. Fan, and J. Sun, “PVN3D: A deep point-wise 3D keypoints voting network for 6DoF pose estimation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 11 632–11 641.
  * [4] W. Gao and R. Tedrake, “kPAM-SC: Generalizable manipulation planning using keypoint affordance and shape completion,” _arXiv preprint arXiv:1909.06980_ , 2019.
  * [5] B. Calli, A. Singh, A. Walsman, S. Srinivasa, P. Abbeel, and A. M. Dollar, “The ycb object and model set: Towards common benchmarks for manipulation research,” in _2015 international conference on advanced robotics (ICAR)_. IEEE, 2015, pp. 510–517.
  * [6] Y. Xiang, T. Schmidt, V. Narayanan, and D. Fox, “PoseCNN: A convolutional neural network for 6D object pose estimation in cluttered scenes,” in _Proceedings of Robotics: Science and Systems_ , 2018.
  * [7] T. Grenzdörffer, M. Günther, and J. Hertzberg, “YCB-M: a multi-camera RGB-D dataset for object recognition and 6DoF pose estimation,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 3650–3656.
  * [8] J. Tremblay, T. To, B. Sundaralingam, Y. Xiang, D. Fox, and S. Birchfield, “Deep object pose estimation for semantic robotic grasping of household objects,” _arXiv preprint arXiv:1809.10790_ , 2018.
  * [9] D. Stumpf, S. Krauß, G. Reis, O. Wasenmüller, and D. Stricker, “SALT: A semi-automatic labeling tool for RGB-D video sequences,” in _Proceedings of the 16th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications - Volume 4: VISAPP_ , 2021, pp. 595–603.
  * [10] C. Rother, V. Kolmogorov, and A. Blake, ““GrabCut” - Interactive foreground extraction using iterated graph cuts,” _ACM Transactions on Graphics (TOG)_ , vol. 23, no. 3, pp. 309–314, 2004.
  * [11] T. Simon, H. Joo, I. Matthews, and Y. Sheikh, “Hand keypoint detection in single images using multiview bootstrapping,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2017, pp. 1145–1153.
  * [12] Y. Yao, Y. Jafarian, and H. S. Park, “MONET: Multiview semi-supervised keypoint detection via epipolar divergence,” in _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , 2019, pp. 753–762.
  * [13] B. Tekin, S. N. Sinha, and P. Fua, “Real-time seamless single shot 6D object pose prediction,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 292–301.
  * [14] K. Park, T. Patten, and M. Vincze, “Pix2Pose: Pixel-wise coordinate regression of objects for 6D pose estimation,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2019, pp. 7668–7677.
  * [15] S. Zakharov, I. Shugurov, and S. Ilic, “DPOD: 6D pose object detector and refiner,” in _Proceedings of the IEEE International Conference on Computer Vision_ , 2019, pp. 1941–1950.
  * [16] S. Peng, Y. Liu, Q. Huang, X. Zhou, and H. Bao, “PVNet: Pixel-wise voting network for 6DoF pose estimation,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 4561–4570.
  * [17] Y. Hu, J. Hugonot, P. Fua, and M. Salzmann, “Segmentation-driven 6D object pose estimation,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 3385–3394.
  * [18] S. Suwajanakorn, N. Snavely, J. J. Tompson, and M. Norouzi, “Discovery of latent 3D keypoints via end-to-end geometric reasoning,” in _Advances in Neural Information Processing Systems_ , 2018, pp. 2059–2070.
  * [19] C. Wang, R. Martín-Martín, D. Xu, J. Lv, C. Lu, L. Fei-Fei, S. Savarese, and Y. Zhu, “6-PACK: Category-level 6D pose tracker with anchor-based keypoints,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 10 059–10 066.
  * [20] M. Vecerik, J.-B. Regli, O. Sushkov, D. Barker, R. Pevceviciute, T. Rothörl, C. Schuster, R. Hadsell, L. Agapito, and J. Scholz, “S3K: Self-supervised semantic keypoints for robotic manipulation via multi-view consistency,” _arXiv preprint arXiv:2009.14711_ , 2020.
  * [21] H. Wang, S. Sridhar, J. Huang, J. Valentin, S. Song, and L. J. Guibas, “Normalized object coordinate space for category-level 6D object pose and size estimation,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 2642–2651.
  * [22] L. Manuelli, W. Gao, P. Florence, and R. Tedrake, “kPAM: Keypoint affordances for category-level robotic manipulation,” in _International Symposium on Robotics Research_ , 2019.
  * [23] X. Sun, B. Xiao, F. Wei, S. Liang, and Y. Wei, “Integral human pose regression,” in _Proceedings of the European Conference on Computer Vision (ECCV)_ , 2018, pp. 529–545.
  * [24] L. Oth, P. Furgale, L. Kneip, and R. Siegwart, “Rolling shutter camera calibration,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2013, pp. 1360–1367.
  * [25] F. Furrer, M. Fehr, T. Novkovic, H. Sommer, I. Gilitschenski, and R. Siegwart, _Evaluation of Combined Time-Offset Estimation and Hand-Eye Calibration on Robotic Datasets_. Cham: Springer International Publishing, 2017.
  * [26] A. M. Andrew, “Multiple view geometry in computer vision,” _Kybernetes_ , 2001\.
  * [27] H. Law, Y. Teng, O. Russakovsky, and J. Deng, “CornerNet-Lite: Efficient keypoint based object detection,” _arXiv preprint arXiv:1904.08900_ , 2019\.
  * [28] T.-Y. Lin, M. Maire, S. Belongie, J. Hays, P. Perona, D. Ramanan, P. Dollár, and C. L. Zitnick, “Microsoft COCO: Common objects in context,” in _European Conference on Computer Vision_. Springer, 2014, pp. 740–755.
  * [29] H. Yang, C. Doran, and J.-J. Slotine, “Dynamical pose estimation,” _arXiv preprint arXiv:2103.06182_ , 2021.
  * [30] F. Tschopp, M. Riner, M. Fehr, L. Bernreiter, F. Furrer, T. Novkovic, A. Pfrunder, C. Cadena, R. Siegwart, and J. Nieto, “VersaVIS—an open versatile multi-camera visual-inertial sensor suite,” _Sensors_ , vol. 20, no. 5, p. 1439, 2020.