# James-Stein Estimation in Guassian Metrology : V2

Wilfred Salmon  DAMTP, Centre for Mathematical Sciences, University of
Cambridge, United Kingdom  Hitachi Cambridge Laboratory, J. J. Thomson Avenue,
CB3 0HE, Cambridge, United Kingdom  Sergii Strelchuk  DAMTP, Centre for
Mathematical Sciences, University of Cambridge, United Kingdom  Department of
Computer Science, University of Warwick, United Kingdom  David R. M.
Arvidsson-Shukur  Hitachi Cambridge Laboratory, J. J. Thomson Avenue, CB3 0HE,
Cambridge, United Kingdom

# James-Stein Estimation in Quantum Gaussian Sensing

Wilfred Salmon  DAMTP, Centre for Mathematical Sciences, University of
Cambridge, United Kingdom  Hitachi Cambridge Laboratory, J. J. Thomson Avenue,
CB3 0HE, Cambridge, United Kingdom  Sergii Strelchuk  DAMTP, Centre for
Mathematical Sciences, University of Cambridge, United Kingdom  Department of
Computer Science, University of Warwick, United Kingdom  David R. M.
Arvidsson-Shukur  Hitachi Cambridge Laboratory, J. J. Thomson Avenue, CB3 0HE,
Cambridge, United Kingdom

###### Abstract

The James-Stein estimator is a biased estimator—for a finite number of samples
its expected value is not the true mean. The maximum-likelihood estimator
(MLE), is unbiased and asymptotically optimal. Yet, when estimating the mean
of $3$ or more normally-distributed random variables, the James-Stein
estimator has a smaller total (expected) error than the MLE. We introduce the
James-Stein estimator to the field of quantum metrology, from both the
frequentist and Bayesian perspectives. We characterise the effect of quantum
phenomena on the James-Stein estimator through the lens of quantum Gaussian
sensing, the task of estimating the mean of an unknown multivariate quantum
Gaussian state. We find that noiseless entanglement or coherence improves
performance of the James-Stein estimator, but diminishes its advantage over
the MLE. In the presence of noise, the James-Stein advantage is restored.
Quantum effects can also boost the James-Stein advantage. We demonstrate this
by investigating multivariate postselective metrology (generalised weak-value
amplification), a strategy that uses quantum effects to measure parameters
with imperfect detectors. Simply by post-processing measured data differently,
our techniques reduce errors in quantum experiments.

##  I Introduction

The James-Stein estimator underlies one of the most counterintuitive phenomena
in classical statistics [1, 2]. Suppose one wishes to estimate 3 or more
independent normally distributed quantities, minimising the total mean-squared
error. The most common strategy for such a task is the maximum-likelihood
estimator (MLE), which treats each quantity separately. However, Stein [1]
showed the surprising result that it is always better to combine the data for
each quantity, although the underlying quantities are completely independent.
Such estimation strategies are now referred to as James-Stein estimators [2].
We refer to the advantage of the James-Stein estimator over the MLE as the
James-Stein advantage.

In this Article, we consider the application of James-Stein estimation to
parameters encoded in quantum systems. We focus on the particular case of
Gaussian shift models [3], which are the quantum generalisations of a Gaussian
distribution with an unknown mean. The central-limit theorem [4] ensures that
the mean of many identical, independent measurements of any parameterised
state is approximately Guassian. Moreover, quantum local asymptotic normality
[3] implies that many copies of any parameterised quantum state may be locally
approximated as a Gaussian shift model. Therefore, the study of Gaussian shift
models is ubiquitous in quantum metrology. We analyse the James-Stein
advantage in commonly studied areas of quantum metrology: noiseless metrology,
noisy metrology, Bayesian metrology and metrology with postselection (all
defined below). We find that boosting metrology by noiseless techniques,
employing entanglement or coherence, generally leads to a decrease in the
James-Stein advantage. In practice, noise precludes the achievement of
theoretically optimal quantum improvements [5]. We show that in the presence
of noise, the James-Stein advantage can be considerable. We demonstrate these
results in the frequentist and Bayesian regimes. We also consider other
experimental limitations, such as imperfect detectors. In such scenarios, it
is advantageous to limit the number of measurements of quantum states via
postselective filters (see below). These filters harness use quantum effects
[6] to compress the information from many copies of a parameterised quantum
state, into (arbitrarily) fewer copies. By utilising such compression, one
requires few measurements of quantum states to estimate an unknown parameter.
We construct an explicit iterative scheme for postselection of Gaussian
states. Analysing this scheme, we find that the James-Stein advantage can be
increased by postselective filtering. Moreover, postselective filtering can be
improved by James-Stein estimation.

##  II Preliminaries

###  II.1 Frequentist parameter estimation

Consider an unknown $n$-dimensional parameter $\theta\in\Theta$, where
$\Theta\subseteq\mathbb{R}^{n}$ is a known parameter space. To estimate
$\theta$, one can sample a random variable $X$ whose distribution $P_{\theta}$
has a probability density function (pdf) $f(x|\theta)$. If $X$ takes the value
$x$, one denotes the corresponding estimate of $\theta$ by
$\hat{\theta}^{\text{}}(x)$. The function $\hat{\theta}^{\text{}}$ is called
an estimator. The performance of an estimator is quantified by its risk:

$R(\hat{\theta}^{\text{}},\theta)=\mathbb{E}_{X\sim
P_{\theta}}\left[\norm{\hat{\theta}^{\text{}}(X)-\theta}^{2}\right].$ (1)

Using the risk, one can quantify the advantage of one estimator over another
at a fixed value of the parameter $\theta$:

$\text{AD}(\hat{\theta}^{\text{1}},\hat{\theta}^{\text{2}},\theta)=\frac{R(\hat{\theta}^{\text{2}},\theta)}{R(\hat{\theta}^{\text{1}},\theta)}.$
(2)

If
$\text{AD}(\hat{\theta}^{\text{1}},\hat{\theta}^{\text{2}},\theta^{\prime})>1$,
it is beneficial to use $\hat{\theta}^{\text{1}}$ over
$\hat{\theta}^{\text{2}}$ when $\theta=\theta^{\prime}$.

In the case of $N$ separate observations $X=(X_{1},\dots,X_{N})$, one
considers a family of estimators $\hat{\theta}^{\text{}}_{N}$; one for each
value of $N$. The sequence $R(\hat{\theta}^{\text{}}_{N},\theta)$ quantifies
how well an estimation strategy $\hat{\theta}^{\text{}}_{N}$ performs with
increasing resources. The most commonly studied regime in parameter estimation
is the asymptotic limit, where $\theta$ is fixed and $N\to\infty$ [7, 8, 9].
In rough terms, the “best” possible estimator $\hat{\theta}^{\text{OPT}}_{N}$
has risk scaling as

$NR(\hat{\theta}^{\text{OPT}}_{N},\theta)\to\text{Tr}[F(\theta)^{-1}],$ (3)

where $F(\theta)$ is the Fisher information matrix [10], given by

$\displaystyle F(\theta)_{ij}$ $\displaystyle=$
$\displaystyle\mathbb{E}_{X\sim
P_{\theta}}\left[\frac{\partial}{\partial\theta_{i}}\log
f(X|\theta)\frac{\partial}{\partial\theta_{j}}\log f(X|\theta)\right].$ (4)

In the quantum setting [11], a parameter $\theta$ is encoded in some quantum
state $\rho(\theta)$. Often, one has $N$ copies of the same state
$\sigma(\theta)$, in which case $\rho(\theta)=\sigma(\theta)^{\otimes N}$. One
can extract information by measuring $\rho(\theta)$. This involves choosing a
POVM [12], $M=\\{M_{k}\\}$, where $M_{k}\geq 0$, and
$\sum_{k}M_{k}=\mathbbm{1}$. The measurement and state induce a probability
distribution parametrised by $\theta$:

$p_{k}^{(M)}(\theta)=\text{Tr}[\rho(\theta)M_{k}].$ (5)

Once a measurement is chosen, the quantum-parameter-estimation problem has
been reduced to a classical problem (where $X$ takes values in the set of
possible measurement outcomes). The quantum analogue of equation (3) is given
by

$NR(\hat{\theta}^{\text{qOPT}}_{N},\theta)\to C^{H}(\theta),$ (6)

where $C^{H}$ is called the Holevo-Cramer-Rao bound; see Ref. [13] for
details. Note that the optimal estimation strategy
$\hat{\theta}^{\text{qOPT}}_{N}$ includes optimisation over both measurements
and estimators.

###  II.2 The James-Stein Estimator

A particularly important parameter-estimation problem concerns data that is
normally distributed around an unknown parameter:
$Z\mid\theta\sim\mathcal{N}(\theta,\Sigma)$111Here, the notation $Z|\theta$
denotes a random variable $Z$ for a given value of $\theta$. Moreover, $\sim$
means “distributed as”.. Here, $\mathcal{N}(\theta,\Sigma)$ denotes the
(multivaritate) Gaussian (normal) distribution [14], with mean
$\theta\in\mathbb{R}^{n}$ and covariance matrix $\Sigma$. For simplicity we
assume that $\Sigma$ is known, though this assumption can be dropped [15].

For normally distributed data, $Z\mid\theta\sim\mathcal{N}(\theta,\Sigma)$, an
intuitive estimator is the MLE, which aligns the estimate with the observed
data: $\hat{\theta}^{\text{MLE}}(z)=z$. The MLE’s risk is
$R(\hat{\theta}^{\text{MLE}},\theta)=\text{Tr}(\Sigma)=\text{Tr}[F(\theta)^{-1}]$.
Below, we compare the MLE with two types of James-Stein estimators [16]. The
first is parameterised by a fixed vector $\nu\in\mathbb{R}^{n}$ and is given
by

$\hat{\theta}^{\text{{\hbox{\nu}}-JS}}(z)=z-\frac{(n-2)\Sigma^{-1}}{(z-\nu)^{T}\Sigma^{-2}(z-\nu)}(z-\nu),$
(7)

