diff options
| -rw-r--r-- | bloch_messiah.m | 85 | ||||
| -rw-r--r-- | polardecomp.m | 56 | ||||
| -rw-r--r-- | python_originals/bloch_messiah.py | 96 | ||||
| -rw-r--r-- | python_originals/polar.py | 105 | ||||
| -rw-r--r-- | python_originals/takagi.py | 75 | ||||
| -rw-r--r-- | takagi.m | 77 |
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 + |