Figure 6.24: Fitting a convex function to given data

% Section 6.5.5
% Boyd & Vandenberghe "Convex Optimization"
% Original by Lieven Vandenberghe
% Adapted for CVX by Argyris Zymnis - 11/27/2005
%
% Here we find the convex function f that best fits
% some given data in the least squares sense.
% To do this we solve
%     minimize    ||yns - yhat||_2
%     subject to  yhat(j) >= yhat(i) + g(i)*(u(j) - u(i)), for all i,j

clear

% Noise level in percent and random seed.
rand('state',29);
noiseint=.05;

% Generate the data set
u = [0:0.04:2]';
m=length(u);
y = 5*(u-1).^4 + .6*(u-1).^2 + 0.5*u;
v1=u>=.2;
v2=u<=.6;
v3=v1.*v2;
dipvec=((v3.*u-.4*ones(1,size(v3,2))).^(2)).*v3;
y=y+40*(dipvec-((.2))^2*v3);

% add perturbation and plots the input data
randf=noiseint*(rand(m,1)-.5);
yns=y+norm(y)*(randf);
figure
plot(u,yns,'o');

% min. ||yns-yhat||_2
% s.t. yhat(j) >= yhat(i) + g(i)*(u(j) - u(i)), for all i,j
cvx_begin
    variables yhat(m) g(m)
    minimize(norm(yns-yhat))
    subject to
        yhat*ones(1,m) >= ones(m,1)*yhat' + (ones(m,1)*g').*(u*ones(1,m)-ones(m,1)*u');
cvx_end

nopts =1000;
t = linspace(0,2,nopts);
f = max(yhat(:,ones(1,nopts)) + ...
      g(:,ones(1,nopts)).*(t(ones(m,1),:)-u(:,ones(1,nopts))));
plot(u,yns,'o',t,f,'-');
axis off
%print -deps interpol_convex_function2.eps
 
Calling SDPT3: 2602 variables, 103 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints = 103
 dim. of socp   var  = 52,   num. of socp blk  =  1
 dim. of linear var  = 2548
 dim. of free   var  =  2 *** 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|1.1e+04|1.4e+02|1.3e+07| 0.000000e+00| 0:0:00| chol  1  1 
 1|0.162|0.575|9.1e+03|5.9e+01|1.1e+07|-1.361373e+02| 0:0:00| chol  1  1 
 2|0.600|0.110|3.6e+03|5.2e+01|6.2e+06|-1.784641e+02| 0:0:00| chol  1  1 
 3|0.569|0.598|1.6e+03|2.1e+01|3.5e+06|-2.706056e+02| 0:0:00| chol  1  1 
 4|0.788|0.860|3.3e+02|3.0e+00|7.7e+05|-1.886178e+02| 0:0:00| chol  1  1 
 5|0.717|0.450|9.4e+01|1.6e+00|2.4e+05|-1.823767e+02| 0:0:00| chol  1  1 
 6|0.811|0.794|1.8e+01|3.4e-01|4.9e+04|-1.687413e+02| 0:0:00| chol  1  1 
 7|0.923|0.560|1.4e+00|1.5e-01|4.3e+03|-1.718447e+02| 0:0:00| chol  1  1 
 8|0.937|0.966|8.6e-02|5.2e-03|5.1e+02|-1.407962e+02| 0:0:01| chol  2  1 
 9|0.575|0.118|3.6e-02|2.2e-02|3.8e+02|-1.333830e+02| 0:0:01| chol  1  1 
10|0.499|0.634|1.8e-02|1.5e-02|2.2e+02|-8.031300e+01| 0:0:01| chol  1  1 
11|0.466|0.442|9.7e-03|1.2e-02|1.7e+02|-6.527428e+01| 0:0:01| chol  1  1 
12|0.346|0.915|6.4e-03|3.0e-03|1.2e+02|-4.489533e+01| 0:0:01| chol  1  1 
13|0.851|0.803|9.5e-04|1.9e-03|8.1e+01|-2.824693e+01| 0:0:01| chol  1  1 
14|0.493|0.440|4.8e-04|1.2e-03|7.1e+01|-2.376340e+01| 0:0:01| chol  1  1 
15|0.672|0.224|1.6e-04|1.1e-03|6.3e+01|-2.021869e+01| 0:0:01| chol  2  2 
16|0.779|0.246|3.5e-05|8.2e-04|5.6e+01|-1.721019e+01| 0:0:01| chol  2  2 
17|0.148|0.270|3.0e-05|6.1e-04|4.8e+01|-1.355218e+01| 0:0:01| chol  1  1 
18|0.881|0.110|3.5e-06|5.5e-04|4.5e+01|-1.267635e+01| 0:0:01| chol  2  1 
19|0.076|0.423|3.3e-06|3.2e-04|3.6e+01|-7.935498e+00| 0:0:01| chol  1  1 
20|0.997|0.103|2.0e-07|2.9e-04|3.2e+01|-7.829282e+00| 0:0:01| chol  1  1 
21|0.209|0.458|3.6e-07|1.5e-04|2.5e+01|-4.387416e+00| 0:0:01| chol  1  1 
22|1.000|0.138|3.4e-07|1.3e-04|2.1e+01|-5.002619e+00| 0:0:01| chol  1  1 
23|1.000|0.486|2.5e-07|6.9e-05|1.5e+01|-3.747187e+00| 0:0:01| chol  2  2 
24|1.000|0.374|1.4e-07|4.3e-05|9.3e+00|-3.840043e+00| 0:0:01| chol  2  2 
25|0.924|0.263|1.3e-07|3.2e-05|7.7e+00|-3.432751e+00| 0:0:01| chol  2  2 
26|1.000|0.382|7.2e-08|2.0e-05|5.0e+00|-3.194082e+00| 0:0:02| chol  2  2 
27|1.000|0.306|5.2e-08|1.4e-05|3.7e+00|-2.948767e+00| 0:0:02| chol  2  1 
28|1.000|0.365|3.1e-08|8.7e-06|2.4e+00|-2.758707e+00| 0:0:02| chol  2  1 
29|1.000|0.326|2.5e-08|5.9e-06|1.7e+00|-2.626543e+00| 0:0:02| chol  2  1 
30|1.000|0.359|1.3e-08|3.8e-06|1.1e+00|-2.524879e+00| 0:0:02| chol  2  2 
31|1.000|0.354|7.6e-09|2.4e-06|7.3e-01|-2.454370e+00| 0:0:02| chol  2  2 
32|1.000|0.412|4.0e-09|1.4e-06|4.2e-01|-2.393915e+00| 0:0:02| chol  2  1 
33|1.000|0.410|4.2e-09|8.5e-07|2.4e-01|-2.351423e+00| 0:0:02| chol  1  2 
34|1.000|0.535|1.5e-09|3.0e-06|1.1e-01|-2.312188e+00| 0:0:02| chol  1  1 
35|1.000|0.897|1.1e-10|2.2e-06|1.3e-02|-2.274944e+00| 0:0:02| chol  1  1 
36|0.533|0.870|1.3e-10|2.8e-07|7.1e-03|-2.273228e+00| 0:0:02| chol  1  1 
37|0.907|0.946|8.0e-11|1.5e-07|2.8e-03|-2.273868e+00| 0:0:02| chol  1  1 
38|0.910|0.976|7.9e-12|5.9e-08|2.1e-04|-2.274201e+00| 0:0:02| chol  1  1 
39|0.930|0.976|1.8e-11|4.4e-09|3.2e-05|-2.274258e+00| 0:0:02| chol  1  1 
40|0.927|0.969|1.7e-11|6.7e-10|6.6e-06|-2.274265e+00| 0:0:02| chol  1  1 
41|1.000|0.982|7.5e-12|1.4e-10|4.3e-07|-2.274266e+00| 0:0:02| chol  1  1 
42|1.000|0.925|4.1e-12|3.8e-11|3.8e-08|-2.274266e+00| 0:0:02|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 42
 primal objective value = -2.27426581e+00
 dual   objective value = -2.27426585e+00
 gap := trace(XZ)       = 3.83e-08
 relative gap           = 6.91e-09
 actual relative gap    = 6.73e-09
 rel. primal infeas     = 4.05e-12
 rel. dual   infeas     = 3.76e-11
 norm(X), norm(y), norm(Z) = 2.7e+00, 8.0e+05, 6.6e+06
 norm(A), norm(b), norm(C) = 8.4e+01, 2.0e+00, 1.9e+01
 Total CPU time (secs)  = 2.4  
 CPU time per iteration = 0.1  
 termination code       =  0
 DIMACS: 4.1e-12  0.0e+00  8.7e-11  0.0e+00  6.7e-09  6.9e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +2.27427