Learning the quantum algorithm for state overlap

The hyperparameters are the quantum resources of the circuit. At the input, the resources are the number of ancilla qubits and data qubits that store the input data for the computation. At the output, the resources are the locations of the measurements (see figure 1). As an example, in the Swap Test for single-qubit states, we are allowed access to one ancilla qubit and two data qubits at the input, and we can measure only the ancilla qubit at the output.

The input data may be classical or quantum, depending on the computation of interest. In the case of state overlap, the input data are quantum states and hence no encoding is necessary. However, for completeness, we note that our approach also applies to classical inputs, in which case the encoding (i.e. storing the classical data in the quantum state of the data qubits) can be treated as a hyperparameter that one fixes while optimizing the algorithm.

Our approach searches for an optimal algorithm, where we consider the algorithm to be a quantum gate sequence with associated classical post-processing. We parameterize (and hence optimize over) both the gate sequence and the post-processing.

Let us first consider the gate sequence. We define a gate set ${ \mathcal A }=\{{A}_{j}(\theta )\}$. Here, each gate A_j is either a one-qubit or two-qubit gate and may also have an internal continuous parameter θ. Hence, ${ \mathcal A }$ is a discrete set, but each element of ${ \mathcal A }$ may have a continuous parameter associated with it. The precise choice of ${ \mathcal A }$ depends on which hardware one is considering. For example, the connectivity differs between IBM and Rigetti hardware, and the former employs $\mathrm{CNOT}$ gates while the latter employs controlled-Z gates. For IBM's five-qubit computer 'ibmqx4' we can write out the gate set as

$$\begin{eqnarray}\begin{array}{rcl}{{ \mathcal A }}_{\mathrm{ibmqx}4} & = & \{{{\rm{CNOT}}}^{10},{{\rm{CNOT}}}^{20},{{\rm{CNOT}}}^{21},{{\rm{CNOT}}}^{32},\\ & & {{\rm{CNOT}}}^{24},{{\rm{CNOT}}}^{34},{U}^{0}(\theta ),{U}^{1}(\theta ),{U}^{2}(\theta ),\\ & & {U}^{3}(\theta ),{U}^{4}(\theta )\},\end{array}\end{eqnarray} \tag{ 1 }$$

where U^j(θ) is an arbitrary gate on qubit j and CNOT^jk is a CNOT from control qubit j to target qubit k. Angles θ in equation (1) may be encoding multiple parameters. In this article, we treat all one-qubit gates equally in the sense that all one-qubit gates are equally complex to implement, although our approach could easily be generalized to account for different complexities for different one-qubit gates.

We consider a generic sequence of d gates

$$\begin{eqnarray}&&{G}_{\vec{k}}(\vec{\theta })={A}_{{k}_{d}}({\theta }_{d})\cdots {A}_{{k}_{2}}({\theta }_{2}){A}_{{k}_{1}}({\theta }_{1}),\end{eqnarray} \tag{ 2 }$$

where $\vec{k}=({k}_{1},\,\ldots ,\,{k}_{d})$ is the vector of indices describing which gates are employed in the gate sequence and $\vec{\theta }=({\theta }_{1},\ldots ,{\theta }_{d})$ is the vector of continuous parameters associated with these gates.

The measurement results give rise to an outcome probability vector $\vec{p}=({p}_{1},\ldots ,{p}_{l},\,\ldots )$. The desired output might be one of these probabilities p_l, or it might be some simple function of these probabilities. Hence, we allow for some simple classical post-processing of $\vec{p}$ in order to reveal the desired output. While there is enormous freedom in applying a function to $\vec{p}$, we consider a simple linear combination of probabilities:

$$\begin{eqnarray}&&g(\vec{p})=\vec{c}\cdot \vec{p}=\displaystyle \sum _{l}{c}_{l}{p}_{l},\end{eqnarray} \tag{ 3 }$$

where $\vec{c}$ is a vector of coefficients whose elements are chosen according to c_l ∈ {−1, 0, 1}. This post-processing is sufficient for the application in this paper (state overlap), although other applications may require a more general form of post-processing. Note that in our approach it is enough to consider measurements in the computational basis, as any change of the measurement basis can be incorporated into the gate sequence in equation (2). In particular, this implies that equation (3) is general enough to cover the expectation values of all Pauli product operators.

