---
abstract: 'This document contains additional plots proving the validity of the results in the main paper and plots showing intermediate measurements used for the adaptive schedule. Furthermore different parametrization for nonlinear annealing schedules are tested.'
author:
- Daniel Herr
- Ethan Brown
- Bettina Heim
- Mario Könz
- Guglielmo Mazzola
- Matthias Troyer
bibliography:
- 'biblio.bib'
title: 'Supplementary Material for: Optimizing schedules for Quantum Annealing'
---

Methods used in the main paper {#methods-used-in-the-main-paper .unnumbered}
------------------------------

In the following a more thorough description of algorithms used to simulate both Classical annealing (CA) and Quantum annealing (QA) is given.

### Classical Annealing {#classical-annealing .unnumbered}

In CA [@kirkpatrick1983] the system is first initialized to a random configuration at a high temperature and then gradually cooled. This allows the system to escape from local minima and to relax towards lower energy configurations [@sa_and_sa_conv_cond]. During each step the Metropolis-Monte-Carlo algorithm [@Metropolis] allows re-equilibration. Each step of this algorithm proposes a change to the current configuration. The change is always accepted if favorable, whereas if the cost increases it will only be accepted by a probability dependent on the temperature of the system.

### Quantum Annealing {#quantum-annealing .unnumbered}

During quantum annealing (QA) [@Ray1989; @Finnila1994; @Kadowaki1998; @idea_of_qa; @qareview] the adiabatic theorem enables the system to evolve from a trivial ground state to the solution of the complex problem $\mathcal{H}_P$, if the change in the Hamiltonian is sufficiently slow. For this one can state a new system Hamiltonian $$\mathcal{H}=s \mathcal{H}_P +  \left(1-s\right) \mathcal{H}_D
    \label{eq:generalhamil}$$ that introduces the control parameter $s \in \left[0 ,1\right]$ which enables a transition from an initial system $\mathcal{H}_D$ whose ground state can be easily obtained to the problem Hamiltonian. In contrast to CA, the system is held at a constant inverse temperature $\beta=1/T$ while the control parameter is slowly increased to 1. Thus, the QA schedule takes the system from a preparable state to a state where quantum fluctuation are suppressed and the system remains frozen in its configuration. The problem Hamiltonian encodes the classical problem that is about to be solved. In QA classical values of $s_i$ are replaced by the Pauli spin-z operators $\sigma_i^z$. And the driver Hamiltonian $\mathcal{H}_D$ introduces randomness to the quantum and is given by $$\mathcal{H}_D = \Gamma_0 \sum_i \sigma^x_i.$$ Here $\sigma_i^x$ is the Pauli spin-x operator for the $i$th spin and enable transitions between $\uparrow$ and $\downarrow$ for a nonzero transverse field $\Gamma_0$.

To simulate QA we used a discrete time SQA algorithm. We use the path integral Monte Carlo (PIMC) method to map the 3 dimensional system into $3+1$ dimensions by the introduction of the imaginary time dimension [@suzuki_orig]. This extra dimension is discretized with $M$ Trotter slices, which are copies of the classical Ising system coupled to each other by a value dependent on $\Gamma$ and $\beta$. As the transverse field decreases the coupling between the Trotter slices gets stronger. In the classical limit $\Gamma = 0$ the quantum fluctuation are suppressed and the system remains frozen in its configuration. The simulation using a finite number of Trotter slices $M$ is called discrete time SQA (DT-SQA) and comes with a discretization error of $O(\beta^3/M^2)$. The physical limit corresponds to the continuous time limit $M\rightarrow \infty$. To compare the computational effort between QA and CA the number of Monte Carlo sweeps of SQA and CA are compared. In SQA a single update consists of a Swendsen-Wang [@swendsenwang] cluster move in imaginary time direction only opposed to the single spin flip CA performs. This has proven to be a reliable comparison for the individual runtime [@Heimb15].

Since recent work by some of the authors indicates that using a PIMC technique with open boundary conditions in imaginary time is the best possible SQA algorithm and may better reproduce the scaling with system’s size of a coherent quantum annealer [@mazzola2017quantum; @Guglielmo], the simulations were run accordingly. The results serve as a best case estimate for the performance of a physical quantum annealer, e.g. the D-Wave 2, since the coherence length in such devices may or may not extend across the entire system. Furthermore all simulations were run in the physical limit where the number of Trotter slices is required to be high enough to ensure convergence. This leads to $M=1024$.

Observables in PIMC {#observables-in-pimc .unnumbered}
-------------------

To evaluate the expectation value of the QA counterpart of the specific heat $C_q$, we need to compute the value of $\left<\sigma_x\right>$. A derivation on how to calculate this observable using the PIMC method can be found in [@Krzakala08] and will be reproduced in the following.\
Generally an observable can be evaluated as: $$\left<O\right> = \frac{1}{Z} \text{tr}\left( O e^{-\beta \mathcal{H}}\right)$$ Thus the expectation value of $\left<\sigma_x\right>$ will evaluate to the following: $$\left<\sigma_i\right>=\frac{\sum_{\sigma_i} \bra{\sigma_i} \sigma_i^x e^{-\beta\mathcal{H}}\ket{\sigma_i}}{\sum_{\sigma_i} \bra{\sigma_i} e^{-\beta \mathcal{H}} \ket{\sigma_i}}$$ Introducing the Trotter decomposition $e^{\beta\mathcal{H}} = {\left( e^{\frac{\beta}{M} \mathcal{H}}\right)}^M$ with $\tau = \frac{\beta}{M}$ and adding identity operators will give $$\left<\sigma_i\right>\Big|_k = \frac{\sum_{\sigma_i} \bra{\sigma_i} \sigma_i^x e^{-\beta\mathcal{H}}\ket{\sigma_i}}{\sum_{\sigma_i} \bra{\sigma_i} e^{-\beta \mathcal{H}} \ket{\sigma_i}}.$$ But there is freedom to choose for which Trotter slice to evaluate $\sigma_i^x$. So one can also take the average of all these choices: $$\left<\sigma_i\right>= \frac{1}{M} \sum_{k=1}^{M} \left<\sigma_i\right>\Big|_k$$ Now with the introduction of a small error of order $O(\tau^2)$ one can split $e^{-\tau (\mathcal{H}_p + \Gamma\sum_i\sigma_i^x)}= e^{-\tau \mathcal{H}_p} e^{-\tau \Gamma \sum_i \sigma_i^x}$ and then find for the expectation value $$\left<\sigma_i^x\right> = \frac{1}{M}\sum_{k=1}^{M} \left( \frac{\bra{\sigma_i^{k+1}} \sigma_i^x e^{-\tau \Gamma \sum_i \sigma_i^x }\ket{\sigma_i^k}}{\bra{\sigma_i^{k+1}} e^{-\tau \Gamma \sum_i \sigma_i}\ket{\sigma_i^k}} \right)$$ Using $\bra{\uparrow} e^{a\sigma^x}\ket{\uparrow} = \cosh(a)$ and $\bra{\uparrow} e^{a\sigma^x}\ket{\downarrow} = \sinh(a)$ one can obtain the expectation value of the Pauli spin x operator: $$\left<\sigma_i^x\right>= \frac{1}{M} \sum_{k=1}^{M} {\tanh\left(-\tau \Gamma \right)}^{s_i^k s_i^{k+1}}$$

![\[fig:sigmax\] The expectation value of $\sigma_x$ operator depending on $\Gamma$](sigmax.pdf){width="\linewidth"}

In Figure \[fig:sigmax\] the expectation value of $\sigma_x$ is plotted for different values of the transverse field. For the more specific schedule relying on the Hamiltonian $$\mathcal{H} = \mathcal{H}_p + \Gamma(t) \sum_{i=0}^N \sigma_x^i$$ one can get an expectation value for $C_q$ by the following formula: $$C_q = \frac{1}{\beta} \frac{d^2}{d\Gamma^2} \log(Z)= \beta\Gamma\left(1 - \left<\sigma_x\right>^2\right)$$ Finally one should note that any constant prefactor in the measurement of $C_q$ will cancel further on when $\lambda$ is determined by fixing the number of MCS.

Schedules {#schedules .unnumbered}
---------

For optimized parameters an exponential schedule and a mixture of both classical and quantum annealing, where the transverse field was linearly decreased while the inverse temperature was linearly increased, were tested. The performance of these schedules is plotted in Figure \[fig:nonlin\_perf\]. Even the hybrid between CA and QA did not show a considerable performance increase, as can be seen in Figure \[fig:BG\_sched\]. The additional parameters that need to be optimized only seem to result in marginal changes. The residual energy seems behave similarly to SQA runs with fixed high temperature at few MCS, while still obtaining a slightly lower state which low temperature runs possess for large numbers of MCS.

![\[fig:nonlin\_perf\]Performance of the different schedules for optimized start and end parameters](Ising_final.pdf){width="\linewidth"}

![\[fig:convergence\]Change in number of trotter slice to check if DT-SQA already converged](convergence.pdf){width="\linewidth"}

![\[fig:BG\_sched\]Median residual energy for a schedule with linear change in both $\beta$ and $\Gamma$](BG_results.pdf){width="\linewidth"}

Notes on the adaptive schedule {#notes-on-the-adaptive-schedule .unnumbered}
------------------------------

The adaptive schedule presented in this paper is derived using the system Hamiltonian given in Eq. \[eq:generalhamil\]. Yet, a more common annealing procedure [@Heimb15; @santoro2] is when the problem Hamiltonian is kept at constant strength while the transverse field is decreased from $\Gamma_0$ to $0$. The simulations performed in this paper use this more specific schedule.

Validity of Results {#validity-of-results .unnumbered}
-------------------

In Figure \[fig:convergence\] one can see that at 1024 Trotter slices a change in discretization has almost no influence on the behavior of QA. Thus one can conclude that the algorithm is sufficiently converged.

Ferromagnet {#ferromagnet .unnumbered}
-----------

Additionally the same analysis was conducted for a 3D Ferromagnet whose degenerate ground state was lifted by a local field in $\sigma_z$ direction. The parameter for the annealing speed is plotted in Figure \[fig:sched\_ferro\] and in Figure \[fig:Ferro\_perf\] for different unoptimized annealing parameters. Again one can see that the adaptive schedule is better than the linear schedule for unoptimized parameters. But we noticed that optimized values, can further improve the performance such that a well optimized starting value beats the unoptimized adaptive schedule in the case of a Ferromagnet. Still this schedule becomes more advantageous the farther from optimized start values the annealing run is.

![\[fig:sched\_ferro\]Measured indicator for the annealing speed of a 3D Ferromagnet with $\beta = 32$](Ferro_CQ.pdf){width="\linewidth"}

![\[fig:Ferro\_perf\]Performance comparison between the adaptive and linear for a 3D Ferromagnet](Ferro_unoptimized.pdf){width="\linewidth"}

Comparison to 2D {#comparison-to-2d .unnumbered}
----------------

Simulations not in the converged limit show the exact same behavior as for 2D Ising spin glasses, cf. Figure \[fig:bettina\_B\] and the paper on 2D Ising spin glasses [@Heimb15]. The same can be observed for the converged case, that is presented in the main paper.

![\[fig:bettina\_B\]temperature behavior of the DT-SQA algorithm without the requirement of convergence to the continuous limit](bettina_B.pdf){width="\linewidth"}
