Figure 6.2: Penalty function approximation

% Section 6.1.2
% Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX Argyris Zymnis - 10/2005
%
% Comparison of the ell1, ell2, deadzone-linear and log-barrier
% penalty functions for the approximation problem:
%       minimize phi(A*x-b),
%
% where phi(x) is the penalty function
% Log-barrier will be implemented in the future version of CVX

% Generate input data
randn('state',0);
m=100; n=30;
A = randn(m,n);
b = randn(m,1);

% ell_1 approximation
% minimize   ||Ax+b||_1
disp('ell-one approximation');
cvx_begin
    variable x1(n)
    minimize(norm(A*x1+b,1))
cvx_end

% ell_2 approximation
% minimize ||Ax+b||_2
disp('ell-2');
x2=-A\b;

% deadzone penalty approximation
% minimize sum(deadzone(Ax+b,0.5))
% deadzone(y,z) = max(abs(y)-z,0)
dz = 0.5;
disp('deadzone penalty');
cvx_begin
    variable xdz(n)
    minimize(sum(max(abs(A*xdz+b)-dz,0)))
cvx_end


% log-barrier penalty approximation
%
% minimize -sum log(1-(ai'*x+bi)^2)

disp('log-barrier')

% parameters for Newton Method & line search
alpha=.01; beta=.5;

% minimize linfty norm to get starting point
cvx_begin
    variable xlb(n)
    minimize norm(A*xlb+b,Inf)
cvx_end
linf = cvx_optval;
A = A/(1.1*linf);
b = b/(1.1*linf);

for iters = 1:50

   yp = 1 - (A*xlb+b);  ym = (A*xlb+b) + 1;
   f = -sum(log(yp)) - sum(log(ym));
   g = A'*(1./yp) - A'*(1./ym);
   H = A'*diag(1./(yp.^2) + 1./(ym.^2))*A;
   v = -H\g;
   fprime = g'*v;
   ntdecr = sqrt(-fprime);
   if (ntdecr < 1e-5), break; end;

   t = 1;
   newx = xlb + t*v;
   while ((min(1-(A*newx +b)) < 0) | (min((A*newx +b)+1) < 0))
       t = beta*t;
       newx = xlb + t*v;
   end;
   newf = -sum(log(1 - (A*newx+b))) - sum(log(1+(A*newx+b)));
   while (newf > f + alpha*t*fprime)
       t = beta*t;
       newx = xlb + t*v;
       newf = -sum(log(1-(A*newx+b))) - sum(log(1+(A*newx+b)));
   end;
   xlb = xlb+t*v;
end


% Plot histogram of residuals

ss = max(abs([A*x1+b; A*x2+b; A*xdz+b;  A*xlb+b]));
tt = -ceil(ss):0.05:ceil(ss);  % sets center for each bin
[N1,hist1] = hist(A*x1+b,tt);
[N2,hist2] = hist(A*x2+b,tt);
[N3,hist3] = hist(A*xdz+b,tt);
[N4,hist4] = hist(A*xlb+b,tt);


range_max=2.0;  rr=-range_max:1e-2:range_max;

figure(1), clf, hold off
subplot(4,1,1),
bar(hist1,N1);
hold on
plot(rr, abs(rr)*40/3, '-');
ylabel('p=1')
axis([-range_max range_max 0 40]);
hold off

subplot(4,1,2),
bar(hist2,N2);
hold on;
plot(rr,2*rr.^2),
ylabel('p=2')
axis([-range_max range_max 0 11]);
hold off

subplot(4,1,3),
bar(hist3,N3);
hold on
plot(rr,30/3*max(0,abs(rr)-dz))
ylabel('Deadzone')
axis([-range_max range_max 0 25]);
hold off

subplot(4,1,4),
bar(hist4,N4);
rr_lb=linspace(-1+(1e-6),1-(1e-6),600);
hold on
plot(rr_lb, -3*log(1-rr_lb.^2),rr,2*rr.^2,'--')
axis([-range_max range_max 0 11]);
ylabel('Log barrier'),
xlabel('r')
hold off
ell-one approximation
 
Calling SDPT3: 230 variables, 100 equality constraints
------------------------------------------------------------

 num. of constraints = 100
 dim. of socp   var  = 200,   num. of socp blk  = 100
 dim. of free   var  = 30 *** convert ublk to lblk
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
    NT      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      mean(obj)   cputime
-------------------------------------------------------------------
 0|0.000|0.000|9.0e-01|6.3e+01|3.2e+05| 1.045672e+02| 0:0:00| chol  1  1 
 1|1.000|0.882|1.0e-05|7.5e+00|2.0e+04| 2.622372e+03| 0:0:00| chol  1  1 
 2|1.000|0.990|4.3e-06|9.5e-02|1.6e+03| 8.016733e+02| 0:0:00| chol  1  1 
 3|0.953|1.000|2.3e-06|3.0e-03|7.5e+01| 6.592898e+01| 0:0:00| chol  1  1 
 4|0.784|0.271|3.6e-05|2.3e-03|3.3e+01| 5.175219e+01| 0:0:00| chol  1  1 
 5|0.920|0.209|1.8e-05|1.8e-03|2.1e+01| 4.961556e+01| 0:0:00| chol  1  1 
 6|0.928|0.439|8.9e-06|1.0e-03|1.1e+01| 5.162648e+01| 0:0:00| chol  1  1 
 7|0.975|0.470|4.8e-06|5.4e-04|5.4e+00| 5.300594e+01| 0:0:00| chol  1  1 
 8|1.000|0.328|3.5e-07|3.6e-04|3.6e+00| 5.365787e+01| 0:0:00| chol  1  1 
 9|1.000|0.474|5.0e-08|1.9e-04|1.9e+00| 5.431941e+01| 0:0:00| chol  1  1 
10|1.000|0.617|1.8e-08|7.3e-05|7.0e-01| 5.480849e+01| 0:0:00| chol  1  1 
11|0.712|0.265|1.9e-08|5.4e-05|5.1e-01| 5.489203e+01| 0:0:01| chol  1  1 
12|1.000|0.216|1.2e-08|4.2e-05|4.1e-01| 5.494415e+01| 0:0:01| chol  1  1 
13|1.000|0.448|6.3e-09|2.3e-05|2.2e-01| 5.502483e+01| 0:0:01| chol  1  1 
14|1.000|0.780|1.4e-09|5.1e-06|4.9e-02| 5.510542e+01| 0:0:01| chol  1  1 
15|0.986|0.975|1.5e-10|1.3e-07|1.2e-03| 5.512833e+01| 0:0:01| chol  1  1 
16|0.989|0.989|1.4e-11|1.6e-06|2.3e-05| 5.512891e+01| 0:0:01| chol  1  1 
17|0.973|0.945|3.9e-13|3.1e-08|1.0e-06| 5.512892e+01| 0:0:01| chol  1  1 
18|0.569|0.944|1.7e-13|1.4e-09|1.0e-07| 5.512892e+01| 0:0:01|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 18
 primal objective value =  5.51289216e+01
 dual   objective value =  5.51289215e+01
 gap := trace(XZ)       = 1.00e-07
 relative gap           = 9.00e-10
 actual relative gap    = 8.41e-10
 rel. primal infeas     = 1.70e-13
 rel. dual   infeas     = 1.39e-09
 norm(X), norm(y), norm(Z) = 1.2e+01, 8.8e+00, 1.3e+01
 norm(A), norm(b), norm(C) = 7.9e+01, 1.0e+01, 1.1e+01
 Total CPU time (secs)  = 0.8  
 CPU time per iteration = 0.0  
 termination code       =  0
 DIMACS: 5.8e-13  0.0e+00  7.6e-09  0.0e+00  8.4e-10  9.0e-10
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +55.1289
ell-2
deadzone penalty
 
Calling SDPT3: 300 variables, 130 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 130
 dim. of socp   var  = 200,   num. of socp blk  = 100
 dim. of linear var  = 100
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
    NT      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      mean(obj)   cputime
-------------------------------------------------------------------
 0|0.000|0.000|9.5e+00|1.4e+01|1.2e+04| 3.535534e+01| 0:0:00| chol  1  1 
 1|1.000|1.000|2.9e-06|1.0e-01|8.6e+02|-4.193082e+02| 0:0:00| chol  1  1 
 2|1.000|0.856|1.0e-07|2.3e-02|1.3e+02|-6.131583e+01| 0:0:00| chol  1  1 
 3|0.783|0.818|3.8e-08|5.0e-03|6.9e+01|-3.789956e+01| 0:0:00| chol  1  1 
 4|0.828|0.790|2.6e-08|1.1e-03|3.2e+01|-2.827240e+01| 0:0:00| chol  1  1 
 5|0.796|0.872|5.3e-09|1.5e-04|1.2e+01|-2.303367e+01| 0:0:00| chol  1  1 
 6|0.902|0.973|5.6e-10|5.1e-06|3.9e+00|-2.164155e+01| 0:0:00| chol  1  1 
 7|0.782|0.804|5.9e-10|1.1e-06|1.6e+00|-2.150812e+01| 0:0:00| chol  1  1 
 8|0.898|0.870|2.5e-10|1.5e-07|3.6e-01|-2.148161e+01| 0:0:00| chol  1  1 
 9|0.738|0.860|4.2e-11|2.2e-08|1.1e-01|-2.146101e+01| 0:0:00| chol  1  1 
10|0.950|0.904|2.1e-12|2.2e-09|8.3e-03|-2.146885e+01| 0:0:00| chol  1  1 
11|0.932|0.980|1.4e-13|5.3e-11|3.8e-04|-2.146814e+01| 0:0:00| chol  1  1 
12|0.988|0.988|2.5e-14|1.6e-12|4.4e-06|-2.146821e+01| 0:0:00| chol  1  1 
13|0.993|0.996|3.0e-14|1.0e-12|5.8e-08|-2.146821e+01| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 13
 primal objective value = -2.14682118e+01
 dual   objective value = -2.14682118e+01
 gap := trace(XZ)       = 5.77e-08
 relative gap           = 1.31e-09
 actual relative gap    = 1.31e-09
 rel. primal infeas     = 3.05e-14
 rel. dual   infeas     = 1.01e-12
 norm(X), norm(y), norm(Z) = 1.1e+01, 4.4e+00, 1.2e+01
 norm(A), norm(b), norm(C) = 5.8e+01, 1.1e+01, 1.1e+01
 Total CPU time (secs)  = 0.5  
 CPU time per iteration = 0.0  
 termination code       =  0
 DIMACS: 1.7e-13  0.0e+00  3.9e-12  0.0e+00  1.3e-09  1.3e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +21.4682
log-barrier
 
Calling SDPT3: 200 variables, 31 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 31
 dim. of socp   var  = 200,   num. of socp blk  = 100
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
    NT      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      mean(obj)   cputime
-------------------------------------------------------------------
 0|0.000|0.000|7.0e+01|1.2e+01|1.8e+03| 0.000000e+00| 0:0:00| chol  1  1 
 1|0.981|1.000|1.3e+00|9.8e-02|4.5e+01|-6.232429e+00| 0:0:00| chol  1  1 
 2|1.000|1.000|9.6e-08|9.8e-03|6.5e+00|-3.345428e+00| 0:0:00| chol  1  1 
 3|1.000|0.799|1.8e-08|2.8e-03|1.5e+00|-1.031687e+00| 0:0:00| chol  1  1 
 4|0.683|0.777|7.6e-09|6.9e-04|6.8e-01|-1.181387e+00| 0:0:00| chol  1  1 
 5|0.900|0.845|8.1e-10|1.2e-04|1.9e-01|-1.190459e+00| 0:0:00| chol  1  1 
 6|0.725|0.701|3.6e-10|3.5e-05|7.7e-02|-1.192951e+00| 0:0:00| chol  1  1 
 7|0.826|0.879|6.6e-11|4.3e-06|2.1e-02|-1.196197e+00| 0:0:00| chol  1  1 
 8|0.618|0.956|2.4e-11|2.0e-07|7.6e-03|-1.198192e+00| 0:0:00| chol  1  1 
 9|0.815|1.000|2.8e-11|9.9e-10|2.3e-03|-1.200331e+00| 0:0:00| chol  1  1 
10|0.930|0.971|1.9e-12|1.3e-10|1.6e-04|-1.201199e+00| 0:0:00| chol  1  1 
11|0.969|0.812|6.0e-14|2.5e-11|8.5e-06|-1.201268e+00| 0:0:00| chol  1  1 
12|0.911|1.000|5.8e-15|1.0e-12|1.1e-06|-1.201270e+00| 0:0:00| chol  1  1 
13|0.597|1.000|8.7e-15|1.0e-12|5.8e-07|-1.201270e+00| 0:0:00| chol  1  1 
14|0.602|1.000|5.4e-15|1.0e-12|3.1e-07|-1.201270e+00| 0:0:00| chol  1  1 
15|0.603|1.000|2.5e-15|1.0e-12|1.7e-07|-1.201270e+00| 0:0:00| chol  1  1 
16|0.604|1.000|1.2e-15|1.0e-12|8.8e-08|-1.201270e+00| 0:0:00| chol  1  1 
17|0.605|1.000|3.6e-15|1.0e-12|4.6e-08|-1.201270e+00| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 17
 primal objective value = -1.20127042e+00
 dual   objective value = -1.20127047e+00
 gap := trace(XZ)       = 4.64e-08
 relative gap           = 1.36e-08
 actual relative gap    = 1.36e-08
 rel. primal infeas     = 3.62e-15
 rel. dual   infeas     = 1.00e-12
 norm(X), norm(y), norm(Z) = 3.1e-01, 1.5e+00, 1.5e+01
 norm(A), norm(b), norm(C) = 5.7e+01, 2.0e+00, 1.0e+01
 Total CPU time (secs)  = 0.4  
 CPU time per iteration = 0.0  
 termination code       =  0
 DIMACS: 3.6e-15  0.0e+00  3.4e-12  0.0e+00  1.4e-08  1.4e-08
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.20127