Query-optimal estimation of unitary channels in diamond distance

TL;DR

This paper presents an optimal algorithm for estimating unknown unitary quantum channels with minimal applications, improving efficiency over previous methods and establishing matching lower bounds.

Contribution

The authors develop a query-optimal algorithm for unitary channel estimation in diamond norm, achieving Heisenberg scaling with minimal resource usage.

Findings

Achieves $ extsf{d}^2/ extsf{ε}$ applications for $ extsf{d}$-dimensional unitaries.

Introduces a bootstrap technique to improve estimation accuracy.

Proves a matching lower bound confirming optimality.

Abstract

We consider process tomography for unitary quantum channels. Given access to an unknown unitary channel acting on a $\textsf{d}$-dimensional qudit, we aim to output a classical description of a unitary that is $\varepsilon$-close to the unknown unitary in diamond norm. We design an algorithm achieving error $\varepsilon$ using $O(\textsf{d}^2/\varepsilon)$ applications of the unknown channel and only one qudit. This improves over prior results, which use $O(\textsf{d}^3/\varepsilon^2)$ [via standard process tomography] or $O(\textsf{d}^{2.5}/\varepsilon)$ [Yang, Renner, and Chiribella, PRL 2020] applications. To show this result, we introduce a simple technique to "bootstrap" an algorithm that can produce constant-error estimates to one that can produce $\varepsilon$-error estimates with the Heisenberg scaling. Finally, we prove a complementary lower bound showing that estimation requires $\Omega(\textsf{d}^2/\varepsilon)$ applications, even with access to the inverse or controlled versions of the unknown unitary. This shows that our algorithm has both optimal query complexity and optimal space complexity.

Equations

$$\mathcal{U}(X):\rho \mapsto X\rho X^{†}.$$

$$\underset{B}{\sup }\underset{{\rho }_{AB}}{\max }\{d(({\mathcal{E}}_1\otimes {\mathbb{I}}_B)({\rho }_{AB}),({\mathcal{E}}_2\otimes {\mathbb{I}}_B)({\rho }_{AB}))\},$$

$${\mathrm{dist}}_{\mathrm{phaseop}}(U,V)=\underset{\varphi \in \mathsf{U}(1)}{\min }\Vert \varphi U-V{\Vert }_{\mathrm{op}}.$$

$$\mathrm{dist}(U,V)=\underset{\text{paths}\gamma }{\inf }\mathrm{len}(\gamma ),\text{where}\mathrm{len}(\gamma )={\int }_0^1\Vert {\gamma }^{\prime }(t)\gamma (t{)}^{-1}{\Vert }_{\mathrm{op}}\mathrm{d}t={\int }_0^1\Vert {\gamma }^{\prime }(t){\Vert }_{\mathrm{op}}\mathrm{d}t.$$

$$\mathrm{dist}(U,\mathsf{U}(1)V)=\underset{\varphi \in \mathsf{U}(1)}{\inf }\mathrm{dist}(U,\varphi V)$$

$${\mathrm{dist}}_{\mathrm{phaseop}}(U,V)\le \mathrm{dist}(U,\mathsf{U}(1)V)\le \frac{\pi }{2}{\mathrm{dist}}_{\mathrm{phaseop}}(U,V)$$

$$\frac{1}{2}\Vert \mathcal{U}(U)-\mathcal{U}(V){\Vert }_{\diamond }\le {\mathrm{dist}}_{\mathrm{phaseop}}(U,V)\le \Vert \mathcal{U}(U)-\mathcal{U}(V){\Vert }_{\diamond }$$

$$\overline{F}({\mathcal{E}}_1,{\mathcal{E}}_2)=1-F(J({\mathcal{E}}_1),J({\mathcal{E}}_2)),$$

$$\overline{F}(\mathcal{U}(U),\mathcal{U}(V))=1-{\left|\frac{1}{\mathsf{d}}\mathrm{Tr}(U^{†}V)\right|}^2.$$

$$\sqrt{1-\sqrt{1-\overline{F}(\mathcal{U}(U),\mathcal{U}(V))}}=\frac{1}{\sqrt{2\mathsf{d}}}\underset{\varphi \in \mathsf{U}(1)}{\min }\Vert U-\varphi V{\Vert }_{\text{F}}$$

$$4\overline{F}(\mathcal{U}(U),\mathcal{U}(V))\le \Vert \mathcal{U}(U)-\mathcal{U}(V){\Vert }_{\diamond }^2\le 2\mathsf{d}\cdot \overline{F}(\mathcal{U}(U),\mathcal{U}(V)).$$

$$\underset{U\notin \mathsf{U}(1)V}{\sup }\frac{\overline{F}(\mathcal{U}(U),\mathcal{U}(V))}{\Vert \mathcal{U}(U)-\mathcal{U}(V){\Vert }_{\diamond }^2}=\frac{1}{4}\text{ if d is even, and}\underset{U\notin \mathsf{U}(1)V}{\sup }\frac{\Vert \mathcal{U}(U)-\mathcal{U}(V){\Vert }_{\diamond }^2}{\overline{F}(\mathcal{U}(U),\mathcal{U}(V))}=2\mathsf{d}.$$

$$E[\mathrm{error}(\mathcal{U}(Z),\mathcal{U}(\hat{\mathbf{Z}}))]<\epsilon .$$

$$\underset{\mathcal{C}}{\mathrm{Pr}}[\mathrm{error}(\mathcal{U}(Z),\mathcal{C})>\epsilon ]\le \eta .$$

$$\mathcal{M}(\mathbf{U})=\sum _X\mathrm{Pr}[\mathbf{U}=X]\mathcal{U}(X).$$

$$\mathrm{error}(\mathcal{U}(Z),\mathcal{M}(\hat{\mathbf{Z}})).$$

$$F(\mathcal{U}(Z),\mathcal{M}(\hat{\mathbf{Z}}))=⟨\Phi |Z^{†}\underset{\hat{\mathbf{Z}}}{E}\hat{\mathbf{Z}}\Phi {\hat{\mathbf{Z}}}^{†}Z|\Phi ⟩=\underset{\hat{\mathbf{Z}}}{E}⟨\Phi |Z^{†}\hat{\mathbf{Z}}\Phi {\hat{\mathbf{Z}}}^{†}Z|\Phi ⟩=\underset{\hat{\mathbf{Z}}}{E}F(\mathcal{U}(Z),\mathcal{U}(\hat{\mathbf{Z}}))$$

$$\overline{F}(\mathcal{U}(Z),\mathcal{M}(\hat{\mathbf{Z}}))=\Vert \mathcal{U}(Z)-\mathcal{M}(\hat{\mathbf{Z}}){\Vert }_{\diamond }.$$

$$\Vert \mathcal{U}(Z)-\mathcal{M}(\hat{\mathbf{Z}}){\Vert }_{\diamond }=\overline{F}(\mathcal{U}(Z),\mathcal{M}(\hat{\mathbf{Z}}))=\underset{\hat{\mathbf{Z}}}{E}\overline{F}(\mathcal{U}(Z),\mathcal{U}(\hat{\mathbf{Z}}))\le \frac{1}{4}\underset{\hat{\mathbf{Z}}}{E}[\Vert \mathcal{U}(Z)-\mathcal{U}(\hat{\mathbf{Z}}){\Vert }_{\diamond }^2].$$

$$|\mathbf{u}⟩=\mathbf{\varphi }\sqrt{1-\mathbf{\epsilon }}|z⟩+\sqrt{\mathbf{\epsilon }}|\mathbf{w}⟩$$

$$\exists \text{ a diagonal unitary Φ such that }\Vert Z\Phi -\mathbf{V}{\Vert }_{\mathrm{op}}\le \epsilon \le \frac{1}{8}.$$

$$\exists \text{ diagonal unitary Φ such that }\Vert Z\Phi -\mathbf{U}{\Vert }_{\mathrm{op}}\le \epsilon /2,$$

$$|{\mathbf{u}}_j⟩={\mathbf{\varphi }}_j\sqrt{1-{\mathbf{\epsilon }}_j}|z_j⟩+\sqrt{{\mathbf{\epsilon }}_j}|{\mathbf{w}}_j⟩$$

$${\bm{U}}-Z\Phi=Z\Phi{\bm{D}}+{\bm{W}}{\bm{E}},\quad\text{where}\quad{\bm{D}}\coloneqq\operatorname{diag}\Big{\lparen}\big{\lparen}\sqrt{1-{\bm{\varepsilon}}_{j}}-1\big{\rparen}_{j}\Big{\rparen},\ {\bm{E}}\coloneqq\operatorname{diag}(\sqrt{{\bm{\varepsilon}}}).$$

$$\Vert \mathbf{D}{\Vert }_{\mathrm{op}}\le \Vert \mathbf{E}{\Vert }_{\mathrm{op}}=\max \{\sqrt{{\mathbf{\epsilon }}_j}\}\le \sqrt{{\epsilon }_0}$$

$$\text{Claim:}\Vert \mathbf{W}{\Vert }_{\mathrm{op}}\le C\text{except with probability at most 1002.}$$

$$\Vert \mathbf{U}-Z\Phi {\Vert }_{\mathrm{op}}\le 1\cdot 1\cdot \sqrt{{\epsilon }_0}+C\cdot \sqrt{{\epsilon }_0}$$

$$|{\mathbf{y}}_j⟩:=\sqrt{{\mathbf{\delta }}_j}{\mathbf{\psi }}_j|z_j⟩+\sqrt{1-{\mathbf{\delta }}_j}|{\mathbf{w}}_j⟩$$

$$\Vert \mathbf{Y}{\Vert }_{\mathrm{op}}\le K\text{ except with probability at most 1001}.$$

$${\bm{Y}}=Z{\bm{\Delta}}_{0}+{\bm{W}}{\bm{\Delta}}_{1},\quad\text{where}\quad{\bm{\Delta}}_{0}\coloneqq\operatorname{diag}\Big{\lparen}\big{\lparen}\sqrt{{\bm{\delta}}_{j}}{\bm{\psi}}_{j}\big{\rparen}_{j}\Big{\rparen},\quad{\bm{\Delta}}_{1}\coloneqq\operatorname{diag}\Big{\lparen}\big{\lparen}\sqrt{1-{\bm{\delta}}_{j}}\big{\rparen}_{j}\Big{\rparen}.$$

$$\Vert \mathbf{W}{\Vert }_{\mathrm{op}}\le (\Vert \mathbf{Y}{\Vert }_{\mathrm{op}}+\Vert Z{\mathbf{\Delta }}_0{\Vert }_{\mathrm{op}})\cdot \Vert {\mathbf{\Delta }}_1^{-1}{\Vert }_{\mathrm{op}}\le (K+1)\cdot (1-{\max }_j\{{\mathbf{\delta }}_j\}{)}^{-1/2}$$

Full text

1. 1 Introduction

2. 1.1 Preliminaries: Distances for unitary (and other) channels

3. 1.2 Unitary estimation and related problems

4. 1.3 Prior work

5. 1.4 Techniques

6. 1.5 Discussion

7. 2 Base tomography using $O({\mathsf{d}}^{2}/\varepsilon^{2})$ applications

8. 2.1 Ancillary results for base tomography

9. 3 Bootstrap of precision to $O({\mathsf{d}}^{2}/\varepsilon)$ applications

10. 3.1 Key lemma: Geometry of unitary groups

11. 3.2 Bootstrap algorithm

12. 4 Lower bounds

13. 4.1 Bound for constant $\varepsilon$

14. 4.2 Bound for general $\varepsilon$

15. A Deferred proofs for Section 1.1

16. B Estimating eigenvalues

17. C Gate-efficient pure state tomography with optimal sample complexity

18. D Compressed proof of Proposition 4.2

1 Introduction

Quantum mechanics models the time evolution of a closed system as a unitary operator on the vector space of states. From the perspective of an experimentalist, this model serves as a mechanism for producing statistical predictions. This operational point of view prompts a fundamental question: when and how one can interact with a system with some unknown dynamics to learn its corresponding unitary operator? This question generally goes under the name of process tomography, and has been studied extensively in a variety of settings depending on chosen figures of merit; see our discussion in Section 1.3 below. The figures of merit may be one or more of the following: the accuracy of an estimate (measured with some metric on the space of quantum channels), the number of times one has to run the unknown dynamics, the complexity of initial states and that of measurements, the size of coherent Hilbert space required, and the complexity of classical processing.

In this paper, we design an algorithm with the following properties. Let $\lVert\cdot\rVert_{\diamond}$ denote the diamond norm, the completely bounded norm on the space of linear operators acting on density matrices. Let ${\mathsf{U}}({\mathsf{d}})$ denote the set of all ${\mathsf{d}}\times{\mathsf{d}}$ unitary matrices. If $X\in{\mathsf{U}}({\mathsf{d}})$ is any unitary operator, we denote by $\mathcal{U}(X)$ the associated unitary channel

[TABLE]

Theorem 1.1** (Upper bound).**

There is a quantum algorithm $\operatorname{\mathcal{A}}$ that, given black-box access to an unknown ${\mathsf{d}}$-dimensional unitary channel $Z\in{\mathsf{U}}({\mathsf{d}})$ and any $\varepsilon>0$, outputs a classical description of a unitary. This output is probabilistic and can be viewed as a ${\mathsf{U}}({\mathsf{d}})$-valued random variable $\widehat{{\bm{Z}}}$. The algorithm satisfies the following properties:

(Query complexity)* $\operatorname{\mathcal{A}}$ queries the black box $O({\mathsf{d}}^{2}/\varepsilon)$ times.* 2. 2.

(Space complexity)* $\operatorname{\mathcal{A}}$ only uses one qudit of dimension ${\mathsf{d}}$. Specifically, it only prepares states of the form $V_{2}(ZV_{1})^{p}V_{0}\ket{0}$, where the positive integer $p$ and unitaries $V_{0},V_{1},V_{2}\in{\mathsf{U}}({\mathsf{d}})$ are adaptively chosen, and measures them in the computational basis.* 3. 3.

(Gate complexity)* If the ${\mathsf{d}}$-dimensional qudit is embedded in ${\mathsf{n}}=\lceil\log_{2}{\mathsf{d}}\rceil$ qubits, then $\operatorname{\mathcal{A}}$ uses $\operatorname{poly}({\mathsf{d}},1/\varepsilon)$ one- and two-qubit quantum gates in total beyond the queries to the black box, and the classical time complexity is $\operatorname{poly}({\mathsf{d}},1/\varepsilon)$.* 4. 4.

(Closeness of output unitary)* The unitary output by $\operatorname{\mathcal{A}}$ is $\varepsilon$-close to the unknown $Z$ in diamond norm with high probability. Specifically, $\operatorname*{\mathbf{E}}[~{}\lVert\mathcal{U}(\widehat{{\bm{Z}}})-\mathcal{U}(Z)\rVert_{\diamond}^{2}~{}]\leq\varepsilon^{2}$, which implies $\Pr[~{}\lVert\mathcal{U}(\widehat{{\bm{Z}}})-\mathcal{U}(Z)\rVert_{\diamond}\leq 3\varepsilon~{}]\geq\tfrac{2}{3}$.* 5. 5.

(Closeness of output channel)* We can also view the output of $\operatorname{\mathcal{A}}$ as a channel. Let $\mathcal{M}$ denote the mixed unitary channel induced by running $\operatorname{\mathcal{A}}$ on $Z$ and then applying its output $\widehat{{\bm{Z}}}$. Then $\left\lVert\mathcal{M}-\mathcal{U}(Z)\right\rVert_{\diamond}\leq\varepsilon^{2}$.*

Complementary to these performance guarantees, we also show the following lower bound, which also holds against algorithms that can use $Z^{\dagger}$ and controlled versions of $Z$ and $Z^{\dagger}$.

Theorem 1.2** (Lower bound).**

Let $\operatorname{\mathcal{A}}$ be an algorithm that, for an unknown ${\mathsf{d}}$-dimensional unitary $Z\in{\mathsf{U}}({\mathsf{d}})$ accessible through black box oracles that implement $Z$, $Z^{\dagger}$, ${\mathrm{c}}Z=\ket{0}\bra{0}\otimes{\mathbb{I}}+\ket{1}\bra{1}\otimes Z$, and ${\mathrm{c}}Z^{\dagger}=\ket{0}\bra{0}\otimes{\mathbb{I}}+\ket{1}\bra{1}\otimes Z^{\dagger}$, can output a classical description of a unitary channel $\mathcal{U}(\widehat{{\bm{Z}}})$ such that $\lVert\mathcal{U}(\widehat{{\bm{Z}}})-\mathcal{U}(Z)\rVert_{\diamond}<\varepsilon<\frac{1}{8}$ with probability $\geq\tfrac{2}{3}$. Then $\operatorname{\mathcal{A}}$ must use $\Omega(\frac{{\mathsf{d}}^{2}}{\varepsilon})$ oracle queries.

Due to subtle differences in assumptions and objectives, we give comparison to the prior art after we properly define mathematical notions. In brief: standard process tomography for general channels [Leu00] uses one qudit but requires a $1/\varepsilon^{2}$ error dependence, which we improve quadratically.555We were not able to locate a $O({\mathsf{d}}^{2}/\varepsilon^{2})$-query process tomography algorithm for diamond-norm distance in the literature; we prove this result in Section 2. Most prior work takes the accuracy measure to be either the one in Item 5 of Theorem 1.1, referred to in this work as the problem of Approximate Storage-And-Retrieval, or the closeness of the output unitary in entanglement infidelity. Prior work of Yang, Renner, and Chiribella [YRC20] gives an algorithm that applies the unknown $Z$ in parallel to achieve the current best guarantees for these problems. In particular, it recovers Theorem 1.1.5 with the same query complexity. However, these problems are easier than that of unitary estimation to diamond-norm distance: norm conversion gives an algorithm for diamond-norm distance with $O({\mathsf{d}}^{2.5}/\varepsilon)$ applications. We improve on this by a factor of $O(\sqrt{{\mathsf{d}}})$, and improve the number of qudits used from $O({\mathsf{d}}^{2.5}/\varepsilon)$ to $1$.

Notation.

Throughout the paper, all $\log$s have base $e=\sum_{k\geq 0}\frac{1}{k!}\approx 2.718$, and $\exp(x)=e^{x}=\sum_{j=0}^{\infty}\frac{x^{j}}{j!}$ for any $x\in{\mathbb{R}}$. The dimension of a qudit is denoted by ${\mathsf{d}}\geq 1$, which we sometimes take to be ${\mathsf{n}}$ qubits, so that ${\mathsf{d}}=2^{\mathsf{n}}$. For matrices, $\lVert\cdot\rVert_{\mathrm{op}}$ is the largest singular value (also known as the operator norm or spectral norm), $\lVert\cdot\rVert_{1}$ is the trace norm (also known as the Schatten 1-norm), and $\lVert\cdot\rVert_{\textup{F}}$ is the Frobenius norm. For a vector $v\in{\mathbb{C}}^{\mathsf{d}}$, $\operatorname{diag}(v)\in{\mathbb{C}}^{{\mathsf{d}}\times{\mathsf{d}}}$ denotes the matrix with $v$ on the diagonal. We use ${\mathbb{I}}$ to denote the identity element of ${\mathsf{U}}({\mathsf{d}})$.

