1 files changed, 77 insertions, 0 deletions
diff --git a/takagi.m b/takagi.m
new file mode 100644
index 0000000..f9375a8
--- /dev/null
+++ b/takagi.m
@@ -0,0 +1,77 @@
+function [rl, U] = takagi(N, tol, rounding)
+    % Autonne-Takagi decomposition of a complex symmetric (not Hermitian!) matrix.
+    %
+    % Args:
+    %   N (complex matrix): square, symmetric matrix N
+    %   tol (double): tolerance for symmetry check (default: 1e-13)
+    %   rounding (integer): decimal places for rounding singular values (default: 13)
+    %
+    % Returns:
+    %   rl (vector): rounded singular values
+    %   U (matrix): Takagi unitary, such that N = U * diag(rl) * U.'
+
+    if nargin < 2
+        tol = 1e-13;
+    end
+    if nargin < 3
+        rounding = 13;
+    end
+
+    [n, m] = size(N);
+
+    if n ~= m
+        error('The input matrix must be square');
+    end
+
+    if norm(N - N.') >= tol
+        error('The input matrix is not symmetric');
+    end
+
+    if all(abs(imag(N(:))) < eps)
+        N = real(N);
+    end
+
+    if all(N(:) == 0)
+        rl = zeros(n, 1);
+        U = eye(n);
+        return;
+    end
+
+    if isreal(N)
+        % If the matrix N is real, use its eigendecomposition
+        [U, L] = eig(N);
+        l = diag(L);
+        vals = abs(l);
+        phases = sqrt(complex(sign(l)));
+        Uc = U * diag(phases);
+
+        [sorted_l, idx] = sort(vals, 'descend');
+        Uc = Uc(:, idx);
+
+        rl = sorted_l;
+        U = Uc;
+        return;
+    end
+
+    [v, l, ws] = svd(N);
+    w = ws';
+    rl = round(diag(l), rounding);
+
+    % Group degenerate singular values
+    [~, ~, ic] = unique(rl);
+    result = accumarray(ic, (1:numel(rl))', [], @(x) {x});
+
+    % Generate the matrices qs of the degenerate subspaces
+    qs = cell(1, length(result));
+    for i = 1:length(result)
+        idx = result{i};
+        vas = v(:, idx);
+        was = w(:, idx);
+        qs{i} = sqrtm(vas' * was);
+    end
+
+    % Construct the Takagi unitary
+    qb = blkdiag(qs{:});
+    U = v * conj(qb);
+end
+