From 6187aede8dc13a4fc4d5fd5463574955d702d7f4 Mon Sep 17 00:00:00 2001 From: "Eugeniy E. Mikhailov" Date: Tue, 19 Nov 2024 13:00:14 -0500 Subject: initial release --- Readme.md | 6 ++++ python_src/williamson.py | 70 ++++++++++++++++++++++++++++++++++++++++++++ python_src/xpxp_to_xxpp.py | 26 +++++++++++++++++ python_src/xxpp_to_xpxp.py | 25 ++++++++++++++++ sympmat.m | 12 ++++++++ williamson.m | 73 ++++++++++++++++++++++++++++++++++++++++++++++ xpxp_to_xxpp.m | 32 ++++++++++++++++++++ xxpp_to_xpxp.m | 30 +++++++++++++++++++ 8 files changed, 274 insertions(+) create mode 100644 Readme.md create mode 100644 python_src/williamson.py create mode 100644 python_src/xpxp_to_xxpp.py create mode 100644 python_src/xxpp_to_xpxp.py create mode 100644 sympmat.m create mode 100644 williamson.m create mode 100644 xpxp_to_xxpp.m create mode 100644 xxpp_to_xpxp.m diff --git a/Readme.md b/Readme.md new file mode 100644 index 0000000..d64b1ef --- /dev/null +++ b/Readme.md @@ -0,0 +1,6 @@ +This is matlab implementation of some function +from python package [Strawberry Fields](https://strawberryfields.ai/) + +The initial conversion was done with Perplexity AI, but +it required a human intervention for xpxp_to_xxpp code. + diff --git a/python_src/williamson.py b/python_src/williamson.py new file mode 100644 index 0000000..19a2f5f --- /dev/null +++ b/python_src/williamson.py @@ -0,0 +1,70 @@ +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 + diff --git a/python_src/xpxp_to_xxpp.py b/python_src/xpxp_to_xxpp.py new file mode 100644 index 0000000..01d0bf2 --- /dev/null +++ b/python_src/xpxp_to_xxpp.py @@ -0,0 +1,26 @@ + +def xpxp_to_xxpp(S): + """Permutes the entries of the input from xpxp ordering to xxpp ordering. + + Args: + S (array): input even dimensional square matrix or vector + + Returns: + (array): permuted matrix or vector + """ + shape = S.shape + n = shape[0] + + if n % 2 != 0: + raise ValueError("The input array is not even-dimensional") + + n = n // 2 + ind = np.arange(2 * n).reshape(-1, 2).T.flatten() + + if len(shape) == 2: + if shape[0] != shape[1]: + raise ValueError("The input matrix is not square") + return S[:, ind][ind] + + return S[ind] + diff --git a/python_src/xxpp_to_xpxp.py b/python_src/xxpp_to_xpxp.py new file mode 100644 index 0000000..940d8b8 --- /dev/null +++ b/python_src/xxpp_to_xpxp.py @@ -0,0 +1,25 @@ +def xxpp_to_xpxp(S): + """Permutes the entries of the input from xxpp ordering to xpxp ordering. + + Args: + S (array): input even dimensional square matrix or array + + Returns: + (array): permuted matrix or array + """ + shape = S.shape + n = shape[0] + + if n % 2 != 0: + raise ValueError("The input array is not even-dimensional") + + n = n // 2 + ind = np.arange(2 * n).reshape(2, -1).T.flatten() + + if len(shape) == 2: + if shape[0] != shape[1]: + raise ValueError("The input matrix is not square") + return S[:, ind][ind] + + return S[ind] + diff --git a/sympmat.m b/sympmat.m new file mode 100644 index 0000000..9eb84db --- /dev/null +++ b/sympmat.m @@ -0,0 +1,12 @@ +function omega = sympmat(n) + % Create a symplectic form matrix + % 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. + I = eye(n); + Z = zeros(n); + omega = [Z I; -I Z]; +end + diff --git a/williamson.m b/williamson.m new file mode 100644 index 0000000..f969ab1 --- /dev/null +++ b/williamson.m @@ -0,0 +1,73 @@ +function [Db, S] = williamson(V, tol) + % Williamson decomposition of positive-definite (real) symmetric matrix. + % + % Args: + % V (matrix): positive definite symmetric (real) matrix + % tol (scalar): tolerance for symmetry check (default: 1e-11) + % + % Returns: + % Db (matrix): diagonal matrix + % S (matrix): symplectic matrix such that V = S' * Db * S + + if nargin < 2 + tol = 1e-11; + end + + [n, m] = size(V); + + if n ~= m + error('The input matrix is not square'); + end + + diffn = norm(V - V'); + if diffn >= tol + error('The input matrix is not symmetric'); + 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); + + vals = eig(V); + if any(vals <= 0) + error('Input matrix is not positive definite'); + end + + Mm12 = real(sqrtm(inv(V))); + r1 = Mm12 * omega * Mm12; + [K, s1] = schur(r1); + + X = [0 1; 1 0]; + I = eye(2); + seq = cell(1, n); + + % Construct a permutation matrix p + for i = 1:n + if s1(2*i-1, 2*i) > 0 + seq{i} = I; + else + seq{i} = X; + end + end + p = blkdiag(seq{:}); + + Kt = K * p; + s1t = p * s1 * p; + s1t + + dd = xpxp_to_xxpp(s1t); + perm_indices = xpxp_to_xxpp(1:2*n); + Ktt = Kt(:, perm_indices); + + Db_diag = [1 ./ diag(dd(1:n, n+1:end)); 1 ./ diag(dd(1:n, n+1:end))]; + Db = diag(Db_diag); + + S = Mm12 * Ktt * sqrtm(Db); + S = inv(S)'; +end + + diff --git a/xpxp_to_xxpp.m b/xpxp_to_xxpp.m new file mode 100644 index 0000000..fc732e2 --- /dev/null +++ b/xpxp_to_xxpp.m @@ -0,0 +1,32 @@ +function S_permuted = xpxp_to_xxpp(S) + % Permutes the entries of the input from xpxp ordering to xxpp ordering. + % + % Args: + % S (array): input even dimensional square matrix or a row vector + % + % Returns: + % S_permuted (array): permuted matrix or vector + + shape = size(S); + if length(shape) > 2 + error('The input matrix is not 2 dimensional'); + end + + n = shape(2); + + if mod(n, 2) ~= 0 + error('The input array is not even-dimensional'); + end + + n = n / 2; + ind = reshape(1:2*n, 2, [])'; + ind = ind(:)'; + + + if shape(1) == shape(2) + S_permuted = S(ind, ind); + else + S_permuted = S(ind); + end +end + diff --git a/xxpp_to_xpxp.m b/xxpp_to_xpxp.m new file mode 100644 index 0000000..d6efbe9 --- /dev/null +++ b/xxpp_to_xpxp.m @@ -0,0 +1,30 @@ +function S_permuted = xxpp_to_xpxp(S) + % Permutes the entries of the input from xxpp ordering to xpxp ordering. + % + % Args: + % S (array): input even dimensional square matrix or array + % + % Returns: + % S_permuted (array): permuted matrix or array + + shape = size(S); + n = shape(1); + + if mod(n, 2) ~= 0 + error('The input array is not even-dimensional'); + end + + n = n / 2; + ind = reshape(1:2*n, 2, [])'; + ind = ind(:); + + if length(shape) == 2 + if shape(1) ~= shape(2) + error('The input matrix is not square'); + end + S_permuted = S(ind, ind); + else + S_permuted = S(ind); + end +end + -- cgit v1.2.3