Tame your data – constrained spline regression

Harald Biller, Continental Automotive Technologies

The Wolfram System’s recent advances in Disciplined Convex Programming
are applied to a multi-objective optimization problem:
find a cubic spline function with small data distance and curvature.
We add positivity, monotonicity, and convexity constraints,
which we express as norm cone constraints.

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

Least squares fits

Get and fit some data

SeedRandom[2341];
data=Table[{x,Exp[2x]-2x-1+RandomReal[{-0.5,0.5}]},{x,-1,1,0.1}];
ListPlot[data]

builtInFit=Fit[data,{1,t,(t^2),(t^3)},t]

Plot[builtInFit,{t,-1,1},Prolog->{PointSize[0.01],Point[data]},PlotRange->{-0.5,4},PlotStyle->Thickness[0.005],ImageSize->1200,BaseStyle->{FontSize->20}]

The same fit from QuadraticOptimization

unconstrainedCubicPolynomial=
Block[
{f,a,x,y,sse,constraints,variables,solution},
f=Function[t,a[1]+a[2]t+a[3](t^2)+a[4](t^3)]; (* ansatz *) 
sse=Sum[x=pt[[1]];y=pt[[2]];(((y-f[x]))^2),{pt,data}]; (* sum of squared errors *) 
constraints={};
variables={a[1],a[2],a[3],a[4]};
solution=Echo @ QuadraticOptimization[sse,constraints,variables];
{f[t],sse}/.solution
]

builtInFit

Cubic spline functions

Recall that a spline function is a piecewise polynomial function which is “as differentiable as reasonable with respect to its degree”.
A cubic spline function has continuous second derivatives.

Clear[unconstrainedCubicSpline]
unconstrainedCubicSpline[n_]:=
Block[
{f,g,a,x,y,sse,variables,solution,internalSmoothnessConstraints},
f=Table[Function[t,Evaluate @ Sum[Indexed[a,{j,k}](t^(k-1)),{k,4}]],{j,n}];
g=Function[t,Evaluate @ Piecewise[Table[{f[[j]][t],-1+(2/n)(j-1)≤t≤-1+(2/n)j},{j,n-1}],f[[n]][t]]];
sse=Sum[x=pt[[1]];y=pt[[2]];(((y-g[x]))^2),{pt,data}];
internalSmoothnessConstraints=
Table[
{f[[j]][-1+(2/n)j]==f[[j+1]][-1+(2/n)j],f[[j]]'[-1+(2/n)j]==f[[j+1]]'[-1+(2/n)j],f[[j]]''[-1+(2/n)j]==f[[j+1]]''[-1+(2/n)j]},
{j,n-1}
];
variables={a∈Matrices[{n,4},Reals]};
solution=QuadraticOptimization[sse,internalSmoothnessConstraints,variables];
{g[t],sse}/.solution
]
unconstrainedCubicSpline[2]

TableForm @ Table[{(2^n),Last @ unconstrainedCubicSpline[(2^n)]},{n,0,4}]

Plot[Evaluate[{First @ unconstrainedCubicPolynomial,First @ unconstrainedCubicSpline[4]}],{t,-1,1},]

Equality constraints can often be eliminated explicitly. For example, use the dedicated BSplineBasis function.

Cubic spline functions with convexity constraints

Clear[convexCubicSpline]
convexCubicSpline[n_]:=
Block[
{f,g,a,x,y,sse,variables,solution,internalSmoothnessConstraints,convexityConstraints},
f=;
g=;
sse=;
internalSmoothnessConstraints=;
convexityConstraints=Append[Table[f[[j]]''[-1+(2/n)j]≥0,{j,n}],f[[1]]''[-1]≥0];
variables=;
solution=
QuadraticOptimization[
sse,
{internalSmoothnessConstraints,convexityConstraints},
variables
];
{g[t],sse}/.solution
]
convexCubicSpline[2]

