% 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: 136 variables, 64 equality constraints
------------------------------------------------------------
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
Split 8 free variables
eqs m = 64, order n = 129, dim = 161, blocks = 17
nnz(A) = 416 + 0, nnz(ADA) = 432, nnz(L) = 248
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            3.52E+01 0.000
  1 :   1.55E+00 1.00E+01 0.000 0.2843 0.9000 0.9000   3.63  1  1  1.3E+00
  2 :   4.23E+00 2.61E+00 0.000 0.2609 0.9000 0.9000   1.46  1  1  4.1E-01
  3 :   4.26E+00 5.75E-01 0.000 0.2203 0.9000 0.9000   1.16  1  1  1.2E-01
  4 :   4.23E+00 1.54E-01 0.000 0.2684 0.9000 0.9000   1.05  1  1  3.6E-02
  5 :   4.22E+00 3.78E-02 0.000 0.2447 0.9000 0.9000   1.02  1  1  9.3E-03
  6 :   4.21E+00 1.58E-03 0.389 0.0419 0.9900 0.6987   1.00  1  1  6.8E-04
  7 :   4.21E+00 3.57E-06 0.397 0.0023 0.9990 0.9990   1.00  1  1  2.7E-06
  8 :   4.21E+00 1.15E-06 0.206 0.3222 0.9000 0.9000   1.00  1  1  8.7E-07
  9 :   4.21E+00 1.96E-07 0.000 0.1700 0.9080 0.9000   1.00  2  1  2.4E-07
 10 :   4.21E+00 2.40E-10 0.000 0.0012 0.8470 0.9000   1.00  2  2  4.0E-08
 11 :   4.21E+00 1.73E-11 0.000 0.0724 0.9529 0.9900   1.00  2  2  2.9E-09

iter seconds digits       c*x               b*y
 11      0.1   Inf  4.2097054126e+00  4.2097054210e+00
|Ax-b| =   3.0e-10, [Ay-c]_+ =   6.7E-09, |x|=  9.4e+00, |y|=  6.0e+00

Detailed timing (sec)
   Pre          IPM          Post
2.000E-02    6.000E-02    0.000E+00    
Max-norms: ||b||=2, ||c|| = 4,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.43169.
------------------------------------------------------------
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.163 0.1084 0.9461 0.9450   1.00  2  2  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: 136 variables, 80 equality constraints
------------------------------------------------------------
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
Split 8 free variables
eqs m = 80, order n = 129, dim = 161, blocks = 17
nnz(A) = 688 + 0, nnz(ADA) = 1360, nnz(L) = 720
 it :     b*y       gap    delta  rate   t/tP*  t/tD*   feas cg cg  prec
  0 :            3.47E+01 0.000
  1 :   5.14E+00 9.84E+00 0.000 0.2839 0.9000 0.9000   3.58  1  1  9.4E-01
  2 :   5.06E+00 1.99E+00 0.000 0.2021 0.9000 0.9000   1.45  1  1  2.7E-01
  3 :   4.35E+00 4.82E-01 0.000 0.2423 0.9000 0.9000   1.16  1  1  9.0E-02
  4 :   4.23E+00 1.34E-01 0.000 0.2779 0.9000 0.9000   1.05  1  1  2.8E-02
  5 :   4.21E+00 3.49E-02 0.000 0.2606 0.9000 0.9000   1.02  1  1  7.9E-03
  6 :   4.21E+00 4.11E-03 0.167 0.1177 0.9000 0.0000   1.01  1  1  1.9E-03
  7 :   4.21E+00 1.13E-04 0.000 0.0276 0.9906 0.9900   1.00  1  1  1.3E-04
  8 :   4.21E+00 2.09E-06 0.000 0.0184 0.9900 0.9900   1.00  1  1  2.3E-06
  9 :   4.21E+00 5.44E-09 0.161 0.0026 0.9000 0.0000   1.00  2  1  9.6E-07
 10 :   4.21E+00 2.67E-10 0.187 0.0491 0.7780 0.9900   1.00  2  2  5.4E-08
 11 :   4.21E+00 3.33E-11 0.143 0.1248 0.9446 0.9450   1.00  2  2  6.8E-09

iter seconds digits       c*x               b*y
 11      0.1   Inf  4.2097055852e+00  4.2097056119e+00
|Ax-b| =   8.9e-10, [Ay-c]_+ =   1.4E-08, |x|=  1.0e+01, |y|=  4.4e+00

Detailed timing (sec)
   Pre          IPM          Post
1.000E-02    7.000E-02    0.000E+00    
Max-norms: ||b||=2, ||c|| = 4,
Cholesky |add|=0, |skip| = 0, ||L.L|| = 1.42968.
------------------------------------------------------------
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.6561   -0.6562   -0.6561
    0.8309    0.8310    0.8309