Probabilistic unitary synthesis with optimal accuracy

In quantum simulation and quantum computation, a global unitary transformation on a many-body quantum system is obtained as a sequence of unitary transformations on a fixed-size system, e.g., those obtained by nearest-neighbor interactions. To guarantee and increase the accuracy of obtaining such transformations, rather than controlling their continuous parameters, each unitary transformation on the fixed-size system is realized as a sequence of gates chosen from a finite gate set $\{g_{i}\}_{i}$, where each $g_{i}$ results in a fixed unitary transformation with negligible error thanks to the sophisticated calibration, quantum error correction (Terhal, 2015) or the nature of the system (Kitaev, 2003). If $\{g_{i}\}_{i}$ is universal, arbitrary unitary transformation can be approximated by a unitary transformation $g_{i_{n}}\circ\cdots\circ g_{i_{2}}\circ g_{i_{1}}$ obtained as a gate sequence for an appropriate choice of gate length $n$ depending on the approximate error one wants to achieve. For a given universal gate set such as the set of the Hadamard, controlled-NOT, and $\pi/8$ gates (Nielsen and Chuang, 2000), an algorithm to find a gate sequence for a given unitary transformation and an approximation error bound is called a unitary synthesis algorithm.

To suppress the effect of decoherence or overhead caused by the fault-tolerant implementation of each gate (Knill et al., 1998; Aharonov and Ben-Or, 2008), various studies (Kitaev et al., 2002; Harrow et al., 2002; Kliuchnikov et al., 2016, 2013; Ross, 2015; Bocharov et al., 2015; Fowler, 2011; Bouland and Giurgica-Tiron, 2021) have proposed unitary synthesis algorithms for minimizing the length of the output gate sequence. Following the celebrated Solovay-Kitaev algorithm (Kitaev et al., 2002), many algorithms are used to find one of the shortest gate sequences that can approximate a target unitary transformation $\Upsilon$ within the desired approximation error. Obviously, the goal can be achieved by brute force search (Fowler, 2011). However, to guarantee their efficiency, many algorithms are designed for synthesizing restricted classes of unitary transformations by using particular gate sets or for finding a nearly shortest gate sequence.

While approximating an $\Upsilon$ by using a single sequence of gates is a natural approach, the advantage of another approach using a probabilistic mixture of unitaries has been demonstrated (Hastings, 2017; Campbell, 2017; Kliuchnikov et al., 2022). Suppose that a synthesis algorithm produces a gate sequence for implementing a unitary transformation in $\{\Upsilon_{\vec{i}}\}=\{g_{i_{n}}\circ\cdots\circ g_{i_{2}}\circ g_{i_{1}}\}_{\vec{i}}$ in accordance with the probability distribution $p(\vec{i})$ to approximate an $\Upsilon$. If the algorithm independently samples $\vec{i}$ for each time the $\Upsilon$ is used in the entire circuit, the physical transformation governed by the randomly executed unitary transformation $\Upsilon_{\vec{i}}$ in accordance with the $p(\vec{i})$ is described by a probabilistic mixture $\sum_{\vec{i}}p(\vec{i})\Upsilon_{\vec{i}}$ of unitaries. In this case, the approximation error should be measured by the distance between the $\Upsilon$ and $\sum_{\vec{i}}p(\vec{i})\Upsilon_{\vec{i}}$.

