Contents
function h = control
Problem data
s = RandStream.create('mt19937ar','seed',0);
RandStream.setDefaultStream(s);
n = 50;
m = 5;
A = randn(m,4*n,8);
b = randn(m,1);
f = @(u,A) max(A*u(:) + b);
CVX
tic
cvx_begin quiet
variables y_aaa(n,4) y_aab(n,4) y_aba(n,4) y_abb(n,4) ...
y_baa(n,4) y_bab(n,4) y_bba(n,4) y_bbb(n,4);
obj = f(y_aaa,A(:,:,1)) + f(y_aab,A(:,:,2)) + f(y_aba,A(:,:,3)) + f(y_abb,A(:,:,4)) + ...
f(y_baa,A(:,:,5)) + f(y_bab,A(:,:,6)) + f(y_bba,A(:,:,7)) + f(y_bbb,A(:,:,8));
minimize(obj);
subject to
y_aaa(:,1) == y_aab(:,1);
y_aab(:,1) == y_aba(:,1);
y_aba(:,1) == y_abb(:,1);
y_abb(:,1) == y_baa(:,1);
y_baa(:,1) == y_bab(:,1);
y_bab(:,1) == y_bba(:,1);
y_bba(:,1) == y_bbb(:,1);
y_aaa(:,2) == y_aab(:,2);
y_aab(:,2) == y_aba(:,2);
y_aba(:,2) == y_abb(:,2);
y_baa(:,2) == y_bab(:,2);
y_bab(:,2) == y_bba(:,2);
y_bba(:,2) == y_bbb(:,2);
y_aaa(:,3) == y_aab(:,3);
y_aba(:,3) == y_abb(:,3);
y_baa(:,3) == y_bab(:,3);
y_bba(:,3) == y_bbb(:,3);
y_aaa <= 1;
y_aab <= 1;
y_aba <= 1;
y_abb <= 1;
y_baa <= 1;
y_bab <= 1;
y_bba <= 1;
y_bbb <= 1;
y_aaa >= -1;
y_aab >= -1;
y_aba >= -1;
y_abb >= -1;
y_baa >= -1;
y_bab >= -1;
y_bba >= -1;
y_bbb >= -1;
cvx_end
toc
h.p_cvx = cvx_optval/8;
cvx_begin quiet
variables w_aaa(n,4) w_aab(n,4) w_aba(n,4) w_abb(n,4) ...
w_baa(n,4) w_bab(n,4) w_bba(n,4) w_bbb(n,4);
obj = f(w_aaa,A(:,:,1)) + f(w_aab,A(:,:,2)) + f(w_aba,A(:,:,3)) + f(w_abb,A(:,:,4)) + ...
f(w_baa,A(:,:,5)) + f(w_bab,A(:,:,6)) + f(w_bba,A(:,:,7)) + f(w_bbb,A(:,:,8));
minimize(obj);
subject to
w_aaa <= 1;
w_aab <= 1;
w_aba <= 1;
w_abb <= 1;
w_baa <= 1;
w_bab <= 1;
w_bba <= 1;
w_bbb <= 1;
w_aaa >= -1;
w_aab >= -1;
w_aba >= -1;
w_abb >= -1;
w_baa >= -1;
w_bab >= -1;
w_bba >= -1;
w_bbb >= -1;
cvx_end
h.p_cvx_pres = cvx_optval/8;
Elapsed time is 0.291739 seconds.
ADMM
MAX_ITER = 100;
ABSTOL = 1e-4;
RELTOL = 1e-2;
lambda = 1;
x_aaa = zeros(n,4); x_aab = zeros(n,4); x_aba = zeros(n,4); x_abb = zeros(n,4);
x_baa = zeros(n,4); x_bab = zeros(n,4); x_bba = zeros(n,4); x_bbb = zeros(n,4);
z_0 = zeros(n,1);
z_a = zeros(n,1); z_b = zeros(n,1);
z_aa = zeros(n,1); z_ab = zeros(n,1); z_ba = zeros(n,1); z_bb = zeros(n,1);
z_aaa = zeros(n,4); z_aab = zeros(n,4); z_aba = zeros(n,4); z_abb = zeros(n,4);
z_baa = zeros(n,4); z_bab = zeros(n,4); z_bba = zeros(n,4); z_bbb = zeros(n,4);
z = [ z_aaa z_aab z_aba z_abb z_baa z_bab z_bba z_bbb ];
u_aaa = zeros(n,4); u_aab = zeros(n,4); u_aba = zeros(n,4); u_abb = zeros(n,4);
u_baa = zeros(n,4); u_bab = zeros(n,4); u_bba = zeros(n,4); u_bbb = zeros(n,4);
tstart = tic;
tx = 0;
for k = 1:MAX_ITER
titer = tic;
x_aaa = prox_f(z_aaa - u_aaa, lambda, A(:,:,1), b);
x_aab = prox_f(z_aab - u_aab, lambda, A(:,:,2), b);
x_aba = prox_f(z_aba - u_aba, lambda, A(:,:,3), b);
x_abb = prox_f(z_abb - u_abb, lambda, A(:,:,4), b);
x_baa = prox_f(z_baa - u_baa, lambda, A(:,:,5), b);
x_bab = prox_f(z_bab - u_bab, lambda, A(:,:,6), b);
x_bba = prox_f(z_bba - u_bba, lambda, A(:,:,7), b);
x_bbb = prox_f(z_bbb - u_bbb, lambda, A(:,:,8), b);
tx = tx + toc(titer);
z_0 = avg(x_aaa(:,1), x_aab(:,1), x_aba(:,1), x_abb(:,1), ...
x_baa(:,1), x_bab(:,1), x_bba(:,1), x_bbb(:,1));
z_a = avg(x_aaa(:,2), x_aab(:,2), x_aba(:,2), x_abb(:,2));
z_b = avg(x_baa(:,2), x_bab(:,2), x_bba(:,2), x_bbb(:,2));
z_aa = avg(x_aaa(:,3), x_aab(:,3));
z_ab = avg(x_aba(:,3), x_abb(:,3));
z_ba = avg(x_baa(:,3), x_bab(:,3));
z_bb = avg(x_bba(:,3), x_bbb(:,3));
z_aaa = [z_0 z_a z_aa x_aaa(:,4)];
z_aab = [z_0 z_a z_aa x_aab(:,4)];
z_aba = [z_0 z_a z_ab x_aba(:,4)];
z_abb = [z_0 z_a z_ab x_abb(:,4)];
z_baa = [z_0 z_b z_ba x_baa(:,4)];
z_bab = [z_0 z_b z_ba x_bab(:,4)];
z_bba = [z_0 z_b z_bb x_bba(:,4)];
z_bbb = [z_0 z_b z_bb x_bbb(:,4)];
u_aaa = u_aaa + x_aaa - z_aaa;
u_aab = u_aab + x_aab - z_aab;
u_aba = u_aba + x_aba - z_aba;
u_abb = u_abb + x_abb - z_abb;
u_baa = u_baa + x_baa - z_baa;
u_bab = u_bab + x_bab - z_bab;
u_bba = u_bba + x_bba - z_bba;
u_bbb = u_bbb + x_bbb - z_bbb;
x = [ x_aaa x_aab x_aba x_abb x_baa x_bab x_bba x_bbb ];
zold = z;
z = [ z_aaa z_aab z_aba z_abb z_baa z_bab z_bba z_bbb ];
u = [ u_aaa u_aab u_aba u_abb u_baa u_bab u_bba u_bbb ];
h.r_norm(k) = norm(x - z,'fro');
h.s_norm(k) = norm(-(z - zold)/lambda,'fro');
h.eps_pri(k) = sqrt(numel(x))*ABSTOL + RELTOL*max(norm(x,'fro'), norm(-z,'fro'));
h.eps_dual(k) = sqrt(numel(x))*ABSTOL + RELTOL*norm(u/lambda,'fro');
h.objval(k) = (1/8)*(f(z_aaa,A(:,:,1)) + f(z_aab,A(:,:,2)) + f(z_aba,A(:,:,3)) + ...
f(z_abb,A(:,:,4)) + f(z_baa,A(:,:,5)) + f(z_bab,A(:,:,6)) + ...
f(z_bba,A(:,:,7)) + f(z_bbb,A(:,:,8)));
if mod(k,10) == 0
fprintf('%3d\t%10.4f\t%10.4f\t%10.4f\t%10.4f\n', k, ...
h.r_norm(k), h.eps_pri(k), h.s_norm(k), h.eps_dual(k));
end
if h.r_norm(k) < h.eps_pri(k) && h.s_norm(k) < h.eps_dual(k)
break;
end
end
tx/k
toc(tstart)
10 0.3626 0.3495 1.3970 0.1377
20 0.1103 0.3718 0.5046 0.1376
30 0.0878 0.3801 0.2681 0.1377
40 0.0399 0.3839 0.1831 0.1377
ans =
0.0058
Elapsed time is 0.299336 seconds.
Plotting
figure(1);
hold on;
plot(h.objval);
plot(repmat(h.p_cvx,k,1), 'k--');
plot(repmat(h.p_cvx_pres,k,1), 'k--');
hold off;
xlim([0 25]);
ylim([-80 -10]);
print -depsc ctrl_progress.eps;
end
function x = avg(varargin)
N = length(varargin);
x = 0;
for k = 1:N
x = x + varargin{k};
end
x = x/N;
end
function x = prox_f(v, lambda, A, b)
[n T] = size(v);
params.A_0 = A(:,1:n);
params.A_1 = A(:,n+1:2*n);
params.A_2 = A(:,2*n+1:3*n);
params.A_3 = A(:,3*n+1:4*n);
params.lambda = lambda;
params.b = b;
params.v_0 = v(:,1);
params.v_1 = v(:,2);
params.v_2 = v(:,3);
params.v_3 = v(:,4);
settings.verbose = 0;
vars = csolve(params, settings);
x = [vars.x_0 vars.x_1 vars.x_2 vars.x_3];
end
function x = prox_f_cvx(v, lambda, f, A)
[n T] = size(v);
cvx_begin quiet
variable x(n,T)
minimize(f(x,A) + (1/(2*lambda))*square_pos(norm(x - v,'fro')))
subject to
x <= 1;
x >= -1;
cvx_end
end