% Example 6.6: Comparison of worst-case robust, Tikhonov, and nominal least squares % Section 6.4.2, Figure 6.16 % Boyd & Vandenberghe "Convex Optimization" % Original by Lieven Vandenberghe % Adapted for CVX Argyris Zymnis - 11/27/05 % (a figure is generated) % % Consider the least-squares problem: % minimize ||(A0 + u1*A1 + u2*A2)x - b||_2 % where u = [u1 u2]' is an uncertain parameter and ||u||_2 <= 1 % Three approximate solutions are found: % 1- nominal optimal (i.e. letting u=0) % 2- Tikhonov Regularized Solution: % minimize ||A0*x - b||_2 + delta*||x||_2 % for some delta (in this case we set delta = 0.1) % 3- worst-case robust approximation: % minimize sup{||u||_2 <= 1} ||(A0 + u1*A1 + u2*A2)x - b||_2) % (reduces to solving an SDP, see pages 323-324 in the book) 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); % we consider LS problems || (A0 + u1*A1 + u2*A2) x - b|| % where ||u|| leq rho % Nominal Solution xnom = A0\b; % Tikhonov Regularized Solution delta = .1; xtych = [A0; sqrt(delta)*eye(n)] \ [b; zeros(n,1)]; % Robust Least Squares solution cvx_begin sdp quiet 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 % Generate Random Trials notrials=100000; r = sqrt(rand(notrials,1)); % random on [0,1] with pdf g(r) = 2r; theta = 2*pi*rand(notrials,1); % uniform on [0,2pi] 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; % Plot histograms figure %subplot(211) [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