1.1 Preliminaries: Distances for unitary (and other) channels

We are mainly interested in distance measures between unitary channels. However, for comparison with prior work we need to consider distances between general channels. For more information on the distances discussed herein, the reader may consult, e.g., [GLN05, YF17]. We defer proofs to Appendix A.

A standard way to measure distances between channels $\mathcal{E}_{1},\mathcal{E}_{2}$ operating on a register “$A$” is to choose a distance $d({\cdot},{\cdot})$ on mixed quantum states and then consider

[TABLE]

where $B$ is an ancilla register and $\rho_{AB}$ is a mixed state on both registers. When $d$ is trace distance, this leads to the diamond norm:

Definition 1.3**.**

When the state metric in Equation 2 is trace distance, $d(\rho_{1},\rho_{2})={\tfrac{1}{2}}\lVert\rho_{1}-\rho_{2}\rVert_{1}$, we obtain the diamond-norm distance, denoted $\operatorname{dist}_{\diamond}(\mathcal{E}_{1},\mathcal{E}_{2})={\tfrac{1}{2}}\lVert\mathcal{E}_{1}-\mathcal{E}_{2}\rVert_{\diamond}$.

The diamond-norm distance is a natural choice for measuring worst-case error. It has direct operational meaning: if a channel sitting in a quantum circuit is replaced by another channel that is $\varepsilon$-different in diamond-norm, then the output distribution of the algorithm is at most $\varepsilon$-different from the original in total variation distance. Between unitary channels, the diamond-norm distance is equivalent (up to constants) to other familiar metrics on unitary matrices. We will switch between these metrics throughout the paper depending on which is convenient. One such metric is operator norm distance up to phase, which we use in Section 2.

Definition 1.4**.**

For unitaries $U,V\in{\mathsf{U}}({\mathsf{d}})$ we define

[TABLE]

We must take distance up to phase because a unitary channel only specifies a corresponding unitary up to global phase; for example, $U$ and $-U$ perform identical operations on quantum states. Another equivalent metric is the intrinsic metric on the projective Lie group ${\mathsf{PU}}({\mathsf{d}})={\mathsf{U}}({\mathsf{d}})/{{\mathsf{U}}(1)}$, which we use in Section 3.1.

Definition 1.5**.**

We consider the length of a path on the Lie group ${\mathsf{U}}({\mathsf{d}})$ with respect to the operator norm on the tangent space ${\mathfrak{u}}({\mathsf{d}})$. For unitaries $U,V\in{\mathsf{U}}({\mathsf{d}})$, we define $\operatorname{dist}(U,V)$ to be the length of a shortest smooth path $\gamma:[0,1]\to{\mathsf{U}}({\mathsf{d}})$ between them ($\gamma(0)=U,\gamma(1)=V$). As an equation,

[TABLE]

We further define the corresponding metric on the projective group ${\mathsf{PU}}({\mathsf{d}})$ via this metric.

[TABLE]

An alternative way to define $\operatorname{dist}(U,V)$ is that it is the unitarily invariant metric satisfying $\operatorname{dist}({\mathbb{I}},\operatorname{diag}(e^{\mathrm{i}\phi_{1}},\ldots,e^{\mathrm{i}\phi_{\mathsf{d}}}))=\max_{k\in[{\mathsf{d}}]}\lvert\phi_{k}\rvert$ for $\phi_{1},\ldots,\phi_{\mathsf{d}}\in[-\pi,\pi)$.666 Taking this definition as a starting point, one may find it involved to prove triangle inequality, which is trivial under the definition by the minimization over paths. We will later see that these two definitions coincide.

That is, $\operatorname{dist}(U,V)$ is the largest angle an eigenvalue of $U^{\dagger}V$ forms with $1$.

Proposition 1.6** (Equivalence of diamond norm, operator norm, intrinsic Lie metrics).**

For unitaries $U,V\in{\mathsf{U}}({\mathsf{d}})$, the metrics $\lVert\mathcal{U}(U)-\mathcal{U}(V)\rVert_{\diamond}$, $\operatorname{dist}_{\mathrm{phaseop}}(U,V)$, and $\operatorname{dist}(U,{\mathsf{U}}(1)V)$ are equivalent up to constants. Specifically, the following inequalities hold.

[TABLE]

Other work on unitary estimation and “Storage-and-Retrieval” problem, defined below, considers an accuracy measure that is not equivalent: the entanglement (in)fidelity. This is the accuracy measure that naturally arises when one estimates a channel by performing state tomography on its Choi state. Another way to recover this norm is to consider the channel norm induced by Bures distance (square-root of infidelity), i.e. taking $d$ to be Bures distance in Equation 2, but rather than maximizing over all $\rho_{AB}$ one simply fixes $\rho_{AB}$ to be the maximally entangled state.777When $\mathcal{E}_{1}$ is unitary, the entanglement infidelity is the same as the so-called average gate infidelity up to the small constant factor of ${\mathsf{d}}/({\mathsf{d}}+1)$ [HHH99]. Without this restriction of fixing $\rho_{AB}$, the resulting distance is an alternative distance studied in [YF17]. It is termed the minimum gate fidelity when $\mathcal{E}_{1}$ is unitary.

Definition 1.7**.**

The entanglement infidelity between ${\mathsf{d}}$-dimensional channels $\mathcal{E}_{1},\mathcal{E}_{2}$ is defined to be

[TABLE]

where $F$ is squared quantum state fidelity $F(\rho,\sigma)=(\lVert\sqrt{\rho}\sqrt{\sigma}\rVert_{1})^{2}$ and $J(\mathcal{E})$ is the Choi state $(\mathcal{E}_{A}\otimes{\mathbb{I}}_{B})(\tfrac{1}{{\mathsf{d}}}\sum_{i,j}\ket{ii}\bra{jj})$ of $\mathcal{E}$.

The notation $\overline{F}({\cdot},{\cdot})$ is not standard. In case $\mathcal{E}_{1}=\mathcal{U}(U)$ and $\mathcal{E}_{2}=\mathcal{U}(V)$ for $U,V\in{\mathsf{U}}({\mathsf{d}})$, we have a simpler formula

[TABLE]

Remark 1.8**.**

$\overline{F}(\cdot,\cdot)$* is not a metric on ${\mathsf{PU}}({\mathsf{d}})={\mathsf{U}}({\mathsf{d}})/({{\mathsf{U}}(1)}{\mathbb{I}})$, but it is related to one:*

[TABLE]

is a metric (which is roughly $\sqrt{\overline{F}/2}$ for small $\overline{F}$). The Frobenius norm on the right-hand side can be thought of as the two-norm distance between the Choi states of $U$ and $V$, treated as pure state vectors, not density matrices, minimized over global, unphysical phase factors.

An elementary argument shows that estimating a unitary channel to entanglement infidelity $\varepsilon$ only implies estimating it to $\sqrt{{\mathsf{d}}\varepsilon}$ error in diamond-norm distance. The factor of $\sqrt{{\mathsf{d}}}$ reflects that entanglement infidelity is an average-case accuracy measure, whereas diamond norm is a worst-case accuracy measure. For example, the entanglement infidelity between the multiply controlled NOT gate $C^{k}X$ (e.g., the usual CNOT if $k=1$ or Toffoli if $k=2$) and the identity is $O(2^{-k})$ for large $k$, but the diamond-norm distance between the two is $1$.

Proposition 1.9**.**

For all ${\mathsf{d}}$ and any unitaries $U,V\in{\mathsf{U}}({\mathsf{d}})$,

[TABLE]

The left-hand side is sharp for ${\mathsf{d}}$ even and the right-hand side is sharp for all ${\mathsf{d}}\geq 2$:

[TABLE]

1.2 Unitary estimation and related problems

We are interested in the following task:

Definition 1.10**.**

The estimation task for unitary channels is the following: One is given black-box access to a unitary channel by $Z\in{\mathsf{U}}({\mathsf{d}})$. After applying the unitary channel some $Q$ times, the algorithm should output a classical description888Throughout this paper, for ease of exposition we will be treating these descriptions as if they have infinite precision. To formalize our results when numbers can only be specified up to a machine precision $\varepsilon_{\mathrm{mach}}$, one could discretize the space of unit vectors and unitary matrices with $O(\varepsilon_{\mathrm{mach}})$-nets. Then, a classical description of a unitary matrix or a unit vector is a list of numbers defining a matrix or vector that need not exactly satisfy the desired constraints, but is $\varepsilon_{\mathrm{mach}}$-close to an element of the net. When we need to implement a unitary from its classical description, we implement this element. This incurs an error of $\operatorname{poly}(d,1/\varepsilon)\varepsilon_{\mathrm{mach}}$ throughout the argument, which can be made to be $O(\varepsilon)$ by incurring $O(\log(d/\varepsilon))$ overhead in classical time complexity, beyond the factors of $O(\log(1/\varepsilon_{\mathrm{mach}}))$ from incorporating finite precision arithmetic (which we ignore in our discussions of time complexity). This issue only affects constant factors in our proofs in a very minor way, and we trust the diligent reader to take note of the required changes. of a unitary $\widehat{{\bm{Z}}}\in{\mathsf{U}}({\mathsf{d}})$. The algorithm may be probabilistic and hence $\widehat{{\bm{Z}}}$ should be thought of as a random variable. The goal is to achieve

[TABLE]

Here, $\operatorname{error}({\cdot},{\cdot})$ is some accuracy measure between unitary channels.

Note here we are only using the error measure between two unitary channels.

Remark 1.11**.**

More generally, one can introduce another parameter $\eta$ and have the goal be that $\operatorname{error}(\mathcal{U}(Z),\mathcal{U}(\widehat{{\bm{Z}}}))\leq\varepsilon$ except with probability at most $\eta$. As long as $\operatorname{error}({\cdot},{\cdot})$ is a metric, this can be achieved with $O(Q\log(1/\eta))$ applications of $Z$ as follows: First, Markov’s inequality implies that $\Pr[\operatorname{error}(\mathcal{U}(Z),\mathcal{U}(\widehat{{\bm{Z}}}))>3\varepsilon]\leq 1/3$. Then, one can reduce to a general $\eta>0$ at the expense of an $O(\log(1/\eta))$ factor using the “median trick”; see Proposition 2.4.

In order to more carefully quantify the resources needed to solve the unitary estimation task, we generally consider hybrid classical/quantum algorithms. These mix classical computation and measurement-outcome-processing with quantum state preparation, applications of $Z$ and other quantum operations, and measurements. Besides the query complexity $Q$, there are some additional resources one wishes to minimize, one of which is the following.

Definition 1.12**.**

In the task of estimating $Z\in{\mathsf{U}}({\mathsf{d}})$, assume that ${\mathsf{d}}=2^{{\mathsf{n}}}$ for some number of qubits ${\mathsf{n}}$. The space overhead of an estimation algorithm is the total number of qubits beyond the minimum, ${\mathsf{n}}$, used in the course of the algorithm.

There are several tasks closely related to unitary estimation in the literature, going under the names of “universal programming”, “unitary cloning”, “reference frame transmission”, and “unitary learning”. For comparison with previous work, it suffices for us to discuss the problem commonly known as “learning” a unitary channel. However, since “learning” and “estimation” are often regarded as synonymous in casual speech, we will follow [SBZ19] and use the terminology “Storage-And-Retrieval”, to avoid confusion.

Definition 1.13**.**

The SAR task for unitary channels (Storage-And-Retrieval) involves developing two quantum algorithms [BCD*+*10]. The “storage” algorithm $\mathcal{S}$ is given black-box access to a unitary $Z\in{\mathsf{U}}({\mathsf{d}})$; after applying the unitary some $Q$ times, $\mathcal{S}$ should output a possibly mixed quantum state $\psi$. The “retrieval” algorithm $\mathcal{R}$ takes as input $\psi$ and implements a possibly nonunitary qudit channel $\bm{{\mathcal{C}}}$, which may be random with some distribution. The goal is to achieve the following:

[TABLE]

Note that here $\operatorname{error}({\cdot},{\cdot})$ is measuring the error between a unitary channel and a possibly nonunitary channel. The storage complexity of the algorithm is the number of qubits used for $\psi$.

Remark 1.14**.**

The case of minimizing $\eta$ while insisting on $\varepsilon=0$ is known as Probabilistic SAR and is studied in, e.g., [SBZ19]. In this case, the channel $\bm{{\mathcal{C}}}$ is necessarily unitary. The output of the storage algorithm in [SBZ19] is always pure. On the other hand, for comparison to prior work, we focus on the case of $\eta=0$, which may be termed Approximate SAR. This latter problem is the one commonly known as “learning” a unitary.

We can perform this task with a unitary estimation algorithm by simply taking the estimate’s classical description and then synthesizing and applying it. We provide some notation formalizing this.

Definition 1.15**.**

If ${\bm{U}}$ is a ${\mathsf{U}}({\mathsf{d}})$-valued classical random variable, we denote by $\mathcal{M}({\bm{U}})$ the mixed-unitary channel associated to ${\bm{U}}$

[TABLE]

Note that if ${\bm{U}}$ is a unitary-valued random variable, the notation $\mathcal{U}({\bm{U}})$ stands for a unitary-channel-valued random variable; that is, $\mathcal{U}({\bm{U}})$ is not the same thing as $\mathcal{M}({\bm{U}})$. We will now compare unitary estimation and SAR. An instructive example to keep in mind is the case where an algorithm $\operatorname{\mathcal{A}}$ takes a black box for the identity channel $Z={\mathbb{I}}$ and outputs the matrix $\widehat{{\bm{Z}}}$, which is $[\begin{smallmatrix}1\\ &e^{\mathrm{i}\varepsilon}\end{smallmatrix}]$ or $[\begin{smallmatrix}1\\ &e^{-\mathrm{i}\varepsilon}\end{smallmatrix}]$ with probability $\frac{1}{2}$. One can verify that $\operatorname*{\mathbf{E}}[\lVert\mathcal{U}(Z)-\mathcal{U}(\widehat{{\bm{Z}}})\rVert_{\diamond}]=\Theta(\varepsilon)$ and $\lVert\mathcal{U}(Z)-\mathcal{M}(\widehat{{\bm{Z}}})\rVert_{\diamond}=\Theta(\varepsilon^{2})$, so by the following remark, on this input $Z={\mathbb{I}}$ with diamond-norm distance as the error measure, $\operatorname{\mathcal{A}}$ performs unitary estimation to error $\Theta(\varepsilon)$ and Approximate SAR to error $\Theta(\varepsilon^{2})$.

Remark 1.16**.**

A unitary estimation algorithm $\operatorname{\mathcal{A}}$ may be turned into an SAR algorithm with the same $\varepsilon$ and $\eta$ parameters: the storage algorithm $\mathcal{S}$ is $\operatorname{\mathcal{A}}$, with $\psi$ being the mixed state of orthogonal quantum states each of which is classical and encodes an output $\widehat{{\bm{Z}}}$ of $\operatorname{\mathcal{A}}$; and the retrieval algorithm $\mathcal{R}$ simply synthesizes and applies $\bm{{\mathcal{C}}}=\mathcal{U}(\widehat{{\bm{Z}}})$ from its classical description. Here, the channel $\bm{{\mathcal{C}}}$ is unitary.

More subtly, we can also turn $\operatorname{\mathcal{A}}$ into an Approximate SAR algorithm as follows. We let the storage algorithm the same as before, but let the retrieval algorithm implement a mixed unitary ${\mathcal{C}}=\mathcal{M}(\widehat{{\bm{Z}}})$ with probability one. This reduction is identical to the previous one, with the difference being that we treat the randomness of $\widehat{{\bm{Z}}}$ as part of the channel being applied. The accuracy of this Approximate SAR is given by

[TABLE]

This error measure is different in general from the error measure in the unitary estimation algorithm’s promise, $\operatorname*{\mathbf{E}}_{\widehat{{\bm{Z}}}}\operatorname{error}(\mathcal{U}(Z),\mathcal{U}(\widehat{{\bm{Z}}}))\leq\varepsilon$.

If $\operatorname{error}$ is diamond-norm distance, then $\operatorname{error}(\mathcal{U}(Z),\mathcal{M}(\widehat{{\bm{Z}}}))\leq\operatorname*{\mathbf{E}}_{\widehat{{\bm{Z}}}}\operatorname{error}(\mathcal{U}(Z),\mathcal{U}(\widehat{{\bm{Z}}}))\leq\varepsilon$ by the convexity of the norm, so in this setting Approximate SAR reduces to unitary estimation. However, $\operatorname{error}(\mathcal{U}(Z),\mathcal{M}(\widehat{{\bm{Z}}}))$ can in fact be made much smaller, to $O(\varepsilon^{2})$; see Proposition 1.18. We can think about Approximate SAR as, roughly, the version of unitary estimation where we only need to be right “in expectation”. In particular, though the output of the unitary estimation algorithm $\mathcal{U}(\widehat{{\bm{Z}}})$ is explicitly given, the output of the converted Approximate SAR algorithm $\mathcal{M}(\widehat{{\bm{Z}}})$ is not, so we might not know what it is.

Remark 1.17**.**

(Continued from Remark 1.16) If $\operatorname{error}=\overline{F}$, the entanglement infidelity, the accuracy of a unitary estimation algorithm is equal to that of its converted Approximate SAR algorithm. (Note that $\overline{F}$ is not a metric.) To see this, recall that the squared quantum state fidelity between a pure state $\psi$ and a mixed state $\rho$ is linear in $\rho$: $\lVert\sqrt{\psi}\sqrt{\rho}\rVert_{1}^{2}=\bra{\psi}\rho\ket{\psi}$. Then, with a maximally entangled state $\Phi=\ket{\Phi}\bra{\Phi}$, we have

[TABLE]

where we omitted subsystem indices and identity tensor factors. Hence, $\overline{F}(\mathcal{U}(Z),\mathcal{M}(\widehat{{\bm{Z}}}))=\operatorname*{\mathbf{E}}\overline{F}(\mathcal{U}(Z),\mathcal{U}(\widehat{{\bm{Z}}}))$.

One can use this fact about entanglement infidelity to show that the variance of a unitary estimation algorithm bounds the error of its converted Approximate SAR algorithm. In some sense, this gives a quadratically better version of Remark 1.16, since $\varepsilon$-error unitary estimation can give $\varepsilon^{2}$-error Approximate SAR.

