Quantum annealing speedup of embedded problems via suppression of Griffiths singularities

For this study we set out to explore strategies for setting parameters for embeddings in coherent quantum annealers, which replace logical qubits with blocks of $M$ ferromagnetically coupled physical qubits. This is done primarily to increase the connectivity of qubit interaction graph [to $M(c-2)+2$, where $c$ is the degree of connectivity graph of physical hardware]. The extreme example achieves all-to-all connectivity at the cost of quadratic reduction of the number of logical qubits [to $O\bigl{(}\sqrt{N}\bigr{)}$] as seen in Refs. Choi, 2008, 2011. Another important application is extending the range of Ising couplings, which is more relevant for the present study of a linear chain ($c=2$). Lastly, using block spins has been suggested as a form of error-correction Pudenz et al. (2014); *pudenz2015quantum; Vinci et al. (2015); Matsuura et al. (2017).

As we introduce $M$ ancillary spins for each logical variable (see Fig. 1), the classical Hamiltonian reads

For simplicity, the ferromagnetic couplings $K_{i}$ within a subchain are taken to be uniform but we allow variations from one logical qubit to another. To ensure that these ferromagnetic links are never broken in any local minimum we insist that

Since the model retains its 1D character, existing analytical techniques are still applicable and numerical studies can explore a regime where annealing becomes intractable for large $N$.

Specifically, we observe that for sufficiently small values of the transverse field $\tilde{\Gamma}<\min_{i}K_{i}$ we may perform the real-space renormalization procedure due to D. S. Fisher Fisher (1992); *fisher1995critical. Each subchain now behaves as a spin subjected to a transverse field with the renormalized value

(In a reversal of notation used in Ref. Fisher, 1992; *fisher1995critical, we write $\tilde{\Gamma}=\frac{s}{1-s}$ to represent the bare value of the transverse field and reserve non-accented letters $\Gamma_{i}$ to denote the renormalized fields.) The effective model is still an Ising chain

albeit with the local transverse field that can potentially have spatial variations.

The embedding procedure thus unlocks an important resource since many existing implementations either lack such local control entirely or have limited dynamic range. It allows us to restore the polynomial scaling of annealing complexity, balancing out disorder fluctuations to synchronize the phase transition across entire chain. We make the following ansatz for embedding parameters

for some rescaling constant $C.$ At the edges we can use, e.g., $K_{0}=CJ_{1}^{-1/(M+1)}$ and $K_{N-1}=CJ_{N-1}^{-1/(M+1)}$.

With this choice, renormalized local fields in the bulk become $\Gamma_{i}=\frac{\tilde{\Gamma}^{M}}{C^{M-1}}\sqrt{J_{i}J_{i+1}}$ so that fluctuations of the value of the critical field $\Delta\tilde{\Gamma}_{c}$ for large subchains now scale as $1/N$.

In Fig. 2 we compare the median minimum gap for different parameter setting prescriptions, showing that the balanced choice in Eq. (8) recovers a polynomial gap. It is imperative to discuss how the parameters of our numerical study were chosen. The black solid line plots the estimate of the median value (taken over disorder distribution) of the minimum gap for strong disorder: $J_{i}$ are taken to be uniformly distributed in the interval [0;1]. The data clearly cannot be fit by any power law (which would have been a straight line on a log-log plot).

It is evident from Eq. (8) that allowing arbitrarily small values of $J_{i}$ would in turn lead to arbitrarily large values of $K_{i}$. Moreover, small values of $J_{\text{min}}$ are not expected to appear naturally for Ising representations of practical problems Lucas (2014) and would make the problem particularly prone to misspecification errors in hardware annealers.

To mitigate this problem, we choose $J_{i}\in[1/M;1]$, where $M$ is the number of spins in a block. Fig. 2 uses colored dotted line to plot gaps of the logical problem for such disorder distributions. Dashed lines show the gaps for an embedded problem with a canonical choice of embedding parameters ($K_{i}=1$). Balanced choice of embedding parameters from Eq. (8) with $C=1$ would have yielded values $K_{i}$ in a range $[1;M^{1/(M-1)}]$. Current practice for existing hardware implementations of quantum annealing is to use largest possible values of couplings $J_{i}$ in order to minimize misspecification error, hence it is natural to assume that $J=1$ already represents the maximum programmable value. In order to accommodate the balanced parameter setting, all ferromagnetic couplings must be scaled down by a factor of $M^{1/(M-1)}$, i.e. $J_{i}\in[M^{-M/(M-1)};M^{-1/(M-1)}]$ and $C=M^{-M/(M-1)^{2}}$. This rescaling step had been made for the data plotted in Fig. 2 using solid colored lines. The motivation for this peculiar choice of disorder distribution is that it maximizes the range $J_{\text{max}}-J_{\text{min}}$ of inter-chain couplings after the rescaling.