In summary, the free parameters that we optimize over (while fixing the hyperparameters) are the gate sequence vector $\vec{k}$, the continuous parameter vector $\vec{\theta }$, and the post-processing vector $\vec{c}$. For a given set of resources, these three vectors define the quantum algorithm, which we denote ${Q}_{\vec{m}}$, where $\vec{m}=(\vec{k},\vec{\theta },\vec{c})$ is the concatenated vector.

Optimizing these parameters involves defining and minimizing a cost function. The cost quantifies the discrepancy between the desired output and the actual output for a given training data set.

Suppose we want to find the algorithm that computes the function $x\to f(x)$. We generate data of the form

$$\begin{eqnarray}&&T={\{({x}^{(i)},f({x}^{(i)}))\}}_{i=1}^{2N}.\end{eqnarray} \tag{ 4 }$$

Half of this data is used for training the algorithm, i.e. optimizing the cost function. The other half is used for testing, i.e. evaluating the algorithm's performance. The training data must be sufficiently general to cover the space of possible inputs. An estimate of the amount of training data needed for state overlap is $N\approx {2}^{2{n}_{D}}$, where n_D is the number of data qubits. This can be seen by noting that our algorithm (which includes both the gate sequence and the post-processing) acts as a linear map from the data qubits' density operator space, which has dimension ${2}^{2{n}_{D}}$, to the output which is just a number and hence has dimension one. So our algorithm is basically a $1\times {2}^{2{n}_{D}}$ matrix, and an estimate of the number of constraints (and hence the number of training data points) needed to fix the algorithm's parameters is ${2}^{2{n}_{D}}$.

For example, when training the algorithm that computes overlap, x⁽ⁱ⁾ = (ρ⁽ⁱ⁾, σ⁽ⁱ⁾) consists of two quantum states ρ⁽ⁱ⁾ and σ⁽ⁱ⁾, and $f({x}^{(i)})=\mathrm{Tr}({\rho }^{(i)}{\sigma }^{(i)})$ quantifies their overlap. One can show that any algorithm that computes pure-state overlap also computes mixed-state overlap. Hence, we generate our training data by randomly choosing pure states according to the Haar measure.

Next we define a cost function. For algorithm ${Q}_{\vec{m}}$, the cost is

$$\begin{eqnarray}&&{C}_{\vec{m}}=\displaystyle \sum _{i=1}^{N}{(f({x}^{(i)})-{y}_{\vec{m}}^{(i)})}^{2}.\end{eqnarray} \tag{ 5 }$$

The cost quantifies the difference between the ideal output f(x⁽ⁱ⁾) and the actual output ${y}_{\vec{m}}^{(i)}$ for each training data point. The actual output can be written as

$$\begin{eqnarray}&&{y}_{\vec{m}}^{(i)}={y}_{(\vec{k},\vec{\theta },\vec{c})}^{(i)}=\vec{c}\cdot {\vec{p}}_{(\vec{k},\vec{\theta })}^{(i)},\end{eqnarray} \tag{ 6 }$$

where $\vec{c}$ is the post-processing vector and ${\vec{p}}_{(\vec{k},\vec{\theta })}^{(i)}$ is the outcome probability vector for input x⁽ⁱ⁾. For example, in the Swap Test, the outcome probability vector corresponds to the ancilla qubit's measurement in the Z basis. Choosing $\vec{c}=(1,-1)$ ensures that ${y}_{\vec{m}}^{(i)}$ is the expectation value of the Pauli Z operator.

For a fixed circuit gate count d, we search over the algorithm space to minimize the cost, as discussed below. We consider various d, incrementing from small to large values. When an exact algorithm exists, we typically are able to minimize the cost. That is, we can find a ${Q}_{\vec{m}}$ with ${C}_{\vec{m}}\approx 0$, for $d\geqslant {d}_{\min }$, where d_min is the minimum number of gates needed to minimize the cost (see figure 3 for example plots of final cost versus d). Note that some elements of the gate set in equation (1) commute with each other. As a consequence, there are typically many ${Q}_{\vec{m}}$ that give zero cost for $d\geqslant {d}_{\min }$. This freedom is used to simplify the algorithm at the end of the cost optimization. So, in the main results section, we present our simplest representation of such algorithms.