where $n\geq 3$. $\hat{\theta}^{\text{$\nu$-JS}}(z)$ shrinks the Gaussian
random variable $Z$ towards the vector $\nu$. The risk of
$\hat{\theta}^{\text{$\nu$-JS}}(z)$ is

$\displaystyle R(\hat{\theta}^{\text{$\nu$-JS}},\;$
$\displaystyle\theta)=\text{Tr}(\Sigma)$
$\displaystyle-(n-2)^{2}\mathbb{E}\left[\frac{1}{(Z-\nu)^{T}\Sigma^{-2}(Z-\nu)}\right],$
(8)

where the expectation is taken over
$Z\mid\theta\sim\mathcal{N}(\theta,\Sigma)$ [16]. The second term is strictly
negative. Thus,
$\text{AD}(\hat{\theta}^{\text{JS}},\hat{\theta}^{\text{MLE}},\theta)>1$ for
any value of $\theta$. The advantage of the James-Stein estimator is biggest
when $\theta$ is close to $\nu$ and $\Sigma$ is “large”. Below, we fix
$\nu=0$, and let $\hat{\theta}^{\text{JS}}\equiv\hat{\theta}^{\text{0-JS}}$.
If one expects $Z$ to be close to $\nu\neq 0$ (and hence that
$\hat{\theta}^{\text{$\nu$-JS}}$ will perform better than
$\hat{\theta}^{\text{JS}}$), one can simply translate the origin to $\nu$.

Instead of shrinking towards a fixed vector, the second type of James-Stein
estimator shrinks $Z$ towards its mean value. For $z\in\mathbb{R}^{n}$, let
$z_{m}\equiv(1/n)(\sum_{i}z_{i})(1,\dots,1)$ be a vector whose entries are all
the mean value of $z$. For $n\geq 4$, the second type of James-Stein estimator
is then given by

$\hat{\theta}^{\text{mJS}}(z)=z-\frac{(n-3)\Sigma^{-1}}{(z-z_{m})^{T}\Sigma^{-2}(z-z_{m})}(z-z_{m}),$
(9)

with risk

$\displaystyle R($
$\displaystyle\hat{\theta}^{\text{mJS}},\theta)=\text{Tr}(\Sigma)$
$\displaystyle-(n-3)^{2}\mathbb{E}\left[\frac{1}{(Z-Z_{m})^{T}\Sigma^{-2}(Z-Z_{m})}\right].$
(10)

Here, the expectation is evaluated over
$Z\mid\theta\sim\mathcal{N}(\theta,\Sigma)$ [16]. The advantage of this James-
Stein estimator is bigger the more isotropic $\theta$ is (and the “larger”
$\Sigma$ is). Even if there are no correlations between $Z$’s components,
i.e., $\Sigma$ is diagonal, both $\hat{\theta}^{\text{JS}}$ and
$\hat{\theta}^{\text{mJS}}$ combine the components of $Z$ in a non-trivial
manner.

Since the risks of $\hat{\theta}^{\text{JS}}$ and $\hat{\theta}^{\text{mJS}}$
have similar forms, it is usually sufficient to consider
$\hat{\theta}^{\text{JS}}$; $\hat{\theta}^{\text{mJS}}$ will behave similarly.
An exception is given by postselective metrology (section VI), where we show
that only $\hat{\theta}^{\text{mJS}}$ is relevant. Otherwise, one should
choose $\hat{\theta}^{\text{JS}}$ or $\hat{\theta}^{\text{mJS}}$ depending on
the expected properties of $\theta$: if $\theta$ is small, use
$\hat{\theta}^{\text{JS}}$; if $\theta$ is isotropic, use
$\hat{\theta}^{\text{mJS}}$. If no information about $\theta$ is known, one
can use either estimator, or a convex combination of them.

If one receives independent identically distributed normally distributed
copies $Z_{1},\dots Z_{N}$, then the sample mean $\bar{Z}_{N}$ is also
normally distributed: $\bar{Z}_{N}\sim\mathcal{N}(\theta,\Sigma/N)$. Thus one
can define $\hat{\theta}^{\text{MLE}}_{N},\hat{\theta}^{\text{JS}}_{N}$ and
$\hat{\theta}^{\text{mJS}}_{N}$: one simply applies each estimator to the
sample mean. For the James-Stein estimators, this involves modifying $\Sigma$
by a factor of $1/N$ in equations (7) and (9). One can see that

$\text{AD}(\hat{\theta}^{\text{MLE}}_{N},\hat{\theta}^{\text{JS}}_{N},\theta)=1-\Theta(1/N),$
(11)

for both types of James-Stein estimator222We say that $f(x)=\Theta[g(x)]$ if
$f$ and $g$ have the same asymptotic scaling, i.e. there exist constants
$C,D\geq 0$ such that for $x$ sufficiently large $Cg(x)\leq f(x)\leq Dg(x)$..
Thus, the James-Stein advantage diminishes with more data. The MLE and James-
Stein estimators all saturate equation (3): they are asymptotically optimal.

###  II.3 Gaussian states

The most natural application of the James-Stein estimator to quantum metrology
arises when measurement outcomes are normally distributed. Thus, we consider
Gaussian states [17]—states which are characterised by a “mean” and
“variance”, like a Gaussian distribution. More precisely, suppose that the
Hilbert space $\mathcal{H}$ is a tensor product of $\ell$ single-mode Fock
spaces. Examples include the Hilbert space of $\ell/3$ 3-dimensional
particles, or a system with $q$ optical modes [17]. There are $\ell$ pairs of
canonical variables $Q_{i},P_{i},\;i=1,\dots,\ell$ satisfying the canonical
commutation relations

$[Q_{i},Q_{j}]=[P_{i},P_{j}]=0,\quad[Q_{i},P_{j}]=i\delta_{ij}.$ (12)

We group the canonical variables into a single vector
$R=(Q_{1},P_{1},\dots,Q_{\ell},P_{\ell})$. For a state $\rho$ and vector
$\xi\in\mathbb{R}^{2\ell}$ we define $\rho$’s characteristic function [17] by

$\chi_{\rho}(\xi)=\text{Tr}(\rho e^{i\xi^{T}R}).$ (13)

As in classical probability theory, $\rho$ is completely characterised by
$\chi_{\rho}(\xi)$ [18]. A state $\rho$ is called Gaussian [3] if there exists
some vector $r\in\mathbb{R}^{2\ell}$ and some positive semi-definite matrix
$A\geq 0$ such that

$\chi_{\rho}(\xi)=e^{i\xi^{T}r}e^{-\frac{1}{2}(\xi^{T}A\xi)}.$ (14)

This is exactly the characteristic function of a classical Gaussian
probability distribution [14], where $r$ is the mean of the distribution, and
$A$ is its covariance matrix.

##  III Unknown Gaussian Channel

In this Section, we consider the application of the James-Stein estimator to
estimating the mean of a Gaussian state, encoded by some quantum process. We
assume that unknown parameters are encoded by a unitary $U(\theta)$ [19]. In
analogy with the number of classical samples from a probability distribution,
we define $N$ as the number of applications of $U(\theta)$. Thus, $N$
quantifies experimental resources.

Let $\rho_{0}$ be a Gaussian state with mean zero and covariance matrix $A$,
and let $U(\theta)=e^{i(\Omega\theta)^{T}R}$, where

$\Omega=\begin{pmatrix}\omega&&\\\ &\ddots&\\\ &&\omega\end{pmatrix}\text{,
and }\omega=\begin{pmatrix}0&1\\\ -1&0\end{pmatrix}.$ (15)

Then, we find that $\rho(\theta)=U(\theta)\rho_{0}U(\theta)^{\dagger}$ is a
Guassian state with mean $\theta\in\mathbb{R}^{2\ell}$ and covariance matrix
$A$.

As $Q$ and $P$ are incompatible, we cannot simultaneously measure them. In
Ref. [3], it was shown that this problem can be circumvented with the use of
an ancilla Gaussian state $\tilde{\rho}$. The ancilla has the same Hilbert
space as $\rho(\theta)$, with canonical variables
$\tilde{Q}_{i},\tilde{P}_{i},\;i=1,\dots,\ell$. We prepare $\tilde{\rho}$ with
mean 0 and some covariance matrix $\tilde{A}$. Then, we measure the state
$\rho(\theta)\otimes\tilde{\rho}$ with respect to the commuting operators
$Q_{i}+\tilde{P}_{i},P_{i}+\tilde{Q_{i}},\;i=1,\dots,\ell$.

We define $\Sigma=A+\tilde{A}$; note that the Heisenberg uncertainty relation
fundamentally lower bounds how small $\Sigma$ can be (see Ref. [3] for
details). In Appendix A, we show, via a calculation using the Wigner function,
that if $Z$ is a random variable corresponding to the outcome of the
aforementioned measurement process, then
$Z\mid\theta\sim\mathcal{N}(\theta,\Sigma)$. Since $Z$ has a Gaussian
distribution we may use the James-Stein estimator.

If we are given $N$ copies of the unitary $U(\theta)$, a “classical”
estimation strategy is given by preparing $N$ copies of $\rho(\theta)$ and
measuring each separately [20]. One finds that the sample mean has
distribution $\bar{Z}_{N}\sim\mathcal{N}(\theta,\Sigma/N)$ and thus that

$R(\hat{\theta}^{\text{MLE}}_{N},\theta)=\text{Tr}(\Sigma)/N.$ (16)

The James-Stein risk scales as in equation (11).

