% 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
%           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)
w >= 0;
cvx_end

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 SDPT3 4.0: 152 variables, 64 equality constraints
------------------------------------------------------------

num. of constraints = 64
dim. of sdp    var  = 32,   num. of sdp  blk  = 16
dim. of socp   var  = 32,   num. of socp blk  = 16
dim. of linear var  = 64
dim. of free   var  =  8 *** convert ublk to lblk
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version  predcorr  gam  expon  scale_data
HKM      1      0.000   1        0
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
0|0.000|0.000|7.7e+00|9.7e+00|4.5e+04| 1.752634e+03  0.000000e+00| 0:0:00| chol  1  1
1|1.000|0.742|7.3e-06|2.6e+00|8.7e+03| 2.323152e+03 -4.834414e+02| 0:0:00| chol  1  1
2|1.000|0.975|4.5e-06|6.9e-02|1.2e+03| 9.909793e+02 -1.527538e+02| 0:0:00| chol  1  1
3|0.896|0.911|1.5e-06|6.6e-03|1.6e+02| 8.750550e+01 -6.991137e+01| 0:0:00| chol  1  1
4|1.000|0.084|1.3e-06|7.1e-03|1.4e+02| 6.982174e+01 -6.345068e+01| 0:0:00| chol  1  1
5|0.832|0.451|3.1e-07|3.9e-03|8.5e+01| 4.154794e+01 -4.168044e+01| 0:0:00| chol  1  1
6|1.000|0.380|1.4e-07|2.4e-03|5.0e+01| 2.295562e+01 -2.589417e+01| 0:0:00| chol  1  1
7|1.000|0.600|3.5e-08|9.7e-04|1.7e+01| 8.164432e+00 -8.294511e+00| 0:0:00| chol  1  1
8|0.934|0.775|1.1e-08|2.2e-04|3.5e+00| 4.943187e+00  1.451602e+00| 0:0:00| chol  1  1
9|1.000|0.131|6.3e-10|2.4e-04|2.9e+00| 4.665170e+00  1.782327e+00| 0:0:00| chol  1  1
10|1.000|0.803|4.2e-10|4.7e-05|6.0e-01| 4.297611e+00  3.697458e+00| 0:0:00| chol  1  1
11|0.993|0.909|1.0e-10|4.3e-06|5.0e-02| 4.212306e+00  4.162719e+00| 0:0:00| chol  1  1
12|0.976|0.981|1.0e-11|7.9e-08|9.4e-04| 4.209768e+00  4.208832e+00| 0:0:00| chol  1  1
13|0.987|0.988|1.6e-11|7.3e-07|2.3e-05| 4.209706e+00  4.209695e+00| 0:0:00| chol  1  1
14|1.000|0.986|2.3e-13|1.8e-08|6.3e-07| 4.209705e+00  4.209705e+00| 0:0:00| chol  1  1
15|1.000|0.987|1.0e-14|4.9e-10|1.9e-08| 4.209705e+00  4.209705e+00| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations   = 15
primal objective value =  4.20970522e+00
dual   objective value =  4.20970521e+00
gap := trace(XZ)       = 1.86e-08
relative gap           = 1.98e-09
actual relative gap    = 1.18e-09
rel. primal infeas (scaled problem)   = 1.02e-14
rel. dual     "        "       "      = 4.89e-10
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual     "        "       "      = 0.00e+00
norm(X), norm(y), norm(Z) = 8.8e+00, 6.0e+00, 1.6e+01
norm(A), norm(b), norm(C) = 1.9e+01, 1.1e+01, 1.8e+01
Total CPU time (secs)  = 0.33
CPU time per iteration = 0.02
termination code       =  0
DIMACS: 3.6e-14  0.0e+00  1.8e-09  0.0e+00  1.2e-09  2.0e-09
-------------------------------------------------------------------

------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +4.20971

Computing the solution of the least-squares problem with variable weights...

Calling SDPT3 4.0: 304 variables, 40 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

