Bell decomp (#584)

* first commit for Bell decomposition

* reck notebook

* Reck seems to be working :D

* clements nearly working

* fix odd external phase index error

* restore the matrix

* put codes in decomposition_test.py and put the demo in the notebook Decomposition_Demo

* sympy example for the beam splitter

* fixed my name spelling XD

* recompose now works in order gates are applied

* added fig4 step of algorithm

* use arctan2 to remove numpy overflow warnings on matrices with 0s in them

* add some notations

* visualize the circuit oK. NEXT STEP: write down the operator sMZI + show the Unitary of the program

* new modifications on the compiler and ops

* reconstruct the matrix is okkkk

* added new compact decompositions

* fix bugs to make decomposition run with engine

* remove notebook from pull request

* fix whitespace problem in code

* check with CodeFactor

* test random unitary on a squeezed state

* test decomposition then recompose

* added recompose files, fix triangular decomposition for matrix with 0s

* lots of tests got modified by black

* disable pylint warning

* try to undo the black file changes

* Update doc/requirements.txt

Co-authored-by: Nicolas Quesada <nicolas@xanadu.ai>

* Update doc/requirements.txt

Co-authored-by: Nicolas Quesada <nicolas@xanadu.ai>

* undo black on test file

* add test for sMZgate on gaussian unitary compiler

* adding back some new tests which i forgot to add when i cleaned up before

* update changelog

* add PR link

* bracket

* run black on decompositions.py

* docstring correction

* added some more context to changelog and added an example

* updated changelog example to remove dependency on recompose function

* move recompose functions to test file

* Update strawberryfields/ops.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

* move command code for interferometer to an independent private function

* Update strawberryfields/decompositions.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

* change Exception to ValueError

* fixed cmd functions by adding reg to args

* tidy some doc strings, delete recompose function which i thought i deleted earlier

* tweaked docstrings, changed decomposition failure error from exception to assert

* fixed assert params

* docstring math fix

* Update .github/CHANGELOG.md

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update .github/CHANGELOG.md

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update .github/CHANGELOG.md

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update strawberryfields/decompositions.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update strawberryfields/decompositions.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update strawberryfields/decompositions.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update strawberryfields/decompositions.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update tests/frontend/test_decompositions.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

* fix the codecov/pylint warning for ops.py

* Update strawberryfields/decompositions.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

* Update strawberryfields/decompositions.py

Co-authored-by: Josh Izaac <josh146@gmail.com>

Co-authored-by: Yuan Yao <yuanyao@YuandeMacBook-Pro.local>
Co-authored-by: Nicolas Quesada <nicolas@xanadu.ai>
Co-authored-by: Josh Izaac <josh146@gmail.com>
Co-authored-by: Yuan <sylvieyao502@gmail.com>

.github/CHANGELOG.md

@@ -2,6 +2,43 @@
 
 <h3>New features since last release</h3>
 
+* Compact decompositions as described in <https://arxiv.org/abs/2104.07561>,
+ (``rectangular_compact`` and ``triangular_compact``) are now available in
+ the ``sf.decompositions`` module, and as options in the ``Interferometer`` operation.
+ [(#584)](https://github.com/XanaduAI/strawberryfields/pull/584)
+
+ This decomposition allows for lower depth photonic circuits in physical devices by applying two 
+ independent phase shifts in parallel inside each Mach-Zehnder interferometer.
+ ``rectangular_compact`` reduces the layers of phase shifters from 2N+1 to N+2 
+ for an N mode interferometer when compared to e.g. ``rectangular_MZ``.
+
+ Example:
+
+ ```python 
+  import numpy as np 
+  from strawberryfields import Program
+  from strawberryfields.ops import Interferometer
+  from scipy.stats import unitary_group
+
+  M = 10
+
+  # generate a 10x10 Haar random unitary
+  U = unitary_group.rvs(M)
+
+  prog = Program(M)
+
+  with prog.context as q:
+      Interferometer(U, mesh='rectangular_compact') | q
+
+  # check that applied unitary is correct 
+  compiled_circuit = prog.compile(compiler="gaussian_unitary")
+  commands = compiled_circuit.circuit
+  S = commands[0].op.p[0] # symplectic transformation
+  Uout = S[:M,:M] + 1j * S[M:,:M] # unitary transformation
+
+  print(np.allclose(U, Uout))
+ ```
+
 <h3>Improvements</h3>
 
 * Cleanup `backends/tfbackend/ops.py` to reduce line count, clarify function
@@ -29,7 +66,7 @@
 
 This release contains contributions from (in alphabetical order):
 
-Aaron Robertson, Jeremy Swinarton, Antal Száva.
+Jake Bulmer, Aaron Robertson, Jeremy Swinarton, Antal Száva, Yuan Yao.
 
 # Release 0.18.0 (current release)
 

strawberryfields/compilers/fock.py

@@ -64,6 +64,7 @@ class Fock(Compiler):
         "CXgate": {},
         "CZgate": {},
         "MZgate": {},
+        "sMZgate": {},
         "Xgate": {},
         "Zgate": {},
         "Fouriergate": {},

strawberryfields/compilers/gaussian.py

@@ -59,6 +59,7 @@ class Gaussian(Compiler):
         "CXgate": {},
         "CZgate": {},
         "MZgate": {},
+        "sMZgate": {},
         "Xgate": {},
         "Zgate": {},
         "Fouriergate": {},

strawberryfields/compilers/gaussian_unitary.py

@@ -83,6 +83,7 @@ class GaussianUnitary(Compiler):
         "Rgate",
         # multi mode gates
         "MZgate",
+        "sMZgate",
         "BSgate",
         "S2gate",
         "Interferometer",  # Note that interferometer is accepted as a primitive
@@ -175,6 +176,17 @@ def compile(self, seq, registers):
                         [dict_indices[modes[0]], dict_indices[modes[1]]],
                         nmodes,
                     )
+                elif name == "sMZgate":
+                    exp_sigma = np.exp(1j * (params[0] + params[1]) / 2)
+                    delta = (params[0] - params[1]) / 2
+                    U = exp_sigma * np.array(
+                        [[np.sin(delta), np.cos(delta)], [np.cos(delta), -np.sin(delta)]]
+                    )
+                    S = expand(
+                        interferometer(U),
+                        [dict_indices[modes[0]], dict_indices[modes[1]]],
+                        nmodes,
+                    )
                 Snet = S @ Snet
                 rnet = S @ rnet
 

strawberryfields/decompositions.py

@@ -17,6 +17,7 @@
 """
 
 from itertools import groupby
+from collections import defaultdict
 
 import numpy as np
 from scipy.linalg import block_diag, sqrtm, polar, schur
@@ -632,6 +633,276 @@ def triangular(V, tol=1e-11):
     return list(reversed(tlist)), np.diag(localV), None
 
 
+def M(n, sigma, delta, m):
+    r"""The symmetric Mach Zehnder interferometer matrix. (Eq 1 of the paper (arXiv:2104.0756).)
+
+    Args:
+        n (int): the starting mode of sMZI
+        sigma (complex): parameter of the sMZI :math:`\frac{(\theta_1+\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        delta (complex): parameter of the sMZI :math:`\frac{(\theta_1-\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        m (int): the length of the unitary matrix to be decomposed
+
+    Returns:
+        array[complex,complex]: the sMZI matrix between n-th and (n+1)-th mode
+    """
+    mat = np.identity(m, dtype=np.complex128)
+    mat[n, n] = np.exp(1j * sigma) * np.sin(delta)
+    mat[n, n + 1] = np.exp(1j * sigma) * np.cos(delta)
+    mat[n + 1, n] = np.exp(1j * sigma) * np.cos(delta)
+    mat[n + 1, n + 1] = -np.exp(1j * sigma) * np.sin(delta)
+    return mat
+
+
+def P(j, phi, m):
+    r"""The phase shifter matrix. (Eq 2 of the paper (arXiv:2104.0756).)
+
+    Args:
+        j (int): the starting mode of phase-shifter
+        phi (complex): parameter of the phase-shifter
+        m (int): the length of the unitary matrix to be decomposed
+
+    Returns:
+        array[complex,complex]: the phase-shifter matrix on the j-th mode
+    """
+    mat = np.identity(m, dtype=np.complex128)
+    mat[j, j] = np.exp(1j * phi)
+    return mat
+
+
+def triangular_compact(U, rtol=1e-12, atol=1e-12):
+    r"""Triangular decomposition of a unitary matrix with sMZIs and phase-shifters, as given in FIG. 2 and "The Reck Scheme" section of (arXiv:2104.0756).
+
+    Args:
+        U (array): unitary matrix
+
+    Returns:
+        dict: A dictionary containing the following items:
+
+        * ``m``: the length of the matrix
+        * ``phi_ins``: parameter of the phase-shifter at the beginning of the mode
+        * ``sigmas``: parameter of the sMZI :math:`\frac{(\theta_1+\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        * ``deltas``: parameter of the sMZI :math:`\frac{(\theta_1-\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        * ``zetas``: parameter of the phase-shifter at the end of the mode
+    """
+
+    if not U.shape[0] == U.shape[1]:
+        raise ValueError("Matrix is not square")
+
+    if not np.allclose(U @ U.conj().T, np.eye(U.shape[0]), rtol=rtol, atol=atol):
+        raise ValueError("The input matrix is not unitary")
+
+    V = U.conj()
+    m = U.shape[0]
+
+    phases = dict()
+    phases["m"] = m
+    phases["phi_ins"] = dict()  # mode : phi
+    phases["deltas"] = dict()  # (mode, layer) : delta
+    phases["sigmas"] = dict()  # (mode, layer) : sigma
+    phases["zetas"] = dict()  # mode : zeta
+
+    for j in range(m - 1):
+        x = m - 1
+        y = j
+        phi_j = -np.angle(V[x, y + 1]) + np.angle(V[x, y])
+        Pj = P(j + 1, phi_j, m)
+        phases["phi_ins"][j] = phi_j
+        V = V @ Pj
+        for k in range(j + 1):
+            n = j - k
+            if V[x, y] == 0:
+                delta = 0.5 * np.pi
+            else:
+                delta = np.arctan2(-abs(V[x, y + 1]), abs(V[x, y]))
+            V_temp = V @ M(n, 0, delta, m)
+            sigma = np.angle(V_temp[x - 1, y - 1]) - np.angle(V_temp[x - 1, y])
+            phases["deltas"][n, k] = delta
+            phases["sigmas"][n, k] = sigma
+            V = V @ M(n, sigma, delta, m)
+            x -= 1
+            y -= 1
+
+    # these next two lines are just to remove a global phase
+    zeta = -np.angle(V[0, 0])
+    phases["zetas"][0] = zeta
+    V = V @ P(0, zeta, m)
+
+    for j in range(1, m):
+        zeta = np.angle(V[0, 0]) - np.angle(V[j, j])
+        phases["zetas"][j] = zeta
+        V = V @ P(j, zeta, m)
+
+    assert np.allclose(V, np.eye(m), rtol=rtol, atol=atol), "decomposition failed"
+
+    return phases
+
+
+def _rectangular_compact_init(
+    U, rtol=1e-12, atol=1e-12
+):  # pylint: disable=too-many-statements, too-many-branches
+    r"""Rectangular decomposition of a unitary with sMZIs and phase-shifters, as given in FIG. 3 and "The Clements Scheme" section of (arXiv:2104.0756).
+
+    Args:
+        U (array): unitary matrix
+
+    Returns:
+        dict: A dictionary containing the following items:
+
+        * ``m``: the length of the matrix
+        * ``phi_ins``: parameter of the phase-shifter at the beginning of the mode
+        * ``sigmas``: parameter of the sMZI :math:`\frac{(\theta_1+\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        * ``deltas``: parameter of the sMZI :math:`\frac{(\theta_1-\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        * ``zetas``: parameter of the phase-shifter at the middle of the mode
+        * ``phi_outs``: parameter of the phase-shifter at the end of the mode
+
+    """
+
+    V = U.conj()
+    m = U.shape[0]
+
+    phases = dict()
+    phases["m"] = m
+    phases["phi_ins"] = dict()  # mode : phi
+    phases["deltas"] = dict()  # (mode, layer) : delta
+    phases["sigmas"] = dict()  # (mode, layer) : sigma
+    phases["zetas"] = dict()  # mode : zeta
+    phases["phi_outs"] = dict()  # mode : phi
+
+    for j in range(m - 1):
+        if j % 2 == 0:
+            x = m - 1
+            y = j
+            phi_j = np.angle(V[x, y + 1]) - np.angle(V[x, y])  # reversed order from paper
+            V = V @ P(j, phi_j, m)
+            phases["phi_ins"][j] = phi_j
+            for k in range(j + 1):
+                if V[x, y] == 0:
+                    delta = 0.5 * np.pi
+                else:
+                    delta = np.arctan2(-abs(V[x, y + 1]), abs(V[x, y]))
+                n = j - k
+                V_temp = V @ M(n, 0, delta, m)
+                sigma = np.angle(V_temp[x - 1, y - 1]) - np.angle(V_temp[x - 1, y])
+                V = V @ M(n, sigma, delta, m)
+                phases["deltas"][n, k] = delta
+                phases["sigmas"][n, k] = sigma
+                x -= 1
+                y -= 1
+        else:
+            x = m - j - 1
+            y = 0
+            phi_j = np.angle(V[x - 1, y]) - np.angle(V[x, y])
+            V = P(x, phi_j, m) @ V
+            phases["phi_outs"][x] = phi_j
+            for k in range(j + 1):
+                if V[x, y] == 0.0:
+                    delta = 0.5 * np.pi
+                else:
+                    delta = np.arctan2(abs(V[x - 1, y]), abs(V[x, y]))
+                V_temp = M(x - 1, 0, delta, m) @ V
+                n = m + k - j - 2
+                if j != k:
+                    sigma = np.angle(V_temp[x + 1, y + 1]) - np.angle(V_temp[x, y + 1])
+                else:
+                    sigma = 0
+                phases["deltas"][n, m - k - 1] = delta
+                phases["sigmas"][n, m - k - 1] = sigma
+                V = M(n, sigma, delta, m) @ V
+                x += 1
+                y += 1
+
+    # these next two lines are just to remove a global phase
+    zeta = -np.angle(V[0, 0])
+    V = V @ P(0, zeta, m)
+    phases["zetas"][0] = zeta
+
+    for j in range(1, m):
+        zeta = np.angle(V[0, 0]) - np.angle(V[j, j])
+        V = V @ P(j, zeta, m)
+        phases["zetas"][j] = zeta
+
+    assert np.allclose(V, np.eye(m), rtol=rtol, atol=atol), "decomposition failed"
+
+    return phases
+
+
+def _absorb_zeta(phases):
+    r"""Adjust rectangular decomposition to relocate residual phase-shifters of interferometer to edge-shifters, as given in FIG. 4 and "Relocating residual phase-shifts" section of (arXiv:2104.0756).
+
+    Args:
+        phases (dict): output of _rectangular_compact_init
+
+    Returns:
+        dict: A dictionary containing the following items:
+
+        * ``m``: the length of the matrix
+        * ``phi_ins``: parameter of the phase-shifter at the beginning of the mode
+        * ``sigmas``: parameter of the sMZI :math:`\frac{(\theta_1+\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        * ``deltas``: parameter of the sMZI :math:`\frac{(\theta_1-\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        * ``phi_edges``: parameters of the edge phase shifters
+        * ``phi_outs``: parameter of the phase-shifter at the end of the mode
+
+    """
+    m = phases["m"]
+    new_phases = phases.copy()
+    del new_phases["zetas"]
+    new_phases["phi_edges"] = defaultdict(float)  # (mode, layer) : phi
+
+    if m % 2 == 0:
+        new_phases["phi_outs"][0] = phases["zetas"][0]
+        for j in range(1, m):
+            zeta = phases["zetas"][j]
+            layer = m - j
+            for mode in range(j, m - 1, 2):
+                new_phases["sigmas"][mode, layer] += zeta
+            for mode in range(j + 1, m - 1, 2):
+                new_phases["sigmas"][mode, layer - 1] -= zeta
+            if layer % 2 == 1:
+                new_phases["phi_edges"][m - 1, layer] += zeta
+            else:
+                new_phases["phi_edges"][m - 1, layer - 1] -= zeta
+    else:
+        for j in range(m):
+            zeta = phases["zetas"][j]
+            layer = m - j - 1
+            for mode in range(j, m - 1, 2):
+                new_phases["sigmas"][mode, layer] += zeta
+            for mode in range(j + 1, m - 1, 2):
+                new_phases["sigmas"][mode, layer - 1] -= zeta
+            if layer % 2 == 0:
+                new_phases["phi_edges"][m - 1, layer] += zeta
+            else:
+                new_phases["phi_edges"][m - 1, layer - 1] -= zeta
+    return new_phases
+
+
+def rectangular_compact(U, rtol=1e-12, atol=1e-12):
+    r"""Rectangular decomposition of a unitary with sMZIs and phase-shifters, as given in FIG. 3+4 and "The Clements Scheme" section of (arXiv:2104.0756).
+
+    Args:
+        U (array): unitary matrix
+
+    Returns:
+        dict: A dictionary containing the following items:
+
+        * ``m``: the length of the matrix
+        * ``phi_ins``: parameter of the phase-shifter at the beginning of the mode
+        * ``sigmas``: parameter of the sMZI :math:`\frac{(\theta_1+\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        * ``deltas``: parameter of the sMZI :math:`\frac{(\theta_1-\theta_2)}{2}`, where :math:`\theta_{1,2}` are the values of the two internal phase-shifts of sMZI
+        * ``phi_edges``: parameters of the edge phase shifters
+        * ``phi_outs``: parameter of the phase-shifter at the end of the mode
+    """
+
+    if not U.shape[0] == U.shape[1]:
+        raise ValueError("Matrix is not square")
+
+    if not np.allclose(U @ U.conj().T, np.eye(U.shape[0]), rtol=rtol, atol=atol):
+        raise ValueError("The input matrix is not unitary")
+
+    phases_temp = _rectangular_compact_init(U, rtol=rtol, atol=atol)
+    return _absorb_zeta(phases_temp)
+
+
 def williamson(V, tol=1e-11):
     r"""Williamson decomposition of positive-definite (real) symmetric matrix.