QPU Solvers: Configuration

When preparing your problem for submission, consider the analog nature of quantum computers, and take measures to increase the likelihood of finding good solutions.

This chapter presents guidance on the following topics:

Read-Anneal Cycles

The Solver Parameters section of the Getting Started with the D-Wave System guide explains why you should always request multiple read-anneal cycles for problems you submit to a D-Wave QPU solver.

As a guideline, equalizing your annealing and reading times is cost-effective. Increasing the number of reads can give better solutions, but improvements diminish over some problem-dependent number. If you were sampling independently from a normal distribution on energies, the sample minimum would decrease as the logarithm of the number of samples. If you were sampling from a lower-bounded distribution (e.g. containing a minimum-energy ground state), taking more samples beyond some number no longer improves the result.

Increasing both annealing time and num_reads increases your probability of success for a problem submission. The optimal combination of annealing time and number of reads to produce the best solution in a fixed time budget is problem dependent, and requires some experimentation.

Example

This example submits the same BQM, randomly generated using the dimod ran_r() class and identically embedded in the same QPU, first with num_read=10000 and an anneal time of 1 ms and then with num_reads=1000 and an anneal time of 50ms. The percentage of best solutions are compared.

>>> from dwave.system import DWaveSampler, FixedEmbeddingComposite, EmbeddingComposite
>>> qpu = DWaveSampler()
>>> import dimod
>>> import dwave.inspector
...
>>> bqm = dimod.generators.ran_r(5, 25)
>>> sampleset_1 = EmbeddingComposite(qpu).sample(bqm,
...                                              return_embedding=True,
...                                              answer_mode="raw,
...                                              num_reads=10000,
...                                              annealing_time=1)
>>> embedding = sampleset_1.info["embedding_context"]["embedding"]
>>> sampleset_50 = FixedEmbeddingComposite(qpu, embedding).sample(bqm,
...                                                               answer_mode="raw,
...                                                               num_reads=1000,
...                                                               annealing_time=50)
>>> print("Best solutions are {}% of samples.".format(
...       len(sampleset_1.lowest(atol=0.5).record.energy)/100))
Best solutions are 1.31% of samples.
>>> print("Best solutions are {}% of samples.".format(
...       len(sampleset_50.lowest(atol=0.5).record.energy)/10))
Best solutions are 4.0% of samples.

Figure 59 shows energy histograms from this example’s two submissions.

image

Fig. 59 Energy histograms for 10,000 samples with an anneal time of 1 ms (left) versus 1,000 samples with an anneal time of 50ms (right).

For this particular example, it makes sense to increase the anneal time at the expense of fewer anneal cycles (number of reads). That might not be true of other problems.

Spin-Reversal (Gauge) Transforms

Coupling \(J_{i,j}\) adds a small bias to qubits \(i\) and \(j\) due to leakage. This can become significant for chained qubits: see integrated control errors (ICE) 1. Additionally, qubits are biased to some small degree in one direction or another due to QPU imperfections.

Applying a spin-reversal transform can improve results by reducing the impact of such unintended biases. A spin-reversal transform does not alter the Ising problem; the transform simply amounts to reinterpreting spin up as spin down, and visa-versa, for a particular spin.

  • Changing too few spins leaves most errors unchanged, and therefore has little effect.
  • Changing too many spins means that most couplers connect spins that are both transformed, thus \(J_{i,j}\) does not change sign. As a result, some systematic errors associated with the couplers are unaffected.

Specify the number of spin-reversal transforms using the num_spin_reversal_transforms solver parameter. Note that increasing this parameter increases the total run time of the problem.

Example

This example solves Ocean’s AND gate example using illustratively long chains for two of the variables (for reference, you can embed this AND gate onto three qubits on a QPU). The first submission does not use spin-reversal transforms while the second does. An AND gate has four feasible states: \(x1, x2, out\) should take values \(000, 010, 100, 111\). The example prints the percentage of samples found for each of the feasible states of all lowest-energy samples with unbroken chains (typically this example also produces a small number of solutions with broken chains). Ideally, for a perfectly balanced QPU, feasible states would be found in equal numbers: [25 25 25 25] percent.

Note

The qubits selected below for chains are available on the particular Advantage QPU used for the example. Select a suitable embedding for the QPU you run examples on.