Plot[Evaluate[{First @ unconstrainedCubicPolynomial,First @ unconstrainedCubicSpline[4],First @ convexCubicSpline[4]}],{t,-1,1},]

Non-negative spline functions

Example: quadratic polynomials

Quadratic polynomials which are positive on the interval [-1,1] form a convex cone – a “non-standard” cone.

Reduce[ForAll[x,-1≤x≤1,((a-c)/2)+b x+((a+c)/2)(x^2)>0],{a,b,c}]

RegionPlot3D[a>0&&-a<b<a&&c<(Sqrt[(a^2)-(b^2)]),{a,-0.2,1},{b,-1,1},{c,-1,1},BoxRatios->Automatic,AxesLabel->{"a","b","c"},BaseStyle->20,ImageSize->800,ViewPoint->{-3,-0.8,1.5}]

Some helpful theory: non-negative polynomials on the interval [-1,1]

Theorem (Pólya and Szegő 1976, attributed to Lukács; B.)

Let f be a real polynomial of degree n such that f(x)≥0 for all x∈[-1,1]. Then there exist real polynomials (Subscript[g, 1]), (Subscript[g, 2]), (Subscript[g, 3]), and (Subscript[g, 4]) such that

f(X)=(1+X)(Subscript[g, 1])(((X))^2)+(1-X)(Subscript[g, 2])(((X))^2)
	          =(Subscript[g, 3])(((X))^2)+(1-(X^2))(Subscript[g, 4])(((X))^2)

and such that all summands have degree at most n+1.

Corollary

(a)	If n=2m+1 then there exist positive semidefinite symmetric matrices (Subscript[A, 1]) and (Subscript[A, 2]) (even tridiagonal) such that

f(X)=(1+X)(1,X,(X^2),…,(X^m))(Subscript[A, 1])(({{1}, {X}, {⋮}, {(X^m)}}))+(1-X)(1,X,(X^2),…,(X^m))(Subscript[A, 2])(({{1}, {X}, {⋮}, {(X^m)}})).

Thus the set of such f has a linear parameterization from (ℝ^(2m+1))×(ℝ^(2m+1))≅(ℝ^(2n)) which is characterized by a linear matrix inequality.

Proof: choose (Subscript[c, 1]),(Subscript[c, 2])∈(ℝ^(m+1)) such that (Subscript[g, k])(X)=(1,X,(X^2),…,(X^m))(Subscript[c, k]), and set (Subscript[A, k])=(Subscript[c, k])(Subscript[c, k])ᵀ. (Make tridiagonal keeping anti-diagonal sums.)

(b) 	If n=2m then there exist positive semidefinite symmetric matrices (Subscript[A, 3]) and (Subscript[A, 4]) (even tridiagonal) such that

f(X)=(1,X,(X^2),…,(X^m))(Subscript[A, 3])(({{1}, {X}, {⋮}, {(X^m)}}))+(1-(X^2))(1,X,(X^2),…,(X^(m-1)))(Subscript[A, 4])(({{1}, {X}, {⋮}, {(X^(m-1))}})).

Thus the set of such f has a linear parameterization from (ℝ^(2m+1))×(ℝ^(2m-1))≅(ℝ^(2n)) which is characterized by a linear matrix inequality.

Cubic spline functions with convexity and positivity constraints

