Use Boolean Hamiltonian Gate for QAOA

examples/examples_test.py

@@ -97,7 +97,7 @@ def test_example_heatmaps():
 
 
 def test_example_runs_qaoa():
-    examples.qaoa.main(repetitions=10, maxiter=5)
+    examples.qaoa.main(repetitions=1, maxiter=1, p=1)
 
 
 def test_example_runs_quantum_teleportation():

examples/qaoa.py

@@ -1,155 +1,162 @@
 """Runs the Quantum Approximate Optimization Algorithm on Max-Cut.
 
-=== EXAMPLE OUTPUT ===
+=== EXAMPLE CIRCUIT ===
+  x0    x1    x2     x3
+  │     │     │      │
+  H     H     H      H
+  │     │     │      │
+  │     │     Y^-0.5 Y^-0.5
+  │     │     │      │
+  │     @─────@      │
+  │     │     │      │
+  │     │     Y^0.5  │
+  │     │     │      │
+  │     │     Rz(0)  │
+  │     │     │      │
+  │     │     Y^-0.5 │
+  │     │     │      │
+  │     @─────@      │
+  │     │     │      │
+┌╴│     │     │      │     ╶┐
+│ │     │     Y^0.5  │      │
+│ │     @─────┼──────@      │
+└╴│     │     │      │     ╶┘
+  │     │     │      │
+  │     │     Y^-0.5 Y^0.5
+  │     │     │      │
+  @─────┼─────@      Rz(0)
+  │     │     │      │
+  │     │     Y^0.5  Y^-0.5
+  │     │     │      │
+┌╴│     │     │      │     ╶┐
+│ │     │     Rz(0)  │      │
+│ │     @─────┼──────@      │
+└╴│     │     │      │     ╶┘
+  │     │     │      │
+  │     Rx(0) Y^-0.5 Y^0.5
+  │     │     │      │
+  @─────┼─────@      Y^-0.5
+  │     │     │      │
+┌╴│     │     │      │     ╶┐
+│ │     │     Y^0.5  │      │
+│ @─────┼─────┼──────@      │
+└╴│     │     │      │     ╶┘
+  │     │     │      │
+  │     │     Rx(0)  Y^0.5
+  │     │     │      │
+  │     │     │      Rz(0)
+  │     │     │      │
+  │     │     │      Y^-0.5
+  │     │     │      │
+  @─────┼─────┼──────@
+  │     │     │      │
+  Rx(0) │     │      Y^0.5
+  │     │     │      │
+  │     │     │      Rx(0)
+  │     │     │      │
+  M─────M─────M──────M('m')
+  │     │     │      │
 
-Example QAOA circuit:
-  0           1           2           3           4           5
-  │           │           │           │           │           │
-  H           H           H           H           H           H
-  │           │           │           │           │           │
-  ZZ──────────ZZ^(-4/13)  │           │           │           │
-┌ │           │           │           │           │           │           ┐
-│ ZZ──────────┼───────────ZZ^(-4/13)  │           │           │           │
-│ │           ZZ──────────┼───────────ZZ^(-4/13)  │           │           │
-└ │           │           │           │           │           │           ┘
-┌ │           │           │           │           │           │           ┐
-│ ZZ──────────┼───────────┼───────────┼───────────ZZ^(-4/13)  │           │
-│ │           ZZ──────────┼───────────┼───────────┼───────────ZZ^(-4/13)  │
-└ │           │           │           │           │           │           ┘
-  Rx(0.151π)  Rx(0.151π)  ZZ──────────┼───────────ZZ^(-4/13)  │
-  │           │           │           │           │           │
-  ZZ──────────ZZ^-0.941   ZZ──────────┼───────────┼───────────ZZ^(-4/13)
-  │           │           │           ZZ──────────ZZ^(-4/13)  │
-┌ │           │           │           │           │           │           ┐
-│ │           │           Rx(0.151π)  ZZ──────────┼───────────ZZ^(-4/13)  │
-│ │           │           │           │           Rx(0.151π)  │           │
-└ │           │           │           │           │           │           ┘
-  ZZ──────────┼───────────ZZ^-0.941   Rx(0.151π)  │           Rx(0.151π)
-┌ │           │           │           │           │           │           ┐
-│ ZZ──────────┼───────────┼───────────┼───────────ZZ^-0.941   │           │
-│ │           ZZ──────────┼───────────ZZ^-0.941   │           │           │
-└ │           │           │           │           │           │           ┘
-  Rx(-0.448π) ZZ──────────┼───────────┼───────────┼───────────ZZ^-0.941
-  │           │           ZZ──────────┼───────────ZZ^-0.941   │
-  │           │           │           │           │           │
-  │           Rx(-0.448π) ZZ──────────┼───────────┼───────────ZZ^-0.941
-  │           │           │           ZZ──────────ZZ^-0.941   │
-┌ │           │           │           │           │           │           ┐
-│ │           │           Rx(-0.448π) ZZ──────────┼───────────ZZ^-0.941   │
-│ │           │           │           │           Rx(-0.448π) │           │
-└ │           │           │           │           │           │           ┘
-  │           │           │           Rx(-0.448π) │           Rx(-0.448π)
-  │           │           │           │           │           │
-  M('m')──────M───────────M───────────M───────────M───────────M
-  │           │           │           │           │           │
-Optimizing objective function ...
-The largest cut value found was 7.
-The largest possible cut has size 7.
-The approximation ratio achieved is 1.0.
+=== EXAMPLE OUTPUT ===
+Brute force max cut: 16
+μ=7.80 max=12
+μ=7.40 max=8
+μ=8.60 max=10
+μ=8.00 max=10
+μ=8.80 max=12
+μ=8.20 max=12
+μ=10.00 max=16
+μ=8.00 max=16
+μ=4.40 max=12
+μ=5.80 max=12
+
+   Return from subroutine COBYLA because the MAXFUN limit has been reached.
+
+   NFVALS =   10   F =-1.000000E+01    MAXCV = 0.000000E+00
+   X = 0.000000E+00   1.000000E+00   0.000000E+00   1.000000E+00   0.000000E+00
+       1.000000E+00   0.000000E+00   0.000000E+00   0.000000E+00   0.000000E+00
 """