Proposition 1.18** ([L]emma 2).**

*Yang2020]

If a unitary estimation algorithm $\operatorname{\mathcal{A}}$ with an output distribution $\varphi(\widehat{{\bm{Z}}}|Z){\mathrm{d}}\widehat{{\bm{Z}}}$ is unitarily covariant, i.e., $\varphi(A\widehat{{\bm{Z}}}B|AZB)=\varphi(\widehat{{\bm{Z}}}|Z)$ for all $A,B\in{\mathsf{U}}({\mathsf{d}})$, then*

[TABLE]

Therefore, by Remarks 1.17 and 1.9,

[TABLE]

1.3 Prior work

As far as we are aware, ours is the first work to study unitary estimation with the more stringent diamond-norm distance as its accuracy measure. We survey the existing literature on related topics, and then compare it to our result.

Standard process tomography.

Standard quantum process tomography [NC12, Chapter 8.4.2] solves the task of general-channel estimation problem by preparing a basis of quantum states, passing them through the channel, and performing state tomography on the results. One has to take care to analyze how the error bounds for state tomography affect the error bounds for channel estimation; see, e.g., [LSS*+*20] estimation of minimum gate fidelity for bounds of the form $Q\leq\operatorname{poly}({\mathsf{d}})/\varepsilon^{2}$ for general-channel estimation with respect to minimum fidelity. For the special case of unitary channels one can get improved bounds, as one only needs to work with $O({\mathsf{d}})$ pure-state estimation tasks. Naively analyzing this strategy yields an $Q\leq O({\mathsf{d}}^{3}/\varepsilon^{2})$ sample upper bound, as an $\varepsilon$ error for each state estimate can compound to a $O(\sqrt{{\mathsf{d}}}\varepsilon)$ error in diamond-norm distance for the resulting channel. But as we show in Section 2, with care this method can be used to obtain $Q\leq O({\mathsf{d}}^{2}/\varepsilon^{2})$ query complexity for unitary estimation with respect to diamond-norm distance. An advantage of this method is that it uses zero space overhead.

Another well-known approach to quantum process tomography is the ancilla-assisted method dating back to Leung [Leu00]. Here one prepares the maximally entangled state in ${\mathsf{d}}^{2}$ dimensions, passes the first half through the channel, and uses state tomography. Again, if the channel is unitary one can use pure state tomography, and with this approach it is not hard to deduce that $Q\leq O({\mathsf{d}}^{2}/\delta)$ queries suffices to obtain unitary estimation with respect to entanglement infidelity $\delta$. Besides using $O(\log{\mathsf{d}})$ space overhead, this approach also only gives an $Q\leq O({\mathsf{d}}^{3}/\varepsilon^{2})$ query complexity bound for diamond-norm distance $\varepsilon$, via the relations Propositions 1.6 and 1.9. However, like with standard process tomography, with a tighter analysis, one may be able to give an $Q\leq O({\mathsf{d}}^{2}/\varepsilon^{2})$ query complexity bound.

Unitary estimation.

The specific task of unitary estimation was perhaps first studied by Acín, Jané, and Vidal [AJV01], where representation theory was used to determine an optimal algorithmic strategy with respect to entanglement infidelity under the assumption that the unknown $Z$ is applied in parallel to one half of a bipartite system. Asymptotic analysis of the error was not given, however. Later, Peres and Scudo [PS02] gave an alternate method establishing that $Q\leq O(1/\sqrt{\delta})$ queries suffice to obtain entanglement infidelity $\delta$ in the case ${\mathsf{d}}=2$. Then Bagan, Baig, and Muñoz-Tapia [BBMT04a] established the same scaling for the method from [AJV01]; and, they [BBMT04b] and Chiribella, D’Ariano, Perinotti, and Sacchi [CDPS04] did the same for a similar method that didn’t require entangled measurements on both parts of the bipartite system. See also independent work of Hayashi [Hay06]. Chiribella, D’Ariano, and Sacchi [CDS05] again showed a more general optimality result for entanglement infidelity, implying that $Q\geq\Omega(1/\sqrt{\delta})$ is a lower bound for this accuracy measure when ${\mathsf{d}}=2$; however, the asymptotic dependence for ${\mathsf{d}}>2$ was not established. Later, Kahn [Kah07] showed that the optimal scaling for general ${\mathsf{d}}$ is of the form $Q=\Theta(f({\mathsf{d}})/\sqrt{\delta})$ for some function $f({\mathsf{d}})$, but was unable to asymptotically analyze it.

Finally, in 2020, Yang, Renner, and Chiribella [YRC20] were able to analyze the optimal unitary estimation algorithm for entanglement infidelity and showed that it achieves error $\delta$ with $Q\leq O({\mathsf{d}}^{2}/\sqrt{\delta})$ queries. Moreover, they showed Proposition 1.18, and since their algorithm is unitarily covariant, this implies an Approximate SAR algorithm with $Q\leq O({\mathsf{d}}^{2}/\sqrt{\varepsilon})$ with respect to diamond-norm error $\varepsilon$. We remark that this optimal unitary estimation algorithm from [YRC20] applies the unknown unitary $Z$ in parallel and hence has space overhead on the order of the query complexity $\Theta({\mathsf{d}}^{2}\log({\mathsf{d}})/\sqrt{\varepsilon})$.

Comparison to [YRC20].

Our main result recovers the unitary estimation (with respect to entanglement infidelity) and Approximate SAR (with respect to diamond-norm distance) results of Yang, Renner, and Chiribella [YRC20]. First, by norm conversion Proposition 1.9, our algorithm with $Q\leq O({\mathsf{d}}^{2}/\varepsilon)$ achieves the guarantee $\operatorname*{\mathbf{E}}[~{}\overline{F}(\mathcal{U}(\widehat{{\bm{Z}}}),\mathcal{U}(Z))~{}]\leq\operatorname*{\mathbf{E}}[~{}\lVert\mathcal{U}(\widehat{{\bm{Z}}})-\mathcal{U}(Z)\rVert_{\diamond}^{2}~{}]\leq\varepsilon^{2}$, giving the unitary estimation guarantee. Second, since any algorithm can be made unitarily covariant, using the same reduction of Proposition 1.18, we achieve an Approximate SAR algorithm with $Q\leq O({\mathsf{d}}^{2}/\sqrt{\varepsilon})$ with respect to diamond-norm error $\varepsilon$.

Neither of these results from [YRC20] imply our result of unitary estimation to diamond norm error. Through norm conversion Proposition 1.9, a $O({\mathsf{d}}^{2}/\sqrt{\delta})$-query unitary estimation algorithm with respect to entanglement infidelity $\delta$ implies a $O({\mathsf{d}}^{2.5}/\varepsilon)$-query algorithm with respect to diamond-norm distance $\varepsilon$. This conversion is tight, even if the unitary estimation algorithm is unitarily covariant. Also, as discussed in Remark 1.16, the output of an approximate SAR algorithm is correct in expectation, but any individual output need not be close in diamond norm.

Finally, our algorithms have space overhead of zero, improving over [YRC20] and making our algorithm significantly closer to practical. Note that the main figure of merit in [YRC20], program cost, denotes the size of the output of the storage algorithm, and is different from the space referred to here, which is the space complexity of the storage algorithm (Definition 1.13).

Comparison to [vACGN23].

Van Apeldoorn, Cornelissen, Gilyén, and Nannicini shows that, given a unitary $Z\in{\mathsf{U}}({\mathsf{d}})$ implementing an unknown state $Z\ket{0}=\ket{z}\in\mathbb{C}^{\mathsf{d}}$, one can compute an estimate $\ket{\widehat{{\bm{z}}}}$ which is $\varepsilon$-close in Euclidean norm with $\geq 1-\delta$ probability [vACGN23, Theorem 23]. Their algorithm uses $O(\frac{{\mathsf{d}}}{\varepsilon}\log\frac{{\mathsf{d}}}{\delta})$ applications of the controlled black box unitary ${\mathrm{c}}Z$ and its inverse ${\mathrm{c}}Z^{\dagger}$, along with $O({\mathsf{d}}\log\frac{{\mathsf{d}}}{\varepsilon})$ qubits. This implies an algorithm for unitary estimation, by using this state tomography algorithm for standard process tomography (see Propositions 2.3 and 2.2) which uses $O(\frac{{\mathsf{d}}^{2}}{\varepsilon}\log\frac{{\mathsf{d}}}{\delta})$ applications of ${\mathrm{c}}Z$ and ${\mathrm{c}}Z^{\dagger}$ and $\Theta({\mathsf{d}}\log\frac{{\mathsf{d}}}{\varepsilon})$ space overhead. This algorithm is based on quantum singular value transformation, so the use of ${\mathrm{c}}Z$ and ${\mathrm{c}}Z^{\dagger}$ and the nonzero space overhead appear to be inherent limitations of their approach. Similarly, other standard primitives in quantum algorithms like amplitude estimation, which we might ordinarily reach for when aiming for a quadratic improvement in error, have similar limitations in terms of requiring stronger access to $Z$ or space overhead. Our algorithm is more direct and so does not lose anything in query complexity or space complexity, and only uses queries to $Z$, and not ${\mathrm{c}}Z$ or ${\mathrm{c}}Z^{\dagger}$. The gate complexity of both algorithms is similar.

Comparison of techniques to prior work.

Existing algorithms for estimating a parametrized class of gates with Heisenberg scaling, like metrology with a GHZ state [JWD*+*08] and robust phase estimation [KLY15], proceed by applying the black box many times in parallel or in series (respectively). Our algorithm uses the same principles (see the warmup in Fig. 1), but as discussed in the next section, naive generalizations of these approaches cannot work. We introduce a novel “shift to identity” step to avoid areas in the space of unitary channels where applying the black box in series fails; this adaptivity allows us to extend robust phase estimation, which can estimate a parameter on the complex unit circle, to work for elements in a much trickier space.

The process tomography protocol most similar to the algorithm presented in our work is that of gate set tomography [NGR*+*21], which can achieve Heisenberg scaling, $\Theta(1/\varepsilon)$, in regimes of practical interest by interleaving long sequences of gates. This is qualitatively similar to what our algorithm does, though we are not aware of work giving theoretical bounds on such protocols.

Prior work on lower bounds.

Lower bounds on channel estimation and discrimination are just as well-studied as upper bounds. Some prior work has tried to characterize the optimal strategy for various problems, such as showing they are sequential or parallel [DFY07, BMQ22]. As mentioned previously, though these results describe optimal strategies for these problems, it is unclear if they imply lower bounds better than $\Omega({\mathsf{d}}^{2}+\frac{1}{\varepsilon})$. Up to log factors, this lower bound is an immediate consequence of parameter counting (e.g. [NC12, Chapter 8.4.2], [BKD14], which one can make rigorous with a Holevo bound [CvDNT99]) and the optimality of Heisenberg scaling. However, these two arguments do not combine naturally to get the $\Omega({\mathsf{d}}^{2}/\varepsilon)$ bound we would expect.

To our knowledge, we give the first jointly optimal lower bound, in the sense that for any decaying function ${\mathsf{d}}\mapsto\varepsilon=\varphi({\mathsf{d}})$ our lower bound $\Omega({\mathsf{d}}^{2}/\varphi({\mathsf{d}}))$ matches the query complexity of our algorithm up to a universal constant. Our lower bound combines prior work on unitary channel discrimination [BMQ22] with a reduction technique from quantum query complexity [CGM*+*09, BCC*+*17]. Though quantum query complexity focuses on diagonal unitaries, we observe that roughly the same technique can be applied to this general unitary estimation setting with some modification.

1.4 Techniques

An algorithm with suboptimal $\varepsilon$ scaling.

A standard idea in process tomography is to estimate $Z$ by doing state tomography on $Z|0\rangle$, $Z|1\rangle$, and so on, and then somehow collate these estimates into a full estimate of $Z$. Since each state tomography requires $O({\mathsf{d}}/\varepsilon^{2})$ samples, this whole procedure can be done with $O({\mathsf{d}}^{2}/\varepsilon^{2})$ applications of $Z$ (Theorem 2.1). There are two minor issues to address when formalizing this. First, an $\varepsilon$ error in each state tomography could cascade to a $\varepsilon\sqrt{{\mathsf{d}}}$ error in the final estimate of $Z$. We address this by ensuring that the state tomography algorithm produces Haar-random error, so that with high probability the errors do not compound. Second, we can only estimate each column up to a phase, so we need to additionally deduce the relative phases between columns. We address this by showing that learning the columns up to phase of both $Z$ and $ZF$, $F$ being the discrete Fourier transform, suffices to deduce these relative phases (Proposition 2.3).

Achieving the optimal Heisenberg scaling.

We wish to improve the $O(1/\varepsilon^{2})$ dependence of process tomography to $O(1/\varepsilon)$. Our strategy will be to use process tomography as a subroutine, but boost the error some other way. Specifically, we give a bootstrapping procedure (Algorithm 1) that estimates a unitary to $\varepsilon$ error by making calls to a “base” algorithm that can only output unitary estimates to error $\frac{1}{200}$, paying an $1/\varepsilon$ factor of overhead in the query complexity (Theorem 3.3). Using the $O({\mathsf{d}}^{2}/\varepsilon_{0}^{2})$ unitary estimation algorithm with $\varepsilon_{0}=\frac{1}{200}$, as the base, this gives the desired $O({\mathsf{d}}^{2}/\varepsilon)$ query complexity.

A good first try is to get constant-error estimates for $Z^{2^{j}}$ for $j$ going from 1 to $\lceil\log_{2}(1/\varepsilon)\rceil$, with the hope that the estimate of $Z^{2^{j}}$ will refine the estimate of $Z$ to $2^{-j}$ error. Taking powers of an unknown gate with exponentially increasing degree is an existing method in quantum metrology for achieving the Heisenberg limit. In particular, we can view this as a non-coherent version of phase estimation [KLY15], where the application of $Z^{2^{j}}$ to an eigenvector of $Z$ extracts the $j$th bit of its corresponding eigenvalue. This argument can be shown to work in simple cases, such as when $Z$ is a rotation in two dimensions (Fig. 1), but if all the eigenvalues of $Z$ are $\pm 1$, then any power is either ${\mathbb{I}}$ or $Z$ so it is impossible to gain any better information about the eigenvectors of $Z$ from its powers.

Nonetheless, we find that if the powers of $Z$ are always close to the identity for all sufficiently high exponents (and hence no eigenvalue is close to $-1$), then this idea works out. Specifically, we show in Lemma 3.1 that if unitaries $U$ and $V$ are $\alpha$-close in diamond-norm distance, then $U^{1/p}$ and $V^{1/p}$ are $50\alpha/p$-close, provided that $U$ and $V$ are $0.01$-close to the identity. We can use this lemma to bootstrap a constant-error estimate to an $\varepsilon$-error one, provided we always apply the base algorithm to matrix powers that are close to the identity. Thus, we always bring the unknown unitary close to the identity: instead of running the base unitary estimation algorithm on $Z^{2^{j}}$, we run it on $(ZV_{j}^{\dagger})^{2^{j}}$, where $V_{j}$ is our best estimate to $Z$ so far in the algorithm. This recenters the unitary at the identity so that we can power it up even further. We note here that, upon formulating the right technical lemma to use (Lemma 3.1), the analysis of the resulting algorithm (Algorithm 1) becomes surprisingly simple.

Lower bound.

To achieve our lower bound of $\Omega({\mathsf{d}}^{2}/\varepsilon)$, we consider a hard instance of the unitary estimation problem, in which one is asked to identify one of $\exp(\Omega({\mathsf{d}}^{2}))$ candidate unitary channels that are $\varepsilon$-apart from one another and are all $O(\varepsilon)$-close to the identity in the diamond norm. Specifically, we choose the ensemble to be $\varepsilon$th-power of a net of reflections.

Simple arguments give a lower bound of $\Omega({\mathsf{d}}^{2}+\frac{1}{\varepsilon})$, with ${\mathsf{d}}^{2}$ being a generic lower bound for identifying one of $\exp(\Omega({\mathsf{d}}^{2}))$ candidate unitary channels [BMQ22] and $\frac{1}{\varepsilon}$ being the number of applications of a channel necessary to discriminate two channels that are $\varepsilon$-close in diamond-norm distance. To improve this, we wish to argue that, since a unitary $Z$ taken from the hard instance applies an $\varepsilon$th power of a reflection, the task of distinguishing $Z$ is a factor of $1/\varepsilon$ harder than the task of distinguishing the net of reflections. Intuitively, this holds because $Z$ is essentially the identity for all but an $\varepsilon$ fraction of the time it is applied.

This heuristic picture can be made precise by showing a reduction from “fractional-query” algorithms to “discrete-query” algorithms [CGM*+*09, BCC*+*17]. This reduction converts a circuit that calls $Z$, the $\varepsilon$th power of a reflection, into one that calls that reflection conditioned on some ancilla. So, hardness of a problem using queries to the reflection lifts to hardness of the problem using queries to $Z$. Translating to our setting, directly applying this reduction shows that a $\Omega(d^{2})$ to distinguishing a net of reflections translates to a $\Omega(\frac{{\mathsf{d}}^{2}}{\varepsilon}\frac{\log\log(1/\varepsilon)}{\log(1/\varepsilon)})$ lower bound for the hard instance. This reduction is tight, though, so removing the $\frac{\log\log(1/\varepsilon)}{\log(1/\varepsilon)}$ term requires additional insight.

We consider again the reduction of [CGM*+*09, BCC*+*17] converting a fractional-query circuit to a discrete-query circuit. This reduction adds one ancilla qubit per oracle call to $Z$ and the joint state of all the ancillas is a superposition of some bitstrings. The average weight of the bitstring is an effective number of the oracle calls to the reflection, which is smaller than the naive query complexity for $Z$ by a factor of $\varepsilon\ll 1$. Without much change in the output distribution of the overall algorithm, one can modify the algorithm to monitor the weight of the ancilla bitstring and actually reduce the number of oracle calls by the factor of $\varepsilon$. This reduction however comes at a price: the modified algorithm now requires postselection whose success probability is exponentially small in the query number. Addressing this small success probability is where one picks up the sub-logarithmic term. However, for this particular lower bound we need not address it: an exponentially small success probability is still hard to achieve for unitary discrimination. Bavaresco, Murao, and Quintino [BMQ22] gives a useful relation between the query complexity of a unitary discrimination problem and success probability, and we observe that small success probability is still meaningful as long as the algorithm is appreciably better than merely guessing the answer. This gives the optimal $\Omega({\mathsf{d}}^{2}/\varepsilon)$ lower bound.

1.5 Discussion

Improving gate complexity.

