Example 6.8: Spline fitting
n=4;
m=40;
randn('state',0);
u = linspace(-1,1,m);
v = 1./(5+40*u.^2) + 0.1*u.^3 + 0.01*randn(1,m);
a = -1/3; b = 1/3;
u1 = u(find(u<a)); m1 = length(u1);
u2 = u(find((u >= a) & (u<b))); m2 = length(u2);
u3 = u(find((u >= b))); m3 = length(u3);
A1 = vander(u1'); A1 = fliplr(A1(:,m1-n+[1:n]));
A2 = vander(u2'); A2 = fliplr(A2(:,m2-n+[1:n]));
A3 = vander(u3'); A3 = fliplr(A3(:,m3-n+[1:n]));
fprintf(1,'Computing splines in the case of L2-norm...');
cvx_begin
variables x1(n) x2(n) x3(n)
minimize ( norm( [A1*x1;A2*x2;A3*x3] - v') )
[1 a a^2 a^3]*x1 == [1 a a^2 a^3]*x2;
[0 1 2*a 3*a^2]*x1 == [0 1 2*a 3*a^2]*x2;
[0 0 2 6*a ]*x1 == [0 0 2 6*a ]*x2;
[1 b b^2 b^3]*x2 == [1 b b^2 b^3]*x3;
[0 1 2*b 3*b^2]*x2 == [0 1 2*b 3*b^2]*x3;
[0 0 2 6*b ]*x2 == [0 0 2 6*b ]*x3;
cvx_end
fprintf(1,'Done! \n');
fprintf(1,'Computing splines in the case of Linfty-norm...');
cvx_begin
variables xl1(n) xl2(n) xl3(n)
minimize ( norm( [A1*xl1;A2*xl2;A3*xl3] - v', inf) )
[1 a a^2 a^3]*xl1 == [1 a a^2 a^3]*xl2;
[0 1 2*a 3*a^2]*xl1 == [0 1 2*a 3*a^2]*xl2;
[0 0 2 6*a ]*xl1 == [0 0 2 6*a ]*xl2;
[1 b b^2 b^3]*xl2 == [1 b b^2 b^3]*xl3;
[0 1 2*b 3*b^2]*xl2 == [0 1 2*b 3*b^2]*xl3;
[0 0 2 6*b ]*xl2 == [0 0 2 6*b ]*xl3;
cvx_end
fprintf(1,'Done! \n');
u1s = linspace(-1.0,a,1000)';
p1 = x1(1) + x1(2)*u1s + x1(3)*u1s.^2 + x1(4).*u1s.^3;
p1l1 = xl1(1) + xl1(2)*u1s + xl1(3)*u1s.^2 + xl1(4).*u1s.^3;
u2s = linspace(a,b,1000)';
p2 = x2(1) + x2(2)*u2s + x2(3)*u2s.^2 + x2(4).*u2s.^3;
p2l1 = xl2(1) + xl2(2)*u2s + xl2(3)*u2s.^2 + xl2(4).*u2s.^3;
u3s = linspace(b,1.0,1000)';
p3 = x3(1) + x3(2)*u3s + x3(3)*u3s.^2 + x3(4).*u3s.^3;
p3l1 = xl3(1) + xl3(2)*u3s + xl3(3)*u3s.^2 + xl3(4).*u3s.^3;
us = [u1s;u2s;u3s];
p = [p1;p2;p3];
pl = [p1l1;p2l1;p3l1];
d = plot(us,p,'b-',u,v,'go', us,pl,'r--',...
[-1 -1], [-0.1 0.25], 'k--', [1 1], [-0.1 0.25], 'k--', ...
[a a], [-0.1 0.25], 'k--', [b b], [-0.1 0.25], 'k--');
title('Approximation using 2 cubic splines');
xlabel('u');
ylabel('f(u)');
legend('L_2 norm','Data points','L_{\infty} norm', 'Location','Best');
Computing splines in the case of L2-norm...
Calling SDPT3: 47 variables, 13 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------
num. of constraints = 13
dim. of socp var = 41, num. of socp blk = 1
dim. of free var = 6 *** 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|2.7e+00|3.7e+01|5.3e+03| 0.000000e+00| 0:0:00| chol 1 1
1|1.000|0.992|5.5e-06|4.0e-01|1.6e+01|-2.240770e+00| 0:0:00| chol 1 1
2|1.000|0.964|5.9e-07|2.4e-02|2.9e-01|-8.901858e-02| 0:0:00| chol 1 1
3|0.940|0.914|2.5e-07|3.0e-03|1.6e-02|-1.172493e-01| 0:0:00| chol 1 1
4|0.986|0.986|1.3e-07|1.4e-04|2.1e-04|-1.165695e-01| 0:0:00| chol 1 1
5|0.989|0.988|7.2e-09|8.8e-06|4.1e-06|-1.165984e-01| 0:0:00| chol 1 1
6|0.992|0.985|1.8e-09|1.4e-07|6.4e-08|-1.166033e-01| 0:0:00| chol 1 1
7|1.000|0.961|1.3e-11|4.7e-09|1.4e-09|-1.166034e-01| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations = 7
primal objective value = -1.16603357e-01
dual objective value = -1.16603352e-01
gap := trace(XZ) = 1.38e-09
relative gap = 1.12e-09
actual relative gap = -4.42e-09
rel. primal infeas = 1.27e-11
rel. dual infeas = 4.74e-09
norm(X), norm(y), norm(Z) = 1.5e+00, 1.9e+00, 1.6e-01
norm(A), norm(b), norm(C) = 1.3e+01, 2.0e+00, 1.7e+00
Total CPU time (secs) = 0.2
CPU time per iteration = 0.0
termination code = 0
DIMACS: 1.3e-11 0.0e+00 6.8e-09 0.0e+00 -4.4e-09 1.1e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.116603
Done!
Computing splines in the case of Linfty-norm...
Calling SDPT3: 86 variables, 13 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------
num. of constraints = 13
dim. of socp var = 80, num. of socp blk = 40
dim. of free var = 6 *** 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|2.8e+01|4.5e+01|5.6e+03| 0.000000e+00| 0:0:00| chol 1 1
1|1.000|0.990|1.2e-05|5.3e-01|1.8e+01|-3.488167e+00| 0:0:00| chol 1 1
2|1.000|0.983|1.4e-05|1.8e-02|1.8e-01|-5.808630e-02| 0:0:00| chol 1 1
3|1.000|0.801|1.6e-06|4.5e-03|3.8e-02|-2.326026e-02| 0:0:00| chol 1 1
4|0.530|0.772|2.7e-06|1.1e-03|2.1e-02|-3.011469e-02| 0:0:00| chol 1 1
5|0.982|0.572|4.5e-08|4.8e-04|6.6e-03|-3.299069e-02| 0:0:00| chol 1 1
6|1.000|0.809|1.3e-08|1.1e-04|9.7e-04|-3.230400e-02| 0:0:00| chol 1 1
7|0.985|0.981|1.5e-09|1.1e-05|1.9e-05|-3.204242e-02| 0:0:00| chol 1 1
8|0.990|0.989|7.1e-11|2.1e-07|3.3e-07|-3.203834e-02| 0:0:00| chol 1 1
9|0.990|0.989|3.5e-12|3.6e-09|5.1e-09|-3.203833e-02| 0:0:00|
stop: max(relative gap, infeasibilities) < 1.49e-08
-------------------------------------------------------------------
number of iterations = 9
primal objective value = -3.20383318e-02
dual objective value = -3.20383363e-02
gap := trace(XZ) = 5.14e-09
relative gap = 4.83e-09
actual relative gap = 4.22e-09
rel. primal infeas = 3.52e-12
rel. dual infeas = 3.60e-09
norm(X), norm(y), norm(Z) = 6.8e-01, 2.4e+00, 2.4e-01
norm(A), norm(b), norm(C) = 1.5e+01, 2.0e+00, 1.7e+00
Total CPU time (secs) = 0.2
CPU time per iteration = 0.0
termination code = 0
DIMACS: 3.5e-12 0.0e+00 5.1e-09 0.0e+00 4.2e-09 4.8e-09
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +0.0320383
Done!