# Bayesian Phase Estimation via Active Learning

Yuxiang Qiu  Guangdong Provincial Key Laboratory of Quantum Metrology and
Sensing $\&$ School of Physics and Astronomy, Sun Yat-Sen University (Zhuhai
Campus), Zhuhai 519082, China  State Key Laboratory of Optoelectronic
Materials and Technologies, Sun Yat-Sen University (Guangzhou Campus),
Guangzhou 510275, China  Min Zhuang  Guangdong Provincial Key Laboratory of
Quantum Metrology and Sensing $\&$ School of Physics and Astronomy, Sun Yat-
Sen University (Zhuhai Campus), Zhuhai 519082, China  State Key Laboratory of
Optoelectronic Materials and Technologies, Sun Yat-Sen University (Guangzhou
Campus), Guangzhou 510275, China  Jiahao Huang  Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Guangdong Provincial Key
Laboratory of Quantum Metrology and Sensing $\&$ School of Physics and
Astronomy, Sun Yat-Sen University (Zhuhai Campus), Zhuhai 519082, China  State
Key Laboratory of Optoelectronic Materials and Technologies, Sun Yat-Sen
University (Guangzhou Campus), Guangzhou 510275, China  Chaohong Lee  Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>Guangdong Provincial Key
Laboratory of Quantum Metrology and Sensing $\&$ School of Physics and
Astronomy, Sun Yat-Sen University (Zhuhai Campus), Zhuhai 519082, China  State
Key Laboratory of Optoelectronic Materials and Technologies, Sun Yat-Sen
University (Guangzhou Campus), Guangzhou 510275, China

###### Abstract

Bayesian estimation approaches, which are capable of combining the information
of experimental data from different likelihood functions to achieve high
precisions, have been widely used in phase estimation via introducing a
controllable auxiliary phase. Here, we present a non-adaptive Bayesian phase
estimation (BPE) algorithms with an ingenious update rule of the auxiliary
phase designed via active learning. Unlike adaptive BPE algorithms, the
auxiliary phase in our algorithm is determined by a pre-established update
rule with simple statistical analysis of a small batch of data, instead of
complex calculations in every update trails. As the number of measurements for
a same amount of Bayesian updates is significantly reduced via active
learning, our algorithm can work as efficient as adaptive ones and shares the
advantages (such as wide dynamic range and perfect noise robustness) of non-
adaptive ones. Our algorithm is of promising applications in various practical
quantum sensors such as atomic clocks and quantum magnetometers.

††preprint: AIP/123-QED

##  I introduction

Quantum phase estimation is at the core of precision measurement and sensing
[1, 2, 3, 4, 5, 6]. The estimation of an unknown phase via interferometric
techniques are widely used in quantum sensors [7, 8, 9, 10, 11] such as atomic
clocks, magnetometers, and gravimeters. Generally, there are two different
approaches to accomplish the phase estimation problem: frequentist and
Bayesian. Frequentist and Bayesian phase estimation strategies lead to
conceptually different information on the estimated parameters and their
uncertainties according to the results of measurements [12]. Compared to
conventional frequentist estimation approaches [13, 14], the Bayesian approach
is capable to obtain information from every single measurement output. The
Bayesian approach makes use of the Bayes’ theorem to update the posterior
probability, which describes the current knowledge about the random variable
based on the available measurement results. This allows the Bayesian approach
to provide statistical information for any number of measurements.

Bayesian phase estimation (BPE) is known to be particularly efficient and
versatile. In recent, BPE protocol becomes a good choice on account of its
ability to reduce the measurement repeats needed while preserving the
robustness against noises [15, 16, 17, 18]. It is of great value in practical
application where only a limited number of measurements are available [19,
20]. Generally, BPE algorithms can be classified into adaptive (online) [21,
22, 23, 15, 24, 25] and non-adaptive (offline) algorithms [26, 27, 28, 29,
30]. In adaptive BPE algorithms, the auxiliary phases (or other controlled
quantities) are calculated during the process of experiments by taking into
account the previous measurement data. While in non-adaptive BPE algorithms,
the auxiliary phases for each Bayesian update step is pre-determined in
advance. Obviously, adaptive algorithms inevitably require laborious
calculations and operations in order to find step-wise optimized auxiliary
phases. Differently, the auxiliary phases for non-adaptive ones are pre-
determined [28]. Besides, non-adaptive algorithms generally require lower
measurement fidelity, and show better dynamic range and greater consistency in
sensitivity [30]. However, non-adaptive algorithms may require more
measurement times so that the total experimental duration may be much longer.
To benefit from non-adaptive BPE algorithms, can we effectively reduce the
required measurement times?

Machine learning, which involves various algorithms and modeling tools for
data processing, has been widely used in the fields of quantum science and
technologies [31]. It provides a powerful tool for understanding and
exploiting quantum effects, such as classifying many-body quantum phases [32,
33], speeding up many-body quantum simulations [34, 35] and improving the
performances of quantum sensors [17, 24, 36]. In particular, active learning,
which involves human intervention in data preprocessing, is a promising
technique to solve the time- or resources-demanding problems [37, 38]. The
main idea in active learning is that, if a learning algorithm may choose the
data that worth to learn from, it can perform better than traditional methods
with substantially less data. Here, in order to enjoy the advantages of non-
adaptive algorithms and meanwhile reduce the required measurement times, we
propose a BPE algorithm via active learning, which can select significant data
in a pre-learning process to increase the efficiency. Our algorithm is not
only as maneuverable as the non-adaptive algorithms, but also can provide a
desirable measurement precision with a reduced number of measurement times.

In this article, we show how to combine non-adaptive BPE algorithm with active
learning to estimate a phase with reduced measurement times. Compared with
conventional BPE algorithms, our algorithm can save up to $85\%$ measurement
times. Our numerical simulations show that the performances of the error and
the uncertainty versus the Bayesian update times keep the same level achieved
by conventional BPE algorithms. The uncertainty may reach the Ghosh bound of
BPE, which scales as the standard quantum limit. The reduction of the required
measurement times and the pre-determination of the variation of auxiliary
phases makes our algorithm as efficient as the adaptive ones [24]. Moreover,
our algorithm shows good dynamic range and robustness against noises.

Figure 1:  The schematic diagrams for (a) the conventional BPE algorithm and
(b) our BPE algorithm via active learning (labeled as BPE-AL). (a) In the
$n$-th update, the auxiliary phase $\Phi_{n}$ is chosen as $2\pi n/k$, and
then a measurement is implemented. The measurement outcome data $u_{n}$ is fed
for Bayesian update to transform the prior probability into the posterior
probability. Then, the posterior is used as the prior for the next Bayesian
update. For the conventional BPE, the measurement times equals the Bayesian
update times $N$. (b) By using the active learning based on entropy-based
sampling, our BPE-AL only need to select the significant auxiliary phases
$\Phi_{n}$ with $n\in\\{\tilde{n}\\}$ for measurement. Thus, the measurement
times can be reduced to $k^{\prime}N/k$, where $k^{\prime}$ can be much
smaller than $k$.  Figure 2:  (a) The sketch of active learning process in our
BPE-AL algorithm. From $M$ independent repeats of $k$ pre-estimation, the
means $\bar{u}_{n}$ and its stds $\sigma_{n}$ of the resultant data are
calculated. By using the entropy-based sampling, $k^{\prime}$ data with
significant stds are picked out, as indexed by $\tilde{n}_{1}$ to
$\tilde{n}_{k^{\prime}}$ in the orange box. (b) Results of the active learning
process with $k=40$ and $M=20$. Blue points and errorbars correspond to
$\bar{u}_{n}$ and $\sigma_{n}$, respectively. The data with significant stds
are highlighted by orange color, see the orange dashed rectangles. (c) The
mean errors obtained from $50$ independent simulations. Blue points are
results of the conventional BPE algorithm, while orange diamonds stand for the
results of our BPE-AL algorithm with $k=40$ and $k^{\prime}=6$. (d) The
associated stds and uncertainties are plotted by points and dashed lines
respectively. The uncertainties are calculated by the stds of the $N$-th
posterior probability. (e) Scalings of mean errors with our BPE-AL for
$k=40,100,200$ (respectively represented by orange, green and red diamonds)
versus the measurement times $N_{\mathrm{meas}}\approx k^{\prime}N/k$. (f) The
measurement times versus the Bayesian update times $N$. Blue line corresponds
the conventional BPE, while orange, green and red lines respectively
correspond to our BPE-AL with $k=40$, 100 and 200.

##  II Algorithm

###  II.1 General Procedure of Bayesian Phase Estimation

In a general quantum phase estimation procedure, a probe state $|\psi\rangle$
undergoes a transformation to an output state $|\psi(\phi)\rangle$ that
depends on an unknown phase $\phi$. The goal is to estimate $\phi$ according
to the measurement results on $|\psi(\phi)\rangle$. To perform the Bayesian
phase estimation, one may introduce an auxiliary phase $\Phi$ to adjust the
probability distribution. Thus, a measurement outcome $u$ occurs with
probability
$p(u|\phi,\Phi)=\langle\psi(\phi,\Phi)|\hat{\Lambda}_{u}|\psi(\phi,\Phi)\rangle$,
where $\psi(\phi,\Phi)$ is the output state and $\hat{\Lambda}_{u}$ is a POVM
operator [39]. In most situations, the probability $p(u|\phi,\Phi)$ is a
periodic function of the unknown phase $\phi$.

