Efficient quantum circuits for diagonal unitaries without ancillas

The Paley-ordered Walsh functions are defined on the continuous interval $0\leqslant x<1$ as [13]

$$\begin{eqnarray}&&{{w}_{j}}\left( x \right)={{\left( -1 \right)}^{\underset{n}{\overset{i=1}{\mathop{\mathop{\sum }^{}}}}\,{{j}_{i}}{{x}_{i}}}},\end{eqnarray} \tag{ 1 }$$

for integer $j=0,1,2,\ldots ,\infty$. They form a complete and orthonormal set, $\mathop{\int }^{}_{0}^{1}{{w}_{j}}\left( x \right){{w}_{l}}\left( x \right)dx={{\delta }_{jl}}$. This definition may be extended to the entire real line by periodic repetition. Here ${{j}_{i}}$ is the ith bit in the binary expansion, $j=\underset{n}{\overset{i=1}{\mathop{\mathop{\sum }^{}}}}\,\;{{j}_{i}}{{2}^{\left( i-1 \right)}}$, and ${{x}_{i}}$ is the ith bit in the dyadic expansion, $x=\underset{\infty }{\overset{i=1}{\mathop{\mathop{\sum }^{}}}}\,{{x}_{i}}/{{2}^{i}}$, where n is the index of the most significant non-zero bit of j. In standard binary notation, therefore, we have $j=\left( {{j}_{n}}{{j}_{n-1}}\ldots {{j}_{1}} \right)$ and $x=\left( {{x}_{1}}{{x}_{2}}\ldots {{x}_{n}} \right)$, where the most significant bit is on the left.

The ${{w}_{j}}$ with indices that are powers of two, $j={{2}^{n}}$, $n=1,\ldots ,\infty$ are square waves known as Rademacher functions. The Rademacher function of order n is denoted ${{R}_{n}}\left( x \right)={{\left( -1 \right)}^{{{x}_{n}}}}$. The first eight Walsh functions are plotted in figure 1. Rademacher functions are in red.

Zoom In Zoom Out Reset image size

Like trigonometric functions, Walsh functions can be used as a basis for orthonormal expansion. For discretely sampled functions this is accomplished by a discrete Walsh–Fourier transform. For arbitrary n, let us discretize the interval $\left[ 0,1 \right)$ into $N={{2}^{n}}$ points, ${{x}_{k}}=k/N$, $k=0,\ldots ,N-1$. We define discrete Walsh functions as ${{w}_{jk}}={{w}_{j}}\left( {{x}_{k}} \right)$. In terms of the bits of j, k, and x, we have

$$\begin{eqnarray}&&{{w}_{jk}}={{\left( -1 \right)}^{\underset{n}{\overset{i=1}{\mathop{\mathop{\sum }^{}}}}\,{{j}_{i}}{{k}_{i}}}}={{\left( -1 \right)}^{\underset{n}{\overset{i=1}{\mathop{\mathop{\sum }^{}}}}\,{{j}_{i}}{{x}_{i}}}},\end{eqnarray} \tag{ 2 }$$

where ${{k}_{i}}$ is the ith bit in the dyadic expansion, $k=\underset{n}{\overset{i=1}{\mathop{\mathop{\sum }^{}}}}\,{{k}_{i}}{{2}^{\left( n-i \right)}}$. The second equality shows that the functional form is the same whether x is continuous or discrete, the only difference being the number of bits in the expansion of x. This makes Walsh series useful for representing discretely sampled continuous functions.

The orthonormality and completeness properties in the discrete case are $\frac{1}{N}\underset{N-1}{\overset{k=0}{\mathop{\mathop{\sum }^{}}}}\,{{w}_{jk}}{{w}_{lk}}={{\delta }_{jl}}$ and $\frac{1}{N}\underset{N-1}{\overset{j=0}{\mathop{\mathop{\sum }^{}}}}\,{{w}_{jk}}{{w}_{jl}}={{\delta }_{kl}}$, respectively. The discrete Walsh–Fourier transform ${{a}_{j}}$ of a function ${{f}_{k}}=f\left( {{x}_{k}} \right)$ is

$$\begin{eqnarray}&&{{a}_{j}}=\frac{1}{N}\underset{N-1}{\overset{k=0}{\mathop{\mathop{\sum }^{}}}}\,{{f}_{k}}{{w}_{jk}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{f}_{k}}=\underset{N-1}{\overset{j=0}{\mathop{\mathop{\sum }^{}}}}\,{{a}_{j}}{{w}_{jk}}.\end{eqnarray} \tag{ 4 }$$

