The Effect of Penalty Factors of Constrained Hamiltonians on the Eigenspectrum in Quantum Annealing

Quantum annealing is a meta-heuristic most commonly known for solving optimization and decision problems [19, 21, 28]. Although this meta-heuristic can also be simulated classically, it has been implemented in quantum hardware by companies such as D-Wave Systems. Those quantum annealers are designed to minimize a spin glass system, described by an Ising Hamiltonian in the following form:where \(\sigma _z^{(i)}\) is the Pauli z-matrix operating on qubit i, \(h_{i}\) is the independent energy or bias of qubit \(i,\) and \(J_{ij}\) are the interaction energies or couplings of qubits i and j.

Within the fundamental process of quantum annealing, an initial Hamiltonian \(\mathcal {H}_{I}\) with an easy-to-prepare minimal energy configuration (or ground state) is physically interpolated to a problem Hamiltonian \(\mathcal {H}_{P}\) whose minimal energy configuration is sought (see Equation (2)). The minimal energy configuration of the problem Hamiltonian corresponds to the best solution of the defined problem. The physical principle on which the D-Wave computation process is based on can be described by a time-dependent Hamiltonian as follows

\(\mathcal {A}(t)\) and \(\mathcal {B}(t)\) are the anneal functions of D-Wave machines, with \(\mathcal {A}(t)\) stating the tunneling energy and \(\mathcal {B}(t)\) being the energy of the problem Hamiltonian at time t in units of joules. The anneal functions must satisfy \({\mathcal {B}(t = 0) = 0}\) and \({\mathcal {A}(t = \tau) = 0}\), with \(\tau\) being the total evolution time. As the state evolution changes from \(t=0\) to \(t=\tau\), the annealing process, described by \(\mathcal {H}(t),\) leads to the final form of the Hamiltonian corresponding to the objective Ising problem that needs to be minimized. Therefore, the ground state of the initial Hamiltonian \(\mathcal {H}(0) = \mathcal {H}_I\) evolves to the ground state of the problem Hamiltonian \({\mathcal {H}(\tau) = \mathcal {H}_P}\). The measurements performed at time \(\tau\) deliver low energy states of the Ising Hamiltonian as stated in Equation (1).

According to the adiabatic theorem [2], if this process is executed sufficiently slow and smooth (i.e., \(\tau\) is large) and the coherence is preserved long enough, the probability to acquire the ground state of the problem Hamiltonian is close to 1 [1]. However, since no real-world computation can run in perfect isolation, the annealing process can suffer from non-adiabatic effects (i.e., thermal fluctuations), which can lead the system to jump from the ground state to an excited state. The minimum distance between the ground state and the first excited state—the one with the lowest energy apart from the ground state—throughout any point in the anneal process is called the minimum spectral gap \(g_\text{min}\) of \(\mathcal {H}(t)\) and is defined aswhere \(E_j(t)\) is the energy of any excited state and \(E_0(t)\) is the energy of the ground state at time t [17]. By computing all those energy states and their corresponding eigenenergies, one can analyze the eigenspectrum of a (relatively small) Hamiltonian and assess its minimum spectral gap (Figure 1 presents an example of an eigenspectrum). However, every problem that one can specify has a different Hamiltonian and therefore a different corresponding eigenspectrum. According to D-Wave Systems, the most difficult problems in terms of quantum annealing are generally those with the smallest spectral gaps [11]. For completeness, it should be noted that there is an alternative formulation to the Ising spin glass system that is used frequently. The so-called Quadratic Unconstrained Binary Optimization (QUBO) formulation is mathematically equivalent to the Ising model and replaces each Pauli z-operator \(\sigma _z^{(i)}\) with a Boolean variable \(x_i\); the conversion is as simple as setting \(\sigma _z^{(i)} \rightarrow 2x_i - 1\) [3, 29]. The D-Wave Systems annealer is also able to minimize the functional form of the QUBO formulation \(x^TQx\), with \(x \in \lbrace 0,1\rbrace ^n\) being a vector of size n of binary variables and \(Q \in \mathbb {R}^{n \times n}\) being a symmetric \(n \times n\) real-valued matrix describing the interactions between the variables. Given matrix Q, the annealing process tries to find binary variable assignments x that minimize the objective function.

In this section, the six investigated problem Hamiltonians with constraints are stated. The constraints are weighted with the penalty factors to ensure valid solutions for the corresponding problem.