Zoom In Zoom Out Reset image size

The cost in equation (5) is a function of several parameters that can be divided into two groups: discrete and continuous. Discrete parameters are those which describe the circuit topology and post-processing of the algorithm. These are the gate sequence vector $\vec{k}$ and the post-processing vector $\vec{c}$. The angles $\vec{\theta }$ are treated as continuous parameters. They define all gates that depend on a parameter. For IBM and Rigetti architectures considered here, angles $\vec{\theta }$ specify all one-qubit gates present in the algorithm. Only the total number of gates d is fixed during optimization, which means that while the length of $\vec{k}$ does not change, the number of elements of $\vec{\theta }$ may vary as the optimization proceeds.

The optimization is performed in iterations until the cost reaches a (possibly local) minimum. Figure 4 shows a schematic description of a single iteration of the optimization algorithm. Each iteration begins with an attempt to modify $\vec{k}$ and $\vec{c}$. While modifying $\vec{k}$, we consider random updates that may involve an arbitrary number of gates. However, updates affecting a smaller number of gates are more probable. In this process, we may change the position or support of a given gate or change its type, e.g. from a one-qubit gate to a CNOT. The update is constrained to result in an algorithm that cannot be easily shortened. For example, the gate sequence that involves two one-qubit gates next to each other is not allowed, since those gates can be combined to a single one-qubit gate. This is a desired feature, as we optimize with a fixed total number of gates. Similarly, we randomly modify $\vec{c}$ giving more preference to changes affecting fewer measurements.

Zoom In Zoom Out Reset image size

Every change in $\vec{k}$ or $\vec{c}$ is followed by reoptimization of continuous parameters $\vec{\theta }$. This is an important step, as changing the gate sequence or post-processing function alone, without reoptimizing the gates' internal parameters $\vec{\theta }$, will most likely cause the cost to increase significantly, effectively suppressing any update of $\vec{k}$ or $\vec{c}$. The continuous part of the optimization is done in a sweeping fashion in which all one-qubit gates are updated sequentially. In this approach, at a given time, a single one-qubit gate is updated while all remaining gates are fixed. After the best one-qubit gate (the one that minimizes the cost) is identified, the optimization algorithm moves to the next one-qubit gate. We allow for randomly changing the order of updating one-qubit gates as a means to avoid local minima. We use a steepest descent method to optimize single one-qubit gates. Note that an arbitrary one-qubit gate can be described (up to a global phase, that does not affect the algorithm) by three real parameters. That is, the steepest descent method mentioned above operates in three-dimensional space. The continuous part of the optimization is repeated, until convergence of the cost function is achieved.

Once the continuous optimization has converged, we compare the final cost C in a given iteration with the current best one C_best. If the cost C is lower than the current best, the new discrete parameters $\vec{k}$ and $\vec{c}$ are accepted. If it is larger, the change is accepted with probability exponentially decreasing in the difference C − C_best following the simulated annealing method.

Every few iterations we check whether the current gate sequence ${G}_{\vec{k}}$ can be compressed. This goes well beyond the simple checks following the update of vector $\vec{k}$ described above. Here, we are trying to find a subsequence of ${G}_{\vec{k}}$ that can be nontrivially rewritten using the same or a smaller number of gates. If such a subsequence is found, we modify ${G}_{\vec{k}}$ accordingly, as this may lead to shortening the gate sequence without increasing the cost. Since the total number of gates is fixed, such compression results in the ability to add gates to the sequence. If that is the case, we insert one-qubit identity gates and reoptimize their continuous parameters as described above. To check if a given subsequence can be rewritten we recursively use the same approach that we use for the full algorithm, which is essentially described in figure 4 except in this case we do not consider the post-processing vector $\vec{c}$.

We remark that the cost function may be difficult to optimize primarily due to many low lying local minima. Thus, it is important to develop techniques to increase the chances of avoiding them. We found it particularly useful to compress the gate sequence periodically, as random updates to vectors $\vec{k}$ and $\vec{c}$ tend to produce local minima that usually include redundant subsequences. As described above, we have developed automated tools to remove such subsequences, which usually allows us to escape local minima.

