Identity Check Problem for Shallow Quantum Circuits

A fundamental task in the analysis of quantum algorithms and devices is to determine whether a given quantum circuit $U$ implements the intended unitary transformation $V$. In practice, exact implementation is rarely possible. Common sources of errors include:

1. 1.

   Hardware noise: Imperfect control and decoherence during execution.

2. 2.

   Compilation error: Approximations introduced when mapping the target circuit into a device’s native gate set.

3. 3.

   Algorithmic error: Inherent approximations in the algorithm itself. For example, Hamiltonian simulation aims to implement the time evolution $e^{iHt}$ of a quantum system with Hamiltonian $H$. A common method, Trotterization approximates this evolution by sequentially applying simpler unitary operations corresponding to terms in $H$, set by parameters such as the number of Trotter steps. Adjusting these parameters trades off between simulation accuracy and resource cost.

Efficiently estimating how close the implemented circuit $U$ is to the ideal unitary $V$ is therefore crucial both for validating quantum algorithms and for tuning algorithmic or device parameters to improve overall performance. For any unitarily invariant norm $\Vert U-V\Vert =\Vert U^{†}V-I\Vert$, so this task is equivalent to estimating the distance from $U^{†}V$ to the identity. This latter formulation is known as the identity check problem.

Unfortunately, estimating this distance between $n$-qubit unitaries is computationally difficult. Rosgen and Watrous showed [11, 10] that estimating the distance between two shallow (with depth logarithmic in $n$) quantum circuits allowing mixed states is PSPACE-hard. This essentially rules out efficient classical or quantum algorithms. Likewise, Janzing, Wocjan, and Beth established QMA-hardness of estimating the distance between two unitary circuits [5]. Ji and Wu [6] strengthened this by showing that this problem remains QMA-hard even if the circuits are constant-depth with only one-dimensional qubit connectivity. This may come as a surprise since one-dimensional shallow circuits are easy to simulate classically using Matrix Product States [17].

It is important that the no-go results stated above hold only if the distance between quantum circuits has to be estimated with a small additive error scaling inverse polynomially with the number of qubits $n$. Is it possible that some less stringent approximation of the distance can be computed efficiently? In this work, we show that the answer is YES and report linear-time classical algorithms approximating the diamond-norm and the operator-norm distances between constant-depth geometrically-local quantum circuits with a constant multiplicative error. Such approximation may be good enough for practical purposes. Note that an estimate of the distance with a constant multiplicative error is informative regardless of how small the distance is. For example, our algorithm can efficiently approximate the distance even if the latter is exponentially small in $n$. This would be impossible for an algorithm that achieves an additive error approximation scaling inverse polynomially with $n$.

Our identity check algorithm borrows many ideas from the recent breakthrough work by Huang, Liu, et al. [3] on learning shallow quantum circuits. The main ingredients of our algorithm, described below, are bounds on the diamond-norm distance $\delta (U)$ that depend on the norm of commutators between $U$ and certain observables composed of SWAP gates. These bounds and their proof are largely based on Ref. [3].

Consider $2n$ qubits labeled by integers $1,\mathrm{\dots },2n$. Let $W_i$ be the SWAP gate applied to qubits $i$ and $i+n$. Given a subset $A\subseteq [n]$, define a $2n$-qubit operator

By definition, $W_A$ acts non-trivially on $2|A|$ qubits.

The quantity $\gamma$ defined in Eq. (6) or its rescaled version $1.16\gamma$ will be the desired estimator of the distance $\delta (U)$. In the next section we show how to choose a partition $[n]=A_1\mathrm{\dots }{A}_m$ with $m=D+1$ parts such that each subset $A_j$ is a union of well-separated hypercubes of linear size $O(hD)$ and all commutators $W_{A_j}(U\otimes I)W_{A_j}(U^{†}\otimes I)$ that appear in Eq. (6) are tensor products of local commutators supported on individual hypercubes. Our construction is based on Ref. [19] which introduced so-called reclusive partitions of the $D$-dimensional Euclidean space. The key ingredient of our algorithm is an additivity lemma stated in Section 4. This lemma expresses the norm of commutators $\Vert W_{A_j}(U\otimes I)W_{A_j}(U^{†}\otimes I)-I\otimes I\Vert$ in terms of the norm of analogous local commutators supported on individual hypercubes. Each local commutator acts on a subset of at most $O{(hD)}^D$ qubits and its eigenvalues can be computed by the exact diagonalization. The additivity lemma then provides a linear time algorithm for computing the norm of global commutators $\Vert W_{A_j}(U\otimes I)W_{A_j}(U^{†}\otimes I)-I\otimes I\Vert$ which is all we need to compute the estimator $\gamma$ defined in Lemma 4.

