Exercise 4.5: Show the equivalence of 3 convex problem formations

% From Boyd & Vandenberghe, "Convex Optimization"
% Joëlle Skaf - 08/17/05
%
% Shows the equivalence of the following 3 problems:
% 1) Robust least-squares problem
%           minimize    sum_{i=1}^{m} phi(a_i'*x - bi)
%    where phi(u) = u^2             for |u| <= M
%                   M(2|u| - M)     for |u| >  M
% 2) Least-squares with variable weights
%           minimize    sum_{i=1}^{m} (a_i'*x - bi)^2/(w_i+1) + M^2*1'*w
%               s.t.    w >= 0
% 3) Quadratic program
%           minimize    sum_{i=1}^{m} (u_i^2 + 2*M*v_i)
%               s.t.    -u - v <= Ax - b <= u + v
%                       0 <= u <= M*1
%                       v >= 0

% Generate input data
randn('state',0);
m = 16; n = 8;
A = randn(m,n);
b = randn(m,1);
M = 2;

% (a) robust least-squares problem
disp('Computing the solution of the robust least-squares problem...');
cvx_begin
    variable x1(n)
    minimize( sum(huber(A*x1-b,M)) )
cvx_end

% (b)least-squares problem with variable weights
disp('Computing the solution of the least-squares problem with variable weights...');
cvx_begin
    variable x2(n)
    variable w(m)
    minimize( sum(quad_over_lin(diag(A*x2-b),w'+1)) + M^2*ones(1,m)*w)
    w >= 0;
cvx_end

% (c) quadratic program
disp('Computing the solution of the quadratic program...');
cvx_begin
    variable x3(n)
    variable u(m)
    variable v(m)
    minimize( sum(square(u) +  2*M*v) )
    A*x3 - b <= u + v;
    A*x3 - b >= -u - v;
    u >= 0;
    u <= M;
    v >= 0;
cvx_end

% Display results
disp('------------------------------------------------------------------------');
disp('The optimal solutions for problem formulations 1, 2 and 3 are given');
disp('respectively as follows (per column): ');
[x1 x2 x3]
Computing the solution of the robust least-squares problem...
 
Calling sedumi: 112 variables, 56 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 56, order n = 97, dim = 129, blocks = 33
nnz(A) = 256 + 0, nnz(ADA) = 976, nnz(L) = 516
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            1.92E+01 0.000
  1 :   1.26E+01 4.67E+00 0.000 0.2436 0.9000 0.9000   1.65  1  1  1.1E+00
  2 :  -1.66E+00 1.05E+00 0.000 0.2243 0.9000 0.9000   1.26  1  1  1.4E+00
  3 :  -3.68E+00 2.21E-01 0.000 0.2109 0.9000 0.9000   1.07  1  1  1.6E-01
  4 :  -4.08E+00 5.49E-02 0.000 0.2483 0.9000 0.9000   1.02  1  1  3.7E-02
  5 :  -4.18E+00 1.30E-02 0.000 0.2364 0.9000 0.9000   1.01  1  1  8.5E-03
  6 :  -4.21E+00 8.27E-04 0.000 0.0637 0.9900 0.9900   1.00  1  1  5.4E-04
  7 :  -4.21E+00 1.04E-08 0.376 0.0000 0.9977 0.9990   1.00  1  1  9.1E-07
  8 :  -4.21E+00 2.08E-09 0.076 0.2012 0.9000 0.9000   1.00  1  1  1.8E-07
  9 :  -4.21E+00 2.56E-10 0.173 0.1229 0.9450 0.9450   1.00  1  1  2.2E-08
 10 :  -4.21E+00 2.41E-11 0.356 0.0942 0.9900 0.9900   1.00  2  2  2.1E-09

iter seconds digits       c*x               b*y
 10      0.1   Inf -4.2097051512e+00 -4.2097051484e+00
|Ax-b| =   1.3e-09, [Ay-c]_+ =   1.8E-09, |x|=  1.5e+01, |y|=  3.3e+00

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    6.000E-02    0.000E+00    
Max-norms: ||b||=1, ||c|| = 2,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 2.94583.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +4.20971
Computing the solution of the least-squares problem with variable weights...
 
Calling sedumi: 304 variables, 40 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 40, order n = 49, dim = 305, blocks = 17
nnz(A) = 192 + 0, nnz(ADA) = 640, nnz(L) = 340
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            4.38E-01 0.000
  1 :  -8.64E+00 1.00E-01 0.000 0.2287 0.9000 0.9000   2.39  1  1  6.8E-01
  2 :  -2.01E+01 3.89E-03 0.000 0.0388 0.9900 0.9900   1.18  1  1  2.4E-02
  3 :  -2.02E+01 2.92E-06 0.000 0.0008 0.9999 0.9999   1.01  1  1  1.8E-05
  4 :  -2.02E+01 5.86E-07 0.042 0.2011 0.9000 0.9000   1.00  1  1  3.7E-06
  5 :  -2.02E+01 1.12E-07 0.000 0.1905 0.9000 0.9000   1.00  1  1  7.0E-07
  6 :  -2.02E+01 1.91E-09 0.000 0.0171 0.8869 0.9000   1.00  1  1  1.2E-07
  7 :  -2.02E+01 2.07E-10 0.162 0.1086 0.9461 0.9450   1.00  1  1  1.3E-08

iter seconds digits       c*x               b*y
  7      0.1   Inf -2.0209705011e+01 -2.0209704929e+01
|Ax-b| =   3.3e-08, [Ay-c]_+ =   8.0E-09, |x|=  1.7e+01, |y|=  2.9e+00

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    5.000E-02    0.000E+00    
Max-norms: ||b||=3, ||c|| = 1.488490e+00,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.6581.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +4.2097
Computing the solution of the quadratic program...
 
Calling sedumi: 128 variables, 56 equality constraints
   For improved efficiency, sedumi is solving the dual problem.
------------------------------------------------------------
SeDuMi 1.21 by AdvOL, 2005-2008 and Jos F. Sturm, 1998-2003.
Alg = 2: xz-corrector, Adaptive Step-Differentiation, theta = 0.250, beta = 0.500
eqs m = 56, order n = 113, dim = 145, blocks = 17
nnz(A) = 400 + 0, nnz(ADA) = 688, nnz(L) = 372
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            4.06E+01 0.000
  1 :   1.88E+00 1.16E+01 0.000 0.2860 0.9000 0.9000   3.61  1  1  1.5E+00
  2 :  -4.11E+00 2.33E+00 0.000 0.2002 0.9000 0.9000   1.45  1  1  5.8E-01
  3 :  -4.15E+00 5.60E-01 0.000 0.2407 0.9000 0.9000   1.16  1  1  1.2E-01
  4 :  -4.19E+00 1.55E-01 0.000 0.2774 0.9000 0.9000   1.05  1  1  3.4E-02
  5 :  -4.20E+00 4.02E-02 0.000 0.2590 0.9000 0.9000   1.02  1  1  9.0E-03
  6 :  -4.21E+00 5.47E-03 0.000 0.1360 0.9046 0.9000   1.01  1  1  1.4E-03
  7 :  -4.21E+00 4.12E-05 0.000 0.0075 0.9901 0.9900   1.00  1  1  1.6E-05
  8 :  -4.21E+00 8.54E-07 0.167 0.0207 0.9900 0.9900   1.00  1  1  3.3E-07
  9 :  -4.21E+00 1.68E-08 0.000 0.0197 0.9900 0.9900   1.00  1  1  7.4E-09

iter seconds digits       c*x               b*y
  9      0.0   7.9 -4.2097051377e+00 -4.2097051852e+00
|Ax-b| =   1.4e-08, [Ay-c]_+ =   1.4E-09, |x|=  1.4e+01, |y|=  2.9e+00

Detailed timing (sec)
   Pre          IPM          Post
0.000E+00    4.000E-02    0.000E+00    
Max-norms: ||b||=4, ||c|| = 2,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 2.20232.
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +4.20971
------------------------------------------------------------------------
The optimal solutions for problem formulations 1, 2 and 3 are given
respectively as follows (per column): 

ans =

    0.3888    0.3888    0.3888
    0.1262    0.1262    0.1262
   -0.3337   -0.3337   -0.3337
    0.1326    0.1326    0.1326
    0.5500    0.5500    0.5500
    0.3526    0.3526    0.3526
   -0.6562   -0.6562   -0.6562
    0.8309    0.8310    0.8309