\section{Introduction}
Quantum computing is much in the news these days, with the recent achievement of quantum supremacy by the Google device \cite{Naturegoogle} and the renewed interest in building a quantum computer (see, for example, Refs.~\cite{Nature1,Nature2}). The main obstacle to realizing a quantum computer of a useful scale and accuracy is the noise that threatens to destroy the quantum features that give quantum computers their power. A key focus of any implementation of a quantum information processing device is hence the characterization, and subsequent control and mitigation, of the noise that unavoidably accompanies the operation of a quantum gate.

An often used measure of the quality of a quantum gate implementation is the fidelity of the gate, i.e., a quantification of how close the action of the actual gate is to the theoretical ideal. Efficient and easy-to-implement procedures such as randomized benchmarking \cite{RB1,RB2,RB3} can offer information on the \emph{average} fidelity of a gate, averaged over input states according to some distribution of interest (e.g., Haar-distributed pure states for randomized benchmarking). What is often more telling about the gate performance, however, is the \emph{minimum} fidelity, i.e., the fidelity of the gate operation, minimized over all possible input states. This is a state-distribution-independent quantity, and gives a minimum guarantee for the quality of the gate, which is important to assure that the quantum computer functions correctly in all scenarios.

Conventionally, the minimum fidelity is estimated by first a full quantum process tomography (QPT) of the noisy gate operation, followed by a numerical minimization, using the estimated process matrix, to find the minimum fidelity over all possible input states. This is a potentially costly procedure: For a $d$-dimensional system, QPT requires an accurate estimation of $d^4$ real parameters, all to yield one number, i.e., the minimum fidelity. Resource estimates of the number of uses of the channel required for an accurate full QPT can be deduced from bounds for state tomography \cite{Kueng2014,SampleComplexity}, giving $O(d^6/\epsilon'^2)$, where $\epsilon'$ is the accuracy of the reconstruction, measured by a distance between quantum processes.

One expects to be able to do better if the desire is only to estimate the minimum fidelity to some desired accuracy $\epsilon$, without the full reconstruction of the process matrix of the noisy gate. In this work, we propose a direct route to estimating the minimum fidelity, without the use of QPT. It combines a numerical gradient-free minimization (GFM) algorithm with the direct target fidelity estimation scheme of Refs.~\cite{DF1,DF2} to perform a hybrid numerical-experimental descent of the fidelity function to the minimum. We compare this GFM scheme with an alternative route, which we name \emph{QPT fidelity estimation}. QPT fidelity estimation does go through QPT, but rather than employing a process distance as the figure of merit as done in standard QPT, we make direct use of the minimum fidelity as the stopping rule. In this sense, this second route can also be regarded as a direct estimation of the minimum fidelity and offers fair comparison with the GFM scheme.

Resource estimates suggest that QPT fidelity estimation needs $O(d^8/\epsilon^2)$ uses of the gate to achieve an estimate of the minimum fidelity to accuracy $\epsilon$. This follows from an argument (see Sec.~\ref{sec:theory}) that relates the desired minimum fidelity accuracy of $\epsilon$ to the Choi-state trace-distance accuracy of $\epsilon'$ used in Refs.~\cite{Kueng2014,SampleComplexity}. This $O(d^8)$ uses of the gate, for fixed $\epsilon$, is a prohibitively high resource cost. The GFM approach, in contrast, is expected (see Sec.~\ref{sec:theory}) to require much milder $O(d^4)$ gate uses. It is thus, from this perspective, a much more efficient procedure for estimating the minimum fidelity.

For specific classes of noisy gates, however, our numerical results suggest a somewhat different conclusion. We numerically simulated both GFM and QPT fidelity estimation procedures for two natural classes of noisy gates: noisy gates random in the Hilbert--Schmidt sense, and gates with random Pauli and amplitude-damping noise. For both classes, our GFM scheme performs close to (though better than) the resource estimates; the QPT fidelity estimation scheme, surprisingly, took only $O(d^4)$ gate uses, much fewer than the $O(d^8)$ prediction. The performance of the two schemes is hence comparable for these two classes of noisy gates, with only slightly better scaling for the GFM scheme. This better performance for QPT is especially unexpected as our simulations do not follow the optimal procedure behind the theoretical bounds obtained in Refs.~\cite{Kueng2014,SampleComplexity}. This reminds us that, for specific classes of noisy gates, the resource estimates, which account even for worst-case scenarios, may not provide a reliable gauge of typical performance.

\makeatletter\global\advance\@colroom-10pt\relax\set@vsize\makeatother

Nevertheless, while there is little room for improvement in the resource scaling for QPT fidelity estimation, the scaling for the GFM scheme is limited largely by the efficiency of the numerical GFM algorithm, a subject of intense study in the field of numerical optimization. Our work also emphasizes the importance of making direct use of the quantity of interest, in this case the minimum fidelity, in the measurement procedure, rather than a secondary quantity such as a process distance. Both schemes explored here, the GFM scheme as well as the QPT fidelity estimation, can be considered when looking for methods of direct estimation of the minimum gate fidelity, with the GFM offering potential for further reduction in resource cost, while the QPT fidelity estimation offers the advantage of a measurement setup familiar from standard QPT.

Whether the observation of Ref.~\cite{Huang+2:2020}, namely, that very
  few measurements, chosen at random, may be sufficient for the estimation of
  particular properties, has a bearing on estimating the minimum gate fidelity
  is currently unknown, and this deserves to be explored.
  The matter is, however, not within the scope of this work.

Below, we begin in Sec.~\ref{sec:directEst} with a description of our GFM algorithm, also providing a reminder of the direct target fidelity estimation scheme following the analysis of Ref.~\cite{DF1}. We then explain the QPT fidelity estimation procedure in Sec.~\ref{sec:QPT}. Section \ref{sec:perf} gives the resource estimates and numerical simulations that compare the performance of the two schemes. We conclude in Sec.~\ref{sec:conc}.

\section{Direct estimation with gradient-free minimization}\label{sec:directEst}
For a known ideal gate $\mathcal{G}$, acting on an $n$-qubit system, and its noisy implementation $\widetilde\cG$, we are interested in how close $\widetilde\cG$ is to $\mathcal{G}$. One way to quantify this is to compare the fidelity between the output state of $\widetilde\cG$ with the ideal output state after $\mathcal{G}$. The minimum fidelity, minimized over all input states, is defined as
\begin{equation}\label{eq:Fmin}
F_{\min}\equiv\min_{\psi}F\bigl(\mathcal{G}(\psi),\widetilde\cG(\psi)\bigr)=\min_{\psi}F\bigl(\psi,\mathcal{E}(\psi)\bigr).
\end{equation}
Here, $\psi\equiv |\psi\rangle\langle\psi|$ is a pure state on the $d(=\!2^n)$-dimensional Hilbert space of the quantum system, $F(|\phi\rangle\langle\phi|,\rho)\equiv \langle\phi|\rho|\phi\rangle$ denotes the squared fidelity between the pure state $|\phi\rangle$ and the (possibly mixed) state $\rho$, and $\mathcal{E}\equiv \mathcal{G}^{-1} \circ\widetilde\cG$ is the noise process that describes the imperfections in the gate implementation $\widetilde\cG$. $\mathcal{G}^{-1}$ is the inverse of the ideal unitary gate $\mathcal{G}(\cdot)=U(\cdot)U^\dagger$, for unitary $U$, i.e., $\mathcal{G}^{-1}(\cdot)=U^\dagger(\cdot) U$. We assume that $\widetilde\cG$ is a completely positive (CP) and trace-preserving (TP) map---also referred to as a quantum channel---and hence so is $\mathcal{E}$. Even though the minimization in Eq.~\eqref{eq:Fmin} appears to only be over pure states $\psi$, $F_{\min}$ is, in fact, the minimum fidelity over all states, pure or mixed, as the concavity of $F$ ensures that the minimum is attained on a pure state.

The goal here is to estimate $F_{\min}$ without first estimating the full process matrix of the noisy gate $\widetilde\cG$, or, equivalently, of the noise channel $\mathcal{E}$. We assume the following experimental capabilities: (i) We can prepare any input state of our choice; (ii) we can send that input state through the noisy gate and access the output state; (iii) we can perform product-Pauli measurements on the output state. The noisy gate is regarded here as a black box $\widetilde\cG$ in the laboratory that takes the input $\psi$ and gives back the output $\widetilde\cG(\psi)$.

We need two additional ingredients. The first is the technique of direct estimation of target fidelity (DTFE) invented in Ref.~\cite{DF1}, allowing the estimation of the fidelity of an $n$-qubit state $\rho$ with some target pure state $|\psi\rangle$ without full state tomography. One needs only make probabilistic product-Pauli measurements, according to a distribution determined by the target pure state (see Sec.~\ref{sec:DTFE} for more details).

The second ingredient is a numerical method for gradient-free minimization (GFM), implementable on a classical computer. A GFM method finds a local minimum of a function in the case where function values, but not the gradient values, are easily available as inputs to the algorithm. A common situation is one where the function itself cannot be written down explicitly, but is accessible only through a numerical procedure. The gradient of the function is hence also not available as a function that can be written down explicitly, and methods of numerically approximating the gradient typically do not work well or are prohibitively expensive to evaluate.
GFM methods incorporate the gradient estimation with the minimization, by choosing trial points in the domain space, deducing some gradient information (often only a rough estimate) from the function values evaluated at those points, taking a small step in a direction expected to lower the function value according to that gradient information, and then repeating the process with new trial points at the new location. With standard regularity criteria on the function, one eventually arrives at a local minimum, and a repetition of the procedure sufficiently many times gives a good chance of finding the global minimum.

We now put the ingredients together for our GFM scheme of direct estimation of the minimum fidelity. Again, the goal is to estimate $F_{\min}$ without full knowledge of $\widetilde\cG$, apart from access to it as an input-output black box in the laboratory. As is typical in such problems, the requirement is to estimate $F_{\min}$ to within a target accuracy of $\epsilon$ from the true value with high probability, i.e., $|\widehatF_{\min}-F_{\min}|\leq \epsilon$ for some small $\epsilon>0$, where $\widehatF_{\min}$ is the estimate while $F_{\min}$ is the (unknown) true value. The function we are minimizing here is $f(\psi)\equiv F\bigl(\psi,\mathcal{E}(\psi)\bigr)$, over the domain of pure states $\psi$. A full description of $f$ and its gradient for minimization using gradient-descent methods requires knowledge of $\mathcal{E}$, which we do not possess. However, we can evaluate the function value in the laboratory by preparing $\psi$, feeding it into the black box $\widetilde\cG$, and then estimating $f(\psi)$ by carrying out the DTFE scheme in the laboratory, treating $\mathcal{G}(\psi)$ as the target (pure) state and $\widetilde\cG(\psi)$ as the state to be compared with the target. In this way, we have access, from measurements in the laboratory on the black box $\widetilde\cG$, to the function values $f(\psi)$.

Our scheme then proceeds as follows: We run, on a classical computer, a GFM algorithm to minimize $f(\psi)$. At each iterative step, the GFM suggests a set of trial states $\psi$, for which it requires the values $f(\psi)$. These values are obtained from the experiment by preparing those $\psi$ states and then performing the DTFE algorithm in the laboratory for each $\psi$. We iterate the GFM algorithm until a stopping rule is satisfied to confirm the attainment of a local minimum. The whole procedure is repeated sufficiently many times to find the global minimum with a specified (high) probability. The stopping rule is carefully tuned, as we will describe below, to attain the desired $\epsilon$ accuracy. The scheme is summarized in Fig. \ref{fig:scheme}.

\makeatletter\global\advance\@colroom10pt\relax\set@vsize\makeatother

\begin{figure}[t]
\centering
\includegraphics[width=\columnwidth]{DEOMGF-Fig1.pdf}
\caption{\label{fig:scheme}Schematic depictions of the two approaches explored in our paper: GFM scheme and QPT fidelity estimation. Both approaches involve interwoven experimental and numerical steps. In the GFM approach, the experimental---via the direct target fidelity estimation procedure---and the numerical---via gradient-free descent---steps alternate seamlessly in a hybrid minimization algorithm, until convergence to the minimum fidelity $\widehatF_{\min}$. In QPT fidelity estimation, the experimental reconstruction of the process matrix occurs first, followed by a numerical minimization to obtain $\widehatF_{\min}$, and this is repeated, if needed, with more uses of the gate until the estimated accuracy of $F_{\min}$ meets our set target.}
\end{figure}

In the following two sections, we describe, in greater detail, first the DTFE scheme of Ref.~\cite{DF1} and then the GFM algorithm employed in our numerical examples.

\subsection{Direct estimation of the target fidelity}\label{sec:DTFE}
We describe here the basic ideas of the DTFE scheme of Ref.~\cite{DF1} needed to understand our work. The analysis and formulas here are taken from that reference, and we refer the reader to the original paper for further details. The target fidelity, i.e., the quantity of interest in Ref.~\cite{DF1}, is the fidelity between a particular state $\rho$, imagined to be the actual state produced by a source, and a pure state $|\psi\rangle$, taken to be the target state the source is designed to produce in the ideal situation. The target fidelity hence quantifies how close the source is to that ideal. That the target state is pure is important technically for the linearity, in $\rho$ and $\psi$, of the squared fidelity used as the logical basis of the scheme; a pure target state is anyway the often-encountered situation in many quantum information processing tasks. For our purpose here, it suffices as well, as our minimization is over pure input states only.

We write the square of the target fidelity, for a $d=2^n$-dimensional situation, as
\begin{equation}\label{eq:fid1}
F(\psi,\rho)=\mathrm{tr}(\psi\rho)=\sum_kx_k(\psi)x_k(\rho),
\end{equation}
where $x_k(\,\cdot\,)\equiv \frac{1}{\sqrt d}\mathrm{tr}\bigl\{W_k(\,\cdot\,)\bigr\}$, with $W_k$ an $n$-qubit product Pauli operator, so that the set $\bigl\{\frac{1}{\sqrt d}W_k\bigr\}_{k=1}^{d^2}$ is an orthonormal operator basis. $x_k(\tau)$ for any state $\tau$ is then the (real) coefficient of $\frac{1}{\sqrt d}W_k$ when writing $\tau$ as a linear combination of elements of that Pauli operator basis. The sum over $k$ in Eq.~\eqref{eq:fid1} above is understood to be only over those $k$ values with $x_k(\psi)\neq 0$. In our present situation, the $x_k(\psi)$s are known, while the $x_k(\rho)$s are not.

The fidelity can be rewritten as
\begin{equation} \label{eq:fid2}
F(\psi,\rho)=\sum_k p_kX_k,
\end{equation}
with $p_k\equiv [x_k(\psi)]^2$ and $X\equiv\frac{x_k(\rho)}{x_k(\psi)}$. Note that $\{p_k\}$ is a probability distribution: $\sum_kp_k=F(\psi,\psi)=1$, and $p_k>0$. Let us define a random variable $X$ which takes value $X_k$ with probability $p_k$. Then, observe that $F(\psi,\rho)=\langle X\rangle$, the expectation value of $X$, which depends on the unknown $\rho$. Now, for each $k$, $x_k(\rho)$ can be estimated from the experiment by measuring $W_k$, a product-Pauli measurement, on $\rho$. Then, the target fidelity, now understood to be equal to $\langle X\rangle$, can be estimated by repeated trials where $W_k$ is chosen as the measurement to be performed on $\rho$ with probability $p_k$.

How good is this estimate? Suppose we do $h\equiv 1/(\eta\delta^2)$ trials, with $k_1,k_2,\ldots, k_h$ chosen according to the distribution $\{p_k\}$, and so obtain estimates for $X_{k_1},X_{k_2},\ldots,X_{k_h}$. Then, $Y\equiv \frac{1}{h}\sum_{\ell=1}^hX_{k_\ell}$ satisfies the inequality: $\textrm{Probability}{\left(|Y-\langle X\rangle|\geq\eta\right)}\leq \delta$.
We cannot, however, know $X_{k_\ell}$ precisely with only a finite number of copies of $\rho$. To estimate $X_{k_\ell}$, it suffices to measure $W_{k_\ell}$ on $t_{k_\ell}$ copies of $\rho$, where
\begin{equation}
t_{k_\ell}\equiv \frac{2\log (2/\delta)}{dp_{k_\ell}h\eta^2}.
\end{equation}
Then, $Y$, now built from estimates of $X_{k_\ell}$ using the above number of copies, satisfies instead the inequality
\begin{equation}
\textrm{Probability}{\left(|Y-\langle X\rangle|\geq2\eta\right)}\leq 2\delta.
\end{equation}
The expected total number of copies of $\rho$ required, to attain an estimate of the target fidelity $2\eta$ away from the true value is then
\begin{equation}
h\sum_kp_kt_k= \sum_k\frac{2\log (2/\delta)}{d \eta^2}\leq 2\log(2/\delta)\frac{d}{\eta^2},
\end{equation}
giving a scaling of $O(d/\eta^2)$ copies, for fixed $\delta$. A more accurate estimate of the expected number of copies is given in Ref.~\cite{DF1} by taking into account that $t_{k_\ell}$ has to be an integer, but that does not change the $O(d/\eta^2)$ scaling with $d$.

\subsection{Gradient-free minimization with CMA-ES}\label{CMAES}
Our GFM direct estimation scheme accommodates the use of any GFM algorithm. The efficiency of the GFM algorithm is of crucial importance in the present context, but other considerations such as ease of coding and number of tuning parameters can also affect one's choice.
There are many known GFM algorithms, including the downhill simplex method \cite{Simplex,ComplexitySimplex}, the directional direct search \cite{dds1,dds2,dds3}, the stochastic method \cite{GF5,SG1}, and the covariance matrix adaptation evolution strategy \cite{CMAES1,CMAES2,CMAES3,CMAES4} (CMA-ES), among others. The downhill simplex method is arguably the most well-studied one, but may not perform well in high-dimensional problems \cite{Torczon1989}. Trying out a few methods on our problem, we found CMA-ES to work well, converging more quickly than the downhill simplex method in our problem instances, with more stable performance than the directional direct search, and having fewer tuning parameters than the stochastic method. We thus focus on CMA-ES in our numerical examples below. We emphasize that a different user can choose a different GFM algorithm and observe a different efficiency performance; therein lies the potential for improvement beyond what we report here, with more efficient GFM methods.

Here, we provide some pertinent details of the CMA-ES algorithm, referring the reader to the original papers (Refs.~\cite{CMAES1,CMAES2,CMAES3,CMAES4}) for further explanation. We pay special attention to how we choose the stopping criteria that terminates the CMA-ES algorithm, as they are crucial for attaining the desired accuracy for $F_{\min}$.

To implement the CMA-ES algorithm, we parameterize the domain space, the $d$-dimensional Hilbert space of pure ($n$-qubit) states, as $\psi=\ell \ell^\dagger/\mathrm{tr}(\ell\ell^\dagger)$, where $\ell$ is the complex column vector
\begin{equation}
\ell=(\ell_1+\mathrm{i} \ell_2, \ell_3+\mathrm{i} \ell_4,\ldots, \ell_{2d-1}+\mathrm{i} \ell_{2d})^\mathrm{T}.
\end{equation}
The domain space is $2d$-dimensional, but there is a single extra parameter---the trace of $\psi$, constrained to be 1---when regarded as a parameterization for $\psi$. This extra parameter does no harm to the minimization, and, in fact, we observe this $2d$-parameter approach to work better for CMA-ES in our numerical examples than an alternative $(2d-1)$-dimensional parameterization with spherical coordinates.

In the $k$th iterative step of the CMA-ES algorithm, a set of $\lambda$ points in the domain space, $L^{(k)}\equiv{\left\{\ell_a^{(k)}\equiv\bigl(\ell_{a,1}^{(k)},\ell_{a,2}^{(k)},\cdots,\ell_{a,2d}^{(k)}\bigr)^\mathrm{T}\right\}}_{a=1}^\lambda$, is drawn from the normal distribution,
\begin{equation}
\ell_a^{(k)} \sim m^{(k)} + \sigma^{(k)}\mathcal{N}{\bigl(0,C^{(k)}\bigr)}, \quad \text{for} \; a=1,...,\lambda,
\end{equation}
where $m^{(k)}$, $\sigma^{(k)}$, and $C^{(k)}$ are the mean, step-size, and covariance matrix, respectively, for the $k$th step. At the initialization step, $m^{(1)}$ is set as a Haar-random pure state, $C^{(1)}$ is set equal to the identity matrix, and $\sigma^{(1)}$ is set to be 0.3. $m^{(1)}$ is set to be a Haar-random pure state to reflect our lack of knowledge of the state which attains the minimum value of fidelity. We observed empirically that the performance remains roughly the same even if we increase or decrease $\sigma^{(1)}$ by an order of magnitude. As long as $\sigma^{(1)}$ is not set to be too large or too small, the performance will not be affected too much. For the case with many local minima, choosing a larger value of $\sigma^{(1)}$ might increase the chance of finding the global minimum. The fidelity value $f_a^{(k)}\equiv f\bigl(\psi_a^{(k)}\bigr)$, for $\psi_a^{(k)}$ built from $\ell_a^{(k)}$, is estimated using the DTFE scheme for each $\ell_a^{(k)}\in L^{(k)}$. The points are then ranked according to the $f_a^{(k)}$ values, with $\ell_{a:\lambda}^{(k)}$ having the $a$th smallest value of $f$, i.e., $f\bigl(\psi_{1:\lambda}^{(k)}\bigr)\leq f\bigl(\psi_{2:\lambda}^{(k)}\bigr)\leq\cdots\leq f\bigl(\psi_{\lambda:\lambda}^{(k)}\bigr)$. The mean for the normal distribution is then updated by the weighted average,
\begin{equation}\label{eq:update}
m^{(k+1)} = \sum\limits^{\mu}_{a=1}w_a \ell^{(k)}_{a:\lambda},
\end{equation}
where $w_k>0$ and $\mu\leq\lambda$. The step-size and covariance matrix are also updated according to rules based on the ranking of the $f$ values. Further details of the algorithm, as well as the appropriate choice of parameters, can be found in Ref.~\cite{Hansen2016}.

In effect, the updates of $m^{(k)}$, $\sigma^{(k)}$, and $C^{(k)}$ move the ``region of interest'' within the domain space in the direction of smaller $f$ values, in correspondence with our goal of minimizing $f$. In the next iterative step, we draw the domain points $L^{(k+1)}$ from the updated region of interest and continue the move towards small $f$ values. This continues until the stopping criteria (more on these below) are met. Our estimate for $F_{\min}$, for one run of CMA-ES, is then the fidelity for the best (smallest $f$) point found in the final iterative step, $\ell^{(\ldots)}_{1:\lambda}$. Each run of the CMA-ES algorithm returns an estimate of a local minimum value for $f$; the entire procedure is repeated sufficiently many times, with different initial conditions, to have confidence that one of those local minima is the global minimum. We estimate the number of repeats needed, for a $95\%$ chance of obtaining the global minimum, using random channels in our numerical simulation (see the numerical examples in Sec.~\ref{Examples}).

To determine when to stop the CMA-ES descent, we calculate two quantities in each iterative step. The first is the range of fidelities for the current set of sample points $L^{(k)}$,
\begin{equation}
\Delta f^{(k)} \equiv \max\limits_a \{f^{(k)}_a\}-\min\limits_a \{f^{(k)}_a\}.
\end{equation}
The second is
\begin{equation}
\nabla f^{(k)} \equiv \frac{1}{\lambda}\sum\limits_{a=1}^\lambda \frac{|f^{(k)}_c - f^{(k)}_a|}{|\ell^{(k)}_c - \ell^{(k)}_a|},
\end{equation}
where $\ell^{(k)}_c$ is the centroid of $L^{(k)}$, and $f^{(k)}_c$ is the corresponding fidelity. $\nabla f^{(k)}$ can be taken as an estimation of the magnitude of the gradient at $\ell^{(k)}_c$.
The stopping criteria are then imposed adaptively to accommodate the variety of possible behaviors of the descent function $f(\psi)$ for different noisy channels. We set an initial threshold $g^{(1)}$ for the gradient estimate $\nabla f^{(1)}$. The target $F_{\min}$ accuracy is $\epsilon$, as before. In subsequent iterative steps, if $\nabla f^{(k)}> g^{(k)}$, we retain the gradient threshold level, i.e., set $g^{(k+1)}=g^{(k)}$, and continue with the next iterative step. If, instead, $\nabla f^{(k)}<g^{(k)}$ but $\Delta f^{(k)}>\epsilon$, we halve the gradient threshold, i.e., set $g^{(k+1)} = g^{(k)}/2$, before taking the next iterative step. If $\nabla f^{(k)}<g^{(k)}$ and $\Delta f^{(k)}<\epsilon$, the algorithm is terminated. A schematic of the stopping criteria is given in Fig.~\ref{fig:GFMStop}.

\begin{figure}[t]
\centering
\includegraphics[width=0.65\columnwidth]{DEOMGF-Fig2.pdf}
\caption{ \label{fig:GFMStop} Diagrammatic representation of the stopping criteria used in the GFM scheme. See main text for a description of the various quantities.}
\end{figure}

The rationale for the adaptive stopping criteria is as follows. The algorithm should terminate when it gets to a stationary point, a minimum, where both the gradient and the range of fidelity values in the set of sample  points are small. That the range is small is demanded by our set target accuracy; that the gradient is small enough is needed to ensure that we have arrived at the minimum value. If we know the $f$ function, requiring the gradient alone to be small enough would also guarantee that the range is small, as the sample points are picked from a local region around the current position. However, since we do not have information about $f(\psi)$, we also do not know how small a gradient threshold is needed for sufficient convergence in the range. Instead, the gradient threshold has to be adjusted adaptively according to the noise channel at hand. We make use of the range of the fidelity values to judge whether the current gradient threshold is sufficiently stringent to ensure that we have arrived at the minimum point. If the gradient is below the set threshold, but the range is still beyond $\epsilon$, this indicates that our gradient criterion is simply not stringent enough; $g^{(k)}$ is thus halved and the iteration continues.

Note that the value of initial gradient threshold $g^{(1)}$ has to be pre-chosen according to the expected channel distribution and the target $\epsilon$ value, to prevent early termination of the algorithm. In our numerical examples, the value of $g^{(1)}$ was chosen, for each distribution of $n$-qubit channels, so that the probability (over the channel distribution) that $|\widehat{F_\text{min}}-F_\text{min}|\leq \epsilon$ is around $0.95$, with $\epsilon=0.01$, with the $0.95$ judged from numerical simulations using a number of trial true channels.

There is one more important detail that has to be addressed: the accuracy of the DTFE procedure needed for the CMA-ES algorithm, as quantified by the $\delta$ and $\eta$ parameters of Sec.~\ref{sec:DTFE}. The DTFE, with a given number of copies of the noisy state, estimates $f(\psi)$ only to a certain accuracy. If the DTFE accuracy is set too high, the resource cost of our GFM direct estimation scheme will become very high; if it is set too low, the CMA-ES algorithm may not converge to the correct minimum value, if it converges at all. CMA-ES, like many evolutionary algorithms of a similar flavor, is relatively robust to noise in function evaluations \cite{Arnold2002}, but proper handling of the accuracy of the function evaluation is still needed.  In our numerical examples, we set $\delta=0.05$ in our DTFE sub-routine and update the needed $\eta$ value in the $k$th iteration of the CMA-ES according to the uncertainty handling algorithm \cite{uncertaintyhandle,uncertaintyhandle2} described next.

The CMA-ES update rule of Eq.~\eqref{eq:update} is based on the ranking of the $f$ values. This means that noisy function evaluation will not affect it as long as the inaccuracies are not large enough to change the $f$ ranking. One can judge whether the function evaluation is sufficiently accurate by checking for rank changes after reevaluations of the function. In our specific context, reevaluation of $f(\psi)$ refers to repeating the DTFE procedure with the same parameters, to get a second estimate of the value of $f(\psi)$ for a given target state $\psi$. The inherent randomness of the DTFE procedure will yield different values, with larger variations from fewer number of uses of the channels, and hence more inaccurate evaluation of $f$.
We follow the uncertainty handling algorithm proposed in \cite{uncertaintyhandle,uncertaintyhandle2}. At each iterative step $k$, we reevaluate $f(\psi_a)$ for each $\psi_a$ built from $\ell_a^{(k)}\in L^{(k)}$, calculate the rank changes after the reevaluations, compute a measure of the uncertainty level, and adjust the accuracy of the DTFE procedure as needed. Specifically, we carry out the following steps:

\begin{center}
\textbf{Uncertainty Handling Algorithm}
\end{center}
\begin{enumerate}
\item For each $L^{(k)}$, compute $f$ for each $\ell_a^{(k)}\in L^{(k)}$ twice. Denote the two values obtained as $f_a^{(k)}$ and $\widetilde{f}_a^{(k)}$.
\item For each $a=1,...,\lambda$, compute the rank change $\Delta R_a^{(k)} = |R(f_a^{(k)})-R(\widetilde{f}_a^{(k)})|-1$, where $R(f_a^{(k)})$ and $R(\widetilde{f}_a^{(k)})$ are the ranks of $f_a^{(k)}$ and $\widetilde{f}_a^{(k)}$, respectively, in the combined set $\{f_a^{(k)},\widetilde{f}_a^{(k)}\}_{a=1}^\lambda$.
\item Compute the uncertainty level,
\begin{align}
\qquad s^{(k)}\!&\equiv\frac{1}{\lambda}\sum_{a=1}^\lambda \Bigl[2\Delta R_a^{(k)}\!-\!\pi_\theta\Bigl(\!R\bigl(\widetilde{f}_a^{(k)}\bigr)\!-\mathrm{H}\bigl(\widetilde{f}_a^{(k)}\!\!-\!f_a^{(k)}\bigr)\!\Bigr)\nonumber\\
&\quad\qquad -\pi_\theta\Bigl(R\bigl(f_a^{(k)}\bigr)-\mathrm{H}\bigl(f_a^{(k)}-\widetilde{f}_a^{(k)}\bigr)\Bigr)\Bigr],
\end{align}
where $\mathrm{H}(\cdot)$ is the step function and $\pi_\theta(R)$ is the $(50\theta)$-th percentile of all possible rank changes (given by the set $\{|1-R|,|2-R|,...,|2\lambda-1-R|\}$) for a given rank $R$.
\item If $s^{(k)}>0$, increase the accuracy in DTFE in the $(k+1)$-th iteration by setting $\eta^{(k+1)} = \alpha\eta^{(k)}$, with $0<\alpha<1$. If $s^{(k)}<0$, set $\eta^{(k+1)} = \eta^{(k)}/\alpha$.
\item For $a=1,...,\lambda$, set $f_a^{(k)} = \tfrac{1}{2}\bigl(f_a^{(k)}+\widetilde{f}_a^{(k)}\bigr)$.
\end{enumerate}
In our numerical examples, $\theta=0.7$ and $\alpha=\frac{1}{\sqrt{2}}$ were found to be good choices.

\section{QPT Fidelity Estimation}\label{sec:QPT}
The performance of our GFM direct estimation scheme has to be compared with the standard alternative of quantum process tomography (QPT), whereby the full process matrix of the unknown channel $\mathcal{E}$ is first estimated, and then the $F_{\min}$ value is numerically computed from the obtained process matrix. For completeness, we remind the reader of a few aspects of QPT important for our work; of course, many textbooks and papers are available on the subject (see, for example, Ref.~\cite{LNP649}).

QPT attempts to reconstruct the full description of an unknown quantum channel (or process) $\mathcal{E}$ from a finite number of uses of the channel. A chosen set of states $\{\rho_k\}$ is sent in as inputs to the channel, and state tomography is done on the outputs of the channel, using a measurement (a positive operator-valued measure, or POVM) $\{\Pi_\ell\}$. The probability of getting outcome $\Pi_\ell$ if the input state $\rho_k$ was sent is given by the Born rule,
\begin{equation}\label{eq:BornRule}
p_{k\ell}=\mathrm{tr}\bigl(\Pi_\ell\mathcal{E}(\rho_k)\bigr).
\end{equation}
One estimates the values of the $p_{k\ell}$s from the experiment using the channel $N_{k}$ times for each input state $\rho_k$, amounting to a total of $N=\sum_{k}N_k$ uses of the channel. The sets $\{\rho_k\}$ and $\{\Pi_\ell\}$ are chosen to be informationally complete, i.e., there is a one-to-one mapping between the $p_{k\ell}$s and the quantum channel $\mathcal{E}$. Once we have the estimate for $\{p_{k\ell}\}$, we apply the mapping to get an estimate for $\mathcal{E}$.

In our numerical examples below, we used a particular choice of $\{\rho_k\}$ and $\{\Pi_\ell\}$, namely, the product tetrahedron states and the product tetrahedron measurement. The set of $n$-qubit product tetrahedron states is the pure states $\psi_{k_1,k_2,\ldots,k_n}\equiv\psi_{k_1}\otimes\psi_{k_2}\ldots\otimes\psi_{k_n}$, with each $\psi_k$ a single-qubit state written in the Bloch-sphere representation as
\begin{equation}
\psi_k=\tfrac{1}{2}(\mathbbm{1}+\bm{a}_k\cdot\bm{\sigma}),\qquad k=1,2,3,\textrm{and }4,
\end{equation}
where the four $\bm{a}_k$s are three-dimensional unit vectors subtending a tetrahedron in the Bloch sphere \cite{Tetrahedronpaper}, and $\bm{\sigma}$ is the vector of Pauli operators. The product tetrahedron measurement is also defined in terms of these tetrahedron states: $\Pi_{\ell_1,\ell_2,\ldots, \ell_n}\equiv \Pi_{\ell_1}\otimes\Pi_{\ell_2}\ldots\otimes\Pi_{\ell_n}$, with each $\Pi_\ell\equiv \frac{1}{2}\psi_\ell$, for $\ell=1,2,3,$ and 4. The tetrahedron states and measurements form an informationally complete QPT scheme. For $n$ qubits, the scheme requires $4^{2n}$ different settings of input state and output measurement. Owing to the tetrahedron geometry, this scheme uses the minimal number of different input states and measurement outcomes in a symmetric way.

The data from the QPT experiment consist of a sequence of detector clicks, which we summarize into $\{n_{k\ell}\}$, with $n_{k\ell}$ being the number of clicks in the detector for $\Pi_\ell$ for input state $\rho_k$, where $k,\ell=1,2,\ldots,4^n$. The total number of uses of the channel is $N= \sum_{k\ell}n_{k\ell}$. To reconstruct the channel $\mathcal{E}$ from the data, we first do linear inversion, i.e., solving for $\mathcal{E}$ (using the Choi-state representation) by replacing, in the Born rule given by Eq.~\eqref{eq:BornRule}, $p_{k\ell}$ with the relative frequency $n_{k\ell}/N$. This does not guarantee that the $\mathcal{E}$ obtained from the linear inversion will be a valid, i.e., completely positive, channel. Complete positivity is then enforced by projecting the solution from linear inversion onto the nearest CPTP map using the algorithm of Ref.~\cite{CPTPalgorithm}, thereby giving us an estimate $\widehat \mathcal{E}$ of the channel. The minimum fidelity $F_{\min}$, our quantity of interest, is then estimated by numerically minimizing $F(\psi,\widehat{\mathcal{E}}(\psi))$ over $\psi$, with standard conjugate-gradient methods.

We want to obtain an estimate of $F_{\min}$ that is within $\epsilon$ of the true value. The accuracy of $F_{\min}$ is controlled by the accuracy of the estimate $\widehat{\mathcal{E}}$ --- a better estimate of $\mathcal{E}$, obtained with a larger $N$, will give a more accurate estimate of $F_{\min}$. \textit{A priori}, we do not know the $N$ needed to attain the specified accuracy on $F_{\min}$. That depends on the unknown channel $\mathcal{E}$. For a fair comparison with our GFM scheme, we need $N$ just large enough so that $F_{\min}$ attains the desired accuracy. We thus increase $N$ slowly and look for convergent behavior. See Fig.~\ref{fig:scheme} for the summary of the scheme.

We begin with a small number of uses of the channel $N_1$ to estimate $\mathcal{E}$ and then $F_{\min}$. In the next iteration, we double the number of uses by measuring a further $N_1$ uses, combining the obtained data with that from the previous iteration to obtain a second estimate of $F_{\min}$. In subsequent iterations, the amount of data used to estimate $F_{\min}$ is doubled each time by doubling the number of uses of the channel. We continue the iteration, obtaining (hopefully) more and more accurate estimates of $F_{\min}$ until both stopping criteria are met. The first requires that the error bar $\widehat\sigma^{(k)}$ ---estimated through bootstrapping of the obtained data---of $F_{\min}$ is smaller than $2\epsilon$. The second criterion demands that the change of $F_{\min}$ from the previous iteration (with half the number of channel uses) is smaller than some threshold $g_\mathrm{thres}$. A schematic of the stopping criteria is given in Fig.~\ref{fig:QPTStop}.

\begin{figure}[t]
\centering
\includegraphics[width=0.65\columnwidth]{DEOMGF-Fig3.pdf}
\caption{ \label{fig:QPTStop} Diagrammatic representation of the stopping criteria used in the QPT fidelity estimation. See main text for a description of the various quantities.}
\end{figure}

Note that despite the similarity to our GFM approach where the procedure was also terminated according to range and gradient criteria, the reason why a gradient criterion is needed here for QPT is rather different. The $F_{\min}$ estimation from QPT was observed in our numerical investigations to be biased, usually underestimating the $F_{\min}$ unless $N$ is large enough. Similar to how we chose our $g^{(1)}$ values for the GFM scheme, suitable values of $g_\mathrm{thres}$, different for different number of qubits and $\epsilon$ values, were chosen by running tests for 100 randomly chosen trial channels and setting $g_\mathrm{thres}$ such that the $F_{\min}$ value is within the desired accuracy for 95\% of the test channels.

\section{Performance of our scheme}\label{sec:perf}

How well does our GFM scheme perform, in terms of the number of uses of the channel, compared with the alternative using QPT fidelity estimation? We are interested, in particular, with the scaling of the resource cost with the number of qubits, for fixed target precision $\epsilon$. The resource cost of QPT is well known to scale poorly with the dimension; does our scheme do better in that respect?

Intuitively, one would expect our GFM scheme to outperform QPT fidelity estimation. Our scheme directly estimates the one quantity of interest, $F_{\min}$, while QPT first estimates the full process matrix, which contains all information about the channel, and then estimates $F_{\min}$, discarding the rest of the gathered information. Our minimalistic approach should thus win over QPT. Indeed, resource estimates suggest this, as we explain below. However, our numerical examples, which take into full account the complexities of both schemes, as well as the particular class of channels considered, tell a somewhat different story. Below, we first discuss resource estimates to give some indication of potential performance and then describe the numerical comparison obtained from simulations of two specific classes of noise channels.

\subsection{Resource estimates}\label{sec:theory}

We can gain some insights into the possible performance by considering known theoretical bounds on the various components of our scheme as well as on QPT. We begin with QPT. A general scaling law for the resource cost of quantum \emph{state} tomography was derived in Refs.~\cite{Kueng2014,SampleComplexity}: $O(dr^2/\epsilon'^2)$ number of copies are needed to estimate the density operator $\rho$ of a $d$-dimensional system to accuracy $\epsilon'$, as measured by the trace-distance deviation from the true state. $r$ is the rank of $\rho$. The channel-state duality via the Choi--Jamio{\l}kowski isomorphism \cite{choipaper,Jamio} allows us to apply this result directly to QPT. Specifically, a $d$-dimensional quantum channel corresponds to a $d^2$-dimensional Choi state. Further assuming that the Choi state is full rank (supposing we have no reason to assume otherwise), we see that the bound for the resource cost of QPT becomes $O(d^6/\epsilon'^2)$ uses of the quantum channel, to estimate the Choi state within $\epsilon'$ \mbox{(trace)} distance from that of the true channel.

We are, however, interested in the accuracy of $F_{\min}$, not in the Choi-state trace distance. Accurate reconstruction of the Choi state of the channel, of course, assures that the estimate of $F_{\min}$ will be accurate as well.
Specifically, we can show a relationship between the deviation in minimum fidelity and the Choi-state trace distance (see Appendix~\ref{sec:app1}),
\begin{equation}\label{eq:fidTr}
|F_{\min}'-F_{\min}|\leq 2d\Vert\rho_{\mathcal{E}'}-\rho_\mathcal{E}\Vert_\mathrm{tr},
\end{equation}
for the Choi states $\rho_{\mathcal{E}}$ and $\rho_{\mathcal{E}'}$ of two arbitrary channels $\mathcal{E}$ and $\mathcal{E}'$, respectively, where $F_{\min}$ and $F_{\min}'$ are the minimum fidelity of the channels $\mathcal{E}$ and $\mathcal{E}'$, respectively, and $\Vert M\Vert_\mathrm{tr}\equiv \frac{1}{2}\mathrm{tr}\bigl(\sqrt{M^\dagger M}\bigr)$ denotes the trace norm. The factor of $d$ in the first term on the right-hand side of Eq.~\eqref{eq:fidTr} is unavoidable; in fact, there exist (see Appendix~\ref{sec:app1}) one-parameter families of channels $\mathcal{E}$ and $\mathcal{E}'$ such that $|F_{\min}'-F_{\min}|= d\Vert\rho_{\mathcal{E}'}-\rho_\mathcal{E}\Vert_\mathrm{tr}$.

To attain accuracy $\epsilon$ for $F_{\min}$, i.e., $|\widehat{F_{\min}}-F_{\min}|\leq \epsilon$, Eq.~\eqref{eq:fidTr} indicates that it is sufficient that the trace-distance accuracy of the QPT scheme satisfies $\Vert\widehat{\rho_\mathcal{E}}-\rho_\mathcal{E}\Vert_\mathrm{tr} \leq \epsilon'$, where $\epsilon'=\epsilon/(2d)$ and $\widehat{\rho_\mathcal{E}}$ is the reconstructed Choi state from QPT. This yields a resource cost of $O(d^6/(\epsilon/d)^2)=O(d^8/\epsilon^2)$, i.e., $O(d^8)$ uses of the channel for fixed $\epsilon$, an altogether prohibitive scaling.

Next, we look at the resource estimates for our GFM scheme. Our two main ingredients, the DTFE scheme and the chosen GFM algorithm, directly determine the resource cost. The efficiency of the chosen GFM algorithm determines the number of function evaluations needed to arrive at the minimum point; the efficiency of the DTFE scheme determines how many uses of the channel are needed for each function evaluation to obtain an estimate of the required accuracy.
For the latter, Ref.~\cite{DF1} gives $O(d/\eta^2)$ for estimating the target fidelity to an accuracy $\eta$. We cannot, however, directly use this estimate. As explained earlier, in our GFM scheme, $\eta$ is not a fixed value, but is a quantity that is adjusted as the GFM iteration proceeds. We thus also need to consider possible dimensional dependence of $\eta$. In all our numerical examples, we observe that $\eta\sim 1/d$, for the $\eta$ values used in the course of the iterations, for $F_{\min}$ accuracy $\epsilon=0.01$. Using this empirical estimate, we thus have that each use of the DTFE subroutine requires $O(d^3)$ uses of the channel.

For a GFM algorithm, there are generally two aspects to consider when estimating the resource cost: (i) the number of function evaluations needed in the course of the descent to a local minimum; (ii) the number of repeats of the descent to obtain the global minimum with high confidence. In our numerical examples fixed at $\epsilon=0.01$, we observe that only a constant, $d$-independent number of repeats was needed to arrive at the global minimum with high confidence~\cite{FN}. We thus disregard aspect (ii) in our considerations of the resource cost. For aspect (i), the CMA-ES algorithm has been observed empirically \cite{CMAES3} to require $O(d)$ to $O(d^2)$ function evaluations in total, to attain a fixed accuracy of the extremum function value. Combining this with the fact that each function evaluation invokes the DTFE scheme that requires $O(d^3)$ uses of the channel, we see that our scheme requires $O(d^4)$ to $O(d^5)$ channel uses, a possible improvement over QPT fidelity estimation.

\subsection{Numerical examples}\label{Examples}
While resource estimates can provide initial clues to the performance of a scheme, a more accurate test comes from numerical simulations of the scheme. Numerical tests, in particular, are able to give indications of variations in performance for different classes of channels, an aspect often not easily captured in a theoretical analysis.
In this section, we present our numerical comparisons of the performance of our GFM scheme with that of QPT fidelity estimation to estimate the minimum fidelity for channels on one to five qubits.

\begin{figure}
\centering
\includegraphics[width=\columnwidth]{DEOMGF-Fig4.pdf}
\caption{\label{fig:numerics}The total number of channel uses for our GFM scheme and for QPT fidelity estimation for $n$ qubits and fixed target $F_{\min}$ accuracy of $\epsilon=0.01$. Two classes of random channels were used: (a) Class-HS, and (b) Class-PA (see main text). Each histogram is generated from 100 randomly chosen channels, with the same 100 channels used for both schemes. The medians are indicated by the arrows; the gray bars mark the height of 10 counts. The numbers on the right (top number for GFM, bottom for QPT) of the histograms give the number of channels, out of 100, for which the estimated $F_{\min}$ is within the desired accuracy. The stopping criteria are adjusted so that these numbers are around 95. The scaling of the number of channel uses with the dimension $d=2^n$ is observed, using the medians of the histograms, as (a) Class-HS: $O(d^{4.2})$ for QPT and $O(d^{2.7})$ for GFM; and (b) Class-PA: $O(d^{3.8})$ for QPT and $O(d^{3.4})$ for GFM.}
\end{figure}

The test channels used in our numerical simulations are drawn from two commonly encountered classes of random channels. It is possible to construct specific channels where one of the schemes outperforms the other. However, these examples are artificial and seldom relevant in practice. The first class, which we label as Class-HS, comprises channels chosen randomly according to the Hilbert--Schmidt measure on the space of quantum channels \cite{randomChannels}. To sample from this class, Haar-random unitary matrices $U$ are chosen on a $d^2$-dimensional Hilbert space, considered as the space of the system with an equal-dimensional ancilla A, and the channel is defined as $\mathcal{E}(\cdot)\equiv \mathrm{tr}_\mathrm{A}\bigl(U(\cdot\otimes\psi_\mathrm{A})U^\dagger\bigr)$, for $\psi_\mathrm{A}$ a pure state on A. This first class of channels is a well-studied class used in many discussions of generic properties of quantum channels.

The second class, which we label as Class-PA, imitates noise channels observed in many quantum experiments today: $n$-qubit weak noise channels composed of a random Pauli channel followed by a random amplitude-damping channel. Here, the Pauli channel is $\mathcal{E}_\mathrm{Pauli}(\cdot)\equiv\sum_{k=1}^{d^2}p_k W_k\rho W_k$, for $d\equiv 2^n$, with $W_k$ a tensor-product Pauli operator and $\{p_k\}$ a probability distribution. The random Pauli channel is generated by drawing numbers $q_2,q_3,\ldots, q_{d^2}$ uniformly from a $(d^2-2)$-simplex, and then setting $p_k=uq_k$, for $k=2,3,\ldots,d^2$, with $u$ uniformly randomly chosen from the range $0$ to $0.1$, and $p_1\equiv 1-\sum_{k=2}^{d^2}p_k$. The upper limit of $0.1$ on $u$ ensures that we have a weak (i.e., close to the identity) Pauli channel. The amplitude-damping channel is $\mathcal{E}_\mathrm{AD}(\cdot)\equiv \sum_{j=0}^{d-1} E_j (\cdot)E_j^\dagger$, with
\begin{align}
E_0&\equiv \sum_{j=0}^{d-1}\sqrt{1-\gamma_j}|\psi_j\rangle\langle\psi_j|,\nonumber\\
\textrm{and} \quad E_j&\equiv \sqrt{\gamma_j} |\psi_0\rangle\langle \psi_j|, \quad j=1,2,\ldots,d-1,
\end{align}
where $\gamma_0= 0$, $\gamma_1= 0.1$, $\gamma_2,\gamma_3,\ldots \gamma_{d-1}$ are uniformly randomly chosen from the range $0$ and $0.1$, and $\{\psi_j\}_{j=0}^{d-1}$ is a Haar-random basis for the system Hilbert space. The amplitude-damping channel models population decay of the $n$-qubit system towards some $|\psi_0\rangle$ state. Physically, $|\psi_0\rangle$ is typically the energetic ground state of the system, which may or may not be aligned with the chosen Pauli axes directions.

In each numerical experiment, we randomly (from the chosen class of channels) choose an $n$-qubit channel and estimate $F_{\min}$ using our GFM scheme and, separately, using QPT fidelity estimation. The target accuracy for $F_{\min}$ is set to $\epsilon=0.01$. The experiment is repeated 100 times (i.e., 100 different channels) for a fixed channel class and fixed $n$. The tuning parameters, namely, the parameters in the CMA-ES algorithm and the DTFE procedure, and those in the stopping criteria for the GFM scheme and the QPT fidelity estimation scheme are prechosen for each channel class and $n$ (see Sec.~\ref{sec:directEst}) and fixed throughout the 100 experiments. The histograms for the experiments are given in Fig.~\ref{fig:numerics}.

The numerical experiments tell a rather different story than the resource estimates of the previous section. The most striking feature of Fig.~\ref{fig:numerics} is that QPT fidelity estimation requires fewer uses of the channel than our GFM scheme, at least for up to the tested five-qubit situation. In terms of scaling with the dimension of the system, the number of uses of the channel for our GFM scheme is about $O(d^3)$, not far from the resource estimates; that for QPT is, however, closer to $O(d^4)$, rather than the $O(d^8)$ behavior of the resource estimates. Overall, the numerics suggest that our GFM scheme shows some advantage over QPT in terms of scaling with the size of the system, with a larger advantage for Class-HS than for Class-PA (see Fig.~\ref{fig:numerics}). However, when the system size is small, i.e., the case of practical interest in the near future, QPT outperforms our scheme in actual number of uses of the channel.

One possible reason behind the significant difference between the observed numerical scaling and the resource estimates on the number of channel uses is that the inequality relation in Eq.~\eqref{eq:fidTr} merely provides an upper bound on the $F_{\min}$ deviation. Given a true channel $\mathcal{E}$ and its estimate $\widehat\mathcal{E}$ from QPT, the deviation in minimum fidelity of $\mathcal{E}$ and $\widehat\mathcal{E}$ is typically much smaller than the upper bound indicated in Eq.~\eqref{eq:fidTr}. This is particularly true for the two classes of channels in our numerical examples, as we observed empirically. Therefore, the use of Eq.~\eqref{eq:fidTr} grossly overestimates how stringent we need to be in the trace-distance deviation to achieve a desired fidelity deviation. This highlights the importance of directly incorporating the desired figure of merit---the minimum fidelity in our case---into the QPT scheme, using it within the stopping criteria, rather than imposing a target accuracy on the deviation of the full process matrix.

Note that the spread in the number of channel uses visible in Fig.~\ref{fig:numerics} is not solely due to our random choice of channels. The fluctuations in the data (and hence the final estimated $F_{\min}$ values) for each run of both schemes also contribute to the observed spread. The histograms shown in Fig.~\ref{fig:numerics} show only one run per scheme per random channel. Additional numerical studies indicate that if each procedure were repeated 100 times on the \emph{same} channel, the spread in the number of channel uses would be roughly as large as what is seen in Fig.~\ref{fig:numerics}. Also, note that the gaps in the histograms for the QPT fidelity estimation in Fig.~\ref{fig:numerics} (blue bars) are due to our choices of $N_1$, the base number of channel uses, subsequently doubled in each round of the iterative procedure.

A further remark concerns the resource scaling of our GFM scheme for Class-HS. Observe in Fig.~\ref{fig:numerics}(a) that we need $O(d^{2.7})$ channel uses for our scheme, a more favorable scaling than the $O(d^{4.2})$ of QPT fidelity estimation. In fact, we suspect a further slowdown in the increase in channel uses with $d$ beyond five qubits for our scheme, gaining further advantage over QPT fidelity estimation. This is because, numerically, we observe that for a typical channel from Class-HS, the range of fidelity values over all pure states shrinks logarithmically as $d$ increases. This means that for large $d$ and fixed target accuracy for the estimation of $F_{\min}$, we only need to pick \emph{any} pure state $\psi$, estimate $f(\psi)$ using the DTFE scheme, and that is already close to the true $F_{\min}$ value, even without further minimization using the GFM algorithm. We thus expect the number of uses of the channel for our GFM scheme to be much reduced in that case. Note that this reduction in the fidelity range is not observed for Class-PA.

As the classical computational resource remains one of the limiting factors of tomography, it is also important to compare the two schemes in this aspect. The classical computational resource cost of the GFM scheme is estimated to be $O(d^5)$. For the QPT fidelity estimation scheme, the estimated classical computational resource cost is at least $O(d^6)$. See Appendix~\ref{sec:app2} for the discussion. In this aspect, the GFM scheme has better performance compared with the QPT fidelity estimation scheme.

\section{Conclusion}\label{sec:conc}

We have explored two different direct schemes for estimating the minimum gate fidelity. The GFM scheme presents an interesting hybrid application of numerical GFM algorithms together with the experimental procedure of direct target fidelity estimation. We compared this to the alternative approach of QPT fidelity estimation, a key difference of which from standard QPT procedures is the direct incorporation of the quantity of interest as the stopping criterion in the iteration. Resource estimates suggest an extremely high cost for QPT fidelity estimation, with a scaling for $O(d^8)$ gate uses, compared with $O(d^4)$ for the GFM approach. This large difference in performance, however, was not seen for the numerical tests carried out on specific classes of noise channels, with both schemes showing a gate-use scaling closer to $O(d^4)$. This reminds us of the need to examine the performance of every procedure for the particular context at hand, in addition to general results that apply in all situations.

In practice, one could consider using either scheme. QPT fidelity estimation has the advantage that the QPT experiment can be easier to perform, using familiar measurement setups. Its performance, however, appears limited to $O(d^4)$ gate uses (for our specific numerical examples) to achieve an estimate of the minimum fidelity of a specified accuracy. On the other hand, while the current performance in our numerical examples is also $O(d^4)$, the GFM approach presents a potential for future improvement with an increase in the efficiency of the numerical GFM scheme used. The GFM approach is also more efficient in terms of classical computational resource cost. The downside, though, is that the GFM scheme, with the need for preparing arbitrary input states suggested by the GFM iteration,  can be more difficult to implement experimentally. That may, nevertheless, be well worth the effort if the number of gate uses can be significantly lowered for large $d$.

\vspace*{\fill}

\medskip
\begin{acknowledgements}
This work is supported by the Ministry of Education, Singapore (through Grant No. MOE2016-T2-1-130). The Centre for Quantum Technologies is a Research Centre of Excellence funded by the Ministry of Education and the National Research Foundation of Singapore.
\end{acknowledgements}

\vspace*{\fill}


\section{}

\section{Introduction}

The coronavirus SARS-CoV-2, which causes the disease COVID-19, has already changed the lives of billions of people globally. Many people have been ordered to stay in their homes, and economies worldwide have largely come to a halt. Thousands have died and at least several million have been infected. In the United States, social distancing measures began being implemented in mid-March, as states such as Washington, California, and New York saw a sharp rise in the number of hospitalizations attributable to the virus. Government response has been hindered by insufficient supply of materials needed for testing, and by the large proportion of infected individuals who are asymptomatic and therefore are unlikely to seek testing even if it is available. This has made it very difficult to know the true size of the infected population. Because effective surveillance has not been possible, policymakers have instead turned to social distancing policies as the only available tool to slow the spread of the virus.

Here, we seek to address two key questions: (1) how many people are actually infected or have ever been infected with SARS-CoV-2; and (2) Are the social distancing measures currently in place effective at suppressing the virus -- that is, can they be expected to lead to a decrease in the number of infections. The main challenge in performing such an analysis is data availability and quality. Ideally, we could use data on the number of confirmed cases to understand the prevalence of the disease and assess policy measures. However, for the reasons outlined above, case counts are wildly unreliable and mostly reflect the availability of testing and proportion of individuals who are asymptomatic rather than the actual size of the infected population. A better data source would be hospitalizations, but even these data are not universally available across different U.S. states, and in some states, such as California, the policy is that even hospitalized patients should not be tested unless the outcome of the test would change the clinical management of the disease. Differences in testing policies for hospitalized patients across states make these data unreliable. Thus, here we build a model of the spread of COVID-19 fit only on data that has reasonable chance of being accurate:  data on deaths attributable to COVID-19. Even these data are imperfect. For example, New York City recently revised its death count upwards by roughly 30 percent by including victims who were COVID-19 suspected but not confirmed. However, patients with severe disease are the most likely to be tested, and therefore the death count data are likely the most reliable data source we have that contains any information about the rate of infection.

Much of the epidemiological literature focuses on fitting or calibrating ordinary differential equation (ODE) models to observed data \cite{hethcote2000mathematics}, and this underpins our approach as well. Ideally, these data would be counts of infected and recovered individuals over time. In order to fit such a model to death data, modifications to the standard approach must be made. The observed deaths need to be linked to the underlying state variables of the ODE model via a sensible likelihood. We do this by means of a distribution for the time from infection to death, and an assumption about the infection fatality rate (IFR) -- the proportion of infected individuals who will eventually die of COVID-19. We base these values on clinical data from Wuhan, China, the Diamond Princess cruise ship, and South Korea. Using only death data, these parameters and the parameters of the underlying ODE model are not separately identifiable. For this reason, the time to death distribution and IFR are assumptions of our model. There exist external estimates of the time to death distribution based on high quality data. For the IFR, precisely because of the difficulty outlined above in establishing the true number of infections, there remains considerable uncertainty about its true value.  Because the IFR is such an important assumption of our model, we perform analysis for several values of the IFR. Conditional on these assumptions and a sampling model for deaths given infections, we fit the parameters of an underlying ODE model of epidemic dynamics to the observed death data.

The remainder of this paper is organized as follows. In section 2 we review related work. In section 3 we introduce our model and explain how the model can be used to infer the number of infections and the effect of social distancing measures. In section 4 we describe an MCMC algorithm for fitting our model. In section 5 we give results. Section 6 concludes and offers some insight on policy options.

\section{Related Work}

Several previous studies have attempted to model the dynamics of the pandemic in various geographical locations, with varying goals. Several of these studies have provided some analysis or estimate of the amount by which confirmed cases undercount the true number of infections. For example, \citet{li2020substantial} propose that in the first month of the epidemic in China, 82-90 percent of infections were undocumented. \citet{riou2020adjusted} use a SEIR model, an epidemiological ODE model, and calibrate their model to the time series of reported deaths and reported infections. By modeling the underreporting of symptomatic cases, and by assuming that approximately half of infections lead to symptomatic cases, they estimate the infected population in Hubei, finding that approximately 30\% of infections were documented. \citet{perkins2020estimating} estimate directly that in the US, more than 90\% of infections have been undocumented by tests using Chinese data and initial reports in the US.

\citet{ferguson2020impact} model the effect of transmission between susceptible and infectious individuals using a microsimulation model built on synthetic populations designed to mimic the populations of the United Kingdom and United States. They assume a fixed time-to-onset and a range of $R_0$ values from 2.0-2.6, and they assume symptomatic cases to be 50\% more infectious than asymptomatic cases. Having calibrated their model to the cumulative number of deaths, they estimate deaths and hospital loads under different non-pharmaceutical interventions (NPI) involving social distancing and isolation. This analysis was arguably the most influential in triggering the adoption of social distancing policies in the United States, as it indicated that without social distancing measures in place, there would likely be around 2.2 million deaths, and hospitals nationwide would be completely overwhelmed. They also predicted that social distancing measures could reduce the number of deaths substantially and prevent exceedance of hospital resources, but only if they were dynamically turned ``on and off'' by triggering mechanisms based on the current number of COVID-19 patients. The report envisioned social distancing remaining in place for roughly 18 months until a vaccine is likely to be available.

The \href{https://penn-chime.phl.io/}{CHIME app} \cite{weissman2020locally} is an online tool created by researchers at the University of Pennsylvania to help hospitals anticipate the number of incoming COVID-19 patients and their needs. The CHIME model uses the current number of COVID-19 hospitalizations to ``back out'' the total number of cases based on a user provided hospitalization rate conditional on infection. Similar to our work, their model does not rely on the case counts. They make forward projections for the number of hospital admissions, ICU admissions, and ventilators needed over the coming weeks. They allow the user to specify the parameters of their underlying epidemiological model as inputs in terms of the doubling time for the infected population.

The model given in \citet{murrayteam2020forecast} has a goal similar to CHIME (hospital use planning). This model (the IHME model) takes a different approach by fitting parametric curves to observed cumulative death rates. They use a hierarchical model on the parameters of the parametric curve. The model essentially projects that the future course of the cumulative death rate curve in the United States will follow a path similar to that observed in other locations that are farther along in the course of their epidemics. There is no underlying model of epidemic dynamics, but there is a sampling model for death rates, which allows them to give confidence intervals. The \emph{New York Times} \href{https://www.nytimes.com/interactive/2020/03/25/opinion/coronavirus-trump-reopen-america.html}{online tool} allows the user to specify inputs to understand how those inputs affect likely infections, hospital loads, and deaths; infections are a side-effect of the rest of the model.

The recent work of \citet{flaxman2020estimating} takes a similar approach to ours to evaluate the effect of social distancing orders across countries in Europe. They define a likelihood for death-only data and incorporate mitigation efforts in the model. They build a hierarchical model to estimate the number of infections and effects of mitigation across countries in Europe. This work is the most closely related to ours, and the overall approach is quite similar, but their analysis is done for Europe whereas here we consider only the United States. \citet{lewnard2020incidence} also evaluates the effect of mitigation efforts. They observe reductions in estimates of the effective reproduction number for patients in three hospital systems in Northern California, Southern California, and Washington State as a consequence of the implementation of non-pharmaceutical interventions, like social distancing.

\section{Model}

Let $\nu_t$ be the number of \emph{new} infections on day $t$ of the epidemic, and let $p$ denote the infection fatality rate, i.e. the probability of death given infection. We denote the day of the first infection by $T_0$. Let $\theta = \{\theta_s : s = 0, 1, ...., m\}$ be the set of probabilities defining the discrete time-to-death distribution, where $\theta_s$ denote the probability that, for those who die, death from COVID-19 occurs $s$ days after the initial infection.  Let $X(t,t')$ denote the number of individuals newly infected on day $t$ who die on day $t'$. Our death model is
\begin{equs}
X(t,t') \mid p, \theta \sim \Pois(p \nu_t \theta_{(t'-t)}).
\end{equs}
The observed deaths on day $r$ are thus given by
\begin{equs}
D(r) = \sum_{t=1}^r X(t,r),
\end{equs}
the sum over all previous days of the number of individuals infected on that day who went on to die on day $r$. This has marginal distribution
\begin{equs} \label{eq:Ddist}
D(r) \sim \Pois\left(p \sum_{t=1}^r \nu_t \theta_{(r-t)}\right).
\end{equs}

This equation defines our likelihood. The use of a Poisson distribution in specifying our model may seem unnatural compared to the specification $X(t, t') \mid p, \theta \sim \Binom(\nu_t, p\theta_{(t-t')})$. The Poisson specification allows $\nu_t$ to take real values as opposed to integer values, as would be required for a Binomial distribution. This allows us to use simpler, deterministic models for the underlying epidemiological curves defining the $\nu_t$s, which simplifies computation. Fortunately in cases where $p\theta_s$ is small and $\nu$ is large (precisely the situation in which we find ourselves after the very early days of the epidemic), the $\Pois(p\theta_s\nu)$ distribution is a good approximation to $\Binom(\nu, p\theta_s)$.

The observed number of deaths $D$ are linked to a compartmental epidemiological model via the total number of newly infected individuals on day $t$, $\nu_t$. We use an SIR model of epidemic dynamics, with state evolution given by the following ODE
\begin{equs} \label{eq:sir}
\frac{\partial {S}_t}{\partial t} &= \begin{cases} -\beta S_t I_t N^{-1} & t < T_1  \\ -\phi \beta S_t I_t N^{-1} & t \ge T_1 \end{cases} \\
\frac{\partial I_t}{\partial t} &= \begin{cases} \beta S_t I_t N^{-1} - \gamma I_t & t < T_1 \\ \phi \beta S_t I_t N^{-1} - \gamma I_t & t \ge T_1 \end{cases} \\
\frac{\partial R_t}{\partial t} &= \gamma I_t,
\end{equs}
where $S_t$ is the number of susceptible individuals at time $t$, $I_t$ is the number of infected individuals at time $t$, $N$ is the total population size, $R_t$ is the number of removed (recovered or deceased) individuals at time $t$, and $T_1$ is the time at which social distancing went into effect. This set of equations defines a Susceptible-Infected-Removed (SIR) model. Prior to mitigation policies, it is a standard SIR model with parameters $\{\beta, \gamma\}$; post-mitigation, it takes  parameters $\{\phi \beta, \gamma\}$. 1-$\phi$ represents the percent reduction in the rate of infection that occurred following implementation of mitigation policies.

In a SIR model with parameters $\{\beta, \gamma\}$, infected individuals are considered contagious, and the mean time between infection and removal (the end of the contagious period) is given by $\gamma^{-1}$.  $R_0$ is given by $\beta \gamma^{-1}$ -- this is the number people who are infected by a single infected individual in a population in which essentially everyone is susceptible. The $\nu_t$ variables on which our likelihood is conditioned can be extracted directly from this model. The number of new infections that occurred during day $t$ is simply the difference in the size of the $S$ compartment at time $t$ and $t+1$. This is because under this model, individuals may only leave the susceptible compartment by becoming infected, so the new number of infections that occurred in any time period is precisely the reduction in the number susceptible in that same period.

Figure \ref{fig:realization} shows one realization of our modified SIR model. The top panel shows one draw of the daily number of deaths. Notice that this is not a smooth curve. When fitting our model, these $D$s will be our data, and will be the only thing we observe directly. The other panels are what we will infer from the $D$s. The middle panel of the figure shows the $\nu$s, the number of daily new infections. The series of $D$s has roughly the same shape as the $\nu$s, though it lags it by about 25 days. This lag is due to the shape of the time from infection to death distribution $\theta$, which we describe below. The bottom panel gives the standard view of the SIR model, showing the total number of susceptible, infected, and removed at any given time point.

\begin{figure}[h]
\centering
\includegraphics[width=5in]{example_realization.png}
\caption{\label{fig:realization}One realization of our modified SIR models}
\end{figure}

Due to identifiability constraints, we do not estimate all of parameters of the model. Rather, we fix those for which there exists high quality information on reasonable values that is applicable in the context studied here. By and large, this includes those parameters pertaining more to the biology of the disease than to the social dynamics of its spread.\footnote{It is impossible to make a clean distinction between these two types of parameters. For example, if the rate of spread allows for the number of people requiring care to overwhelm existing healthcare infrastructures, the case fatality rate will increase, as those who require medical support to survive but cannot get it will die.} These include $\theta$ and $p$.

For $\theta$, we draw from two studies. \citet{zhou2020clinical} reports that in Wuhan, China, the time from symptom onset to death had a median time of 18.5 days with an interquartile range of 15 to 22 days. A Poisson-Gamma distribution with parameters $\{\alpha, \beta\} = \{27.75, 1.5\}$ matches the reported quantiles well. \citet{lauer2020incubation} estimate the incubation period of the disease, also using data from China. They report a median incubation time of 5.1 days with an estimated 97.5th percentile of 11.5 days and a 99th percentile of 14 days. The quantiles of a Poisson-Gamma distribution with parameters $\{\alpha, \beta\} = \{5.5, 1.1\}$ match these reported quantiles well. We calculate the distribution of the total time from infection to death by generating 100,000 samples from the described distribution of the incubation period and the described distribution of the time from symptom onset to death. The time to death is the sum of these two numbers. We truncate the maximum time to death from infection to be the 99th percentile of the generated samples. This results in a time to death distribution shown in Figure \ref{fig:time-to-death}.

\begin{figure}[h]
\centering
\includegraphics[width=3in]{theta.png}
\caption{\label{fig:time-to-death} Conditional on death occurring, the probability of death on each day following infection.}
\end{figure}

To set $p$ we rely on several external data points. \citet{russell2020estimating} use data from individuals on the Diamond Princess cruise ship to estimate an infection fatality rate of 1.2\% (95\% CI: 0.38\%-2.7\%) after adjusting for delays between infection confirmation and death. The ship was a closed population: we know who was on the ship and therefore who to test, so we have confidence in the denominator (all individuals infected with COVID-19 were identified because all people on the ship could be tested). It is possible that some negative tests were false negatives, which would lower the estimate of the infection fatality rate by increasing the denominator. The tests used for these individuals were based on quantitative polymerase chain reaction (qPCR), which has a much higher potential for false negatives than for false positives, so we do not expect false positives were a significant factor in these estimates.

Another external information source is reported case fatality rates in countries that have done aggressive testing and contact tracing, including testing asymptomatic individuals. In these places, the case numbers may approach the true number of infections, so case fatality information in these locations provide a reasonably tight upper bound on infection fatality rates. Here, we look to South Korea, which has done the most testing per capita. As of 25 March 2020 there were 9,137 confirmed positive cases and 126 deaths in South Korea.\footnote{\url{https://www.cdc.go.kr/board/board.es?mid=a30402000000&bid=0030}} This gives a crude case fatality rate of 1.4\%. \cite{JordanCFR} use  a method for adjusting crude case fatality rates to account for the fact that some open cases may still result in death. They estimate an adjusted case fatality rate in South Korea of about 0.29 \%.

Finally, \citep{li2020substantial} estimates the total size of the infected population during the epidemic in China . Their estimates of the total number of infected implies an adjusted infection fatality rate of about 0.4 percent (see the note between Li et al. (2020) co-author J. Shaman \href{https://statmodeling.stat.columbia.edu/2020/03/07/coronavirus-age-specific-fatality-ratio-estimated-using-stan/}{here}). \cite{wu2020estimating} estimate a symptomatic case fatality of 1.4\% in Wuhan, China. Because not all cases are symptomatic, this implies a lower IFR.

In each of these cases, differences in underlying demographics, comorbities, and other risk factors could limit the applicability of these estimates to the United States at large. However, such factors are difficult to adjust for, since data on risk factors is preliminary and limited. We believe that the comparison is still useful in providing a range of rough ballpark estimates. Thus, consistent with all of these estimates, and because the infection fatality rate is arguably the most important parameter in our model that cannot be learned from the data, we consider three different cases that cover the spectrum of proposed values: 0.5 percent, 1.0 percent, and 1.5 percent, though we focus on the 1.0 percent case in some figures.

The parameters we estimate in our model are $\gamma_r$, $\beta$, and $T_0$. We base our prior on the infectious period on previously published studies. \citet{ferguson2020impact}  give a mean generation time ($\gamma^{-1}$ in the SIR model) of 6.4 days; \citet{prem2020effect} assume an infectious period of 3 or 7 days in their simulations with an incubation period of about 6 days. While these aren't directly comparable due to differences in the model, we use these as rough guidance. We choose the prior $\gamma^{-1} \sim N_{[3.4, 9.4]}(6.4, 1.5)$, a truncated normal distribution corresponding to an infectious period that ranges from about 3 to 9 days. This is a fairly tight prior, consistent with other models for this disease. To set the prior on $\beta \mid \gamma$, we use the reported 95 percent confidence interval in \cite{li2020early} to establish rough upper and lower bounds on the $R_0$, and put $\beta \mid \gamma \sim N_{[\gamma, 4\gamma]}(2.5\gamma, 1.5\gamma)$, encompassing possible $R_0$ values of between 1 and 4 for all of the allowed values of $\gamma$.  We place a uniform prior on $T_0$ between 1 January and 20 February.
We also place a Uniform(0.01,0.99) prior on $\phi$ allowing it to encompass the possibility that mitigation policies almost completely stopped the spread of the disease ($\phi \approx 0$) and that the policy led to no reduction at all ($\phi \approx 1$). This prior rules out the possibility that the mitigation policies, in fact, increased $R_0$.

\section{Computation}\label{sec:computation}
We do computation by MCMC using the adaptive Metropolis algorithm (see \cite{haario2001adaptive}). The algorithm produces samples from the Bayesian posterior distribution of the parameters of our model, $\eta =\{ \beta, \gamma, T_0, \phi \}$. Specifically, we update $\eta$  by proposing a new set of parameters, $\eta^*$ from $N(\eta,\Sigma)$, with the time-inhomogeneous covariance $\Sigma$ computed using the method of \cite{haario2001adaptive}. We then calculate the  set of $\nu^*_t$s implied by $\eta^*$ by simulating from the SIR model.  We accept or reject the proposed $\eta^*$ using Metropolis-Hastings, with target density proportional to
\begin{equs} \label{eq:LogPost}
\ell(\nu(\eta), \eta) = \sum_{t=1}^T p(D(t) \mid \nu) + \log(\pi(\eta)),
\end{equs}
where $p(D(t) \mid \nu)$ is the Poisson probability mass function implied by \eqref{eq:Ddist}, and $\pi()$ represents the prior density. Notice that because $\nu$ is a deterministic function of $\eta$, both terms in \eqref{eq:LogPost} depend on the values of these parameters.

We run for 50,000 iterations (100,000 for the state of New York), begin adaptation after 10,000 iterations (20,000 for New York), and use 25,000 iterations of burn-in (50,000 for New York). Trace plots are shown in the appendix. R code for all of the analysis here is available at \url{https://github.com/jamesjohndrow/Covid19-modeling}.

\section{Results}
We fit our model to the daily number of deaths for several states in the United States: California, Florida, New York, and Washington. We use the state-level data compiled by the \textit{New York Times}.\footnote{\url{https://github.com/nytimes/covid-19-data}}. The last date in our data is 17 April 2020. We take the date of implementation of social distancing to be the first day on which restaurants and schools were both closed statewide, as recorded in a GitHub repository maintained by researchers at the University of Washington\footnote{\url{https://github.com/COVID19StatePolicy/SocialDistancing}}. This definition is of course somewhat arbitrary, but given limitations of the available data, it is difficult to conceive of a richer model that would allow the effects of social distancing to phase in gradually without making similarly strong assumptions about how much each type of measure is ``worth'' compared to the eventual statewide lockdowns that were implemented everywhere.

Figure \ref{fig:fit} shows model fit for each of the states for  $p = 0.01$. The gray bands show point-wise 95\% posterior predictive intervals; the black line shows the values in our data. Despite the low-dimensionality of the free parameters in our model, our model fits the data reasonably well. Occasionally, the true values stray outside of the intervals. This suggests that perhaps a negative binomial or other over-dispersed count distribution be considered for future models, though of course one expects to occasionally see values outside of a pointwise 95\% interval. We also believe that some of the extreme deviations from the trend of the model are recording artifacts. For example, the large spike immediately followed by a dip in deaths evident in the New York figure in early April does not appear in data provided by the city of New York \footnote{see the ``Daily Counts" chart here: \url{https://www1.nyc.gov/site/doh/covid/covid-19-data.page}}. In fact, during the time period in which the unexpectedly low count occurs at about 280 in our data, New York City alone recorded nearly 400 deaths. It may be the case that some of the deaths allocated to the previous day in our data, where we see an unexpectedly large count, actually occurred on the following day, where our data have a surprisingly low count. This smoother New York City data, however, does not include state-wide counts, making it unsuitable for our state-level analysis. For the sake of consistency and replicability, we use the same data source for all states.

\begin{figure}[h]
\includegraphics[width=3in]{data_fit_CA_0_01.png}
\includegraphics[width=3in]{data_fit_FL_0_01.png}
\includegraphics[width=3in]{data_fit_NY_0_01.png}
\includegraphics[width=3in]{data_fit_WA_0_01.png}
\caption{Posterior predictive analysis  for the trajectory of deaths in California, Florida, New York, and Washington. 95\% pointwise posterior predictive intervals shown in gray. Data is shown in black. \label{fig:fit}}
\end{figure}

Table \ref{tab:R0} shows posterior means and 95\% posterior credible intervals of pre-intervention $R_0$ for each state and for all considered values of $p$. These estimates are not especially sensitive to $p$ (good news since this is the parameter about which we are most uncertain), and all fall within the range between about 1.7 and 2.5.

\input{state_R0.tex} \label{tab:R0}

One quantity estimated by our model is the cumulative number of SARS-CoV-2 infections over time. Comparing these estimates with the reported number of confirmed cases by state allows us to estimate the extent by which confirmed cases of COVID-19 undercount the number of infections. We define the undercount at each time point to be our estimated number of cumulative infections as of that time point divided by the cumulative number of confirmed cases at that time. That it, the undercount is the multiplicative factor by which the recorded number of confirmed cases under-estimates the true number. Figure \ref{fig:undercount} shows the undercount for each state across time. Naturally, when $p$ is smaller, and therefore the number of deaths is a smaller proportion of the true number of infections, we estimate the undercount to be greater. Under the range of values for $p$ considered, we find that the undercount was somewhere between 10 and 60-fold across states as of early April. This is similar to the estimate of 50-85 fold in early April obtained by \citet{Bendavid2020.04.14.20062463}, which was specific to Santa Clara county, California. As testing has increased, the undercount factor has reduced to around five to 25. Table \ref{tab:undercount} shows the undercount and associated 95\% posterior credible intervals as of the last day that appears in our data, 17 April 2020.

\begin{figure}[h]
\includegraphics[width=5in]{state_undercount_daily.png}
\caption{ \label{fig:undercount} Estimates of the undercount factor for each state by day. }
\end{figure}

\input{state_undercount.tex}

Social distancing policies began coming into effect in parts of the United States around 15 March, and thus as of \today, enough time has now passed since the implementation of social distancing in at least some states that the effects will be visible in the deaths data. Recalling the time from infection to death distribution from Figure \ref{fig:time-to-death}, changes in the dynamics of new infections should begin to be visible in death data around two weeks following the onset of the policy change, and should be mostly visible by around four weeks. This makes right now an appropriate time to attempt to estimate the effect of social distancing on infection dynamics. In the SIR model, the number of infected individuals will grow when the rate of new infections is higher than the rate of removal and will decrease otherwise. In other words, the rates are equal when the time derivative of $I_t$ is zero, i.e.
\begin{equs}
\phi \beta S_t I_t N^{-1} &= \gamma I_t \\
\frac{\phi \beta S_t}{N \gamma} &= 1 \equiv R_t
\end{equs}
We thus focus on estimating the quantity $R_t$ as of \today. When $R_t>1$, the number of infected individuals is still growing, while if $R_t < 1$, it is declining and the virus is being suppressed. It is important to note that, because this quantity depends on the current state of $S_t$, which is decreasing in time, it is possible for the current measures to fail to suppress the virus today yet be sufficient to suppress the virus at some later time because $S_t$ will have declined.

Table \ref{tab:RT1} shows estimates of the posterior mean of $R_{T_1}$ along with 95 percent posterior credible intervals. Results are reported for four states and for three values of $p$: 0.005, 0.01, and 0.015. It is clear that this quantity is not very sensitive to the assumption about $p$, which is encouraging since this is the major assumption of our model. Based on these results, the current measures are probably sufficient to suppress the virus in New York and Washington, possibly sufficient in Florida, and probably insufficient in California.

\begin{footnotesize}
\input{state_RT1.tex}
\end{footnotesize}

We  complement our estimates of $R_t$ with short term forecasts for each state. We do not produce longer term forecasts because the significant uncertainty regarding the future course of policy regarding social distancing in each state makes such forecasts highly uncertain and to some extent misleading; however, we do provide some tentative analysis about whether whether the current measures have been sufficient to suppress the spread of the virus and by how much. Three week projections for New York, California, Washington, and Florida for $p= 0.01$ are shown in Figure \ref{fig:forecasts}. Gray areas indicate pointwise 95\% posterior predictive intervals. Consistent with the posterior estimates of $R_{T_1}$, deaths in New York and Washington are projected to decline from recent peaks. In California, a slow increase in the number of deaths is predicted by our model. Florida appears likely to see growth in the number of cases, though this growth may be very slow, and may appear more like a plateau than obvious growth.

\begin{figure}[h]
\includegraphics[width=3in]{CA_Death_forecast.png}
\includegraphics[width=3in]{FL_Death_forecast.png}
\includegraphics[width=3in]{NY_Death_forecast.png}
\includegraphics[width=3in]{WA_Death_forecast.png}
\caption{History with three week projections based on posterior mean (dark line), with posterior 95 percent pointwise credible intervals (gray region), based on $p=0.01$. \label{fig:forecasts}}
\end{figure}

\section{Conclusion}

We have built a model for the transmission of SARS-CoV-2 using only information which has a reasonable chance of being measured correctly:  the observed number of daily deaths, timing of containment measures, and information on the clinical progression of the disease. Even having endeavored to build a model around the best available information, it is likely that some of this information is imperfect. For example, anecdotal reports from New York and elsewhere suggest that not all deaths from COVID-19 are being accurately reported, so even deaths attributable to COVID-19 are likely not perfectly measured. However, we believe death data to be more accurate than either cases or hospitalizations, which are not even available for some states. Information on the clinical progression of the disease, including the infection fatality rate, will likely improve over time, as more data on this emerging disease becomes available and is analyzed. Because of the importance of the infection fatality rate in our model and the substantial uncertainty about its value, we consider three different scenarios which roughly cover the range of possibilities that have been proposed. Analysis for $p = 0.01$ is shown here; detailed analysis of alternative values of $p=0.005$ and $p=0.015$ are given in supplementary material.

In contrast to models that use a wider range of information that may be less precisely measured, we believe the main strengths of our approach are that it minimizes assumptions and prioritizes simplicity, interpretability, and parameter identifiability. Our model also has a proper likelihood and is fit to data in the Bayesian paradigm, rather than relying on ad hoc calibration of model parameters to produce trajectories resembling the observed data. This allows us to formally account for uncertainty in all of our estimates via reporting of posterior credible intervals. Finally, our model is underpinned by an SIR model of infection dynamics to link observed deaths to the underlying unobserved infections. It would not be difficult with our approach to substitute some other model in place of the SIR. In particular, any other compartmental model, such as a Susceptible-Exposed-Infected-Removed (SEIR) model could be used, or the model could be elaborated to incorporate more change points in the parameters of the compartmental model to account for fine grained policy analysis. However, the existing data does not seem adequate to burden the model with so many parameters at this time.

Using our model, we estimate that official case counts substantially undercount of the number of infections.  This is not a surprise. Despite recent increases in testing capabilities, our analysis suggests that we continue to undercount the number of infections by a factor of around five to 25. This estimate is, of course, highly dependent on our assumption about $p$. Our model also suggests wide variability in the initial transmissibility of the disease prior to intervention policies as well as wide variability in the effects of those policies. For example, while our model suggests that New York's epidemic has peaked and will continue to see a decrease in the number of deaths, it indicates that California may not yet have peaked, and in Florida a slow increase in the number of deaths seems likely. All of these projections are conditioned on the current policies staying in place and their effects on transmission remaining constant. Tightening or relaxing mitigation policies, or ``quarantine fatigue'' leading to a decrease in compliance by the public, would impact the trajectory of the spread, and this is not accounted for in our model. Finally, we are just entering the time period where the effects of mitigation are becoming visible in the death data. That leaves us effectively with little data to estimate the impact of the mitigation efforts. This makes our estimates sensitive to perturbations due to recording errors or delays (perhaps caused by holiday weekends). As time progresses, the effect of mitigation policies will be observable in a larger proportion of the death data, improving the reliability of the estimates. We anticipate updating the analysis as more data become available, and continually working to enrich the methodology as data quality and availability improves sufficiently to support more complex modeling.

\section{Acknowledgements}

The authors thank Dr Megan Price and Tarak Shah for their comments and suggestions.  We thank Alexander D'Amour for finding a mis-stated equation in an earlier version. This work was supported by grants to the Human Rights Data Analysis Group by the John D. and Catherine T. MacArthur Foundation and the Oak Foundation.