The next lemma shows that estimation of the operator-norm distance can be reduced to estimation of diamond-norm distance given any point in the eigenvalue polygon of $U$.

In the rest of this section we prove Lemma 4 and 5.

Given a quantum circuit $U$ acting on $n$ qubits, the lightcone $ℒ(j)$ of a qubit $j\in [n]$ is defined as the set of all output qubits $i\in [n]$ that can be reached by moving through the circuit diagram forward in time starting from the input qubit $j$. For example, if $U$ is a one-dimensional circuit of depth $h$ then

For any subset of qubits $S\subseteq [n]$ let $ℒ(S)$ be the lightcone of $S$ defined as

We say that a subset of qubits $S$ is the support of an operator $O$ and write $S=\mathrm{supp}(O)$ if $O$ acts trivially on all qubits $j\notin S$. By definition,

for any operator $O$. Furthermore, $UO{U}^{†}=U_{loc}O{U}_{loc}^{†}$, where $U_{loc}$ is a “localized” circuit obtained from $U$ by removing all gates acting on qubits outside of the lightcone $ℒ(\mathrm{supp}(O))$.

Two subsets of qubits $S_1$ and $S_2$ are said to be lightcone separated if $ℒ(S_1)\cap ℒ(S_2)=\mathrm{\varnothing }$. If $O_1$ and $O_2$ are operators supported on $S_1$ and $S_2$ then $U{O}_1{O}_2{U}^{†}$ is a product of operators $U{O}_1{U}^{†}$ and $U{O}_2{U}^{†}$ with disjoint supports.

Suppose now that $n$ qubits are located at cells of a $D$-dimensional rectangular array. We shall consider partitions of the array into $D$-dimensional cubes known as reclusive partitions [19]. The linear size of each cube in the partition will be chosen as

where $h$ is the depth of $U$.

Figure 2 shows examples of 1D and 2D reclusive partitions, see Ref. [19] for the 3D example. We defer the proof of Lemma 6 to Appendix A since it is a simple rephrasing of the results established in [19]. By construction, each cube in the partition contains at most $L^D$ qubits (cubes located near the boundary of the array may be truncated) and any pair of cubes from the same set $A_j$ is lightcone separated due to Eq. (23). Write

where ${\mathrm{\ell }}_j$ is the number of cubes in $A_j$ and $A_{j,p}$ denotes the $p$-th cube in $A_j$. By constriction, we have

Since the lightcone of a cube with a linear size $L$ can be enclosed by a cube of linear size $L+2h$, the number of qubits contained in any lightcone $ℒ(A_{j,p})$ is bounded as

Here we used Eq. (23).

Consider the diamond-norm distance $\delta (U)$ and specialize the commutator-based bound of Lemma 4 to the reclusive partition $[n]=A_1\mathrm{\dots }{A}_{D+1}$. By definition,

Lightcone separation of cubes $A_{j,p}$ implies that operators $(U\otimes I)W_{A_{j,p}}(U^{†}\otimes I)$ acts on pairwise disjoint subsets of qubits. Thus

where we defined commutators

The above shows that $K_{j,p}$ are operators acting on pairwise disjoint subsets of qubits (for a fixed $j$). Let $U_{j,p}$ be a “localized” circuit obtained from $U$ by replacing all gates acting on at least one qubit outside of the lightcone $ℒ(A_{j,p})$ with the identity. Then $U_{j,p}$ acts non-trivially only on the lightcone $ℒ(A_{j,p})$ and

The support of $K_{j,p}$ includes all qubits in the left $n$-qubit register contained in $ℒ(A_{j,p})$ as well as all qubits in the right $n$-qubit register contained in $A_{j,p}$. Thus

Eigenvalues of a unitary operator acting on $m$ qubits can be computed in time $O(2^{3m})$ by the exact diagonalization of a unitary $2^m\times 2^m$ matrix. Thus one can compute all eigenvalues of the commutator $K_{j,p}$ in time

In the next section we show that the norm

that appears in the bound of Lemma 4 is a simple function of eigenvalues of individual commutators $K_{j,p}$.