Campbell (Campbell, 2017) and Vadym et al. (Kliuchnikov et al., 2022) constructed algorithms to compute a probability distribution $\{p(\vec{i})\}_{\vec{i}}$ for a given $\Upsilon$ and a set $\{\Upsilon_{\vec{i}}\}_{\vec{i}}$ of unitary transformations implemented as a gate sequence such that the approximation error of $\sum_{\vec{i}}p(\vec{i})\Upsilon_{\vec{i}}$ against $\Upsilon$ is almost quadratically better than that of a single optimal unitary transformation in $\{\Upsilon_{\vec{i}}\}_{\vec{i}}$. More precisely, $\left\|\Upsilon-\sum_{\vec{i}}p(\vec{i})\Upsilon_{\vec{i}}\right\|_{\diamond}=O(\epsilon^{2})$ for the worst approximation error $\epsilon=\max_{\Upsilon}\min_{\vec{i}}\frac{1}{2}\left\|\Upsilon-\Upsilon_{\vec{i}}\right\|_{\diamond}$ caused by deterministic synthesis, where $\left\|\mathcal{A}-\mathcal{B}\right\|_{\diamond}$ is the diamond norm (Watrous, 2018; Kitaev et al., 2002). This also indicates that probabilistically executing $\Upsilon_{\vec{i}}$ in accordance with $p(\vec{i})$ can further reduce the length of the shortest gate sequence without increasing the approximation error (if one measures the error by using the above diamond norm) (Campbell, 2017). In general, probabilistic synthesis consists of two procedures; (i) computing a probability distribution $\{p(\vec{i})\}_{\vec{i}}$ and sampling a description $\vec{i}$ of a gate sequence with a classical computer, and (ii) implementing the sampled gate sequence on a quantum computer. In contrast to deterministic synthesis, procedure (ii) may require a quantum computer to rearrange a gate sequence each time it realizes a target unitary transformation. However, such rearrangeability is usually assumed as a standard functionality of a quantum computer.

However, procedure (i) should be meticulously designed to construct a practical synthesis algorithm. This is because the number of possible gate sequences grows exponentially with respect to the length $n$ of the sequence, resulting in a large degree of freedom in choosing $\{p(\vec{i})\}_{\vec{i}}$. While a probabilistic synthesis algorithm with guaranteed time complexity exists for single-qubit unitary transformations that correspond to axial rotations (Kliuchnikov et al., 2022), no such algorithm was known even for general single-qubit unitary transformations. Furthermore, a fundamental question remained open regarding the optimality of existing synthesis algorithms in comparison to the minimum approximation error $\min_{p}\left\|\Upsilon-\sum_{\vec{i}}p(\vec{i})\Upsilon_{\vec{i}}\right\|_{\diamond}$. Minimax optimization makes it difficult to investigate the minimum approximation error from an analytical perspective except for a few specific $\Upsilon$ and sets $\{\Upsilon_{\vec{i}}\}_{\vec{i}}$ (Sacchi and Sacchi, 2017).

In this section, we summarize basic notations used throughout the paper. Note that we consider only finite-dimensional Hilbert spaces. In particular, two-dimensional Hilbert space $\mathbb{C}^{2}$ is called a qubit. The $\mathbf{L}\left(\mathcal{H}\right)$ and $\mathbf{Pos}\left(\mathcal{H}\right)$ represent the set of linear operators and positive semidefinite operators on Hilbert space $\mathcal{H}$, respectively. $\mathbb{I}\in\mathbf{Pos}\left(\mathcal{H}\right)$ represents the identity operator, and we sometimes use the subscript to specify the system where $\mathbb{I}$ acts as $\mathbb{I}_{\mathcal{H}}$. For Hermitian operators $A$ and $B$ on $\mathcal{H}$, $A\geq B$ represents $A-B\in\mathbf{Pos}\left(\mathcal{H}\right)$, and $A>B$ represents $A-B$ is positive definite. The $\mathbf{S}\left(\mathcal{H}\right):=\left\{\rho\in\mathbf{Pos}\left(\mathcal{H}\right):\text{tr}\left[\rho\right]=1\right\}$ and $\mathbf{P}\left(\mathcal{H}\right):=\left\{\rho\in\mathbf{S}\left(\mathcal{H}\right):\text{tr}\left[\rho^{2}\right]=1\right\}$ represent the set of quantum states and that of pure states, respectively. Pure state $\phi\in\mathbf{P}\left(\mathcal{H}\right)$ is sometimes alternatively represented by complex unit vector $|{\phi}\rangle\in\mathcal{H}$ satisfying $\phi=|{\phi}\rangle\langle{\phi}|$. Any physical transformation of the quantum state can be represented by a completely positive and trace preserving (CPTP) linear mapping $\Gamma:\mathbf{L}\left(\mathcal{H}_{1}\right)\rightarrow\mathbf{L}\left(\mathcal{H}_{2}\right)$. There exists one-to-one correspondence between a linear mapping $\Xi:\mathbf{L}\left(\mathcal{H}_{1}\right)\rightarrow\mathbf{L}\left(\mathcal{H}_{2}\right)$ and its Choi-Jamiołkowski operator $J(\Xi):=\sum_{i,j}|{i}\rangle\langle{j}|\otimes\Xi(|{i}\rangle\langle{j}|)\in\mathbf{L}\left(\mathcal{H}_{1}\otimes\mathcal{H}_{2}\right)$.

