-
Notifications
You must be signed in to change notification settings - Fork 10
/
matrix_sqrt.m
118 lines (96 loc) · 3.07 KB
/
matrix_sqrt.m
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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
% Set parameters for testing
n = 512;
tau = 1;
numPoints = 1000;
isFloat = true;
isGPU = true;
maxIter = 10;
runIter = 10;
fprintf('Experiment:\n')
fprintf('>>>> n %i, tau %.1f, numPoints %i, float %i, gpu %i, maxIter %i, runIter %i\n', ...
n, tau, numPoints, isFloat, isGPU, maxIter, runIter);
% Generate points randomly to create a PSD matrix
X = randn(numPoints, n);
A = (X'*X)/size(X,1) + tau*eye(n);
dldz = rand(n);
dldz = 0.5*(dldz + dldz');
% Calculate quantities with full precision on the CPU
sqrtA_true = sqrtm(A);
dlda_true = lyap2(sqrtA_true, -dldz);
if isFloat
A = single(A);
dldz = single(dldz);
fprintf('>>>> floating point: single\n');
else
fprintf('>>>> floating point: double\n');
end
if isGPU
A = gpuArray(A);
dldz = gpuArray(dldz);
fprintf('>>>> GPU: true\n');
else
fprintf('>>>> GPU: false\n');
end
% Evaluate forward models
% SVD
fprintf('\nForward time:\n');
tic
for i = 1:runIter
[sqrtA_svd, cache] = sqrt_forward(A, 'svd');
end
fprintf(' %fs SVD\n', toc);
% Denman Beavers
tic
for i = 1:runIter
[sqrtA_db, ~] = sqrt_forward(A, 'db', maxIter);
end
fprintf(' %fs Denman-Beavers (%i iter)\n', toc, maxIter);
% Newton Schulz
tic
for i = 1:runIter
[sqrtA_ns, ~] = sqrt_forward(A, 'ns', maxIter);
end
fprintf(' %fs Newton-Schulz (%i iter)\n', toc, maxIter);
% Evalute backward models
fprintf('\nBackward time:\n');
tic
for i = 1:runIter
dlda_mb = sqrt_backward(A, dldz, cache, 'matrix-backprop', maxIter);
end
fprintf(' %fs Matrix-backprop\n', toc);
tic
for i = 1:runIter
dlda_svd = sqrt_backward(A, dldz, cache, 'lyap-svd', maxIter);
end
fprintf(' %fs Lyapunov SVD\n', toc);
tic
for i = 1:runIter
dlda_ns = sqrt_backward(sqrtA_ns, dldz, cache, 'lyap-ns', maxIter);
end
fprintf(' %fs Lyapunov Newton-Schulz (%i iter)\n', toc, maxIter);
if isGPU
sqrtA_svd = double(gather(sqrtA_svd));
sqrtA_db = double(gather(sqrtA_db));
sqrtA_ns = double(gather(sqrtA_ns));
dlda_mb = double(gather(dlda_mb));
dlda_svd = double(gather(dlda_svd));
dlda_ns = double(gather(dlda_ns));
end
matrix_norm = sqrt(sum(sum(A.*A)));
matrix_norm_sqrt = sqrt(sum(sum(sqrtA_true.*sqrtA_true)));
% Forward errors
fprintf('\nForward error:\n');
error_fwd = sqrt(sum(sum((sqrtA_svd - sqrtA_true).*(sqrtA_svd - sqrtA_true))))/matrix_norm_sqrt;
fprintf(' %d SVD\n', error_fwd);
error_fwd = sqrt(sum(sum((sqrtA_db - sqrtA_true).*(sqrtA_db - sqrtA_true))))/matrix_norm_sqrt;
fprintf(' %d Denman-Beavers (%i iter)\n', error_fwd, maxIter);
error_fwd = sqrt(sum(sum((sqrtA_ns - sqrtA_true).*(sqrtA_ns - sqrtA_true))))/matrix_norm_sqrt;
fprintf(' %d Newton-Schulz (%i iter)\n', error_fwd, maxIter);
% Backward errors
fprintf('\nBackward error:\n');
error_bkwd = sqrt(sum(sum((dlda_mb - dlda_true).*(dlda_mb - dlda_true))))/matrix_norm;
fprintf(' %d Matrix-backprop\n', error_bkwd);
error_bkwd = sqrt(sum(sum((dlda_svd - dlda_true).*(dlda_svd - dlda_true))))/matrix_norm;
fprintf(' %d Lyapunov SVD\n', error_bkwd);
error_bkwd = sqrt(sum(sum((dlda_ns - dlda_true).*(dlda_ns - dlda_true))))/matrix_norm;
fprintf(' %d Lyapunov Newton-Schulz (%i iter)\n', error_bkwd, maxIter);