MATLAB Answers

Constrained Polynomial Regression

140 views (last 30 days)
Jared
Jared on 6 Mar 2012
Answered: Kin Sung Chan on 30 Sep 2019
Hey all,
I'm sure someone will have a quick answer for my probably dumb question. Performing polynomial least squares regression of a set of [x,y] data. Using polyfit, polyval, and corrcoef, I solved for the standard regression.
My problem is when I need to solve for a constrained polynomial such as: f(x)=a3*x.^3+a1*x+a0
So what am I missing that leaves me with no understanding how to solve this problem.

  1 Comment

Michael Nemesure
Michael Nemesure on 15 Mar 2017
One way to achieve a polynomial fit with some coefficients constrained is to use the psedo-inverse pinv on an appropriately modified Vandermonde matrix. An example: for a 4th order fit to n data points (x,y) with linear and quadratic coefficents fixed at p1 and p2, compute
pinv([ones(n,1) zeros(n,1) zeros(n,1) X.^3 X.^4])* (Y - p1*X - p2*X.^2)
The resulting column vector will have 0's in the constrained locations to be replaced by p1 and p2.

Sign in to comment.

Accepted Answer

Geoff
Geoff on 6 Mar 2012
So you mean that you want to constrain a cubic polynomial regression such that the x^2 term is zero?
For this I would use fsolve if you have the optimization toolbox.
If you don't, there is always the Nelder-Mead algorithm implemented as the function fminsearch.
You'd go something like this:
constraint = [1,0,1,1]; % Ignore possible x^2 term
polyFunc = @(p) polyval(p.*constraint,x);
objectiveFunc = @(p) (y - polyFunc(p)).^2;
p0 = [1, 0, 2, 3]; % It pays to have a realistic initial guess
p = fminsearch( objectiveFunc, p0 );
Instead of using polyval you could always just write the polynomial formula verbatim. However, this way you get a result that is consistent with other polynomial functions.
Warning: Nelder-Mead can be a hairy beast. Read the documentation, play around, and find ways to tame it.

  2 Comments

Geoff
Geoff on 6 Mar 2012
Sorry, my mistake!
The objective function should be the SUM of squares.
objectiveFunc = @(p) sum((y - polyFunc(p)).^2);
Also, just because you constrained x^2 term to zero and supplied zero starting value for it, doesn't mean the algorithm won't output a non-zero value. So you might want to multiply the result by the constraint again:
p = p .* constraint;
Jared
Jared on 6 Mar 2012
This way seemed to work like a dream and it uses methods that I can somewhat understand.

Sign in to comment.

More Answers (6)

Richard Willey
Richard Willey on 6 Mar 2012
Curve Fitting Toolbox allows you to specify constraints for individual regression coefficients.
With this said and done, could you please elaborate why you need to constrain a regression coefficient to be precisely equal to zero rather than simply leaving that term out of your model?

  0 Comments

Sign in to comment.


Mark Selby
Mark Selby on 6 Mar 2012
The simple first principles way is as follows...
x=(1:10)'; % invent some data for prediction
y = x.^3 + 2; % and generate some "experimental" data
% create a matrix that contains all you predictors
% note just include the terms you want
vm=[x.^3 ones(size(x))]; % the vandermode matrix
% then solve the least squares problem.
beta = vm\y %this is essentially what polyfit does.
Alternatively download "polyfitn" from the FEX.... it does what it says on the tin AND your problem with a nice interface and some help.

  4 Comments

Show 1 older comment
Sean de Wolski
Sean de Wolski on 6 Mar 2012
http://en.wikipedia.org/wiki/Vandermonde_matrix
You create a vandermonde matrix for your above expression (y = a3.*x^3+a1*x-200) such that v, contains the xpressions to solve for) and v*x=y;
V = [x(:).^3, x(:), ones(size(x(:)))]
Jared
Jared on 6 Mar 2012
Nevermind. I'm dumb and apparently knew what you meant. I just called it a Z matrix all the time instead of vandermode.
Sean de Wolski
Sean de Wolski on 6 Mar 2012
Totally okay. And this question has yielded lots of good edumacational information, so +1.

Sign in to comment.


Jared
Jared on 6 Mar 2012
I am an engineering student learning numerical methods and matlab for the first time. I understand how to do standard polynomial regression however I do not know how to just leave the term out of the model and still solve for the coefficients.
The long way seems to be set the derivatives of the sum of the squares of the residuals to zero with respect to each of the unknown coefficients and I would end up with a column of zeros in my coefficient matrix but I figured there had to be an easier way.

  0 Comments

Sign in to comment.


Richard Willey
Richard Willey on 6 Mar 2012
The 12a release of Statistics Toolbox has some very nice new features for regression analysis. There is a new function named LinearModel for linear regression. The function implements a scripting syntax called "Wilkinson's Notation" which is designed for implementing custom linear models. (Wilkinson's Notation" is pretty standard amongst statisticians)
% Create a data set
X = linspace(1, 10, 100);
X = X';
Y = 3*X.^3 + 5*X + 7;
% Fit a linear regression to the model
myFit = LinearModel.fit(X,Y, 'y~x1:x1:x1 + x1')
plot(myFit)
methods(myFit)
The text string 'y~x1:x1:x1 + x1' translates as y is modeled as function of X^3, X, and a constant. The model excludes the term X^2 which is equivalent to setting the coefficient equal to zero.

  0 Comments

Sign in to comment.


Roland
Roland on 3 Jul 2019
For others looking into polynomial regression with constraints:
A general framework for constraint polynomial approximation (value constraints, derivative constraints, coefficient constraints) can be found in a paper of me and my colleagues:
  • O’Leary, P., Ritt, R. and Harker, M. (2019) ‘Constrained Polynomial Approximation for Inverse Problems in Engineering’, in Wahab, M. A. (ed.) Proceedings of the 1st International Conference on Numerical Modelling in Engineering. Springer Singapore, pp. 225–244. doi: 10.1007/978-981-13-2273-0_19.

  0 Comments

Sign in to comment.


Kin Sung Chan
Kin Sung Chan on 30 Sep 2019
One way is to use the fittype and fix the lower and upper bound. (x,y) are some vectors that you want to fit.
The good thing about this method is you can fit with whatever equation you like.
s = fitoptions('Method','NonlinearLeastSquares',...
'Lower',[10 -10 ],...
'Upper',[10 10 ],...
'MaxIter',10000,...
'TolX',1E-8,...
'Startpoint',[10,0]);
f = fittype('a1*x+b','options',s);
ffit = fit(x, y, f)

  0 Comments

Sign in to comment.


Translated by