In Ref. [3], it was shown that upon optimisation of $\tilde{A}$,
$\text{Tr}(\Sigma)=C^{H}(\theta)$. Thus, the aforementioned measurement
combined with the MLE saturates equation (6) and is asymptotically optimal.
However, since
$\text{AD}(\hat{\theta}^{\text{JS}},\hat{\theta}^{\text{MLE}},\theta)>1$ we
note that $R(\hat{\theta}^{\text{JS}}_{N},\theta)<C^{H}(\theta)/N$. Thus, for
any finite $N$, the James-Stein estimator has risk smaller than the asymptotic
minimum, given in equation (6). The James-Stein advantage must vanish as
$N\to\infty$ [as required by equation (6)] but can be substantial for finite
$N$; see Fig. 1 (orange curve).

A common “quantum” strategy [21, 22, 19, 20, 23] is to apply $U(\theta)$
sequentially $N$ times to produce the state
$\rho_{N}(\theta)=U^{N}(\theta)\rho_{0}{U^{N}(\theta)^{\dagger}}$. Now,
$\rho_{N}(\theta)$ is a Gaussian state with covariance matrix $A$ and mean
$N\theta$. Thus, letting $Z^{(N)}$ denote the outcome of measuring the mean of
$\rho_{N}(\theta)$, we have
$Z^{(N)}/N\mid\theta\sim\mathcal{N}(\theta,\Sigma/N^{2})$.

We denote by $\hat{\theta}^{\text{qMLE}}_{N}$ and
$\hat{\theta}^{\text{qJS}}_{N}$ the MLE and James-Stein estimators applied to
$Z^{(N)}/N$, respectively. We find that

$\text{AD}(\hat{\theta}^{\text{qMLE}}_{N},\hat{\theta}^{\text{MLE}}_{N},\theta)=N.$
(17)

Using the unitary sequentially gives a factor of $N$ improvement over using it
on separate copies of $\rho_{0}$. Similarly,
$\text{AD}(\hat{\theta}^{\text{qJS}}_{N},\hat{\theta}^{\text{JS}}_{N},\theta)=\Theta(N)$,
so that the James-Stein estimator is also improved. However, there is a
corresponding decrease in the James-Stein advantage compared to equation (11):

$\text{AD}(\hat{\theta}^{\text{qMLE}}_{N},\hat{\theta}^{\text{qJS}}_{N},\theta)=1-\Theta(1/N^{2}).$
(18)

In summary, whilst the James-Stein estimator still has an advantage over the
MLE, this advantage is suppressed by coherent quantum effects, in the form of
sequential unitary applications. See the dotted curve in Fig. 1. The use of an
entangled probe state will show similar behaviours.

##  IV Gaussian Channels with Noise

In Section III, we assumed perfect (noiseless) unitary evolution. In this
Section, we lift this assumption and consider a common noise model. We show
that the James-Stein and quantum-metrological advantages can be simultaneously
large.

Suppose that there is coherent noise, such that we apply $U(\phi)$ instead of
$U(\theta)$, where $\phi$ fluctuates randomly around $\theta$. This is
modelled by the channel

$\Lambda_{\theta}(\sigma)=\int d\phi f_{\theta}(\phi)U(\phi)\sigma
U(\phi)^{\dagger},$ (19)

where $f_{\theta}$ is the pdf of a normal distribution with mean $\theta$ and
covariance matrix $\Delta$. One can check that $\Lambda_{\theta}$ maps
Gaussian states to Gaussian states: if $\rho$ is a Gaussian state with mean
$r$ and covariance $A$, then $\Lambda_{\theta}(\rho)$ is a Gaussian state with
mean $r+\theta$ and covariance $\Delta+A$. Channels that map Gaussian states
to Gaussian states are also called Gaussian [24].

We let $\bar{W}_{N}$ denote the sample mean of measuring $N$ copies of
$\Lambda_{\theta}(\rho_{0})$, and $W^{(N)}$ correspond to measuring a single
copy of $\Lambda_{\theta}^{N}(\rho_{0})$. We find that

$\displaystyle\bar{W}_{N}\mid\theta$
$\displaystyle\sim\mathcal{N}(\theta,(\Sigma+\Delta)/N)$ (20) $\displaystyle
W^{(N)}/N\mid\theta$
$\displaystyle\sim\mathcal{N}(\theta,\Sigma/N^{2}+\Delta/N).$ (21)

We label the estimation strategies $\hat{\theta}^{\text{MLE}}$,
$\hat{\theta}^{\text{JS}}$, $\hat{\theta}^{\text{qMLE}}$ and
$\hat{\theta}^{\text{qJS}}$ applied to the (noisy) $W$ variables with a “w” in
superscript. We find that the unbounded asymptotic quantum-metrological
advantage of equation (17) is reduced to a constant factor:

$\text{AD}(\hat{\theta}^{\text{wqMLE}}_{N},\hat{\theta}^{\text{wMLE}}_{N},\theta)=1+\frac{\text{Tr}(\Sigma)}{\text{Tr}(\Delta)}-\Theta(1/N).$
(22)

Nevertheless, we stress that it can still be advantageous to sequentially
apply $\Lambda_{\theta}$. Depending on the relative size of $\Sigma$ and
$\Delta$, this advantage may be considerable. Moreover, we find that

$\text{AD}(\hat{\theta}^{\text{wqMLE}}_{N},\hat{\theta}^{\text{wqJS}}_{N},\theta)=1-\Theta(1/N).$
(23)

Comparing equations (18) and (23), we see that the James-Stein advantage is
larger in the presence of noise. This is demonstrated in Fig. 1 for a fixed
$\theta\in\mathbb{R}^{4}$. The blue curve
$\text{AD}(\hat{\theta}^{\text{wqJS}}_{N},\hat{\theta}^{\text{wJS}}_{N},\theta)$
shows that sequentially applying $U(\theta)$ leads to a significant advantage,
even in the presence of noise. When $N=50$, we find that
$\text{AD}(\hat{\theta}^{\text{wqJS}}_{N},\hat{\theta}^{\text{wJS}}_{N},\theta)\approx
1.8$; a quantum metrological improvement of roughly $80\%$. From the other
curves, we see that the James-Stein advantage falls off less rapidly when
$U(\theta)$ is applied sequentially with noise than without noise. When
$N=25$, the James-Stein advantage is about 10% when sequentially applying
$U(\theta)$ with noise; it is 1% when sequentially applying $U(\theta)$
without noise. We emphasise that these James-Stein advantages are “free”. No
extra experimental resources are required, simply an alternative post-
processing of the measured data.

Figure 1: Comparison of various estimation strategies when
$\theta=(0.5,-0.2,0.3,0.1)\in\mathbb{R}^{4}$, and
$\Sigma,\Delta=4\times\mathbbm{1}$. The James-Stein advantage with sequential
unitary applications in the presence of noise (green curve) may be large. This
contrasts to the case of sequential applications with no noise (dotted red
line), where the James-Stein advantage quickly becomes negligible. The James-
Stein advantage is even larger when applying $\Lambda_{\theta}$ on separate
copies of $\rho(\theta)$ (orange curve), but in this case one loses the
quantum advantage from sequential applications (demonstrated by the blue
curve).

##  V Bayesian Gaussians

In Bayesian parameter estimation, one assumes that there is some prior
distribution $\pi$ on the parameter space $\Theta$ that encodes an a priori
belief about the parameter to be estimated [4]. With respect to this belief,
one can consider the Bayes risk of an estimator $\hat{\theta}^{\text{}}$:

$R_{\pi}(\hat{\theta}^{\text{}})=\mathbb{E}_{\vartheta\sim\pi}[R(\hat{\theta}^{\text{}},\vartheta)]=\int
R(\hat{\theta}^{\text{}},\theta)d\pi(\theta).$ (24)

An estimator $\hat{\theta}^{\text{}}$ is called Bayesian if it minimises the
Bayes risk, as defined in equation (24).

We consider the case of a normally distributed prior
$\theta\sim\mathcal{N}(\theta_{0},\Xi)$. To capture the different cases of
Sections III and IV, we consider a random variable

$Z_{N}\mid\theta\sim\mathcal{N}(\theta,\Gamma_{N}),$ (25)

where $\Gamma_{N}$ is some full-rank covariance matrix that depends on $N$. We
denote the Bayesian, MLE and James-Stein estimators by
$\hat{\theta}^{\text{bMLE}}$ and $\hat{\theta}^{\text{bJS}}$, respectively. In
Appendix B.1, we show that the Bayes estimator for equation (25) is

$\hat{\theta}^{\text{B}}(z)=(\Gamma_{N}^{-1}+\Xi^{-1})^{-1}(\Gamma_{N}^{-1}z+\Xi^{-1}\theta_{0}).$
(26)

We also calculate the Bayes risk of the MLE and James-Stein estimators (see
appendix B.2); the results are summarised in Table 1.

$\hat{\theta}^{\text{}}$ |  $R_{\pi}(\hat{\theta}^{\text{}})$
---|---
$\hat{\theta}^{\text{bMLE}}$ |  $\text{Tr}(\Gamma_{N})$
$\hat{\theta}^{\text{bJS}}$ |  $\text{Tr}(\Gamma_{N})-(2q-2)^{2}\mathbb{E}\left[\frac{1}{{Y_{N}}^{T}\Gamma_{N}^{-2}{Y_{N}}}\right]$
$\hat{\theta}^{\text{B}}$ |  $\text{Tr}(\Gamma_{N})-\text{Tr}\left[\Gamma_{N}^{2}(\Gamma_{N}+\Xi)^{-1}\right]$

Table 1:  Comparison of the Bayes risk of the maximum-likelihood, James-Stein
and Bayes estimators. Here, $Y_{N}\sim\mathcal{N}(\theta_{0},\Gamma_{N}+\Xi)$.

By the definition of the Bayes estimator, and since the James-Stein estimator
performs better everywhere than the MLE [see equation (11)], we have that

$R_{\pi}(\hat{\theta}^{\text{B}})<R_{\pi}(\hat{\theta}^{\text{bJS}})<R_{\pi}(\hat{\theta}^{\text{bMLE}}).$
(27)

