-
Notifications
You must be signed in to change notification settings - Fork 0
/
main.m
92 lines (85 loc) · 2.23 KB
/
main.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
clear;
N = [50 64 82 100 128 150 200 256];
% N = [100, 256, 512, 1024, 2048, 4096, 8142, 10000];
smooth = [1, 5, 10];
% smooth = [1];
lsm = length(smooth);
algs = cell(lsm+1, 1); algs{1} = 'APG';
for i = 1 : lsm
algs{i+1} = sprintf('MG-%d', smooth(i));
end
choice = [1:4];
infObj = 0;
tol = 1e-15;
verbose = 1;
linespec = ['o', '+', '*', '.', 'x'];
options = optimoptions('quadprog', 'Display', 'off', 'Algorithm', 'interior-point-convex', ...
'MaxIterations', 10, 'OptimalityTolerance', tol, 'StepTolerance', tol*0.01, 'LinearSolver', 'sparse');
for k = choice
n = N(k);
h = 1 / (n + 1);
Q0 = gallery('poisson', n) / h^2;
L0 = 8 / h^2;
% L0 = svds(Q0,1);
mu0 = svds(Q0, 1, 'smallest');
conv_fact = 1 - mu0 / L0;
phi = max(sin((1 : n)*3*pi/n), 0);
phi = vec(phi' * phi);
p0 = Q0 * phi;
x_ini = rand(n^2, 1);
% L = floor(2*log2(n)) - 1;
x = cell(1+lsm, 1);
hist = cell(1+lsm, 1);
disp('APG');
tic;
[x{1}, hist{1}] = apg(Q0, p0, L0, x_ini, tol, verbose);
toc;
infObj = min(infObj, min(hist{1}.F));
levels = [1];
llv = length(levels);
for i = 1 : lsm
s = smooth(i);
for j = 1 : llv
L = levels(j);
fprintf('\nMGProx-%d with level %d\n', s, L);
t0 = tic;
[x{i*j+1}, hist{i*j+1}] = mgproxL(Q0, p0, L0, x_ini, tol, L, s, verbose);
toc(t0);
infObj = min(infObj, min(hist{i*j+1}.F));
end
end
end
return;
figure(1);
for i = 1 : lsm + 1
plot((hist{i}.F-infObj)/hist{i}.F(1), linespec(i));
hold on;
end
hold off;
title('Function value');
xlabel('Iteration k');
ylabel('$(F_k-F_{\mbox{min}})/F_{\mbox{ini}}$', 'interpreter', 'latex');
legend(algs);
set(gca, 'YScale', 'log');
figure(2);
for i = 1 : lsm + 1
plot(hist{i}.G/hist{i}.G(1), linespec(i));
hold on;
end
hold off;
title('Norm of the proximal gradient mapping');
xlabel('Iteration k');
ylabel('$\Vert G_k \Vert/\Vert G_{\mbox{ini}} \Vert$', 'interpreter', 'latex');
legend(algs);
set(gca, 'YScale', 'log');
% tic;
% cvx_begin quiet
% variable y(N^2)
% minimize (0.5*quad_form(y,Q0) + p0'*y)
% subject to
% y >= 0;
% cvx_end
% toc
% function f = obj(q, p, x)
% f = 0.5*x'*q*x + p'*x;
% end