>>> from dwave.system import DWaveSampler, FixedEmbeddingComposite
>>> qpu = DWaveSampler()
>>> Q = {('x1', 'x2'): 1, ('x1', 'z'): -2, ('x2', 'z'): -2, ('z', 'z'): 3}
>>> embedding = {'x1': [2146, 2131, 2145, 2147, 3161, 3176, 3191, 3206, 3221,
...              3236,  3281, 3296, 3311, 3326], 'x2': [3251, 2071, 2086, 2101,
...              2116, 2161, 2176, 2191, 2206, 2221, 2236, 3250, 3252], 'z': [3266]}
...
>>> sampleset = FixedEmbeddingComposite(qpu, embedding).sample_qubo(Q, num_reads=5000)
>>> print("Feasible: {} percent. Time: {} ms".format(
...       sampleset.record.num_occurrences[:4]/50,
...       sampleset.info["timing"]["qpu_access_time"]))
Feasible: [16.9  17.02  8.76 56.84] percent. Time: 446041 ms
...
>>> sampleset_srt = FixedEmbeddingComposite(qpu, embedding).sample_qubo(Q,
...                       num_reads=5000,  num_spin_reversal_transforms=100)
>>> print("Feasible: {} percent. Time: {} ms".format(
...       sampleset_srt.record.num_occurrences[:4]/50,
...       sampleset_srt.info["timing"]["qpu_access_time"]))
Feasible: [18.76 25.14 29.6  25.48] percent. Time: 3023974 ms

Note that the submission using spin reversals produced more balanced solutions (the four feasible configurations for an AND gate are closer to being 25% of the lowest-energy samples with unbroken chains). Note too that the QPU access time increased from about half a second to about three seconds.

Further Information

Postprocessing

Postprocessing optimization and sampling algorithms provide local improvements with minimal overhead to solutions obtained from the quantum processing unit (QPU).

Ocean software provides postprocessing tools, and you can optionally run postprocessing online on D-Wave 2000Q and earlier systems.

Example: Broken-Chain Fixing

By default, Ocean embedding composites such as the EmbeddingComposite class fix broken chains.

This three-variable example ferromagnetically couples variable a, represented by a two-qubit chain, to two variables, b and c, that have opposing biases and are represented by one qubit each. Setting a chain strength that is smaller than the ferromagnetic coupling makes it likely for the chain to break.

Figure 60 shows the problem graph and its embedding in an Advantage QPU.

image

Fig. 60 The problem graph (left) and a particular embedding on a QPU (right), with a broken chain, as displayed by the problem inspector.

The first submission uses Ocean’s default postprocessing of chains to set a value for variable a; the second submission discards samples with broken chains.

Note

The qubits selected below are available on the particular Advantage QPU used for the example. Select a suitable embedding for the QPU you run examples on.

>>> from dwave.system import DWaveSampler, FixedEmbeddingComposite
>>> from dwave.embedding import chain_breaks
...
>>> qpu = DWaveSampler(solver={'topology__type': 'pegasus'})
>>> embedding={'a': [4755, 99], 'b': [69], 'c': [4785]}
...
>>> sampleset = FixedEmbeddingComposite(qpu, embedding=embedding).sample_ising(
...                                     {'b': +1, 'c': -1}, {'ab': -1, 'ac': -1},
...                                     chain_strength=0.8,
...                                     num_reads=1000)
>>> print(sampleset)
   a  b  c energy num_oc. chain_b.
0 +1 -1 +1   -2.0     672 0.333333
1 -1 -1 +1   -2.0      49      0.0
2 +1 +1 +1   -2.0     118      0.0
3 -1 -1 -1   -2.0      74      0.0
4 +1 -1 +1   -2.0      87      0.0
['SPIN', 5 rows, 1000 samples, 3 variables]
...
>>> sampleset = FixedEmbeddingComposite(qpu, embedding=embedding).sample_ising(
...                                     {'b': +1, 'c': -1}, {'ab': -1, 'ac': -1},
...                                     chain_strength=0.8,
...                                     num_reads=1000,
...                                     chain_break_method=chain_breaks.discard)
>>> print(sampleset)
   a  b  c energy num_oc. chain_.
0 -1 -1 +1   -2.0      60     0.0
1 +1 +1 +1   -2.0      79     0.0
2 -1 -1 -1   -2.0     142     0.0
3 +1 -1 +1   -2.0      77     0.0
['SPIN', 4 rows, 358 samples, 3 variables]

D-Wave 2000Q and Earlier Systems

For online postprocessing (postprocessing executed on the remote solver), you can choose:

  • No postprocessing (the default).
  • Optimization postprocessing to obtain a set of samples with lowest energy on a graph \(G\).
  • Sampling postprocessing to obtain a set of samples that correspond to a target Boltzmann distribution with inverse temperature \(\beta\) defined on the logical graph \(G\).

Postprocessing does not affect throughput—it works in tandem with the system and adds little overhead time (except for the last batch). However, if you have fewer than about 50 reads, postprocessing will likely not batch the outputs, so you do get overhead.

Further Information

Imprecision of Biases

Ising problems with high-precision parameters (\(h_i\) and \(J_{i,j}\)) present a challenge for quantum computers due to the finite precision available on \(\vc{h}\) and \(\vc{J}\). A problem may have lowest energy states that are sensitive to small variations in \(h\) or \(J\) while also requiring a large range of \(h\) or \(J\) values or high penalty values to enforce constraints on chains of qubits.

