add out_flow_constraint function and private function, add unit tests (…

…#1220)

* start qaoa cycles module, add mapping between edges and wires, add tests

* update formatting and imports

* update qaoa documentation, change cycles.py to cycle.py

* add private helper functions for squaring Hamiltonian terms and collecting duplicate terms, add unit tests

* add out_flow_constraint function and private function, add unit tests

* Apply black to tests

* update WIP entry in CHANGELOG.md

* Apply suggestions from code review

Co-authored-by: Tom Bromley <49409390+trbromley@users.noreply.github.com>

* remove _collect_duplicates

* add assertions for operators in unit tests

* fix docstring spacing, remove type hint

* apply black

* remove _collect_duplicates()

* apply black

* Line widths

* typo

* Apply suggestions from code review

Co-authored-by: Tom Bromley <49409390+trbromley@users.noreply.github.com>

* add WIP entry to CHANGELOG.md

* retrigger checks

* apply black

Co-authored-by: trbromley <brotho02@gmail.com>
Co-authored-by: Tom Bromley <49409390+trbromley@users.noreply.github.com>

.github/CHANGELOG.md

@@ -61,8 +61,9 @@
   [(#1209)](https://github.com/PennyLaneAI/pennylane/pull/1209)
   [(#1251)](https://github.com/PennyLaneAI/pennylane/pull/1251)
   [(#1213)](https://github.com/PennyLaneAI/pennylane/pull/1213)
+  [(#1220)](https://github.com/PennyLaneAI/pennylane/pull/1220)
   [(#1214)](https://github.com/PennyLaneAI/pennylane/pull/1214)
-
+  
 * Adds `QubitCarry` and `QubitSum` operations for basic arithmetic.
   [(#1169)](https://github.com/PennyLaneAI/pennylane/pull/1169)
 

pennylane/qaoa/cycle.py

@@ -308,6 +308,52 @@ def _square_hamiltonian_terms(
     return squared_coeffs, squared_ops
 
 
+def out_flow_constraint(graph: nx.DiGraph) -> qml.Hamiltonian:
+    r"""Calculates the `out flow constraint <https://1qbit.com/whitepaper/arbitrage/>`__
+    Hamiltonian for the maximum-weighted cycle problem.
+
+    Given a subset of edges in a directed graph, the out-flow constraint imposes that at most one
+    edge can leave any given node, i.e., for all :math:`i`:
+
+    .. math:: \sum_{j,(i,j)\in E}x_{ij} \leq 1,
+
+    where :math:`E` are the edges of the graph and :math:`x_{ij}` is a binary number that selects
+    whether to include the edge :math:`(i, j)`.
+
+    A set of edges satisfies the out-flow constraint whenever the following Hamiltonian is minimized:
+
+    .. math::
+
+        \sum_{i\in V}\left(d_{i}^{out}(d_{i}^{out} - 2)\mathbb{I}
+        - 2(d_{i}^{out}-1)\sum_{j,(i,j)\in E}\hat{Z}_{ij} +
+        \left( \sum_{j,(i,j)\in E}\hat{Z}_{ij} \right)^{2}\right)
+
+
+    where :math:`V` are the graph vertices, :math:`d_{i}^{\rm out}` is the outdegree of node
+    :math:`i`, and :math:`Z_{ij}` is a qubit Pauli-Z matrix acting
+    upon the qubit specified by the pair :math:`(i, j)`. Mapping from edges to wires can be achieved
+    using :func:`~.edges_to_wires`.
+
+    Args:
+        graph (nx.DiGraph): the directed graph specifying possible edges
+
+    Returns:
+        qml.Hamiltonian: the out flow constraint Hamiltonian
+
+    Raises:
+        ValueError: if the input graph is not directed
+    """
+    if not hasattr(graph, "out_edges"):
+        raise ValueError("Input graph must be directed")
+
+    hamiltonian = qml.Hamiltonian([], [])
+
+    for node in graph.nodes:
+        hamiltonian += _inner_out_flow_constraint_hamiltonian(graph, node)
+
+    return hamiltonian
+
+
 def net_flow_constraint(graph: nx.DiGraph) -> qml.Hamiltonian:
     r"""Calculates the `net flow constraint <https://doi.org/10.1080/0020739X.2010.526248>`__
     Hamiltonian for the maximum-weighted cycle problem.
@@ -332,6 +378,7 @@ def net_flow_constraint(graph: nx.DiGraph) -> qml.Hamiltonian:
     Pauli-Z matrix acting upon the wire specified by the pair :math:`(i, j)`. Mapping from edges to
     wires can be achieved using :func:`~.edges_to_wires`.
 
+
     Args:
         graph (nx.DiGraph): the directed graph specifying possible edges
 
@@ -352,9 +399,55 @@ def net_flow_constraint(graph: nx.DiGraph) -> qml.Hamiltonian:
     return hamiltonian
 
 
+def _inner_out_flow_constraint_hamiltonian(graph: nx.DiGraph, node) -> qml.Hamiltonian:
+    r"""Calculates the inner portion of the Hamiltonian in :func:`out_flow_constraint`.
+    For a given :math:`i`, this function returns:
+
+    .. math::
+
+        d_{i}^{out}(d_{i}^{out} - 2)\mathbb{I}
+        - 2(d_{i}^{out}-1)\sum_{j,(i,j)\in E}\hat{Z}_{ij} +
+        ( \sum_{j,(i,j)\in E}\hat{Z}_{ij}) )^{2}
+
+    Args:
+        graph (nx.DiGraph): the directed graph specifying possible edges
+        node: a fixed node
+
+    Returns:
+        qml.Hamiltonian: The inner part of the out-flow constraint Hamiltonian.
+    """
+    coeffs = []
+    ops = []
+
+    edges_to_qubits = edges_to_wires(graph)
+    out_edges = graph.out_edges(node)
+    d = len(out_edges)
+
+    for edge in out_edges:
+        wire = (edges_to_qubits[edge],)
+        coeffs.append(1)
+        ops.append(qml.PauliZ(wire))
+
+    coeffs, ops = _square_hamiltonian_terms(coeffs, ops)
+
+    for edge in out_edges:
+        wire = (edges_to_qubits[edge],)
+        coeffs.append(-2 * (d - 1))
+        ops.append(qml.PauliZ(wire))
+
+    coeffs.append(d * (d - 2))
+    ops.append(qml.Identity(0))
+
+    H = qml.Hamiltonian(coeffs, ops)
+    H.simplify()
+
+    return H
+
+
 def _inner_net_flow_constraint_hamiltonian(graph: nx.DiGraph, node) -> qml.Hamiltonian:
     r"""Calculates the squared inner portion of the Hamiltonian in :func:`net_flow_constraint`.
 
+
     For a given :math:`i`, this function returns:
 
     .. math::

tests/test_qaoa.py

@@ -31,6 +31,8 @@
     _square_hamiltonian_terms,
     cycle_mixer,
     _partial_cycle_mixer,
+    out_flow_constraint,
+    _inner_out_flow_constraint_hamiltonian,
 )
 from scipy.linalg import expm
 from scipy.sparse import csc_matrix, kron
@@ -813,9 +815,9 @@ def test_cycle_mixer(self):
                 flows[start] += 1
                 flows[end] -= 1
 
-            # A bitstring is valid if the net flow is zero and we aren't the empty set or the set of all
-            # edges. Note that the max out-flow constraint is not imposed, which means we can pass
-            # through nodes more than once
+            # A bitstring is valid if the net flow is zero and we aren't the empty set or the set
+            # of all edges. Note that the max out-flow constraint is not imposed, which means we can
+            # pass through nodes more than once
             if sum(np.abs(flows)) == 0 and 0 < len(edges) < n_wires:
                 valid_bitstrings_indx.append(indx)
             else:
@@ -838,7 +840,8 @@ def test_cycle_mixer(self):
         # Now consider a unitary generated by the Hamiltonian
         h_matrix_e = expm(1j * h_matrix)
 
-        # We expect non-zero transitions among the set of valid bitstrings, and no transitions outside
+        # We expect non-zero transitions among the set of valid bitstrings, and no transitions
+        # outside
         for indx in valid_bitstrings_indx:
             column = h_matrix_e[:, indx]
             destination_indxs = np.argwhere(column != 0).flatten().tolist()
@@ -1064,6 +1067,27 @@ def test_square_hamiltonian_terms(self):
             ]
         )
 
+    def test_inner_out_flow_constraint_hamiltonian(self):
+        """Test if the _inner_out_flow_constraint_hamiltonian function returns the expected result
+        on a manually-calculated example of a 3-node complete digraph relative to the 0 node"""
+
+        g = nx.complete_graph(3).to_directed()
+        h = _inner_out_flow_constraint_hamiltonian(g, 0)
+
+        expected_ops = [
+            qml.Identity(0),
+            qml.PauliZ(0) @ qml.PauliZ(1),
+            qml.PauliZ(0),
+            qml.PauliZ(1),
+        ]
+
+        expected_coeffs = [2, 2, -2, -2]
+
+        assert expected_coeffs == h.coeffs
+        for i, expected_op in enumerate(expected_ops):
+            assert str(h.ops[i]) == str(expected_op)
+        assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
+
     def test_inner_net_flow_constraint_hamiltonian(self):
         """Test if the _inner_net_flow_constraint_hamiltonian function returns the expected result on a manually-calculated
         example of a 3-node complete digraph relative to the 0 node"""
@@ -1086,6 +1110,22 @@ def test_inner_net_flow_constraint_hamiltonian(self):
             assert str(h.ops[i]) == str(expected_op)
         assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
 
+    def test_inner_out_flow_constraint_hamiltonian_non_complete(self):
+        """Test if the _inner_out_flow_constraint_hamiltonian function returns the expected result
+        on a manually-calculated example of a 3-node complete digraph relative to the 0 node, with
+        the (0, 1) edge removed"""
+        g = nx.complete_graph(3).to_directed()
+        g.remove_edge(0, 1)
+        h = _inner_out_flow_constraint_hamiltonian(g, 0)
+
+        expected_ops = [qml.PauliZ(wires=[0])]
+        expected_coeffs = [0]
+
+        assert expected_coeffs == h.coeffs
+        for i, expected_op in enumerate(expected_ops):
+            assert str(h.ops[i]) == str(expected_op)
+        assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
+
     def test_inner_net_flow_constraint_hamiltonian_non_complete(self):
         """Test if the _inner_net_flow_constraint_hamiltonian function returns the expected result on a manually-calculated
         example of a 3-node complete digraph relative to the 0 node, with the (1, 0) edge removed"""
@@ -1109,6 +1149,54 @@ def test_inner_net_flow_constraint_hamiltonian_non_complete(self):
             assert str(h.ops[i]) == str(expected_op)
         assert all([op.wires == exp.wires for op, exp in zip(h.ops, expected_ops)])
 
+    def test_out_flow_constraint(self):
+        """Test the out-flow constraint Hamiltonian is minimised by states that correspond to
+        subgraphs that only ever have 0 or 1 edge leaving each node
+        """
+        g = nx.complete_graph(3).to_directed()
+        h = out_flow_constraint(g)
+        m = wires_to_edges(g)
+        wires = len(g.edges)
+
+        # We use PL to find the energies corresponding to each possible bitstring
+        dev = qml.device("default.qubit", wires=wires)
+
+        def states(basis_state, **kwargs):
+            qml.BasisState(basis_state, wires=range(wires))
+
+        cost = qml.ExpvalCost(states, h, dev, optimize=True)
+
+        # Calculate the set of all bitstrings
+        bitstrings = itertools.product([0, 1], repeat=wires)
+
+        # Calculate the corresponding energies
+        energies_bitstrings = ((cost(bitstring).numpy(), bitstring) for bitstring in bitstrings)
+
+        for energy, bs in energies_bitstrings:
+
+            # convert binary string to wires then wires to edges
+            wires_ = tuple(i for i, s in enumerate(bs) if s != 0)
+            edges = tuple(m[w] for w in wires_)
+
+            # find the number of edges leaving each node
+            num_edges_leaving_node = {node: 0 for node in g.nodes}
+            for e in edges:
+                num_edges_leaving_node[e[0]] += 1
+
+            # check that if the max number of edges is <=1 it corresponds to a state that minimizes
+            # the out_flow_constraint Hamiltonian
+            if max(num_edges_leaving_node.values()) > 1:
+                assert energy > min(energies_bitstrings)[0]
+            elif max(num_edges_leaving_node.values()) <= 1:
+                assert energy == min(energies_bitstrings)[0]
+
+    def test_out_flow_constraint_undirected_raises_error(self):
+        """Test `out_flow_constraint` raises ValueError if input graph is not directed"""
+        g = nx.complete_graph(3)  # undirected graph
+
+        with pytest.raises(ValueError):
+            h = out_flow_constraint(g)
+
     def test_net_flow_constraint(self):
         """Test if the net_flow_constraint Hamiltonian is minimized by states that correspond to a
         collection of edges with zero flow"""