Rescaling of couplings parameters downward handicaps quantum annealing performance. Using rescaled $J^{\prime}=\lambda J$ we can show that $\Delta E_{c}^{\prime}\cdot\Delta s_{c}^{\prime}=(\lambda^{2}s_{c}/s^{\prime}_{c})\Delta E_{c}\cdot\Delta s_{c}$. The scaling factor can be rewritten as $(1-s_{c})\lambda^{3}+s_{c}\lambda^{2}$, which is less than unity if $\lambda<1$. Numerical data for the gap is consistent with this qualitative behavior. More rigorous analysis of the low-energy spectrum along the lines described in Appendix A also predicts that the gap for canonical embedding is smaller than that of the logical problem with the same disorder by $O(M)$ factor, also in agreement with the numerical data.

We encountered a problem specific to instances using balanced embedding, likely due to a quirk in its spectrum: numerical diagonalization could only performed for a smaller range of $N$ as the limits of machine precision were reached. It may be possible to overcome this obstacle using multiprecision arithmetic to find the roots of characteristic polynomial for the tridiagonal matrix. Fortunately, the median gap fits polynomial scaling perfectly (straight line on a log-log plot), which gives us sufficient confidence to extrapolate the data for large $N$.

The one-dimensional Ising chain is the simplest example of an integrable quantum model. A rotation ($X_{i}\rightarrow-Z_{i}$, $Z_{i}\rightarrow X_{i}$) followed by a Jordan-Wigner transformation expresses the Hamiltonian as a quadratic form in fermionic operators. The solution of time-dependent Bogolyubov-de Gennes (BdG) equations is used to write these operators in Heisenberg representation (which is a linear combination of Schrödinger operators since $H$ is quadratic Dziarmaga (2006); Caneva et al. (2007)). For instance, using Majorana representation (see e.g. Ref. Fradkin, 2013) ($\chi_{2i}=X_{i}\prod_{k=0}^{i-1}Z_{k}$ and $\chi_{2i+1}=Y_{i}\prod_{k=0}^{i-1}Z_{k}$) we can write

Here $\mathbf{S}$ is special orthogonal $2n\times 2n$ matrix, whereas $\mathbf{M}$ is skew-symmetric and also happens to be tridiagonal with the elements on the upper diagonal $M_{2i,2i+1}=1-s$ and $M_{2i-1,2i}=sJ_{i-1}$. This special structure further aids numerics; on top of logarithmic reduction in complexity afforded by the free-fermion mapping.

Eigenvalues of $\mathbf{M}$ (e.g. for ordering $\lambda_{2k}=+\mathrm{i}\epsilon_{k}$, $\lambda_{2k+1}=-\mathrm{i}\epsilon_{k}$ with $0<\epsilon_{0}<\cdots<\epsilon_{N-1}$) completely determine the spectrum of the instantaneous Hamiltonian $H(t)$ as a sum of single-particle excitations: $E_{\{n_{k}\}}=2\sum_{k}n_{k}\epsilon_{k}$ where $n_{k}\in\{0,1\}$. Specifically, the gap between the ground state and the second excited state (which is second lowest in energy among those in the symmetric subspace) is $\Delta E=2(\epsilon_{1}+\epsilon_{2})$. This value enters the adiabatic condition used to estimate the runtime (see Appendix A).

To better quantify the performance we concentrate on a TTS metric of Eq. (2). The probability to find the system in its ground state at the end of the annealing cycle at $t=T$ is

where $b_{k}$ and $b_{k}^{\dagger}$ are fermionic quasiparticle operators in Dirac representation that diagonalize the final classical Hamiltonian $H_{C}$. Expectation value can be computed at time $t=0$ provided that Heisenberg representation for $b_{k}(t)$ and $b_{k}^{\dagger}(t)$ is used. The initial state satisfies $a_{i}|0\rangle=0$ where $a_{i}$ and $a_{i}^{\dagger}$ are the Dirac fermion operators that diagonalize the transverse field Hamiltonian $H_{D}$.