Let us now discuss the scaling of the approach described above. The optimization requires the cost to be evaluated multiple times during every iteration. As part of computing the cost, one has to evaluate ${y}_{\vec{m}}^{(i)}$ in equation (5) for each training data point, which necessarily scales exponentially with the number of qubits on a classical computer. However, it can be outsourced to a quantum computer. Such a hybrid algorithm will efficiently compute the contribution to the cost from a single element of a training data set, although the resulting cost will reflect the quantum hardware's noise. In this work, we evaluate the cost on a classical computer, as we are mainly interested in the discovery of theoretical algorithms without device-specific noise considerations.

Another aspect of the algorithm scaling is training data. In general, its size will scale exponentially with the number of data qubits. However, we would like to stress that this fact does not jeopardize our approach since we numerically obtain solutions (algorithm instances) only for a small number of data qubits. Those optimization problems require training data that is still manageable. Algorithm instances are then used to manually recognize the pattern and generalize to arbitrary system size.

Finally, the search space defined by $\vec{k}$ is exponentially large in the number of gates. This makes it impossible to systematically check all possibilities in the search for an optimal algorithm. On the other hand, the heuristic approach described above seems to be capable of finding the solution efficiently.

For a fixed problem size, we minimize the cost. If the cost goes to zero (which we define as a cost less than 10⁻⁶), we say we have an algorithm instance. In particular, this corresponds to fixing the size of the data and hence fixing n_D, the number of data qubits. To study the generalization of the algorithm, we grow the size of the problem by increasing n_D. In some cases, one may also need to increase the number of ancilla qubits, n_A, and/or the number of measurements in order to minimize the cost.

This gives us a set of algorithm instances for various problem sizes. An important challenge is to abstract a general algorithm from these instances. This challenge is particularly difficult because one can typically only find algorithm instances for small problem sizes. This is due to the fact that the search space for vectors $\vec{k}$ grows rapidly with problem size, namely as ${n}_{T}^{2d}$, where n_T = n_D + n_A is the total number of qubits and d is the circuit gate count.

In this work, we were able to manually recognize the pattern by which the algorithm generalizes to arbitrary problem size by inspecting the various algorithm instances. In future work, we will explore automated methods for recognizing the general algorithm.

Our main results are short-depth algorithms for quantifying overlap on idealized quantum computing hardware. For the latter, we consider full connectivity, and we allow for arbitrary one-qubit gates as well as $\mathrm{CNOT}$ gates between all of the qubits.

We consider two sets of resources. The first set of resources are identical to those allowed for the Swap Test, i.e. access to one ancilla qubit and two data qubits, as well as one measurement on the ancilla qubit. The cost versus number of gates for these resources is shown in figure 3(A), and we obtained essentially zero cost for d = 8. To understand how the algorithm generalizes, we increase the number of qubits in ρ and σ to n = 2, giving a minimum gate count of d = 14, as shown in figure 3(B). As discussed below this generalizes to an algorithm (shown in figure 5) that we refer to as our ABA.

Zoom In Zoom Out Reset image size

The second set of resources we consider allows for measurements on all of the qubits. For these additional resources, figure 3(C) shows that zero cost is obtained for d = 2. To recognize how this algorithm generalizes, we increase the number of qubits to n = 2, giving a minimum gate count of d = 4, as shown in figure 3(D). The surprising result is that the ancilla qubit is not used at all in this algorithm, even though we train the algorithm in the presence of an ancilla. This allows us to display the resulting general algorithm, our BBA, in figure 6 without the ancilla qubit.

Zoom In Zoom Out Reset image size

In both cases discussed above, we managed to discover the general (valid for arbitrary problem size) form of the algorithm from its two smallest instances. We expect that in other applications, the general form of the algorithm may be harder to find and more sophisticated tools will have to be developed.

Figure 5(A) shows the ABA for one-qubit states ρ and σ. The unitary U in this circuit is $U={T}^{\dagger }H$. This circuit employs 4 CNOT gates and 4 one-qubit gates for a total of 8 gates. It uses a simple post-processing vector $\vec{c}=(1,-1)$ that amounts to measuring the Pauli Z operator on the ancilla qubit, which is the same observable measured in the Swap Test. Not only does this circuit have a lower gate count than typical implementations of the Swap Test (see e.g. the circuit in figure 1(B)), but actually it implements a completely different unitary.