In this section, we describe the procedure of generating the data and computing the necessary information—that is, penalty factor ratio, the minimum spectral gap, and its location in time w.r.t. the anneal path—as input for the ML analysis methods. As mentioned in Section 2.2, each constrained Hamiltonian has a certain valid half-open interval for its penalty factor A. Although we fixed the penalty factor \(B=1\), the sampling interval of penalty factor A was computed according to the problem-specific constraint (see Equations (5), (7)–(15)).

Using the MECP as an example, the sampling interval for A was set to \(\left[B\cdot n+0.1, B\cdot n+5.0\right]\), with n being the number of sets of the problem instance. The MECP penalty factor ratio was calculated via \(\frac{B\cdot n}{A}\). Note that we restricted the in general half-open interval of A to \(\left[B\cdot n+0.1, B\cdot n+5.0\right]\). Although one could theoretically use a larger sampling range, one has to consider that D-Wave Systems’ auto-scaling feature scales every Ising model weight and bias to a given hardware solver dependent range of \(h_{i}\) and \(J_{ij}\) [12]. This means that the largest value in the Ising model (probably affected by the constraint penalty factor A) is scaled to the upper bound of the range and everything else proportionally smaller. In conclusion, to get good and meaningful solutions, one should set the Ising model coefficients of the constraints large enough to ensure valid solutions but small enough to maintain the importance of the individual Ising terms.

For each problem instance, we draw 50 evenly spaced penalty factors for A from the interval to see the rate of change of the minimum spectral gap and its location in the anneal path. The minimum spectral gap of each constraint Hamiltonian with its unique sampled penalty factors was then calculated according to D-Wave Systems’ anneal functions \(\mathcal {A}(t)\) and \(\mathcal {B}(t)\) of their Advantage 4.1 system [10]. Note that we also incorporated D-Wave Systems’ auto-scaling feature during the generation of the data, before computing the minimum spectral gap, so that the experiments match the behavior of the D-Wave hardware [12].

In Figure 2, the corresponding preliminary raw plots of 25 random instances for each problem class are visualized, to understand the following experiments in Section 5. However, for analyzing the trends of the minimum spectral gap, its location in the anneal path, and the penalty factor ratio with ML methods, we generated 1,000 unique problem instances per problem class. Note that we restricted the size of the instances to eight variables, respectively logical qubits, since the numerical computation of the eigenspectrum of each instance is quite time consuming and increases with the number of variables.

Note that to generate such large datasets of Ising problem instances and their minimum spectral gaps for training ML methods, we computed the eigenspectra of the comparatively smaller logical Ising problem instances, which differ to some extent from the hardware embedded ones. We address this topic in Appendix A.

Since it is hard to visually find trends in the 1,000 generated problem instances of each problem class, DBSCAN was used to cluster the generated trajectories as described in Section 4. In Figure 3, the DBSCAN clustering results of the minimum spectral gap and its location over the set of penalty factor ratios are plotted. The axis labeling is exactly as in Figure 2. In each subplot, the mean trajectory and the \(95\%\) confidence interval, as error bars, of each found cluster are plotted. In addition, the cluster label is annotated in text form. Table 1 shows the associated best found DBSCAN hyperparameters and the achieved PCC for each problem class, which was used as metric for the clustering.

In general, the clustering of MVCP, MCP, and SPP are similar in that they all exhibit multiple low-lying clusters w.r.t. the size of the minimum spectral gap (trajectories at the very bottom of each plot) and a location of around 60 to 90 in the discretized anneal path ranging from 0 to 100 (see the y-axis). Additionally, they all have a cluster growing in terms of the size of the minimum spectral gap with increasing penalty factor ratio and a comparatively early gap location of 10 to 20 (cf. Figure 3(b), (c), and (e)). This means that in theory the problem instances of these clusters can be optimized with a certain penalty factor ratio w.r.t. the minimal spectral gap size. Moreover, we assume that the early spectral gap of those instances might be favorable, since the coherence time of NISQ computers, like D-Wave Systems, is still limited. With an early minimum spectral gap, the probability to jump from the ground to an exited state due to natural quantum decoherence might be decreased. Note that for MCP, there is one cluster (label \(-1\)) of problem instances with a comparatively large minimum spectral gap trajectory. Investigation showed that those MCP instances represent fully connected graphs/cliques, which is a trivial MCP and results in an omitted constraint, respectively penalty factor \(A,\) in the Ising formulation (see Equation (10)). Therefore, it cannot be optimized with a certain penalty factor ratio and is visualized as a horizontal trajectory in Figure 3(c).