The trace distance $\left\|\rho-\sigma\right\|_{\text{tr}}$ of two quantum states $\rho,\sigma\in\mathbf{S}\left(\mathcal{H}\right)$ is defined as $\left\|M\right\|_{\text{tr}}:=\frac{1}{2}\text{tr}\left[\sqrt{MM^{\dagger}}\right]$ for $M\in\mathbf{L}\left(\mathcal{H}\right)$. It represents the maximum total variation distance between probability distributions obtained from measurements performed on two quantum states. A similar notion measuring the distinguishability of $\rho$ and $\sigma$ is the fidelity function, defined by $F\left(\rho,\sigma\right):=\max\text{tr}\left[\Phi^{\rho}\Phi^{\sigma}\right]$, where $\Phi^{\rho}\in\mathbf{P}\left(\mathcal{H}\otimes\mathcal{H}^{\prime}\right)$ is a purification of $\rho$, i.e., $\rho=\text{tr}_{\mathcal{H}^{\prime}}\left[\Phi^{\rho}\right]$, and the maximization is taken over all the purifications. Fuchs-van de Graaf inequalities (Fuchs and van de Graaf, 1999) provide relationships between the two measures with respect to the distinguishability as follows:

holds for any state $\rho,\sigma\in\mathbf{S}\left(\mathcal{H}\right)$, where the equality of the right inequality holds when $\rho$ and $\sigma$ are pure.

The distance measuring the distinguishability of two CPTP mappings $\mathcal{A},\mathcal{B}:\mathbf{L}\left(\mathcal{H}_{1}\right)\rightarrow\mathbf{L}\left(\mathcal{H}_{2}\right)$ corresponding to the trace distance is the diamond norm $\left\|\mathcal{A}-\mathcal{B}\right\|_{\diamond}$ defined by $\frac{1}{2}\left\|\mathcal{A}-\mathcal{B}\right\|_{\diamond}:=\max_{\Phi\in\mathbf{P}\left(\mathcal{H}_{1}\otimes\mathcal{H}_{3}\right)}\left\|((\mathcal{A}-\mathcal{B})\otimes id)(\Phi)\right\|_{\text{tr}}$, where $id$ represents the identity mapping acting on $\mathcal{H}_{3}$.

Let $\Xi:\mathbf{L}\left(\mathcal{H}_{1}\right)\rightarrow\mathbf{L}\left(\mathcal{H}_{2}\right)$ be a linear Hermitian-preserving mapping and $A$ and $B$ be Hermitian operators on $\mathcal{H}_{1}$ and $\mathcal{H}_{2}$, respectively. SDP is an optimization problem formally defined with a triple $(\Xi,A,B)$ as follows (Watrous, 2018):

where $\Xi^{\dagger}:\mathbf{L}\left(\mathcal{H}_{2}\right)\rightarrow\mathbf{L}\left(\mathcal{H}_{1}\right)$ is the adjoint of $\Xi$, defined as the linear mapping satisfying $\text{tr}\left[Y^{\dagger}\Xi(X)\right]=\text{tr}\left[(\Xi^{\dagger}(Y))^{\dagger}X\right]$ for all $X\in\mathbf{L}\left(\mathcal{H}_{1}\right)$ and $Y\in\mathbf{L}\left(\mathcal{H}_{2}\right)$. We can easily verify that the solution to the primal problem is smaller than or equal to that of the dual problem. The situation when the two solutions coincide is called a strong duality. Slater’s theorem states that the strong duality holds if either of the following conditions holds:

1. (1)

   The solution to the primal problem is finite, and there exists a Hermitian operator $Y$ on $\mathcal{H}_{2}$ such that $\Xi^{\dagger}(Y)>A$.

