1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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)

$\endgroup$
  • $\begingroup$ I've used the idea here (but coded in Python with cvxopt) and it works a treat. Thanks Johan. $\endgroup$ – n00b Dec 19 '13 at 20:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.