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
|
def takagi(N, tol=1e-13, rounding=13):
r"""Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
Note that singular values of N are considered equal if they are equal after np.round(values, tol).
See :cite:`cariolaro2016` and references therein for a derivation.
Args:
N (array[complex]): square, symmetric matrix N
rounding (int): the number of decimal places to use when rounding the singular values of N
tol (float): the tolerance used when checking if the input matrix is symmetric: :math:`|N-N^T| <` tol
Returns:
tuple[array, array]: (rl, U), where rl are the (rounded) singular values,
and U is the Takagi unitary, such that :math:`N = U \diag(rl) U^T`.
"""
(n, m) = N.shape
if n != m:
raise ValueError("The input matrix must be square")
if np.linalg.norm(N - np.transpose(N)) >= tol:
raise ValueError("The input matrix is not symmetric")
N = np.real_if_close(N)
if np.allclose(N, 0):
return np.zeros(n), np.eye(n)
if np.isrealobj(N):
# If the matrix N is real one can be more clever and use its eigendecomposition
l, U = np.linalg.eigh(N)
vals = np.abs(l) # These are the Takagi eigenvalues
phases = np.sqrt(np.complex128([1 if i > 0 else -1 for i in l]))
Uc = U @ np.diag(phases) # One needs to readjust the phases
list_vals = [(vals[i], i) for i in range(len(vals))]
list_vals.sort(reverse=True)
sorted_l, permutation = zip(*list_vals)
permutation = np.array(permutation)
Uc = Uc[:, permutation]
# And also rearrange the unitary and values so that they are decreasingly ordered
return np.array(sorted_l), Uc
v, l, ws = np.linalg.svd(N)
w = np.transpose(np.conjugate(ws))
rl = np.round(l, rounding)
# Generate list with degenerancies
result = []
for k, g in groupby(rl):
result.append(list(g))
# Generate lists containing the columns that correspond to degenerancies
kk = 0
for k in result:
for ind, j in enumerate(k): # pylint: disable=unused-variable
k[ind] = kk
kk = kk + 1
# Generate the lists with the degenerate column subspaces
vas = []
was = []
for i in result:
vas.append(v[:, i])
was.append(w[:, i])
# Generate the matrices qs of the degenerate subspaces
qs = []
for i in range(len(result)):
qs.append(sqrtm(np.transpose(vas[i]) @ was[i]))
# Construct the Takagi unitary
qb = block_diag(*qs)
U = v @ np.conj(qb)
return rl, U
|