In this section we show how to compute the norm of commutators that appear in Lemma 4. First, let us introduce some terminology. Let $S^1=\{z\in ℂℂ:|z|=1\}$ be the unit circle. If $U$ is a unitary operator, let $\mathrm{𝖾𝗂𝗀}(U)\subseteq S^1$ be the set of eigenvalues of $U$ (ignoring multiplicities). Consider $2n$ qubits, a subset $A\subseteq [n]$, and a SWAP operator $W_A={\prod }_{i\in A}{W}_i$ where $W_i$ is the SWAP gate acting on qubits $i$ and $i+n$. Consider a commutator

We claim that $\mathrm{𝖾𝗂𝗀}(K_A)=\mathrm{𝖾𝗂𝗀}(K_A^{†})$. Indeed, $K_A^{†}=W_A{K}_A{W}_A$. Since $W_A$ is both unitary and hermitian, conjugation by $W_A$ does not change the eigenvalue spectrum. Thus eigenvalues of $K_A$ have a form $e^{\pm i\phi }$ with $0\le \phi \le \pi$. For each $\phi$ one can choose both positive and negative sign in the exponent. Define a function $\theta$ that maps subsets of qubits $A\subseteq [n]$ to real numbers in the interval $[0,\pi ]$ such that

Note that $e^{i\theta (A)}$ is the unique eigenvalue of $K_A$ with the maximum distance from $1$ and a non-negative imaginary part. Accordingly,

We shall need the following simple fact.

Combining all above ingredients we arrive at the following algorithm for the $D$-dimensional identity check problem. We first consider the case when the input circuit $U$ is sufficiently close to the identity such that . Below we assume that a reclusive partition $[n]=A_1\mathrm{\dots }{A}_{D+1}$ of the $D$-dimensional qubit array has been already computed, see Appendix A for details. We claim that the following algorithm outputs an estimator $\gamma$ satisfying $\delta (U)\le \gamma \le (D+1)\delta (U)$.

Indeed, if line 9 is never reached, Corollary 10 of the additivity lemma imply that the output of the algorithm coincides with the quantity $\gamma$ defined in Lemma 4 specialized to the reclusive partition. In this case correctness of the algorithm follows directly from Lemma 4. Otherwise, the algorithm outputs $\gamma =2$, while Corollary 10 implies that $\delta (U)\ge \sqrt{2}$. In this case $\gamma =2$ satisfies the bounds $\delta (U)\le \gamma \le (D+1)\delta (U)$ for $D\ge 1$. We claim that the algorithm runs in time $O(n{2}^{12{(2hD)}^D})$. Indeed, the total number of cubes $A_{j,p}$ is $O(n)$. Computing the function $\theta (A_{j,p})$ at line 7 requires eigenvalues of a unitary operator $K_{A_{j,p}}$ acting on at most $4{(2hD)}^D$ qubits, as discussed in Section 3. This computation takes time $O(2^{12{(2hD)}^D})$. Hence the total runtime is $O(n{2}^{12{(2hD)}^D})$.

Next consider the general case when it is possible that $\delta (U)=2$. Define our estimator of $\delta (U)$ as $1.16\gamma$, where $\gamma$ is the output of Algorithm 1. We claim that

If the algorithm never reaches line 9 then its output coincides with the quantity $\gamma$ defined in Lemma 4 and Eq. (33) follows directly from Lemma 4, see Eq. (7). Otherwise, if the algorithm reaches line 9, it outputs $\gamma =2$ while $\delta (U)\ge \sqrt{2}$ due to Corollary 10 of the additivity lemma. In this case the first inequality in Eq. (33) follows from $\delta (U)\le 2$ and the second inequality becomes $2\le (D+1)\delta (U)$ which is true for any $D\ge 1$ since $\delta (U)\ge \sqrt{2}$. The runtime analysis is the same as before.

Since the runtime scales exponentially with the size of cubes $A_{j,p}$, one may wish to choose a partition with smaller cubes even if this negatively impacts the approximation quality. As an extreme case, one can choose each cube $A_{j,p}$ as a single qubit. However ensuring the lightcone separation between cubes in the same subset $A_j$ would require $\approx {(4h+1)}^D$ subsets $A_j$ instead of $D+1$ subsets ¹¹1Since any qubit is lightcone separated from all but at most ${(1+4h)}^D$ other qubits, Vizing’s theorem implies that qubits can be partitioned into at most $1+{(1+4h)}^D$ lightcone separated subsets.. Accordingly, the approximation ratio would become $\alpha =\mathrm{\Omega }({(4h+1)}^D)$ instead of $\alpha =D+1$.