The simplest and most widely used example of quantum phase estimation is the
Ramsey or Mach-Zehnder interferometry with individual two-mode particles [40,
41, 42, 43, 44]. The output state can be expressed as
$|\psi(\phi,\Phi)\rangle=(e^{-i\Phi}|0\rangle+e^{-i\phi}|1\rangle)/\sqrt{2}$
[16, 45, 46, 47, 48, 24, 49, 50], where $|0\rangle$ and $|1\rangle$
respectively represent the two modes. The auxiliary phase $\Phi$ can be well
controlled as desired [51, 28, 15, 24, 52]. After the operation of
recombination, one can finally get the probability of finding the state in
$|u\rangle$ (whose measurement outcome $u$ equals 0 or 1) as [26, 16, 24, 45]:

$p(u|\phi,\Phi)=\frac{1}{2}\left[1+(-1)^{u}\cos{(\phi-\Phi)}\right].$ (1)

In single-shot measurements, a binary outcome $u\in\\{0,1\\}$ can always be
obtained. For an ensemble of $R$ particles ($R>1$), the binary measurement
outcome can be obtained via the threshold method [45, 29, 53]. The threshold
measurement only needs to tell the particle occupation on which state is
quantitatively surpass the particles on the other state, instead of preparing
a system available for single shot measurements. If $r$ particles are found in
$|1\rangle$, then the resultant outcome $u$ is given by

