feat(decompositions): Clements decomposition

The Clements decomposition has been implemented in this patch
(truthfully, it has mostly been copied from my thesis work).

`scipy` had to be added as a development dependency for the random
unitary matrix generation for the tests.

`BaseOperator` has been extended with the `is_unitary` method to verify
that the current object is unitary.

Several tests have been written to tests
`piquasso.decompositions.clements` under `tests.test_clements`.

All the fixtures have been placed under `tests.conftest`: this way we
can share the fixtures between all tests in the `tests` module.

References: Budapest-Quantum-Computing-Group#6

piquasso/decompositions/clements.py

@@ -0,0 +1,199 @@
+#
+# Copyright (C) 2020 by TODO - All rights reserved.
+#
+
+"""
+Implementation for the Clements decomposition
+
+References
+----------
+.. [1] William R. Clements, Peter C. Humphreys, Benjamin J. Metcalf,
+    W. Steven Kolthammer, Ian A. Walmsley, "An Optimal Design for
+    Universal Multiport Interferometers", arXiv:1603.08788.
+"""
+
+import numpy as np
+
+from piquasso.operator import BaseOperator
+
+
+class T(BaseOperator):
+    """
+    A definition for the representation
+    of the beamsplitter.
+    """
+
+    def __new__(cls, operation, d):
+        """Produces the beamsplitter matrix.
+
+        The matrix is automatically embedded in a `d` times `d` identity
+        matrix, which can be readily applied during decomposition.
+
+        Args:
+            operation (dict): A dict containing the angle parameters and the
+                modes on which the beamsplitter gate are applied.
+            d (int): The total number of modes.
+
+        Returns:
+            T: The beamsplitter matrix.
+        """
+        theta, phi = operation["params"]
+        i, j = operation["modes"]
+
+        matrix = np.array(
+            [
+                [np.exp(1j * phi) * np.cos(theta), -np.sin(theta)],
+                [np.exp(1j * phi) * np.sin(theta), np.cos(theta)],
+            ]
+        )
+
+        self = np.asarray(np.identity(d, dtype=complex))
+
+        self[i, i] = matrix[0, 0]
+        self[i, j] = matrix[0, 1]
+        self[j, i] = matrix[1, 0]
+        self[j, j] = matrix[1, 1]
+
+        return self.view(BaseOperator)
+
+    @classmethod
+    def transposed(cls, operation, d):
+        """Produces the transposed beamsplitter matrix.
+
+        Args:
+            operation (dict): A dict containing the angle parameters and the
+                modes on which the beamsplitter gate are applied.
+            d (int): The total number of modes.
+
+        Returns:
+            T: The transposed beamsplitter matrix.
+        """
+
+        theta, phi = operation["params"]
+
+        return np.transpose(
+            cls({"params": (theta, -phi), "modes": operation["modes"]}, d=d)
+        )
+
+    @classmethod
+    def i(cls, *args, **kwargs):
+        """Shorthand for :meth:`transposed`.
+
+        The inverse of the matrix equals the transpose in this case.
+
+        Returns:
+            T: The transposed beamsplitter matrix.
+        """
+        return cls.transposed(*args, **kwargs)
+
+
+class Clements:
+    def __init__(self, U, decompose=True):
+        """
+        Args:
+            U (numpy.ndarray): The (square) unitary matrix to be decomposed.
+            decompose (bool): Optional, if `True`, the decomposition is
+                automatically calculated. Defaults to `True`.
+        """
+        self.U = U
+        self.d = U.shape[0]
+        self.inverse_operations = []
+        self.direct_operations = []
+        self.diagonals = None
+
+        if decompose:
+            self.decompose()
+
+    def decompose(self):
+        """
+        Decomposes the specified unitary matrix by application of beamsplitters
+        prescribed by the decomposition.
+        """
+
+        for column in reversed(range(0, self.d - 1)):
+            if column % 2 == 0:
+                self.apply_direct_beamsplitters(column)
+            else:
+                self.apply_inverse_beamsplitters(column)
+
+        self.diagonals = np.diag(self.U)
+
+    def apply_direct_beamsplitters(self, column):
+        """
+        Calculates the direct beamsplitters for a given column `column`, and
+        applies it to `U`.
+
+        Args:
+            column (int): The current column.
+        """
+
+        for j in range(self.d - 1 - column):
+            operation = self.eliminate_lower_offdiagonal(column + j + 1, j)
+            self.direct_operations.append(operation)
+
+            beamsplitter = T(operation, d=self.d)
+
+            self.U = beamsplitter @ self.U
+
+    def apply_inverse_beamsplitters(self, column):
+        """
+        Calculates the inverse beamsplitters for a given column `column`, and
+        applies it to `U`.
+
+        Args:
+            column (int): The current column.
+        """
+
+        for j in reversed(range(self.d - 1 - column)):
+            operation = self.eliminate_upper_offdiagonal(column + j + 1, j)
+            self.inverse_operations.append(operation)
+
+            beamsplitter = T.i(operation, d=self.d)
+            self.U = self.U @ beamsplitter
+
+    def eliminate_lower_offdiagonal(self, i, j):
+        """
+        Calculates the parameters required to eliminate the lower triangular
+        element `i`, `j` of `U` using `T`.
+        """
+        r = -self.U[i, j] / self.U[i - 1, j]
+        theta = np.arctan(np.abs(r))
+        phi = np.angle(r)
+
+        return {
+            "modes": (i - 1, i),
+            "params": (theta, phi),
+        }
+
+    def eliminate_upper_offdiagonal(self, i, j):
+        """
+        Calculates the parameters required to eliminate the upper triangular
+        `i`, `j` of `U` using `T.transposed`.
+        """
+        r = self.U[i, j] / self.U[i, j + 1]
+        theta = np.arctan(np.abs(r))
+        phi = np.angle(r)
+
+        return {
+            "modes": (j, j + 1),
+            "params": (theta, phi),
+        }
+
+    @staticmethod
+    def from_decomposition(decomposition):
+        """
+        Creates the unitary operator from the Clements operations.
+        """
+        U = np.identity(decomposition.d, dtype=complex)
+
+        for operation in decomposition.inverse_operations:
+            beamsplitter = T(operation, d=decomposition.d)
+            U = beamsplitter @ U
+
+        U = np.diag(decomposition.diagonals) @ U
+
+        for operation in reversed(decomposition.direct_operations):
+            beamsplitter = T.i(operation, d=decomposition.d)
+            U = beamsplitter @ U
+
+        return BaseOperator(U)