Thus, the Bayes estimator $\hat{\theta}^{\text{B}}$ is preferable. However, in
practice, one may be unable to use $\hat{\theta}^{\text{B}}$, since it depends
explicitly on the covariance matrix $\Xi$, which may be unknown at the point
of experiment. Additionally, there is no hope to estimate $\Xi$; $\theta$ is
sampled according to $\pi$, and subsequently $\theta$ is encoded into $N$
copies of $U(\theta)$. Thus, the $N$ copies of $U(\theta)$ yield no
information about $\Xi$. Consequently, the James-Stein estimator may be the
best available alternative to the Bayes estimator.

If $\Gamma_{N}\to 0$ as $N\to\infty$, then the Bayesian advantage over the MLE
is second order in $\Gamma_{N}$:
$R_{\pi}(\hat{\theta}^{\text{B}})=R_{\pi}(\hat{\theta}^{\text{bMLE}})-O(\Gamma_{N}^{2})$
. Thus,
$R_{\pi}(\hat{\theta}^{\text{B}})/R_{\pi}(\hat{\theta}^{\text{bMLE}})=1-O(\Gamma_{N})\to
1$. Moreover, the faster, $\Gamma_{N}$ converges to 0, the faster
$R_{\pi}(\hat{\theta}^{\text{B}})/R_{\pi}(\hat{\theta}^{\text{bMLE}})$ will
converge to 1. Therefore, sequential applications of $U(\theta)$ will move the
MLE risk, and hence by equation (27), the James-Stein risk, closer to the
optimal Bayes estimator’s risk. That is, the use of quantum-metrological
techniques can move the James-Stein risk closer to the unattainable minimal
Bayes risk.

##  VI Postselected Guassian states

In this Section, we consider the application of postselection (defined below)
[25, 26, 27, 6] to the estimation of the mean of a Gaussian state. We also
investigate how postselection affects the James-Stein advantage. For
notational simplicity, we assume that the covariance matrix
$\Sigma\propto\mathbbm{1}$ and that our state is pure. Additionally, by the
ancilla trick described in Section III, we may assume that $\theta$ is encoded
solely in the positional degree of freedom of our state. That is,

$\ket{\psi_{\theta}}=C^{n/2}\int d^{n}xe^{-{\norm{x-\theta}^{2}/B}}\ket{x},$
(28)

for some known $n,B>0$, some unknown $\theta\in\mathbb{R}^{n}$, and
$C=\sqrt{2/\pi B}$. For example, this is the state of $n\geq 4$ (continuous-
variable) non-entangled, quantum sensors, each estimating some parameter
$\theta_{i}\in\mathbb{R}$. Given $N$ copies of $\ket{\psi_{\theta}}$, as
detailed in Section III, one could measure the position of each state and take
the sample mean to estimate $\theta$ with the corresponding risk scaling as
$1/N$.

Many experiments suffer from limitations such as detector saturation [28, 29,
30, 31] (if the intensity of photons arriving at a photon detector is too
high, the detector performance degrades) and dead-times [32, 33, 34] (a
detector must reset for a brief time between particle observations).
Sometimes, one can only measure $M$ of the $N$ states $\ket{\psi_{\theta}}$
produced per unit time. It may be that $M\ll N$, so that the rate at which one
can measure information is much less than the rate that one can produce it.

If one were to simply discard a random fraction $(1-M/N)$ of the copies of
$\ket{\psi_{\theta}}$, the risk of our estimation strategies would scale as
$1/M$—worse than $1/N$. However, if one knows that $\theta$ is a small
displacement from some known value $\theta_{0}$, i.e.
$\delta=\theta-\theta_{0}$ satisfies $\norm{\delta}^{2}\ll M/N$, then the
information carried by $N$ probe states can be losslessly compressed down to
$M\ll N$ probe states [25]. This is similar in spirit to weak-value
amplification techniques [35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47].

Instead of discarding states at random, one can utilise a more clever filter
described by a 2-outcome POVM, with operators $F$ and $\mathbbm{1}-F$. The
outcome corresponding to $F$ (alt. $\mathbbm{1}-F$) has the probe pass (alt.
not pass) the filter. Let $t\in[0,1]$ be a transmission parameter and take

$F=\mathbbm{1}-(1-t^{2})\outerproduct{\psi_{\theta_{0}}}{\psi_{\theta_{0}}}.$
(29)

$F$ transmits all states perpendicular to $\ket{\psi_{\theta_{0}}}$, but only
allows $\ket{\psi_{\theta_{0}}}$ through with probability $t^{2}$. As in [25],
we implement $F$ with the natural Kraus operator

$K=\mathbbm{1}-(1-t)\outerproduct{\psi_{\theta_{0}}}{\psi_{\theta_{0}}}.$ (30)

If a state passes the filter, it will be in the state
$K\ket{\psi_{\theta}}/\norm{K\ket{\psi_{\theta}}}$. In Appendix C.1, we show
that if $t\gg\norm{\delta}$, then the post-filter state of the particle is

$\frac{K\ket{\psi_{\theta}}}{\norm{K\ket{\psi_{\theta}}}}\approx C^{n/2}\int
d^{n}x\;e^{-\norm{x-\theta_{0}-\delta/t}^{2}/B}\ket{x}$ (31)

Moreover, we show the probability of passing the filter is approximately
$t^{2}$.

If we measure the probe’s position after postselection, we observe a random
variable $X$, approximately normally distributed as

$X\mid\theta\sim\mathcal{N}(\theta_{0}+\delta/t,[B/4]\mathbbm{1}).$ (32)

Due to the factor of $1/t$, $X$ is more sensitive in changes to $\theta$ than
position measurements of $\ket{\psi_{\theta}}$. Thus
$Y=t(X-\theta_{0})+\theta_{0}$ is distributed as

$Y\mid\theta\sim\mathcal{N}(\theta,t^{2}[B/4]\mathbbm{1}).$ (33)

It takes (on average) $1/t^{2}$ copies of $\ket{\psi_{\theta}}$ until a state
passes the filter, but the resulting state has a variance that is $t^{2}$
smaller than before. Thus, the information in $N$ copies of
$\ket{\psi_{\theta}}$ is losslessly compressed into $M\approx t^{2}N$ copies
of $K\ket{\psi_{\theta}}/\norm{K\ket{\psi_{\theta}}}$. Such lossless
compression has been shown to require genuine quantum effects [6, 25, 46].

We proceed to investigate the advantage of the James-Stein estimator in
experiments with post-selective filters. We model the aforementioned situation
in which one is limited by a detector. If we can only measure $M$ states per
unit time, we send $N=M/t^{2}$ copies of $\ket{\psi_{\theta}}$ through the
filter. In this way, we ensure that (on average) $M$ states arrive at the
detector per unit time. Thus, we can make $t$ arbitrarily small (in the
approximation $M\ll N$).

We start with with no prior knowledge about $\theta$ and hence no initial
guess $\theta_{0}$. This means that $\hat{\theta}^{\text{JS}}$ is a sub-
optimal choice of James-Stein estimator, since we do not expect $\theta$ to be
close to any a priori known value. However, if the probes are measuring
similar quantities, one could expect that the $\theta_{i}$’s would be similar
and thus that $\hat{\theta}^{\text{mJS}}$ would still perform well. Thus, we
only consider $\hat{\theta}^{\text{mJS}}$ in this section, and not
$\hat{\theta}^{\text{JS}}$. Because we require that $t\gg\norm{\delta}$, one
cannot decrease $t$ until one has a good estimate $\theta_{0}$ of $\theta$. We
run an iterative strategy (see Appendix C.2 for the algorithm and Appendix C.3
for details on the numerics), where one estimates $\delta$ (using sample
variance) and then decreases $t$ when one is confident that $\delta$ is
sufficiently small. We consider using the MLE or James-Stein estimator to
generate our estimates $\theta_{0}$ of $\theta$, and denote these strategies
by $\hat{\theta}^{\text{pMLE}}$ and $\hat{\theta}^{\text{pmJS}}$,
respectively. As a baseline, we compare this to measuring $N$ copies of
$\ket{\psi_{\theta}}$ with no postselection ($t=1$), and using either the MLE
$\hat{\theta}^{\text{MLE}}$ or James-Stein estimator
$\hat{\theta}^{\text{mJS}}$.

In Fig. 2(a), we fix $\theta\in\mathbb{R}^{4}$. The blue curve shows that,
when using the MLE, postselection gives an advantage. After some initial few
measurements, this advantage increases linearly with the number of
measurements. This is due to the increased sensitivity by $1/t$ in the
postselected states, which leads to a more accurate estimate of $\theta$. The
orange curve shows that the postselected James-Stein estimator also
outperforms the non-postselected MLE. The green (alt. red) curve shows the
James-Stein advantage without (alt. with) postselection. We see that using
postselection the James-Stein advantage is initially increased, but eventually
decays to below the non-postselected James-Stein advantage. This is because
using $\hat{\theta}^{\text{mJS}}$ allows one to postselect more quickly, so
that the initial advantage increases. However, for larger values of $N$, $t$
shrinks, so the sample variance (and hence the James-Stein advantage)
decreases. This behaviour is further analysed below. For $N<40$, postselection
boosts the James-Stein advantage.

As discussed above, we expect $\hat{\theta}^{\text{mJS}}$ to perform better
when $\theta$ is more “isotropic”. Letting $\bar{\theta}=\sum\theta_{i}/n$, we
quantify $\theta$’s isotropy by

$v(\theta)=\sum(\theta_{i}-\bar{\theta})^{2}$ (34)

To quantify the effect of postselection on the James-Stein advantage, we
consider the ratio

$\text{PAD}(N,\theta)=\frac{\text{AD}(\hat{\theta}^{\text{pmJS}}_{N},\hat{\theta}^{\text{pMLE}}_{N},\theta)}{\text{AD}(\hat{\theta}^{\text{mJS}}_{N},\hat{\theta}^{\text{MLE}}_{N},\theta)}.$
(35)

