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
|
function [rl, U] = takagi(N, tol)
% 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)
%
% Returns:
% rl (vector): rounded singular values
% U (matrix): Takagi unitary, such that N = U * diag(rl) * U.'
if nargin < 2
tol = 1e-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 = diag(l);
% 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
|
