Clear Filters
Clear Filters

Merging piecewise polynomial (pp) structures

10 views (last 30 days)
Michal
Michal on 24 Mar 2023
Edited: Bruno Luong on 30 Mar 2023
Is there any method how to merge two or more "pp" (peicewise polynomial) structures in a case of consecutive "pp" structures with some range overlap to get one "pp" structure describing merged function?

Answers (3)

Bruno Luong
Bruno Luong on 28 Mar 2023
Edited: Bruno Luong on 28 Mar 2023
% Example as in John's answer
breaks = [4 10]; coefs = [ 0 0 1 -2 53];
PP1 = mkpp(breaks,coefs);
breaks = [8 15]; coefs = [ -1 6 1 4 77];
PP2 = mkpp(breaks,coefs);
% overlap interval
oI = [8 10];
li = diff(oI); % its length
f = @(x) ppval(PP1,x);
g = @(x) ppval(PP2,x);
% define smooth transition weight function
sigfun = @(x) x.^2.*(3-2*x);
w2 = @(x) sigfun((x-oI(1))/li);
w1 = @(x) 1-w2(x);
% weighted sum of f and g, applied only in oI
h = @(x) w1(x).*f(x) + w2(x).*g(x);
%% Fit on 3 intervals
arg = {5,[],[],struct('KnotRemoval', 'yes','Animation',0)}; % EDIT first argument 4 => 5
xo = linspace(oI(1),oI(2),100);
ppo = BSFK(xo,h(xo),arg{:});
x1 = linspace(PP1.breaks(1),oI(1));
ppf = BSFK(x1,f(x1),arg{:}); % we can also use PP1 instead
x2 = linspace(oI(end),PP2.breaks(end));
ppg = BSFK(x2,g(x2),arg{:});
% Construct ppmerge structure
ppmerge = ppf;
ppmerge.breaks = [ppf.breaks(1:end-1),ppo.breaks(1:end-1),ppg.breaks(1:end)];
ppmerge.coefs = [ppf.coefs;
ppo.coefs;
ppg.coefs];
ppmerge.pieces = size(ppmerge.coefs,1);
% Graphical check
figure
hold on
xf = linspace(PP1.breaks(1),PP1.breaks(end));
xg = linspace(PP2.breaks(1),PP2.breaks(end));
plot(xf,f(xf),'b+');
plot(xg,g(xg),'rx');
xm = linspace(ppmerge.breaks(1),ppmerge.breaks(end));
plot(xm,ppval(ppmerge,xm),'g','LineWidth',2);
legend('f','g','merge')
A closer look in the overlap interval, it seems the transition is smooth as expected
  11 Comments
Michal
Michal on 30 Mar 2023
OK, you probably mean the primary sampling of the signals. Am I right?
Bruno Luong
Bruno Luong on 30 Mar 2023
Edited: Bruno Luong on 30 Mar 2023
What I means in "subdivisions of f/g subintervals inside overlap region" is that
in the region of [8,10] in the example is the "overlap region".
In this region f has multiple break points, they define a set of "subintervals of f inside overlap region". Each of the subinterval of f is a polynomial of order k-1, since f is pp. So if you select k sample points in each sub intervals (eg. equidistance including the two end break points, thus the word "subdivision") then f is fully characterized by the values computed at this subdivision samples.
That define the sample points for f. Then do the same procedure for g. Take the union of two sets of sampling points, this gives us a way to select the points in the overlap region, and I'm sure it is large enough to recover accurately the approximation of
h = w1*f + w2*g.
on the overlap region by fitting. It is better than selecting emprirical 100 points as I did in my code and it can capture the sharp corner where it needs, etc...

Sign in to comment.


Matthieu
Matthieu on 24 Mar 2023
Are you trying to 'concatenate' your polynomial structures ?
I am not an expert, although I know those are defined by 1 vector and a matrix, breaks & coefs. break defines the range of action of the coefficients of your polynomial, stored in coefs.
You can get those vectors with:
breaks = [0 4 10 15];
coefs = [0 1 -1 1 1; 0 0 1 -2 53; -1 6 1 4 77];
pp = mkpp(breaks,coefs) ;
pp.breaks
ans = 1×4
0 4 10 15
pp.coefs
ans = 3×5
0 1 -1 1 1 0 0 1 -2 53 -1 6 1 4 77
A way of merging those would then be to concatenate the breaks and coefs of all piecewise polynomials you want to merge, as below :
breaks2 = [15 20 22 28];
coefs2 = [0 0 1 2 3; 0 0 0 0 5; -0.1 1 0 0 8] ;
breaks_merged = horzcat(breaks,breaks2(2:end));
coefs_merged = vertcat(coefs,coefs2);
pp_merged = mkpp(breaks_merged,coefs_merged) ;
t = -5:0.1:33 ;
plot(t,ppval(pp_merged,t))
Is this what you needed ?
  4 Comments
