draft of bloch_messiah

3 files changed, 276 insertions, 0 deletions

diff --git a/python_originals/bloch_messiah.py b/python_originals/bloch_messiah.py
new file mode 100644
index 0000000..a4cc310
--- /dev/null
+++ b/python_originals/bloch_messiah.py
@@ -0,0 +1,96 @@
+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
+
+
diff --git a/python_originals/polar.py b/python_originals/polar.py
new file mode 100644
index 0000000..69b6a35
--- /dev/null
+++ b/python_originals/polar.py
@@ -0,0 +1,105 @@
+def polar(a, side="right"):
+    """
+    Compute the polar decomposition.
+
+    Returns the factors of the polar decomposition [1]_ `u` and `p` such
+    that ``a = up`` (if `side` is "right") or ``a = pu`` (if `side` is
+    "left"), where `p` is positive semidefinite. Depending on the shape
+    of `a`, either the rows or columns of `u` are orthonormal. When `a`
+    is a square array, `u` is a square unitary array. When `a` is not
+    square, the "canonical polar decomposition" [2]_ is computed.
+
+    Parameters
+    ----------
+    a : (m, n) array_like
+        The array to be factored.
+    side : {'left', 'right'}, optional
+        Determines whether a right or left polar decomposition is computed.
+        If `side` is "right", then ``a = up``.  If `side` is "left",  then
+        ``a = pu``.  The default is "right".
+
+    Returns
+    -------
+    u : (m, n) ndarray
+        If `a` is square, then `u` is unitary. If m > n, then the columns
+        of `a` are orthonormal, and if m < n, then the rows of `u` are
+        orthonormal.
+    p : ndarray
+        `p` is Hermitian positive semidefinite. If `a` is nonsingular, `p`
+        is positive definite. The shape of `p` is (n, n) or (m, m), depending
+        on whether `side` is "right" or "left", respectively.
+
+    References
+    ----------
+    .. [1] R. A. Horn and C. R. Johnson, "Matrix Analysis", Cambridge
+           University Press, 1985.
+    .. [2] N. J. Higham, "Functions of Matrices: Theory and Computation",
+           SIAM, 2008.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from scipy.linalg import polar
+    >>> a = np.array([[1, -1], [2, 4]])
+    >>> u, p = polar(a)
+    >>> u
+    array([[ 0.85749293, -0.51449576],
+           [ 0.51449576,  0.85749293]])
+    >>> p
+    array([[ 1.88648444,  1.2004901 ],
+           [ 1.2004901 ,  3.94446746]])
+
+    A non-square example, with m < n:
+
+    >>> b = np.array([[0.5, 1, 2], [1.5, 3, 4]])
+    >>> u, p = polar(b)
+    >>> u
+    array([[-0.21196618, -0.42393237,  0.88054056],
+           [ 0.39378971,  0.78757942,  0.4739708 ]])
+    >>> p
+    array([[ 0.48470147,  0.96940295,  1.15122648],
+           [ 0.96940295,  1.9388059 ,  2.30245295],
+           [ 1.15122648,  2.30245295,  3.65696431]])
+    >>> u.dot(p)   # Verify the decomposition.
+    array([[ 0.5,  1. ,  2. ],
+           [ 1.5,  3. ,  4. ]])
+    >>> u.dot(u.T)   # The rows of u are orthonormal.
+    array([[  1.00000000e+00,  -2.07353665e-17],
+           [ -2.07353665e-17,   1.00000000e+00]])
+
+    Another non-square example, with m > n:
+
+    >>> c = b.T
+    >>> u, p = polar(c)
+    >>> u
+    array([[-0.21196618,  0.39378971],
+           [-0.42393237,  0.78757942],
+           [ 0.88054056,  0.4739708 ]])
+    >>> p
+    array([[ 1.23116567,  1.93241587],
+           [ 1.93241587,  4.84930602]])
+    >>> u.dot(p)   # Verify the decomposition.
+    array([[ 0.5,  1.5],
+           [ 1. ,  3. ],
+           [ 2. ,  4. ]])
+    >>> u.T.dot(u)  # The columns of u are orthonormal.
+    array([[  1.00000000e+00,  -1.26363763e-16],
+           [ -1.26363763e-16,   1.00000000e+00]])
+
+    """
+    if side not in ['right', 'left']:
+        raise ValueError("`side` must be either 'right' or 'left'")
+    a = np.asarray(a)
+    if a.ndim != 2:
+        raise ValueError("`a` must be a 2-D array.")
+
+    w, s, vh = svd(a, full_matrices=False)
+    u = w.dot(vh)
+    if side == 'right':
+        # a = up
+        p = (vh.T.conj() * s).dot(vh)
+    else:
+        # a = pu
+        p = (w * s).dot(w.T.conj())
+    return u, p
+
diff --git a/python_originals/takagi.py b/python_originals/takagi.py
new file mode 100644
index 0000000..4f3fbd9
--- /dev/null
+++ b/python_originals/takagi.py
@@ -0,0 +1,75 @@
+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
+