From 4f84d4b2a1687afd9376dae629ea95f33defa893 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Zolt=C3=A1n=20Kolarovszki?= Date: Tue, 28 Jul 2020 08:28:14 +0200 Subject: [PATCH] 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: #6 --- piquasso/decompositions/__init__.py | 0 piquasso/decompositions/clements.py | 199 ++++++++++++++++++++++++++++ piquasso/operator.py | 12 ++ poetry.lock | 2 +- pyproject.toml | 1 + tests/conftest.py | 28 ++++ tests/test_clements.py | 120 +++++++++++++++++ tests/test_program.py | 12 -- 8 files changed, 361 insertions(+), 13 deletions(-) create mode 100644 piquasso/decompositions/__init__.py create mode 100644 piquasso/decompositions/clements.py create mode 100644 tests/conftest.py create mode 100644 tests/test_clements.py diff --git a/piquasso/decompositions/__init__.py b/piquasso/decompositions/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/piquasso/decompositions/clements.py b/piquasso/decompositions/clements.py new file mode 100644 index 00000000..acc353c2 --- /dev/null +++ b/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) diff --git a/piquasso/operator.py b/piquasso/operator.py index ecdb379a..1d8472de 100644 --- a/piquasso/operator.py +++ b/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() diff --git a/poetry.lock b/poetry.lock index 5f74eacf..93b0a08b 100644 --- a/poetry.lock +++ b/poetry.lock @@ -538,7 +538,7 @@ docs = ["sphinx", "jaraco.packaging (>=3.2)", "rst.linker (>=1.9)"] testing = ["jaraco.itertools", "func-timeout"] [metadata] -content-hash = "698cf07668f8b1dd2b2399d516518bd5a2785aa34a850028f5cd403a59c6494f" +content-hash = "6fba3709cb624fbee0a438cbe1c0795b898231cad249b93a64ed92863e8a74fc" python-versions = "^3.6" [metadata.files] diff --git a/pyproject.toml b/pyproject.toml index b17da912..a006663d 100644 --- a/pyproject.toml +++ b/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"] diff --git a/tests/conftest.py b/tests/conftest.py new file mode 100644 index 00000000..fa9dbdeb --- /dev/null +++ b/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 diff --git a/tests/test_clements.py b/tests/test_clements.py new file mode 100644 index 00000000..faabb35c --- /dev/null +++ b/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() diff --git a/tests/test_program.py b/tests/test_program.py index 51b84021..ae5c5205 100644 --- a/tests/test_program.py +++ b/tests/test_program.py @@ -2,27 +2,15 @@ # Copyright (C) 2020 by TODO - All rights reserved. # -import pytest import numpy as np from piquasso.program import Program from piquasso.context import Context from piquasso.gates import B from piquasso.mode import Q -from piquasso.state import State from piquasso.backend import FockBackend -@pytest.fixture -def dummy_state(): - return State.from_state_vector([1, 0]) - - -@pytest.fixture(scope="session") -def tolerance(): - return 10E-9 - - def test_current_program_in_program_context(dummy_state): program = Program(state=dummy_state, backend=FockBackend)