aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Readme.md6
-rw-r--r--python_src/williamson.py70
-rw-r--r--python_src/xpxp_to_xxpp.py26
-rw-r--r--python_src/xxpp_to_xpxp.py25
-rw-r--r--sympmat.m12
-rw-r--r--williamson.m73
-rw-r--r--xpxp_to_xxpp.m32
-rw-r--r--xxpp_to_xpxp.m30
8 files changed, 274 insertions, 0 deletions
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
+
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-25 06:45:31 +0000