To complete the analogy with Fourier series, we recall that orthonormal functions arise as the irreducible representations of symmetry groups [17]. For trigonometric functions, the relevant group is that of translations. For Walsh functions up to order ${{2}^{n}}$, it is the group $\mathbb{Z}_{2}^{\otimes n}$, which is formed by a basis for diagonal operators on n qubits. These are the Walsh operators introduced below.

The state of an n-qubit register in a quantum computer is typically expressed as a superposition, $\left| \psi \right\rangle=\underset{N-1}{\overset{k=0}{\mathop{\mathop{\sum }^{}}}}\,{{c}_{k}}\left| k \right\rangle$, of $N={{2}^{n}}$ states in the computational basis [6], defined as

$$\begin{eqnarray}&&\left| k \right\rangle=\left| {{k}_{1}},\ldots ,{{k}_{n}} \right\rangle .\end{eqnarray} \tag{ 5 }$$

Here $k=0,1,\ldots ,N-1$ is represented as an n-bit dyadic expansion, as above, $k=\underset{n}{\overset{i=1}{\mathop{\mathop{\sum }^{}}}}\,{{k}_{i}}{{2}^{\left( n-i \right)}}$. The bits ${{k}_{i}}=0$ or 1 denote the state of the ith qubit. A unitary operator $\hat{U}\,={{e}^{i\hat{f}\,}}$ that is diagonal in this basis is given in terms of its eigenvalues as $\hat{f}\,\left| k \right\rangle ={{f}_{k}}\left| k \right\rangle$. Functions f(x) of a continuous variable, $x\in \left[ 0,L \right)$, may be represented in this way if they are discretely sampled. Here we will assume a constant sampling interval, and define the sampling (grid) points as ${{x}_{k}}=kL/N$, so that ${{f}_{k}} \equiv f\left( {{x}_{k}} \right)$. To simplify the discussion, we use units such that L = 1 to restrict the variable ${{x}_{k}}$ to the interval $\left[ 0,1 \right)$ where Walsh functions are defined. The results for general L are obtained by replacing w(x) by $w\left( x/L \right)$. We will also use the notation $\left| k \right\rangle$, $\left| {{x}_{k}} \right\rangle$, and $\left| x \right\rangle$ interchangably, dropping the subscript k on x when there is no loss of clarity.

Let ${{\hat{Z}\,}_{i}}$ denote the Pauli Z operator acting on the ith qubit, ${{\hat{Z}\,}_{i}}\left| {{k}_{1}},\ldots ,{{k}_{n}} \right\rangle={{\left( -1 \right)}^{{{k}_{i}}}}\left| {{k}_{1}},\ldots ,{{k}_{n}} \right\rangle$. We define the Walsh operator of order j on n qubits as

$$\begin{eqnarray}&&{{\hat{w}\,}_{j}}=\otimes _{i=1}^{n}{{\left( {{\hat{Z}\,}_{i}} \right)}^{{{j}_{i}}}}={{\left( {{\hat{Z}\,}_{1}} \right)}^{{{j}_{1}}}}\otimes {{\left( {{\hat{Z}\,}_{2}} \right)}^{{{j}_{2}}}}\otimes \ldots \otimes {{\left( {{\hat{Z}\,}_{n}} \right)}^{{{j}_{n}}}},\end{eqnarray} \tag{ 6 }$$

where $j=1,\ldots ,{{2}^{n}}$, and ${{j}_{i}}$ is the $i\text{th}$ bit in the binary expansion, $j=\underset{n}{\overset{i=1}{\mathop{\mathop{\sum }^{}}}}\,\;{{j}_{i}}{{2}^{\left( i-1 \right)}}$. Powers of ${{\hat{Z}\,}_{i}}$ are defined as ${{\left( {{\hat{Z}\,}_{i}} \right)}^{1}} \equiv {{\hat{Z}\,}_{i}}$ and ${{\left( {{\hat{Z}\,}_{i}} \right)}^{0}} \equiv \hat{1}\,$. The set of all Walsh operators $j=1,\ldots ,{{2}^{n}}$ forms a basis for diagonal operators on n qubits, given by all possible tensor products of single-qubit Pauli Z gates. Their eigenvalues in the computational basis $\left| x \right\rangle$, $x\in \left[ 0,1 \right)$, are Walsh functions with index j and independent variable x: ${{\hat{w}\,}_{j}}\left| x \right\rangle=\otimes _{i=1}^{n}{{\left( {{\hat{Z}\,}_{i}} \right)}^{{{j}_{i}}}}\left| k \right\rangle\;=\mathop{\prod }^{}_{i=1}^{n}{{\left( -1 \right)}^{{{j}_{i}}{{k}_{i}}}}\left| {{k}_{1}},\ldots ,{{k}_{n}} \right\rangle={{w}_{jk}}\left| k \right\rangle={{w}_{j}}\left( x \right)\left| x \right\rangle$. For the case of Rademacher functions, this relationship was pointed out by Sornborger [18], who observed that the eigenvalue of a single Pauli Z gate acting on the ith qubit in (5) is a binary-valued function of x with period $1/{{2}^{\left( i-1 \right)}}$.

