Example 6.6: Comparison of worst-case robust, Tikhonov, and nominal least squares
clear
cvx_quiet(false);
m = 50;
n = 20;
randn('state',0);
rand('state',0);
A0 = randn(m,n);
[U,S,V] = svd(A0);
S= diag(fliplr(logspace(-0.7,1,n)));
A0 = U(:,1:n)*S*V';
A1 = randn(m,n); A1 = A1/norm(A1);
A2 = randn(m,n); A2 = A2/norm(A2);
Aperb0 = [A1;A2];
p = 2;
b = U(:,1:n)*randn(n,1) + .1*randn(m,1);
xnom = A0\b;
delta = .1;
xtych = [A0; sqrt(delta)*eye(n)] \ [b; zeros(n,1)];
cvx_begin sdp
variables t lambda xrob(n)
minimize(t+lambda)
subject to
[eye(m) A1*xrob A2*xrob A0*xrob-b; ...
[A1*xrob A2*xrob]' lambda*eye(2) zeros(2,1); ...
[A0*xrob-b]' zeros(1,2) t] >= 0;
cvx_end
notrials=100000;
r = sqrt(rand(notrials,1));
theta = 2*pi*rand(notrials,1);
v = [r.*cos(theta) r.*sin(theta)];
ls_res = zeros(1,notrials);
rob2_res = zeros(1,notrials);
rob_res = zeros(1,notrials);
tych_res = zeros(1,notrials);
for i =1:notrials
A = A0 + v(i,1)*A1 + v(i,2)*A2;
ls_res(i) = norm(A*xnom-b);
rob_res(i) = norm(A*xrob-b);
tych_res(i) = norm(A*xtych-b);
end;
figure
[N1, hist1] = hist(ls_res,[min(ls_res):.1:max(ls_res)]);
freq1 = N1/notrials;
[N2, hist2] = hist(rob_res,hist1);
freq2 = N2/notrials;
[N3, hist3] = hist(tych_res,hist1);
freq3 = N3/notrials;
h = bar(hist3,freq3);
text(3, 0.07, 'Tikhonov');
set(h,'FaceColor',0.90*[1 1 1]);
hold on
h = bar(hist2,freq2);
text(4.2, 0.05, 'Nominal');
set(h,'FaceColor',0.80*[1 1 1]);
h = bar(hist2,freq2);
set(h,'FaceColor','none');
text(2.6, 0.2, 'Robust LS');
h = bar(hist3,freq3);
set(h,'FaceColor','none');
h = bar(hist1,freq1);
set(h,'FaceColor','none');
xlabel('||(A0 + u1*A1 + u2*A2)*x - b||_2')
ylabel('Frequency')
hold off
Calling SDPT3: 1431 variables, 22 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------
num. of constraints = 22
dim. of sdp var = 53, num. of sdp blk = 1
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version predcorr gam expon scale_data
HKM 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap mean(obj) cputime
-------------------------------------------------------------------
0|0.000|0.000|4.9e+01|6.7e+00|3.0e+04| 1.325000e+03| 0:0:00| chol 1 1
1|0.958|1.000|2.0e+00|6.9e-02|2.4e+03| 1.030273e+03| 0:0:00| chol 1 1
2|0.976|1.000|4.9e-02|6.9e-03|7.6e+01| 2.153922e+01| 0:0:00| chol 1 1
3|1.000|1.000|2.7e-07|1.1e-02|2.1e+01| 2.747313e+00| 0:0:00| chol 1 1
4|0.965|0.877|1.7e-08|1.4e-03|8.9e-01|-6.388350e+00| 0:0:00| chol 1 1
5|0.994|0.935|3.0e-09|9.5e-05|1.5e-01|-6.689767e+00| 0:0:00| chol 1 1
6|0.965|0.983|1.8e-09|2.3e-06|5.8e-03|-6.756737e+00| 0:0:00| chol 1 1
7|0.934|1.000|1.1e-09|6.9e-08|4.0e-04|-6.759313e+00| 0:0:01| chol 1 1
8|0.858|1.000|6.5e-09|2.3e-10|1.2e-04|-6.759442e+00| 0:0:01| chol 1 1
9|0.970|0.986|1.9e-09|3.4e-10|4.4e-06|-6.759498e+00| 0:0:01| chol 2 2
10|1.000|1.000|3.4e-12|3.9e-10|5.8e-07|-6.759499e+00| 0:0:01| chol 1 2
11|1.000|1.000|1.4e-12|1.0e-12|1.6e-08|-6.759500e+00| 0:0:01|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations = 11
primal objective value = -6.75949959e+00
dual objective value = -6.75949960e+00
gap := trace(XZ) = 1.63e-08
relative gap = 1.13e-09
actual relative gap = 1.12e-09
rel. primal infeas = 1.36e-12
rel. dual infeas = 1.00e-12
norm(X), norm(y), norm(Z) = 6.3e+00, 5.6e+00, 9.9e+00
norm(A), norm(b), norm(C) = 2.6e+01, 2.4e+00, 1.1e+01
Total CPU time (secs) = 0.8
CPU time per iteration = 0.1
termination code = 0
DIMACS: 1.6e-12 0.0e+00 4.3e-12 0.0e+00 1.1e-09 1.1e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +6.7595