Since the Hamiltonian is quadratic, the expectation value (10) is computed by performing the sum over all possible pairings and can be expressed as

where the only non-zero elements of $\mathbf{A}$, $\mathbf{B}$ are $A_{2k,2k+1}=B_{2k-1,2k}=B_{0,2N-1}=1$ and $A_{2k+1,2k}=B_{2k,2k-1}=B_{2N-1,0}=-1$ (see Appendix B).

Investigating the complexity using the TTS metric is somewhat challenging numerically. It entails integrating systems of differential equations of more than $10^{4}$ variables for multiple choices of annealing times (in excess of $10^{7}$ in dimensionless units). Numerical error accumulates proportionally to the evolution time, and is further amplified when we evaluate the determinant (11). After multiple tests, we implemented an integrator based on Cayley transform and Magnus expansion Blanes et al. (2000); *blanes2002high up to 8th order. However, we use Padé approximation of the exponential to obtain a straightforward mapping to Cayley transform and exploit the tridiagonal structure of matrix $\mathbf{M}$ to attain best performance. Compared with traditional Runge-Kutta approaches used for these types of problems Mbeng et al. (2019), our method maintains orthogonality of $\mathbf{S}$ at all update steps which significantly improves stability for long integration times (see Appendix C). For a given timestep our implementation outperforms comparable methods such as Dormand-Prince algorithms that are the method of choice in celestial mechanics Dormand and Prince (1980), but also maintains excellent precision for large step sizes $\Delta t$.

Note that since optimal annealing time in Eq. (2) cannot be known in advance, it seems unfair to perform optimization of for each individual instance. We use ensemble-wide $T_{\text{opt}}(N)$, which is determined through minimization, as a function of annealing time, of the median TTS, from a sample of random instances of the same size. This choice is similar to the practice of benchmarking of application problems.

We plot the median time-to-solution for random instances of increasing size in Fig. 3. As before, colored dotted lines refer to logical problem with $J_{i}\in[1/M;1]$, dashed lines correspond to canonical embedding parameters and solid lines present results for balanced embedding parameters. Finally, black solid line refers to a logical problem with strong disorder $J_{i}\in[0;1]$, where non-polynomial behavior is most pronounced. We observe that for the largest sizes balanced embedding clearly outperforms the canonical one ($K_{i}=1$).

We have investigated quantum annealing of embeddings of a simple problem, where logical qubits were replaced by ferromagnetic chains of Ising spins. A novel ansatz for a balanced choice of coupling parameters based on renormalization group intuition results in an exponential improvement in annealing times. This is corroborated numerically, using time-to-solution metric of complexity. For large sizes the balanced choice significantly outperforms both canonical parameter setting and the quantum annealing of the logical problem using the same distribution of disorder (but with no embedding).

The protocol proposed in this letter succeeds because embedding parameters are set in such a fashion that spatially separate regions achieve criticality simultaneously. A situation were multiple domains oriented in random directions are created is thus avoided. Remarkably, this can be accomplished without local control of the transverse field. Complementary approaches described in literature rely on “growing” the domain and require significant modifications to quantum annealing protocol using time-dependent local control of the transverse field Rams et al. (2016).

An important venue of research is the generalization of this idea to higher-dimensional problems. Dearth of exactly solvable models suggests using approximate methods to find embedding parameters that achieve synchronization of local phase transitions. We expect to be able to suppress Griffiths singularities relieving the phase transition bottleneck of annealing complexity. However, Ising models in higher dimensions exhibit frustration, which dominates complexity for very large instances Knysh (2016) (crafted problems in 1D can also exhibit frustration bottleneck Roberts et al. (2020)). We are cautiously optimistic that this general technique can be useful for intermediate-range problems that could be implemented in highly coherent quantum annealers Novikov et al. (2018). Another research direction is the extension of this work to open systems where the annealing dynamics is incoherent, driven by coupling with an external reservoir Smelyanskiy et al. (2017); Mishra et al. (2018). This would open up the opportunity of testing the balanced embedding approach on D-Wave machines as well Bando et al. .