We give an algorithm for unitary channel estimation that is query-efficient and space-efficient. However, we still need to run quantum circuits with $\operatorname{poly}({\mathsf{d}},1/\varepsilon)$ depth, in addition to the oracle calls. If this could be improved to $\operatorname{poly}(\log({\mathsf{d}}),1/\varepsilon)$, it would make this algorithm significantly more practical. However, the depth is bottlenecked by the “shifting to identity” step, which requires a unitary synthesis and so is high depth. As discussed in Section 1.4, this step appears to be necessary to avoid the hard case of unitaries with $-1$ eigenvalues.

A simpler question is whether the $O({\mathsf{d}}^{2}/\varepsilon^{2})$ algorithm can be made gate-efficient. In view of this goal, we give a gate-efficient version of pure state tomography with optimal sample complexity in Appendix C; prior gate-efficient algorithms lose log factors in the sample complexity. This should be able to give a $O({\mathsf{d}}^{2}/\varepsilon^{2})$ query complexity algorithm with the desired gate complexity and optimal space complexity, though we do not prove this.

Estimating the eigenvalues of a unitary channel.

We show in Appendix B that shifting to the identity is not necessary when we merely want to learn the eigenvalues of the unitary without the eigenvectors. In this Appendix, we describe how to use a (gate-efficient) control-free version of phase estimation to achieve this result (though notably, the estimates we achieve do not distinguish between eigenvalues of different multiplicity). This algorithm appears to be folklore. If one can reduce unitary estimation to $\operatorname{poly}({\mathsf{d}})$ instances of eigenvalue estimation, then this result would give gate-efficient unitary channel tomography with query complexity $\operatorname{poly}({\mathsf{d}})/\varepsilon$.

Generalizing beyond unitary channels.

Finally, we note that our results do not extend to general channels. For example, the channel that destroys the input and outputs the outcome of a biased coin requires $\Omega(1/\varepsilon^{2})$ queries to estimate. However, we leave open the question of analyzing this algorithm’s tolerance to noise, and the related question of whether $o(1/\varepsilon^{2})$ process tomography is possible for “close-to-unitary” channels.

2 Base tomography using $O({\mathsf{d}}^{2}/\varepsilon^{2})$ applications

In this section we show a base tomography algorithm that has a quadratic dependence on the desired precision. We will use the “operator norm up to phase” distance for this section, defined in Definition 1.4, which is equivalent to diamond norm up to constants.

Theorem 2.1**.**

There is a tomography algorithm that, given access to an unknown unitary $Z\in{\mathsf{U}}({\mathsf{d}})$, as well as parameters $\varepsilon,\eta>0$, applies $Z$ at most $O({\mathsf{d}}^{2}/\varepsilon^{2})\cdot\log(1/\eta)$ times and outputs an estimate ${\bm{V}}$ satisfying $\operatorname{dist}_{\mathrm{phaseop}}(Z,{\bm{V}})\leq\varepsilon$ except with probability at most $\eta$. The algorithm uses only a ${\mathsf{d}}$-dimensional Hilbert space (in particular, only ${\mathsf{n}}$ qubits when ${\mathsf{d}}=2^{{\mathsf{n}}}$).

In fact, by standard methods (explicitly shown in Proposition 2.4 below) it suffices to prove this theorem for a fixed confidence value less than $\frac{1}{2}$, say $\eta=\frac{1}{3}$.

Our algorithm will essentially learn the unknown $Z$ column by column, using the below pure state tomography result, Proposition 2.2. This result was essentially previously known, but we will take care of a few minor details in Section 2.1.

Proposition 2.2**.**

There is a pure state tomography algorithm with the following behavior. Given access to copies of a pure state $\ket{z}\in{\mathbb{C}}^{{\mathsf{d}}}$, it sequentially and nonadaptively makes von Neumann measurements on $O({\mathsf{d}}/\varepsilon_{0})$ copies of $\ket{z}$ (using only ${\mathsf{d}}$-dimensional Hilbert space). Then, after classically collating and processing the measurement outcomes, it outputs (a classical description of) an estimate pure state

[TABLE]

such that: (i) ${\bm{\phi}}$ is a complex phase; (ii) the infidelity ${\bm{\varepsilon}}$ is at most $\varepsilon_{0}$ except with probability at most $\exp(-5{\mathsf{d}})\leq\frac{1}{100{\mathsf{d}}}$; (iii) the vector $\ket{{\bm{w}}}$ is distributed Haar-randomly999Perfect Haar-randomness is only possible by making idealized assumptions about the algorithm’s hardware; such technical issues of algorithmic complexity are deferred to Footnote 8. among all states orthogonal to $\ket{z}$.

This column-by-column technique has the minor downside that each column estimate can be off by a different complex phase ${\bm{\phi}}$. However, there are a few simple ways to work around this flaw; in particular, the following result (proven in Section 2.1 below) gives a completely black-box method:

Proposition 2.3**.**

Let $\operatorname{\mathcal{L}}$ be a tomography algorithm as in Theorem 2.1, except with the following weaker guarantee about ${\bm{V}}$:

[TABLE]

Then there is a tomography algorithm (using $\operatorname{\mathcal{L}}$ twice) that achieves Theorem 2.1, with $25\varepsilon$ and $2\eta$ in place of $\varepsilon$ and $\eta$, and $O({\mathsf{d}}^{3})$ additional classical time complexity.

Putting all of the above together, we can establish our base tomography result:

Proof of Theorem 2.1.

From the preceding discussion, it suffices to obtain a unitary tomography routine achieving Equation 21 with failure probability at most $\eta=\frac{1}{6}$. We apply the pure state tomography routine from Proposition 2.2, with $\varepsilon_{0}=\Omega(\varepsilon^{2})$ to be chosen later, on each of $\ket{z_{1}},\dots,\ket{z_{{\mathsf{d}}}}$, where $|z_{j}\rangle=Z|j\rangle$. This indeed uses $Z$ at most $O({\mathsf{d}}^{2}/\varepsilon^{2})$ times, within a Hilbert space of dimension only ${\mathsf{d}}$, and produces estimates $\ket{{\bm{u}}_{1}},\dots,\ket{{\bm{u}}_{{\mathsf{d}}}}$. Let ${\bm{U}}\in{\mathbb{C}}^{{\mathsf{d}}\times{\mathsf{d}}}$ denote the (possibly nonunitary) matrix with the $\ket{{\bm{u}}_{j}}$’s as columns; we will show that

[TABLE]

except with probability at most $\frac{1}{50}$. Then if we output any unitary ${\bm{V}}$ satisfying $\lVert{\bm{U}}-{\bm{V}}\rVert_{\mathrm{op}}\leq\varepsilon/2$ (for example, ${\bm{V}}={\bm{X}}{\bm{Y}}^{\dagger}$, where ${\bm{U}}={\bm{X}}{\bm{\Sigma}}{\bm{Y}}^{\dagger}$ is a singular value decomposition), then ${\bm{V}}$ will satisfy Equation 21, as desired.

To establish Equation 22, write

[TABLE]

for $j\in[d]$ as in Equation 20. Then with $\Phi=\operatorname{diag}({\bm{\phi}}_{1},\dots,{\bm{\phi}}_{d})$ and ${\bm{W}}$ the matrix with $\ket{{\bm{w}}_{j}}$’s as columns, we have

[TABLE]

We have

[TABLE]

except with probability at most $\frac{1}{100}$, by a union bound. We will also shortly show:

[TABLE]

(Here $C$ is a universal constant.) Combining the above, we conclude

[TABLE]

except with probability at most $\frac{1}{100}+\frac{2}{100}<\frac{1}{6}$, and this norm bound is at most $\varepsilon/2$ (as needed for Equation 22) provided we take the constant in $\varepsilon_{0}=\Omega(\varepsilon^{2})$ small enough.

It remains to prove the claim from Equation 26. Recall that ${\bm{W}}$ has independent unit columns $\ket{{\bm{w}}_{j}}$, with $\ket{{\bm{w}}_{j}}$ Haar-random orthogonal to $\ket{z_{j}}$. Introduce i.i.d. real random variables ${\bm{\delta}}_{1},\dots,{\bm{\delta}}_{{\mathsf{d}}}$, where ${\bm{\delta}}_{j}$ is distributed as $\lvert\braket{1}{{\bm{x}}}\rvert^{2}$ for $\ket{{\bm{x}}}$ a Haar-random unit vector in ${\mathbb{C}}^{{\mathsf{d}}}$. If we further introduce i.i.d. uniformly random complex phases ${\bm{\psi}}_{1},\dots,{\bm{\psi}}_{{\mathsf{d}}}$, then the unit vectors

[TABLE]

are in fact Haar-random and independent. Letting ${\bm{Y}}$ denote the matrix with the $\ket{{\bm{y}}_{j}}$’s as columns, it is a standard fact in random matrix theory101010For example, if ${\bm{Y}}$’s columns were independent Haar-random unit vectors in ${\mathbb{R}}^{{\mathsf{d}}}$ (as opposed to ${\mathbb{C}}^{{\mathsf{d}}}$) then [Ver18, Theorems 3.4.6, 4.6.1] would directly yield that $\lVert{\bm{Y}}\rVert_{\mathrm{op}}\leq K$ except with probability at most $2\exp(-cK^{2}{\mathsf{d}})$, for some constant $c>0$. The generalization to the complex case is very minor. that, for some universal constant $K$, we have

[TABLE]

We can now rewrite Equation 28 as

[TABLE]

from which we can conclude

[TABLE]

except with probability at most $\frac{1}{100}$. Thus to complete the proof of the claim in Equation 26, it suffices to show

[TABLE]

for some constant $K^{\prime}$. Note that for each constant value of ${\mathsf{d}}$, the random variable ${\bm{\delta}}_{j}$ has a continuous probability density on $[0,1]$; from this observation, it’s easy to deduce that it suffices to prove Equation 32 under the assumption ${\mathsf{d}}\geq{\mathsf{d}}_{0}$ for some constant ${\mathsf{d}}_{0}$. But this is easy: ${\bm{\delta}}_{j}^{2}$ has mean $1/{\mathsf{d}}$ and is sub-exponentially distributed with parameter $\Theta(1/{\mathsf{d}})$ [Ver18, Proposition 2.7.1, Lemma 2.7.6, Theorem 3.4.6]; hence $\Pr[{\bm{\delta}}_{j}^{2}\geq 1/2]\leq\exp(-\Omega({\mathsf{d}}))$ and—with a union bound—this is more than sufficient for Equation 32, once ${\mathsf{d}}$ is sufficiently large. ∎

2.1 Ancillary results for base tomography

We begin by giving a proof of the pure state tomography result we needed:

Proof of Proposition 2.2.

This result was essentially proven in [CL14, KRT17]. To be precise, we will refer to the analysis from [GKKT20, Theorem 2]. The result therein is exactly what we need, except for the following distinctions:

• •

Rather than making von Neumann measurements, [GKKT20] refer to performing the “uniform POVM” on each copy of $\ket{z}$; this is the continuous POVM with elements labeled by unit vectors $\ket{v}\in{\mathbb{C}}^{{\mathsf{d}}}$, in which the $\ket{v}$-element has density ${\mathsf{d}}\cdot\ket{v}\!\bra{v}$ with respect to Haar measure on $\ket{v}$. However, this is mathematically equivalent to first using classical randomness to choose a Haar-random ${\bm{V}}\sim{\mathsf{U}}({\mathsf{d}})$, and then projectively measuring in the basis of ${\bm{V}}$’s columns.

• •

The classical post-processing algorithms in [KRT17, GKKT20] do not necessarily output a pure (rank-one) hypothesis; they output a possibly mixed state ${\bm{\rho}}$, counting it as a success (in the case of [GKKT20]) when $\lVert|z\rangle\langle z|-{\bm{\rho}}\rVert_{1}\leq\sqrt{\varepsilon_{0}}$. This also means they do not explicitly confirm condition (iii) in Proposition 2.2, concerning the Haar-randomness of $\ket{{\bm{w}}}$.

But inspection of the actual algorithm in [GKKT20] shows that this second issue is easily fixed. The algorithm first forms ${\bm{L}}=({\mathsf{d}}+1)\mathrm{avg}\{\ket{{\bm{v}}_{j}}\!\bra{{\bm{v}}_{j}}\}-{\mathbb{I}}$, where the $\ket{{\bm{v}}_{j}}$’s are the measurement outcomes. With $O({\mathsf{d}}/\varepsilon_{0})$ measurements, this matrix is shown to satisfy $\lVert|z\rangle\langle z|-{\bm{L}}\rVert_{\mathrm{op}}\leq\frac{1}{4}\varepsilon_{0}$ except with probability at most $\exp(-5{\mathsf{d}})$. The authors of [GKKT20] then “round” ${\bm{L}}$ to a quantum state ${\bm{\rho}}$ by first diagonalizing it as ${\bm{W}}\operatorname{diag}({\bm{\lambda}}){\bm{W}}^{\dagger}$ for ${\bm{W}}\in{\mathsf{U}}({\mathsf{d}})$ and ${\bm{\lambda}}\in{\mathbb{C}}^{{\mathsf{d}}}$, and then taking ${\bm{\rho}}={\bm{W}}\operatorname{diag}({\bm{\lambda}}^{\prime}){\bm{W}}^{\dagger}$, where ${\bm{\lambda}}^{\prime}$ is the nearest probability vector to ${\bm{\lambda}}$. With this adjustment, they show that $\lVert\ket{u}\!\bra{u}-{\bm{\rho}}\rVert_{1}\leq\sqrt{\varepsilon_{0}}$ as needed.

Note that this last inequality implies that the closest rank-$1$ matrix ${\bm{M}}$ to ${\bm{\rho}}$ must satisfy $\lVert{\bm{M}}-{\bm{\rho}}\rVert_{1}\leq\sqrt{\varepsilon_{0}}$. On the other hand, it is well known that ${\bm{M}}$ is simply given by ${\bm{W}}\operatorname{diag}({\bm{\lambda}}^{\prime\prime}){\bm{W}}^{\dagger}$, where ${\bm{\lambda}}^{\prime\prime}$ is formed from ${\bm{\lambda}}^{\prime}$ by zeroing out all entries except the largest. It follows that the zeroed-out entries sum to at most $\sqrt{\varepsilon_{0}}$, and hence $\operatorname*{Tr}({\bm{M}})\geq 1-\sqrt{\varepsilon_{0}}$. If we now form $\widetilde{{\bm{M}}}=\frac{{\bm{M}}}{\operatorname*{Tr}({\bm{M}})}$, then $\widetilde{{\bm{M}}}$ is a rank-$1$ state with $\lVert\widetilde{{\bm{M}}}-{\bm{M}}\rVert_{1}\leq O(\sqrt{\varepsilon_{0}})$ and hence $\lVert\widetilde{{\bm{M}}}-{\bm{\rho}}\rVert_{1}\leq O(\sqrt{\varepsilon_{0}})$ and $\lVert\ket{u}\!\bra{u}-\widetilde{{\bm{M}}}\rVert_{1}\leq O(\sqrt{\varepsilon_{0}})$. Thus we may express $\widetilde{{\bm{M}}}=\ket{\widehat{{\bm{u}}}}\!\bra{\widehat{{\bm{u}}}}$ and output $\ket{\widehat{{\bm{u}}}}$ (after slightly adjusting the constant on $\varepsilon_{0}$).

It only remains to observe that the process of forming ${\bm{L}},{\bm{\rho}},{\bm{M}},\widetilde{{\bm{M}}},\ket{\widehat{{\bm{u}}}}$ is completely symmetric with respect to the subspace orthogonal to $\ket{u}$, and hence the “error vector” $\ket{{\bm{w}}}$ is indeed distributed Haar-randomly. ∎

Next, we give the algorithm for Proposition 2.3, which fixes the column phases for our unitary tomography algorithm:

Proof of Proposition 2.3.

Let $F\in{\mathsf{U}}({\mathsf{d}})$ denote the discrete Fourier transform, so $\braket{a}{F}{b}=(1/\sqrt{{\mathsf{d}}})\exp(-2\pi\mathrm{i}ab/{\mathsf{d}})$. (The only property we use about $F$ is that each entry has the same magnitude; if ${\mathsf{d}}=2^{\mathsf{n}}$, one may use the Hadamard transform $H^{\otimes{\mathsf{n}}}$ instead.) Given $Z$, our algorithm applies $\operatorname{\mathcal{L}}$ once to $Z$ and once to $ZF^{\dagger}$; call the results ${\bm{V}},{\bm{G}}\in{\mathsf{U}}({\mathsf{d}})$, respectively. Except with probability at most $2\eta$, we get that there exist diagonal unitaries ${\bm{\Phi}}_{\bm{V}}$ and ${\bm{\Phi}}_{\bm{G}}$ such that

[TABLE]

Given this, we conclude that ${\bm{G}}^{\dagger}{\bm{V}}$ is close to $F$, up to phases on rows and columns.

[TABLE]

We can then deduce ${\bm{\Phi}}_{\bm{V}}$ by computing ${\bm{G}}^{\dagger}{\bm{V}}$ in $O({\mathsf{d}}^{3})$ time, and then essentially reading off the relative column phases. First, notice that by pigeonhole principle applied to Eq. 34,

[TABLE]

Let $C(a,b)$ denote the inequality in Eq. 35. Since $\braket{a}{{\bm{\Phi}}_{\bm{G}}^{\dagger}F{\bm{\Phi}}_{\bm{V}}}{b}=\overline{{\bm{\phi}}}_{\bm{G}}^{(a)}{\bm{\phi}}_{{\bm{V}}}^{(b)}\braket{a}{F}{b}$ has magnitude $\frac{1}{\sqrt{{\mathsf{d}}}}$, we can define $P\in\mathbb{C}^{{\mathsf{d}}\times{\mathsf{d}}}$ to be $\braket{a}{P}{b}=\braket{a}{{\bm{G}}^{\dagger}{\bm{V}}}{b}/\braket{a}{F}{b}$, and it follows that

[TABLE]

By Eq. 35, $C(a,1)$ and $C(a,b)$ hold for at least half of $a$’s. Let ${\bm{\psi}}_{b}$ denote the coordinate-wise (real and imaginary) median of $\{\braket{a}{P}{b}/\braket{a}{P}{1}\}_{a\in[{\mathsf{d}}]}$. Then by Eq. 37,

[TABLE]

Let ${\bm{\Psi}}$ be the matrix with ${\bm{\psi}}_{b}$ on the diagonal. Then $\lVert{\bm{\phi}}_{\bm{V}}^{(1)}{\bm{\Psi}}-{\bm{\Phi}}_{\bm{V}}\rVert_{\mathrm{op}}\leq 24\varepsilon$. Combining this with $\lVert Z{\bm{\Phi}}_{\bm{V}}-{\bm{V}}\rVert_{\mathrm{op}}\leq\varepsilon$ from Equation 33 easily yields

[TABLE]