2. (2)

   The solution to the dual problem is finite, and there exists a positive definite operator $X$ on $\mathcal{H}_{1}$ such that $\Xi(X)=B$.

For a metric space $(X,d)$ and two subsets $S,T\subseteq X$, $S$ is called an $\epsilon$-covering of $T$ if $\sup_{t\in T}\inf_{s\in S}d(s,t)\leq\epsilon$. In this article, we basically assume that $X$ is the set of CPTP mappings, the metric is defined as $d(\mathcal{A},\mathcal{B})=\frac{1}{2}\left\|\mathcal{A}-\mathcal{B}\right\|_{\diamond}$, $S$ is a finite set of unitary transformations and $T$ is a subset of unitary transformations such as a $2\epsilon$-ball $\left\{\Upsilon^{\prime}:\frac{1}{2}\left\|\Upsilon^{\prime}-\Upsilon\right\|_{\diamond}\leq 2\epsilon\right\}$ around a unitary transformation $\Upsilon:\mathbf{L}\left(\mathcal{H}\right)\rightarrow\mathbf{L}\left(\mathcal{H}\right)$.

In this section, we construct an SDP for computing the optimal probability distribution that minimizes the diamond norm between the target CPTP mapping $\mathcal{A}$ and a probabilistic mixture of CPTP mappings $\{\mathcal{B}_{x}\}_{x}$. We can compute the optimal probability distribution in probabilistic unitary synthesis by solving this SDP by restricting $\mathcal{A}$ and $\{\mathcal{B}_{x}\}_{x}$ as unitary transformations. We also mention the relationship between our SDP and the algorithm proposed by Campbell (Campbell, 2017).

For a given $\Upsilon:\mathbf{L}\left(\mathcal{H}\right)\rightarrow\mathbf{L}\left(\mathcal{H}\right)$ and a given set $\{\Upsilon_{x}:\mathbf{L}\left(\mathcal{H}\right)\rightarrow\mathbf{L}\left(\mathcal{H}\right)\}_{x\in X}$ of unitary transformations implemented as a gate sequence, which forms an $\epsilon$-covering the set of unitary transformations with sufficiently small $\epsilon$, “the convex hull finding algorithm” proposed by Campbell (Campbell, 2017) can find a probability distribution $\{\tilde{p}(x)\}_{x\in\tilde{X}}$ such that $\sum_{x\in\tilde{X}}\tilde{p}(x)H_{x}=0$, where $\Upsilon_{x}(\rho)=\Upsilon(e^{iH_{x}}\rho e^{-iH_{x}})$ and $H_{x}=O(\epsilon)$ for all $x\in\tilde{X}\subseteq X$. Note that $M=O(\epsilon)$ represents $\left\|M\right\|_{\infty}=O(\epsilon)$ as $\epsilon\rightarrow 0$ for a linear operator $M\in\mathbf{L}\left(\mathcal{H}\right)$ depending on $\epsilon$. By using the dual problem in Proposition 3.1, we can verify that the distance $\epsilon$, which is achievable by a deterministic unitary synthesis finding the closest $\Upsilon_{x}$ to approximate $\Upsilon$, can be improved into $O(\epsilon^{2})$ by mixing unitaries in accordance with the probability distribution $\{\tilde{p}(x)\}_{x\in\tilde{X}}$ as follows. First, by using the dual problem of the SDP to compute the diamond norm between two CPTP mappings, we obtain

where $\mathcal{H}^{\prime}$ represents the the output system of $\Upsilon$, which is isomorphic to $\mathcal{H}$. Second, by using the Taylor expansions $e^{iH_{x}}=\mathbb{I}+iH_{x}+R_{x}$, where $R_{x}=O(\epsilon^{2})$, we obtain