MECP, KP, and BILP all show different clustering. As already seen in the preliminary raw data, MECP has instances with different trends in the clustering (cf. Figures 2(a) and 3(a)). Note that there is one cluster (label \(-1\)) with no slope and three clusters with upward (label 0 and 1) and downward (label 2) trends with an increasing penalty factor ratio. KP shows two clusters with both having a small upward trend with an increasing penalty factor ratio (cf. Figure 3(d)). Within the BILP clustering, all found clusters remain with the same minimum spectral gap size despite varying penalty factor ratios. Three clusters (labels \(-1\), 7, and 8) have a relatively small minimum spectral gap and a late location (80–100) in the anneal path, whereas the remaining clusters gained larger minimum spectral gaps within an early stage (20–40) in the anneal path (cf. Figure 3(f)). Investigation showed that BILP instances do not differ in the on- and off-diagonal values in the Ising formulation, respectively \(h_i\) and \(J_{ij}\) values of the Ising model. Therefore, no change in the ratio of penalty factors is possible (see Equation (15)).

In the next step, we investigate whether an NN is able to predict the penalty factor ratio associated with the largest minimum spectral gap, given the problem instance as input. Figure 4 shows the resulting RMSE and \(R^2\) coefficients while training our model, including the \(95\%\) confidence interval over 10 runs. Trivially, for KP, MVCP, MCP, and SPP, the \(R^2\) coefficient reaches around \(0.981 \pm 0.002\) to \(1.0 \pm 0.0\) at the end of training, since for these problems the best penalty ratio to choose is the maximal one (i.e., where factor A is set to the lowest possible value in our setting). Additionally, the RMSE converges to 0.0, which tells us that the distance between the predicted penalty factors made by our regression model and the actual best penalty factors in our setting is very small or null. Interestingly, training a good NN model for MECP and BILP is more difficult. The best achieved \(R^2\) coefficient was \(0.918 \pm 0.048\) and \(0.959 \pm 0.012\), respectively. This also reflects in a comparatively worse RMSE of \(0.091 \pm 0.048\) and \(0.072 \pm 0.001\). A possible reason for that can be observed in the corresponding clustering (see Figure 3). Whereas KP, MVCP, MCP, and SPP all have no or upward trends with an increasing penalty factor ratio, MECP correspondingly contains upward and downward trends, which leads to different penalty factors being best. We presume that this leads to a slightly worse performance in predicting the best penalty factors in our setting. Regarding BILP, no trajectory cluster has a slope, which leads to all penalty factor ratios being equally good w.r.t. the size of the minimum spectral gap. We assume that this results in a similar behavior to the MECP for the NN model.

Subsequently, we test whether choosing the best penalty factor ratio in our predefined specific range influences the approximation ratio on the D-Wave Advantage 4.1 system and whether using a NN for predicting the penalty factors is beneficial compared to choosing them arbitrarily from the valid penalty factor range.

For each given problem class, a set of 100 random problem instances is generated and a penalty ratio is inferred from the NN model. Next, the Ising model of those 100 problem instances is formulated (using the predicted penalty ratio) and run on the D-Wave machine. In Figure 5, the resulting approximation ratio of the model is shown as the red line including the \(95\%\) confidence interval over the 100 problem instances. Next, an experiment is started:. For 50 iterations, a random penalty ratio is sampled from the predefined range and the Ising model is formulated and sent to the D-Wave machine. The running best approximation ratio is kept (green line) with a \(95\%\) confidence interval. Note that even though the approximation ratio was calculated as stated Section 4.2.3, we used the average approximation ratio of the model as the baseline and therefore set it to 1.0, whereas the approximation ratio of the randomly sampled penalty factors was set relative to that of the model, to better visualize at what point the random sampling process reaches or surpasses the baseline.