$u=\left\\{\begin{array}[]{l}0,\;\quad\;r\leq R/2,\\\
1,\;\quad\;r>R/2.\end{array}\right.$ (2)

This refers to the “majority voting”, which is widely achieved by analysis of
spin-dependent photo-luminescence data intensity [45, 54] or voltage signals,
or equivalently by repetitive measurements upon a single-spin system, which is
also usually implemented in real experiments.

The BPE algorithms involve a set of measurements and updating the prior
probability distribution according to the Bayes’ theorem [55, 56]. Given the
first $n-1$ outcomes $\vec{u}_{n-1}=(u_{1},u_{2},\ldots,u_{n-1})$, the
posterior probability distribution

$p(\phi|{\vec{u}}_{n})=\mathcal{N}p(u_{n}|\phi,\Phi_{n})p(\phi|\vec{u}_{n-1}),$
(3)

where $p(\phi|\vec{u}_{n-1})$ is the prior probability,
$p(u_{n}|\phi,\Phi_{n})$ is the likelihood and
$\mathcal{N}=[\int{p(u|\phi,\Phi)p(\phi)}]^{-1}$ is a normalization factor.

The auxiliary phase for the $n$-th iteration $\Phi_{n}$ can be designed by
non-adaptive algorithms. For the likelihood function given by Eq.(1), the
variation of the auxiliary phase can be designed with equal steps [27, 28,
30], i.e.,

$\Phi_{n}=n\cdot\frac{2\pi}{k},$ (4)

where $k$ denotes the number of auxiliary phases in a period. This updating
rule of $\Phi_{n}$, which uses pre-established measurement settings or
policies [57, 58, 52], avoids expensive calculations and over-fast
experimental rate [16]. The value of $k$ can be properly chosen within the
permission of the experimental adjusting precision.

As shown in FIG.1(a), the phase estimation procedure for each Bayesian update
can be described as follows.

  * •

Step 1: Given a prior probability distribution $p(\phi|\vec{u}_{n})$ [The
initial $p(\phi|u_{0})$ can be set as a uniform distribution over interval
$[0,2\pi)$].

  * •

Step 2: Perform the experiment with auxiliary phase $\Phi_{n}$ according to
Eq.(4), and record the measurement outcome $u_{n}$.

  * •

Step 3: Update the posterior probability distribution according to Eq.(3).

  * •

Step 4: Evaluate expectation $\phi_{\mathrm{est}}=\int\phi
p(\phi|\vec{u}_{n})d\phi$ and associated uncertainty
$\Delta\phi_{\mathrm{est}}=\sqrt{\int\phi^{2}p(\phi|\vec{u}_{n})d\phi-\phi_{\mathrm{est}}^{2}}$
via the posterior probability $p(\phi|\vec{u}_{n})$ [18, 24, 19, 20].

For the next update, the current posterior probability distribution
$p(\phi|\vec{u}_{n})$ is regarded as a new prior probability. Then return to
Step 1 to start a new cycle until $n=N$.

This iteration makes BPE approaches different from traditional frequentist
approaches [12], and allows better efficiency and noise robustness [15]. The
key issue is the updating rule of auxiliary phase $\Phi$ because the algorithm
cannot get any new useful information when the value of $\Phi$ is fixed.
Adaptive algorithms make it efficiently feasible by designing particular
updating rules or policies of $\Phi$ [22, 57, 30, 16, 24, 52]. While our
algorithm is based upon the above non-adaptive procedure, as shown in the
following, it utilizes the ideas of active learning to reduce the actual
measurement times.

###  II.2 Bayesian Phase Estimation via Active Learning

As a concept from machine learning, active learning [42, 59] involves a
learning algorithm that can choose the data it worth to learn from. It can
perform better than traditional methods with substantially less data. Here we
adopt the so-called _entropy-based sampling_ [60] from active learning to
select the measurement data that significantly affect the results. The entropy
of a discrete probability distribution is defined as [61]
$H=-\sum_{i=1}^{n}\mathrm{P}(x_{i})\log{\mathrm{P}(x_{i})}$, where ${x_{i}}$
refers to all possible values of a random variable $x$, $P(x_{i})$ is the
probability it occurs. The data with larger entropy generally contain more
information [62, 18]. Here, in our algorithm, we introduction a learning
process to select the specific data with large entropy.

In the learning process, we repeat the measurements $M$ times, where each
repeat contains $k$ times of the Bayesian update (3) with $\Phi_{n}=2\pi n/k$
and $n=\\{1,2,\cdots,k\\}$. The measurement times $Mk$ in the learning process
is usually much smaller than the Bayesian update times $N$ in the whole
experiment. The measurement data
$\\{u_{1}^{(m)},u_{2}^{(m)},\ldots,u_{k}^{(m)}\\}_{m=1,2,\ldots,M}$ are
collected and shown in FIG. 1 (b). For the data obtained from the $n$-th
update using $\Phi_{n}$, the associated entropy is given by
$H_{n}=-\sum_{u=0,1}\mathrm{P}_{n}(u)\log{\mathrm{P}_{n}(u)}$. From this
definition, one can find that when $P_{n}(0)=P_{n}(1)=1/2$, the entropy
reaches its maximum. While when $\\{P_{n}(0)~{}=~{}0,~{}P_{n}(1)~{}=~{}1\\}$
or $\\{P_{n}(0)~{}=~{}1,~{}P_{n}(1)~{}=~{}0\\}$, the entropy equals zero. For
the binary measurement outcomes, the standard deviation (std) has similar
property of the entropy. Therefore selecting the phase $\Phi_{n}$ bringing
significant std is equivalent to picking out the data with large entropy. The
std can be calculated by the primal definition:

$\sigma_{n}=\sqrt{\frac{1}{M}\sum_{m=1}^{M}(u_{n}^{(m)}-\bar{u}_{n})^{2}},$
(5)

where the arithmetic mean $\bar{u}_{n}=\frac{1}{M}\sum_{m=1}^{M}u_{n}^{(m)}$.
Our algorithm aims to search for those $\Phi_{n}$ associated with
$\sigma_{n}>0.01$, as shown in FIG. 2 (a).

The above active learning process will make the BPE procedure more efficient.
As shown in FIG. 2 (b), through the active learning, the informative data (the
orange points) with $\sigma_{n}>0.01$ (which are labelled as
$\\{\tilde{n}_{1},\tilde{n}_{2},...,\tilde{n}_{k^{\prime}}\\}$) are selected.
The other data with small std (the blue points) are regarded as uninformative.
Owing to the periodicity, we can deduce that the measurement data obtained
after a Bayesian update using $\Phi_{\tilde{n}}$ are more informative than the
others. Therefore the indices can be explicitly written as

$\\{\tilde{n}\\}=\\{ck+\tilde{n}_{1},ck+\tilde{n}_{2},\ldots,ck+\tilde{n}_{k^{\prime}}\\}_{c=1,2,...}.$
(6)

As shown in the light orange box in FIG. 1 (b), we only perform real
measurements for these $\Phi_{\tilde{n}}$ and record the measurement data
$u_{\tilde{n}}$. While for other $\Phi_{n}$ with $n\notin\\{\tilde{n}\\}$, we
do not perform real measurements and the corresponding outcomes are given as
$0$ or $1$ according to the results of learning.

By including the entropy-based sampling operation, the procedure of our BPE-AL
algorithm is as following.

  * •

Step 0: Acquire the informative set $\\{\tilde{n}\\}$ via the entropy-based
sampling method.

  * •

Step 1: Given a prior probability distribution $p(\phi|\vec{u}_{n})$ [The
initial $p(\phi|u_{0})$ can be set as a uniform distribution over interval
$[0,2\pi)$].

  * •

Step 2: Obtain the data $u_{n}$ via a real measurement with the auxiliary
phase $\Phi_{n}$ only if $n\in\\{\tilde{n}\\}$. Otherwise the measurement is
canceled and the data is given by the results of learning.

  * •

Step 3: Update the posterior probability distribution according to Eq.(3).

  * •

Step 4: Evaluate the expectation $\phi_{\mathrm{est}}$ and the associated
uncertainty $\Delta\phi_{\mathrm{est}}$.

The difference between BPE and BPE-AL only occurs in the Step 2, where we add
a conditional statement to decide whether a real measurement needs to be
implemented. Through picking out the most informative data with limited
measurement times, the BPE-AL algorithm is capable to perform the procedure
more fast. In this way, for the same Bayesian update times, the required
measurement times can be observably smaller than the one for the conventional
BPE. Thus, our algorithm equips the advantages of non-adaptive algorithms and
meanwhile dramatically reduces the real measurement times. In the following,
we will show the performance analysis of our BPE-AL.

##  III Performance analysis

###  III.1 Measurement Times and Measurement Precision

We discuss the measurement times at first. For the learning process, we set
$M=20$, $k=40$, and the estimated phase $\phi=2.7624$ for demonstration. As
shown in FIG.2 (b), the mean $\bar{u}_{n}$ and its std $\sigma_{n}$ are marked
by points and errorbars, respectively. Particularly, the informative data
whose $\sigma_{n}>0.01$ are highlighted by orange color. We find that, in a
period of $\Phi_{n}$, only $k^{\prime}=6$ informative data are selected from
$k=40$ Bayesian updates. Thus, after the learning process, for every period
(which includes $k=40$ Bayesian updates), only $k^{\prime}=6$ real
measurements are needed to performed. For the conventional BPE, the
measurement times equals the Bayesian update times, $N_{\mathrm{meas}}=N$. By
comparison, our BPE-AL algorithm can reduce the measurement times to
$N_{\mathrm{meas}}=Mk+k^{\prime}(N-Mk)/k$. Roughly, if the Bayesian update
times is sufficiently large $N\gg Mk$, the measurement times can be reduced to
$N_{\mathrm{meas}}\approx k^{\prime}N/k$. Thus, the measurement times
performed in experiments can be greatly reduced. For example, for $k=40$,
$k^{\prime}=6$, the measurement times of BPE-AL is reduced to $15\%$. Despite
the measurement times decreases substantially, the Bayesian update times
remains the same, which guarantee the measurement precision.

The comparison of measurement precisions with BPE and BPE-AL are shown in
FIG.2 (c) and (d). We evaluate the absolute error $|\phi_{\mathrm{est}}-\phi|$
and the uncertainty $\Delta\phi_{\mathrm{est}}$. To avoid possible numerical
errors, the phase $\phi$ is discretized to $10^{6}$ points within a period
$[0,2\pi)$. Based upon the results of $50$ independent simulation, we compare
the mean errors $\frac{1}{50}\sum_{k=1}^{50}|\phi_{\mathrm{est}}^{(k)}-\phi|$
and the corresponding stds
$\sqrt{\frac{1}{50}\sum_{k=1}^{50}\left[|\phi_{\mathrm{est}}^{(k)}-\phi|-\left(\frac{1}{50}\sum_{k=1}^{50}|\phi_{\mathrm{est}}^{(k)}-\phi|\right)\right]^{2}}$
for the conventional BPE and our BPE-AL. The mean errors and uncertainties of
our BPE-AL are all as good as the ones of the conventional BPE. This indicates
that our BPE-AL is valid for phase estimation with high-precision as the BPE
does.

Further, we give the simulation results via our BPE-AL with different $k$, as
shown in FIG.2 (e). There are almost no difference between the performances
for different values of $k$ because the decreasing ratios of measurement times
$k^{\prime}/k$ are similar. Under the condition of similar measurement
precision, the measurement times required by our BPE-AL is much smaller
compared with the one of the conventional BPE. The measurement times versus
the Bayesian update times $N$ via BPE and BPE-AL are shown in FIG.2 (f). The
reduction of measurement times can save much experimental time.

###  III.2 Dynamic Range

The dynamic range [10, 29] of a physical quantity is its maximum value that
can be detected with high precision. The ambiguity-free dynamic range of a
phase is fundamentally limited to $2\pi$. Usually, for adaptive BPE
algorithms, the cosine likelihood function (1) can cause inadequately
convolution of the probability distribution functions that may result in bad
performance [15] for $\phi$ near $0$ or $2\pi$. Luckily, our BPE-AL algorithm
preserves the dynamic range of non-adaptive BPE algorithms [30], which works
well for all phases in the interval $[0,2\pi)$. We perform our BPE-AL
algorithm to estimate different unknown phases $\phi$ in the interval
$[0,2\pi)$ and obtain their mean errors and uncertainties from $50$
independent simulations, see FIG. 3(a). For all phases in $[0,2\pi)$, their
mean errors and uncertainties are similar without significant changes. This
indicates that our BPE-AL algorithm can provide effective estimation for a
unknown phase $\phi$ over the whole range of $[0,2\pi)$.

Figure 3:  (a) Mean errors and stds of $50$ independent simulations versus the
estimated phases $\phi$ within $(0,2\pi)$. Blue and orange points with
errorbars correspond to our BPE-AL algorithm with $N/k=100$ and $N/k=1000$,
respectively. (b) Mean errors and stds for $\phi=2.7624$ versus the particle
number $R$ (from $20$ to $200$). The orange dotted dashed line stands for the
measurement times $k^{\prime}$ in each $k=40$ iterations. (c) The squared
uncertainties $\Delta^{2}\phi_{\mathrm{est}}$ with our BPE-AL algorithm (blue
points) versus the measurement times $k^{\prime}N/k$, while green points are
squared uncertainties obtained by GO algorithm. The Ghosh bound is drawn with
orange dashed line, and the standard quantum limit (SQL) is plotted as a black
dashed line. (d) Mean errors of GO algorithm and our BPE-AL algorithm are
shown as green and red dotted lines respectively. Results with our BPE-AL
algorithm are reported in every $k^{\prime}=6$ points to keep the same
consumption of particle number.

###  III.3 Particle Number

The particle number $R$ may influence the binary outcomes through the
threshold criteria (2) and so that it will have an influence on the final
results. Generally, a larger particle number $R$ will result in more precise
binary outcomes. The results from different particle numbers $R$ are shown in
FIG. 3(b) for $N/k=1000$. It is shown that our BPE-AL algorithm is robust
against a lack of particles in each iteration, whose effect can be easily
compensated by increase the number of $k^{\prime}$. These results also suggest
that $R=100$ particles are sufficient to obtain an estimator with good
precision and uncertainty. While increasing $R$ can also decrease $k^{\prime}$
and thus reduce the total measurement times, but it is not recommended
$k^{\prime}<6$ when $k=40$ because we find that a too small $k^{\prime}$ makes
the estimation unstable, which is harmful to the precision.

###  III.4 Uncertainty Bound

Below we discuss the precision bound of our BPE-AL algorithm. Instead of the
Cramer-Rao lower bound often used in frequentist approaches, the precision
bound of our algorithm can be given by the Ghosh bound [63], which is proposed
specifically for Bayesian approaches [12]. The Ghosh bound for an estimated
phase $\phi$ is defined as the reciprocal of the Fisher information [64, 65,
19] of a posterior probability distribution,

$\Delta^{2}\phi_{\rm
GB}(\vec{u}_{n})=\left(\int_{a}^{b}\mathrm{d}\phi\frac{1}{p(\phi|\vec{u}_{n})}\left(\frac{\mathrm{d}p(\phi|\vec{u}_{n})}{\mathrm{d}\phi}\right)^{2}\right)^{-1}.$
(7)

Thus, the Ghosh bound requires
$\Delta^{2}\phi_{\mathrm{est}}(\vec{u}_{n})\geqslant\Delta^{2}\phi_{\rm
GB}(\vec{u}_{n})$, where the estimation uncertainty
$\Delta^{2}\phi_{\mathrm{est}}(\vec{u}_{n})$ is calculated from the variance
of $p(\phi|\vec{u}_{n})$.

The calculation of the Ghosh bound can be implemented iteratively with $n$.
Here, we apply central-limit theorem to approximate the posterior distribution
to a Gaussian distribution to simplify the calculations in Eq.(7). To be
specific, a more strict binomial distribution likelihood function is used to
describe $r$ of $R$ particles found in state $|1\rangle$ as [45]:

$p(r|\phi,\Phi)=\mathrm{C}_{R}^{r}p(1|\phi,\Phi)^{r}p(0|\phi,\Phi)^{R-r},$ (8)

where the two probabilities come from the likelihood function Eq.(1). This
binomial distribution can be replaced by a Gaussian distribution when the
total atom number $R$ is large enough [45]:

$p(r|\phi,\Phi)\approx\frac{1}{\sqrt{2\pi}\sigma}\mathrm{exp}\left[-\frac{{\left(r-Rp(1|\phi,\Phi)\right)}^{2}}{2\sigma^{2}}\right],$
(9)

where $\sigma^{2}=Rp(1|\phi,\Phi)p(0|\phi,\Phi)\approx r(R-r)/R$. When $n>10$
the posterior probability can be approximated as a Gaussian function with
negligible errors [16, 24, 44]. Considering a posterior probability in the
form of Eq. (9), the reciprocal of posterior probabilities in Eq.(7) can be
simplified. The result is shown in FIG.3 (b) as blue points indicating that
the uncertainty of the estimators with our BPE-AL obeys the Ghosh bound pretty
well.

Finally, a simple comparison is made between our BPE-AL algorithm and a
successful adaptive algorithm named the Gaussian optimization (GO) [24]. The
update rule of $\Phi$ in the GO algorithm is designed to keep the variance of
the summation of two likelihood functions Eq.(1) taking its minimum in every
update, resulting in optimal Bayesian updates and estimation efficiency. The
difference is that the update rule of $\Phi$ in GO algorithm is only design
for phase estimation tasks where the binary data $u$ is provided by the
single-particle system (scilicet the $R=1$ case). For this reason in the
comparison every result from GO are simulated by taking the average of $100$
single-paticle measurements. After that the mean errors and stds of $50$
independent simulations are recorded, as well as the uncertainties calculated
in the same way as above. By contrast results from each iteration of our BPE-
AL algorithm are provided in once where the particle number is set as $R=100$.

The comparison results are shown in FIG. 3 (c,d), where the x-axis is marked
by the total measurement times $k^{\prime}N/k$. The uncertainties of GO
algorithm in (c) obeys standard quantum limit (SQL) while above the Ghosh
bound evaluated with our BPE-AL algorithm. In (d) for comparison of estimation
precision, our BPE-AL algorithm reports mean values every $k^{\prime}=6$
points to keep the same consumption of particles. With the same measurement
times, our BPE-AL algorithm can outperform the GO algorithm.

###  III.5 Robustness against Noises

At last, we consider two kinds of noises often occurring in BPE experiments:
depolarization noise and phase noise [24, 52]. Concretely speaking the
depolarization noise caused by errors in measuring apparatus will result in
omitted photon counts, and the phase noise caused by phase fluctuations will
result in random errors in the adjustment of the auxiliary phase $\Phi$. Here
the depolarization and phase noises are simulated by adding Gaussian white
noises on the origin values in Eq.(2) and Eq.(3), which are respectively
characterized by the noise strength parameters $q_{d}$ and $q_{p}$. The values
of readout population number $r$ and auxiliary phase $\Phi$ are then
respectively changed to $r^{\prime}$ and $\Phi^{\prime}$ as:

$\begin{split}r^{\prime}&=(1+\kappa_{d})r,\\\
\Phi^{\prime}&=\Phi+\kappa_{p}\pi,\\\ \end{split}$ (10)

where $\kappa_{d}$ and $\kappa_{p}$ respectively obey the Gaussian
distributions as: $\kappa_{d}\sim\mathcal{N}(0,q_{d}^{2})$ and
$\kappa_{p}\sim\mathcal{N}(0,q_{p}^{2})$.

Simulation results with our BPE-AL algorithm in the presence of noises are
shown in FIG. 4. The measurement times needed in each period $k^{\prime}$ and
the variations of mean errors and stds with noise strengths are given. As the
noise strengths $\kappa_{d}$ and $\kappa_{p}$ increase, the correspondingly
measurement times used in each period $k^{\prime}$ have to be increased.
However, the mean errors and standard deviations do not change a lot when the
strength of noises goes up. This indicates that out BPE-AL algorithm has the
ability to compensate the impact of noises by just increasing the measurement
times. Thus in the presence of noises, it is necessary to select more
informative data to make a desirable estimation via our BPE-AL algorithm.

Figure 4:  The performance in the presence of noises. The measurement times
needed in each period $k^{\prime}$ versus the noise strength under (a)
depolarization noises and (b) phase noises. Mean errors and stds with our BPE-
AL algorithm under (c) depolarization noises and (d) phase noises. Here,
$N/k=1000$ with $k=40$, and the estimated phase is $\phi=2.7624$.

##  IV Discussions

We show how active learning can be used to implement a Bayesian phase
estimation algorithm. Our algorithm involves an active learning process that
can choose the relatively informative data and guide us to perform the
Bayesian phase estimation with substantially less measurement times. For an
example of two-level Ramsey interferometry, we adopt the entropy-based
sampling method from active learning to find the significant auxiliary phases
$\Phi$. After the sampling process, for every $k$ Bayesian updates, only
$k^{\prime}$ (which is much smaller than $k$) measurements are needed to
performed. Thus, in a process of large amounts of Bayesian updates, the real
measurement times needed to perform in experiments can be reduced from $N$ to
$k^{\prime}N/k$ compared with the conventional BPE. In the noiseless case, the
typical measurement times with our algorithm can be $15\%$ of the one with
conventional BPE. While in the presence of noises, by suitably increasing the
measurement times, the performance of our algorithm can remain the same level.
Excepting for reducing the measurement times, our algorithm also has the
advantages of non-adaptive phase estimation algorithms such as perfect dynamic
range and easy accessibility.

The entropy-based sampling method is used in the dataset generation part [59]
of BPE, meaning that reformative designs for BPE are also compatible to our
algorithm. For example, the Markov-chain Monte Carlo (MMC) or particle filter
method [66, 17, 67] and particle guess heuristic (PGH) method [68] can be
combined in, for the reason that we did no changes to the posterior update
part of general BPE algorithms. The particle filter method can be applied by
simply shrinking the distance between the $k^{\prime}$ points in each period
according to the narrowness of posterior functions, which on the other hand
requires better phase adjusting precision. The PGH method can be applied
originally by adjusting the phase accumulation time to achieve better
precision, as it is done in Ref. [16, 15].

In realistic experiments, our algorithm is a promising alternative approach in
real Bayesian parameter estimation tasks to simplify the experiment procedures
such as atom clocks [69] and quantum magnetometers [29], by replacing the
operations of realtime objective function maximization [28, 10, 43, 18, 70,
25], or providing an inherent adequate dynamic range without requirement of
restarts setting in the middle of the estimations [16, 15]. In addition, the
active learning procedure in our algorithm is not relevant to the specific
form of the likelihood function, suggesting that our algorithm can also work
in situation where modeled effect of decoherence time $T_{2}$ are taken into
account [16, 29, 15, 44]. Furthermore, in varying parameters estimation cases
the changing policy of $\Phi$ found by the active learning method can be
optimized simultaneously through reinforcement learning. Our algorithm
provides a promising way for implementing efficient Bayesian phase estimation
in various practical sensors.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(12025509, 11874434), the Key-Area Research and Development Program of
GuangDong Province (2019B030330001), and the Science and Technology Program of
Guangzhou (201904020024). M. Z. is partially supported by the National Natural
Science Foundation of China (12047563). J. H. is partially supported by the
Guangzhou Science and Technology Projects (202002030459).

## data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## References

  * Giovannetti, Lloyd, and Maccone [2006] V. Giovannetti, S. Lloyd, and L. Maccone, “Quantum metrology,” Phys. Rev. Lett. 96, 010401 (2006).
  * Giovannetti, Lloyd, and MacCone [2011] V. Giovannetti, S. Lloyd, and L. MacCone, “Advances in quantum metrology,” Nature Photonics 5, 222–229 (2011).
  * Gross [2012] C. Gross, “Spin squeezing, entanglement and quantum metrology with Bose-Einstein condensates,” Journal of Physics B: Atomic, Molecular and Optical Physics 45 (2012), 10.1088/0953-4075/45/10/103001, 1203.5359 .
  * Degen, Reinhard, and Cappellaro [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, “Quantum sensing,” Rev. Mod. Phys. 89, 035002 (2017).
  * Braun  _et al._ [2018] D. Braun, G. Adesso, F. Benatti, R. Floreanini, U. Marzolino, M. W. Mitchell, and S. Pirandola, “Quantum-enhanced measurements without entanglement,” Rev. Mod. Phys. 90, 35006 (2018).
  * Pezzè  _et al._ [2018] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, “Quantum metrology with nonclassical states of atomic ensembles,” Rev. Mod. Phys. 90 (2018), 10.1103/RevModPhys.90.035005, 1609.01609 .
  * Lane, Braunstein, and Caves [1993] A. S. Lane, S. L. Braunstein, and C. M. Caves, “Maximum-likelihood statistics of multiple quantum phase measurements,” Phys. Rev. A 47, 1667–1696 (1993).
  * Tsang [2012] M. Tsang, “Ziv-zakai error bounds for quantum parameter estimation,” Phys. Rev. Lett. 108, 230401 (2012).
  * Giovannetti and Maccone [2012] V. Giovannetti and L. Maccone, “Sub-heisenberg estimation strategies are ineffective,” Phys. Rev. Lett. 108, 210404 (2012).
  * Waldherr  _et al._ [2012] G. Waldherr, J. Beck, P. Neumann, R. S. Said, M. Nitsche, M. L. Markham, D. J. Twitchen, J. Twamley, F. Jelezko, and J. Wrachtrup, “High-dynamic-range magnetometry with a single nuclear spin in diamond,” Nature Nanotechnology 7, 105–108 (2012).
  * Pezzè, Hyllus, and Smerzi [2015] L. Pezzè, P. Hyllus, and A. Smerzi, “Phase-sensitivity bounds for two-mode interferometers,” Phys. Rev. A 91, 032103 (2015).
  * Li  _et al._ [2018] Y. Li, L. Pezzè, M. Gessner, Z. Ren, W. Li, and A. Smerzi, “Frequentist and Bayesian quantum phase estimation,” Entropy 20, 25–31 (2018).
  * Kay [1993] S. Kay,  _Fundamentals of Statistical Signal Processing: Estimation Theory_ (Prentice Hall: Upper Saddle River, NJ. USA, 1993).
  * Lehmann [1998] G. Lehmann, E.L.; Casella,  _Theory of Point Estimation_ (Springer: Berlin, Germany, 1998).
  * Paesani  _et al._ [2017] S. Paesani, A. A. Gentile, R. Santagati, J. Wang, N. Wiebe, D. P. Tew, J. L. O’Brien, and M. G. Thompson, “Experimental bayesian quantum phase estimation on a silicon photonic chip,” Phys. Rev. Lett. 118, 100503 (2017).
  * Wiebe and Granade [2016] N. Wiebe and C. Granade, “Efficient bayesian phase estimation,” Phys. Rev. Lett. 117, 010503 (2016).
  * Wang  _et al._ [2017] J. Wang, S. Paesani, R. Santagati, S. Knauer, A. A. Gentile, N. Wiebe, M. Petruzzella, J. L. O’brien, J. G. Rarity, A. Laing, and M. G. Thompson, “Experimental quantum Hamiltonian learning,” Nature Physics 13, 551–555 (2017), 1703.05402 .
  * Ruster  _et al._ [2017] T. Ruster, H. Kaufmann, M. A. Luda, V. Kaushal, C. T. Schmiegelow, F. Schmidt-Kaler, and U. G. Poschinger, “Entanglement-based dc magnetometry with separated ions,” Phys. Rev. X 7, 031050 (2017).
  * Rubio, Knott, and Dunningham [2018] J. Rubio, P. Knott, and J. Dunningham, “Non-asymptotic analysis of quantum metrology protocols beyond the cramér–rao bound,” Journal of Physics Communications 2, 015027 (2018).
  * Rubio and Dunningham [2019] J. Rubio and J. Dunningham, “Quantum metrology in the presence of limited data,” New Journal of Physics 21, 043037 (2019).
  * Wiseman [1995] H. M. Wiseman, “Adaptive phase measurements of optical modes: Going beyond the marginal $q$ distribution,” Phys. Rev. Lett. 75, 4587–4590 (1995).
  * Berry and Wiseman [2002] D. W. Berry and H. M. Wiseman, “Adaptive quantum measurements of a continuously varying phase,” Phys. Rev. A 65, 043803 (2002).
  * Armen  _et al._ [2002] M. A. Armen, J. K. Au, J. K. Stockton, A. C. Doherty, and H. Mabuchi, “Adaptive homodyne measurement of optical phase,” Phys. Rev. Lett. 89, 133602 (2002).
  * Lumino  _et al._ [2018] A. Lumino, E. Polino, A. S. Rab, G. Milani, N. Spagnolo, N. Wiebe, and F. Sciarrino, “Experimental phase estimation enhanced by machine learning,” Phys. Rev. Applied 10, 044033 (2018).
  * Dimario and Becerra [2020] M. T. Dimario and F. E. Becerra, “Single-Shot Non-Gaussian Measurements for Optical Phase Estimation,” Phys. Rev. Lett. 125, 120505 (2020).
  * Higgins  _et al._ [2009] B. L. Higgins, D. W. Berry, S. D. Bartlett, M. W. Mitchell, H. M. Wiseman, and G. J. Pryde, “Demonstrating heisenberg-limited unambiguous phase estimation without adaptive measurements,” New Journal of Physics 11, 073023 (2009).
  * Berry and Wiseman [2000] D. W. Berry and H. M. Wiseman, “Optimal states and almost optimal adaptive measurements for quantum interferometry,” Phys. Rev. Lett. 85, 5098–5101 (2000).
  * Said, Berry, and Twamley [2011] R. S. Said, D. W. Berry, and J. Twamley, “Nanoscale magnetometry using a single-spin system in diamond,” Physical Review B - Condensed Matter and Materials Physics 83, 1–7 (2011), 1103.4816 .
  * Nusran, Momeen, and Dutt [2012] N. M. Nusran, M. U. Momeen, and M. V. Dutt, “High-dynamic-range magnetometry with a single electronic spin in diamond,” Nature Nanotechnology 7, 109–113 (2012).
  * Nusran and Dutt [2014] N. M. Nusran and M. V. G. Dutt, “Optimizing phase-estimation algorithms for diamond spin magnetometry,” Phys. Rev. B 90, 024422 (2014).
  * Carleo  _et al._ [2019a] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, “Machine learning and the physical sciences,” Rev. Mod. Phys. 91, 045002 (2019a).
  * Carrasquilla and Melko [2017] J. Carrasquilla and R. G. Melko, “Machine learning phases of matter,” Nature Physics 13, 431–434 (2017).
  * van Nieuwenburg, Liu, and Huber [2017] E. P. L. van Nieuwenburg, Y.-H. Liu, and S. D. Huber, “Learning phase transitions by confusion,” Nature Physics 13, 435–439 (2017).
  * Shen, Liu, and Fu [2018] H. Shen, J. Liu, and L. Fu, “Self-learning monte carlo with deep neural networks,” Phys. Rev. B 97, 205140 (2018).
  * Liu  _et al._ [2017] J. Liu, Y. Qi, Z. Y. Meng, and L. Fu, “Self-learning monte carlo method,” Phys. Rev. B 95, 041101 (2017).
  * Schuff, Fiderer, and Braun [2020] J. Schuff, L. J. Fiderer, and D. Braun, “Improving the dynamics of quantum sensors with reinforcement learning,” New Journal of Physics 22, 035001 (2020).
  * Ding  _et al._ [2020] Y. Ding, J. D. Martín-Guerrero, M. Sanz, R. Magdalena-Benedicto, X. Chen, and E. Solano, “Retrieving quantum information with active learning,” Phys. Rev. Lett. 124, 140504 (2020).
  * Yao  _et al._ [2020] J. Yao, Y. Wu, J. Koo, B. Yan, and H. Zhai, “Active learning algorithm for computational physics,” Phys. Rev. Research 2, 013287 (2020).
  * Wiseman and Milburn [2009] H. M. Wiseman and G. J. Milburn,  _Quantum measurement and control_ (Cambridge university press, 2009).
  * Ramsey [1986] N. Ramsey,  _Molecular Beams_ (Oxford University Press, 1986).
  * Lee  _et al._ [2012] C. Lee, J. Huang, H. Deng, H. Dai, and J. Xu, “Nonlinear quantum interferometry with Bose condensed atoms,” Frontiers of Physics 7, 109–130 (2012).
  * Sattles [2012] B. Sattles,  _Active learning: Synthesis Lectures on Artificial Intelligence and Machine Learning_ (Carnegie Mellon University, 2012).
  * Bonato  _et al._ [2016] C. Bonato, M. S. Blok, H. T. Dinani, D. W. Berry, M. L. Markham, D. J. Twitchen, and R. Hanson, “Optimized quantum sensing with a single electron spin using real-time adaptive measurements,” Nature Nanotechnology 11, 247–252 (2016).
  * Danilin  _et al._ [2018] S. Danilin, A. V. Lebedev, A. Vepsäläinen, G. B. Lesovik, G. Blatter, and G. S. Paraoanu, “Quantum-enhanced magnetometry by phase estimation algorithms with a single artificial atom,” npj Quantum Information 4 (2018), 10.1038/s41534-018-0078-y.
  * Dinani  _et al._ [2019] H. T. Dinani, D. W. Berry, R. Gonzalez, J. R. Maze, and C. Bonato, “Bayesian estimation for quantum sensing in the absence of single-shot detection,” Phys. Rev. B 99, 125413 (2019).
  * Holland and Burnett [1993] M. J. Holland and K. Burnett, “Interferometric detection of optical phase shifts at the heisenberg limit,” Phys. Rev. Lett. 71, 1355–1358 (1993).
  * Lee [2006] C. Lee, “Adiabatic mach-zehnder interferometry on a quantized bose-josephson junction,” Phys. Rev. Lett. 97, 150402 (2006).
  * Cronin, Schmiedmayer, and Pritchard [2009] A. D. Cronin, J. Schmiedmayer, and D. E. Pritchard, “Optics and interferometry with atoms and molecules,” Rev. Mod. Phys. 81, 1051–1129 (2009).
  * Zheng  _et al._ [2020] K. Zheng, M. Mi, B. Wang, L. Xu, L. Hu, S. Liu, Y. Lou, J. Jing, and L. Zhang, “Quantum-enhanced stochastic phase estimation with the su(1,1) interferometer,” Photon. Res. 8, 1653–1661 (2020).
  * Ou and Li [2020] Z. Y. Ou and X. Li, “Quantum su(1,1) interferometers: Basic principles and applications,” APL Photonics 5, 080902 (2020), https://doi.org/10.1063/5.0004873 .
  * Hradil  _et al._ [1996] Z. Hradil, R. Myška, J. Peřina, M. Zawisky, Y. Hasegawa, and H. Rauch, “Quantum phase in interferometry,” Phys. Rev. Lett. 76, 4295–4298 (1996).
  * Rambhatla  _et al._ [2020] K. Rambhatla, S. E. D’Aurelio, M. Valeri, E. Polino, N. Spagnolo, and F. Sciarrino, “Adaptive phase estimation through a genetic algorithm,” Physical Review Research 2, 1–12 (2020).
  * D’Anjou and Coish [2014] B. D’Anjou and W. A. Coish, “Soft decoding of a qubit readout apparatus,” Physical Review Letters 113, 1–5 (2014).
  * Gupta, Hacquebard, and Childress [2016] A. Gupta, L. Hacquebard, and L. Childress, “Efficient signal processing for time-resolved fluorescence detection of nitrogen-vacancy spins in diamond,” J. Opt. Soc. Am. B 33, B28–B34 (2016).
  * Linden, Dose, and Toussaint [2014] W. v. d. Linden, V. Dose, and U. v. Toussaint,  _Bayesian Probability Theory: Applications in the Physical Sciences_ (Cambridge University Press, 2014).
  * Spagnolo  _et al._ [2019] N. Spagnolo, A. Lumino, E. Polino, A. S. Rab, N. Wiebe, and F. Sciarrino, “Machine learning for quantum metrology,” Proceedings 12, 28 (2019).
  * Hentschel and Sanders [2010] A. Hentschel and B. C. Sanders, “Machine learning for precise quantum measurement,” Phys. Rev. Lett. 104, 063603 (2010).
  * Hentschel and Sanders [2011] A. Hentschel and B. C. Sanders, “Efficient algorithm for optimizing adaptive quantum metrology processes,” Physical Review Letters 107, 1–5 (2011).
  * Carleo  _et al._ [2019b] G. Carleo, I. Cirac, K. Cranmer, L. Daudet, M. Schuld, N. Tishby, L. Vogt-Maranto, and L. Zdeborová, “Machine learning and the physical sciences,” Rev. Mod. Phys. 91, 045002 (2019b).
  * Li, He, and Li [2020] L. Li, H. He, and J. Li, “Entropy-based Sampling Approaches for Multi-Class Imbalanced Problems,” IEEE Transactions on Knowledge and Data Engineering 32, 2159–2170 (2020).
  * Shannon [2001] C. E. Shannon, “A mathematical theory of communication,” ACM SIGMOBILE Mobile Comput. Commun. Rev. 5, 3–55 (2001).
  * Gray [2011] R. M. Gray,  _Entropy and Information Theory_ (Springer, 2011).
  * Ghosh [1993] M. Ghosh, “Cramér-Rao bounds for posterior variances,” Statistics & Probability Letters 17, 173–178 (1993).
  * A [1982] H. A,  _Probabilistic and Statistical Aspects of Quantum Theory_ (NorthHolland: Amsterdam, 1982).
  * Paris [2009] M. G. Paris, “Quantum estimation for quantum technology,” International Journal of Quantum Information 7, 125–137 (2009), 0804.2981 .
  * Granade  _et al._ [2012] C. E. Granade, C. Ferrie, N. Wiebe, and D. G. Cory, “Robust online hamiltonian learning,” New Journal of Physics 14, 103013 (2012).
  * Puebla  _et al._ [2020] R. Puebla, Y. Ban, J. F. Haase, M. B. Plenio, M. Paternostro, and J. Casanova, “Versatile atomic magnetometry assisted by bayesian inference,” arXiv (2020), arXiv:2003.02151 .
  * Wiebe  _et al._ [2014] N. Wiebe, C. Granade, C. Ferrie, and D. G. Cory, “Hamiltonian learning and certification using quantum resources,” Phys. Rev. Lett. 112, 190501 (2014).
  * Hume, Rosenband, and Wineland [2007] D. B. Hume, T. Rosenband, and D. J. Wineland, “High-fidelity adaptive qubit detection through repetitive quantum nondemolition measurements,” Phys. Rev. Lett. 99, 1–4 (2007).
  * Dushenko, Ambal, and McMichael [2020] S. Dushenko, K. Ambal, and R. D. McMichael, “Sequential Bayesian Experiment Design for Optically Detected Magnetic Resonance of Nitrogen-Vacancy Centers,” Phys. Rev. Applied 14, 1 (2020).


# A Machine Learning Model for Nowcasting Epidemic Incidence

Saumya Yashmohini Sahai<EMAIL_ADDRESS>Saket Gurukar  Wasiur R. KhudaBukhsh
Srinivasan Parthasarathy  Grzegorz A. Rempała<EMAIL_ADDRESS>Department of
Computer Science and Engineering
The Ohio State University  Division of Biostatistics and Mathematical
Biosciences Institute
The Ohio State University

###### Abstract

Due to delay in reporting, the daily national and statewide COVID-19 incidence
counts are often unreliable and need to be estimated from recent data. This
process is known in economics as nowcasting. We describe in this paper a
simple random forest statistical model for nowcasting the COVID - 19 daily new
infection counts based on historic data along with a set of simple covariates,
such as the currently reported infection counts, day of the week, and time
since first reporting. We apply the model to adjust the daily infection counts
in Ohio, and show that the predictions from this simple data-driven method
compare favorably both in quality and computational burden to those obtained
from the state-of-the-art hierarchical Bayesian model employing a complex
statistical algorithm.

###### keywords:

Nowcasting, Backfilling, COVID-19 Incidence, Random Forest

###### MSC:

[2010] 00-01, 99-00

††journal: Journal of LaTeX Templates

##  1 Introduction

The SARS-CoV-2 virus, first observed in the United States (USA) in January
2020 [1, 2, 3], is highly contagious [4] and has spread in both urban and
rural regions [5, 6] of the USA. To gauge and combat the SARS-CoV-2 spread,
governments and health organizations have set up public information systems
such as COVID-19 dashboards [7, 8, 9, 10]. These dashboards are useful to
brief the public [8] about the current state of COVID-19 in specific regions,
make data-driven public health decisions [10], and improve transparency in
governance [11]. Many of these dashboards show the number of daily new
infections (daily incidence), where the infection count on a particular date
refers to the number of people who started experiencing disease symptoms on
that date (i.e., the onset date of illness). Whereas reporting onset dates is
very useful from the viewpoint of contact tracing and disease spread
monitoring, it is also challenging due to unavoidable delays. [12, 13, 7].
These delays are often due to the time-lags between experiencing initial
symptoms and seeking care, receiving testing results, and updating the
statewide records [14, 15]. As a consequence, the incidence reporting based on
onset counts leads to under-counting of the present and most recent cases.
Dashboards often explicitly warn about this problem [16, 17]. Figure 1 shows
one such example from COVID-19 Dashboard maintained by the Ohio Department of
Health (ODH) [7] where the region of possible under-reporting is marked with a
grey rectangle.

The incomplete current count data poses huge challenges for both local and
national healthcare policymakers as they strive to make difficult public
health decisions (e.g., introduce lockdowns, curfews, evaluate vaccination
effects, etc) in real time to limit the spread of the virus. The use of
statistical methods to moderate the effects of incomplete data could help
reduce uncertainty in public health decision-making during the COVID-19
pandemic and increase public awareness of the most recent disease trends.

Figure 1: Daywise COVID-19 cases in Ohio, as on 12-01-2020. The shaded area –
comprising of 21 days – is the preliminary case data and is likely under-
reported to the ODH due to delayed reporting.

While forecasting COVID-19 cases is typically concerned with predicting the
future burden of the epidemic, nowcasting [18, 19, 20] addresses the problem
of delayed reporting and focuses on the estimation of current case counts from
not-too-distant historic data. Given the under-reported infection data for a
particular date, the nowcasting models estimate the total number of current
infections for that date, which will be reported eventually. In the
literature, there exist several sophisticated statistical methods for
addressing the issue of nowcasting for COVID-19. For instance, Wu et al. [21]
nowcast the probable size of the COVID-19 outbreak in Wuhan, China. The
authors estimate the basic reproduction number $R_{0}$ from their proposed
non-homogeneous counting process modeling the exported number of international
cases from Wuhan and the global human mobility data from/to Wuhan. The authors
then used the estimated $R_{0}$ in the Susceptible-Exposed-Infected-Recovered
or SEIR model [22] for nowcasting and forecasting the outbreak’s size. The
nowcasting problem for delayed reporting of COVID-19 cases is also addressed
by Silva et al. [23] and Greene et al. [13] using Bayesian smoothing approach
[19] where the authors model the delayed number of reported cases with their
proposed Markov counting processes.

In this paper, we propose a simple yet efficient machine learning model that
addresses the problem of nowcasting in a way that is easily understood by non-
experts and therefore suitable for presenting to public health decision-
makers. The only data our proposed model requires can be readily collected
from publicly available dashboards. Despite its simplicity, the model is seen
to predict, with high accuracy (measured with the typical regression-style
$R^{2}$ value), the number of people who start experiencing COVID-19 symptoms
on a particular date. We also show that our proposed model outperforms the
state-of-the-art hierarchical Bayesian model [24] in terms of nowcasting
accuracy while being also approximately 72000x faster. Our model predictions
can also be utilized as input to other forecasting models, for instance, the
ones created for ODH [25] that forecast the future number of infections and
subsequent hospital burden in Ohio.

##  2 Materials and Methods

###  2.1 Data processing

To perform our analysis, we used the public data available at ODH COVID-19
dashboard111https://coronavirus.ohio.gov/wps/portal/gov/covid-19/dashboards/overview,
which is updated daily. It provides the daily partial incidence count, that
is, the count of all individuals $i_{td}$ reported on a given day $t$ to be
confirmed COVID-19 cases with the day of onset $d$ where $d\leq t$. For our
analysis we aggregated cases by the onset date to get the state-level
progression of the onset reporting. Accordingly, the infection count $I_{TD}$
on a specific day $T$ for a given specific onset date $D$ is given by

$I_{TD}=\sum_{d\leq t\leq T}\mathds{1}_{i_{td}},$ (1)

where $\mathds{1}_{i_{td}}$ is the indicator function

$\mathds{1}_{i_{td}}=\begin{cases}1\text{ where }t\leq T,d=D,\\\ 0\text{
otherwise.}\end{cases}$ (2)

Note that for a given $D$, $I_{TD}$ is non-decreasing as a function of $t$
and, assuming that it is also bounded, it has a limit as $T\to\infty$. This is
illustrated in Figure 2 where we see that over the course of 52 days $I_{TD}$
becomes approximately a constant.

Figure 2: Progression of infection count $I_{td}$ for a specific $d$ value
(11-01-2020) over $t$ ranging from 11-01-2020 to 12-22-2020. There is a steep
rise in the infection count in the initial days of data collection as the data
is backfilled, but it gradually stabilizes.

We denote the asymptotic stable value of $I_{TD}$ for an onset date $D$ by

$I^{s}_{D}=\lim_{T\to\infty}I_{TD}$ (3)

and define $F_{TD}$ as the amount of undercounting for a specific $D$ on day
$T$ given by

$F_{TD}=1-\frac{I_{TD}}{I^{s}_{D}}.$ (4)

We may think about $F_{TD}$ as a standardized measure of undercounting that is
also robust to changes in incidence rates during the course of the pandemic.
In what follows, we therefore consider $F_{TD}$ in place of $I_{TD}$. The plot
of $F_{TD}$ as a function of time $T$ is presented in Figure 2. Note that
although in general $F_{TD}\to 0$ as $T\to\infty$, we may see clearly from the
plot that this convergence is not necessarily monotone and that in the fixed
time window $I_{TD}$ only approximately stabilizes as it approaches
$I^{s}_{D}$. In order to improve data stability in the time windows of
interest, we consider the $I_{TD}$ limit to be reached in practice as soon as
$F_{TD}<0.05$. This particular cutoff value was chosen by cross-validation
$[0,0.5]$, as described in Section 2.2.

In order to cross-validate and measure the prediction testing error, data to
be used for nowcasting is split into a training and a validation (testing) set
based on $t$, where all $F_{td}$ with $t<T_{train}$ are in the former and
$t>T_{train}$ are in the latter.

###  2.2 Model

#### Covariates

The model includes the following features to predict the $F_{td}$.

  * 1.

Days since data collection ($\Delta$). For any given infection count
${I_{td}}$ reported on day $t$ with onset date $d$, we define this feature as

$\Delta_{td}=t-d.$ (5)

  * 2.

Day of the week ($\omega_{t}$). This categorical variable denotes the day of
the week for $t$, at which data is being reported, $\omega_{t}\in\\{$Mo, Tu,
We, Th, Fr, Sa, Su$\\}$.

  * 3.

Raw infection count ($I_{td}$). This is the daily partial incidence count for
the pandemic, as described in equation 1.

#### Random forest regression

We train a random forest (RF) regression model [26] on the data partition
defined in section 2.1, to predict $F_{td}$ from the covariates. Formally, we
may write

$F_{td}=f(\Delta_{td},\Delta_{td}^{2},\Delta_{td}^{3},\omega_{i},I_{td}),$ (6)

where $f$ is the RF model.

##  3 Results

#### Goodness of fit

The explained variance ($R^{2}$ value) is used to evaluate the goodness of fit
of the model on both the training data (time window from 10-01-2020 to
11-15-2020) and on the testing data (time window from 11-16-2020 to
12-15-2020). The predictions from the fitted model plotted against the true
values in test data can be seen in Figure 3. The explained variance is $0.99$
on the training data and $0.89$ on the testing data, which shows that the
model’s prediction of $F_{td}$ generalizes well to the unseen data.

Figure 3: Actual vs Predicted $F_{td}$ on the testing dataset. Robust
prediction of $F_{td}$ is crucial for correct prediction of final infection
count $I^{s}_{d}$.

#### Importance of covariates

The relative importance of covariates (the Gini importance or the mean
decrease in impurity) in the fitted model, described in Section 2.2 can be
seen in Table 1. The covariate _days since data collection_ ($\Delta$), along
with its quadratic and cubic transforms turn out to be the most important
features in determining the fraction of missing data $F_{td}$ . The _day of
the week_ $\omega_{i}$ has much less relative importance.

Covariate $\Delta^{2}$ $\Delta^{3}$ $\Delta$ $I$ $\omega_{t}=\text{Th}$
$\omega_{t}=\text{Tu}$ $\omega_{t}=\text{We}$ $\omega_{t}=\text{Fr}$
$\omega_{t}=\text{Mo}$ $\omega_{t}=\text{Sa}$ $\omega_{t}=\text{Su}$
Importance 0.337 0.325 0.311 0.013 0.003 0.003 0.003 0.002 0.001 0.001 0.001

Table 1: Relative (Gini) importance of covariates. Days since data collection
($\Delta$) and its transformations are the most important, with day of the
week $(\omega_{t})$ having the least effect.

#### Prediction of missingness $F_{td}$

Figure 4 shows the prediction of $F_{td}$ for different values of
$\Delta_{td}$. As seen from the plot, the model predictions are close to the
true $F_{td}$ when $\Delta>4$. The good agreement at $\Delta=0$ is trivial, as
at first date of collection, $F_{td}$ is almost always close to $1.0$ and thus
easy to predict. It is also evident that first 3-4 days of data collection
seem to be unreliable in predicting the correct $F_{td}$ and therefore should
be utilized cautiously in the nowcasting predictions.

Figure 4: Predicted missing fraction, $F_{td}$ at various $\Delta_{td}$

#### Actual count prediction

Based on the prediction of $F_{td}$ and the current observed count $I_{td}$,
we use (4), to get the estimate of $I^{s}_{d}$, which is the stable value of
the infection count on day $d$. The typical trends for 4 different days of the
week can be seen in Figure 5. The infection count from the model predicts the
stable value $I^{s}_{d}$ robustly after five days (starting from $\Delta=5$),
and in some cases even earlier. In Figure 5 we may see that irrespective of
the day of the week (Monday, Wednesday, Friday, Sunday), the model is seen to
predict the value of $I^{s}_{d}$ with good accuracy. We may also note that on
Monday and Sunday the model predictions have higher uncertainty likely due to
the effect of weekend test processing slowdown.

Figure 5: Prediction of raw infection count Figure 6: Comparison between the
Bayesian model (BM) and random forest (RF) model from 11-01-2020 to
12-09-2020.  Figure 7: Comparison between the Bayesian model (BM) and random
forest (RF) model from 12-01-2020 to 12-09-2020.

#### Comparison with the Bayesian model

In order to provide some context for assessing the quality of the RF model
predictions, we compare our results with a state-of-the-art hierarchical
Bayesian model proposed recently by Kline et al. [24], which has been used for
the same purpose of nowcasting COVID-19 cases in the state of Ohio. The model,
which we refer to as the Bayesian model (BM) in the following, is more
elaborate than ours as it has also a spatial component. Specifically, it keeps
track of COVID-19 cases over time in different geographical regions (counties
in Ohio). Although in our comparison we aggregate BM spatial counts, for the
sake of completeness we briefly describe here the entire model along with its
spatial component. Denoting by $Y_{i,t}$ the true count of cases in county $i$
with onset date $t$ the BM assumes the following Poisson model for the
dynamics of the disease:

$\displaystyle Y_{i,t}$
$\displaystyle{}\sim\textsc{Poisson}\left(\exp\left(O_{i}+\alpha_{i,t}+X_{t}\eta_{i}\right)\right),$
(7)

where $O_{i}$ is an offset of the logarithm of population of county $i$, the
spatio-temporal random variables $\alpha_{i,t}$ are the latent states of the
process, the design vector $X_{t}$ indicates the day of the week, and the
vector $\eta_{i}$ captures the day of the week effect. It is assumed that
$Y_{i,t}$ is only partially observed for time $t>T_{\text{max}}-D$, where
$T_{\text{max}}$ stands for the last onset date and $D$ (assumed $30$ in [24])
is the maximum reporting delay following onset. BM also uses a semi-local
linear trend model [27] for the spatio-temporal random variables
$\alpha_{i,t}$. Further, the spatial correlation is accounted for using an
intrinsic conditional auto-regressive model. The reporting delay is described
by a Multinomial-Dirichlet model as follows. Denoting by $Z_{i,t,d}$ the count
of cases in county $i$ with onset date $t$, which are observed $d$ days after
$t$, one defines $Z_{i,t}=\left(Z_{i,t,0},Z_{i,t,1},\ldots,Z_{i,t,D}\right)$.
Then, the Multinomial-Dirichlet model prescribes

$\displaystyle Z_{i,t}$
$\displaystyle\sim\textsc{Multinomial}\left(p_{i,t},Y_{i,t}\right),$
$\displaystyle p_{i,t}$
$\displaystyle\sim\textsc{GeneralizedDirichlet}\left(a_{i,t},b_{i,t}\right),$

where the vectors $a_{i,t}$ and $b_{i,t}$ are described in terms of mean and
dispersion parameters [28]. The choice of a Generalized Dirichlet distribution
allows for modeling potential overdispersion in $p_{i,t}$ (see [28]).
Moreover, it leads to a convenient Beta-Binomial conditional distribution
representation for the components $Z_{i,t,d}$. For the purpose of Bayesian
analysis, the authors specify normally distributed priors for the parameters
and use the R package nimble to perform a Markov chain Monte Carlo (MCMC)
algorithm. The authors report a run time of approximately 20 hours for
$30,000$ iterations.

To compare the two models, we calculate the $L_{2}$ distance between the
predictions made by the RF and the Bayesian model, respectively and the actual
known stable values in the Ohio COVID-19 daily counts dataset. We report the
ratio of the two $L_{2}$ distance values as a measure of relative closeness of
the models to the true (stable) data value for days $T-10$ to $T$ and $T-10$
to $T-5$, where $T$ is the last available date in the data. The results are
presented in Table 2. As can be seen in the table, the predictions by the
random forest model are relatively closer to the true values than those
generated using the Bayesian model estimates. The ratio is smaller in the full
10 day window, indicating that the RF model makes better predictions than BM
for days that are close to data collection.

| $T-10:T$ | $T-10:T-5$
---|---|---
$RF/BM$ | 0.565 | 0.726
Table 2: The ratio of $L_{2}$ norm of nowcasted predictions from the Bayesian
model (BM) and random forest model (RF) from the true stable values at two
different time instants $T$ and $T-5$. The ratio values below one indicate
that in both cases the RF model performs better than BM.

##  4 Summary and Discussion

We presented here a simple method for nowcasting COVID-19 cases from historic
data on daily incidence of new cases, as measured by the onset of symptoms.
Such type of data is now widely available for all states in the USA as well as
for most countries in the world. When the need to take immediate decisions on
governance or policy arises, nowcasting can be a useful tool in providing more
accurate estimates about disease incidence and spread. Specifically, our
proposed nowcasting algorithm uses a random forest (RF) regression methodology
and leverages covariates that are based on day of the week, the number of days
passed since first data collection and total incidence so far.

The proposed algorithm is both conceptually simple and computationally
efficient. Our results also suggest that it compares favorably with a much
more elaborate Bayesian model. We have illustrated the application of our
approach on publicly available data from COVID-19 daily onsets in Ohio, as
available from the state’s COVID-19 interactive dashboard. We observed that
the model is able to predict the final incidence for a day, within 3 to 4 days
of data collection. We also find that the number of days passed since first
data collection, along with its transformations (or derivatives), are the most
important covariates in predicting the final incidence.

In order to make our RF method predictions broadly available to interested
researchers and practitioners, we have created a publicly available and
accessible interactive notebook (see below). As described in the repository,
the notebook allows one to use our algorithm to nowcast current COVID-19 onset
occurrences, based on any user-provided historic data supplied in appropriate
format.

The problem of nowcasting historic data is an important one, specially during
the current COVID-19 pandemic, when delays in reporting can snowball into sub-
optimal policies and actions, that can cost lives and create unnecessary
societal burden. Our proposed method allows both general public and health
providers to carefully monitor the pandemic trends and make informed
decisions. The ideas we presented while focused on COVID-19 can be broadly
applicable to similar public health problems in the future.

## Software Availability

The interactive notebook for performing the nowcasting using the random forest
approach described in the paper, along with installation instructions, is
freely available at https://github.com/sahaisaumya/nowcasting

## Acknowledgements

This research was partially funded by NSF grants DMS-1853587 and DMS-2027001
to GAR. The work of WKB was supported by the President’s Postdoctoral Scholars
Program (PPSP) of the Ohio State University. We would like to thank Harley
Vossler for providing helpful feedback on the interactive notebook.

## References

  * [1] T. Bedford, A. L. Greninger, P. Roychoudhury, L. M. Starita, M. Famulare, M.-L. Huang, A. Nalla, G. Pepper, A. Reinhardt, H. Xie, et al., Cryptic transmission of sars-cov-2 in washington state, Science 370 (6516) (2020) 571–575.
  * [2] Centers for Disease Control and Prevention, First Travel-related Case of 2019 Novel Coronavirus Detected in United States, https://www.cdc.gov/media/releases/2020/p0121-novel-coronavirus-travel-case.html (2020).
  * [3] J. R. Fauver, M. E. Petrone, E. B. Hodcroft, K. Shioda, H. Y. Ehrlich, A. G. Watts, C. B. Vogels, A. F. Brito, T. Alpert, A. Muyombwe, et al., Coast-to-coast spread of sars-cov-2 during the early epidemic in the united states, Cell 181 (5) (2020) 990–996.
  * [4] Transmission of SARS-CoV-2: implications for infection prevention precautions, COVID-19 Data Dashboard, https://www.who.int/news-room/commentaries/detail/transmission-of-sars-cov-2-implications-for-infection-prevention-precautions (2020).
  * [5] R. Paul, A. A. Arif, O. Adeyemi, S. Ghosh, D. Han, Progression of covid-19 from urban to rural areas in the united states: a spatiotemporal analysis of prevalence rates, The Journal of Rural Health 36 (4) (2020) 591–601.
  * [6] J. T. Mueller, K. McConnell, P. B. Burow, K. Pofahl, A. A. Merdjanoff, J. Farrell, Impacts of the covid-19 pandemic on rural america, Proceedings of the National Academy of Sciences 118 (1).
  * [7] Ohio Department of Health, COVID-19 Dashboard, https://coronavirus.ohio.gov/wps/portal/gov/covid-19/dashboards (2021).
  * [8] NYC Health, COVID-19: Data, https://www1.nyc.gov/site/doh/covid/covid-19-data.page (2021).
  * [9] California All, Tracking COVID-19 in California, https://covid19.ca.gov/state-dashboard/ (2021).
  * [10] Utah Department of Health, Phased Guidelines for the General Public and Businesses to Maximize Public Health and Economic Reactivation Version 4.5, https://coronavirus-download.utah.gov/Health/Phased_Health_Guidelines_V4.5.3_05262020.pdf (2021).
  * [11] L. Fell, Trust and covid-19: Implications for interpersonal, workplace, institutional, and information-based trust, Digital Government: Research and Practice 2 (1) (2020) 1–5.
  * [12] J. E. Harris, Overcoming reporting delays is critical to timely epidemic monitoring: The case of covid-19 in new york city, MedRxiv.
  * [13] S. K. Greene, S. F. McGough, G. M. Culp, L. E. Graf, M. Lipsitch, N. A. Menzies, R. Kahn, Nowcasting for real-time covid-19 tracking in new york city: An evaluation using reportable disease data from early in the pandemic, JMIR public health and surveillance 7 (1) (2021) e25538.
  * [14] The Wall Street Journal, Covid-19 Data Reporting System Gets Off to Rocky Start, https://www.wsj.com/articles/covid-19-data-reporting-system-gets-off-to-rocky-start-11597178974 (2020).
  * [15] Governer of Ohio, COVID-19 Update: Antigen Testing, K-12 Education Update, DataOhio Portal, https://governor.ohio.gov/wps/portal/gov/governor/media/news-and-media/covid19-update-12072020 (2020).
  * [16] World Health Organization, WHO Coronavirus Disease (COVID-19) Dashboard, https://covid19.who.int/ (2021).
  * [17] Washington State Department of Health, COVID-19 Data Dashboard, https://www.doh.wa.gov/Emergencies/COVID19/DataDashboard (2021).
  * [18] J. van de Kassteele, P. H. Eilers, J. Wallinga, Nowcasting the number of new symptomatic cases during infectious disease outbreaks using constrained p-spline smoothing, Epidemiology (Cambridge, Mass.) 30 (5) (2019) 737.
  * [19] S. F. McGough, M. A. Johansson, M. Lipsitch, N. A. Menzies, Nowcasting by bayesian smoothing: A flexible, generalizable model for real-time epidemic tracking, PLoS computational biology 16 (4) (2020) e1007735.
  * [20] J. Lawless, Adjustments for reporting delays and the prediction of occurred but not reported events, Canadian Journal of Statistics 22 (1) (1994) 15–31.
  * [21] J. T. Wu, K. Leung, G. M. Leung, Nowcasting and forecasting the potential domestic and international spread of the 2019-ncov outbreak originating in wuhan, china: a modelling study, The Lancet 395 (10225) (2020) 689–697.
  * [22] J. L. Aron, I. B. Schwartz, Seasonality and period-doubling bifurcations in an epidemic model, Journal of theoretical biology 110 (4) (1984) 665–679.
  * [23] A. A. M. d. Silva, L. G. Lima-Neto, L. M. M. d. Costa, M. L. B. M. Bragança, A. K. D. Barros Filho, B. B. Wittlin, B. F. d. Souza, B. L. C. A. d. Oliveira, C. A. d. Carvalho, E. B. A. F. Thomaz, et al., Population-based seroprevalence of sars-cov-2 and the herd immunity threshold in maranhão, Revista de Saúde Pública 54 (2020) 131.
  * [24] D. Kline, A. Hyder, E. Liu, M. Rayo, S. Malloy, E. Root, A Bayesian spatio-temporal nowcasting model for public health decision-making and surveillance, arXiv preprint arXiv:2102.04544.
  * [25] Infectious Disease Institute (IDI) COVID-19 Response Modeling Team at The Ohio State University, Predicting COVID-19 Cases and Subsequent Hospital Burden in Ohio, https://idi.osu.edu/assets/pdfs/covid_response_white_paper.pdf (2021).
  * [26] L. B. Statistics, L. Breiman, Random forests, in: Machine Learning, 2001, pp. 5–32.
  * [27] K. H. Brodersen, F. Gallusser, J. Koehler, N. Remy, S. L. Scott, Inferring causal impact using Bayesian structural time-series models, The Annals of Applied Statistics 9 (1) (2015) 247 – 274.
  * [28] O. Stoner, T. Economou, Multivariate hierarchical frameworks for modeling delayed reporting in count data, Biometrics 76 (3) (2020) 789–798.