and so our algorithm may output ${\bm{V}}{\bm{\Psi}}^{\dagger}$. ∎

Finally, we give the (completely standard) trick for boosting confidence:

Proposition 2.4**.**

Let $\mathcal{L}$ be a learning algorithm for objects in a metric space with efficiently computable distance $\mathrm{distance}({\cdot},{\cdot})$. Assume that on input $Z$, the algorithm $\mathcal{L}$ outputs ${\bm{V}}$ satisfying $\mathrm{distance}(Z,{\bm{V}})\leq\varepsilon$ except with probability at most $.49$. Then there is an algorithm $\mathcal{L}^{\prime}$ that takes an additional input $\eta>0$, uses $\mathcal{L}$ at most $O(\log(1/\eta))$ times, and guarantees $\mathrm{distance}(Z,{\bm{V}})\leq 3\varepsilon$ except with probability at most $\eta$.

Proof.

For $T=O(\log(1/\eta))$ times, run $\mathcal{L}$ independently, obtaining ${\bm{V}}_{1},\dots,{\bm{V}}_{T}$. By a standard Chernoff bound, except with probability at most $\eta$, there are at least $.505T$ “good” estimates, where we say ${\bm{V}}_{j}$ is “good” if $\operatorname{dist}(Z,{\bm{V}}_{j})\leq\varepsilon$. By the triangle inequality, every good estimate is also “central”, where we say estimate ${\bm{V}}_{j_{0}}$ is “central” if it has the following property: $\operatorname{dist}({\bm{V}}_{j_{0}},{\bm{V}}_{k})\leq 2\varepsilon$ for at least $.505T$ estimates ${\bm{V}}_{k}$. Let $\mathcal{L}$ now select and output any central estimate ${\bm{V}}_{j^{*}}$; the method of brute-force checking central-ness has time complexity $O(T^{2})$ times the cost of computing a distance, which is $O({\mathsf{d}}^{3})$ for ${\mathsf{d}}\times{\mathsf{d}}$ matrices. Since $.505T+.505T>T$, the Pigeonhole Principle implies that at least one of the estimates ${\bm{V}}_{k}$ for which $\operatorname{dist}({\bm{V}}_{j^{*}},{\bm{V}}_{k})\leq 2\varepsilon$ is also good. Thus the triangle inequality implies $\operatorname{dist}({\bm{V}}_{j^{*}},Z)\leq 3\varepsilon$, as desired. ∎

3 Bootstrap of precision to $O({\mathsf{d}}^{2}/\varepsilon)$ applications

3.1 Key lemma: Geometry of unitary groups

The purpose of this subsection is to recall some well-known metric notions on unitary groups and projective unitary groups; see [Sza97] and [KSV02, §8.3.3]. Only Lemma 3.1 will be used in later sections. Readers who are familiar with the intrinsic metric induced by operator norm may quickly proceed to the next subsection.

Consider the operator norm on the Lie algebra ${\mathfrak{u}}({\mathsf{d}})$ of all ${\mathsf{d}}$-by-${\mathsf{d}}$ antihermitian matrices. By demanding left- and right-invariance we obtain a metric on a Lie group ${\mathsf{U}}({\mathsf{d}})$, as defined in Definition 1.5. Recalling this definition, the length of a smooth path $\gamma:[0,1]\to{\mathsf{U}}({\mathsf{d}})$ is given by

[TABLE]

where $\gamma^{\prime}(t)\gamma(t)^{-1}\in{\mathfrak{u}}({\mathsf{d}})$. The distance between two points of ${\mathsf{U}}({\mathsf{d}})$ is the infimum of the lengths of all smooth paths connecting the two:

[TABLE]

This definition makes it obvious that $\operatorname{dist}$ is a metric, obeying the triangle inequality. Though $\operatorname{dist}$ is equivalent to the “extrinsic” metric $\lVert\cdot\rVert_{\mathrm{op}}$, we will use the intrinsic metric in this section because it leads us to think in terms of the Lie algebra.

For any $U\in{\mathsf{U}}({\mathsf{d}})$ and any real number $r>0$, we define an open metric ball of radius $r$ centered at $U$ by

[TABLE]

As the notation suggests, due to the right invariance of the metric, the ball ${\mathcal{B}}(r)U$ is a shift of ${\mathcal{B}}(r){\mathbb{I}}$. We will write ${\mathcal{B}}(r)$ for ${\mathcal{B}}(r){\mathbb{I}}$.

Every unitary quantum channel is defined by a unitary, but specifies the unitary only up to a global phase factor. This motivates us to consider projective unitary groups ${\mathsf{U}}({\mathsf{d}})/{{\mathsf{U}}(1)}$ where ${{\mathsf{U}}(1)}$ is the center of ${\mathsf{U}}({\mathsf{d}})$ consisting of phase factors. The dimension of a matrix in the center ${{\mathsf{U}}(1)}$ is implicit. The metric $\operatorname{dist}$ induces a metric on the projective unitary group for which we use the same notation. In analogy with ${\mathcal{B}}$, we write

[TABLE]

Strictly speaking, this is not a metric ball of the projective unitary group since ${\mathcal{PB}}(r)$ is a subset of ${\mathsf{U}}({\mathsf{d}})$, but this will hardly cause any confusion below.

We consider the fractional power of a unitary in a small neighborhood of ${\mathbb{I}}$ for any real $r>0$. The small neighborhood is actually ${\mathcal{B}}(\pi)$ and we define

[TABLE]

for any $X\in{\mathfrak{u}}({\mathsf{d}})$ with $\lVert X\rVert_{\mathrm{op}}<\pi$, which we remark is a strict inequality. This is a proper definition because the exponential map is injective in the corresponding neighborhood of $0\in{\mathfrak{u}}({\mathsf{d}})$; see Lemma 3.2 below. We will use the following lemma in our accuracy boosting algorithm in the next subsection.

Lemma 3.1**.**

For any $U,V\in{\mathcal{PB}}(\tfrac{1}{3\pi})$ and $p\geq 1$,

[TABLE]

The proof of this lemma will use the following.

Lemma 3.2** (Following Eq. 5, Lemma 3, and Lemma 4 of [Sza97]).**

For any $U\in{\mathsf{U}}({\mathsf{d}})$, there is $X\in{\mathfrak{u}}({\mathsf{d}})$ with $\lVert X\rVert_{\mathrm{op}}\leq\pi$ such that $U=e^{X}$. If $e^{X}=e^{Y}$ with $\lVert X\rVert_{\mathrm{op}},\lVert Y\rVert_{\mathrm{op}}<\pi$, then $X=Y$. Further, the following hold:

1. (a)

For any $X\in{\mathfrak{u}}({\mathsf{d}})$ such that $\lVert X\rVert_{\mathrm{op}}\leq\pi$, $\operatorname{dist}({\mathbb{I}},e^{X})=\lVert X\rVert_{\mathrm{op}}$. 2. (b)

For any $U,V\in{\mathsf{U}}({\mathsf{d}})$, $\lVert U-V\rVert_{\mathrm{op}}\leq\operatorname{dist}(U,V)\leq\frac{\pi}{2}\lVert U-V\rVert_{\mathrm{op}}$. 3. (c)

For any $X,Y\in{\mathfrak{u}}({\mathsf{d}})$, $\lVert e^{X}-e^{Y}\rVert_{\mathrm{op}}\leq\lVert X-Y\rVert_{\mathrm{op}}$. 4. (d)

For any $X,Y\in{\mathfrak{u}}({\mathsf{d}})$ such that $\lVert X\rVert_{\mathrm{op}},\lVert Y\rVert_{\mathrm{op}}\leq\frac{1}{\pi}$, $\lVert X-Y\rVert_{\mathrm{op}}\leq\pi\lVert e^{X}-e^{Y}\rVert_{\mathrm{op}}$.

We did not optimize the constant $\pi$ in Item 3.2(d).

Proof of Lemma 3.1.

We first prove the non-projective version of the statement: consider $U,V\in{\mathcal{B}}(1/\pi)$ and $p\geq 1$. Write $U=e^{X}$ and $V=e^{Y}$ with $\lVert X\rVert_{\mathrm{op}},\lVert Y\rVert_{\mathrm{op}}<1/\pi$. Then

[TABLE]

The inequalities follow from Item 3.2(b), Item 3.2(c), Item 3.2(d), and Item 3.2(b), respectively.

Now, to prove the lemma, without loss of generality let $U,V$ be their representatives in ${\mathcal{B}}(\tfrac{1}{3\pi})$, and let $e^{\mathrm{i}\theta}$ be the global phase minimizing $\operatorname{dist}(U,V{{\mathsf{U}}(1)})$. Then

[TABLE]

So, $Ue^{-\mathrm{i}\theta/2},Ve^{\mathrm{i}\theta/2}\in{\mathcal{B}}(1/\pi)$, so we can use Eq. 46.

[TABLE]

Proof of Lemma 3.2.

The eigenvalues of a unitary are on the complex unit circle, each of which can be specified uniquely by an angle in $[-\pi,\pi)$. So, the exponential map ${\mathfrak{u}}({\mathsf{d}})\ni X\mapsto e^{X}\in{\mathsf{U}}({\mathsf{d}})$ is injective on the domain where $\lVert X\rVert_{\mathrm{op}}<\pi$. This proves the first two statements.

Item 3.2(a): Let $\theta_{j}\in[-\pi,\pi)$ be the eigenvalues of $-\mathrm{i}X$. Without loss of generality, suppose $\ell=\lvert\theta_{1}\rvert\geq\lvert\theta_{j}\rvert$ for all $j$. Then the path $\gamma$ taking $t\mapsto\exp(\mathrm{i}t\operatorname{diag}(\theta_{1},\ldots,\theta_{\mathsf{d}}))$ goes from ${\mathbb{I}}$ to $e^{X}$ and has length $\ell$. This shows that $\operatorname{dist}({\mathbb{I}},e^{X})\leq\ell$.

To show that $\operatorname{dist}({\mathbb{I}},e^{X})\geq\ell$, consider a smooth path $\gamma:[0,1]\to{\mathsf{U}}({\mathsf{d}})$ joining ${\mathbb{I}}$ and $e^{X}$. The first column of a unitary matrix is a ${\mathsf{d}}$-complex-dimensional unit vector, which corresponds to a point in the standard unit sphere $S^{2{\mathsf{d}}-1}\subset{\mathbb{C}}^{{\mathsf{d}}}$. So, the first column of the matrix $\gamma(t)$, denoted $[\gamma(t)]_{1}$, defines a smooth curve $\xi:[0,1]\to S^{2d-1}$. Since $\lVert A\rVert_{\mathrm{op}}=\sup_{v:\lVert v\rVert_{2}=1}\lVert Av\rVert_{2}\geq\lVert A_{1}\rVert_{2}$ for any matrix $A$, we conclude that the length of $\gamma$ is at least the length of the path $\xi$ under the standard Euclidean metric:

[TABLE]

The path $\xi$ connects the two points $(1,0,\ldots,0),(e^{\mathrm{i}\theta_{1}},0,\ldots,0)\in S^{2d-1}$. It is well known that the length of $\xi$ is at least the (smaller) angle between $1$ and $e^{\mathrm{i}\theta_{1}}$ in the complex plane, which is $\ell$.111111 Since the standard sphere is a closed Riemannian manifold, one can appeal to the Hopf–Rinow theorem to obtain a geodesic realizing the distance between any pair of points, and characterize geodesics $g$ by the geodesic equation $\nabla_{\dot{g}}\dot{g}=0$, to conclude that a path of minimum length between any pair of points must be on a great circle.

Item 3.2(b): By unitary invariance of the two metrics, it suffices to prove the statement for $V={\mathbb{I}}$. For $\theta\in(-\pi,\pi]$, we observe that $\lvert 1-e^{\mathrm{i}\theta}\rvert^{2}=2-2\cos\theta=4\sin^{2}\frac{\theta}{2}$, implying that $\lvert 1-e^{\mathrm{i}\theta}\rvert=\lvert 2\sin(\tfrac{\theta}{2})\rvert\leq\lvert\theta\rvert\leq\lvert\pi\sin(\tfrac{\theta}{2})\rvert=\frac{\pi}{2}\lvert 1-e^{\mathrm{i}\theta}\rvert.$ It follows by Item 3.2(a) that $\lVert{\mathbb{I}}-U\rVert_{\mathrm{op}}\leq\operatorname{dist}({\mathbb{I}},U)\leq\frac{\pi}{2}\lVert{\mathbb{I}}-U\rVert_{\mathrm{op}}$.

Item 3.2(c):

[TABLE]

Item 3.2(d): Note that for any $X,Y\in{\mathfrak{u}}({\mathsf{d}})$, we have

[TABLE]

This is a well-known inequality; see e.g. Eq. (143) in arXiv version of [CST*+*21]. To prove it, notice that $\lVert e^{X}Ye^{-X}-Y\rVert_{\mathrm{op}}=\lVert\int_{0}^{1}{\mathrm{d}}s\,\partial_{s}[e^{sX}Ye^{-sX}]\rVert_{\mathrm{op}}\leq\int_{0}^{1}{\mathrm{d}}s\,\lVert e^{sX}(XY-YX)e^{-sX}\rVert_{\mathrm{op}}=\lVert[X,Y]\rVert_{\mathrm{op}}$, so therefore,

[TABLE]

Then,

[TABLE]

Rearranging, we complete the proof. ∎

3.2 Bootstrap algorithm

Using Lemma 3.1, we can show that we can bootstrap a base tomography algorithm that achieves constant error to an algorithm that gets $\varepsilon^{-1}\log(1/\eta)$ query complexity. Following Section 3.1, for this section we use the distance metric between two unitaries $U,V\in{\mathsf{U}}({\mathsf{d}})$ of $\operatorname{dist}(U,V{\mathsf{U}}(1))$ (Definition 1.5). By Proposition 1.6, this is equivalent to diamond-norm distance up to universal constants.

Theorem 3.3**.**

Suppose we have an oracle capable of applying an unknown unitary channel $Z\in{\mathsf{U}}({\mathsf{d}})$. Further suppose we have an algorithm $\operatorname{\mathcal{A}}$ that, given such an oracle, can output a unitary ${\bm{U}}$ such that $\operatorname{dist}(Z,{\bm{U}}{{\mathsf{U}}(1)})\leq\varepsilon_{0}\leq\frac{1}{600}$ with probability $>0.51$. Then, given error parameters $\varepsilon,\eta\in(0,1)$, Algorithm 1 outputs a unitary ${\bm{U}}$ such that

1. (a)

$\operatorname{dist}(Z,{\bm{U}}{{\mathsf{U}}(1)})\leq\varepsilon$* with probability $\geq 1-\eta$;* 2. (b)

and $\operatorname*{\mathbf{E}}[\operatorname{dist}(Z,{\bm{U}}{{\mathsf{U}}(1)})^{2}]\leq(1+32\eta)\varepsilon^{2}$.

Algorithm 1* has the further property that, if $\operatorname{\mathcal{A}}$ uses $Q$ queries to $Z$, then Algorithm 1 requires only $O(\frac{Q}{\varepsilon}\log\frac{1}{\eta})$ queries to $Z$.*

By plugging in the base tomography algorithm from Theorem 2.1 into this bootstrap algorithm, we obtain our main result Theorem 1.1. We restate the theorem now.

See 1.1

Proof of Theorem 1.1.

It is clear that our bootstrap algorithm needs as many quantum registers as the base tomography algorithm does. The base tomography algorithm combines $O({\mathsf{d}})$ pure state tomography outcomes, each of which uses projective measurements in computational basis on a pure state of form $A(ZV^{\dagger})^{p}B\ket{0}$ where $A,B,V$ are unitaries known at the moment of the state tomography. See Theorem 2.1.

The query complexity and the accuracy guarantees follow from Theorem 3.3, with Theorem 2.1 as the base algorithm.

The number of one- and two-qubit gates used in the overall algorithm beyond the oracle calls, is determined by the complexity of implementing interspersing unitaries $A,B,V$ in the preparation of states $A(ZV^{\dagger})^{p}B\ket{0}$. Those unitaries are on a ${\mathsf{d}}$-dimensional qudit, so it can be compiled to accuracy $\alpha$ in operator norm by $\operatorname{poly}({\mathsf{d}},\log(1/\alpha))$ elementary gates [KSV02]. We only need $\alpha=1/\operatorname{poly}(d/\varepsilon)$. The algorithm prepares $\operatorname{poly}({\mathsf{d}}/\varepsilon)$ states of form $A(ZV^{\dagger})^{p}B\ket{0}$ in total.

The classical time complexity is bounded by those of matrix multiplications, matrix diagonalizations, and finding diagonal phase factors in Proposition 2.3. All of these take $\operatorname{poly}({\mathsf{d}}/\varepsilon)$ time.

The bound on $\lVert\mathcal{M}-\mathcal{U}(Z)\rVert_{\diamond}$ follows from Theorem 3.3, which shows that Algorithm 1 gives $\operatorname*{\mathbf{E}}_{\bm{X}}\lVert{\bm{X}}\rVert_{\mathrm{op}}^{2}\leq O(\varepsilon^{2})$, and Proposition 1.18. Rescaling $\varepsilon$ if necessary, we complete the proof. ∎

The bootstrap algorithm (Algorithm 1) is essentially four lines long: we begin with $V_{0}={\mathbb{I}}$ as an initial estimate for $Z$, and with every iteration, we refine this estimate. In iteration $j$, we take the current estimate ${\bm{V}}_{j}$, and consider the “residual” $Z{\bm{V}}_{j}^{\dagger}$: the residual is close to the identity when ${\bm{V}}_{j}$ is close to $Z$. We run $\operatorname{\mathcal{A}}$ to get an estimate ${\bm{U}}_{j}$ of $(Z{\bm{V}}_{j}^{\dagger})^{p_{j}}$. Though this estimate ${\bm{U}}_{j}$ is only good to constant error, by Lemma 3.1, its $p_{j}$th root will be $(1/p_{j})$-close to $Z{\bm{V}}_{j}^{\dagger}$. So, we can use this root to form the next estimate ${\bm{V}}_{j+1}={\bm{U}}_{j}^{1/p_{j}}{\bm{V}}_{j}$, which is $(1/p_{j})$-close to $Z$. After $\log_{2}(1/\varepsilon)$ iterations, this estimate is good enough to output. We have just given nearly a full proof of the first part of Theorem 3.3. We swept two details under the rug: first, we need that $(Z{\bm{V}}_{j}^{\dagger})^{p_{j}}$ is close enough to the identity to satisfy the conditions for Lemma 3.1, which follows from induction; and we need to account for failure. We deal with these details below.

Proof of Theorem 3.3.

