How to recreate a piecewise polynomial from polyfit to values in each interval?

I need to recreate the piecewise polynomial coefficients from values of the evaluated piecewise polynomial. I want to create an equivalent pp struct.

The pp has been stored as values sampled to enable a unique fit of each polynomial segment (i.e. values at the break times, and 4 internal points for each 6th order polynomial segment: 6 values, 6 polynomial coefficients to solve).

I thought I could use polyfit() on each break interval of samples to get each segment's polynomial coefficients, then use mkpp(). Scaled and centred, or not, the polyfit() coefficients in pp form do not reproduce the original samples, but evaluating polyval(polyfit()) for each segment does. I'm not clear why there's a difference.

What have I missed? Can it be done without the Curve Fitting Tool box?

An example for a single 6th order polynomial, split into 3 segments shows the discrepancy:

% Define polynomial coefs, x and y
p = [1, 2, 3, 4, 5, 6];
x = -5:10;
y = polyval(p, x);

% Define a pp form with knots at -5, 0, 5, 10
d = 1;
knotInterval = 5;
breaks = -5:knotInterval:10;
pieces = length(breaks)-1;

coefs = nan(pieces*d, length(p));
coefsScaled = nan(pieces*d, length(p));
mu = nan(pieces*d, 2);
yPolyval = nan(pieces, knotInterval+1);
% For the 6 (x,y) values in each pp break interval, fit an order 6
% polynomial, centered and scaled polynomial, and evaluate the standard
% fitted polynomial.
for nPiece = 1:pieces
    idxPiece = x >= breaks(nPiece) & x <= breaks(nPiece+1);
    coefs(nPiece, :) = polyfit(x(idxPiece), y(idxPiece), length(p)-1);
    [coefsScaled(nPiece, :), ~, mu(nPiece, :)] = polyfit(x(idxPiece), y(idxPiece), length(p)-1);
    yPolyval(nPiece,:) = polyval(coefs(nPiece, :), x(idxPiece));
end

% Build pp for the coefs and scaled coefs
pp = mkpp(breaks, coefs, d);
ppScaled = mkpp(breaks, coefsScaled, d);

% Plot the original data, and the fitted pp and polynomials
xPolyval = [-5:0; 0:5; 5:10];
plot(x, y, 'x-', x, ppval(pp, x), 'o--',  x, ppval(ppScaled, x), '^--', ...
    xPolyval(1, :), yPolyval(1, :), 'sq-.', xPolyval(2, :), yPolyval(2, :), 'sq-.', ...
    xPolyval(3, :), yPolyval(3, :), 'sq-.')
ylim(1e5*[-1, 3])
legend({'True data', 'pp', 'ppScaled', ...
    'Recovered polyval 1', 'Recovered polyval 2', 'Recovered polyval 3'})

Answer 1

[dpb]

OK, I did have time to go read ppval; indeed, it does presume coefficients are evaluated from each section as the origin, not in absolute coordinates--inside the function is

...
[b,c,l,o]=unmkpp(pp);                     % unpacks pp-->[bkpt,coefs,pieces,order]
[~,index] = histc(xs,[-inf,b(2:l),inf]);  % get which section each input xq goes with
...
% now go to local coordinates ...         % subtracts beginning breakpoint x from input
xs = xs-b(index);                         % vector of x so it thinks each is at origin
...

and goes on to evaluate the functional from there.

To make it work as you've defined the coefficients and coordinates, you have to evaluate over a modified xq vector that counteracts this internal shifting...

[~,ix]=histc(x,[-inf pp.breaks(2:pp.pieces) inf]); % get the equivalent index vector
xq=x+pp.breaks(ix);                                % "real" x plus breakpoint by section

with that then,

plot(x,y,'*-')
hold on
plot(x,ppval(pp,xq),'r-')

returns the expected result

Alternatively, compute the coefficients for each region based on the breakpoint for the region so each section origin is the breakpoint for the section.

Have to read the details of the documentation very carefully; the examples don't show that detail at all clearly. The description for the coefs though does spell it out; of course I didn't read that until just now after all the above! :)

coefs -- Poynomial coefficients

Polynomial coefficients, specified as an L-by-k matrix with the ith row coefs(i,:) containing
the local coefficients of an order k polynomial on the ith interval, [breaks(i), breaks(i+1)].
In other words, the polynomial is
coefs(i,1)*(X-breaks(i))^(k-1) +  coefs(i,2)*(X-breaks(i))^(k-2) + ... + coefs(i,k-1)*(X-breaks(i)) + coefs(i,k).

[Will]

Many thanks for this, it was a bit too subtle for me!

Your suggestion to compute the coefficients in each segment with an x origin at the break point is the best way to go I think, and the unscaled and uncentred coefficients should be used.

I tried simply evaluating at shifted x points as you showed, but this doesn't hold for all cases - like when the polynomial is badly conditioned. This is true for my real case, but fitting the coefficients in a shifted interval gives well conditioned polynomials and gets me within 1e-11 of the true values.

So explicitly, the loop in my example should be:

for nPiece = 1:pieces
    idxPiece = x >= breaks(nPiece) & x <= breaks(nPiece+1);
    coefs(nPiece, :) = polyfit(x(idxPiece)-breaks(nPiece), ...
        y(idxPiece), length(p)-1);
end

[dpb]

Yes, agreed. For higher-order polynomials conditioning is key.

It would take a modified ppval to use the centered form of coefficients in polyfit/val but that wouldn't be too difficult to do; it is actually very straight-ahead code once one looks at it; it clearly hasn't been infected by the later penchant to use all the fancy constructs and internals. :)

It might well be worth looking at it and "rolling your own" version if you're doing this regularly and with large x values in order to avoid numerical issues.

You could still use the structure form; there would be no issues I see with adding another field to the pp structure as flag; there's an internal logic branch now to treat returns from interp1 vis a vis mkpp and friends as far as the output orientation of results if the input xq values are array rather than vector; adding a custom field for 'centered' and the necessary statistics wouldn't interfere with that. You'd probably want to incorporate into a custom mkpp as well for completeness.

While old functions, seems like a worthwhile enhancement request on its own in general, not just for this type of use.