-
 import itertools
+from typing import List
 
 import numpy as np
 import networkx
 import scipy.optimize
+from sympy.parsing.sympy_parser import parse_expr
 
 import cirq
 
 
-def main(repetitions=1000, maxiter=50):
-    # Set problem parameters
-    n = 6
-    p = 2
+def brute_force(graph, n):
+    bitstrings = np.array(list(itertools.product(range(2), repeat=n)))
+    mat = networkx.adjacency_matrix(graph, nodelist=sorted(graph.nodes))
+    vecs = (-1) ** bitstrings
+    vals = 0.5 * np.sum(vecs * (mat @ vecs.T).T, axis=-1)
+    vals = 0.5 * (graph.size() - vals)
+    return max(np.round(vals))
 
-    # Generate a random 3-regular graph on n nodes
-    graph = networkx.random_regular_graph(3, n)
 
-    # Make qubits
-    qubits = cirq.LineQubit.range(n)
+def qaoa(booleans: List[str], repetitions: int, maxiter: int, p: int):
+    """Run the QAOA optimization for a list of Boolean expressions.
 
-    # Print an example circuit
-    betas = np.random.uniform(-np.pi, np.pi, size=p)
-    gammas = np.random.uniform(-np.pi, np.pi, size=p)
-    circuit = qaoa_max_cut_circuit(qubits, betas, gammas, graph)
-    print('Example QAOA circuit:')
-    print(circuit.to_text_diagram(transpose=True))
+    Args:
+        booleans: A list of Boolean expressions (we want as many of them to be true as possible).
+        repetitions: The number of times to repeat the measurements.
+        maxiter: The number of iterations of the optimizer.
+        p: The number of times to repeat the Hamiltonian gate.
+    """
+    boolean_exprs = [parse_expr(boolean) for boolean in booleans]
+    param_names = cirq.parameter_names(boolean_exprs)
+    qubits = [cirq.NamedQubit(name) for name in param_names]
 
-    # Create variables to store the largest cut and cut value found
-    largest_cut_found = None
-    largest_cut_value_found = 0
+    def f(x):
+        # Build the circuit.
+        circuit = cirq.Circuit()
+        circuit.append(cirq.H.on_each(*qubits))
 
-    # Initialize simulator
-    simulator = cirq.Simulator()
+        for i in range(p):
+            hamiltonian_gate = cirq.BooleanHamiltonian(
+                {q.name: q for q in qubits}, booleans, 2.0 * x[p + i]
+            )
+            circuit.append(hamiltonian_gate)
+            circuit.append(cirq.rx(2.0 * x[i]).on_each(*qubits))
 