If $\text{PAD}>1$, then postselection has increased the James-Stein advantage
(compared to the non-postselected case). In Fig. 2(b), we plot
$\text{PAD}(N,\theta)$ against $N$ for 4 different values of
$\theta\in\mathbb{R}^{4}$, each with a different $v(\theta)$. As described
above, the James-Stein estimator allows stronger postselection for smaller
$N$, so that $\text{PAD}>1$. However, as $N$ increases, the $1/t$ factor in
sensitivity increases $v(\theta)$ by a factor of $1/t^{2}$, which degrades the
James-Stein advantage. Thus, PAD eventually decreases. These effects are
magnified when $v(\theta)$ is smaller.

##  VII Conclusion

We have introduced the James-Stein estimator to quantum metrology. We have
focused on Gaussian shift estimation, but our techniques generalise to most
common parameterised quantum states. We have shown that across a wide array of
quantum-metrology protocols, the James-Stein estimator outperforms the MLE.
This means that by simply changing post-processing techniques, one can
increase experimental performance, without requiring any additional
experimental resources. Our results highlight the non-trivial relationship
between the asymptotic and non-asymptotic regimes of (quantum) metrology [48,
49]; strategies that are asymptotically optimal may be sub-optimal with finite
resources. We conclude that the James-Stein estimator can be a useful tool in
quantum metrology. Moreover, quantum phenomena have a non-trivial relationship
with the James-Stein advantage; they may increase it (sec. VI), diminish it
(sec. III), or neither (sec. IV).

## Acknowledgements

The authors wish to thank P. Karlsson for prompting them to investigate the
James-Stein estimator in a quantum setting. Further, the authors thank N.
Mertig, F. Venn, C. Long and J. Smith for their useful discussions and
comments. W.S. was supported by the EPSRC and Hitachi. S.S. acknowledges
support from the Royal Society University Research Fellowship. D.R.M.A.-S. was
supported by Girton College.

(a)

(b)

Figure 2: (a) Comparison of various postselected and non-postselected
estimation strategies when $\theta=(1,-2,3,1.5)/30$, $B=1$. (b) Comparison of
the the affect of postselection on the James-Stein advantage for various
$\theta$, $B=1$. The small kink in the blue line is due to computational shot
noise.

## References

  * Stein [1956] C. Stein, Proceedings of the Third Berkeley Symposium on Mathematical Statistics and Probability 1, 197 (1956).
  * Stein and James [1961] C. Stein and W. James, Proceedings of the Fourth Berkeley Symposium on Mathematical Statistics and Probability 1, 361 (1961).
  * Demkowicz-Dobrzański  _et al._ [2020] R. Demkowicz-Dobrzański, W. Górecki, and M. Guţă, Journal of Physics A: Mathematical and Theoretical 53, 363001 (2020).
  * Wasserman [2004] L. Wasserman,  _All of Statistics_ (Springer, New York, NY, 2004) p. 202.
  * Demkowicz-Dobrzański  _et al._ [2012] R. Demkowicz-Dobrzański, J. Kołodyński, and M. Guţă, Nature communications 3, 1063 (2012).
  * Arvidsson-Shukur  _et al._ [2020] D. R. M. Arvidsson-Shukur, N. Yunger Halpern, H. V. Lepage, A. A. Lasek, C. H. W. Barnes, and S. Lloyd, Nat. Commun. 11, 3775 (2020).
  * Lehmann and Casella [1998] E. L. Lehmann and G. Casella,  _Theory of Point Estimation_ (Springer, New York, NY, 1998).
  * Braunstein and Caves [1994] S. L. Braunstein and C. M. Caves, Phys. Rev. Lett. 72, 3439 (1994).
  * Albarelli  _et al._ [2020] F. Albarelli, M. Barbieri, M. Genoni, and I. Gianani, Physics Letters A 384, 126311 (2020).
  * Vaart [1998] A. W. v. d. Vaart, Efficiency of estimators, in  _Asymptotic Statistics_, Cambridge Series in Statistical and Probabilistic Mathematics (Cambridge University Press, 1998) p. 108–124.
  * Helstrom [1969] C. W. Helstrom, Journal of Statistical Physics 1, 231 (1969).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, 2010).
  * Yang  _et al._ [2019] Y. Yang, G. Chiribella, and M. Hayashi, Communications in Mathematical Physics 368, 223 (2019).
  * van der Vaart [2000] A. W. van der Vaart,  _Asymptotic Statistics_ (Cambridge University Press, 2000).
  * Chételat and Wells [2012] D. Chételat and M. T. Wells, The Annals of Statistics 40, 3137 (2012).
  * Li and Bhoj [1988] T. F. Li and D. S. Bhoj, Scandinavian Journal of Statistics 15, 33 (1988).
  * Weedbrook  _et al._ [2012] C. Weedbrook, S. Pirandola, R. García-Patrón, N. J. Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys. 84, 621 (2012).
  * Englert and Wódkiewicz [2003] B.-G. Englert and K. Wódkiewicz, International Journal of Quantum Information 1, 153 (2003).
  * Giovannetti  _et al._ [2006] V. Giovannetti, S. Lloyd, and L. Maccone, Phys. Rev. Lett. 96, 010401 (2006).
  * Giovannetti  _et al._ [2011] V. Giovannetti, S. Lloyd, and L. Maccone, Nature Photonics 5, 222 (2011).
  * Smith  _et al._ [2022a] J. G. Smith, C. H. W. Barnes, and D. R. M. Arvidsson-Shukur, Phys. Rev. A 106, 062615 (2022a).
  * Smith  _et al._ [2022b] J. G. Smith, C. H. W. Barnes, and D. R. M. Arvidsson-Shukur, Phys. Rev. A 106, 062615 (2022b).
  * Ji  _et al._ [2008] Z. Ji, G. Wang, R. Duan, Y. Feng, and M. Ying, IEEE Transactions on Information Theory 54, 5172 (2008).
  * Eisert and Wolf [2007] J. Eisert and M. M. Wolf, Quantum Information with Continous Variables of Atoms and Light 10.48550/ARXIV.QUANT-PH/0505151 (2007).
  * Jenne and Arvidsson-Shukur [2022] J. H. Jenne and D. R. M. Arvidsson-Shukur, Unbounded and lossless compression of multiparameter quantum information (2022).
  * Salvati  _et al._ [2023] F. Salvati, W. Salmon, C. H. Barnes, and D. R. Arvidsson-Shukur, arXiv preprint arXiv:2307.08648 (2023).
  * Lupu-Gladstein  _et al._ [2022a] N. Lupu-Gladstein, Y. B. Yilmaz, D. R. Arvidsson-Shukur, A. Brodutch, A. O. Pang, A. M. Steinberg, and N. Y. Halpern, Physical Review Letters 128, 220504 (2022a).
  * Luu and Jiang [2006] J. X. Luu and L. A. Jiang, Applied optics 45, 3798 (2006).
  * Wang  _et al._ [2010] J. Wang, T. Liu, S. Jiao, R. Chen, Q. Zhou, K. K. Shung, L. V. Wang, and H. F. Zhang, Journal of biomedical optics 15, 021317 (2010).
  * Frehlich and Ochs [1990] R. G. Frehlich and G. R. Ochs, Applied optics 29, 548 (1990).
  * Lupu-Gladstein  _et al._ [2022b] N. Lupu-Gladstein, Y. B. Yilmaz, D. R. M. Arvidsson-Shukur, A. Brodutch, A. O. T. Pang, A. M. Steinberg, and N. Y. Halpern, Phys. Rev. Lett. 128, 220504 (2022b).
  * Hachisu and Ido [2015] H. Hachisu and T. Ido, Japanese Journal of Applied Physics 54, 112401 (2015).
  * Wahl  _et al._ [2020] M. Wahl, T. Röhlicke, S. Kulisch, S. Rohilla, B. Krämer, and A. C. Hocke, Review of Scientific Instruments 91 (2020).
  * Rohde and Ralph [2006] P. P. Rohde and T. C. Ralph, Journal of Modern Optics 53, 1589 (2006).
  * Hosten and Kwiat [2008] O. Hosten and P. Kwiat, Science 319, 787 (2008).
  * Dixon  _et al._ [2009] P. B. Dixon, D. J. Starling, A. N. Jordan, and J. C. Howell, Phys. Rev. Lett. 102, 173601 (2009).
  * Turner  _et al._ [2011] M. D. Turner, C. A. Hagedorn, S. Schlamminger, and J. H. Gundlach, Opt. Lett. 36, 1479 (2011).
  * Pfeifer and Fischer [2011] M. Pfeifer and P. Fischer, Opt. Express 19, 16508 (2011).
  * Starling  _et al._ [2010a] D. J. Starling, P. B. Dixon, A. N. Jordan, and J. C. Howell, Phys. Rev. A 82, 063822 (2010a).
  * Starling  _et al._ [2010b] D. J. Starling, P. B. Dixon, N. S. Williams, A. N. Jordan, and J. C. Howell, Phys. Rev. A 82, 011802 (2010b).
  * Xu  _et al._ [2013] X.-Y. Xu, Y. Kedem, K. Sun, L. Vaidman, C.-F. Li, and G.-C. Guo, Phys. Rev. Lett. 111, 033604 (2013).
  * Magaña-Loaiza  _et al._ [2014] O. S. Magaña-Loaiza, M. Mirhosseini, B. Rodenburg, and R. W. Boyd, Phys. Rev. Lett. 112, 200401 (2014).
  * Strübi and Bruder [2013] G. Strübi and C. Bruder, Phys. Rev. Lett. 110, 083605 (2013).
  * Viza  _et al._ [2013] G. I. Viza, J. Martínez-Rincón, G. A. Howland, H. Frostig, I. Shomroni, B. Dayan, and J. C. Howell, Opt. Lett. 38, 2949 (2013).
  * Egan and Stone [2012] P. Egan and J. A. Stone, Opt. Lett. 37, 4991 (2012).
  * Pusey [2014] M. F. Pusey, Phys. Rev. Lett. 113, 200401 (2014).
  * Kunjwal  _et al._ [2019] R. Kunjwal, M. Lostaglio, and M. F. Pusey, Physical Review A 100, 042116 (2019).
  * Salmon  _et al._ [2023] W. Salmon, S. Strelchuk, and D. Arvidsson-Shukur, Quantum 7, 998 (2023).
  * Meyer  _et al._ [2023] J. J. Meyer, S. Khatri, D. S. França, J. Eisert, and P. Faist, arXiv preprint arXiv:2307.06370 (2023).
  * Tatarskiĭ [1983] V. I. Tatarskiĭ, Soviet Physics Uspekhi 26, 311 (1983).
  * Banerjee  _et al._ [2005] A. Banerjee, X. Guo, and H. Wang, IEEE Transactions on Information Theory 51, 2664 (2005).
  * Casella  _et al._ [2004] G. Casella, C. P. Robert, and M. T. Wells, Lecture Notes-Monograph Series , 342 (2004).

