Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
find roots of a spline

Subject: find roots of a spline

From: Andres

Date: 22 Apr, 2010 12:49:08

Message: 1 of 4

Hi,

I'd like to speed up my code for the following problem:

input: two vectors x, y (numel(x) = numel(y) = ~5e5)
output: roots of the cubic spline through (x,y)

I've written a function 'rootspline' to get the roots of the piecewise polynomial form of the cubic spline 'pp' (see below), so I run

pp = spline(x,y);
xo = rootspline(pp);

compare example 3 in the rootspline doc.
Any ideas how to improve rootspline, its subfunction rootsbnd, or how to do it in a different way?
Looking at the profiler results, it seems that calling rootsbnd for each piece of the spline hurts most...

Thanks for your thoughts!
Andres

>>
function zo = rootspline(pp)

% rootspline find roots of a spline
%
% xo = rootspline(pp) returns the unique roots of the piecewise polynomial
% contained in pp. Typically, pp was constructed using the functions
% interp1, pchip, spline, or the spline utility mkpp.
%
% rootspline uses the roots function, so the input must not contain NaN or
% Inf.
%
% EXAMPLES:
%
% % choose source data:
% % ----------------- 1 -----------------
% x = -2:2;
% y = [ 4 1 0 -1 -4 ];
% % ----------------- 2 -----------------
% x = [-2,-1,1,2];
% y = [ 4 1 1 4];
% % ----------------- 3 -----------------
% x = unique(rand(21,1)*20-10);
% y = rand(size(x))-0.5;
% % -------------------------------------
%
% % calculate spline
% pp = spline(x,y);
% xo = rootspline(pp);
% % plot
% xx = linspace(min(x),max(x),50*numel(x));
% figure
% plot(xx,ppval(pp,xx))
% hold on, grid on
% plot(x,y,'k*')
% plot(xo,0*xo,'ro')
% legend({'spline','source data','roots'},'location','best')
% xlabel('x'), ylabel('y'), title('rootspline example')
%
% Compare examples 1), 2) and note the roundoff susceptibility!
%
%
% See also
% roots, spline, interp1, pchip, mkpp

x = pp.breaks(:);
y = [pp.coefs(:,end);ppval(pp,x(end))];


% preallocate
zo = zeros(pp.pieces,1);
% initialize roots counter
cnt = 0;
% step through all pieces of the spline and find roots
for idx = 1:pp.pieces
    r = rootsbnd(pp.coefs(idx,:),pp.breaks([idx,idx+1]));
    zo(cnt+1:cnt+numel(r)) = r;
    cnt = cnt + numel(r);
end
zo(cnt+1:end) = [];

% get unique roots (zo is already sorted)
if numel(zo)>1
    zo = zo([true;diff(zo)~=0]);
end

end

function r = rootsbnd(p,bnd)

% rootsbnd find roots of a polynomial piece inside an interval
%
% r = rootsbnd(p,bnd) returns the roots of the polynomial represented by
% the coefficient vector p (e.g. p = pp.coefs, pp produced by spline)
% inside the interval bnd = [x1,x2].

r = roots(p);
% sort out non-real and outlying roots
idc = imag(r)==0;
idc = idc & (r>=0 & r<=diff(bnd));
r = r(idc);
% sort multiple roots
if numel(r)>1
    r = sort(r);
end
% shift to interval
r = r + bnd(1);

end
<<

Subject: find roots of a spline

From: John D'Errico

Date: 22 Apr, 2010 13:27:21

Message: 2 of 4

"Andres" <rantore@werb.deNoRs> wrote in message <hqpgk4$l4q$1@fred.mathworks.com>...
> Hi,
>
> I'd like to speed up my code for the following problem:
>
> input: two vectors x, y (numel(x) = numel(y) = ~5e5)
> output: roots of the cubic spline through (x,y)
>
> I've written a function 'rootspline' to get the roots of the piecewise polynomial form of the cubic spline 'pp' (see below), so I run
>
> pp = spline(x,y);
> xo = rootspline(pp);

Some years ago, I wrote an inverse tool for cubic splines.

Rather than using roots, use the explicit formula for the
roots of a cubic polynomial, given the coefficients. This
can be done in a fully vectorized form, so all the intervals
can be processed at once.

Other tricks you can do are to break up knot intervals
into intervals such that the spline is strictly monotone
over each interval. That will allow you to only apply
a roots operation to the intervals where you KNOW a
root exists already.

John

Subject: find roots of a spline

From: Andres

Date: 22 Apr, 2010 14:49:21

Message: 3 of 4

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <hqpirp$11p$1@fred.mathworks.com>...

> Rather than using roots, use the explicit formula for the
> roots of a cubic polynomial, given the coefficients. This
> can be done in a fully vectorized form, so all the intervals
> can be processed at once.
>
> Other tricks you can do are to break up knot intervals
> into intervals such that the spline is strictly monotone
> over each interval. That will allow you to only apply
> a roots operation to the intervals where you KNOW a
> root exists already.
>
> John

Good ideas. One could rely on roots.m if the order of the spline is > 3 (or even 4), and use the explicit formulae otherwise (if the spline is not necessarily cubic).
Of course some special cases (e.g. y=0 inside an interval of x) need to be considered for a general solution as well.
Thanks!

Subject: find roots of a spline

From: John D'Errico

Date: 22 Apr, 2010 16:24:20

Message: 4 of 4

"Andres" <rantore@werb.deNoRs> wrote in message <hqpnlh$605$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <hqpirp$11p$1@fred.mathworks.com>...
>
> > Rather than using roots, use the explicit formula for the
> > roots of a cubic polynomial, given the coefficients. This
> > can be done in a fully vectorized form, so all the intervals
> > can be processed at once.
> >
> > Other tricks you can do are to break up knot intervals
> > into intervals such that the spline is strictly monotone
> > over each interval. That will allow you to only apply
> > a roots operation to the intervals where you KNOW a
> > root exists already.
> >
> > John
>
> Good ideas. One could rely on roots.m if the order of the spline is > 3 (or even 4), and use the explicit formulae otherwise (if the spline is not necessarily cubic).
> Of course some special cases (e.g. y=0 inside an interval of x) need to be considered for a general solution as well.
> Thanks!

You need to check the results carefully in any such
computation.

For example, what if the computation yields a slightly
complex result?

What if the root found is just eps outside the interval
of interest?

What if the function is not truly cubic over the interval?
That is, suppose the cubic coefficient is insignificant
on a given interval? It will generate a spurious root,
but where?

The most careful code will worry about all of these
things.

John

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us