The locations of the $\hat{Z}\,$ operators in ${{\hat{w}\,}_{j}}$ correspond to the positions of the 1's in the bit-reversed binary string for j. For example, the Walsh operator with j = 6 on n = 3 qubits is ${{\hat{w}\,}_{6}}=1\otimes \hat{Z}\,\otimes \hat{Z}\,$, since j = 6 in binary is $\left( {{j}_{3}}{{j}_{2}}{{j}_{1}} \right)=\left( 110 \right)$. The gate representation of ${{w}_{6}}$ is shown in figure 2. The general Walsh operator requires O(n) gates for its implementation: a single Z gate and up to $2n$ controlled NOTs.

Zoom In Zoom Out Reset image size

Using (4), any diagonal operator on n qubits may be expanded as a sum of $N={{2}^{n}}$ Walsh operators, $\hat{f}\,=\underset{N-1}{\overset{j=0}{\mathop{\mathop{\sum }^{}}}}\,{{a}_{j}}{{\hat{w}\,}_{j}}$. Walsh operators commute. Therefore any diagonal unitary may be written as a product of exponentials of Walsh operators,

$$\begin{eqnarray}&&\hat{U}\,={{e}^{i\hat{f}\,}}=\underset{N-1}{\overset{j=0}{\mathop{\mathop{\prod }^{}}}}\,{{e}^{i{{a}_{j}}{{\hat{w}\,}_{j}}}}.\end{eqnarray} \tag{ 7 }$$

Each term in the product, ${{\hat{U}\,}_{j}}={{e}^{i{{a}_{j}}{{\hat{w}\,}_{j}}}}$, is of the form $\exp \left( -i\frac{{{\theta }_{j}}}{2}{{\otimes }_{i}}{{\left( {{\hat{Z}\,}_{i}} \right)}^{{{j}_{i}}}} \right)$, where ${{\theta }_{j}}=-2{{a}_{j}}$. Hence the circuit for ${{\hat{U}\,}_{j}}$ is identical to that for ${{\hat{w}\,}_{j}}$, except the Z-gate is replaced by a Z-rotation, ${{R}_{z}}\left( -2{{a}_{j}} \right)$, where ${{R}_{z}}\left( \theta \right) \equiv {{e}^{-iZ\theta /2}}$. The circuit for $\hat{U}\,$ is given by successively applying the circuits for ${{\hat{U}\,}_{j}}$.

Figure 3 shows two equivalent ways of implementing one such term, specifically ${{\hat{U}\,}_{7}}$. As seen in this figure, the gate configuration is not unique. We adopt the convention in figure 3(b) where the CNOTs are always targeted on the highest order qubit possible. Then a precise rule for constructing the circuit for ${{\hat{U}\,}_{j}}$ can be given in terms of the binary expansion of j: a rotation gate, ${{R}_{z}}\left( -2{{a}_{j}} \right)$, is placed on the qubit corresponding to the most significant non-zero bit (MSB) of j. Then CNOTs are placed on either side, targeted on the same qubit as the rotation gate, and controlled on the qubits corresponding to the 1's other than the MSB in the binary expansion of j. This rule will be used in the next section to construct an optimal circuit for $\hat{U}\,$.

Zoom In Zoom Out Reset image size

Equation (7) easily generalizes to more than one dimension. For a d-dimensional system represented by d registers of n qubits each, the single Walsh operators will be replaced by tensor products of up to d Walsh operators over the different registers. The exact number depends on the number of variables in the function f. For applications to quantum simulation, this does not significantly increase the gate complexity as interaction potentials are generally few-body. Since products of Walsh operators are also Walsh operators, the expansions have the same form as (7) with $N={{2}^{dn}}$.

The utility of (7) is that it relates the circuit depth of $\hat{U}\,$ to the number of coefficients in the Walsh series for f. If some of these coefficients are zero or may be neglected, this leads to a reduction in the circuit depth for implementing $\hat{U}\,$. We will examine such cases below, but first we discuss methods to find optimal-depth circuits given a Walsh series for f, as well as to calculate the circuit depth based on the number of non-zero coefficients ${{a}_{j}}$.