Using the trick in Proposition 2.4, we assume that our oracle $\operatorname{\mathcal{A}}$ can boost its success probability to $\geq 1-\gamma$ with only $O(\log\frac{1}{\gamma})$ overhead to output a unitary that’s correct up to $c=\frac{1}{200}$ error.

We analyze the output of Algorithm 1 for arbitrary $T>0$, first assuming that Algorithm 1 is successful for all $j$. Let

[TABLE]

be the error after iteration $j-1$ of the while-loop. We will prove ${\bm{\delta}}_{k}<2^{-k-4}$. For iteration [math], we have $p_{0}=1$ so by our assumption on $\operatorname{\mathcal{A}}$,

[TABLE]

We proceed by induction in $k$. First, notice that the input and output of the call to $\operatorname{\mathcal{A}}$ on the $k$th iteration are in ${\mathcal{PB}}(\frac{1}{3\pi})$:

[TABLE]

So, we can apply Lemma 3.1 to conclude that

[TABLE]

Therefore, the final accuracy guarantee is

[TABLE]

The query complexity is sum of those at Algorithm 1 for all $j$: For any given $j$, this call queries the oracle for $(Z{\bm{V}}_{j}^{\dagger})^{p_{j}}$ for $O(Q\log\frac{1}{\eta_{j}})$ times, which can be performed using $O(Qp_{j}\log\frac{1}{\eta_{j}})$ queries to $Z$. Summing over all the iterations, we get the desired total query complexity,

[TABLE]

Notice that $\eta_{j}$ is chosen so the query complexity does not incur the $\log\log(1/\varepsilon)$ factor that comes from a naive union bound. The failure probability of Algorithm 1 for each $j$ is $\eta_{j}$, so the probability that any of them fail is bounded by $\eta$, proving Item 3.3(a):

[TABLE]