For MVCP, MCP, and SPP, the NN model performs comparatively very well. In the case of MVCP, it takes 43 iterations for the random process to reach the model on average. In the case of SPP, it takes about 50 iterations to nearly reach the model, whereas in the case of MCP, the random process never reaches the model in the 50 samples that we generated. Of course, these three problems that work well are also the ones for which inference is easy (it consists of predicting the highest penalty value ratio). Regarding the other problem classes (MECP, KP, BILP), it takes on average just 5 iterations for the random process to outperform the model. However, it should be noted that the confidence interval of the model performance in these cases completely overlaps with the approximation ratio of the random process. Possible reasons can be found in the clustering results. Even though the trends of the minimum spectral gap trajectories in the KP clustering are similar to MVCP, MCP, and SPP, the performance of the predicted optimal penalty factors is comparatively bad on real hardware (cf. Figures 3 and 5). We believe this is due to the small minimum spectral gap interval of KP problems, which is in the range of 0.0 to 0.05 (cf. Figure 2(d)). We assume that those small differences in the minimum spectral gap size of KP Hamiltonians have no effect when executing it on D-Wave’s hardware. Regarding MECP, a similar problem can be observed. Although the minimum spectral gap interval is by a factor of 10 larger than the one of KP (cf. Figure 2(a)), the difference of the gap size between the worst and best penalty factor ratio of a problem instance is very small (low slope), which leads to the same assumption as for KP. Since in the BILP clustering no slope is present, no optimization can be achieved here in theory, but also practically, on real hardware.

Since in this experiment the running best approximation ratio is kept for the random process, we now compare the average performance of the NN model prediction against the average of all 50 random samples from the iterative process of the 100 problem instances per problem class, as presented in the previous paragraph. As a metric, the approximation ratio is used again. Table 2 shows the corresponding results. It is clear that the model on average outperforms at least \(12.9\% \pm 7.7\%\) and up to \(167.1\% \pm 29.8\%\) against the random process. Thus, from this perspective, the model should always be used for the six problems we analyzed here.

Furthermore, an interesting observation in terms of the number of optimal solutions could be made. Using the DBSCAN clustering that we found previously and grouping the clusters with no slope around the gap size of 0.0 (cf. Figure 3) and correlating it to each problem instance having either exactly one optimal solution or multiple, the PCC shows in general very high correlations. As seen in Table 3, MCP, MVCP, SPP, and BILP all have extremely strong (either positive or negative) correlations, with MECP still having a rather highly negative correlation. Only KP is an exception, with a comparatively lower correlation with the clustering. It is quite interesting that having exactly one optimal solution leads to minimum spectral gap trajectories that can be optimized well, whereas having two or more optimal solutions results in a flat minimum spectral gap trajectory with a gap size of around 0.0. This is the case for all problems with high correlations. It should be noted that correlating the clustering with the exact number of optimal solutions (e.g., with exactly three optimal solutions) reduces the overall PCC for all problems. This in turn shows that the exact number of optimal solutions is comparatively rather insignificant for the correlation.

In addition, we investigated the sparseness of the problem instances, to check if the trajectory clusters with their different trends also correlate with the sparseness of the problems. The sparseness of an problem instance (Ising matrix) was calculated as \(Sparseness=\frac{\#Elements-\#NonZeroElements}{\#Elements}\). Even though sparse problems are easier to embed into the sparse connectivity of D-Wave hardware and should therefore be easier to solve, since smaller (embedded) problems are produced, the minimum spectral gap trajectories and their trends did not correlate with the sparseness of the problems. The PCC was around 0 for each problem class. In general, each cluster (with upward, downward, and also the flat trends) of the problem classes contained problem instances with different sparseness.

Last, we investigate the minimum spectral gap trajectories when scaling the problem instances up and down in terms of variables. Since for all previous experiments the problem instances were restricted to 8 variables (due to the time-consuming computation of the eigenspectra), we now look into instances with 6, 8, and 10 variables for only the comparatively well-performing problem classes MVCP, MCP, and SPP. Note that for problem instances of size 6 and 8, we generated 1,000 instances per problem class, and for size 10, we generated only 500 instances per problem class. Figure 6 shows the composed DBSCAN clustering of problem instances of different size. Clusters of same problem size are identically colored. The achieved PCC (with optimized DBSCAN hyperparameters) was \(-0.900\), \(-0.992\), and \(-0.901\) for MVCP, MCP, and SPP, respectively. One can clearly see that the same trajectory trends occur regardless of the problem size. We therefore assume that the observed trends are size independent and will also hold for other problem sizes.