Fit polynomial to the maximums of a function

Question

I want to fit a polynomial to noisy data so that the approximated polynomial is always >= the original data. For example:

x = linspace (-2, 6);
y = (x-2).^2 + 1 + 2 * randn (size (x));

function ret = delta (P, x, y)
  yP = polyval (P, x);
  d = yP - y;
  d (d < 0) *= 1000;
  ret = sumsq (d);
endfunction

P0 = polyfit (x, y, 2);
f = @(P) delta (P, x, y);
[P, FVAL] = sqp (P0, f)

xi = linspace (min(x), max(x), 100);
yi = polyval (P, xi);

plot(x, y, xi, yi);
grid on

Is there a better way/method which also works with higher order polynomials?

The easies way would be to just use polyfit and then calculate max(y-yi) and add this as offset but this wouldn't be optimal...

EDIT: I want to use GNU OCtave but added "matlab" as tag because the language and functions are similiar.

EDIT: based on tvo's answer and real data:

x = [10 20 30 40 50 60 80 100];
y = [0.2372, 0.1312, 0.0936, 0.0805, 0.0614, 0.0512, 0.0554, 0.1407];

function ret = delta (P, x, y)
  ret = sumsq (polyval (P, x) - y);
endfunction

f = @(P) delta (P, x, y);
h = @(P) polyval(P, x) - y;

P0 = polyfit (x, y, 3);
[P] = sqp (P0, f, [], h)

xi  = linspace (min(x), max(x));
yi  = polyval (P0, xi);
yio = polyval (P, xi);

plot(x, y, xi, yi, ";initial guess;", xi, yio, ";optimized;");
grid on

but as you can see, the optimized and evaluated poly has points < the orginal point which isn't allowed.

Answer 1

your method seems fine, and I see no reason it can't be used for higher-order polynomials actually. Please explain why if you think it can't be used.

You are using Octave's 'sqp' solver. Documentation is here: http://www.gnu.org/software/octave/doc/v4.0.1/Nonlinear-Programming.html

You may want to avoid multiplying the error by an arbitrary number (1000 in your example) when it is negative. This may fail for different data sets, especially if they are larger, i.e. more data points.

You could try to use the non-linear inequality constraint options that Octave's 'sqp' offers, i.e. the h(x)>=0 (see doc).

As objective function phi you could use a square norm error, as in your example, and add constraints of the form h(x)>=0 for every data point. Note that 'x' would be the polynomial coefficients you want to fit and h(x) is the polynomial evaluated at a specific data point.

For example:

phi = @(P) delta_mod (P, x, y); % mod: Don't increase the importance of negative residuals!!
h = @(P) polyval(P, x1) - y1;

Psol = sqp(P0, phi, [], h);

Notice that the constraint function 'h' ensures that the polynomial will lie above (x1,y1), and the objective function 'phi' will try to keep it as close as possible. You can extend 'h' to contain one constraint for every data point in your set (see doc).