Likewise, we expect that the runtime can be improved at the cost of a worse approximation ratio $\alpha$ by computing the norm of commutators $K_{A_{j,p}}-I$ using a randomized version of the power method [8]. It is known that this method can approximate the operator norm of a matrix of size $2^m\times 2^m$ with a multiplicative error $1+ϵ$ using $O(m/ϵ)$ matrix-vector multiplications [8]. In our case, $K_{A_{j,p}}$ is specified by a quantum circuit acting on $m=4{(2hD)}^D$ qubits with $poly(m)$ gates, see Section 3. Thus one can implement matrix-vector multiplication for the matrix $K_{A_{j,p}}-I$ in time $poly(m)2^m$. Accordingly, the power method runs in time $poly(m)2^m/ϵ$, whereas the exact diagonalization of $K_{A_{j,p}}-I$ requires time $\mathrm{\Omega }(2^{3m})$.

In this section, we implement the algorithm described in Section 5 to approximate the distance between identity and a constant-depth circuit $U$ of up to 100 qubits. We consider $U=U_1{U}_2^{†}$, where $U_1,U_2$ are two different unitaries that both approximate the time evolution $e^{-iH\tau }$ of $n$ qubits under the one-dimensional XY model:

In the limit of small $\tau$, $U=U_1{U}_2^{†}\approx I$ approximate a forward evolution followed by a backward evolution under the same Hamiltonian. Explicitly, $U_1$ and $U_2$ are the first-order Trotter approximations with, respectively, an odd-even ordering and an X-Y ordering:

We note that $X_j{X}_{j+1}$ and $Y_j{Y}_{j+1}$ are both antisymmetric under the unitary conjugation by the staggered Pauli string $X_1{Y}_2{X}_3{Y}_4\mathrm{\dots }$. Therefore, the eigenvalues of $U$ comes in complex conjugate pairs which results in a simple relationship between the diamond-norm and the operator-norm distances. Namely, a simple algebra shows that $\delta (U)=2\sin (\phi )$, where $\phi \in [0,\pi /2)$ is defined by $\Vert U-I\Vert =|e^{i\phi }-1|$. In addition, using a well-known mapping from the XY model to free fermions [15], we can compute this distance exactly, providing a benchmark for our algorithm.

In Fig. 3, we compare the exact distance $\delta (U)$ against the bounds presented in Lemma 4 for up to 100 qubits at $\tau =0.01$. For the one-dimensional qubit array, the bounds simplify to $\delta (U)\le \gamma \le 2\delta (U)$, where

Here, $A_1$ and $A_2$ are the qubit partitions illustrated in Fig. 2 with $L=4$. The lightcone separated construction of $A_j$ and the additivity lemma allow us to efficiently compute the commutator $\Vert W_{A_j}(U\otimes I)W_{A_j}(U^{†}\otimes I)-I\Vert$ for each $j$. In particular, computing the bounds reduces to finding eigenvalues of operators that are each supported on at most 12 qubits. Additionally, due to the translational invariance of the unitary $U$, only $O(1)$ such operators are unique, making the complexity of our algorithm independent of the system size.

Both bounds correctly capture the linear dependence of the Trotter error on the system size $n$, with the upper bound $\gamma$ approaching the exact $\delta (U)$ in the limit of large $n$. We note that $\Vert U-I\Vert$ and, thus, $\delta (U)$ can also be estimated by finding the maximum eigenvalue of the Hamiltionian $H_U\equiv {(U-I)}^{†}(U-I)$. Writing this Hamiltonian as a matrix product operator on a one-dimensional lattice, one can efficiently find a lower bound to its maximum eigenvalue using an algorithm based on the density matrix renormalization group (DMRG). While DMRG does not have a performance guarantee, we find that it produces lower bounds to within $3\times {10}^{-7}$ of the exact $\delta (U)$ in this example, providing a complementary approach to our algorithm in one dimension. DRMG simulations were performed using the matrix product representation library for Python $\mathrm{𝗆𝗉𝗇𝗎𝗆}$ [14] with MPS bond dimension $\chi =20$ and two DMRG sweeps in $\mathrm{𝗆𝗉𝗇𝗎𝗆}.\mathrm{𝗅𝗂𝗇𝖺𝗅𝗀}.\mathrm{𝖾𝗂𝗀}$.