Merging piecewise polynomial (pp) structures

Question

Michal 2023-3-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1934684-merging-piecewise-polynomial-pp-structures

编辑： Bruno Luong 2023-3-30

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bruno Luong 2023-3-28

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1934684-merging-piecewise-polynomial-pp-structures#answer_1202824

编辑：Bruno Luong 2023-3-28

在 MATLAB Online 中打开

Using BSFK from this FEX https://www.mathworks.com/matlabcentral/answers/1810010-free-knot-spline-approximation-bsfk-problem?s_tid=ans_lp_feed_leaf

% 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 个评论
显示 9更早的评论隐藏 9更早的评论

Bruno Luong 2023-3-28

编辑：Bruno Luong 2023-3-29

在 MATLAB Online 中打开

Here is a modified version of the code that recycle PP1 and PP2 to build ppf and ppg rather than fit.

% 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 = [PP2.breaks(1) PP1.breaks(end)]; % [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); % sigmoid (0,1) -> (0,1)
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 overlap interval
n = size(PP1.coefs,2);
arg = {n,[],[],struct('KnotRemoval', 'yes','Animation',0)};
xo = linspace(oI(1),oI(2),100);
ppo = BSFK(xo,h(xo),arg{:});
% Recycle PP1 and PP2 on smaller intervals distinct from overlap interval
ppf = pptruncate(PP1, -Inf, oI(1));
ppg = pptruncate(PP2, oI(2), +Inf);
% Construct ppmerge structure
breaks = [ppf.breaks(1:end-1),ppo.breaks(1:end-1),ppg.breaks(1:end)];
coefs = [ppf.coefs;
         ppo.coefs;
         ppg.coefs];
ppmerge = mkpp(breaks,coefs);
% 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');
%% Truncate a pp on a new interval [left,right]
function pp = pptruncate(pp, left, right)
breaks = pp.breaks;
coefs = pp.coefs;
b = breaks < left;
if b(1)
    n = size(coefs,2);
    i = find(b,1,'last');
    breaks(1:i-1) = [];
    coefs(1:i-1,:) = [];
    if length(breaks) >= 2
        if breaks(2) > left % check if interval not reduced to a point
            % shift the first polynomial to the new most left break point == left
            % using polynomial fit (rather than conv, see John d'Errico's shiftPolySeg())
            x = linspace(left-breaks(1), breaks(2)-breaks(1), n);
            y = polyval(coefs(1,:), x);
            coefs(1,:) = polyfit(x-x(1),y,n-1);
            breaks(1) = left;
        else
            % Interval reduced t a point, remove it
            breaks(1) = [];
            coefs(1,:) = [];
        end
    end
end
b = breaks > right;
if ~isempty(b) && b(end)
    i = find(b,1,'first');
    breaks(i+1:end) = [];
    coefs(i+1:end,:) = [];
    breaks(end) = right;
end
pp = mkpp(breaks,coefs);
end

EDIT: modification to make pptruncate more robust

Michal 2023-3-30

OK, you probably mean the primary sampling of the signals. Am I right?

Bruno Luong 2023-3-30

编辑：Bruno Luong 2023-3-30

在 MATLAB Online 中打开

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...

请先登录，再进行评论。

Answer 2

Matthieu 2023-3-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1934684-merging-piecewise-polynomial-pp-structures#answer_1200294

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

John D'Errico 2023-3-24

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 2023-3-27

编辑：Michal 2023-3-27

@Bruno Luong @John D'Errico

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.

请先登录，再进行评论。

Answer 3

John D'Errico 2023-3-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1934684-merging-piecewise-polynomial-pp-structures#answer_1200399

编辑：John D'Errico 2023-3-24

在 MATLAB Online 中打开

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 个评论
显示 24更早的评论隐藏 24更早的评论

John D'Errico 2023-3-27

在 MATLAB Online 中打开

Sorry. What I wrote works for linear segments. But it needs to be fixed for higher degree segments.

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);
fnplt(PP1)
hold on
fnplt(PP2)
allbreaks = sort([PP1.breaks,PP2.breaks]);
S1 = [0,PP1.coefs];

But S3 needs to be corrected for a higher degree segment. It was fine for a linear segment as I did it.

S3 = [0,shiftPolySeg(PP2.coefs,diff(allbreaks(2:3)))];
S2a = shiftPolySeg(PP1.coefs,diff(allbreaks(1:2)));
S2b = PP2.coefs;

The above is fine. But there was also another problem in my code, necessary when the overlap interval is not exactly ONE unit in length. The linear combination between the two pieces needs to be:

intervallen = diff(allbreaks([2,3]));

S2 = conv([-1/intervallen,1],S2a) + conv([1/intervallen 0],S2b);

PPmerged = mkpp(allbreaks,[S1;S2;S3])

PPmerged = struct with fields:

form: 'pp' breaks: [4 8 10 15] coefs: [3×6 double] pieces: 3 order: 6 dim: 1

fnplt(PPmerged)

function Ps = shiftPolySeg(P,b)

% Given a PP segment on the interval [0,a], find the coefficients of the

% same poly on the interval [b,a], then adjusted so the new polynomial

% lives on the interval [0,a-b].

Pdeg = numel(P)-1;

Ps = [zeros(1,Pdeg),P(end)];

T = 1;

for n = 1:Pdeg

T = conv(T,[1,b]);

Ps = Ps + [zeros(1,Pdeg - n),P(Pdeg - n + 1)*T];

end

Michal 2023-3-28

编辑：Michal 2023-3-28

John OK, the main reason why I asked for this specific algorithm is based on overall processing scenario:

Data are created in the fixed time processing window (typically N ~ 1000 - 10000 equidistant samples length)
The data are smoothed (de-noised) by application of free-knot splines approximation procedure. Most of free-knot spline algorithms has typical numerical complexity O(N^2) or O(N^3) (depends on used optimization algorithm) ... so, for increasing processing window length (N) is free-knot spline algorithm progresivelly more time consuming. So, practically is optimal processing window length N ~ 1000-3000 samples
Typical overlap of the two consecutive processing windows is about 50-300 samples, to be able properly process any irregular behavior of signal near to the processing windows bounds. The transiotion (overlap) region containts, in general, several break-points!!!
Smooted data from the step 2 are represented as PPwindow (piecewise polynomial) structure.
Current PPwindow is then incrementally merged with previous archived PP structure PPmerge = mergepp(PPmerge,PPwindow), where function "mergepp" represents the algorithm I am looking for.
Overlap - Transition region is processed via weighted sum algorithm (preserving continuity or (optionally) smoothness at the end points)
Finally, the structure PPmerge containts smoothed signal for many processing steps in very suitable (compressed) form to very effectivelly evaluate any part of smoothed signal.

This scenario could be possibly processed by Bruno's steps, too. But I am still not sure, how to exactly implement this method. What is your opinion about that?

Bruno Luong 2023-3-28

@Michal reply in my answer

Michal 2023-3-29

@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.

请先登录，再进行评论。

Merging piecewise polynomial (pp) structures

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（3 个）

11 个评论
显示 9更早的评论隐藏 9更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

26 个评论
显示 24更早的评论隐藏 24更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Merging piecewise polynomial (pp) structures

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（3 个）

11 个评论 显示 9更早的评论隐藏 9更早的评论

4 个评论 显示 2更早的评论隐藏 2更早的评论

26 个评论 显示 24更早的评论隐藏 24更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

11 个评论
显示 9更早的评论隐藏 9更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

26 个评论
显示 24更早的评论隐藏 24更早的评论