These are typically quantitative optimization problems rather than problems of a purely combinatorial nature (such as finding a subgraph with certain properties), where the number and connectivity of the qubits is more important than the weights, and problems for which near-optimal solutions are unacceptable. The solution’s quality depends on slight differences, at low-energy regions of the solution space, of the problem Hamiltonian as delivered to the QPU from its specification.

Example: Limiting Biases with Embedding

You can improve results by minimizing the range of on-QPU \(J\) or \(h\) values through embeddings.

For example, if a problem variable \(s_i\), which has the largest parameter value \(h_i\), is represented by qubits \(q_i^1, q_i^2\), and \(q_i^3\) having the same value in any feasible solution, \(h_i\) can be shared across the three qubits; i.e., \(h_i s_i \rightarrow (h_i/3)(q_i^1+q_i^2+q_i^3)\), reducing \(h_i\) by a factor of 3. In a similar way, coupling parameters \(J_{i,j}\) may also be shared.

In any embedding there may be multiple edges between chains of qubits representing problem variables. You can enhance precision (at the cost of using extra qubits) by sharing the edge weight across these edges.

>>> from dwave.system import DWaveSampler, EmbeddingComposite, FixedEmbeddingComposite
>>> import networkx as nx
>>> import dimod
>>> import random
>>> import dwave.inspector
...
>>> # Create a 5-variable problem with one outsized bias
>>> G = nx.generators.small.bull_graph()
>>> for edge in G.edges:
...     G.edges[edge]['quadratic'] = random.choice([1,-1])
>>> for node in range(max(G.nodes)):
...    G.nodes[node]['linear'] = random.choice([0.1,-0.1])
>>> G.nodes[max(G.nodes)]['linear'] = 10
>>> bqm = dimod.from_networkx_graph(G,
...                                 vartype='SPIN',
...                                 node_attribute_name ='linear',
...                                 edge_attribute_name='quadratic')
>>> # Submit the problem to a QPU solver
>>> qpu = DWaveSampler(solver={'topology__type': 'pegasus'})
>>> sampleset = EmbeddingComposite(qpu).sample(bqm, num_reads=1000, return_embedding=True)

Figure 63 shows the embedded problem with the large-biased variable represented by qubit 3999.

image

Fig. 62 An embedded problem with one large-biased variable.

>>> embedding = dict(sampleset.info["embedding_context"]["embedding"])
>>> embedding[4]
(3999,)
>>> embedding[4] = [3999, 1715, 1730, 1745]
>>> sampleset = FixedEmbeddingComposite(qpu, embedding).sample(bqm, num_reads=1000)

Figure 63 shows the embedded problem with the large-biased variable represented by four chained qubits.

image

Fig. 63 A large-biased variable represented by four chained qubits.

Example: Limiting Biases by Simplifying the Problem

In problems with interaction \(h_i s_i\), where \(h_i>0\) is much larger than all other problem parameters, it is likely that in low-energy states, \(s_i=-1\) (\(2h\) lower in energy than \(s_i=+1\)). Generally, you may be able to identify, in polynomial time, a subset of variables that always take the same value in the ground state. You can then eliminate such variables from the problem.

Consider preprocessing problems to determine whether certain variable values can be inferred. There is little overhead in attempting to simplify every problem before sending it to the QPU.

The code below preprocesses the problem of the previous section, which has a single outsized value for variable 4.

>>> from dwave.system import DWaveSampler, EmbeddingComposite
>>> import networkx as nx
>>> import dimod
>>> import random
>>> from dimod.reference.composites import RoofDualityComposite
...
>>> # Create a 5-variable problem with one outsized bias
>>> G = nx.generators.small.bull_graph()
>>> for edge in G.edges:
...     G.edges[edge]['quadratic'] = random.choice([1,-1])
>>> for node in range(max(G.nodes)):
...    G.nodes[node]['linear'] = random.choice([0.1,-0.1])
>>> G.nodes[max(G.nodes)]['linear'] = 10
>>> bqm = dimod.from_networkx_graph(G,
...                                 vartype='SPIN',
...                                 node_attribute_name ='linear',
...                                 edge_attribute_name='quadratic')
>>> # Preprocess and submit to a QPU solver
>>> sampler_pp = RoofDualityComposite(EmbeddingComposite(DWaveSampler()))
>>> sampleset = sampler_pp.sample(bqm, num_reads=1000, return_embedding=True)
...
>>> # Check which variables were fixed in preprocessing
>>> print(set(bqm.variables) -
...       set(sampleset.info["embedding_context"]["embedding"].keys()))
{4}

The problem submitted to the QPU has had the value of variable 4 fixed by the RoofDualityComposite composite.

Further Information

Problem Scale

In general, use the full range of \(h\) and \(J\) values available for the QPU when submitting a problem.

Ocean software’s default enabling of the auto_scale solver parameter automatically scales problems to make maximum use of the available ranges. However, if you are using the extended_j_range parameter, autoscaling is disabled.

Example

