How to sample 5000 variables with EmbeddingComposite.sample_qubo()?

[DevelopDaily]

https://github.com/dwavesystems/dwave-system/blob/32e9065cddeb123106b43d947877043c1a2ccc78/dwave/system/composites/embedding.py#L47

This code produces a correct solution. It is basically a very simple Max-Cut Problem of 256 vertices with a known solution.

num_vars = 256

graph_list = list()
for i in range(num_vars-1):
    graph_list.append((i, i+1))
    
G = networkx.Graph()
G.add_edges_from(graph_list)

Q = defaultdict(int)
for i, j in G.edges:
    Q[(i,i)]+= -1
    Q[(j,j)]+= -1
    Q[(i,j)]+= 2
    
sampler = EmbeddingComposite(DWaveSampler() )    
response = sampler.sample_qubo(Q, chain_strength=4, num_reads=10000)

Unfortunately, however, as soon as more vertices, say, num_vars = 1024, are used, the program will never produce the theoretically correct result, which is supposed to have num_vars - 1 cut edges, and num_vars/2 vertices in each set.

Could you please help me reach the goal of as many vertices as possible, say, 5000? After all, it is claimed that D-Wave hardware is capable of handling more than 5000 variables.

Thanks.

[arcondello]

Hi @DevelopDaily, this is an interesting problem. You're essentially testing the maximum antiferromagnetic chain length. it's definitely not true that the hardware can solve all problems sent to it, and in some sense this is a particularly hard problem for the hardware relative to classically.

To explain why, it's probably best to start with the best way to solve it classically, which is to start at one end of the path and then simply alternate values along the chain. This can be solved in O(n) time.

However, classical BQM samplers like dwave-neal and dwave-tabu don't know this trick. They randomly select a starting position and test various changes until they find a good solution. So imagine you have a chain like 0-1-0-1-1-0-1-0-0-1-0-1 you can see that this chain is made up of four valid "sub chains" but the middle subchain is flipped. So in order for a classical solver to "fix it", it would need to flip the entire sub chain, which is difficult for it to do.

Now, on the quantum computer, we have a semi-similar situation. Different qubits can freeze out at different times. So part way through the anneal you can get something like 0-1-s-s-1-0-1-s-s-s-0-1 where s represents a variable still in superposition. The different frozen-out sub-chains can be valid among themselves, but not valid in relation to each other. Because frozen out qubits can only flip thermally, it is difficult to "fix".

I would be curious for your problem, if you're getting runs of valid chains in the solution. Maybe something like

def runs(sampleset):
     sample = sampleset.first.sample
 
     if not sample:
         return []
     chunks = [1]
     current_value = sample[0]
     for v in range(len(sample)):
         if sample[v] == current_value:  # not a cut, so make a new chunk
             chunks.append(1)
         else:
             current_value = sample[v]
             chunks[-1] += 1
     return chunks

would tell you how long each of the valid segments are. Note that the above function assumes that your path graph is labelled [0,n).

Also, FWIW, a simpler way to formulate this problem is

G = nx.path_graph(num_variables)

h = {v: 0 for v in G.nodes}
J = {(u, v): -1 for (u, v) in G.edges}

sampler = EmbeddingComposite(DWaveSampler() )    
sampleset = sampler.sample_ising(h, J, chain_strength=4, num_reads=10000)

[boothby]

It's worth noting that finding a path of length N in an N-node graph is generally quite difficult (this also known as the Hamiltonian Path Problem). I would expect minorminer to have a very hard time with this -- I got lucky and found a path of length 4000 on my first try, but 5000 is going to be very difficult for minorminer. If you're really curious, the first place to look would probably be longest path heuristics -- for example SageMath [1] has both heuristic (fast) and exact (slow) algorithms. Alternatively, you can try a subgraph solver such as Glasgow [2].

However, as Alex notes, this class of problem will be quite difficult for the QPU as well -- running such long paths may be of dubious value.

[1] https://doc.sagemath.org/html/en/reference/graphs/sage/graphs/generic_graph.html#sage.graphs.generic_graph.GenericGraph.longest_path [2] https://github.com/ciaranm/glasgow-subgraph-solver

[DevelopDaily]

Using your function, I print out the lowest energy and the chunks. Here are the results:

10 independent runs of: sampler.sample_qubo(Q, chain_strength=4, num_reads=10)

energy:  -1011.0   chunks: [1, 49, 12, 85, 69, 54, 32, 32, 230, 108, 197, 53, 49, 54]
energy:  -1011.0   chunks: [1, 160, 52, 21, 9, 142, 211, 11, 132, 92, 28, 30, 51, 85]
energy:  -1009.0   chunks: [1, 95, 54, 46, 12, 30, 14, 122, 98, 142, 197, 30, 105, 20, 28, 31]
energy:  -1011.0   chunks: [1, 58, 148, 142, 104, 101, 65, 43, 70, 3, 42, 73, 24, 151]
energy:  -1012.0   chunks: [1, 63, 37, 12, 141, 147, 50, 10, 59, 49, 132, 264, 60]
energy:  -1011.0   chunks: [1, 27, 29, 63, 67, 16, 185, 38, 97, 145, 46, 206, 58, 47]
energy:  -1009.0   chunks: [1, 28, 30, 56, 98, 56, 85, 37, 71, 68, 74, 187, 10, 150, 25, 49]
energy:  -1013.0   chunks: [1, 10, 88, 144, 309, 52, 155, 92, 51, 21, 92, 10]
energy:  -1009.0   chunks: [1, 67, 273, 44, 62, 6, 1, 16, 80, 33, 55, 54, 37, 50, 59, 187]
energy:  -1009.0   chunks: [1, 36, 92, 182, 52, 223, 25, 87, 24, 12, 40, 74, 33, 70, 24, 50]

2 independent runs of: sampler.sample_qubo(Q, chain_strength=4, num_reads=100)

energy:  -1014.0   chunks: [1, 63, 47, 57, 207, 72, 19, 309, 129, 82, 39]
energy:  -1015.0   chunks: [1, 68, 174, 236, 133, 114, 44, 96, 88, 71]

2 independent runs of: sampler.sample_qubo(Q, chain_strength=4, num_reads=10000)

energy:  -1018.0   chunks: [1, 163, 276, 106, 405, 44, 30]
energy:  -1017.0   chunks: [1, 125, 148, 92, 46, 239, 155, 219]

It is pleasantly surprising that the energy goes lower as the num_reads goes higher. The D-Wave API does not allow me to go beyond 10,000. Would you tweak the num_reads internally to, say, 1,000,000? The annealing may reach the theoretical lowest energy. How do you think?

By the way, thanks for your advice on the ising formulation, which I haven't learnt yet. Somehow, your particular ising model does not produce the same result as my model intends to. My graph basically looks like a long string of 1024 pearls.

[arcondello]

See also https://github.com/dwavesystems/dwave-system/issues/106