Tame Your Data—Constrained Spline Regression

Hello everybody, and welcome to my presentation about constraint spline regression. It's about disciplined convex programming. That's convex optimization with specific constraints. Every constraint is that some standard cone, one of a small number of specific cones, must contain an affine function of the variables. Such standard cones are all vector components are non negative, the picture on the left, or the ice cream cone with respect to the Euclidean norm in n dimensions, the picture on the right or the cone of semidefinite matrices. These are the optimization problems which couldn't today be solved in an accurate, fast and robust way. Thankfully, the Wolfram system makes this available. We apply this to cubic spline functions, piecewise cubic polynomial functions with continuous second derivatives. We approximate a data series, just a number of points in the plane. We may also minimize the integral curvature. And then we add constraints such as the spline function should be positive or increasing or convex left turned. And the main point of this talk is the bottom line of the abstract. These constraints turn out to be standard cone constraints. And so Mathematica's framework in disciplined convex optimization solves this problem. I work for Continental, the global automotive supplier, a large and colorful company with all these friendly people and with these basic figures. Last year we had just under 40 billion euro sales, just under 200,000 employees all over the world. About a third of our sales comes from tires for passenger cars and trucks, which you may know. Almost half comes from mechatronic parts for passenger cars, such as the electronic brake control systems for which I work. And the rest is industrial rubber, such as conveyor belts. I work, as I said, for breaks in the business area of safety and motion vision. Zero is what we're aiming at. Zero fatalities, zero injuries, one day, zero road accidents. And we support that by manifold products, some of which you see here, for instance, electronic brake control systems. What I'm specializing in, that's, that's down here. And one slide to present myself. My background is mathematics and physics, studied in various German places and in London. Got my PhD from Stuttgart. Then I joined Continental almost 20 years ago in the development of electronic brake systems. First I did component simulation, that is PDE modeling. Then I moved on to system simulation, that is modeling with ordinary differential equations. And that gave me a specific view on system design as well. So that this kind of system modeling work has also influenced the break systems which Continental has brought out over this period in a significant way. Now let's come to the data fitting problem. Let's first just get some data, add some noise to this function here. And look at these 21 data points here. And we can fit these with a linear combination of these basis functions, I.e. by a cubic polynomial or whatever we choose. So let's look at this cubic polynomial fit, and that already fits the data quite nicely. And why would I want to fit these data? As always, I want to evaluate between points, I want to extrapolate, and I also want a smooth approximation. I want to get the derivative. For instance, in that talk just an hour ago, from finite elements to system modeler, the data points would have been displacements and corresponding forces. And that was the problem, to fit a function and to obtain the derivative, because a spring stiffness, the derivative of force versus displacement, was what they were really interested in. So fitting is really ubiquitous in modeling. And this polynomial fit in that example is really not bad, but you want maybe want something more flexible, more general. And so instead of the built in fit function, we're moving on to quadratic optimization, a more general, more flexible function. So for that we write out the cubic polynomial as an ansatz function. Then we have the sum of squared errors as an objective. So we sum over all the points, call the first coordinate X and the second Y, and evaluate Y minus F of X squared. We don't have any constraints in this first example. And the variables which we want to find in the end are the polynomial coefficients a1 to a4. And then we ask Mathematica to do the quadratic optimization, where the objective is the sum of squared errors. The constraints are voids here, and the variables are the polynomial coefficients. And we ask for the expression of the result and also for the sum of squared errors. Here we go. Yeah, hopefully that's the same polynomial as before. Of course it is. Right? Yeah. And then we want to have a more complicated answers, namely a cubic spline, not a single polynomial, but a piecewise polynomial function. And such a piecewise polynomial function is called the spline function if it is as differentiable as reasonable. For instance, if the polynomial degree is 3, we have cubic functions, then we want the second derivative to be continuous. And so we set up things in a similar way. This time it's a mathematical function with an argument reflecting the number of sub intervals. And we don't have one ansatz polynomial, we have n ansatz polynomials for the N sub intervals, So a table of them and the coefficients now form a matrix. And I use this index notation here, where j runs from 1 to n over the intervals and k runs from 1 to 4 over the exponent. And then G, the overall ansatz function is a piecewise function made up of these polynomial pieces on the sub intervals of length two over n. The sum of squared errors is as before. Now we do have constraints, namely internal smoothness. We have these polynomial pieces and at all the internal boundaries between the sub intervals and these are these numbers here, minus 1 plus 2 over n times j, where j ranges from 1 to n minus 1. And so the number here ranges from minus 1 plus 2 over n up to plus 1 minus 2 over n. And at all the boundary points we want the left and right interval polynomial to be equal and the first derivative of that polynomial and the second derivative of that polynomial. These are the internal smoothness constraints. The variables are now a matrix variable and quadratic optimization handles this nicely. We'll just evaluate the example of two sub intervals. But the only interval boundary is at 0, and from minus 1 to 0 we have the first term, and from 0 to plus 1 we have the second term. And internal smoothness in this case means that only the cubic summons can differ. And the more subintervals we take, the smaller is the error which we expect. So let's just set a table of the second part of the result. It runs quickly and it does decrease the error, of course. And yeah, just demonstrate what it looks like for four sub intervals. So the blue was our original polynomial and the yellowish one is the cubic spline. And it's closer to the data. And so in that example it's not so much closer, but that can be an important example advantage. Now you would probably, if you really implemented that, eliminate the equality constraints and use the B spline basis for performance and also matrix condition reasons that would be favorable. But for the purpose of this talk, let's just stick to that easy to read notation which I used before. And another thing which I really want to change is I don't want that wiggly curve here. So this is closer to the data, but I don't want the curvature sign to change from right to left to right to left. Again, I want to add the additional constraint that the solution is convex. So I set up everything as before, but I add convexity constraints. These are very easy because the second derivative is a piecewise linear function. And the convexity constraint is that the second derivative is always positive. But for this it's necessary and sufficient that it is positive on the interval boundaries on the breakpoints, right? And that's what I write out here as a convexity constraint. And as before, I get such a solution which now fits the data better than the polynomial did but it is convex, which is even more than what the polynomial was. And convexity can be very important. For physical reasons or just for smoothness reasons. Maybe I want more. So the green curve is the best we got up to now. It's convex, but it goes into the negative. And I may not want that either for physical reasons. I may want to constrain the function to be non negative. For instance, because it's a spring stiffness, whatever the y value is. I may be convinced that for physical or for modeling reasons it must be non negative. And that takes me to the second part of my talk. Because non negativity. Yeah, that's more complicated when you leave the realm of linear functions, P size linear functions. For example, already for quadratic polynomials, you can hope to work out by hand a condition on the coefficients of this quadratic polynomial here not to be negative on that interval. Actually the Wolfram system can do it for you. So in the degree two case the term is easy enough. For degree three it wouldn't actually work anymore. And if you visualize that, you see that you do get a convex cone with tip at zero here. But it's not a standard cone. So it's not in the framework of disciplined convex optimization. Unless you twist it a bit. Unless you write it as a sum of an ice cream cone and a half line. These are standard cones, but that involves introducing an auxiliary variable, three variables for the ice cream cone and one for the half line. And actually turns out that that generalizes to higher degree polynomials. So we are interested in cubics now, but you can use it up and up. That's a bit of theory which I need some time to dig out, which I even could improve upon a little bit. So when is a polynomial of degree n non negative on that reference interval from minus 1 to 1? So this is definitely a non negative polynomial. You have this square here that's non negative, and you multiply by another typical polynomial which on that interval will arise from 0 to 2, or with this polynomial which will fall from 2 to 0. So this is definitely non negative polynomial on the reference interval. And the theorem is that surprisingly every non negative polynomial can be written in that way. And the two summons can be chosen such that their degree is at most one above the degree of F. And if we don't take 1 +x and 1x as the defining polynomial, so to speak, but 1 and 1 minus x squared, the same theorem holds. We also have this representation and again the same degree bound. Now this is even degree. The summands are Even degree, the sum can be odd degree. If the highest terms cancel, that's what must happen. And these are odd degree summons. And depending on whether the degree of f is odd or even, you use the representation which is more suitable. So for odd degrees it's more economical to use the first representation. And then there's the observation that you can write the first representation which is still in blue here, in this way, where a 1 and a 2 are positive semidefinite matrices. You can even choose tridiagonal matrices. Namely you write the polynomial G1 or G2 in such a way with a coefficient vector ck. And then you see that you have CK times CK transposed here. And when you think about which matrices lead to the same polynomial here you see that you can replace a one with a tridiagonal matrix with the same anti diagonal sums. So this is an M +1 times M +1 matrix, but actually it only has two M +1 free entries. So what we are having here is a parameterization of the space of positive polynomials from R to the 2m plus 1 squared characterized by semidefinite matrix inequalities. And that's what builds the bridge to the optimization functions, as I announced. Right? And a similar statement for the even degree case. So the corresponding function for non negative convex cubic splines would contain the internal smoothness constraints. We also want to keep things convex. But then we write out every sub polynomial in the way we saw in the theorem. Actually we just changed the model polynomials. We shifted the interval, so to speak. But anyway, if these matrices are positive definite, then this will be a positive polynomial on the sub interval number J. And every positive polynomial comes from such a positive definite pair of coefficient matrices. So the positivity constraints here are just that. These coefficients up there now are positive with respect to the cone of positive semidefinite two by two matrices. And yeah, so we have a a few additional variables here. It used to be N4 and now it is N6 and. It works fine Also fast enough, positivity constraint gives a larger sum of squared error. So these columns are like before unconstrained convex constraint and now positive end convex constraint. So the additional constraint enforces the larger sum of squared errors, but it decreases over the number of intervals. Of course. And when you look at the result, we now have the red curve which is convex and non negative, just as we wanted. Let's briefly look at one automotive application at least after all this theory. So this is the world harmonized like light vehicles, test cycle for fuel consumption tests, etcetera that's given as a list of velocities. But if you actually want to drive it, you want the list of accelerations and decelerations. And yeah, so this is a spline fit to the data and you can zoom in and convince yourself of the shape of the decelerations and accelerations. I also like to see deceleration versus speed. So the decelerations are not very high and they only happen at. At low speeds. Right, Right. So that's one example. And now I'd like to summarize. This is my last slide. I'd like to summarize and also to extend the theory a little bit. So what I have shown you is the squares approximation by a cubic polynomial, a cubic spline, then a convex cubic spline, and then a non negative convex cubic spline where the non negativity constraint involved that theory about representation of non negative polynomials in terms of two positive semidefinite matrices. Now monotonicity, that's a sign constraint on the first derivative. So we can impose it on the same way. We would then probably prescribe for sub intervals only whether the approximating spline function is positive or negative, increasing or decreasing, convex or concave. Another interesting thing is that in low dimensions, positive disseminative two by two matrices are actually the ice cream cone from the title picture. So we can use second order cone optimization, which is a bit more specialized and allows us to use more solvers, maybe manipulate the problem in a more explicit way, things like that. We could also replace the sum of squared errors with the maximum approximation error or the sum of absolute errors. It's all straightforward in the same framework. And we can add a second objective, namely the integral curvature. We could further smoothen the result by allowing a larger error and minimizing the integral curvature instead, as I wrote down here. So we would determine the closest approximating function, f 0 first and then we would relax the distance to the data which we have obtained here in the L1, L2 or L infinity sense. Or we would just allow an additional distance to the first fit and minimize the sum, sorry, the integral over the square of the second derivative, the integral curvature measure under this new constraint, and that's an additional smoothing. So all this is now open to us after the development which I showed, which was mainly setting up things suitable for the disciplined convex programming framework in Mathematica, and then finding that bit of theory which allows us to write non negativity of a polynomial in an interval as a semidefinite cone constraint. And that's the end of it. Thank you very much for listening up to now. And now look to the other screen and see the many questions which you have hopefully posed. Yeah, I'm glad to see so many questions. Anton Antonov is asking about design matrices, but I don't really know that. Please get in touch. Then there are a few questions about neural networks or other machine learning, But. So actually I don't know what quantile regression is again from Anton Antonov, So I don't know very much about it. I thought about it before. Gaussian regression is something you hear a lot about. But for a Gaussian process to be strictly positive is not a natural constraint. But there may be a relation. But when you experiment with predict and Gaussian processes that's also much, much slower. So you've seen that this process here is very fast on these small examples and it goes up to definitely a few hundred intervals with the built in solvers and probably a few tens of thousands of intervals with the commercial solvers which are also available through Wolfram. So you can treat large problems in a short time with very elementary and easy to understand mathematical concepts. But other machine learning methods may add a different view here. Yeah, so please, please get in touch and discuss. Right, let me look on my clock. I am afraid I used up my half hour. Ah, an example for non negativity. Well, actually a spring constant is not such a bad example. Or if you have a force versus travel characteristic of a spring, you would maybe want to use that as a travel versus force characteristic. You want to solve for X in terms of Y and then you depend on the function being strictly increasing. Things like that. Right, yeah. So we'll probably need to go. Thank you very much for listening and have a great rest of the conference. Bye bye.