We evaluate the time evolution of a quantum system,

$$\begin{eqnarray}&&\left| \psi \left( t \right) \right\rangle={{e}^{-i\hat{H}\,t}}\left| \psi \left( 0 \right) \right\rangle ,\end{eqnarray} \tag{ 15 }$$

using the real-space, or first quantized, representation of the wavefunction in terms of position eigenstates, $\left| \psi \left( t \right) \right\rangle = \int \left| \mathbf{x} \right\rangle\left\langle\mathbf{x}| \psi \left( t \right) \right\rangle d \mathbf{x}$. For a d-dimensional system ($d=3\;m$ for m particles), $\left| \mathbf{x} \right\rangle=\left| {{x}^{1}} \right\rangle\ldots \left| {{x}^{d}} \right\rangle$, and $d\mathbf{x} \equiv {{d}^{d}}x$. Each ${{x}^{i}}$ is discretized, and represented on the quantum computer in the computational basis as in (5). Using d registers of n qubits each, the basis states corresponding to a grid of ${{2}^{dn}}$ points can be represented.

Equation (15) is evaluated using the first-order Trotter formula [1, 3, 20, 21]. Assuming a time-independent Hamiltonian $\hat{H}\,=K\left( \hat{p}\, \right)+V\left( \hat{x}\, \right)$, where $\left[ \hat{K}\,,\hat{V}\, \right]\ne 0$, the quantum algorithm for evaluating (15) is

$$\begin{eqnarray}&&\left| \psi \left( t \right) \right\rangle={{\left( {{\hat{f}\,}^{\dagger }}\;{{e}^{-i\hat{K}\,\delta t}}\;\hat{f}\,\;{{e}^{-i\hat{V}\,\delta t}} \right)}^{t/\delta t}}\left| \psi \left( 0 \right) \right\rangle,\end{eqnarray} \tag{ 16 }$$

where $t/\delta t$ is an integer called the Trotter number. The $\hat{f}\,$ operators are quantum Fourier transforms (QFTs) and are inserted to diagonalize the kinetic energy operators $\hat{K}\,$. The potential energy $\hat{V}\,$ is already diagonal in the position representation.

It is generally not possible to evaluate the diagonal unitary kinetic and potential propagators in (16) exactly. At the very least, there will be sampling error in going from the continuous $\left| x \right\rangle$ to the discrete $\left| {{x}_{k}} \right\rangle$ representation. This contribution to the total error is in addition to the Trotter error from splitting the propagator into non-commuting parts. Letting $\hat{U}\,$ denote the operator on the right-hand side of (16), the total simulation error satisfies $E\left( \hat{U}\,,{{e}^{-i\hat{H}\,t}} \right) \equiv \| \hat{U}\,-{{e}^{-i\hat{H}\,t}} \| \leqslant \alpha t\delta t+{{E}_{G}}$, where ${{E}_{G}}$ denotes the gate error in evaluating the kinetic and potential propagators, $\alpha t\delta t$ is the first order Trotter error, and $\alpha =\|\left[ \hat{V}\,,\hat{K}\, \right]\|$ is a problem-specific constant.

As we have seen, the gate error in evaluating a diagonal unitary is equal to the absolute error in the exponent. For the potential energy propagator, approximating V(x) with a function ${{V}_{\epsilon }}\left( x \right)$ satisfying $\;{{\sup }_{x}}\;\left| {{V}_{\epsilon }}\left( x \right)-V\left( x \right) \right|\leqslant\epsilon$, results in an error $E\left( {{e}^{-i{{\hat{V}\,}_{\epsilon }}\delta t}},{{e}^{-i\hat{V}\,\delta t}} \right) \leqslant \epsilon \delta t$. Letting ${{\epsilon }_{V}}$ be the error in V(x) and ${{\epsilon }_{K}}$ be the error in K(p), the total gate error for the algorithm satisfies ${{E}_{G}}\leqslant {{\epsilon }_{V}}t+{{\epsilon }_{K}}t$. The total error in evaluating $\hat{U}\,$ therefore satisfies

$$\begin{eqnarray}&&E\left( \hat{U}\,,{{e}^{-i\hat{H}\,t}} \right)\leqslant \alpha t\delta t+{{\epsilon }_{V}}t+{{\epsilon }_{K}}t.\end{eqnarray} \tag{ 17 }$$