##  Appendix A Distribution of Gaussian measurement outcome

In this appendix, we find the distribution of the measurement described in
Section III (and Ref. [3]). We begin by recapping the definition of a Gaussian
state. Let $\mathcal{H}=\mathcal{F}^{\otimes\ell}$ be the tensor product of
$\ell$ single mode Fock spaces. Let
$R=(X_{1},P_{1},\dots,X_{\ell},P_{\ell})^{T}$ be a vector of operators, where
$(Q_{i},P_{i})$ are the position and momentum canonical variables of the $i$th
system, so that

$[Q_{i},P_{j}]=i\delta_{ij}.$ (36)

The characteristic function of a state $\rho$ is the function

$\chi_{\rho}:\mathbb{R}^{2\ell}\to\mathbb{C},\;\xi\mapsto\text{Tr}(\rho
e^{i\xi^{T}R}).$ (37)

A state is called Gaussian if for some $r\in\mathbb{R}^{2\ell}$ and $A$ some
positive semi-definite matrix, we have that

$\chi_{\rho}(\xi)=e^{i\xi^{T}r}e^{-\frac{1}{2}(\xi^{T}A\xi)}.$ (38)

The Wigner function [50] of the state is defined as the Fourier transformation
of the characteristic function, that is

$W_{\rho}:\mathbb{R}^{2\ell}\to\mathbb{C},\;x\mapsto\int\frac{d^{2\ell}\xi}{(2\pi)^{2\ell}}e^{-ix^{T}\xi}\chi_{\rho}(\xi).$
(39)

If $\rho$ is a Gaussian state, then we find that

$W_{\rho}(x)\propto\exp(-\frac{1}{2}(x-r)^{T}A^{-1}(x-r)),$ (40)

i.e. a Gaussian in phase space.

As in Section III, suppose we have a Gaussian state with a known covariance
matrix $A$, but unknown mean $\theta$ that we wish to estimate. That is

$W_{\rho}(x)\propto\exp(-\frac{1}{2}(x-\theta)^{T}A^{-1}(x-\theta)).$ (41)

We cannot simultaneously measure the position and the momentum of the state,
but detailed in Ref. [3], we can circumvent this problem with the use of an
ancilla. The ancilla is in a Gaussian state $\widetilde{\rho}$ which has the
roles of position and momentum inverted. That is to say we define a vector of
operators

$\widetilde{R}=(\widetilde{P}_{1},\widetilde{Q}_{1},\dots,\widetilde{P}_{\ell},\widetilde{Q}_{\ell})^{T},$
(42)

and ask that $\widetilde{\rho}$ satisfies

$\text{Tr}(\widetilde{\rho}e^{i\xi^{T}\widetilde{R}})=e^{i\xi^{T}\widetilde{r}}e^{-\frac{1}{2}\xi^{T}\widetilde{A}^{-1}\xi},$
(43)

for some vector $\widetilde{r}$ and covariance matrix $\widetilde{A}$.

Defining the swap matrix

$S=\begin{pmatrix}s&&\\\ &\ddots&\\\ &&s\end{pmatrix},\quad
s=\begin{pmatrix}0&1\\\ 1&0\end{pmatrix},$ (44)

we find that $\widetilde{\rho}$ is a Gaussian state with mean vector
$S\widetilde{r}$ and covariance matrix $S\widetilde{A}S$.

The joint state $\rho\otimes\widetilde{\rho}$ has Wigner function

$W_{\rho\otimes\widetilde{\rho}}(x)\propto\exp(-\frac{1}{2}(\bar{x}-\bar{\theta})^{T}\bar{A}^{-1}(\bar{x}-\bar{\theta})),$
(45)

where

$\bar{x}=\begin{pmatrix}\mathbbm{1}&\\\
&S\end{pmatrix}x,\quad\bar{\theta}=\begin{pmatrix}\theta\\\
\widetilde{r}\end{pmatrix},\quad\bar{A}=\begin{pmatrix}A&\\\
&\widetilde{A}\end{pmatrix}.$ (46)

Ref. [3] notes that the operators

$\\{Y_{i}=R_{i}+\widetilde{R}_{i}\\}.$ (47)

mutually commute, and thus can be simultaneously measured. In fact, they show
that this measurement scheme is optimal—the Fisher information of the
resulting distribution is maximal.

Consider the probability density $f$ at the measurement outcome
$Q_{i}+\widetilde{P}_{i}=\lambda_{i}$, $P_{i}+\widetilde{Q}_{i}=\mu_{i}$,
denoted $f(\lambda_{1},\mu_{1},\dots,\lambda_{\ell},\mu_{\ell})$. Using
standard properties of the Wigner function (see Ref. [50] for details), we
find

$\displaystyle f(\lambda_{1},\mu_{1},\dots,\lambda_{\ell},\mu_{\ell})$
$\displaystyle=\int_{\mathbb{R}^{\ell}}d^{\ell}x\int_{\mathbb{R}^{\ell}}d^{\ell}\widetilde{x}\;W_{\rho\otimes\widetilde{\rho}}(x_{1},\dots
x_{\ell},\widetilde{x}_{1},\dots\widetilde{x}_{\ell},\mu_{1}-\widetilde{x}_{1},\dots,\mu_{\ell}-\widetilde{x}_{\ell},\lambda_{1}-x_{1},\dots,\lambda_{\ell}-x_{\ell})$
(48)
$\displaystyle\propto\int_{\mathbb{R}^{\ell}}d^{\ell}x\int_{\mathbb{R}^{\ell}}d^{\ell}\widetilde{x}\;\exp(-\frac{1}{2}(\bar{x}-\bar{\theta})^{T}\bar{A}^{-1}(\bar{x}-\bar{\theta})),$
(49)

where

$\bar{x}=(x_{1},\mu_{1}-\widetilde{x}_{1},\dots,x_{\ell},\mu_{\ell}-\widetilde{x}_{\ell},\lambda_{1}-x_{1},\widetilde{x}_{1},\dots,\lambda_{\ell}-x_{\ell},\widetilde{x}_{\ell})^{T}.$
(50)

We wish to extract the $\lambda$ and $\mu$ dependence from this integral. Let
$x_{0}=(x_{1},-\widetilde{x}_{1},\dots,x_{\ell},-\widetilde{x}_{\ell})^{T}$,
$\mu_{0}=(0,\mu_{1},\dots,0,\mu_{\ell})^{T}-\theta$ and
$\lambda_{0}=(\lambda_{1},0,\dots,\lambda_{\ell},0)^{T}-\widetilde{r}$. Using
the definition of $\bar{A}$ and $\bar{\theta}$, we find that

$f(\lambda_{1},\mu_{1},\dots,\lambda_{\ell},\mu_{\ell})\propto\int_{\mathbb{R}^{\ell}}d^{\ell}x\int_{\mathbb{R}^{\ell}}d^{\ell}\widetilde{x}\;\exp(-\frac{1}{2}\left((x_{0}+\mu_{0})^{T}A^{-1}(x_{0}+\mu_{0})+(\lambda_{0}-x_{0})^{T}\widetilde{A}^{-1}(\lambda_{0}-x_{0})\right)).$
(51)

We complete the square in the exponent, i.e. write it in the form

$-\frac{1}{2}\left((x_{0}+r_{0})^{T}M(x_{0}+r_{0})+R\right),$ (52)

for some $r_{0},R$. From inspection, we can see that

$M=A^{-1}+\widetilde{A}^{-1}\quad\text{and}\quad
r_{0}=M^{-1}\left(A^{-1}\mu_{0}-\widetilde{A}^{-1}\lambda_{0}\right).$ (53)

This gives

$\displaystyle R$
$\displaystyle=\mu_{0}^{T}A^{-1}\mu_{0}+\lambda_{0}^{T}\widetilde{A}^{-1}\lambda_{0}-r_{0}^{T}Mr_{0},$
(54)
$\displaystyle=\mu_{0}^{T}M_{1}\mu_{0}+\lambda_{0}^{T}M_{2}\lambda_{0}+\lambda_{0}^{T}M_{3}\mu_{0}+\mu_{0}^{T}M_{4}\lambda_{0},$
(55)

where the matrices $M_{i}$ are given by

$\displaystyle M_{1}$
$\displaystyle=A^{-1}-A^{-1}(A^{-1}+\widetilde{A}^{-1})^{-1}A^{-1},\quad$
$\displaystyle
M_{2}=\widetilde{A}^{-1}-\widetilde{A}^{-1}(A^{-1}+\widetilde{A}^{-1})^{-1}\widetilde{A}^{-1},$
$\displaystyle M_{3}$
$\displaystyle=\widetilde{A}^{-1}(A^{-1}+\widetilde{A}^{-1})^{-1}A^{-1},$
$\displaystyle
M_{4}=A^{-1}(A^{-1}+\widetilde{A}^{-1})^{-1}\widetilde{A}^{-1}.$ (56)

In fact, all of the $M_{i}$ are equal, as we show in the following lemma.

