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
|
function [Db, S] = williamson(V, tol)
% Williamson decomposition of positive-definite (real) symmetric matrix.
%
% Args:
% V (matrix): positive definite symmetric (real) matrix
% tol (scalar): tolerance for symmetry check (default: 1e-11)
%
% Returns:
% Db (matrix): diagonal matrix
% S (matrix): symplectic matrix such that V = S' * Db * S
if nargin < 2
tol = 1e-11;
end
[n, m] = size(V);
if n ~= m
error('The input matrix is not square');
end
diffn = norm(V - V');
if diffn >= tol
error('The input matrix is not symmetric');
end
if mod(n, 2) ~= 0
error('The input matrix must have an even number of rows/columns');
end
n = n / 2;
omega = sympmat(n);
vals = eig(V);
if any(vals <= 0)
error('Input matrix is not positive definite');
end
Mm12 = real(sqrtm(inv(V)));
r1 = Mm12 * omega * Mm12;
[K, s1] = schur(r1);
X = [0 1; 1 0];
I = eye(2);
seq = cell(1, n);
% Construct a permutation matrix p
for i = 1:n
if s1(2*i-1, 2*i) > 0
seq{i} = I;
else
seq{i} = X;
end
end
p = blkdiag(seq{:});
Kt = K * p;
s1t = p * s1 * p;
dd = xpxp_to_xxpp(s1t);
perm_indices = xpxp_to_xxpp(1:2*n);
Ktt = Kt(:, perm_indices);
Db_diag = [1 ./ diag(dd(1:n, n+1:end)); 1 ./ diag(dd(1:n, n+1:end))];
Db = diag(Db_diag);
S = Mm12 * Ktt * sqrtm(Db);
S = inv(S)';
end
%!test
%! V = load('test_data/covmat0.csv');
%! [Db, S] = williamson(V);
%! assert(Db,eye(8), 1e-14)
%!test
%! V = load('test_data/covmat15.csv');
%! D0 = load('test_data/w_diag15.csv');
%! [Db, S] = williamson(V);
%! assert(Db,D0, 1e-9)
|
