Repository: PennyLaneAI/pennylane
Pull request #1251 (2021-04-22 17:10:07)

Cycle constrained mixer hamiltonian
For the constrained QAOA approach to the maximum weighted cycle problem, the complexity of the constraints is moved to the mixer Hamiltonian.

Add cycle mixer Hamiltonian.

This allows a user to easily implement the constrained approach to this problem. 

This builds on #1207 which added the cost Hamiltonian. The constrained approach would utilise both of these Hamiltonians.

codecov[bot]: # [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) Report
> Merging [#1251](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (7baa594) into [master](https://codecov.io/gh/PennyLaneAI/pennylane/commit/b81a579e47db988db251006ad25516a5c7d2a875?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) (b81a579) will **increase** coverage by `0.00%`.
> The diff coverage is `100.00%`.

[![Impacted file tree graph](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251/graphs/tree.svg?width=650&height=150&src=pr&token=9S3YzlUgGk&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)

```diff
@@           Coverage Diff           @@
##           master    #1251   +/-   ##
=======================================
  Coverage   98.14%   98.14%           
=======================================
  Files         146      146           
  Lines       11035    11079   +44     
=======================================
+ Hits        10830    10874   +44     
  Misses        205      205           
```


| [Impacted Files](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI) | Coverage Δ | |
|---|---|---|
| [pennylane/qaoa/cycle.py](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI#diff-cGVubnlsYW5lL3Fhb2EvY3ljbGUucHk=) | `100.00% <100.00%> (ø)` | |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Last update [b81a579...7baa594](https://codecov.io/gh/PennyLaneAI/pennylane/pull/1251?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PennyLaneAI).

github-actions[bot]: Hello. You may have forgotten to update the changelog!
Please edit [.github/CHANGELOG.md](/PennyLaneAI/pennylane/blob/master/.github/CHANGELOG.md) with:
* A one-to-two sentence description of the change. You may include a small working example for new features.
* A link back to this PR.
* Your name (or GitHub username) in the contributors section.

Patch:
diff --git a/.github/CHANGELOG.md b/.github/CHANGELOG.md
index b15ac5e35..e90dabadc 100644
--- a/.github/CHANGELOG.md
+++ b/.github/CHANGELOG.md
@@ -59,6 +59,7 @@
 * TBD: Add `qml.qaoa.cycle` (when structure is finalised)
   [(#1207)](https://github.com/PennyLaneAI/pennylane/pull/1207)
   [(#1209)](https://github.com/PennyLaneAI/pennylane/pull/1209)
+  [(#1251)](https://github.com/PennyLaneAI/pennylane/pull/1251)
 
 * Adds `QubitCarry` and `QubitSum` operations for basic arithmetic.
   [(#1169)](https://github.com/PennyLaneAI/pennylane/pull/1169)
diff --git a/pennylane/qaoa/cycle.py b/pennylane/qaoa/cycle.py
index 2b347b10a..2d783f327 100644
--- a/pennylane/qaoa/cycle.py
+++ b/pennylane/qaoa/cycle.py
@@ -79,6 +79,114 @@ def wires_to_edges(graph: nx.Graph) -> Dict[int, Tuple]:
     return {i: edge for i, edge in enumerate(graph.edges)}
 
 
+def cycle_mixer(graph: nx.DiGraph) -> qml.Hamiltonian:
+    r"""Calculates the cycle-mixer Hamiltonian.
+
+    Following methods outlined `here <https://arxiv.org/pdf/1709.03489.pdf>`__, the
+    cycle-mixer Hamiltonian preserves the set of valid cycles:
+
+    .. math::
+        \frac{1}{4}\sum_{(i, j)\in E}
+        \left(\sum_{k \in V, k\neq i, k\neq j, (i, k) \in E, (k, j) \in E}
+        \left[X_{ij}X_{ik}X_{kj} +Y_{ij}Y_{ik}X_{kj} + Y_{ij}X_{ik}Y_{kj} - X_{ij}Y_{ik}Y_{kj}\right]
+        \right)
+
+    where :math:`E` are the edges of the directed graph. A valid cycle is defined as a subset of
+    edges in :math:`E` such that all of the graph's nodes :math:`V` have zero net flow (see the
+    :func:`~.net_flow_constraint` function).
+
+    **Example**
+
+    >>> import networkx as nx
+    >>> g = nx.complete_graph(3).to_directed()
+    >>> h_m = cycle_mixer(g)
+    >>> print(h_m)
+      (-0.25) [X0 Y1 Y5]
+    + (-0.25) [X1 Y0 Y3]
+    + (-0.25) [X2 Y3 Y4]
+    + (-0.25) [X3 Y2 Y1]
+    + (-0.25) [X4 Y5 Y2]
+    + (-0.25) [X5 Y4 Y0]
+    + (0.25) [X0 X1 X5]
+    + (0.25) [Y0 Y1 X5]
+    + (0.25) [Y0 X1 Y5]
+    + (0.25) [X1 X0 X3]
+    + (0.25) [Y1 Y0 X3]
+    + (0.25) [Y1 X0 Y3]
+    + (0.25) [X2 X3 X4]
+    + (0.25) [Y2 Y3 X4]
+    + (0.25) [Y2 X3 Y4]
+    + (0.25) [X3 X2 X1]
+    + (0.25) [Y3 Y2 X1]
+    + (0.25) [Y3 X2 Y1]
+    + (0.25) [X4 X5 X2]
+    + (0.25) [Y4 Y5 X2]
+    + (0.25) [Y4 X5 Y2]
+    + (0.25) [X5 X4 X0]
+    + (0.25) [Y5 Y4 X0]
+    + (0.25) [Y5 X4 Y0]
+
+    Args:
+        graph (nx.DiGraph): the directed graph specifying possible edges
+
+    Returns:
+        qml.Hamiltonian: the cycle-mixer Hamiltonian
+    """
+    hamiltonian = qml.Hamiltonian([], [])
+
+    for edge in graph.edges:
+        hamiltonian += _partial_cycle_mixer(graph, edge)
+
+    return hamiltonian
+
+
+def _partial_cycle_mixer(graph: nx.DiGraph, edge: Tuple) -> qml.Hamiltonian:
+    r"""Calculates the partial cycle-mixer Hamiltonian for a specific edge.
+
+    For an edge :math:`(i, j)`, this function returns:
+
+    .. math::
+
+        \sum_{k \in V, k\neq i, k\neq j, (i, k) \in E, (k, j) \in E}\left[
+        X_{ij}X_{ik}X_{kj} + Y_{ij}Y_{ik}X_{kj} + Y_{ij}X_{ik}Y_{kj} - X_{ij}Y_{ik}Y_{kj}\right]
+
+    Args:
+        graph (nx.DiGraph): the directed graph specifying possible edges
+        edge (tuple): a fixed edge
+
+    Returns:
+        qml.Hamiltonian: the partial cycle-mixer Hamiltonian
+    """
+    coeffs = []
+    ops = []
+
+    edges_to_qubits = edges_to_wires(graph)
+
+    for node in graph.nodes:
+        out_edge = (edge[0], node)
+        in_edge = (node, edge[1])
+        if node not in edge and out_edge in graph.edges and in_edge in graph.edges:
+            wire = edges_to_qubits[edge]
+            out_wire = edges_to_qubits[out_edge]
+            in_wire = edges_to_qubits[in_edge]
+
+            t = qml.PauliX(wires=wire) @ qml.PauliX(wires=out_wire) @ qml.PauliX(wires=in_wire)
+            ops.append(t)
+
+            t = qml.PauliY(wires=wire) @ qml.PauliY(wires=out_wire) @ qml.PauliX(wires=in_wire)
+            ops.append(t)
+
+            t = qml.PauliY(wires=wire) @ qml.PauliX(wires=out_wire) @ qml.PauliY(wires=in_wire)
+            ops.append(t)
+
+            t = qml.PauliX(wires=wire) @ qml.PauliY(wires=out_wire) @ qml.PauliY(wires=in_wire)
+            ops.append(t)
+
+            coeffs.extend([0.25, 0.25, 0.25, -0.25])
+
+    return qml.Hamiltonian(coeffs, ops)
+
+
 def loss_hamiltonian(graph: nx.Graph) -> qml.Hamiltonian:
     r"""Calculates the loss Hamiltonian for the maximum-weighted cycle problem.
 


Test patch:
diff --git a/tests/test_qaoa.py b/tests/test_qaoa.py
index de170044a..397d8804a 100644
--- a/tests/test_qaoa.py
+++ b/tests/test_qaoa.py
@@ -15,13 +15,23 @@
 Unit tests for the :mod:`pennylane.qaoa` submodule.
 """
 import pytest
+import itertools
 import numpy as np
 import networkx as nx
 import pennylane as qml
 from pennylane import qaoa
 from networkx import Graph
 from pennylane.wires import Wires
-from pennylane.qaoa.cycle import edges_to_wires, wires_to_edges, loss_hamiltonian
+from pennylane.qaoa.cycle import (
+    edges_to_wires,
+    wires_to_edges,
+    loss_hamiltonian,
+    cycle_mixer,
+    _partial_cycle_mixer,
+)
+from scipy.linalg import expm
+from scipy.sparse import csc_matrix, kron
+
 
 #####################################################
 
@@ -41,6 +51,40 @@ def decompose_hamiltonian(hamiltonian):
     return [coeffs, ops, wires]
 
 
+def matrix(hamiltonian: qml.Hamiltonian, n_wires: int) -> csc_matrix:
+    r"""Calculates the matrix representation of an input Hamiltonian in the standard basis.
+
+    Args:
+        hamiltonian (qml.Hamiltonian): the input Hamiltonian
+        n_wires (int): the total number of wires
+
+    Returns:
+        csc_matrix: a sparse matrix representation
+    """
+    ops_matrices = []
+
+    for op in hamiltonian.ops:
+        op_wires = np.array(op.wires.tolist())
+        op_list = op.non_identity_obs if isinstance(op, qml.operation.Tensor) else [op]
+        op_matrices = []
+
+        for wire in range(n_wires):
+            loc = np.argwhere(op_wires == wire).flatten()
+            mat = np.eye(2) if len(loc) == 0 else op_list[loc[0]].matrix
+            mat = csc_matrix(mat)
+            op_matrices.append(mat)
+
+        op_matrix = op_matrices.pop(0)
+
+        for mat in op_matrices:
+            op_matrix = kron(op_matrix, mat)
+
+        ops_matrices.append(op_matrix)
+
+    mat = sum(coeff * op_mat for coeff, op_mat in zip(hamiltonian.coeffs, ops_matrices))
+    return csc_matrix(mat)
+
+
 class TestMixerHamiltonians:
     """Tests that the mixer Hamiltonians are being generated correctly"""
 
@@ -696,6 +740,142 @@ class TestCycles:
 
         assert r == {0: (0, 1), 1: (0, 2), 2: (0, 3), 3: (1, 2), 4: (1, 3), 5: (2, 3), 6: (3, 4)}
 
+    def test_partial_cycle_mixer_complete(self):
+        """Test if the _partial_cycle_mixer function returns the expected Hamiltonian for a fixed
+        example"""
+        g = nx.complete_graph(4).to_directed()
+        edge = (0, 1)
+
+        h = _partial_cycle_mixer(g, edge)
+
+        ops_expected = [
+            qml.PauliX(0) @ qml.PauliX(1) @ qml.PauliX(7),
+            qml.PauliY(0) @ qml.PauliY(1) @ qml.PauliX(7),
+            qml.PauliY(0) @ qml.PauliX(1) @ qml.PauliY(7),
+            qml.PauliX(0) @ qml.PauliY(1) @ qml.PauliY(7),
+            qml.PauliX(0) @ qml.PauliX(2) @ qml.PauliX(10),
+            qml.PauliY(0) @ qml.PauliY(2) @ qml.PauliX(10),
+            qml.PauliY(0) @ qml.PauliX(2) @ qml.PauliY(10),
+            qml.PauliX(0) @ qml.PauliY(2) @ qml.PauliY(10),
+        ]
+        coeffs_expected = [0.25, 0.25, 0.25, -0.25, 0.25, 0.25, 0.25, -0.25]
+
+        assert h.coeffs == coeffs_expected
+        assert all(op.wires == op_e.wires for op, op_e in zip(h.ops, ops_expected))
+        assert all(op.name == op_e.name for op, op_e in zip(h.ops, ops_expected))
+
+    def test_partial_cycle_mixer_incomplete(self):
+        """Test if the _partial_cycle_mixer function returns the expected Hamiltonian for a fixed
+        example"""
+        g = nx.complete_graph(4).to_directed()
+        g.remove_edge(2, 1)  # remove an egde to make graph incomplete
+        edge = (0, 1)
+
+        h = _partial_cycle_mixer(g, edge)
+
+        ops_expected = [
+            qml.PauliX(0) @ qml.PauliX(2) @ qml.PauliX(9),
+            qml.PauliY(0) @ qml.PauliY(2) @ qml.PauliX(9),
+            qml.PauliY(0) @ qml.PauliX(2) @ qml.PauliY(9),
+            qml.PauliX(0) @ qml.PauliY(2) @ qml.PauliY(9),
+        ]
+        coeffs_expected = [0.25, 0.25, 0.25, -0.25]
+
+        assert h.coeffs == coeffs_expected
+        assert all(op.wires == op_e.wires for op, op_e in zip(h.ops, ops_expected))
+        assert all(op.name == op_e.name for op, op_e in zip(h.ops, ops_expected))
+
+    def test_cycle_mixer(self):
+        """Test if the cycle_mixer Hamiltonian maps valid cycles to valid cycles"""
+        n_nodes = 3
+        g = nx.complete_graph(n_nodes).to_directed()
+        m = wires_to_edges(g)
+        n_wires = len(g.edges)
+
+        # Find Hamiltonian and its matrix representation
+        h = cycle_mixer(g)
+        h_matrix = np.real_if_close(matrix(h, n_wires).toarray())
+
+        # Decide which bitstrings are valid and which are invalid
+        valid_bitstrings_indx = []
+        invalid_bitstrings_indx = []
+
+        for indx, bitstring in enumerate(itertools.product([0, 1], repeat=n_wires)):
+            wires = [i for i, bit in enumerate(bitstring) if bit == 1]
+            edges = [m[wire] for wire in wires]
+
+            flows = [0 for i in range(n_nodes)]
+
+            for start, end in edges:
+                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
+            if sum(np.abs(flows)) == 0 and 0 < len(edges) < n_wires:
+                valid_bitstrings_indx.append(indx)
+            else:
+                invalid_bitstrings_indx.append(indx)
+
+        # Check that valid bitstrings map to a subset of the valid bitstrings
+        for indx in valid_bitstrings_indx:
+            column = h_matrix[:, indx]
+            destination_indxs = set(np.argwhere(column != 0).flatten())
+
+            assert destination_indxs.issubset(valid_bitstrings_indx)
+
+        # Check that invalid bitstrings map to a subset of the invalid bitstrings
+        for indx in invalid_bitstrings_indx:
+            column = h_matrix[:, indx]
+            destination_indxs = set(np.argwhere(column != 0).flatten())
+
+            assert destination_indxs.issubset(invalid_bitstrings_indx)
+
+        # 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
+        for indx in valid_bitstrings_indx:
+            column = h_matrix_e[:, indx]
+            destination_indxs = np.argwhere(column != 0).flatten().tolist()
+            assert destination_indxs == valid_bitstrings_indx
+
+        # Check that invalid bitstrings transition within the set of invalid bitstrings
+        for indx in invalid_bitstrings_indx:
+            column = h_matrix_e[:, indx]
+            destination_indxs = set(np.argwhere(column != 0).flatten().tolist())
+            assert destination_indxs.issubset(invalid_bitstrings_indx)
+
+    def test_matrix(self):
+        """Test that the matrix function works as expected on a fixed example"""
+        g = nx.lollipop_graph(3, 1)
+        h = qml.qaoa.bit_flip_mixer(g, 0)
+
+        mat = matrix(h, 4)
+        mat_expected = np.array(
+            [
+                [0, 1, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
+                [1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+                [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0],
+                [0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+                [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
+            ]
+        )
+
+        assert np.allclose(mat.toarray(), mat_expected)
+
     def test_edges_to_wires_directed(self):
         """Test that edges_to_wires returns the correct mapping on a directed graph"""
         g = nx.lollipop_graph(4, 1).to_directed()