Clear[nonNegativeConvexCubicSpline]
nonNegativeConvexCubicSpline[n_]:=
Block[
{f,g,a,α,x,y,sse,variables,solution,internalSmoothnessConstraints,convexityConstraints,positivityConstraints},
f=Table[Function[t,Evaluate[(t-(-1+(j-1)(2/n))){1,t}.(({{(IndexedDefault[a, j,1]), (IndexedDefault[a, j,2])}, {(IndexedDefault[a, j,2]), (IndexedDefault[a, j,3])}})).{1,t}+(-1+j(2/n)-t){1,t}.(({{(IndexedDefault[a, j,4]), (IndexedDefault[a, j,5])}, {(IndexedDefault[a, j,5]), (IndexedDefault[a, j,6])}})).{1,t}]],{j,n}];
g=;
sse=;
internalSmoothnessConstraints=;
convexityConstraints=;
positivityConstraints=Table[{(({{(IndexedDefault[a, j,1]), (IndexedDefault[a, j,2])}, {(IndexedDefault[a, j,2]), (IndexedDefault[a, j,3])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,(({{(IndexedDefault[a, j,4]), (IndexedDefault[a, j,5])}, {(IndexedDefault[a, j,5]), (IndexedDefault[a, j,6])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},{j,n}];
variables={a∈Matrices[{n,6},Reals]};
solution=
SemidefiniteOptimization[
α,
Flatten[{sse≤α,internalSmoothnessConstraints,convexityConstraints,positivityConstraints}],
Append[variables,α]
];
{Expand[g[t]],α}/.solution
]
nonNegativeConvexCubicSpline[2]

TableForm @ Table[{(2^n),Last @ unconstrainedCubicSpline[(2^n)],Last @ convexCubicSpline[(2^n)],Last @ nonNegativeConvexCubicSpline[(2^n)]},{n,0,4}]

With[{n=4},
Plot[
Evaluate[{unconstrainedCubicPolynomial,unconstrainedCubicSpline[n],convexCubicSpline[n],nonNegativeConvexCubicSpline[n]}[[All,1]]],
{t,-1,1},
]
]

Automotive application: the Worldwide Harmonized Light Vehicles Test Cycle (WLTC)

With[
{β={},style={FontSize->40},size=2400,ar1=(1/6),ar2=0.4,pad={{140,40},{100,20}}},
With[
{
v=BSplineFunction[Join[{0},β,{0}],SplineDegree->3,SplineKnots->Range[-3,Length[β]+2]],
a=BSplineFunction[((Append[β,0]-Prepend[β,0])/3.6),SplineDegree->2,SplineKnots->Range[-2,Length[β]+1]]
},
With[
{diag=Plot[v[t],{t,0,1800},Ticks->None,ImageSize->1200,AspectRatio->(1/20),ImagePadding->{{20,20},{5,5}}]},
Manipulate[
TabView[
{
"speed and acceleration versus time"->
Dynamic @ Column[
{
Plot[v[t],{t,int[[1]],int[[2]]},PlotRange->{-10,140},PlotStyle->Thickness[0.0025],PlotRangePadding->None,PlotTheme->"Detailed",FrameLabel->{"time [s]","speed [km/h]"},PlotLegends->None,ImageSize->size,AspectRatio->ar1,ImagePadding->pad,BaseStyle->style
],
Plot[a[t],{t,int[[1]],int[[2]]},PlotRange->{-1.8,1.8},PlotStyle->Thickness[0.0025],PlotRangePadding->None,PlotTheme->"Detailed",FrameLabel->{"time [s]","acceleration [m/(s^2)]"},PlotLegends->None,ImageSize->size,AspectRatio->ar1,ImagePadding->pad,BaseStyle->style
]
}
],
"acceleration versus speed"->
Dynamic @ ParametricPlot[{v[t],a[t]},{t,int[[1]],int[[2]]},PlotRange->{{-10,140},{-1.8,1.8}},PlotPoints->Ceiling[1+((int[[2]]-int[[1]])/3)],PlotStyle->Thickness[0.0025],PlotRangePadding->None,PlotTheme->"Detailed",FrameLabel->{"speed [km/h]","acceleration [m/(s^2)]"},PlotLegends->None,ImageSize->size,AspectRatio->ar2,ImagePadding->pad,BaseStyle->style
]
},
LabelStyle->Bold,ImageSize->{2600,1200},Alignment->Center
],
{{dummy,0,""},diag&},
{{int,{0,1800},"time range [s]"},0,1800,30,ControlType->IntervalSlider,Appearance->"Labeled",Method->"Push",MinIntervalSize->30,ImageSize->{1200,100}},
LabelStyle->style
]
]
]
]

The complete picture

Summary

least squares approximation by: polynomial – cubic spline – convex cubic spline – non-negative convex cubic spline

Related constraints

Monotonicity is a sign constraint on the first derivative, so it can be imposed in the same way.

Prescribe for sub-intervals whether the approximating spline function is
positive or negative, increasing or decreasing, convex or concave.

Positive semidefinite 2×2 matrices form a 3-dimensional quadratic norm cone.
For cubic splines, we can replace SemidefiniteOptimization with SecondOrderConeOptimization.
Advantages: Xpress solver, …

We can minimize the maximum approximation error or the sum of absolute errors instead of the sum of squared errors.

Two objectives: approximation error versus integral curvature

We can obtain a smoother function when we increase the distance between the fit and the data.

More precisely:

determine the closest approximating function (Subscript[f, 0]) as before

allow an additional distance ((ℓ^1), (ℓ^2), or (ℓ^∞)) to the data, or allow an additional distance ((L^2) or (L^∞)) to (Subscript[f, 0])

minimize ∫(((f''))^2) under this new constraint.

The end.

Appendix

Cubic spline functions with convexity and positivity constraints using BSplineBasis

Prepare positivity constraint on [0,1]

Block[
{b,c,g,j,k,x,eq},
g=Echo @ Collect[(1-x){1,x}.(({{(IndexedDefault[c, 1]), ((IndexedDefault[b, 1])/2)}, {((IndexedDefault[b, 1])/2), (IndexedDefault[b, 2])}})).{1,x}+x{1,x}.(({{(IndexedDefault[b, 3]), ((IndexedDefault[b, 4])/2)}, {((IndexedDefault[b, 4])/2), (IndexedDefault[b, 5])}})).{1,x}-(∑_(k=1)^4)(IndexedDefault[c, k])(x^(k-1)),x];
eq=Thread[0==Rest @ CoefficientList[g,x]];
MatrixForm[(({{(IndexedDefault[b, 3]), ((IndexedDefault[b, 4])/2)}, {((IndexedDefault[b, 4])/2), (IndexedDefault[b, 5])}}))/.Simplify @ First @ Solve[eq,{(IndexedDefault[b, 3]),(IndexedDefault[b, 4]),(IndexedDefault[b, 5])}]/.{(IndexedDefault[b, k_])->(IndexedDefault[b, j,k]),(IndexedDefault[c, k_])->(IndexedDefault[c, j,k])}]
]

Implement fitting

Clear[nonNegativeConvexCubicBSpline]
nonNegativeConvexCubicBSpline[n_]:=
Block[
{knots,f,ff,g,h,a,b,c,α,t,x,y,sse,variables,solution,convexityConstraints,positivityConstraints},
knots=Range[-1-(6/n),1+(6/n),(2/n)];
f=Function[t,Evaluate @ Sum[(IndexedDefault[a, j])BSplineBasis[{3,knots},j-1,t],{j,n+3}]];
sse=Sum[x=pt[[1]];y=pt[[2]];(((y-f[x]))^2),{pt,data}]; 
ff=
Table[(IndexedDefault[a, j]),{j,n+3}].
Table[
g=First @ PiecewiseExpand @ BSplineBasis[{3,knots},j-1,t];
Table[Select[g,#[[2]]/.t->((knots[[k+3]]+knots[[k+4]])/2)&],{k,n}]/.{{}->0,{{h_,_}}->h},
{j,n+3}
];
c=Table[PadRight[CoefficientList[ff[[j]]/.t->(1-t)knots[[j+3]]+t knots[[j+4]],t],4],{j,n}];
convexityConstraints=Table[f''[knots[[j+4]]]≥0,{j,0,n}];
positivityConstraints=Table[{(({{(IndexedDefault[c, j,1]), ((IndexedDefault[b, j,1])/2)}, {((IndexedDefault[b, j,1])/2), (IndexedDefault[b, j,2])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0,(({{-(IndexedDefault[b, j, 1])+(IndexedDefault[c, j, 1])+(IndexedDefault[c, j, 2]), (1/2) ((IndexedDefault[b, j, 1])-(IndexedDefault[b, j, 2])+(IndexedDefault[c, j, 3]))}, {(1/2) ((IndexedDefault[b, j, 1])-(IndexedDefault[b, j, 2])+(IndexedDefault[c, j, 3])), (IndexedDefault[b, j, 2])+(IndexedDefault[c, j, 4])}}))(Underscript[\[VectorGreaterEqual], (SemidefiniteConeList[2])])0},{j,n}];
variables={a∈Vectors[n+3,Reals],b∈Matrices[{n,2},Reals]};
solution=
SemidefiniteOptimization[
α,
Flatten[{sse≤α,convexityConstraints,positivityConstraints}],
Append[variables,α]
];
{BSplineFunction[a/.solution,SplineKnots->knots][t],α/.solution}
]
nonNegativeConvexCubicBSpline[6]

TableForm @ Table[{(2^n),Last @ unconstrainedCubicSpline[(2^n)],Last @ convexCubicSpline[(2^n)],Last @ nonNegativeConvexCubicSpline[(2^n)],Last @ nonNegativeConvexCubicBSpline[(2^n)]},{n,0,4}]

With[{n=4},
Plot[
Evaluate[{unconstrainedCubicPolynomial,unconstrainedCubicSpline[n],convexCubicSpline[n],nonNegativeConvexCubicBSpline[n]}[[All,1]]],
{t,-1,1},
]
]

Compare with first implementation

With[
{n=8},
Plot[Evaluate[Last[nonNegativeConvexCubicSpline[n]]-Last[nonNegativeConvexCubicBSpline[n]]],{t,-1,1},]
]

With[{n=8},{RepeatedTiming[nonNegativeConvexCubicSpline[n]],RepeatedTiming[nonNegativeConvexCubicBSpline[n]]}]

Block[{n=8,f1,f2},
f1=First @ nonNegativeConvexCubicSpline[n];f2=First @ nonNegativeConvexCubicBSpline[n];
{RepeatedTiming[Plot[f1,{t,-1,1}]],RepeatedTiming[Plot[f2,{t,-1,1}]]}
]

Code for the cones on the title slide

CellPrint @ ExpressionCell[
DynamicModule[
{vp={1.3,-2.4,2.},vv={0,0,1},size=500,ax={Arrow @ Tube[{{0,0,0},{1.4,0,0}}],Arrow @ Tube[{{0,0,0},{0,1.4,0}}],Arrow @ Tube[{{0,0,0},{0,0,1.4}}]}},
Row[
{
Graphics3D[{Polygon[{Append[Table[{Cos[ϕ],Sin[ϕ],0},{ϕ,0,(π/2),(π/20)}],{0,0,0}],Append[Table[{Cos[ϕ],0,Sin[ϕ]},{ϕ,0,(π/2),(π/20)}],{0,0,0}],Append[Table[{0,Cos[ϕ],Sin[ϕ]},{ϕ,0,(π/2),(π/20)}],{0,0,0}]}],Arrowheads[.08],ax},PlotRange->{{-0.5,1.5},{-0.5,1.5},{-0.5,1.5}},ViewPoint->Dynamic[vp],ViewVertical->Dynamic[vv],ImageSize->size],Graphics3D[{Cone[{{0,0,1},{0,0,0}}],Arrowheads[.06],ax},PlotRange->{{-1.5,1.5},{-1.5,1.5},{-0.5,1.5}},ViewPoint->Dynamic[vp],ViewVertical->Dynamic[vv],ImageSize->size]
},
Background->White
]
],
"Subtitle"
]