From: "John D'Errico" <woodchips@rochester.rr.com>
Path: news.mathworks.com!newsfeed-00.mathworks.com!webx
Newsgroups: comp.soft-sys.matlab
Subject: Re: Constrained least squares/polynomial fitting
Message-ID: <ef29f88.0@webx.raydaftYaTP>
Date: Thu, 23 Feb 2006 19:32:05 -0500
References: <1140728114.186592.85140@j33g2000cwa.googlegroups.com>
Lines: 73
NNTP-Posting-Host: 66.66.30.34
MIME-Version: 1.0
Content-Type: text/plain; charset="ISO-8859-1"
Content-Transfer-Encoding: 8bit
Xref: news.mathworks.com comp.soft-sys.matlab:332785



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