Removes decomps already in thewalrus

.github/CHANGELOG.md

@@ -10,6 +10,8 @@
   instead accessible via the `sf.program_utils.validate_gate_parameters` function.
   [(#720)](https://github.com/XanaduAI/strawberryfields/pull/720)
 
+* Removes redundant decompositions already available in `thewalrus`
+
 <h3>Bug fixes</h3>
 
 <h3>Documentation</h3>
@@ -18,7 +20,7 @@
 
 This release contains contributions from (in alphabetical order):
 
-Theodor Isacsson
+Theodor Isacsson, Nicolas Quesada
 
 
 # Release 0.23.0 (current release)

requirements.txt

@@ -9,7 +9,7 @@ scipy==1.8.0
 sympy==1.10
 tensorflow==2.8.0
 tensorboard==2.8.0
-thewalrus==0.19.0
+git+https://github.com/xanaduai/thewalrus.git@master
 toml==0.10.2
 typing-extensions==4.1.1
 urllib3==1.26.8

strawberryfields/decompositions.py

@@ -16,90 +16,12 @@
 used to perform gate decompositions.
 """
 
-from itertools import groupby
 from collections import defaultdict
 
 import numpy as np
-from scipy.linalg import block_diag, sqrtm, polar, schur
 from thewalrus.quantum import adj_scaling
-from thewalrus.symplectic import sympmat, xpxp_to_xxpp
-
-
-def takagi(N, tol=1e-13, rounding=13):
-    r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
-
-    Note that singular values of N are considered equal if they are equal after np.round(values, tol).
-
-    See :cite:`cariolaro2016` and references therein for a derivation.
-
-    Args:
-        N (array[complex]): square, symmetric matrix N
-        rounding (int): the number of decimal places to use when rounding the singular values of N
-        tol (float): the tolerance used when checking if the input matrix is symmetric: :math:`|N-N^T| <` tol
-
-    Returns:
-        tuple[array, array]: (rl, U), where rl are the (rounded) singular values,
-            and U is the Takagi unitary, such that :math:`N = U \diag(rl) U^T`.
-    """
-    (n, m) = N.shape
-    if n != m:
-        raise ValueError("The input matrix must be square")
-    if np.linalg.norm(N - np.transpose(N)) >= tol:
-        raise ValueError("The input matrix is not symmetric")
-
-    N = np.real_if_close(N)
-
-    if np.allclose(N, 0):
-        return np.zeros(n), np.eye(n)
-
-    if np.isrealobj(N):
-        # If the matrix N is real one can be more clever and use its eigendecomposition
-        l, U = np.linalg.eigh(N)
-        vals = np.abs(l)  # These are the Takagi eigenvalues
-        phases = np.sqrt(np.complex128([1 if i > 0 else -1 for i in l]))
-        Uc = U @ np.diag(phases)  # One needs to readjust the phases
-        list_vals = [(vals[i], i) for i in range(len(vals))]
-        list_vals.sort(reverse=True)
-        sorted_l, permutation = zip(*list_vals)
-        permutation = np.array(permutation)
-        Uc = Uc[:, permutation]
-        # And also rearrange the unitary and values so that they are decreasingly ordered
-        return np.array(sorted_l), Uc
-
-    v, l, ws = np.linalg.svd(N)
-    w = np.transpose(np.conjugate(ws))
-    rl = np.round(l, rounding)
-
-    # Generate list with degenerancies
-    result = []
-    for k, g in groupby(rl):
-        result.append(list(g))
-
-    # Generate lists containing the columns that correspond to degenerancies
-    kk = 0
-    for k in result:
-        for ind, j in enumerate(k):  # pylint: disable=unused-variable
-            k[ind] = kk
-            kk = kk + 1
-
-    # Generate the lists with the degenerate column subspaces
-    vas = []
-    was = []
-    for i in result:
-        vas.append(v[:, i])
-        was.append(w[:, i])
-
-    # Generate the matrices qs of the degenerate subspaces
-    qs = []
-    for i in range(len(result)):
-        qs.append(sqrtm(np.transpose(vas[i]) @ was[i]))
-
-    # Construct the Takagi unitary
-    qb = block_diag(*qs)
-
-    U = v @ np.conj(qb)
-    return rl, U
-
+from thewalrus.symplectic import sympmat
+from thewalrus.decompositions import takagi
 
 def graph_embed_deprecated(A, max_mean_photon=1.0, make_traceless=False, rtol=1e-05, atol=1e-08):
     r"""Embed a graph into a Gaussian state.
@@ -140,7 +62,7 @@ def graph_embed_deprecated(A, max_mean_photon=1.0, make_traceless=False, rtol=1e
     if make_traceless:
         A = A - np.trace(A) * np.identity(n) / n
 
-    s, U = takagi(A, tol=atol)
+    s, U = takagi(A)
     sc = np.sqrt(1.0 + 1.0 / max_mean_photon)
     vals = -np.arctanh(s / (s[0] * sc))
     return vals, U
@@ -184,7 +106,7 @@ def graph_embed(A, mean_photon_per_mode=1.0, make_traceless=False, rtol=1e-05, a
 
     scale = adj_scaling(A, n * mean_photon_per_mode)
     A = scale * A
-    s, U = takagi(A, tol=atol)
+    s, U = takagi(A)
     vals = -np.arctanh(s)
     return vals, U
 
@@ -903,172 +825,6 @@ def rectangular_compact(U, rtol=1e-12, atol=1e-12):
     return _absorb_zeta(phases_temp)
 
 
-def williamson(V, tol=1e-11):
-    r"""Williamson decomposition of positive-definite (real) symmetric matrix.
-
-    See :ref:`williamson`.
-
-    Note that it is assumed that the symplectic form is
-
-    .. math:: \Omega = \begin{bmatrix}0&I\\-I&0\end{bmatrix}
-
-    where :math:`I` is the identity matrix and :math:`0` is the zero matrix.
-
-    See https://math.stackexchange.com/questions/1171842/finding-the-symplectic-matrix-in-williamsons-theorem/2682630#2682630
-
-    Args:
-        V (array[float]): positive definite symmetric (real) matrix
-        tol (float): the tolerance used when checking if the matrix is symmetric: :math:`|V-V^T| \leq` tol
-
-    Returns:
-        tuple[array,array]: ``(Db, S)`` where ``Db`` is a diagonal matrix
-            and ``S`` is a symplectic matrix such that :math:`V = S^T Db S`
-    """
-    (n, m) = V.shape
-
-    if n != m:
-        raise ValueError("The input matrix is not square")
-
-    diffn = np.linalg.norm(V - np.transpose(V))
-
-    if diffn >= tol:
-        raise ValueError("The input matrix is not symmetric")
-
-    if n % 2 != 0:
-        raise ValueError("The input matrix must have an even number of rows/columns")
-
-    n = n // 2
-    omega = sympmat(n)
-    vals = np.linalg.eigvalsh(V)
-
-    for val in vals:
-        if val <= 0:
-            raise ValueError("Input matrix is not positive definite")
-
-    Mm12 = sqrtm(np.linalg.inv(V)).real
-    r1 = Mm12 @ omega @ Mm12
-    s1, K = schur(r1)
-    X = np.array([[0, 1], [1, 0]])
-    I = np.identity(2)
-    seq = []
-
-    # In what follows I construct a permutation matrix p  so that the Schur matrix has
-    # only positive elements above the diagonal
-    # Also the Schur matrix uses the x_1,p_1, ..., x_n,p_n  ordering thus I use rotmat to
-    # go to the ordering x_1, ..., x_n, p_1, ... , p_n
-
-    for i in range(n):
-        if s1[2 * i, 2 * i + 1] > 0:
-            seq.append(I)
-        else:
-            seq.append(X)
-
-    p = block_diag(*seq)
-    Kt = K @ p
-    s1t = p @ s1 @ p
-    dd = xpxp_to_xxpp(s1t)
-    perm_indices = xpxp_to_xxpp(np.arange(2 * n))
-    Ktt = Kt[:, perm_indices]
-    Db = np.diag([1 / dd[i, i + n] for i in range(n)] + [1 / dd[i, i + n] for i in range(n)])
-    S = Mm12 @ Ktt @ sqrtm(Db)
-    return Db, np.linalg.inv(S).T
-
-
-def bloch_messiah(S, tol=1e-10, rounding=9):
-    r"""Bloch-Messiah decomposition of a symplectic matrix.
-
-    See :ref:`bloch_messiah`.
-
-    Decomposes a symplectic matrix into two symplectic unitaries and squeezing transformation.
-    It automatically sorts the squeezers so that they respect the canonical symplectic form.
-
-    Note that it is assumed that the symplectic form is
-
-    .. math:: \Omega = \begin{bmatrix}0&I\\-I&0\end{bmatrix}
-
-    where :math:`I` is the identity matrix and :math:`0` is the zero matrix.
-
-    As in the Takagi decomposition, the singular values of N are considered
-    equal if they are equal after np.round(values, rounding).
-
-    If S is a passive transformation, then return the S as the first passive
-    transformation, and set the the squeezing and second unitary matrices to
-    identity. This choice is not unique.
-
-    For more info see:
-    https://math.stackexchange.com/questions/1886038/finding-euler-decomposition-of-a-symplectic-matrix
-
-    Args:
-        S (array[float]): symplectic matrix
-        tol (float): the tolerance used when checking if the matrix is symplectic:
-            :math:`|S^T\Omega S-\Omega| \leq tol`
-        rounding (int): the number of decimal places to use when rounding the singular values
-
-    Returns:
-        tuple[array]: Returns the tuple ``(ut1, st1, vt1)``. ``ut1`` and ``vt1`` are symplectic orthogonal,
-            and ``st1`` is diagonal and of the form :math:`= \text{diag}(s1,\dots,s_n, 1/s_1,\dots,1/s_n)`
-            such that :math:`S = ut1  st1  v1`
-    """
-    (n, m) = S.shape
-
-    if n != m:
-        raise ValueError("The input matrix is not square")
-    if n % 2 != 0:
-        raise ValueError("The input matrix must have an even number of rows/columns")
-
-    n = n // 2
-    omega = sympmat(n)
-    if np.linalg.norm(np.transpose(S) @ omega @ S - omega) >= tol:
-        raise ValueError("The input matrix is not symplectic")
-
-    if np.linalg.norm(np.transpose(S) @ S - np.eye(2 * n)) >= tol:
-
-        u, sigma = polar(S, side="left")
-        ss, uss = takagi(sigma, tol=tol, rounding=rounding)
-
-        # Apply a permutation matrix so that the squeezers appear in the order
-        # s_1,...,s_n, 1/s_1,...1/s_n
-        perm = np.array(list(range(0, n)) + list(reversed(range(n, 2 * n))))
-
-        pmat = np.identity(2 * n)[perm, :]
-        ut = uss @ pmat
-
-        # Apply a second permutation matrix to permute s
-        # (and their corresonding inverses) to get the canonical symplectic form
-        qomega = np.transpose(ut) @ (omega) @ ut
-        st = pmat @ np.diag(ss) @ pmat
-
-        # Identifying degenerate subspaces
-        result = []
-        for _k, g in groupby(np.round(np.diag(st), rounding)[:n]):
-            result.append(list(g))
-
-        stop_is = list(np.cumsum([len(res) for res in result]))
-        start_is = [0] + stop_is[:-1]
-
-        # Rotation matrices (not permutations) based on svd.
-        # See Appendix B2 of Serafini's book for more details.
-        u_list, v_list = [], []
-
-        for start_i, stop_i in zip(start_is, stop_is):
-            x = qomega[start_i:stop_i, n + start_i : n + stop_i].real
-            u_svd, _s_svd, v_svd = np.linalg.svd(x)
-            u_list = u_list + [u_svd]
-            v_list = v_list + [v_svd.T]
-
-        pmat1 = block_diag(*(u_list + v_list))
-
-        st1 = pmat1.T @ pmat @ np.diag(ss) @ pmat @ pmat1
-        ut1 = uss @ pmat @ pmat1
-        v1 = np.transpose(ut1) @ u
-
-    else:
-        ut1 = S
-        st1 = np.eye(2 * n)
-        v1 = np.eye(2 * n)
-
-    return ut1.real, st1.real, v1.real
-
 
 def sun_compact(U, rtol=1e-12, atol=1e-12):
     r"""Recursive factorization of unitary transfomations.

strawberryfields/ops.py

@@ -26,6 +26,7 @@
 import scipy.special as ssp
 
 from thewalrus.symplectic import xxpp_to_xpxp
+from thewalrus.decompositions import blochmessiah, williamson
 
 import strawberryfields as sf
 import strawberryfields.program_utils as pu
@@ -2933,7 +2934,7 @@ class GaussianTransform(Decomposition):
         transformation (:math:`S = O_1 Z`).
     """
 
-    def __init__(self, S, vacuum=False, tol=1e-10):
+    def __init__(self, S, vacuum=False):
         super().__init__([S])
         self.ns = S.shape[0] // 2
         self.vacuum = (
@@ -2954,7 +2955,7 @@ def __init__(self, S, vacuum=False, tol=1e-10):
             self.U1 = X1 + 1j * P1
         else:
             # transformation is active, do Bloch-Messiah
-            O1, smat, O2 = dec.bloch_messiah(S, tol=tol)
+            O1, smat, O2 = blochmessiah(S)
             X1 = O1[:N, :N]
             P1 = O1[N:, :N]
             X2 = O2[:N, :N]
@@ -3072,7 +3073,7 @@ def __init__(self, V, r=None, decomp=True, tol=1e-6):
         self.p_disp = r[self.ns :]
 
         # needed only if decomposed
-        th, self.S = dec.williamson(V, tol=tol)
+        th, self.S = williamson(V, rtol=tol)
         self.pure = np.abs(np.linalg.det(V) - 1.0) < tol
         self.nbar = 0.5 * (np.diag(th)[: self.ns] - 1.0)
 

strawberryfields/utils/random_numbers_matrices.py

@@ -110,12 +110,10 @@ def random_interferometer(N, real=False):
     Returns:
         array: random :math:`N\times N` unitary distributed with the Haar measure
     """
+    if N == 1:
+        if real:
+            return np.array([[2 * (np.random.binomial(1, 0.5) - 0.5)]])
+        return np.array([[np.exp(1j * 2 * np.pi * np.random.rand())]])
     if real:
-        z = np.random.randn(N, N)
-    else:
-        z = randnc(N, N) / np.sqrt(2.0)
-    q, r = sp.linalg.qr(z)
-    d = np.diagonal(r)
-    ph = d / np.abs(d)
-    U = np.multiply(q, ph, q)
-    return U
+        return sp.stats.ortho_group.rvs(N)
+    return sp.stats.unitary_group.rvs(N)