1 files changed, 75 insertions, 0 deletions
diff --git a/python_originals/takagi.py b/python_originals/takagi.py
new file mode 100644
index 0000000..4f3fbd9
--- /dev/null
+++ b/python_originals/takagi.py
@@ -0,0 +1,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
+