Determine low-order polynomial lower bound

Question

I have a function $f$ I'd like to determine numerically and I have a bunch of $(x, y)$ pairs which approximate the function in the following sense: all of the points satisfy $f(x) \leq y$, most of the points have $f(x) = y$ (at least to the accuracy I need), a few of the points have $y$ significantly bigger than $f(x)$. (The points come from a simulated annealing which occasionally fails.)

I believe that $f$ can be well approximated by a low-order polynomial, so I'd like to find such which lies below all of the points and which hits as many as possible, but does not mind missing a few. That's a rather woolly thing to want, so I guess I should look for the largest which lies below all of the points, where largest would be with respect to the $L^p$ norm, say $p = 1$ of $2$, on an interval containing all of the $x$s.

Any pointers or references on how to find such a polynomial would be much appreciated.

Answer 1

In the language of regression, you want to be robust against outliers. Penalizing the 1-norm of the approximation error is a simple way to achieve this. This penalty, together with the constraint that the approximation satisfies $f(x_i)\leq y_i$ leads to a linear program (can probably be solved more efficient with specialized methods, but this is the obvious model). Proof-of-concept code below in MATLAB (using the toolbox YALMIP (disclaimer, developed by me))

x = (-1:0.1:1)';
y = x.^2+x.^3+.5*x.^4;
y(3) = 3;
y(11)= 7;
plot(x,y,'*')

sdpvar a b c
yApprox = a+b*x+c*x.^2;
Error = norm(y-yApprox,1);
Below = [yApprox <= y];
solvesdp(Below,Error);
hold on;
plot(x,double(yApprox))

Note that the model and data here basically yields the same solution if you switch to a quadratic objective

Error = (y-yApprox)'*(y-yApprox);

The reason is that the upper bound constraint on the approximation protects the solution from being pulled towards the outliers (which typically happens when using quadratic error measures)