Since the diagonal unitaries can be implemented efficiently, and the QFT requires $\text{poly}\left( n \right)$ gates, the entire algorithm is efficient, and requires $O\left( \text{poly}\left( n \right),1/{{\epsilon }_{V}},1/{{\epsilon }_{K}},t/\delta t \right)$ gates.

The parameters $\delta t$, ${{\epsilon }_{V}}$, and ${{\epsilon }_{K}}$ may be varied to obtain the shortest gate sequence for the simulation given a combined total error tolerance. Here, we only consider the problem of finding the shortest gate sequence for a single Trotter step. This corresponds to finding the shortest Walsh series for the approximate potential and kinetic energies given ${{\epsilon }_{V}}$ and ${{\epsilon }_{K}}$.

The Eckart barrier is defined as $A\;\text{sech}\left( a\;x \right)$ [16], and is plotted in figure 6 for A = 1, a = 0.05. Also shown is a plot of a 19-term Walsh series for this potential that is accurate to 10%. This series was constructed by including a subset of coefficients from the full Walsh–Fourier series starting from the largest, then the next largest, etc until the function was reproduced within the required 10% accuracy. We used a ${{2}^{n}}$-term Walsh–Fourier transform with n = 13 to approximate the infinite series. This introduces a discretization error of about $0.1%$. The 19 largest coefficients included in the approximate Walsh series are those with Paley indices: 1, 2, 4, 7, 8, 11, 13, 14, 16, 19, 21, 22, 25, 32, 35, 37, 38, 64, and 67. The smallest power of 2 that is greater or equal to every index is ${{2}^{7}}=128$, which means that 7 qubits are necessary to represent the series.

Zoom In Zoom Out Reset image size

This approach gives the minimal set of Walsh–Fourier coefficients, and is usually very close to the fully optimized solution found when the magnitudes of the coefficients are allowed to vary [14]. We find that only 7 qubits are necessary to represent the potential to 10% accuracy, with the given set of parameters. If $n>7$, only the qubits corresponding to the seven most significant digits in the register will be used. This illustrates the resource savings possible if n is large. Although useful for illustrating the approach, the classical algorithm we described for finding best subset of Walsh–Fourier coefficients to approximate a function is not efficient since it requires first calculating a high-dimensional Walsh–Fourier transform. In fact it is not necessary to do this. For a given k, efficient methods exist for finding the best k-term Walsh–Fourier series approximation to a given function without calculating the entire transform [22].

Had we opted to use a partial Walsh–Fourier series (keeping all ${{2}^{n}}$ coefficients for some integer n) to approximate the Eckart barrier, we would also find that $n7$ is required to obtain better than 10% accuracy. (The discretization error with n = 7 is $7.8%$, and with n = 6 is $15.6%$.) The efficient circuit produced in this way requires a total of ${{2}^{7}}-3=125$ gates, of which ${{2}^{6}}=64$ are rotation gates. (The Eckart barrier is an even function. Therefore half the Walsh coefficients are zero.) In contrast, the truncation described in the previous paragraph gives a circuit with a total of approximately 50 gates, of which 19 are rotation gates. This is more than a factor of two improvement.

We performed numerical simulations of equation (16) for the Eckart barrier with multiple error tolerances on the potential. The wavefunction was initialized to a Gaussian wavepacket traveling towards the barrier. Since there is a known polynomial-time algorithm for the kinetic energy propagator, ${{e}^{i{{\hat{p}\,}^{2}}/2}}$ [23], we evaluated it with maximum resolution. The time-evolution of the wavefunction is shown in figure 7. One can see from these figures that relatively few Walsh functions are needed for an accurate simulation. Even the lowest fidelity simulation reproduces the important features of the quantum scattering problem, including interference fringes. (Although here the fringes are due to periodic boundary conditions and are not physical.)

Zoom In Zoom Out Reset image size

For the present example, (17) drastically overestimates the total error, since it is a bound over all wavefunctions. To quantify the error in the simulation for the particular initial states under consideration, we found it more convenient to use the fidelity, defined as $F=\left|\left\langle\psi\left(t\right)|{{\psi }_{0}}\left( t \right)\right\rangle \right|$, where $\left| {{\psi }_{0}}\left( t \right) \right\rangle$ is the exact final wavefunction. The fidelity is related to the simulation error defined previously as $E={{\sup }_{\left| \psi \right\rangle}}\sqrt{2\left( 1-F \right)}$. As a proxy for $\left| {{\psi }_{0}}\left( t \right) \right\rangle$, we used a 10-qubit simulation with maximum possible resolution (including all Walsh operators) and 1000 time steps. By numerically analyzing the scaling of the error with the number of time steps, we found this number of time steps gives a Trotter error of less than 1%.