1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
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
|