This example uses a single outsized bias to “squash” the scale of biases available for the remaining variables. Figure 64 shows the (native) BQM embedded on a D-Wave 2000Q QPU: all 128 variables have linear coefficients of 0 or 1 except for one variable with a linear coefficient of -100. (In practice your problem may have a minority of variables that have significantly different values from the majority.)

image

Fig. 64 A BQM with a single large-biased qubit that reduces the range of qubit biases available for the remaining linear biases.

>>> import dimod
>>> import networkx as nx
>>> from dwave.system import DWaveSampler
>>> import dwave_networkx as dnx
...
>>> # Create a native C4 problem with one outsized bias
>>> qpu = DWaveSampler(solver={'topology__type': 'chimera'})
>>> c16_working = dnx.chimera_graph(16, node_list=qpu.nodelist, edge_list=qpu.edgelist)
>>> c4_nodes = [v for v, d in c16_working.nodes(data=True) if all(x < 4 for x in d['chimera_index'][:2])]
>>> c4 = c16_working.subgraph(c4_nodes)
>>> bqm = dimod.generators.randint(c4, "SPIN")
>>> bqm.set_linear(list(bqm.linear.keys())[0], -100)
...
>>> # Submit with and without the outsized bias
>>> sampleset = qpu.sample(bqm, num_reads=1000)
>>> sampler_fixed = dimod.composites.FixedVariableComposite(qpu)
>>> sampleset_fixed = sampler_fixed.sample(bqm, fixed_variables={0: 1}, num_reads=1000)

The two figures below show the energies of returned solutions:

  • Figure 65 is the BQM with the outsized bias.
  • Figure 66 is the (updated) BQM without the outsized bias.
image

Fig. 65 Energies of returned solutions for the original BQM with the outsized bias.

image

Fig. 66 Energies of returned solutions for the BQM with the variable that has an outsized bias fixed.

When the original BQM is embedded on the QPU, the problem range is scaled by 50:

>>> print(min(bqm.linear.values())//min(qpu.properties["h_range"]))
50.0

The qubit biases of all but one variable and the coupler strengths are either 0 or 0.02; that is, over 99% of the qubit biases are in just 1% of the available range of \(h\) values ([-2.0, 2.0] for the QPU on which this example was run).

In a fixed BQM, most variables keep the original coefficient values except for variables connected to the fixed variable. For this execution of the example, two connected variables’ coefficients are changed to a value of 2 to account for the fixed variable. The fixed BQM is embedded with scaling of 1:

>>> bqm.fix_variable(0, 1)
>>> set(bqm.quadratic.values()) | set(bqm.linear.values())
{0.0, 1.0, 2.0}
>>> print(max(bqm.linear.values())//max(qpu.properties["h_range"]))
1.0

The qubit biases of most variables and the coupler strengths are either 0 or 1, and two qubits have biases of 2; that is, over 99% of the qubit biases are in 25% of the available range of \(h\) values for this QPU.

For this artificial example, you can easily fix the variables with bias 2 as well. Doing so reduces the values of biases to three:

>>> set(bqm.quadratic.values()) | set(bqm.linear.values())
{-1.0, 0.0, 1.0}

Now all the BQM’s quadratic coefficients of 0 or 1 are in 50% of the available range of \(J\) values for this particular QPU ([-1.0, 1.0]) and the linear coefficients of 0, 1, and -1 are in 50% of the available range of \(h\) values ([-2.0, 2.0]).

Figure 67 is the energies of returned solutions for the BQM with all variables with coeffcients that are not either 0 or 1 fixed.

image

Fig. 67 Energies of returned solutions for the BQM with all variable that do not have values of 0 or 1 fixed.

Annealing Schedule

Note

Variations on the global anneal schedule are not supported on D-Wave 2X and earlier systems. For D-Wave 2000Q and Advantage systems, the maximum number of points in a schedule is system-dependent.

Some types of problems benefit from the introduction of a pause or a quench at some point in the anneal schedule. A pause dwells for some time at a particular anneal fraction; a quench abruptly terminates the anneal within a few hundred nanoseconds of the point specified.

This degree of control over the global annealing schedule also enables closer study the quantum annealing algorithm.

Pause and Quench

A pause can be a useful diagnostic tool for instances with a small perturbative anticrossing. While pauses early or late in the anneal have no effect, a pause near the expected perturbative anticrossing produces a large increase in the ground-state success rate.

If a quench is fast compared to problem dynamics, then the distribution of states returned by the quench can differ significantly from that returned by the standard annealing schedule. The probability of obtaining ground state samples depends on when in the anneal the quench occurs, with later quenches more likely to obtain samples from the ground state.

Supply the scheduling points using the anneal_schedule solver parameter.

Reverse Anneal

Reverse annealing enables the use of quantum annealing as a component in local search algorithms to refine classical states. Examples of using this feature include Quantum Boltzmann sampling, tunneling rate measurements, and relaxation rate measurements.

Examples

Figure 68 shows embedded in a D-Wave 2000Q[1] QPU a 16-qubit system, which was studied in a nature article. It has an energy gap of 4 between the classical ground state and excited states.

[1]This problem can be similarly embedded in an Advantage QPU but the embedding is simpler on a Chimera structure, as shown in the code below.
image

Fig. 68 A 16-qubit system with an energy gap of 4 between the classical ground state and excited states embedded in a D-Wave 2000Q QPU.

The following code shows how varying the anneal schedule can increase the probability of finding ground states. First, the problem is embedded onto a QPU such that each problem qubit is represented by a single qubit on the QPU.

>>> from minorminer import find_embedding
>>> from dwave.system import DWaveSampler, FixedEmbeddingComposite
>>> import numpy as np
...
>>> # Configure the problem structure
>>> h = {0: 1.0, 1: -1.0, 2: -1.0, 3: 1.0, 4: 1.0, 5: -1.0, 6: 0.0, 7: 1.0,
...      8: 1.0, 9: -1.0, 10: -1.0, 11: 1.0, 12: 1.0, 13: 0.0, 14: -1.0, 15: 1.0}
>>> J = {(9, 13): -1, (2, 6): -1, (8, 13): -1, (9, 14): -1, (9, 15): -1,
...      (10, 13): -1, (5, 13): -1, (10, 12): -1, (1, 5): -1, (10, 14): -1,
...      (0, 5): -1, (1, 6): -1, (3, 6): -1, (1, 7): -1, (11, 14): -1,
...      (2, 5): -1, (2, 4): -1, (6, 14): -1}
...
>>> # Find an embedding with a single QPU qubit representing each problem qubit
>>> qpu_chimera = DWaveSampler(solver={'topology__type': 'chimera'})
>>> max_units = 16*16
>>> reads = 1000
>>> nodes = {s: sorted(list(set(list(range(s, s + max_units*8, 8))) &
...                    set(qpu_chimera.nodelist))) for s in h.keys()}
>>> embedding = find_embedding(J.keys(),
...                            qpu_chimera.edgelist,
...                            restrict_chains=nodes)
>>> sampler = FixedEmbeddingComposite(qpu_chimera, embedding)

Print the percentage of ground states for a 100 ms anneal:

>>> sampleset = sampler.sample_ising(h, J, num_reads=reads, answer_mode='raw',
...                                  annealing_time=100)
>>> counts = np.unique(sampleset.record.energy.reshape(reads,1), axis=0,
...                    return_counts=True)[1]
>>> print("{}% of samples were best energy {}.".format(100*counts[0]/sum(counts),
...       sampleset.first.energy))
0.8% of samples were best energy -20.0.

Print the percentage of ground states for an anneal with a 100 ms pause:

>>> anneal_schedule=[[0.0, 0.0], [40.0, 0.4], [140.0, 0.4], [142, 1.0]]
>>> sampleset = sampler.sample_ising(h, J, num_reads=reads, answer_mode='raw',
...                                  anneal_schedule=anneal_schedule)
>>> counts = np.unique(sampleset.record.energy.reshape(reads,1), axis=0,
...                    return_counts=True)[1]
>>> print("{}% of samples were best energy {}.".format(100*counts[0]/sum(counts),
...       sampleset.first.energy))
3.4% of samples were best energy -20.0.

Print the percentage of ground states for a reverse anneal (of almost 100 ms):

>>> reverse_schedule = [[0.0, 1.0], [5, 0.45], [99, 0.45], [100, 1.0]]
>>> initial = dict(zip(sampleset.variables, sampleset.record[int(reads/2)].sample))
>>> reverse_anneal_params = dict(anneal_schedule=reverse_schedule,
...                              initial_state=initial,
...                              reinitialize_state=True)
>>> sampleset = sampler.sample_ising(h, J, num_reads=reads, answer_mode='raw',
...                                  **reverse_anneal_params)
>>> counts = np.unique(sampleset.record.energy.reshape(reads,1), axis=0,
...                    return_counts=True)[1]
>>> print("{}% of samples were best energy {}.".format(100*counts[0]/sum(counts),
...       sampleset.first.energy))
20.8% of samples were best energy -20.0.

Further Information

Anneal Offsets

Note

Anneal offsets are not supported on D-Wave 2X and earlier systems. Before using this feature, query the solver properties using SAPI calls to determine whether it is supported and, if so, to obtain the available tuning ranges per qubit.

Anneal offsets may improve results for problems in which the qubits have irregular dynamics for some easily determined reason; for example, if a qubit’s final value does not affect the energy of the classical state, you can advance it (with a positive offset) to reduce quantum bias in the system.

Anneal offsets can also be useful in embedded problems with varying chain length: longer chains may freeze out earlier than shorter ones, which means that at an intermediate point in the anneal, some variables act as fixed constants while others remain undecided. If, however, you advance the anneal of the qubits in the shorter chains, they freeze out earlier than they otherwise would. The correct offset will synchronize the annealing trajectory of the shorter chains with that of the longer ones.

If you decide that offsetting anneal paths might improve results for a problem, your next task is to determine the optimal value for the qubits you want to offset. As a general rule, if a qubit is expected to be subject to a strong effective field relative to other qubits, delay its anneal with a negative offset. The ideal offset magnitudes are likely to be the subject of trial and error, but expect that the appropriate offsets for two different qubits in the same problem to be within 0.2 normalized offset units of each other.

Supply the array of offsets for the qubits in the system using the anneal_offsets solver parameter with a length equal to the num_qubits property.

Example: 3-Qubit System

This example is a 3-qubit looks at a system that has a ground state, \(1, 1, 1\), separated from its two closest excited states, \(-1, -1, -1\) and \(-1, -1, 1\), by a small energy gap compared to its remaining excited states. These two first excited states have the same energy and differ by a single flip of qubit 2; consequently, the superposition of these two states is dominant early in the anneal. Figure 69 shows the problem and a possible embedding in one particular Advantage QPU.

image

Fig. 69 A three-qubit system with a small energy gap between the ground state and first two excited states.

dimod’s ExactSolver shows the energies of the ground state, first two excited states, and remaining states:

>>> from dimod import ExactSolver
>>> h = {0: 0, 1: 0.9, 2: -1}
>>> J = {(0, 1): -1, (1, 2): -1}
>>> print(ExactSolver().sample_ising(h, J))
   0  1  2 energy num_oc.
5 +1 +1 +1   -2.1       1
0 -1 -1 -1   -1.9       1
7 -1 -1 +1   -1.9       1
4 -1 +1 +1   -0.1       1
1 +1 -1 -1    0.1       1
6 +1 -1 +1    0.1       1
2 +1 +1 -1    1.9       1
3 -1 +1 -1    3.9       1
['SPIN', 8 rows, 8 samples, 3 variables]

The D-Wave system used for this example is an Advantage that has couplers between active qubits 30, 31, and 2940. Select a suitable embedding for the QPU you run examples on.

>>> from dwave.system import FixedEmbeddingComposite, DWaveSampler
>>> qpu = DWaveSampler()
>>> embedding = {0: [31], 1: [30], 2: [2940]}
>>> sampler = FixedEmbeddingComposite(qpu, embedding)
>>> print(qpu.properties['anneal_offset_ranges'][2940])
[-0.7012257815714587, 0.6717794151250857]

For the default anneal offset of qubit 2, this particular run of 1000 samples, successfully returned the problem’s ground state about one third of the time, and likewise each of the two first excited states a third of the time:

>>> sampleset = sampler.sample_ising(h, J, num_reads=1000)
>>> print(sampleset)
   0  1  2 energy num_oc. chain_.
0 +1 +1 +1   -2.1     386     0.0
1 -1 -1 +1   -1.9     276     0.0
2 -1 -1 -1   -1.9     338     0.0
['SPIN', 3 rows, 1000 samples, 3 variables]

Appying a positive offset to qubit 2 causes it to freeze a bit earlier in the anneal than qubits 0 and 1. Consequently, the superposition of the two lowest excited states, \(-1, -1, -1\) and \(-1, -1, 1\), no longer dominates, and the ground state is found much more frequently.

>>> offset = [0]*qpu.properties['num_qubits']
>>> offset[2940]=0.2
>>> sampleset = sampler.sample_ising(h, J, num_reads=1000, anneal_offsets=offset)
>>> print(sampleset)
   0  1  2 energy num_oc. chain_.
0 +1 +1 +1   -2.1     979     0.0
1 -1 -1 +1   -1.9       7     0.0
2 -1 -1 -1   -1.9      13     0.0
3 -1 +1 +1   -0.1       1     0.0
['SPIN', 4 rows, 1000 samples, 3 variables]

Example: 16-Qubit System

The example problem of the Annealing Schedule section improved solutions for a 16-qubit system, shown in Figure 68 embedded in a D-Wave 2000Q QPU, which was studied in a nature article, and has an energy gap of 4 between the classical ground state and excited states. The eight “outer” qubits, which are coupled to only one other qubit, enable single flips that produce dominant superpositions of excited states in the anneal (these are small-gap anticrossings), reducing the likelihood of finding the ground state.

First, the problem is embedded onto a QPU such that each problem qubit is represented by a single qubit on the QPU.

>>> from minorminer import find_embedding
>>> from dwave.system import DWaveSampler, FixedEmbeddingComposite
>>> import numpy as np
...
>>> # Configure the problem structure
>>> h = {0: 1.0, 1: -1.0, 2: -1.0, 3: 1.0, 4: 1.0, 5: -1.0, 6: 0.0, 7: 1.0,
...      8: 1.0, 9: -1.0, 10: -1.0, 11: 1.0, 12: 1.0, 13: 0.0, 14: -1.0, 15: 1.0}
>>> J = {(9, 13): -1, (2, 6): -1, (8, 13): -1, (9, 14): -1, (9, 15): -1,
...      (10, 13): -1, (5, 13): -1, (10, 12): -1, (1, 5): -1, (10, 14): -1,
...      (0, 5): -1, (1, 6): -1, (3, 6): -1, (1, 7): -1, (11, 14): -1,
...      (2, 5): -1, (2, 4): -1, (6, 14): -1}
...
>>> # Find an embedding with a single QPU qubit representing each problem qubit
>>> qpu_chimera = DWaveSampler(solver={'topology__type': 'chimera'})
>>> max_units = 16*16
>>> reads = 5000
>>> nodes = {s: sorted(list(set(list(range(s, s + max_units*8, 8))) &
...                    set(qpu_chimera.nodelist))) for s in h.keys()}
>>> embedding = find_embedding(J.keys(),
...                            qpu_chimera.edgelist,
...                            restrict_chains=nodes)
>>> sampler = FixedEmbeddingComposite(qpu_chimera, embedding)

Figure 70 shows embedded in a D-Wave 2000Q[2] QPU.

[2]This problem can be similarly embedded in an Advantage QPU but the embedding is simpler on a Chimera structure, as shown in the code.
image

Fig. 70 A 16-qubit system with an energy gap of 4 between the classical ground state and excited states embedded in a D-Wave 2000Q QPU. One of the eight “outer” qubits, problem qubit 4, which is embedded here to QPU qubit 172, is highlighted.

The following code samples 5000 times, with an anneal time of 100 ms and no annealing offsets, and prints the percentage of ground states found.

>>> offset = [0]*qpu_chimera.properties['num_qubits']
>>> sampleset = sampler.sample_ising(h, J, num_reads=reads, answer_mode='raw',
...                                  annealing_time=100,
...                                  anneal_offsets=offset)
>>> counts = np.unique(sampleset.record.energy.reshape(reads,1), axis=0,
...                    return_counts=True)[1]
>>> print("{}% of samples were best energy {}.".format(100*counts[0]/sum(counts),
...       sampleset.first.energy))
1.18% of samples were best energy -20.0.

The minimum range for a positive anneal offset for any of the eight “outer” qubits of this example’s particular embedding can be found as below:

>>> ao_range = qpu_chimera.properties["anneal_offset_ranges"]
>>> ao_step = qpu_chimera.properties["anneal_offset_step"]
>>> outer8 = [172, 180, 168, 176, 171, 179, 175, 183]
>>> print(min([ao_range[outer][1] / ao_step for outer in outer8]))
16.12534800734681

The following code samples 5000 times, with an anneal time of 100 ms and annealing offsets set for the eight “outer” qubits, and prints the percentage of ground states found.

>>> for outer in outer8:
...     offset[outer] = 10*ao_step
>>> sampleset = sampler.sample_ising(h, J, num_reads=reads, answer_mode='raw',
...                                  annealing_time=100,
...                                  anneal_offsets=offset)
>>> counts = np.unique(sampleset.record.energy.reshape(reads,1), axis=0,
...                    return_counts=True)[1]
>>> print("{}% of samples were best energy {}.".format(100*counts[0]/sum(counts),
...       sampleset.first.energy))
32.2% of samples were best energy -20.0.

Further Information

Energy Gap

There are strategies for increasing the gap between ground and excited states during the anneal. For example, different choices of constraints when reformulating a CSP as a QUBO affect the gap.

Consider also the differences between maximizing the gap versus creating a uniform gap.

Example

This example formulates a 3-bit parity check[3] as an Ising model embedded in three quarters of the qubits (6 of 8) of a single Chimera unit cell in two ways, the first with an energy gap of 1 and the second with an energy gap of 2.

[3]Feasible states are those states where the number of spin-up values for three Ising variables, \(s_1, s_2, s_3\), is even.

Figure 71 shows penalty model

\[0.5(s_1 + s_2 + s3) -a + 0.5(s_1s_2 + s_1s_3 + s_2s_3) - as_1 - as_2 - as_3\]

with one auxiliary variable, \(a\), and an energy gap of 1 embedded in a D-Wave 2000Q QPU.

image

Fig. 71 3-bit parity check as an Ising model with one auxiliary variable embedded in a D-Wave 2000Q QPU.

Figure 72 shows penalty model

\[-a_1 +a_2 - a_3 + s_1a_1 + s_1a_2 + s_1a_3 - s_2a_1 + s_2a_2 + s_2a_3 + s_3a_1 + s_3a_2 - s_3a_3\]

with three auxiliary variables, \(a_1, a_2, a_3\), and an energy gap of 2 embedded in a D-Wave 2000Q QPU.

image

Fig. 72 3-bit parity check as an Ising model with three auxiliary variables embedded in a D-Wave 2000Q QPU.

The following code runs on the working graph of a particular D-Wave 2000Q QPU; select suitable embeddings for the QPU you run examples on.

>>> from dwave.system import DWaveSampler, FixedEmbeddingComposite
>>> from dimod import ExactSolver
>>> qpu_chimera = DWaveSampler(solver={'topology__type': 'chimera'})

Solve the first model (energy gap of 1):

>>> h = {'s1': 0.5, 's2': 0.5, 's3': 0.5, 'a': -1}
>>> J = {('s1', 's2'): 0.5, ('s1', 's3'): 0.5, ('s2', 's3'): 0.5,
...      ('s1', 'a'): -1, ('s2', 'a'): -1, ('s3', 'a'): -1}
>>> embedding = {'s1': [0, 4], 's2': [1, 5], 's3': [2], 'a': [6]}
>>> sampler = FixedEmbeddingComposite(qpu_chimera, embedding)
>>> sampleset = sampler.sample_ising(h, J, num_reads=1000)
>>> print(sampleset)
    a s1 s2 s3 energy num_oc. chain_.
0  -1 -1 -1 -1   -2.0     313     0.0
1  +1 -1 +1 +1   -2.0     196     0.0
2  +1 +1 -1 +1   -2.0     216     0.0
3  +1 +1 +1 -1   -2.0     254     0.0
4  +1 +1 +1 +1   -1.0       2     0.0
5  -1 -1 -1 +1   -1.0       5     0.0
6  +1 -1 +1 -1   -1.0       2     0.0
7  -1 -1 +1 -1   -1.0       3     0.0
8  -1 +1 -1 -1   -1.0       5     0.0
9  +1 +1 -1 -1   -1.0       3     0.0
10 +1 +1 +1 +1   -1.0       1    0.25
['SPIN', 11 rows, 1000 samples, 4 variables]

Solve the second model (energy gap of 2):

>>> h = {'a1': -1, 'a2': 1, 'a3': -1}
>>> J = {('s1', 'a1'): 1, ('s1', 'a2'): 1, ('s1', 'a3'): 1,
...      ('s2', 'a1'): -1, ('s2', 'a2'): 1, ('s2', 'a3'): 1,
...      ('s3', 'a1'): 1, ('s3', 'a2'): 1, ('s3', 'a3'): -1}
>>> embedding = {'s1': [0], 's2': [1], 's3': [2], 'a1': [4], 'a2': [5], 'a3': [6]}
>>> sampler = FixedEmbeddingComposite(qpu_chimera, embedding)
>>> sampleset = sampler.sample_ising(h, J, num_reads=1000)
>>> print(sampleset)
  a1 a2 a3 s1 s2 s3 energy num_oc. chain_.
0 -1 -1 +1 +1 -1 +1   -6.0     167     0.0
1 +1 -1 +1 -1 +1 +1   -6.0      69     0.0
2 +1 -1 -1 +1 +1 -1   -6.0     445     0.0
3 +1 +1 +1 -1 -1 -1   -6.0     317     0.0
4 +1 +1 +1 -1 -1 +1   -4.0       1     0.0
5 -1 -1 -1 +1 +1 +1   -4.0       1     0.0
['SPIN', 6 rows, 1000 samples, 6 variables]

For this simple example, both formulations produce a high percentage of ground states (states where the number of spin-up values for variables \(s_1, s_2, s_3\) is even). Had this been part of a more complex problem, you might have needed to weigh the benefit of a larger energy gap against other considerations, such as a larger number of ancillary variables.

Further Information

  • [Bia2014] discusses constructing a penalty function for a given constraint with the largest possible gap, subject to bounds on the supported \(h\)s and \(J\)s.
  • [Pud2014] and [Pud2015] discuss error suppression techniques using auxiliary qubits and the energy gap.

Neighbor Interactions

The dynamic range of \(h\) and \(J\) values may be limited by ICE. Instead of finding low-energy states to an optimization problem defined by \(h\) and \(J\), the QPU solves a slightly altered problem that can be modeled as

\[E^{\delta}_{ising} ({\bf s}) = \sum_{i=1}^N (h_i + \delta h_i ) s_i + \sum_{i=1}^N \sum_{j=i+1}^{N} (J_{i,j} + \delta J_{i,j} ) s_i s_j,\]

where the ICE errors \(\delta h_i\) and \(\delta J_{i,j}\) depend on \(h_i\) and on the values of all incident couplers \(J_{i,j}\) and neighbors \(h_j\), as well as their incident couplers \(J_{j,k}\) and next neighbors \(h_k\). For example, if a given problem is specified by \((h_1 = 1 , h_2 = 1, J_{1,2} = -1)\), the QPU might actually solve the problem \((h_1 = 1.01, h_2 = 0.99, J_{1,2} = -1.01)\).

Changing a single parameter in the problem might change all three error terms, altering the problem in different ways.

Further Information

[Har2010] discusses how applied \(h\) bias leaks from spin \(i\) to its neighboring spins.