Let S_ABA denote the Schmidt rank (across the cut between ancilla and the data qubits) of the unitary G_ABA associated with the ABA gate sequence. It can be verified that S_ABA = 3. This means that G_ABA is not locally equivalent to a controlled-SWAP, whose analogously defined Schmidt rank is 2. Thus, the ABA is fundamentally different from the Swap Test: it cannot be obtained from the Swap Test by local operations.

The general form of the ABA is given in figure 5(B). There is a repeating unit, shown in the inset of the figure, that is applied on each pair of qubits composing ρ and σ as well as on the ancilla qubit. This unit has 4 CNOT gates, so the overall algorithm employs $4n$ CNOT gates and $6n+2$ total gates. Hence, the gate count grows linearly with the number of data qubits.

Figure 6(A) shows the BBA for one-qubit states ρ and σ. This circuit employs one CNOT gate followed by one Hadamard gate, with both qubits being measured. It is straightforward to show that this corresponds to a Bell basis measurement. The post-processing is a bit more complicated, with $\vec{c}=(1,1,1,-1)$, which corresponds to summing the probabilities for the 00, 01, and 10 outcomes and subtracting probability of the 11 outcome. The above post-processing is equivalent to measuring the expectation value of a controlled-Z operator.

The generalization of this algorithm is given in figure 6(B). The repeating unit is simply a CNOT and Hadamard, applied on each pair of qubits composing ρ and σ. Furthermore, every qubit is measured at the output. The total number of gates is simply $2n$, and hence grows linearly with the number of qubits. However, more importantly, the CNOT and Hadamard on each qubit pair can be performed in parallel. This crucial fact means that this algorithm has a constant depth, independent of problem size. Namely, the depth is two quantum gates.

On the other hand, the classical post-processing is somewhat complicated, and its complexity scales linearly with the problem size. Namely, the post-processing vector can be written as $\vec{c}={(1,1,1,-1)}^{\otimes n}$, provided that one arranges the qubits in the order P₁Q₁P₂Q₂ .... P_nQ_n, where P₁P₂ ... P_n and Q₁Q₂ ... Q_n are the subsystems composing ρ and σ respectively. The linear scaling of post-processing follows from the fact that one does not explicitly compute $\vec{c}\cdot \vec{p}$ in equation (3). Rather one bins individual measurement outcomes into one of two bins (either the 1 or −1 bin). Here, the bin is determined by first assigning each of the n qubit pairs a value of 1 or −1, based on the associated eigenvalue of the controlled-Z operator, and then multiplying these n values. The overlap $\mathrm{Tr}(\rho \sigma )$ is then given by the weighted average over all outcomes, where the weights correspond to the bin label (either 1 or −1).

Nevertheless, for NTQCs, due to decoherence and gate infidelity, it is better for the classical post-processing to grow linearly in n than for the quantum circuit depth to grow linearly in n. Hence, the BBA seems to be the superior algorithm in that case.

In 2013, Garcia-Escartin and Chamorro-Posada discovered the BBA for computing state overlap [24]. We were unaware of this important result until after our machine-learning approach found our BBA. More generally, it appears that the quantum computing community seems to be unaware of this article, perhaps because the article was presented in the language of quantum optics rather than that of quantum computing. Indeed, the ancilla-based version of the Swap Test, shown in figure 1, continues to be the algorithm employed in the quantum computing literature (e.g. see [27, 33]).

Although our two algorithms look very different, one can actually show a simple equivalence between our ABA and our BBA. One can see this by converting the classical post-processing in the BBA into a quantum gate. In particular, this gate would be a Toffoli gate, controlled by the two data qubits with the target being an ancilla qubit prepared in the $| 0\rangle$ state. Appendix B shows proof of this statement. After inserting the Toffoli gate (see figure 7(B)), one would do a measurement of the Pauli Z observable on the ancilla to decode the state overlap. By replacing the Toffoli gate with its decomposition from [35] and simplifying the resulting circuit, one can then obtain our ABA (see figure 7(C)). In this sense, our ABA is essentially our BBA but with the classical post-processing transformed into Toffoli gates and a measurement on the ancilla. This equivalence is shown in figure 7 for one-qubit states. The generalization to multi-qubit states is straightforward.

Zoom In Zoom Out Reset image size