% 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 4.0: 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  = 2550
*******************************************************************
   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|1.1e+04|5.1e+01|3.4e+07|-2.463801e-08  0.000000e+00| 0:0:00| chol  1  1 
 1|0.023|0.867|1.1e+04|6.9e+00|3.2e+07| 7.051700e+01 -3.679632e+03| 0:0:00| chol  1  1 
 2|0.918|0.974|8.6e+02|2.7e-01|2.7e+06| 1.401794e+01 -2.943432e+03| 0:0:00| chol  1  1 
 3|0.910|0.987|7.8e+01|3.1e-02|2.4e+05| 1.029536e+01 -6.757387e+02| 0:0:00| chol  1  1 
 4|0.961|0.969|3.1e+00|9.2e-03|1.0e+04| 1.085902e+01 -4.145262e+02| 0:0:00| chol  1  1 
 5|0.965|0.985|1.1e-01|9.8e-04|6.9e+02| 1.092364e+01 -3.723299e+02| 0:0:00| chol  1  1 
 6|0.229|0.314|8.2e-02|7.0e-04|5.6e+02| 1.109744e+01 -3.294966e+02| 0:0:00| chol  1  1 
 7|0.219|0.312|6.4e-02|1.7e-02|4.9e+02| 1.114592e+01 -3.045759e+02| 0:0:00| chol  1  1 
 8|0.720|0.287|1.8e-02|2.5e-02|3.4e+02| 1.122282e+01 -2.740088e+02| 0:0:00| chol  1  1 
 9|0.739|0.306|4.7e-03|2.1e-02|2.6e+02| 1.125567e+01 -2.344734e+02| 0:0:00| chol  1  1 
10|0.217|0.246|3.6e-03|1.7e-02|2.3e+02| 1.127324e+01 -2.101866e+02| 0:0:00| chol  1  1 
11|0.187|0.582|3.0e-03|7.7e-03|1.9e+02| 1.132396e+01 -1.716350e+02| 0:0:00| chol  1  1 
12|0.343|1.000|1.9e-03|5.9e-04|1.4e+02| 1.135751e+01 -1.171362e+02| 0:0:00| chol  1  1 
13|0.911|1.000|1.7e-04|3.9e-04|8.5e+01| 1.130069e+01 -7.287643e+01| 0:0:00| chol  1  1 
14|0.574|1.000|7.4e-05|3.5e-05|5.7e+01| 1.105172e+01 -4.493234e+01| 0:0:00| chol  1  1 
15|0.722|1.000|2.1e-05|1.5e-05|3.7e+01| 1.070173e+01 -2.550582e+01| 0:0:00| chol  1  1 
16|0.507|1.000|1.0e-05|4.1e-06|2.9e+01| 9.626004e+00 -1.850222e+01| 0:0:00| chol  1  1 
17|0.798|0.987|2.1e-06|2.1e-06|1.8e+01| 8.940952e+00 -8.677965e+00| 0:0:00| chol  1  1 
18|0.306|1.000|1.4e-06|4.1e-07|2.5e+01| 3.963923e+00 -2.041665e+01| 0:0:00| chol  1  1 
19|0.820|0.965|2.6e-07|3.0e-07|8.7e+00| 3.284605e+00 -5.242031e+00| 0:0:00| chol  1  1 
20|0.777|0.882|5.7e-08|8.7e-08|5.3e+00| 5.203090e-01 -4.702256e+00| 0:0:00| chol  1  1 
21|0.517|0.735|2.8e-08|3.4e-08|3.5e+00|-1.037220e-01 -3.567379e+00| 0:0:00| chol  1  1 
22|0.798|1.000|5.6e-09|5.5e-09|1.6e+00|-1.131497e+00 -2.689229e+00| 0:0:00| chol  1  1 
23|0.608|0.919|2.2e-09|1.6e-09|9.1e-01|-1.575050e+00 -2.467360e+00| 0:0:00| chol  1  1 
24|0.844|1.000|3.4e-10|4.7e-10|3.3e-01|-2.030276e+00 -2.351939e+00| 0:0:00| chol  1  1 
25|0.748|0.880|8.6e-11|2.1e-10|1.2e-01|-2.186506e+00 -2.299985e+00| 0:0:00| chol  1  1 
26|0.842|1.000|1.4e-11|3.6e-10|3.5e-02|-2.247701e+00 -2.281961e+00| 0:0:00| chol  1  1 
27|0.896|1.000|1.4e-12|3.5e-10|4.1e-03|-2.270685e+00 -2.274766e+00| 0:0:00| chol  1  1 
28|0.952|0.961|3.0e-13|3.7e-10|2.0e-04|-2.274089e+00 -2.274289e+00| 0:0:00| chol  1  1 
29|0.922|0.986|2.5e-12|5.4e-10|4.0e-05|-2.274236e+00 -2.274276e+00| 0:0:00| chol  1  1 
30|0.658|0.991|1.3e-12|1.0e-09|1.5e-05|-2.274253e+00 -2.274268e+00| 0:0:00| chol  1  1 
31|0.920|1.000|8.0e-13|1.3e-09|2.0e-06|-2.274264e+00 -2.274266e+00| 0:0:00| chol  1  1 
32|0.807|1.000|1.7e-12|2.1e-09|4.4e-07|-2.274265e+00 -2.274266e+00| 0:0:00| chol  1  1 
33|0.890|1.000|2.6e-11|1.5e-09|8.1e-08|-2.274266e+00 -2.274266e+00| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 33
 primal objective value = -2.27426576e+00
 dual   objective value = -2.27426584e+00
 gap := trace(XZ)       = 8.06e-08
 relative gap           = 1.45e-08
 actual relative gap    = 1.44e-08
 rel. primal infeas (scaled problem)   = 2.65e-11
 rel. dual     "        "       "      = 1.53e-09
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 2.7e+00, 3.2e+08, 2.6e+09
 norm(A), norm(b), norm(C) = 8.4e+01, 2.0e+00, 1.3e+02
 Total CPU time (secs)  = 0.43  
 CPU time per iteration = 0.01  
 termination code       =  0
 DIMACS: 2.6e-11  0.0e+00  2.3e-08  0.0e+00  1.4e-08  1.5e-08
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +2.27427