num. of constraints = 40
dim. of socp   var  = 288,   num. of socp blk  = 16
dim. of linear var  = 16
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version  predcorr  gam  expon  scale_data
NT      1      0.000   1        0
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
0|0.000|0.000|2.9e+00|7.5e+00|2.8e+03| 0.000000e+00  0.000000e+00| 0:0:00| chol  1  1
1|1.000|1.000|1.8e-06|8.9e-02|4.0e+02|-2.456884e+00 -3.955094e+02| 0:0:00| chol  1  1
2|1.000|0.882|2.1e-07|1.8e-02|4.7e+01|-8.411599e+00 -5.393956e+01| 0:0:00| chol  1  1
3|0.840|0.952|8.4e-08|1.7e-03|4.9e+00|-1.781825e+01 -2.254401e+01| 0:0:00| chol  1  1
4|0.977|0.968|2.0e-08|1.4e-04|1.3e-01|-2.015599e+01 -2.027678e+01| 0:0:00| chol  1  1
5|0.965|0.982|3.6e-09|1.1e-05|3.3e-03|-2.020802e+01 -2.021033e+01| 0:0:00| chol  1  1
6|0.889|0.972|1.1e-08|3.2e-07|2.4e-04|-2.020955e+01 -2.020976e+01| 0:0:00| chol  1  1
7|0.897|0.958|3.2e-09|1.4e-08|2.4e-05|-2.020969e+01 -2.020971e+01| 0:0:00| chol  2  2
8|0.895|0.798|7.9e-10|3.5e-09|4.0e-06|-2.020970e+01 -2.020971e+01| 0:0:00| chol  1  1
9|0.830|0.881|2.3e-10|5.8e-10|6.5e-07|-2.020970e+01 -2.020971e+01| 0:0:00| chol  1  2
10|0.947|0.581|1.2e-11|2.9e-10|2.0e-07|-2.020971e+01 -2.020971e+01| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations   = 10
primal objective value = -2.02097052e+01
dual   objective value = -2.02097054e+01
gap := trace(XZ)       = 2.01e-07
relative gap           = 4.86e-09
actual relative gap    = 4.23e-09
rel. primal infeas (scaled problem)   = 1.23e-11
rel. dual     "        "       "      = 2.88e-10
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual     "        "       "      = 0.00e+00
norm(X), norm(y), norm(Z) = 1.7e+01, 2.9e+00, 3.7e+00
norm(A), norm(b), norm(C) = 1.4e+01, 1.5e+01, 6.4e+00
Total CPU time (secs)  = 0.12
CPU time per iteration = 0.01
termination code       =  0
DIMACS: 4.7e-11  0.0e+00  7.4e-10  0.0e+00  4.2e-09  4.9e-09
-------------------------------------------------------------------

------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +4.20971

Computing the solution of the quadratic program...

Calling SDPT3 4.0: 128 variables, 56 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

num. of constraints = 56
dim. of sdp    var  = 32,   num. of sdp  blk  = 16
dim. of linear var  = 80
*******************************************************************
SDPT3: Infeasible path-following algorithms
*******************************************************************
version  predcorr  gam  expon  scale_data
HKM      1      0.000   1        0
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
0|0.000|0.000|1.4e+01|9.3e+00|1.8e+04| 7.324334e+02  0.000000e+00| 0:0:00| chol  1  1
1|1.000|1.000|3.9e-06|9.4e-02|1.3e+03| 5.898749e+02 -6.551406e+02| 0:0:00| chol  1  1
2|0.986|0.923|2.8e-07|1.6e-02|6.6e+01| 2.708421e+01 -3.712863e+01| 0:0:00| chol  1  1
3|1.000|0.829|6.4e-08|3.5e-03|2.2e+01| 6.612363e+00 -1.492046e+01| 0:0:00| chol  1  1
4|1.000|1.000|1.2e-08|9.4e-05|9.5e+00| 8.063808e-01 -8.733407e+00| 0:0:00| chol  1  1
5|0.905|0.908|1.6e-09|1.7e-05|1.1e+00|-3.569766e+00 -4.710522e+00| 0:0:00| chol  1  1
6|1.000|1.000|3.3e-10|9.4e-07|5.2e-01|-3.918831e+00 -4.437390e+00| 0:0:00| chol  1  1
7|0.947|0.946|1.3e-10|1.4e-07|3.2e-02|-4.191363e+00 -4.223574e+00| 0:0:00| chol  1  1
8|0.982|0.986|5.6e-11|1.1e-08|7.9e-04|-4.209226e+00 -4.210015e+00| 0:0:00| chol  1  1
9|0.988|0.988|3.2e-11|1.4e-10|9.4e-06|-4.209699e+00 -4.209709e+00| 0:0:00| chol  1  1
10|0.995|1.000|1.5e-13|6.5e-12|1.4e-07|-4.209705e+00 -4.209705e+00| 0:0:00| chol  1  1
11|0.999|1.000|4.4e-11|1.0e-12|2.4e-09|-4.209705e+00 -4.209705e+00| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations   = 11
primal objective value = -4.20970521e+00
dual   objective value = -4.20970521e+00
gap := trace(XZ)       = 2.36e-09
relative gap           = 2.50e-10
actual relative gap    = 2.20e-10
rel. primal infeas (scaled problem)   = 4.39e-11
rel. dual     "        "       "      = 1.00e-12
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual     "        "       "      = 0.00e+00
norm(X), norm(y), norm(Z) = 1.4e+01, 2.9e+00, 9.4e+00
norm(A), norm(b), norm(C) = 2.0e+01, 1.7e+01, 1.1e+01
Total CPU time (secs)  = 0.17
CPU time per iteration = 0.02
termination code       =  0
DIMACS: 1.5e-10  0.0e+00  3.8e-12  0.0e+00  2.2e-10  2.5e-10
-------------------------------------------------------------------

------------------------------------------------------------
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.8309    0.8309