MATLAB Central Newsreader - Constrained least squares/polynomial fitting

Constrained least squares/polynomial fitting

wegwerp@gmail.com — Thu, 23 Feb 2006 12:55:14 -0500

Hi Group,

I am trying to find an offset between two sets of data points. Both
sets of data (which are measurement data, rougly the same size) could
be fitted with a certain polynomial:

order = 2;
p1 = polyfit(x1,y1,order);
p2 = polyfit(x2,y2,order);
offset = polyval(p2,0) - polyval(p1,0); % or offset = p2(end) - p1(end)

The offset I am looking needs to be evaluated at x=0

So far so good. In my application, however, I have the additional
knowledge that both datasets originate from the same function, apart
from the offset. I thus want to do the polyfitting, with the constraint
that the all the higher order terms are equal for both sets, and I am
looking for the difference in the 0-th order coefficient. For example:

p1 = [a2,a1,a0]
p1 = [b2,b1,b0]

constraint: a2 = b2, a1=b1 and I am looking for b0-a0. It would be nice
to know a1 and a2 too, but that is not neccessary. The polynomial order
will probably be between 1 and 4

Any thoughts or any pointers on how to solve this? I know it should be
straightforward to implement this with fminsearch, but since both
datasets could be large this might be really slow. I do have the
gut-feeling that what I am looking for is a linear problem which could
be solved directly. It has been a long time ago that I had those linear
algebra classes ...

Many thanks in advance,
Bas

Re: Constrained least squares/polynomial fitting

John D'Errico — Thu, 23 Feb 2006 19:32:05 -0500

wegwerp wrote:
>
>
> Hi Group,
>
> I am trying to find an offset between two sets of data points. Both
> sets of data (which are measurement data, rougly the same size)
> could
> be fitted with a certain polynomial:
>
> order = 2;
> p1 = polyfit(x1,y1,order);
> p2 = polyfit(x2,y2,order);
> offset = polyval(p2,0) - polyval(p1,0); % or offset = p2(end) -
> p1(end)
>
> The offset I am looking needs to be evaluated at x=0
>
> So far so good. In my application, however, I have the additional
> knowledge that both datasets originate from the same function,
> apart
> from the offset. I thus want to do the polyfitting, with the
> constraint
> that the all the higher order terms are equal for both sets, and I
> am
> looking for the difference in the 0-th order coefficient. For
> example:
>
> p1 = [a2,a1,a0]
> p1 = [b2,b1,b0]
>
> constraint: a2 = b2, a1=b1 and I am looking for b0-a0. It would be
> nice
> to know a1 and a2 too, but that is not neccessary. The polynomial
> order
> will probably be between 1 and 4
>
> Any thoughts or any pointers on how to solve this? I know it should
> be
> straightforward to implement this with fminsearch, but since both
> datasets could be large this might be really slow. I do have the
> gut-feeling that what I am looking for is a linear problem which
> could
> be solved directly. It has been a long time ago that I had those
> linear
> algebra classes ...

Use indicator variables.

Assume that x1 nd x2 are column vectors, as are
y1 and y2. Also, the order of the poly will be p.

n1 = length(x1);
n2 = length(x2);
M = zeros(n1+n2,p+2);
x12 = [x1;x2;];

M(p,:) = x12;
for i = 2:p
M(p-i+1,:) = x12.*M(p-i+2,:);
end
M(p+1,:) = [ones(n1,1);zeros(n2,1)];
M(p+1,:) = [zeros(n1,1);ones(n2,1)];

poly = M\[y1;y2];

The last two coefficients in poly are the constant
coefficients for the two polynomials for sets 1
and 2 respectively, but the rest of the terms are
used in common.

HTH,
John D'Errico

Re: Constrained least squares/polynomial fitting

wegwerp@gmail.com — Fri, 24 Feb 2006 03:55:30 -0500

Hi John,

thanks a lot, that was exactly what I was looking for. There were some
minor typos in your version, so I attached my working version below for
completeness sake.

Cheers,
Bas

p = 4; %order of polynomial
t1 = t1(:); t2 = t2(:); %make column vector
n1=length(t1); n2=length(t2);
t12 = [t1;t2];

M = zeros(n1+n2,p+2);
M(:,p) = t12;
for i = p-1:-1:1
M(:,i) = t12.*M(:,i+1);
end
M(:,p+1) = [ones(n1,1);zeros(n2,1)];
M(:,p+2) = [zeros(n1,1);ones(n2,1)];

poly = M\[y1(:);y2(:)];
poly = poly(:)'; %make row vector
p1 = poly(1:p+1);
p2 = [poly(1:p),poly(p+2)];
offset = poly(p+2)-poly(p+1);

plot(t1,y1,'.',t1,polyval(p1,t1),t2,y2,'.',t2,polyval(p2,t2))