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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
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
|
