draft of bloch_messiah

6 files changed, 494 insertions, 0 deletions

diff --git a/bloch_messiah.m b/bloch_messiah.m
new file mode 100644
index 0000000..6724a26
--- /dev/null
+++ b/bloch_messiah.m
@@ -0,0 +1,85 @@
+function [ut1, st1, v1] = bloch_messiah(S, tol, rounding)
+    % Bloch-Messiah decomposition of a symplectic matrix.
+    %
+    % Args:
+    %   S (matrix): symplectic matrix
+    %   tol (double): tolerance for symplectic check (default: 1e-10)
+    %   rounding (int): decimal places for rounding singular values (default: 9)
+    %
+    % Returns:
+    %   ut1, st1, v1 (matrices): Decomposition matrices such that S = ut1 * st1 * v1
+
+    if nargin < 2
+        tol = 1e-10;
+    end
+    if nargin < 3
+        rounding = 9;
+    end
+
+    [n, m] = size(S);
+
+    if n ~= m
+        error('The input matrix is not square');
+    end
+
+    if mod(n, 2) ~= 0
+        error('The input matrix must have an even number of rows/columns');
+    end
+
+    n = n / 2;
+    omega = sympmat(n);
+
+    if norm(S' * omega * S - omega) >= tol
+        error('The input matrix is not symplectic');
+    end
+
+    if norm(S' * S - eye(2*n)) >= tol
+        [u, sigma] = polardecomp(S, 'left');
+        [ss, uss] = takagi(sigma, tol, rounding);
+
+        % Apply permutation matrix
+        perm = [1:n, 2*n:-1:n+1];
+        pmat = eye(2*n);
+        pmat = pmat(perm, :);
+
+        ut = uss * pmat;
+
+        % Apply second permutation matrix
+        qomega = ut' * omega * ut;
+        st = pmat * diag(ss) * pmat;
+
+        % Identify degenerate subspaces
+        st_diag = round(diag(st), rounding);
+        [~, ~, ic] = unique(st_diag(1:n));
+        stop_is = cumsum(accumarray(ic, 1));
+        start_is = [0; stop_is(1:end-1)] + 1;
+
+        % Rotation matrices based on SVD
+        u_list = cell(1, length(start_is));
+        v_list = cell(1, length(start_is));
+
+        for i = 1:length(start_is)
+            start_i = start_is(i);
+            stop_i = stop_is(i);
+            x = real(qomega(start_i:stop_i, n+start_i:n+stop_i));
+            [u_svd, ~, v_svd] = svd(x);
+            u_list{i} = u_svd;
+            v_list{i} = v_svd';
+        end
+
+        pmat1 = blkdiag(u_list{:}, v_list{:});
+
+        st1 = pmat1' * pmat * diag(ss) * pmat * pmat1;
+        ut1 = uss * pmat * pmat1;
+        v1 = ut1' * u;
+    else
+        ut1 = S;
+        st1 = eye(2*n);
+        v1 = eye(2*n);
+    end
+
+    ut1 = real(ut1);
+    st1 = real(st1);
+    v1 = real(v1);
+end
+
diff --git a/polardecomp.m b/polardecomp.m
new file mode 100644
index 0000000..698bd90
--- /dev/null
+++ b/polardecomp.m
@@ -0,0 +1,56 @@
+function [u, p] = polardecomp(a, side)
+    % 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: The array to be factored.
+    %   side: (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: 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: `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.
+    
+    if nargin < 2
+        side = 'right';
+    end
+    
+    if ~strcmp(side, 'right') && ~strcmp(side, 'left')
+        error('`side` must be either ''right'' or ''left''');
+    end
+    
+    if ~ismatrix(a)
+        error('`a` must be a 2-D array.');
+    end
+    
+    [w, s, vh] = svd(a, 'econ');
+    u = w * vh';
+    
+    if strcmp(side, 'right')
+        % a = up
+        p = (vh' * diag(diag(s))) * vh;
+    else
+        % a = pu
+        p = (w * diag(diag(s))) * w';
+    end
+end
+
+%!test
+%! a = [1,-1; 2, 4];
+%! [u, p] = polardecomp(a);
+%! u0 = [0.85749293,  -0.51449576 ;  0.51449576,   0.85749293];
+%! p0 = [1.88648444,  1.2004901; 1.2004901,  3.94446746];
+%! assert(u,u0, 1e-8)
+%! assert(p,p0, 1e-8)
+%! assert(u*u', eye(2), 1e-15)
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
+
diff --git a/takagi.m b/takagi.m
new file mode 100644
index 0000000..f9375a8
--- /dev/null
+++ b/takagi.m
@@ -0,0 +1,77 @@
+function [rl, U] = takagi(N, tol, rounding)
+    % Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
+    %
+    % Args:
+    %   N (complex matrix): square, symmetric matrix N
+    %   tol (double): tolerance for symmetry check (default: 1e-13)
+    %   rounding (integer): decimal places for rounding singular values (default: 13)
+    %
+    % Returns:
+    %   rl (vector): rounded singular values
+    %   U (matrix): Takagi unitary, such that N = U * diag(rl) * U.'
+
+    if nargin < 2
+        tol = 1e-13;
+    end
+    if nargin < 3
+        rounding = 13;
+    end
+
+    [n, m] = size(N);
+
+    if n ~= m
+        error('The input matrix must be square');
+    end
+
+    if norm(N - N.') >= tol
+        error('The input matrix is not symmetric');
+    end
+
+    if all(abs(imag(N(:))) < eps)
+        N = real(N);
+    end
+
+    if all(N(:) == 0)
+        rl = zeros(n, 1);
+        U = eye(n);
+        return;
+    end
+
+    if isreal(N)
+        % If the matrix N is real, use its eigendecomposition
+        [U, L] = eig(N);
+        l = diag(L);
+        vals = abs(l);
+        phases = sqrt(complex(sign(l)));
+        Uc = U * diag(phases);
+
+        [sorted_l, idx] = sort(vals, 'descend');
+        Uc = Uc(:, idx);
+
+        rl = sorted_l;
+        U = Uc;
+        return;
+    end
+
+    [v, l, ws] = svd(N);
+    w = ws';
+    rl = round(diag(l), rounding);
+
+    % Group degenerate singular values
+    [~, ~, ic] = unique(rl);
+    result = accumarray(ic, (1:numel(rl))', [], @(x) {x});
+
+    % Generate the matrices qs of the degenerate subspaces
+    qs = cell(1, length(result));
+    for i = 1:length(result)
+        idx = result{i};
+        vas = v(:, idx);
+        was = w(:, idx);
+        qs{i} = sqrtm(vas' * was);
+    end
+
+    % Construct the Takagi unitary
+    qb = blkdiag(qs{:});
+    U = v * conj(qb);
+end
+