######  Lemma A.1.

$M_{1}=M_{2}=M_{3}=M_{4}=(A+\widetilde{A})^{-1}$.

Proof:  By symmetry of the $M_{i}$ it suffices to prove the claim for $M_{1}$
and $M_{3}$. Note

$M_{3}=(A(A^{-1}+\widetilde{A}^{-1})\widetilde{A})^{-1}=(A+\widetilde{A})^{-1}.$
(57)

Additionally,

$\displaystyle M_{1}(A+\widetilde{A})$
$\displaystyle=\mathbbm{1}+A^{-1}\widetilde{A}-A^{-1}((A^{-1}\widetilde{A}+\mathbbm{1})\widetilde{A}^{-1})^{-1}(\mathbbm{1}+A^{-1}\widetilde{A}),$
(58) $\displaystyle=\mathbbm{1}.$ (59)

The result follows. ∎

Combining equations (51), (55) and lemma A.1, we can make a change of
coordinates in the integral to find

$f(z)=f(\lambda_{1},\mu_{1},\dots,\lambda_{\ell},\mu_{\ell})\propto\exp(-\frac{1}{2}\left((\lambda_{0}+\mu_{0})^{T}(A+\widetilde{A})^{-1}(\lambda_{0}+\mu_{0})\right)).$
(60)

Setting $\widetilde{r}$ to zero (as in [3]), we reach

$Z|\theta\sim\mathcal{N}(\theta,A+\widetilde{A}),$ (61)

as claimed in Section III.

##  Appendix B Bayesian Gaussian calculations

In this appendix, we find the Bayes estimator for the measurement process
described in Section V, as well as the Bayes risk of the maximum-likelihood,
James-Stein and Bayes estimators

###  B.1 Finding the Bayes estimator

Recall that we are considering the (classical) Bayesian estimation problem of
$Z_{N}\mid\theta\sim\mathcal{N}(\theta,\Gamma_{N})$ and
$\theta\sim\mathcal{N}(\theta_{0},\Xi)$. Since we are using least squares loss
[51], the Bayes estimator has the simple form

$\hat{\theta}^{\text{B}}(z)=\mathbb{E}[\;\theta\mid z\;].$ (62)

Thus, consider the distribution of $\theta\mid z$.

$f(\theta|z)=\frac{f(z|\theta)f(\theta)}{\int
f(z|\theta)f(\theta)d^{2n}\theta}.$ (63)

Note

$f(z|\theta)f(\theta)\propto
e^{-\frac{1}{2}[(z-\theta)^{T}\Gamma_{N}^{-1}(z-\theta)+(\theta-\theta_{0})^{T}\Xi^{-1}(\theta-\theta_{0})]}\\\
$ (64)

We complete the square in the exponent, i.e. write it in the form

$(\theta-r_{0})^{T}M(\theta-r_{0})+R.$ (65)

From inspection, it is clear that we need

$M=\Gamma_{N}^{-1}+\Xi^{-1},\qquad
r_{0}=M^{-1}(\Gamma_{N}^{-1}z+\Xi^{-1}\theta_{0}).$ (66)

We deduce

$\theta|z\sim\mathcal{N}(r_{0},M^{-1}),$ (67)

so that
$\hat{\theta}^{\text{B}}(z)=(\Gamma_{N}^{-1}+\Xi^{-1})^{-1}(\Gamma_{N}^{-1}z+\Xi^{-1}\theta_{0})$.
We also find that

$R=z^{T}\Gamma_{N}^{-1}z+\theta_{0}^{T}\Xi^{-1}\theta_{0}-r_{0}^{T}Mr_{0}.$
(68)

In analogy with equation (54) we find that

$R=(z-\theta_{0})(\Gamma_{N}+\Xi)^{-1}(z-\theta_{0}),$ (69)

And thus deduce

$Z_{N}\sim\mathcal{N}(\theta_{0},\Gamma_{N}+\Xi).$ (70)

###  B.2 Calculating Bayes risks

We recall the three estimators we are comparing:

$\hat{\theta}^{\text{MLE}}(z)=z,\qquad\hat{\theta}^{\text{JS}}(z)=\left(\mathbbm{1}-(2q-2)\frac{\Sigma^{-1}}{z^{T}\Sigma^{-2}z}\right)z,\qquad\hat{\theta}^{\text{B}}(z)=(\Gamma_{N}^{-1}+\Xi^{-1})^{-1}(\Gamma_{N}^{-1}z+\Xi^{-1}\theta_{0}).$
(71)

Since $R(\hat{\theta}^{\text{MLE}},\theta)=\text{Tr}(\Gamma_{N})$, it is easy
to see $R_{\pi}(\hat{\theta}^{\text{MLE}})=\text{Tr}(\Gamma_{N})$.

From equation (II.2) we know

$R(\hat{\theta}^{\text{JS}},\theta)=\text{Tr}(\Gamma_{N})-(2q-2)^{2}\mathbb{E}\left[\frac{1}{Z_{N}^{T}\Gamma_{N}^{-2}Z_{N}}\middle|\;\theta\;\right],$
(72)

By definition,

$R_{\pi}(\hat{\theta}^{\text{JS}})=\mathbb{E}_{\theta\sim\pi}[R(\hat{\theta}^{\text{JS}},\theta)],$
(73)

so using the tower law of expectation we find

$R_{\pi}(\hat{\theta}^{\text{JS}})=\text{Tr}(\Gamma_{N})-(2q-2)^{2}\mathbb{E}\left[\frac{1}{Z_{N}^{T}\Gamma_{N}^{-2}Z_{N}}\right],$
(74)

where the expectation is now taken over the full distribution of $Z_{N}$ that
we found in equation (70).

It remains to find the Bayes-risk of the Bayes estimator
$\hat{\theta}^{\text{B}}$. Using by the tower-law of expectation,

$\displaystyle
R_{\pi}(\hat{\theta}^{\text{B}})=\mathbb{E}_{Z_{N}}\left(\mathbb{E}_{\theta|Z_{N}}\left[\norm{\theta-\hat{\theta}^{\text{B}}(Z_{N})}^{2}\middle|\;Z_{N}\;\right]\right)$
(75)

Since $\hat{\theta}^{\text{B}}(z)=\mathbb{E}[\;\theta\mid z\;]$ and
$\theta\mid z$ is distributed according to equation (67), the inner
expectation is the expected squared distance of a multivariate Gaussian
distribution from its mean, which is the trace of the covariance matrix:
$\text{Tr}(\Gamma_{N}^{-1}+\Xi^{-1})^{-1}$. The outer expectation is then the
expectation of a constant so that
$R_{\pi}(\hat{\theta}^{\text{B}})=\text{Tr}(\Gamma_{N}^{-1}+\Xi^{-1})^{-1}$.

In order to compare the Bayes risk of the Bayes estimator with the MLE and JS
risks, it is useful to rewrite this risk using lemma A.1. Taking
$A=\Gamma_{N}^{-1}$ and $\widetilde{A}=\Xi^{-1}$. We find

$R_{\pi}(\hat{\theta}^{\text{B}})=\text{Tr}(\Gamma_{N})-\text{Tr}(\Gamma_{N}(\Gamma_{N}+\Xi)^{-1}\Gamma_{N}).$
(76)

Noting that $\Gamma_{N}$ is symmetric, and $\Gamma_{N}+\Xi$ is positive
definite, we deduce that that
$R_{\pi}(\hat{\theta}^{\text{B}})<R_{\pi}(\hat{\theta}^{\text{MLE}})$.

##  Appendix C Postselection

###  C.1 Approximate postselected state

In this appendix, we justify equation (31) in Section VI. First we calculate
the overlap

$\displaystyle\innerproduct{\psi_{\theta_{0}}}{\psi_{\theta}}$
$\displaystyle=C^{n}\int_{-\infty}^{\infty}d^{n}x\;\;e^{-(\norm{x-\theta}^{2}+\norm{x-\theta_{0}}^{2})/B},$
(77)
$\displaystyle=C^{n}\int_{-\infty}^{\infty}d^{n}x\;\;e^{-(2\norm{x-\theta/2-\theta_{0}/2}^{2}+\norm{\theta-\theta_{0}}^{2}/2)/B},$
(78) $\displaystyle=e^{-\norm{\delta}^{2}/2B}.$ (79)

Then

$\displaystyle K\ket{\psi_{\theta}}$
$\displaystyle=\ket{\psi_{\theta}}-(1-t)e^{-\norm{\delta}^{2}/2B}\ket{\psi_{\theta_{0}}},$
(80) $\displaystyle=C^{n/2}\int
d^{n}x\;(e^{-\norm{x-\theta}^{2}/B}-(1-t)e^{-\norm{\delta}^{2}/2B}e^{-\norm{x-\theta_{0}}^{2}/B})\ket{x},$
(81) $\displaystyle=C^{n/2}\int
d^{n}y\;(e^{-\norm{y-\delta}^{2}/B}-(1-t)e^{-\norm{\delta}^{2}/2B}e^{-\norm{y}^{2}/B})\ket{y+\theta_{0}},$
(82) $\displaystyle=C^{n/2}\int
d^{n}y\;e^{-\norm{\delta}^{2}/B}e^{-\norm{y}^{2}/B}[e^{2y^{T}\delta/B}-(1-t)e^{\norm{\delta}^{2}/2B}]\ket{y+\theta_{0}}.$
(83)

We can see from equation (82) that the integrand only has effective support
for $y=O(\delta)\ll 1$. Thus we can Taylor expand in $\delta$, treating $y$ as
$O(1)$.

