Polynomial fit matlab with some constraints on the coefficients

Question

I have data that I should interpolate with a function which must be of the following kind:

f(x) = ax4 + bx2 + c

with a > 0 and b ≤ 0. Unfortunately, MATLAB's polyfit does not allow any constraints on the coefficients of the polynomial. Does anybody know if there is a MATLAB function to do this? Otherwise, how can I implement it?

Thank you very much in advance,

Elisabetta

Answer 1

You can try using fminsearch, fminunc defining your objective function manually.

Alternatively, you can define your problem slightly different:

f(x) = a2x4 - b2x2 + c

Now, the new a and b can be optimized for without constraints, while ensuring that the final a and b you are looking for are positive (negative resp.).

Answer 2

Without constraints, the problem can be written and solved as a simple linear system:

% Your design matrix ([4 2 0] are the powers of the polynomial)
A = bsxfun(@power, your_X_data(:), [4 2 0]);

% Best estimate for the coefficients, [a b c], found by 
% solving A*[a b c]' = y in a least-squares sense
abc = A\your_Y_data(:)

Those constraints will of course automatically be satisfied iff that constrained model indeed underlies your data. For example,

% some example factors
a = +23.9;
b = -15.75;
c = 4;

% Your model
f = @(x, F) F(1)*x.^4 + F(2)*x.^2 + F(3);

% generate some noisy XY data
x = -1:0.01:1;
y = f(x, [a b c]) + randn(size(x));

% Best unconstrained estimate a, b and c from the data
A = bsxfun(@power, x(:), [4 2 0]);
abc = A\y(:);

% Plot results
plot(x,y, 'b'), hold on
plot(x, f(x, abc), 'r')
xlabel('x (nodes)'), ylabel('y (data)')

However, if you impose constraints on data that are not accurately described by that constrained model, things might go wrong:

% Note: same data, but flipped signs 
a = -23.9;
b = +15.75;
c = 4;

f = @(x, F) F(1)*x.^4 + F(2)*x.^2 + F(3);

% generate some noisy XY data
x = -1:0.01:1;
y = f(x, [a b c]) + randn(size(x));

% Estimate a, b and c from the data, Forcing a>0 and b<0
abc = fmincon(@(Y) sum((f(x,Y)-y).^2), [0 0 0], [-1 0 0; 0 +1 0; 0 0 0], zeros(3,1));

% Plot results
plot(x,y, 'b'), hold on
plot(x, f(x, abc), 'r')
xlabel('x (nodes)'), ylabel('y (data)')

(this solution has a == 0, indicative of an incorrect model choice).

If the exact equality of a == 0 is a problem: there is of course no difference if you set a == eps(0). Numerically, this will not be noticeable for real-world data, but it's nonzero nonetheless.

Anyway, I have a suspicion that your model is not well chosen and the constraints are a "fix" to get everything to work, or your data should actually be unbiased/rescaled before trying to make any fit, or that some similar preconditions apply (I've often seen people do this sort of thing, so yes, I'm a bit biased in this respect :).

So...what are the real reasons behind those constraints?