John D'Errico
John D'Errico on 24 Mar 2023
Note that a weighted sum will implicitly increase the order of the polynomial segment in the merged region.
A linear weight will result in a function that is one degree higher in the overlapped interval, but it will also result in a derivative discontinuity at the break. So continuous, but not differentiable. So it should also be possible to employ a nonlinear (actually just a higher polynomial order) weighting scheme, if it were important to acheive differentiability on the overlapped interval too.
Michal
Michal on 27 Mar 2023
Edited: Michal on 27 Mar 2023
My main intention is to perform running (moving window) approximation of continuously measured signals by piecewise polynomials. The crucial task is suitable merging separate "pp" structures corresponding to the separate consecutive processing windows with overlap. Transition domain (overlap) between two consecutive windows can be processed by "weighted sum" method.
The final merged function (pp structure) could be one of possible way how to archive all separately measured windows at once in compressed form for next additional processing.
What I am looking for is the method how to perform this task effectively in MATLAB.

Sign in to comment.


John D'Errico
John D'Errico on 24 Mar 2023
Edited: John D'Errico on 24 Mar 2023
Is there any method to "merge" two PP functions? Well, no. Should there be one? I've been using and writing splines tools for 40+ years, in MATLAB for 35 years or so, and answering questions about them for as long. And I've never seen that request, certainly not as you are asking. You don't write code to do something that nobody will ever want to use.
Could there be a way to do so,? Well, yes. The obvious one is to just expand the list of breaks, then append the pp segments as suggested already.
However, your question is even less, let me say, "expected", because you have an overlap. What would you do in the overlap region? Arbitrarily choose one function of the other? Take the mean? Flip a coin to decide?
Sigh. COULD you do something? Well, I suppose you could do something. I could think of several choices even here. The problem is in the overlap region, you need to resolve the conflct. That is, suppose we have intervals at the ends of domain 1 as [A,B], and at the beginning of domain 2 as [C,D]. I'll assume from what you have said is that A < C < B < D. Otherwise there is no problem. Actually, I can see another issue. What if B < C? That is, the problem is simple, if B == C. Now the direct merger already suggested works. But if B < C, then there is an undefined region. I suppose you could extrapolate the curves, over the hole in the middle between domains. But extrapolating splines is about the most obscene thing you can do. Avoid it at all costs.
Anyway, suppose we have the sceneario where A < C < B < D? We could now adjust the breaks to be exactly that sequence. Where on the interval [A,B), we use PP1 from the lower domain. On the interval [B,D), we use PP2, so the upper segment.
But on the interval [C,B) What choice do we have? Perhaps we might just use the average of the two functions on that interval. But that would now almost certainly introduce a point of discontinuity at the breaks, a bad thing.
One idea might be to use an average that varies along the interval. Essentially, this would force the curve to be a higher order segment. So instead of a cubic segment there, it would now be, a degree 4 segment at least. That would allow the resulting function to be continuous, though not necessarily differentiable. If we are willing to make it a 5th degree segment, that should allow one to construct a "merged" segment that would be both continous AND differentiable. (Just thinking off the cuff there.)
So, doable. But, for example, what would you do here?
PP1 = mkpp([0 2],[1 0]);
PP2 = mkpp([1 3],[-1 2]);
fnplt(PP1)
hold on
fnplt(PP2)
hold off
Lets see. The new pp form would have 4 breaks.
allbreaks = sort([PP1.breaks,PP2.breaks]);
S1 = [0,PP1.coefs];
S3 = [0,PP2.coefs(1),fnval(PP2,allbreaks(3))];
S2a = [PP1.coefs(1),fnval(PP1,allbreaks(2))];
S2b = PP2.coefs;
S2 = conv([-1,1],S2a) + conv([1 0],S2b);
PPmerged = mkpp(allbreaks,[S1;S2;S3]);
fnplt(PPmerged,'g')
hold on
fnplt(PP1,'r')
fnplt(PP2,'b')
hold off
So the "merged function is the same as the old ones, in the intervals [0,1], and [2,3]. nd it is made to be continuous at the overlap points, but differentiability was not enforced. On the sub-interval [1 2], the segment is now quadratic, so one degree higher order. That was necessary to enforce continuity.
As I said, doable.
  26 Comments
Michal
Michal on 29 Mar 2023
@John D'Errico John, any progress with your method based on "shiftPolySeg" function for multi-break overlap region? Is my purpose description clear enough for you? From my point of view, is the effective and robust overlapped PP structures merging method very important for anybody who needs to apply incremental version of piecewise polynomial approximation on running data stream.
I still think, that your approach could be more general and robust than Bruno's method based on re-fit of overlap region. The Bruno's re-fit based method requires additional user defined control parameter (number of weighted sum samples ~ 100), which is not simple to "optimally" estimate. Moreover, the re-fit introduce additional computing overhead and needs to use proper fitting algorithm. Your method works only with PP structures, which is an interesting way.

Sign in to comment.

Categories

Find more on Polynomials in Help Center and File Exchange

Products


Release

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!