$\displaystyle K\ket{\psi_{\theta}}$ $\displaystyle=C^{n/2}\int
d^{3}y\;e^{-\norm{\delta}^{2}/B}e^{-\norm{y}^{2}/B}[t+2y^{T}\delta/B+O(\norm{\delta}^{2})]\ket{y+\theta_{0}},$
(84) $\displaystyle=tC^{n/2}\int
d^{n}y\;e^{-\norm{\delta}^{2}/t^{2}B}e^{-\norm{y}^{2}/B}e^{2y^{T}\delta/tB}[1+O(\norm{\delta}^{2}/t^{2})]\ket{y+\theta_{0}},$
(85) $\displaystyle=tC^{n/2}\int
d^{n}y\;e^{-\norm{y-\delta/t}^{2}/B}[1+O(\norm{\delta}^{2}/t^{2})]\ket{y+\theta_{0}},$
(86) $\displaystyle=tC^{n/2}\int
d^{n}x\;e^{-\norm{x-\theta_{0}-(\theta-\theta_{0})/t}^{2}/B}[1+O(\norm{\delta}^{2}/t^{2})]\ket{x}.$
(87)

So for $\norm{\delta}\ll t$, we find, as claimed in the main text, that

$\mathbb{P}(\text{pass filter})=\norm{K\ket{\psi_{\theta}}}^{2}\approx
t^{2},\qquad\frac{K\ket{\psi_{\theta}}}{\norm{K\ket{\psi_{\theta}}}}\approx
C^{n/2}\int
d^{n}x\;e^{-\norm{x-\theta_{0}-(\theta-\theta_{0})/t}^{2}/B}\ket{x}.$ (88)

###  C.2 Iterative Strategy

In this appendix, we describe our iterative strategy for estimating $\theta$
using postselection. We wish to compare the James-Stein estimator and MLE. In
order to capture both cases, we describe the algorithm in more general terms,
referring only to an estimator $\hat{\theta}^{\text{}}$ for estimating the
mean of a Gaussian distribution. If one wishes to implement a specific
strategy, one can replace $\hat{\theta}^{\text{}}$ by, for example,
$\hat{\theta}^{\text{MLE}}$ or $\hat{\theta}^{\text{mJS}}$.

Initially, we set $\theta_{0}=0$ and $t=1$.

Suppose we are on the $k$th measurement; $\theta_{0}$ and $t$ have been set to
some known values. Recall (from appendix C.1), that the postselected state is
roughly

$\frac{K\ket{\psi_{\theta}}}{\norm{K\ket{\psi_{\theta}}}}\approx C^{n/2}\int
d^{n}x\;e^{-\norm{x-\theta_{0}-(\theta-\theta_{0})/t}^{2}/B}\ket{x}.$ (89)

We measure the position of the postselected state, let the outcome be $X_{k}$.
Rescaling, we take $Y_{k}=\theta_{0}+t(X-\theta_{0})$ so that
$Y\sim\mathcal{N}(\theta,t^{2}(B/4)\mathbbm{1})$. We make an estimate
$\hat{\theta}^{\text{}}_{k}(Y)$ of $\theta$. Our collated estimate (using all
the data so far) is given by $\bar{\theta}_{k}$, the sample mean of
$\hat{\theta}^{\text{}}_{1},\dots,\hat{\theta}^{\text{}}_{k}$. Denote the
sample standard deviation of
$\hat{\theta}^{\text{}}_{1},\dots\hat{\theta}^{\text{}}_{k}$ by
$\hat{\sigma}_{k}$. We estimate $\delta$ by
$\hat{\delta}_{k}=\hat{\sigma}_{k}/\sqrt{k}$.

If $\delta_{k}$ is ever less than $0.3t$, we set $t$ to $3\hat{\delta}_{k}$
and $\theta_{0}$ to $\bar{\theta}_{k}$: our current best estimate of $\theta$.

If the strategy terminates after measurement $K$, we output
$\bar{\theta}_{K}$. Thus, after measurement $k$, the current error of the
estimation strategy is given by $R_{k}=\norm{\bar{\theta}_{k}-\theta}^{2}$. To
estimate the risk of a particular strategy
$\hat{\theta}^{\text{}}=\hat{\theta}^{\text{mJS}}$ or
$\hat{\theta}^{\text{}}=\hat{\theta}^{\text{MLE}}$ at measurement $k$, we
average $R_{k}$ over many individual repetitions.

###  C.3 Exact sampling from the postselected pdf

In order to numerically simulate the postselection, we must sample exactly
from the distribution

$f_{\theta_{0},t}(x)=\frac{|\matrixelement{x}{K}{\psi_{\theta}}|^{2}}{\norm{K\ket{\psi_{\theta}}}^{2}},$
(90)

given by measuring the position of the postselected state.

In order to do this, we use rejection sampling. We give a brief summary of
rejection sampling here, for details see Ref. [52]. Assume there is a target
pdf $f(x)$ that we wish to sample from. Rejection sampling requires a pdf
$g(x)$ that one can sample from, and constant $M$ satisfying $f(x)\leq Mg(x)$.
The rejection sampling algorithm is:

  1. 1.

Generate a candidate $x$ from the pdf $g(x)$.

  2. 2.

Output this candidate with probability $f(x)/Mg(x)$. If the candidate is
rejected, return to step 1.

Note that the chance of acceptance is given by $\int g(x)[f(x)/Mg(x)]dx=1/M$.
Thus, on average, $M$ samples from $g(x)$ are required to generate one sample
from $f(x)$.

To use rejection sampling, we must find a constant $M_{\theta_{0},t}$ and a
pdf $g_{\theta_{0},t}(x)$ that is easy to sample from, such that
$f_{\theta_{0},t}(x)\leq M_{\theta_{0},t}g_{\theta_{0},t}(x)$. First, we
calculate

$\displaystyle\innerproduct{\psi_{\theta}}{\psi_{\theta_{0}}}=e^{-\norm{\delta}^{2}/2B},$
(91)

where $\delta=\theta-\theta_{0}$. Then, recalling
$K=(t-1)\outerproduct{\theta_{\theta_{0}}}{\theta_{\theta_{0}}}+\mathbbm{1}$,
we find

$\displaystyle|\matrixelement{x}{K}{\psi_{\theta}}|^{2}=C^{n}|e^{-\norm{x-\theta}^{2}/B}-(1-t)e^{-\norm{\delta}^{2}/2B}e^{-\norm{x-\theta_{0}}^{2}/B}|^{2}.$
(92)

By expanding the square and integrating each of the resulting Gaussians, we
find that

$\displaystyle\norm{K\ket{\psi_{\theta}}}^{2}$ $\displaystyle=\int
dx|\matrixelement{x}{K}{\psi_{\theta}}|^{2},$ (93)
$\displaystyle=1+(t^{2}-1)e^{-\norm{\delta}^{2}/B}.$ (94)

If $A,B\geq 0$, note that $|A-B|^{2}\leq A^{2}+B^{2}$. Thus, we may bound
equation (92) by

$|\matrixelement{x}{K}{\psi_{\theta}}|^{2}\leq
C^{n}(e^{-2\norm{x-\theta}^{2}/B}+(1-t)^{2}e^{-\norm{\delta}^{2}/B}e^{-2\norm{x-\theta_{0}}^{2}/B})$
(95)

Thus,

$\displaystyle f_{\theta_{0},t}(x)$ $\displaystyle\leq
C^{d}\frac{e^{-2\norm{x-\theta}^{2}/B}+(1-t)^{2}e^{-\norm{\delta}^{2}/B}e^{-2\norm{x-\theta_{0}}^{2}/B}}{1+(t^{2}-1)e^{-\norm{\delta}^{2}/B}},$
(96)
$\displaystyle=\frac{1+(1-t)^{2}e^{-\norm{\delta}^{2}/B}}{1+(t^{2}-1)e^{-\norm{\delta}^{2}/B}}\left(\frac{1}{1+(1-t)^{2}e^{-\norm{\delta}^{2}/B}}C^{d}e^{-2\norm{x-\theta}^{2}/B}+\frac{(1-t)^{2}e^{-\norm{\delta}^{2}/B}}{1+(1-t)^{2}e^{-\norm{\delta}^{2}/B}}C^{n}e^{-2\norm{x-\theta_{0}}^{2}/B}\right)$
(97) $\displaystyle=M_{\theta_{0},t}g_{\theta_{0},t}(x),$ (98)

where

$M_{\theta_{0},t}=\frac{1+(1-t)^{2}e^{-\norm{\delta}^{2}/B}}{1+(t^{2}-1)e^{-\norm{\delta}^{2}/B}},$
(99)

and

$g_{\theta_{0},t}(x)=\frac{1}{1+(1-t)^{2}e^{-\norm{\delta}^{2}/B}}C^{n}e^{-2\norm{x-\theta}^{2}/B}+\frac{(1-t)^{2}e^{-\norm{\delta}^{2}/B}}{1+(1-t)^{2}e^{-\norm{\delta}^{2}/B}}C^{n}e^{-2\norm{x-\theta_{0}}^{2}/B}.$
(100)

Note that $g_{\theta_{0},t}(x)$ is a pdf given by the convex sum of the pdf of
two different normal distributions. Specifically, let

$\phi_{\mu}(x)=C^{n}e^{-2\norm{x-\mu}^{2}/B},$ (101)

denote the pdf of a normal distribution with mean $\mu$ and covariance matrix
$\frac{B}{4}\mathbbm{1}$. Then,

$g_{\theta_{0},t}(x)=p_{\theta_{0},t}\phi_{\theta}(x)+(1-p_{\theta_{0},t})\phi_{\theta_{0}}(x),$
(102)

where

$p_{\theta_{0},t}=\frac{1}{1+(1-t)^{2}e^{-\norm{\delta}^{2}/B}}.$ (103)

Thus it is easy to sample from $g_{\theta_{0},t}(x)$: sample from
$\phi_{\theta}(x)$ with probability $p_{\theta_{0},t}$, otherwise sample from
$\phi_{\theta_{0}}(x)$. Therefore, we can use rejection sampling to sample
from $f_{\theta_{0},t}$. Note that as $t,\delta\to 0,M\to\infty$, so as
postselection becomes stronger, rejection sampling takes longer (on average)
to run.