From 984d56b4d94bf298a197f1567834451acf49cab9 Mon Sep 17 00:00:00 2001 From: "Eugeniy E. Mikhailov" Date: Tue, 19 Nov 2024 13:40:38 -0500 Subject: dir rename --- python_originals/williamson.py | 70 ++++++++++++++++++++++++++++++++++++++++ python_originals/xpxp_to_xxpp.py | 26 +++++++++++++++ python_originals/xxpp_to_xpxp.py | 25 ++++++++++++++ python_src/williamson.py | 70 ---------------------------------------- python_src/xpxp_to_xxpp.py | 26 --------------- python_src/xxpp_to_xpxp.py | 25 -------------- 6 files changed, 121 insertions(+), 121 deletions(-) create mode 100644 python_originals/williamson.py create mode 100644 python_originals/xpxp_to_xxpp.py create mode 100644 python_originals/xxpp_to_xpxp.py delete mode 100644 python_src/williamson.py delete mode 100644 python_src/xpxp_to_xxpp.py delete mode 100644 python_src/xxpp_to_xpxp.py diff --git a/python_originals/williamson.py b/python_originals/williamson.py new file mode 100644 index 0000000..19a2f5f --- /dev/null +++ b/python_originals/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_originals/xpxp_to_xxpp.py b/python_originals/xpxp_to_xxpp.py new file mode 100644 index 0000000..01d0bf2 --- /dev/null +++ b/python_originals/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_originals/xxpp_to_xpxp.py b/python_originals/xxpp_to_xpxp.py new file mode 100644 index 0000000..940d8b8 --- /dev/null +++ b/python_originals/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/python_src/williamson.py b/python_src/williamson.py deleted file mode 100644 index 19a2f5f..0000000 --- a/python_src/williamson.py +++ /dev/null @@ -1,70 +0,0 @@ -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 deleted file mode 100644 index 01d0bf2..0000000 --- a/python_src/xpxp_to_xxpp.py +++ /dev/null @@ -1,26 +0,0 @@ - -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 deleted file mode 100644 index 940d8b8..0000000 --- a/python_src/xxpp_to_xpxp.py +++ /dev/null @@ -1,25 +0,0 @@ -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] - -- cgit v1.2.3