where $P=\sum_{x\in\tilde{X}}\tilde{p}(x)(H_{x}J(id)H_{x}+R_{x}J(id)R_{x}^{\dagger})\in\mathbf{Pos}\left(\mathcal{H}\otimes\mathcal{H}^{\prime}\right)$, $H_{x}$ and $R_{x}$ acts on $\mathcal{H}^{\prime}$ and we use the fact that $|{\tilde{\phi}}\rangle\langle{\tilde{\psi}}|+|{\tilde{\psi}}\rangle\langle{\tilde{\phi}}|\leq\tilde{\phi}+\tilde{\psi}$ with complex vectors $(|{\tilde{\phi}}\rangle,|{\tilde{\psi}}\rangle)=\left(\left\|R_{x}\right\|_{\infty}^{-\frac{1}{2}}\sum_{j}(\mathbb{I}_{\mathcal{H}}\otimes R_{x})|{jj}\rangle,\left\|R_{x}\right\|_{\infty}^{\frac{1}{2}}\sum_{j}|{jj}\rangle\right)$ and $(|{\tilde{\phi}}\rangle,|{\tilde{\psi}}\rangle)=\left(\left\|R_{x}\right\|_{\infty}^{\frac{1}{2}}\left\|H_{x}\right\|_{\infty}^{-\frac{1}{2}}\sum_{j}(\mathbb{I}_{\mathcal{H}}\otimes H_{x})|{jj}\rangle,i\left\|R_{x}\right\|_{\infty}^{-\frac{1}{2}}\left\|H_{x}\right\|_{\infty}^{\frac{1}{2}}\sum_{j}(\mathbb{I}_{\mathcal{H}}\otimes R_{x})|{jj}\rangle\right)$ in the inequality. Third, by letting $S$ in Eq. (8) be R.H.S. of Eq. (10), we obtain

Since the approximation error $\frac{1}{2}\left\|\Upsilon-\sum_{x\in\tilde{X}}\tilde{p}(x)\Upsilon_{x}\right\|_{\diamond}$ is generally worse than the optimal one $\min_{p}\frac{1}{2}\left\|\Upsilon-\sum_{x\in X}p(x)\Upsilon_{x}\right\|_{\diamond}$, we can obtain a better probability distribution and better estimation of the approximation error by numerically solving the SDP shown in Proposition 3.1. The ellipsoid method guarantees that $\{p(x)\}_{x\in X}$ and $r$ in the dual problem such that the difference between $r$ and the optimum dual value is less than $\epsilon$ can be computed in $poly\left(|X|\log\left(\frac{1}{\epsilon}\right)\right)$-time (Lovász, 2003). Note that we assume the dimension of the Hilbert space is constant since the unitary synthesis is usually executed for $\Upsilon$ on a fixed-size system.

This section investigates the relationship between the discrete approximation of unitary transformations and the probabilistic approximation for a general $\Upsilon$ and general set $\{\Upsilon_{x}\}_{x}$. Specifically, we show the tight relationship between $\min_{p}\left\|\Upsilon-\sum_{x}p(x)\Upsilon_{x}\right\|_{\diamond}$ and $\min_{x}\left\|\Upsilon-\Upsilon_{x}\right\|_{\diamond}$, where the former represents the minimum approximation error obtained by probabilistic synthesis and the latter represents that by deterministic synthesis when $\{\Upsilon_{x}\}_{x}$ is a set of unitary transformations implemented as a gate sequence. The first lemma shows the fundamental limitation of probabilistic synthesis, and the second one shows its superiority over deterministic synthesis.

To the best of our knowledge, the dependence of the approximation error obtained by probabilistic synthesis on the dimension of the Hilbert space shown in this theorem has never been found. This dependence is inevitable since we can also show the sharpness of this theorem in Appendix C. More precisely, we can show that for any real number $\epsilon\in(0,1]$, any integer $d\geq 2$ and any $\Upsilon$, there exists $\{\Upsilon_{x}\}_{x\in X}$ achieving the lower bound in Ineq. (14).

In the following lemma, we show the tight upper bound showing that the worst approximation error caused by deterministic synthesis can be reduced by probabilistic synthesis at least quadratically. Our upper bound slightly improves the various existing upper bounds (Hastings, 2017; Campbell, 2017), which have been proven for several classes of target unitary transformations and $\epsilon$-coverings $\{\Upsilon_{x}\}_{x\in X}$ with small $\epsilon$. Using Proposition 5.5, shown in the next section, we can verify that our upper bound is still tight even if we consider the approximation of axial single-qubit unitary transformations.