piquasso/operator.py

@@ -28,3 +28,15 @@ def adjoint(self):
             (BaseOperator): The adjoint of the operator.
         """
         return self.conjugate().transpose().view(self.__class__)
+
+    def is_unitary(self, tol=1e-10):
+        """
+        Args:
+            tol (float, optional): The tolerance for testing the unitarity.
+                Defaults to `1e-10`.
+
+        Returns:
+            bool: `True` if the current object is unitary within the specified
+                tolerance `tol`, else `False`.
+        """
+        return (self @ self.adjoint() - np.identity(self.shape[0]) < tol).all()

pyproject.toml

@@ -14,6 +14,7 @@ pytest = "^5.2"
 strawberryfields = "^0.14.0"
 pennylane = "^0.10.0"
 pennylane-sf = "^0.9.0"
+scipy = "^1.5.2"
 
 [build-system]
 requires = ["poetry>=0.12"]

tests/conftest.py

@@ -0,0 +1,28 @@
+#
+# Copyright (C) 2020 by TODO - All rights reserved.
+#
+
+import pytest
+
+from scipy.stats import unitary_group
+
+from piquasso.operator import BaseOperator
+from piquasso.state import State
+
+
+@pytest.fixture
+def dummy_unitary():
+    def func(d):
+        return BaseOperator(unitary_group.rvs(d))
+
+    return func
+
+
+@pytest.fixture
+def dummy_state():
+    return State.from_state_vector([1, 0])
+
+
+@pytest.fixture(scope="session")
+def tolerance():
+    return 10E-10

tests/test_clements.py

@@ -0,0 +1,120 @@
+#
+# Copyright (C) 2020 by TODO - All rights reserved.
+#
+
+import pytest
+import numpy as np
+
+from piquasso.decompositions.clements import T, Clements
+
+
+def test_T_beamsplitter_is_unitary():
+    theta = np.pi / 3
+    phi = np.pi / 4
+
+    beamsplitter = T({"params": [theta, phi], "modes": [0, 1]}, d=2)
+
+    assert beamsplitter.is_unitary()
+
+
+def test_eliminate_lower_offdiagonal_2_modes(dummy_unitary, tolerance):
+
+    U = dummy_unitary(d=2)
+
+    decomposition = Clements(U, decompose=False)
+
+    operation = decomposition.eliminate_lower_offdiagonal(1, 0)
+
+    beamsplitter = T(operation, 2)
+
+    rotated_U = beamsplitter @ U
+
+    assert np.abs(rotated_U[1, 0]) < tolerance
+
+
+def test_eliminate_lower_offdiagonal_3_modes(dummy_unitary, tolerance):
+
+    U = dummy_unitary(d=3)
+
+    decomposition = Clements(U, decompose=False)
+
+    operation = decomposition.eliminate_lower_offdiagonal(1, 0)
+
+    beamsplitter = T(operation, 3)
+
+    rotated_U = beamsplitter @ U
+
+    assert np.abs(rotated_U[1, 0]) < tolerance
+
+
+def test_eliminate_upper_offdiagonal_2_modes(dummy_unitary, tolerance):
+
+    U = dummy_unitary(d=2)
+
+    decomposition = Clements(U, decompose=False)
+
+    operation = decomposition.eliminate_upper_offdiagonal(1, 0)
+
+    beamsplitter = T.i(operation, 2)
+
+    rotated_U = U @ beamsplitter
+
+    assert np.abs(rotated_U[1, 0]) < tolerance
+
+
+def test_eliminate_upper_offdiagonal_3_modes(dummy_unitary, tolerance):
+
+    U = dummy_unitary(d=3)
+
+    decomposition = Clements(U, decompose=False)
+
+    operation = decomposition.eliminate_upper_offdiagonal(1, 0)
+
+    beamsplitter = T.i(operation, 3)
+
+    rotated_U = U @ beamsplitter
+
+    assert np.abs(rotated_U[1, 0]) < tolerance
+
+
+@pytest.mark.parametrize("n", [2, 3, 4, 5])
+def test_clements_decomposition_on_n_modes(
+            n, dummy_unitary, tolerance
+        ):
+
+    U = dummy_unitary(d=n)
+
+    decomposition = Clements(U)
+
+    diagonalized = decomposition.U
+
+    assert np.abs(diagonalized[0, 1]) < tolerance
+    assert np.abs(diagonalized[1, 0]) < tolerance
+
+    assert (
+        sum(sum(np.abs(diagonalized)))
+        - sum(np.abs(np.diag(diagonalized))) < tolerance
+    ), "Absolute sum of matrix elements should be equal to the "
+    "diagonal elements, if the matrix is properly diagonalized."
+
+
+@pytest.mark.parametrize("n", [2, 3, 4, 5])
+def test_clements_decomposition_and_composition_on_n_modes(
+            n, dummy_unitary, tolerance
+        ):
+
+    U = dummy_unitary(d=n)
+
+    decomposition = Clements(U)
+
+    diagonalized = decomposition.U
+
+    assert (
+        sum(sum(np.abs(diagonalized)))
+        - sum(np.abs(np.diag(diagonalized))) < tolerance
+    ), "Absolute sum of matrix elements should be equal to the "
+    "diagonal elements, if the matrix is properly diagonalized."
+
+    original = Clements.from_decomposition(decomposition)
+
+    assert (U - original < tolerance).all()