Now, we prove Item 3.3(b). Item 3.3(a) only implies a bound of $\operatorname*{\mathbf{E}}[\operatorname{dist}(Z,{\bm{U}}{\mathsf{U}}(1))]\leq\varepsilon^{2}+\frac{\pi}{2}\eta$, since when the algorithm fails, the output unitary could be as bad as possible, $\pi/2$ away from $Z$. We prove now that Algorithm 1 fails gracefully: namely, if the algorithm first fails in iteration $f$, then the output only has error $O(2^{-f}$. Let $f$ be the least index of the iteration where the base tomography algorithm fails, so that $\Pr[f=0]<\eta_{0}$, and $\Pr[f=j\mid f>j-1]<\eta_{j}$ for $j=1,2,3,\ldots$. If the algorithm never fails, we take $f=T+1$. The argument we described above still applies up to index $f$, so that

[TABLE]

We can further conclude that the output of the algorithm will be about as good of an estimate to $Z$ as ${\bm{V}}_{f+1}$, since the size of the adjustment in each iteration decays exponentially regardless of failure.121212 In the following computation, we use that $\operatorname{dist}({\mathbb{I}},{\bm{U}}_{k}^{1/p_{k}}{{\mathsf{U}}(1)})\leq\pi 2^{-k}$, which holds because the eigenvalues of $U^{1/p_{k}}$ are in $(-\pi/p_{k},\pi/p_{k})$ for all $U\in{\mathsf{U}}({\mathsf{d}})$ with eigenvalues in $(-\pi,\pi)$. However, there is a technical issue: Algorithm 1 cannot be performed if ${\bm{U}}_{k}$ has $-1=e^{\mathrm{i}\pi}$ as an eigenvalue, since in this case, the $1/p_{j}$-th root of this arbitrary output is not defined. One can and should ignore such nongeneric behavior, but strictly speaking our algorithm is not defined for such cases. Formally, one should address this by changing the choice of discretization (Footnote 8) to ensure that the algorithm never outputs an estimate with $-1$ as an eigenvalue. This can be done since it only excludes a measure-zero set.

[TABLE]

This suffices to give a bound on the quality of the output $V_{T+1}$ in expectation.

[TABLE]

Remark 3.4**.**

Our bootstrap crucially depends on the continuous group structure of the (projective) unitary group since we take fractional roots of a result from base tomography and the the fractional root needs to in the set of candidate outputs. It is not important whether the set of candidate outputs is a unitary group or a projective unitary group. Conversely, our bootstrap works for any closed Lie subgroup of the (projective) unitary group (e.g., an orthogonal group $O({\mathsf{d}})\subset{\mathsf{U}}({\mathsf{d}})$) and an obvious analog of Theorem 3.3 is proved verbatim. Note that if a subgroup is not closed, then the intrinsic metric $\operatorname{dist}$ and the extrinsic one $\lVert\cdot\rVert_{\mathrm{op}}$ are not equivalent and our bootstrap will not work. Also note that some care would be needed in the base tomography; the Fourier or Hadamard transform in the proof of Proposition 2.3 might not belong to the Lie subgroup in consideration.

4 Lower bounds

The goal of this section is to prove Theorem 1.2. The first subsection will prove a $\Omega({\mathsf{d}}^{2})$ lower bound using existing results on unitary discrimination [BMQ22], which we lift to a $\Omega({\mathsf{d}}^{2}/\varepsilon)$ lower bound in the following section through a variant of an argument proving the equivalence of the fractional and discrete query models [CGM*+*09, BCC*+*17].

4.1 Bound for constant $\varepsilon$

As a warm up, we first prove a query complexity lower bound for constant $\varepsilon$. We begin by constructing a hard instance of the estimation problem, which is a discrete set of unitaries in which each pair of unitaries is at least distance $1/4$ apart. For this instance, the task of estimating each unitary to error $<1/8$ is equivalent to the task of simply identifying the unknown unitary.

Proposition 4.1**.**

There exists a set ${\mathcal{N}}=\{W_{j}\in{\mathsf{U}}({\mathsf{d}}):j=1,\ldots,N\}$ with $N\geq\exp({\mathsf{d}}^{2}/64)$ of Hermitian unitaries (i.e., the eigenvalues are $\pm 1$) such that $\lVert\mathcal{U}(W_{j})-\mathcal{U}(W_{k})\rVert_{\diamond}\geq 1/4$ for any $j\neq k$.

The constants 64 and $1/4$ are not sharp. A construction of some packing net of unitaries goes back at least to [Sza83]. Here we construct a special net of unitaries that share the same real eigenspectrum.

Proof.

Considering $2$-by-$2$ Pauli matrices embedded into a higher dimensional unitary group, we see that the claim is true for ${\mathsf{d}}\leq 8$. Let ${\mathsf{d}}=2r+2>8$ if ${\mathsf{d}}$ is even, or ${\mathsf{d}}=2r+1>8$ if ${\mathsf{d}}$ is odd.

We first show the existence of a set of far apart unitaries in operator norm. Lemma 8 of [HHJ*+*17] shows that there are $\geq\exp(r^{2}/8)$ density matrices of rank $r$ in dimension $2r$ where exactly half the eigenvalues are $1/r$ and the other half are zero, and any two different density matrices in this set are at least $1/4$-apart in trace distance. We may write this set as $\{({\mathbb{I}}_{2r}+V_{j})/(2r):j=1,\ldots,N\}$ with $N\geq\exp(r^{2}/8)$ where $V_{j}\in{\mathsf{U}}(2r)$ is a Hermitian unitary of trace zero. If $j\neq k$, we have

[TABLE]

Therefore, the set $S=\{V_{j}\in{\mathsf{U}}(2r):j=1,\ldots,N\}$ consists of Hermitian unitaries that are far apart in operator norm. We now need to define unitaries that are far apart in the norm used in the statement of the proposition.

Now, embed $S$ into ${\mathsf{U}}({\mathsf{d}})$ as $V_{j}\mapsto W_{j}=V_{j}\oplus{\mathbb{I}}_{b}\in{\mathsf{U}}({\mathsf{d}})$, where $b=1$ or $2$. We claim that $\{W_{j}:j=1,\ldots,N\}$ is a desired set. Consider $W_{j}$ and $W_{k}$ for $j\neq k$. Then, for all $\phi\in{{\mathsf{U}}(1)}$,

[TABLE]

Minimizing over all $\phi$ and using Eq. 7, we get that $\lVert\mathcal{U}(W_{j})-\mathcal{U}(W_{k})\rVert_{\diamond}\geq\operatorname{dist}_{\mathrm{phaseop}}(W_{j},W_{k})\geq\tfrac{1}{4}$. Finally, $N\geq\exp(r^{2}/8)\geq\exp({\mathsf{d}}^{2}/64)$ for ${\mathsf{d}}>8$. ∎

The following result upper bounds the average success probability of any unitary identification algorithm that makes at most $Q$ uses of the unknown unitary:

Proposition 4.2** (Theorem 5 of [BMQ22]).**

Let $U,\ldots,U_{N}\in{\mathsf{U}}({\mathsf{d}})$ by any unitary channels. For any unitary estimation algorithm $\operatorname{\mathcal{A}}$ that queries an input unitary channel $Q$ times and outputs an index $\hat{x}$ given $U_{x}$ with probability $\Pr[\hat{x}|x]$ where $x,\hat{x}\in[N]$, it holds that

[TABLE]

Ref. [BMQ22] contains a proof that relies on results from several references therein. For the readers’ convenience, we present an overview of the proof in Appendix D.

From this, it straightforwardly follows that if we use the set ${\mathcal{N}}$ we just constructed, if $Q$ is much smaller than $d^{2}$, then we cannot solve the identification problem with high probability. But we can also prove something stronger:

Lemma 4.3**.**

If $\frac{1}{N}\binom{Q+{\mathsf{d}}^{2}-1}{Q}\geq\exp(-cQ)$ where $Q\geq 1$, ${\mathsf{d}}\geq 100/c^{\prime}\geq 1$ and $N\geq\exp(c^{\prime}{\mathsf{d}}^{2})$ for some constants $c,c^{\prime}>0$, then $Q\geq c^{\prime\prime}{\mathsf{d}}^{2}$ for some constant $c^{\prime\prime}$ that depends only on $c,c^{\prime}$.

In combination with Proposition 4.2, this immediately implies a $\Omega({\mathsf{d}}^{2})$ query lower bound for any unitary estimation algorithm that estimates to small constant distance with constant success probability, since one can use it to distinguish between the $\exp(\Omega({\mathsf{d}}^{2}))$ unitaries in the net from Proposition 4.1 with constant success probability. We do not use the full power of this lemma to draw this conclusion, but we will need it in its full form in the next subsection.

Proof.

Suppose for a contradiction $1\leq Q\leq\delta{\mathsf{d}}^{2}$ for some $\delta\in(0,1)$. It is easy to show from Stirling’s approximation that $\binom{x}{sx}\leq xe^{xH(s)}$ for all $x\geq 1$ where $H(s)=-s\ln(s)-(1-s)\ln(1-s)$ is the binary entropy function. Then, $e^{c^{\prime}{\mathsf{d}}^{2}-cQ}\leq\binom{Q+{\mathsf{d}}^{2}-1}{Q}\leq e^{H(\frac{\delta}{1+\delta})(1+\delta){\mathsf{d}}^{2}}(1+\delta){\mathsf{d}}^{2}$. Taking logs and dividing by ${\mathsf{d}}^{2}$, we get

[TABLE]

For ${\mathsf{d}}\geq 100/c^{\prime}$, the left-hand side is at least $c^{\prime}/2$. This completes the proof. ∎

4.2 Bound for general $\varepsilon$

We’re now ready to prove the stronger $\varepsilon$-dependent lower bound. Our high-level strategy is as follows. First, we construct a special set of unitaries that are pairwise $\varepsilon$ apart using our constructed set ${\mathcal{N}}$. This constructed set brings all the unitaries in ${\mathcal{N}}$ closer to the identity (and hence each other) by raising each unitary to a small power $\alpha$, which will be of the order of $\varepsilon$. Then assuming there is an $\operatorname{\mathcal{A}}$ that uses $Q$ queries to solve the unitary estimation problem on this set, there exists another quantum algorithm for identifying unitaries from the set ${\mathcal{N}}$ with $O(\varepsilon Q)$ query complexity, but there is a catch: The query complexity stated is in the fractional query model, a model introduced by Cleve, Gottesman, Mosca, Somma, and Yonge-Mallo [CGM*+*09]. (In this model, it is cheaper to query a small power of a unitary compared to querying the unitary itself.) But we want to show our lower bound in the standard model assumed in Theorem 1.2. So we use a modified proof of the equivalence of the fractional and discrete query models from Berry, Childs, Cleve, Kothari, and Somma [BCC*+*17] to give an algorithm using $O(\varepsilon Q)$ queries in the usual sense, except with failure probability $\geq\exp(-50\varepsilon Q)$ instead of $\frac{2}{3}$. This extremely low success probability is not a problem, since Proposition 4.2 and Lemma 4.3 yield strong lower bounds even with small success probability. Thus any such algorithm, which used $O(\varepsilon Q)$ queries, must use $\Omega({\mathsf{d}}^{2})$ queries, which gives $Q=\Omega({\mathsf{d}}^{2}/\varepsilon)$ as desired.

We start by defining some notation for the fractional power of a reflection. Following convention from [BCC*+*17], this is the power in the negative angle direction: with this definition, $(-1)^{\alpha}=e^{-\mathrm{i}\pi\alpha}$.

Definition 4.4**.**

For a reflection $R\in{\mathsf{U}}({\mathsf{d}})$ (i.e., $R^{2}={\mathbb{I}}$) and $\alpha\in[-1,1]$, we define

[TABLE]

We can now state our main conversion lemma, which takes an algorithm that queries a fractional power of $Z$ and converts it to an algorithm that only queries $Z$ and makes fewer queries, at the cost of having low success probability.

Lemma 4.5**.**

Let $\operatorname{\mathcal{C}}\in{\mathsf{U}}(D)$ be implemented by a quantum circuit that uses $Q$ queries to $Z$ or $Z^{\dagger}$, with $Z=R^{\alpha}$ for $R\in{\mathsf{U}}({\mathsf{d}})$ a reflection and $\alpha\in(0,1]$.

Then there is a circuit $\operatorname{\mathcal{C}}_{R}\in{\mathsf{U}}(2^{Q}D)$ with $Q$ ancilla qubits, which uses $50+100\alpha Q$ queries to ${\mathrm{c}}R$ and such that, for all $\ket{\psi}\in\mathbb{C}^{D}$,

[TABLE]

A similar conversion was also stated in [BCC*+*17, Lemma 3.8], but their lemma converts the circuit $\operatorname{\mathcal{C}}$ that uses $Q$ oracle calls to $Z$ and $Z^{\dagger}$ to a circuit that uses $O(\alpha Q\frac{\log(\alpha Q/\xi)}{\log\log(\alpha Q/\xi)})$ oracle calls to ${\mathrm{c}}R$ and incurs $\ell^{2}$-norm error $\xi$ in the output state. Here, $\alpha Q$ is the fractional query complexity of the circuit, since in that model each application of $Z$ and $Z^{\dagger}$ incurs $\alpha$ cost. So their lemma is similar to ours, but incomparable: It can achieve much smaller error, but the query complexity has a log factor which would weaken our lower bound if used directly.

To eliminate this log factor, we will modify their argument and approximate the complete circuit instead of each individual segment. This gives a circuit with $O(\alpha Q)$ calls to ${\mathrm{c}}R$, but with the significantly worse success probability of $0.9\exp(-\pi\alpha Q)$. This exponentially decaying success probability is generally undesirable, but will turn out not to affect the lower bound argument, which merely needs that the success probability of the algorithm is not significantly better than guessing, which succeeds with probability $\frac{1}{N}\leq\exp(-{\mathsf{d}}^{2}/64)$.

Proof.

We closely follow the fractional query to discrete query reduction in [BCC*+*17]. We begin by writing a circuit to implement $Z$ or $Z^{\dagger}$ from an application of ${\mathrm{c}}R$. Consider the following circuit.

[TABLE]

Let us calculate the evolution of the state vector along this circuit:

[TABLE]

Using the circuit Eq. 76 with either $P^{+}_{\alpha}$ or $P^{-}_{\alpha}$, we can apply $R^{\alpha}=Z$ or $R^{-\alpha}=Z^{\dagger}$, respectively, with one use of ${\mathrm{c}}R$ and postselection of the ancilla on $\ket{0}$. This reproves [BCC*+*17, Lemma 3.3] in the slightly more general setting where we may implement $Z^{\dagger}$ and the reflection $R$ is arbitrary, not necessarily diagonal.

Let $\widehat{\operatorname{\mathcal{C}}}$ denote the circuit $\operatorname{\mathcal{C}}$ augmented so that every instance of $Z$ and $Z^{\dagger}$ is replaced with the circuit in Eq. 76 (where we suppose $\alpha>0$) as shown in Fig. 2. By composing Eq. 77, we get that

[TABLE]

where $(\bra{0}^{\otimes Q}\otimes I)\ket{\Phi^{\perp}}=0$, and the amplitude $\nu$ on the desired state satisfies

[TABLE]

Now it seems like we have made no progress, since the augmented algorithm $\widehat{\operatorname{\mathcal{C}}}$ naively uses $Q$ queries to ${\mathrm{c}}R$, which is the same number of queries used by the original algorithm. We would like to reduce this query complexity to $100\alpha Q$. If $\alpha>{\tfrac{1}{2}}$, the $\operatorname{\mathcal{C}}_{R}=\widehat{\operatorname{\mathcal{C}}}$ is what we want. So, assume $0<\alpha\leq{\tfrac{1}{2}}$.

The observation that lets us drastically reduce the query complexity of this algorithm is that the initial state of the ancilla has very low weight on $Q$-bit strings with a large number of $1$’s. In other words, it is mostly supported on low Hamming weight strings. Observe that the state is

[TABLE]

Here, we think of $\alpha$ as being small, making $1-\gamma=\Theta(\alpha)$. Indeed,

[TABLE]

If $\ket{\textup{anc}}$ were replaced by an approximate state $\ket{\widetilde{\textup{anc}}}$ that is a superposition over only the bit strings of Hamming weight $\leq K$, if $K$ is large enough (roughly $O(\alpha Q)$), then the approximate state will be very close to the original state, and hence will not affect the output too much.

Then notice what the resulting circuit looks like: It has $Q$ c$R$ gates, but the control register never holds a string of Hamming weight larger than $K$. In other words, in no branch of the superposition do we ever apply $R$ more than $K$ times, although the circuit formally has $Q$ copies of c$R$. So morally, this circuit should only require $K$ queries to $R$, not $Q$ queries.

The argument for this is given more explicitly in [BCC*+*17, Page 28], but briefly, a bit string in the superposition can be thought of as encoding a partition of the circuit, dictating how long the circuit should be run before a query from $\operatorname{\mathcal{A}}$ is responded with $R$ instead of ${\mathbb{I}}$; in this way, one can run the circuit with $K$ oracle queries, by conditioning the intermediate unitaries on the ancilla bit string.

Let $\ket{\textup{anc}_{\leq K}}$ denote the subnormalized state of $\ket{\textup{anc}}$ obtained by discarding all computational basis components of Hamming weight $>K$. We consider a cutoff $K=k(1-\gamma)Q$, where $(1-\gamma)Q$ is the “mean” bitstring weight, and $k$ is a parameter to be chosen. Then

[TABLE]

where we used a Chernoff bound since the error can be written as the probability that a binomial random variable takes a large value. With this, we can compute the $\ell^{2}$ error incurred by replacing $\ket{\textup{anc}}$ with a normalized state $\ket{\textup{anc}_{\leq K}}/\lVert\ket{\textup{anc}_{\leq K}}\rVert_{2}$.

[TABLE]

We take $K=\lceil 40+40(1-\gamma)Q\rceil$, so $k\geq 40+\frac{40}{(1-\gamma)Q}$, and

[TABLE]

So, let $\widehat{\operatorname{\mathcal{C}}}_{\leq K}$ denote the circuit of $\widehat{\operatorname{\mathcal{C}}}$, modified so that upon being given $\ket{0}^{\otimes Q}$, the circuit prepares the truncated $\ket{\textup{anc}_{\leq K}}/\lVert\ket{\textup{anc}_{\leq K}}\rVert_{2}$ instead of $\ket{\textup{anc}}$. This will be our final choice of circuit $\operatorname{\mathcal{C}}_{R}$, so our goal now is to establish the claimed properties of the resulting circuit

[TABLE]

where $\ket{\Phi^{\perp}}$ is some state satisfying $(\bra{0}^{\otimes Q}\otimes I)\ket{\Phi^{\perp}}=0$.

There is some phase ambiguity in definitions, since $\ket{\widetilde{\mathrm{out}}}$ could have any phase as long as $\nu_{\leq K}$ compensates for it so that the product of the phases is correct. So we will choose $\nu_{\leq K}$ to have the same phase as $\nu$, so that $\lvert\nu_{\leq K}-\nu\rvert=\lvert\lvert\nu_{\leq K}\rvert-\lvert\nu\rvert\rvert$. If the ancilla state changes by $\xi$ in $\ell^{2}$ distance, since the rest of the circuit is unitary, this changes the output state of $\widehat{\operatorname{\mathcal{C}}}$ by $\xi$ error. Using this, we can conclude

[TABLE]

A consequence of this is that $\lvert\nu_{\leq K}\rvert\geq\exp(-\alpha\pi\tfrac{Q}{2})-0.01\exp(-100\alpha\pi Q)\geq 0.99\exp(-\alpha\pi\tfrac{Q}{2})$. Similarly, we can show

[TABLE]

This gives the desired properties. The number of queries made is $K=\lceil 40(1+(1-\gamma)Q)\rceil\leq 50+100\alpha Q$. ∎

We are now ready to prove the final lower bound, which we restate for the reader’s convenience:

See 1.2

Proof.

Consider a quantum circuit $\operatorname{\mathcal{A}}$ with the assumed properties and that uses $Q$ queries to the oracle, where $Q$ is a function of ${\mathsf{d}}$ and $\varepsilon$. Consider the net of reflections ${\mathcal{N}}\subset{\mathsf{U}}({\mathsf{d}})$ as described in Proposition 4.1. We let $Z=R^{\alpha}$. If we run a process tomography algorithm $\operatorname{\mathcal{A}}$ with the properties assumed in the statement on $Z$, then $\operatorname{\mathcal{A}}$ is promised to output an estimate ${\bm{U}}$ such that $\operatorname{dist}({\mathsf{U}}(1){\bm{U}},Z)<\varepsilon$ with probability $\geq 2/3$ using $Q$ queries of ${\mathrm{c}}Z$ and ${\mathrm{c}}Z^{\dagger}$. Since $({\mathrm{c}}R)^{\alpha}={\mathrm{c}}(R^{\alpha})={\mathrm{c}}Z$, these ${\mathrm{c}}Z$ and ${\mathrm{c}}Z^{\dagger}$ queries can be thought of as fractional queries to ${\mathrm{c}}R$. We first show that the assumed output guarantee, run on an element $R\in{\mathcal{N}}$ of the net, suffices to determine which element $R$ is, upon choosing $\alpha=\Theta(\varepsilon)$; then, we show how to perform this algorithm, only using queries to ${\mathrm{c}}{\mathrm{c}}R$, the doubly controlled $R$.

For the first step, suppose we have access to an unknown $R\in{\mathcal{N}}$, and wish to identify which element it is. We can do this by applying $\operatorname{\mathcal{A}}$ to $Z=R^{\alpha}$ for $\frac{1}{\alpha}=\lfloor\frac{1}{8\varepsilon}\rfloor$ (where we assume that $\varepsilon<\frac{1}{8}$ so that $\frac{1}{\alpha}\geq 1$). On a successful run, which occurs with probability $\geq\frac{2}{3}$, the output ${\bm{U}}$ satisfies $\lVert\mathcal{U}({\bm{U}})-\mathcal{U}(Z)\rVert_{\diamond}<\varepsilon$. In this case, we can identify which element in the net it is by computing ${\bm{U}}^{1/\alpha}$, and picking the element of the net it is closest to. The correctness guarantee implies that ${\bm{U}}^{1/\alpha}$ and $R$ are close in diamond norm:

[TABLE]

where above we used the triangle inequality and unitary invariance of the diamond norm. By the property of the net ${\mathcal{N}}$, only one element is ever within $\frac{1}{8}$ of any particular ${\bm{U}}^{1/\alpha}$, so the closest net element to it must be the correct element $R$.

Now, we show that $\operatorname{\mathcal{A}}$ can still be run given appropriate access to $R$ only. Applying Lemma 4.5 on $\operatorname{\mathcal{A}}$ with $Q$ queries to ${\mathrm{c}}Z,{\mathrm{c}}Z^{\dagger}$, we get a circuit $\operatorname{\mathcal{A}}_{R}$ that uses $50+100\alpha Q$ queries to ${\mathrm{c}}{\mathrm{c}}R$. Further, if we run $\operatorname{\mathcal{A}}_{R}$, the ancilla qubits will be measured in $\ket{0}^{\otimes Q}$ with probability $\geq 0.9\exp(-\alpha\pi Q)$, and the resulting pure state will have $\ell^{2}$-norm error $\exp(-99\alpha\pi Q-20)$ from the true output of $\operatorname{\mathcal{A}}$. We will perform $\operatorname{\mathcal{A}}_{R}$, measure all of the ancilla in the computational basis, and if it is $\ket{0}^{\otimes Q}$ (postselection), we output the output of $\operatorname{\mathcal{A}}_{R}$ (which we assume to be the measurement of the rest of the qubits in the computational basis); otherwise, we declare failure. Since $\operatorname{\mathcal{A}}$ produces a correct estimate with probability $\geq\frac{2}{3}$, and since closeness (in $\ell^{2}$ norm) of states $\ket{a},\ket{b}$ implies closeness (in total variation distance) between the probability distributions resulting from measuring in the computational basis, $\operatorname{dist}_{\mathrm{TV}}(\mathbf{a},\,\mathbf{b})\leq 4\lVert\ket{a}-\ket{b}\rVert_{2}$ [BV97, Lemma 3.6], we have

[TABLE]

To summarize, we have an algorithm that, with $50+100\alpha Q$ queries to ${\mathrm{c}}{\mathrm{c}}R$, can output a classical description ${\bm{U}}$ of $Z=R^{\alpha}$ such that $\lVert\mathcal{U}({\bm{U}})-\mathcal{U}(Z)\rVert_{\diamond}<\varepsilon$ with probability $\geq\frac{1}{2}\exp(-\alpha\pi Q)$.

If we are promised that $R\in{\mathcal{N}}$ with ${\mathcal{N}}$ the net from Proposition 4.1, then choosing $\alpha=\lfloor\frac{1}{8\varepsilon}\rfloor^{-1}$, the output of this algorithm can distinguish an element $R$ in the net with probability $\geq\frac{1}{2}\exp(-\alpha\pi Q)$. It can also distinguish elements of the net ${\mathrm{c}}{\mathrm{c}}{\mathcal{N}}=\{{\mathrm{c}}{\mathrm{c}}R\mid R\in{\mathcal{N}}\}\subset{\mathsf{U}}(4{\mathsf{d}})$, since identifying $R$ is equivalent to identifying ${\mathrm{c}}{\mathrm{c}}R$. By Proposition 4.2 and Lemma 4.3, any algorithm that can distinguish between all the elements in ${\mathrm{c}}{\mathrm{c}}{\mathcal{N}}$ with probability $\geq\exp(\delta_{2}-\delta_{1}Q)$, where $\delta_{1}\geq 0,\delta_{2}\in\mathbb{R}$ are constants, must use $Q=\Omega({\mathsf{d}}^{2})$ queries. From this we can conclude the bound:

[TABLE]

Appendix A Deferred proofs for Section 1.1

See 1.6

Proof.

We can get Eq. 6 from minimizing the terms in the inequality Item 3.2(b) over a global phase. So, it suffices to prove Eq. 7, which relates diamond-norm distance to the other distances.

Let $\lambda_{1},\ldots,\lambda_{\mathsf{d}}$ denote the eigenvalues of $U^{\dagger}V$, and let $\sigma$ denote their spread, i.e. the length of the shortest arc (of the complex unit circle) containing the spectrum of $U^{\dagger}V$:

[TABLE]

Then, examining the geometry, we can rewrite all of the distances in terms of this spread $\sigma$.131313 This reflects that all these norms are worst-case norms, in that they in some sense maximize error over all possible directions.

For diamond-norm distance, we get a simple formula using special properties for when the two channels are isometries.

[TABLE]

To get Eq. 93, if the spread is $\geq\pi$, then the convex hull must contain the origin, and so the minimum is [math]; otherwise, the minimum is achieved at the midpoint of the extremal eigenvalues, ${\tfrac{1}{2}}(e^{\mathrm{i}a}+e^{\mathrm{i}b})$ according to the notation in Eq. 92.

We can also see

[TABLE]

Using the following basic inequalities,

[TABLE]

we can conclude from Eqs. 94 and 95 that Eq. 7 holds:

[TABLE]

See 1.9

Proof.

For ${\mathsf{d}}=1$, there is only one unitary channel (the identity channel), so all distances are zero and the equality holds trivially. So, suppose ${\mathsf{d}}\geq 2$ and, as in the proof above, let $\lambda_{1},\ldots,\lambda_{\mathsf{d}}$ denote the eigenvalues of $U^{\dagger}V$. Let $\sigma\in[0,2\pi)$ be the length of a shortest arc on the unit complex circle that covers all the eigenvalues. By rotating if necessary, we may set $\lambda_{1}=e^{\mathrm{i}\sigma/2}$ and $\lambda_{\mathsf{d}}=e^{-\mathrm{i}\sigma/2}$. Recall the following expressions.

[TABLE]

Since $\frac{1}{{\mathsf{d}}}\sum_{i}\lambda_{i}$ is a point in the convex hull, the first two lines here give the left-hand side of the proposition. For the right-hand side, if $\sigma\leq\pi$, then

[TABLE]

If $\sigma>\pi$, then ${\mathsf{d}}>2$ and there must be some eigenvalue, say, $\lambda_{j}=e^{\mathrm{i}\theta}$ where $\theta$ is in the interval $J:=[-\frac{\pi}{2}+\frac{\sigma-\pi}{2},\frac{\pi}{2}-\frac{\sigma-\pi}{2}]$, since otherwise there would be a shorter arc containing all the eigenvalues. It is easy to see by drawing a parallelogram on the complex plane that the magnitude of $e^{\mathrm{i}\sigma/2}+e^{-\mathrm{i}\sigma/2}+e^{\mathrm{i}x}$ for varying $x\in J$ assumes the largest value $1$ when $x$ is either end point of the interval $J$. So,

[TABLE]

Therefore, if we set $\sigma^{\prime}=\min(\sigma,\pi)$ in all cases, we have

[TABLE]

where the second inequality in the second line uses ${\mathsf{d}}\geq 2$ and $\cos(\sigma^{\prime}/2)\geq 0$.

The tightness of the left-hand side is seen by setting $\lambda_{1}=\cdots=\lambda_{{\mathsf{d}}/2}\neq\lambda_{{\mathsf{d}}/2+1}=\cdots=\lambda_{\mathsf{d}}$ where ${\mathsf{d}}$ is even. The tightness of the right-hand side is seen by setting $\lambda_{1}=e^{\mathrm{i}\sigma/2}$, $\lambda_{\mathsf{d}}=e^{-\mathrm{i}\sigma/2}$ with small $\sigma$ and $\lambda_{i}=1$ for all other $i$. ∎

See 1.18

The inequality follows from [YRC20, Lem. 2], but we give another short proof (with worse constants) using the idea of “Mixing Unitary” Lemma of [Has16, Cam17].

Proof of Eq. 19.

Write $Z^{\dagger}{\bm{U}}=e^{\bm{X}}$ where $\lVert{\bm{X}}\rVert_{\mathrm{op}}\leq\pi$. Since the distribution of ${\bm{X}}$ is conjugation invariant,

[TABLE]

Then, for any state $\rho$ that may be entangled with some ancilla, we see, using nested commutators $[{\bm{X}},\rho]_{k}=[{\bm{X}},[{\bm{X}},\rho]_{k-1}]$ for $k>0$ and $[{\bm{X}},\rho]_{0}=\rho$,

[TABLE]

The result follows using norm equivalence (Proposition 1.6). ∎

Appendix B Estimating eigenvalues

In this section we describe an algorithm that, given access to an unknown unitary $Z\in{\mathsf{U}}({\mathsf{d}})$, outputs an approximation to the eigenvalues of $Z$ that is close in “$L_{\infty}$-distance up to a global phase”. Let us make this precise:

Definition B.1**.**

Let $\Theta_{1},\Theta_{2}\subset[0,2\pi)$ be two finite sets of “eigenphases”. We write

[TABLE]

where $d_{\mathrm{H}}$ denotes Hausdorff distance in ${\mathbb{R}}/2\pi{\mathbb{Z}}$. In other words, $d_{\widetilde{\mathrm{H}}}(\Theta_{1},\Theta_{2})\leq\varepsilon$ iff, for some $\tau$, the sets $\Theta_{1}^{\prime}=\Theta_{1}+\tau$ and $\Theta_{2}$ have the following property (when all numbers are taken mod $2\pi$): For every $\theta_{1}\in\Theta^{\prime}_{1}$ there is $\theta_{2}\in\Theta_{2}$ with $\lvert\theta_{1}-\theta_{2}\rvert\leq\varepsilon$, and vice versa (interchanging the roles of $\Theta^{\prime}_{1}$ and $\Theta_{2}$).

Our goal is to prove the following theorem:

Theorem B.2**.**

There is an eigenvalue estimation algorithm that, given $\varepsilon>0$ and access to an unknown unitary $Z\in{\mathsf{U}}({\mathsf{d}})$, applies $Z$ at most $O({\mathsf{d}}/\varepsilon)\cdot\log^{2}{\mathsf{d}}$ times and outputs the classical description of a set $\widehat{\Theta}\subset[0,2\pi)$ such that $d_{\widetilde{\mathrm{H}}}(\Theta,\widehat{\Theta})\leq\varepsilon$ except with probability at most ${\mathsf{d}}^{-100}$, where $\Theta$ is the set of $Z$’s eigenphases. (That is, $\Theta=\{\theta\in[0,2\pi):e^{\mathrm{i}\theta}\text{ is an eigenvalue of }Z\}$.) Moreover, the quantum space requirement for the algorithm is only $2$ qudits plus $1$-qubit, and the gate complexity beyond the uses of $Z$ is only $\operatorname{poly}\log({\mathsf{d}}/\varepsilon)$.

The overall algorithm is similar to Phase Estimation. We begin as follows:

Proposition B.3**.**

There is an algorithm that, given $\varepsilon,\eta>0$, access to an unknown $Z\in{\mathsf{U}}({\mathsf{d}})$, and two eigenvectors $\ket{a}$, $\ket{b}$ with eigenphases $\alpha,\beta\in[0,2\pi)$ (respectively), has the following behavior: The algorithm uses one additional qubit, applies $Z$ at most $O(\log(1/\eta)/\varepsilon)$ times, and outputs an estimate $\theta$ that is within $\pm\varepsilon$ of $\beta-\alpha$, except with probability at most $\eta$. In addition, at the end of the algorithm still holds $\ket{a}\otimes\ket{b}$, unentangled with its results.

Proof.

The algorithm strongly resembles Quantum Phase Estimation. Given $\ket{a}\otimes\ket{b}$ in “registers $1$ and $2$”, suppose we:

• •

Adjoin a qubit $\ket{+}$ in a new “register [math]”.

• •

Performed controlled-SWAP on registers $1,2$, controlled on register [math].

• •

Apply $Z$ to register $1$.

• •

Perform controlled-SWAP again.

• •

Detach the qubit in register [math].

It is easy to calculate that this procedure leaves registers $1$ and $2$ in the state $\ket{a}\otimes\ket{b}$, with the detached qubit in state $\sqrt{1/2}(e^{\mathrm{i}\alpha}\ket{0}+e^{\mathrm{i}\beta}\ket{1})$. Up to a global phase, this is $\sqrt{1/2}(\ket{0}+e^{\mathrm{i}(\beta-\alpha)}\ket{1})$. We can then repeat this procedure but with $Z^{2}$, $Z^{4}$, $Z^{8}$, etc. in place of $Z$, yielding qubits $\sqrt{1/2}(\ket{0}+e^{\mathrm{i}2^{k}(\beta-\alpha)}\ket{1})$ for all $0\leq k<m$, at a cost of $2^{m}-1$ uses of $Z$. This is precisely the scenario arising within textbook Quantum Phase Estimation [CEMM98], prior to its QFT. If we were to finish with QFT, we would essentially complete the proof, except that it would use $m=\log(1/\varepsilon)+O(1)$ qubits to achieve error $\varepsilon$ with probability $2/3$. Instead, if we use Iterative Quantum Phase Estimation [DJSW07], we only need one additional qubit to get the same guarantee. Finally, repeating the whole algorithm $O(\log(1/\eta))$ times and taking the median result completes the proof. ∎

Proposition B.4**.**

There is an algorithm that, given $\varepsilon>0$, access to an unknown $Z\in{\mathsf{U}}({\mathsf{d}})$, and an eigenvector $\ket{a}$, has the following behavior: The algorithm uses one additional qudit and one additional qubit, applies $Z$ at most $O({\mathsf{d}}/\varepsilon)\cdot\log^{2}{\mathsf{d}}$ times, and outputs the classical description of a set $\widehat{\Theta}\subset[0,2\pi)$ such that $d_{\widetilde{\mathrm{H}}}(\Theta,\widehat{\Theta})\leq\varepsilon$ except with probability at most ${\mathsf{d}}^{-100}$, where $\Theta$ is the set of $Z$’s eigenphases.

Proof.

Write $\alpha$ for the eigenphase of $Z$ on $\ket{a}$. The algorithm repeatedly does the following: Adjoin to $\ket{a}$ a second qudit in the maximally mixed state. Then use the routine from Proposition B.3, with $\eta={\mathsf{d}}^{-C}$ for some large constant $C$. The result is that except with probability at most ${\mathsf{d}}^{-C}$, the routine uses $O(1/\varepsilon)\cdot\log d$ applications of $Z$, and ends up holding the following: $\ket{a}\otimes\ket{{\bm{b}}}$ together with an $\varepsilon$-accurate estimate of $\bm{\beta}-\alpha$, where $\ket{{\bm{b}}}$ is a uniformly random eigenvector of $Z$ and $e^{\mathrm{i}\bm{\beta}}$ is the associated eigenvalue. At this point, the difference $\bm{\beta}-\alpha$ may be recorded, and the register containing $\ket{{\bm{b}}}$ discarded.

By repeating this procedure $O({\mathsf{d}}\log{\mathsf{d}})$ times, the Coupon Collector analysis ensures that the algorithm will record $\varepsilon$-accurate values of $\theta-\alpha$ for the at most $d$ distinct values $\theta\in\Theta$. (Except with probability at most ${\mathsf{d}}^{-100}$, having chosen $C$ appropriately.) It follows that the collection $\widehat{\Theta}$ of recorded values satisfies $d_{\widetilde{\mathrm{H}}}(\Theta,\widehat{\Theta})\leq\varepsilon$ (independent of what $\alpha$ is), as required. ∎

Remark B.5**.**

If $\log(1/\varepsilon)\ll\log{\mathsf{d}}$ then the log factors in the preceding analysis may be slightly improved. In this case, since the number of possible values for $\theta-\alpha$ is only $O(1/\varepsilon)$, one can use analysis of Non-Uniform Coupon Collecting to show that only $O({\mathsf{d}}\log 1/\varepsilon)$ repetitions are required. We omit further details.

Finally, Theorem B.2 now follows by applying the algorithm from Proposition B.4 with the maximally mixed state in place of $\ket{a}$.

Appendix C Gate-efficient pure state tomography with optimal sample complexity

Although it doesn’t seem to significantly improve the gate complexity of our algorithm, we record the following result:

Theorem C.1**.**

There is a tomography algorithm for ${\mathsf{d}}$-dimensional pure states, ${\mathsf{d}}=2^{\mathsf{n}}$, with the following behavior. Given $\varepsilon>0$ and $O({\mathsf{d}}/\varepsilon)$ copies of an unknown ${\mathsf{n}}$-qubit pure state $\ket{u}\in{\mathbb{C}}^{{\mathsf{d}}}$, the algorithm sequentially measures the copies using von Neumann measurements and collects the classical results. The measurements are nonadaptively chosen using classical randomness, and each is implemented with gate complexity ${\mathsf{n}}^{O(1)}\cdot\log(1/\varepsilon)$ on the minimal number of qubits, ${\mathsf{n}}$. Finally, after the classical results are processed in $d^{O(1)}$ time, an estimate $\ket{\widehat{{\bm{u}}}}$ is output, and except with probability at most ${\mathsf{d}}^{-100}$ it satisfies $\lvert\braket{u}{\widehat{{\bm{u}}}}\rvert\geq 1-\varepsilon$.

We were unable to find this theorem in the literature, although it follows from relatively standard ideas. It would seem that all previous works on gate-efficient pure state tomography (e.g., [GLF*+*10, CL14, KRT17, GKKT20]) have sample complexity no better than $O({\mathsf{d}}\log({\mathsf{d}})/\varepsilon)$.

Proof.

Recall the standard pure state tomography algorithm from Proposition 2.2 and its analysis from Section 2.1. Slightly more strongly we use that, with $m=O({\mathsf{d}}/\varepsilon)$ copies, it can be made to satisfy ([GKKT20, Theorem 5])

[TABLE]

simultaneously for all constants $t>0$, from which

[TABLE]

easily follows (for some constant $c$). Hence for any positive even integer $k$,

[TABLE]

Recall that

[TABLE]

where ${\bm{V}}_{1},\dots,{\bm{V}}_{m}$ are independent Haar-random unitaries and ${\bm{J}}_{t}$ denotes the column index of the $t$th measurement outcome. For a given $t\in[m]$, let us define a “${\bm{V}}_{t}$-entry” to be an expression of the form $\braket{i}{{\bm{V}}_{t}}{{\bm{J}}_{t}}\braket{{\bm{J}}_{t}}{{\bm{V}}_{t}^{\dagger}}{i^{\prime}}$, where $\ket{i},\ket{i^{\prime}}$ are standard basis vectors.

Now

[TABLE]

which in turn is a sum of $({\mathsf{d}}+1+m{\mathsf{d}}^{2})^{k}$ “monomials”, each of the form

[TABLE]

for a constant $c$ with $|c|\leq 1$ and nonnegative integers $k_{t}$ summing to at most $k$. Note that the expected value of a monomial as in Equation 109 is equal to

[TABLE]

since ${\bm{V}}_{1},\dots,{\bm{V}}_{m}$ are independent.

Let us now introduce independent but pseudorandom unitaries ${\bm{V}}^{\prime}_{1},\dots,{\bm{V}}^{\prime}_{m}$. We write ${\bm{L}}^{\prime}$ and ${\bm{J}}^{\prime}_{t}$ and repeat the above development. Now suppose the following holds:

[TABLE]

Then—using that any ${\bm{V}}_{t}$-entry or ${\bm{V}}^{\prime}_{t}$-entry has magnitude at most $1$—it’s not hard to conclude that

[TABLE]

Thus for some $\eta\leq({\mathsf{d}}/\varepsilon)^{-O(k)}$ the difference is at most $\varepsilon^{k}$ and from Equation 107 we get

[TABLE]

But when $\operatorname*{Tr}((\ket{u}\!\bra{u}-{\bm{L}}^{\prime})^{k})\geq 2{\mathsf{d}}^{101}(c\sqrt{\varepsilon})^{k}$ we have

[TABLE]

provided $k=\Omega(\log d)$.

In summary, provided:

• •

${\bm{V}}^{\prime}_{1},\dots,{\bm{V}}^{\prime}_{m}$ achieve Equation 111 for $k=\Omega(\log d)$ and $\eta=O({\mathsf{d}}/\varepsilon)^{-O(k)}$,

• •

each ${\bm{V}}^{\prime}_{t}$ can be implemented with gate complexity ${\mathsf{n}}^{O(1)}\cdot\log(1/\varepsilon)$,

we get that $\lVert\ket{u}\!\bra{u}-{\bm{L}}^{\prime}\rVert_{\mathrm{op}}\leq O(\sqrt{\varepsilon})$ except with probability at most ${\mathsf{d}}^{-100}$, from which the result follows as in Section 2.1.

Let us now consider a generic product of $k_{t}$ ${\bm{V}}_{t}$-entries. Dropping the $t$ subscript for notational simplicity, it looks like

[TABLE]

Since $\Pr[{\bm{J}}=j\mid{\bm{V}}=V]=\braket{u}{V}{j}\!\braket{j}{V^{\dagger}}{u}$, the expectation of the above quantity is

[TABLE]

Since $\ket{u}$ is a unit vector, it follows that the above expectation changes by no more than $\sqrt{{\mathsf{d}}}\cdot\lambda$ in magnitude if ${\bm{V}}$ is replaced by ${\bm{V}}^{\prime}$ distributed as a quantum $({\mathsf{d}},\lambda,k+1)$-tensor-product-expander (see Definition 1 and equations (4)–(6) of [BaHH16]). By the work of [BaHH16] (see also the latest strengthening from [Haf22]), we know that a $({\mathsf{d}},\lambda,k)$-TPE can be generated by certain simple probability distribution on ${\mathsf{n}}$-qubit unitary circuits with $\mathrm{poly}({\mathsf{n}},k)\cdot\log(1/\lambda)$ gates (provided $k=o({\mathsf{d}})$). With $k=\Theta(\log{\mathsf{d}})=\Theta({\mathsf{n}})$ and $\lambda=O({\mathsf{d}}/\varepsilon)^{O(k)}$, this is ${\mathsf{n}}^{O(1)}\log(1/\varepsilon)$ gates, as needed to complete the proof of the theorem. ∎

Appendix D Compressed proof of Proposition 4.2

Here we give a summary of the proof of [BMQ22, Thm. 5]. The content in this appendix is an excerpt of results in [CDP08, CE16, HHHH10, BMQ22]. Our exposition will omit much of general discussion in those references, but present necessary pieces.

Step 1: Quantum Testers.

It is well known that a quantum channel ${\mathcal{C}}$ (unitary or not) corresponds to a Choi operator $C$, which is nothing but a density operator obtained by applying ${\mathcal{C}}$ to an unnormalized pure density matrix $\sum_{i,i^{\prime}}\ket{ii}\bra{i^{\prime}i^{\prime}}$ that is maximally entangled with an ancilla of dimension equal to that of the input of ${\mathcal{C}}$. The Choi operator of the composition of quantum channels has a formula in terms of the Choi operators of the component channels, using so-called link product denoted by $*$ [CDP08]. Although the link product involves partial transposes (that generally do not preserve positivity) whenever an output of a channel is plugged in to another channel, the result of a link product is always a positive semidefinite (PSD) operator.

Let us call a sequence of quantum channels a quantum network, which is not necessarily a unitary. It has been observed [CDP08] that even if some part of a quantum network is unspecified, one can associate a Choi operator to a quantum network by introducing a Hilbert space for each of input and output of an unspecified channel, giving a blank slot that will be filled by a quantum channel. Given a Choi operator $T_{x}$ for a quantum network with $Q$ slots unfilled and with an outcome $x$ postselected, if we have $Q$ Choi operators $C_{1},C_{2},\ldots,C_{Q}$ of quantum channels to fill the $Q$ slots, we can express the probability of obtaining $x$ as

[TABLE]

where $C=C_{1}\otimes C_{2}\otimes\cdots\otimes C_{Q}$ and $*$ is the link product [CE16, §2.4]. This encompasses the scenario where the $Q$ channels are not parallel but causally ordered since $T_{x}$ can include e.g. identity channels in between the slots. Here, $C^{T}$ is the full, rather than a partial, transpose of $C$. Normally, a link product of two Choi operators takes partial transpose on the input Hilbert space of the causally succeeding factor, but since each channel $C_{j}$ takes input from an open “leg” of $T_{x}$ and some other leg of $T_{x}$ takes input from the output of $C_{j}$, we equivalently take full transpose of $C_{j}$ while keeping $T_{x}$ intact. We may think of this as a consequence of the fact that $T_{x}*C$ is a one-dimensional operator, a number. Note that $C_{j}^{T}$ is a Choi operator of a quantum channel; in particular, if $C_{j}$ is that of a unitary channel, so is $C_{j}^{T}$. The collection $\{T_{x}\}_{x}$ is called a quantum tester [CDP09].

Step 2: Semidefinite programming.

Now we consider the problem in the statement of Proposition 4.2. There are $N$ candidate channels, labeled by $x=1,\ldots,N$. We are given access to $Q$ uses of one of the candidates. Letting a quantum network to output $y$, we have a tester $\{T_{y}~{}|~{}y=1,2,\ldots,N\}$ that has $Q$ blank slots to which we plug in $Q$ identical channels, chosen uniformly at random from the candidates. Given $Q$ uses of a candidate channel labeled by $x$, the probability of outputting $y$ is $\operatorname*{Tr}(T_{y}(V_{x}^{\otimes Q})^{T})$ by Eq. 117 where $V_{x}$ is the corresponding Choi operator. The average probability of outputting a correct label is given by

[TABLE]

Now, each operator $T_{x}$ is a link product of some PSD operators, and hence is PSD. In addition, since probabilities must add up to $1$, we have $\operatorname*{Tr}(\sum_{y}T_{y}(V_{x}^{\otimes Q})^{T})=1$ for each $x$.

Actually, the blank slots of a tester may be filled with arbitrary channels and must still produce a probability distribution. In particular, we must have

[TABLE]

for every $V$ that is the Choi operator of a unitary channel. Here we do not need to take transpose because the transpose of the Choi operator $V(U)$ of a unitary channel $\mathcal{U}(U)$ is still the Choi operator $V(U^{*})$ of the complex conjugate unitary channel $\mathcal{U}(U^{*})$. This condition consists of infinitely many equations; however, one observes [CE16] that this condition is equivalent to demanding $\operatorname*{Tr}(T_{\text{all}}W)=1$ for any affine combination $W$ in the affine span

[TABLE]

This seemingly more complicated condition simplifies the situation because any affine space in a finite dimensional vector space has a finite basis. So, the condition transcribes to a finite number of equations $\operatorname*{Tr}(T_{\text{all}}W_{j})=1$, one for each affine basis element $W_{j}$ of ${\mathcal{W}}$.

A set of PSD operators $T_{y}$ with the affine constraint $\operatorname*{Tr}(T_{\text{all}}W_{j})=1$ may go beyond what a physical quantum network gives and may not qualify as a tester; we do not claim that physically realizable testers are fully characterized by ${\mathcal{W}}$. Nonetheless, we consider a semidefinite program [CE16] specified by $\{V_{x}~{}|~{}x=1,2,\ldots,N\}$:

[TABLE]

Since any physically realizable tester is a feasible solution, the optimal value of this program is an upper bound on $p_{succ}$. Suppose that $\{T_{x}\}_{x}$ is a feasible solution. Observe that

[TABLE]

because, if $W=\sum_{j}\xi_{j}W_{j}$ with $\sum_{j}\xi_{j}=1$, then

[TABLE]

It is an instance of the weak duality of semidefinite programming. The larger the class ${\mathcal{W}}$ is, the stronger our bound $\lambda$ will be.

Step 3: Some representation theory.

By construction, ${\mathcal{W}}$ contains

[TABLE]

where ${\mathrm{d}}U$ is the Haar measure on ${\mathsf{U}}({\mathsf{d}})$ and $V(U)$ is the Choi operator of the unitary channel $\mathcal{U}(U)$. We claim [BMQ22] that for any $T\in{\mathsf{U}}({\mathsf{d}})$,

[TABLE]

which is in the form of Eq. 122 and will therefore complete the proof of Proposition 4.2. Now, Eq. 125 is equivalent to saying [HHHH10] that there is $\gamma=\binom{Q+{\mathsf{d}}^{2}-1}{Q}>0$ with which $\lvert\braket{\phi}{\psi}\rvert^{2}\leq\gamma\operatorname*{\mathbf{E}}_{U}\bra{\phi}R(U)\ket{\psi}\bra{\psi}R(U)^{\dagger}\ket{\phi}$ for any vector $\ket{\phi}$ where $\ket{\psi}$ is a pure unnormalized state corresponding to $V(T)^{\otimes Q}$ and $R(U)=({\mathbb{I}}\otimes U)^{\otimes Q}$ is a unitary representation of a compact Lie group ${\mathsf{U}}({\mathsf{d}})$. It is standard to work in a basis where $R(U)$ is block-diagonal, so $R(U)=\bigoplus_{\mu}R_{\mu}(U)\otimes{\mathbb{I}}_{m(\mu)}$ where $\mu$ ranges over all inequivalent irreps occurring in $R$. The natural number $m(\mu)$ is the multiplicity of $\mu$ within $R$. Correspondingly we have $\ket{\psi}=\bigoplus_{\mu}\ket{\psi_{\mu}}$ and $\ket{\phi}=\bigoplus_{\mu}\ket{\phi_{\mu}}$ where each orthogonal summand, $\ket{\psi_{\mu}}$ or $\ket{\phi_{\mu}}$, can be viewed as an entangled state across the irrep $R_{\mu}$ and the “multiplicity space.” Let $d_{\mu}=\dim R_{\mu}$ and $d^{\prime}_{\mu}=\min(d_{\mu},m(\mu))$. We may write $\ket{\psi_{\mu}}=\sum_{a=1}^{d^{\prime}_{\mu}}\ket{\psi^{R}_{\mu,a}}\ket{\psi^{M}_{\mu,a}}$ where each $\ket{\psi^{R}_{\mu,a}}$ is normalized. By Cauchy–Schwarz,

[TABLE]

On the other hand, The average of conjugation by an irrep $R_{\mu}$ is to trace it out and replace it by the identity, scaled to preserve the trace (Schur’s lemma). Hence,

[TABLE]

Therefore, we have proved $V(T)^{\otimes Q}\preceq\gamma^{\prime}\bar{C}$ [HHHH10].

It remains to show that $\gamma^{\prime}\leq\binom{Q+{\mathsf{d}}^{2}-1}{Q}=\gamma$. Because $d^{\prime}_{\mu}\leq d_{\mu}=\dim R_{\mu}$ by definition, we aim to show

[TABLE]

As noted in [BMQ22], this is a formula appearing in Schur’s thesis [Sch01, Eq. (57)]. We note the following pointers to modern textbooks (e.g. [FH04, App. A]) to prove this identity. Obviously, an irrep appears in $R(U)=({\mathbb{I}}\otimes U)^{\otimes Q}$ if and only if it does in $U^{\otimes Q}$. It is well known that the irreps $R_{\mu}$ of ${\mathsf{U}}({\mathsf{d}})$ in the $Q$-fold tensor representation correspond to Young diagrams $\mu$ with $Q$ boxes and at most ${\mathsf{d}}$ rows (English notation). The character $\operatorname*{Tr}R_{\mu}(U)$ is given by the Schur polynomial $s_{\mu}(u_{1},\ldots,u_{\mathsf{d}})$ of degree $Q$ evaluated at the eigenvalues $u_{k}$ of $U$. An identity that is useful for us is the Cauchy identity

[TABLE]

where $\lambda$ ranges over all Young diagrams with at most ${\mathsf{d}}$ nonzero rows but unlimited number of boxes. Here all the inverse polynomials should be understood as a formal power series. The dimension $d_{\mu}$ is $s_{\mu}(1,1,\ldots,1)$, so if we take the degree $2n$ part of the right-hand side of Eq. 129 and evaluate it at $u_{k}=v_{\ell}=1$ for all $k,\ell$, then we obtain $\sum_{\mu}(\dim R_{\mu})^{2}$. This is equivalent to reading off the coefficient of $t^{Q}$ in the series expansion of $(1-t)^{-{\mathsf{d}}^{2}}$. Applying $(Q!)^{-1}\partial_{t}^{Q}$, the identity Eq. 128 follows.