-    # Define objective function (we'll use the negative expected cut value)
+        circuit.append(cirq.measure(*qubits, key='m'))
 
-    def f(x):
-        # Create circuit
-        betas = x[:p]
-        gammas = x[p:]
-        circuit = qaoa_max_cut_circuit(qubits, betas, gammas, graph)
-        # Sample bitstrings from circuit
-        result = simulator.run(circuit, repetitions=repetitions)
+        # Measure
+        result = cirq.Simulator().run(circuit, repetitions=repetitions)
         bitstrings = result.measurements['m']
-        # Process bitstrings
-        nonlocal largest_cut_found
-        nonlocal largest_cut_value_found
-        values = cut_values(bitstrings, graph)
-        max_value_index = np.argmax(values)
-        max_value = values[max_value_index]
-        if max_value > largest_cut_value_found:
-            largest_cut_value_found = max_value
-            largest_cut_found = bitstrings[max_value_index]
-        mean = np.mean(values)
-        return -mean
-
-    # Pick an initial guess
-    x0 = np.random.uniform(-np.pi, np.pi, size=2 * p)
-
-    # Optimize f
-    print('Optimizing objective function ...')
-    scipy.optimize.minimize(f, x0, method='Nelder-Mead', options={'maxiter': maxiter})
 
-    # Compute best possible cut value via brute force search
-    all_bitstrings = np.array(list(itertools.product(range(2), repeat=n)))
-    all_values = cut_values(all_bitstrings, graph)
-    max_cut_value = np.max(all_values)
+        # Evaluate
+        values = []
+        for rep in range(repetitions):
+            subs = {name: val == 1 for name, val in zip(param_names, bitstrings[rep, :])}
+            values.append(
+                sum(1 if boolean_expr.subs(subs) else 0 for boolean_expr in boolean_exprs)
+            )
 
-    # Print the results
-    print(f'The largest cut value found was {largest_cut_value_found}.')
-    print(f'The largest possible cut has size {max_cut_value}.')
-    print(f'The approximation ratio achieved is {largest_cut_value_found / max_cut_value}.')
+        print('μ=%.2f max=%d' % (np.mean(values), max(values)))
 
+        return -np.mean(values)
 
-def rzz(rads):
-    """Returns a gate with the matrix exp(-i Z⊗Z rads)."""
-    return cirq.ZZPowGate(exponent=2 * rads / np.pi, global_shift=-0.5)
+    x0 = np.zeros(2 * p)
+    scipy.optimize.minimize(f, x0, method='COBYLA', options={'maxiter': maxiter, 'disp': True})
 
 
-def qaoa_max_cut_unitary(qubits, betas, gammas, graph):  # Nodes should be integers
-    for beta, gamma in zip(betas, gammas):
-        yield (rzz(-0.5 * gamma).on(qubits[i], qubits[j]) for i, j in graph.edges)
-        yield cirq.rx(2 * beta).on_each(*qubits)
+def main(repetitions=10, maxiter=250, p=5):
+    # Set problem parameters
+    n = 6
 
+    # Generate a random bipartite graph.
+    graph = networkx.complete_multipartite_graph(n, n)
 
-def qaoa_max_cut_circuit(qubits, betas, gammas, graph):  # Nodes should be integers
-    return cirq.Circuit(
-        # Prepare uniform superposition
-        cirq.H.on_each(*qubits),
-        # Apply QAOA unitary
-        qaoa_max_cut_unitary(qubits, betas, gammas, graph),
-        # Measure
-        cirq.measure(*qubits, key='m'),
-    )
+    # Compute best possible cut value via brute force search
+    print('Brute force max cut: %d' % (brute_force(graph, 2 * n)))
 
+    # Build the boolean expressions
+    booleans = [f"x{i} ^ x{j}" for i, j in graph.edges]
 
-def cut_values(bitstrings, graph):
-    mat = networkx.adjacency_matrix(graph, nodelist=sorted(graph.nodes))
-    vecs = (-1) ** bitstrings
-    vals = 0.5 * np.sum(vecs * (mat @ vecs.T).T, axis=-1)
-    vals = 0.5 * (graph.size() - vals)
-    return vals
+    qaoa(booleans, repetitions=repetitions, maxiter=maxiter, p=p)
 
 
 if __name__ == '__main__':