Google Groups no longer supports new Usenet posts or subscriptions. Historical content remains viewable.
Dismiss

Speed up interpolation

25 views
Skip to first unread message

denes.cs...@gmail.com

unread,
Dec 11, 2008, 8:27:26 AM12/11/08
to
Hi Everyone!

I'm writing an app with some math parts, and I wrote the following fn
for Lagrange interpolation. It's the slowest part of the app, used for
locating peak points of curves. I hope the comments make it clear. Any
suggestions how to make it faster?


;;; The following operator creates the Lagrange interpolation
;;; function for a curve defined by an arbitrary number of points.
;;; If we have 3 points, with coordinates x0 y0, x1 y1, x2 y2, then:
;;; poly0(x) = (x-x1)*(x-x2)/(x0-x1)*(x0-x2)
;;; poly1(x) = (x-x0)*(x-x2)/(x1-x0)*(x1-x2)
;;; poly2(x) = (x-x0)*(x-x1)/(x2-x0)*(x2-x1)
;;; lag(x) = y0*poly0(x) + y1*poly1(x) + y2*poly2(x)

;; LAGRANGE-FN returns a function for Lagrange interpolation on the
;; given set of coordinates. It pre-calculates some staying parts of
;; the polynoms and uses some combination to speed up the 'runtime'
;; part.

(defun lagrange-fn (coords)
(let* ((xs (mapcar #'first coords)) ; x coords (wcont)
(permanent (mapcar #'(lambda (coord xx) ; permanent part of
polynoms:
(/ (second coord) ; Ys / denominators
(reduce #'*
(mapcar #'(lambda (es)
(- xx es))
(remove xx xs)))))
coords xs)))
(labels ((given-x? (x coords) ; if x = the x of
one of the
(if (consp coords) ; coords: return the
(if (= x (first (car coords))) ; corresponding y
(second (car coords)) ; (the following algorythm
(given-x? x (cdr coords))) ; assumes that x - xn
nil))) ; is never zero!
(lambda (x)
(let ((given (given-x? x coords)))
(if given
given
(let* ((combinations
(mapcar #'(lambda (xx) ; calculate all x - xn
(- x xx)) xs)) ; combinations once
(their-sum ; and multiply them:
(reduce #'* combinations))) ; ONE division/polynom
(reduce #'+ (mapcar #'(lambda (permanent-part xx)
(* permanent-part
(/ their-sum xx)))
permanent combinations)))))))))


Thanks
d

Rainer Joswig

unread,
Dec 11, 2008, 9:37:33 AM12/11/08
to
In article
<e3a4921b-de93-42ca...@e1g2000pra.googlegroups.com>,
denes.cs...@gmail.com wrote:

One of the problems is that combinations like
reduce mapcar remove are nice to write down,
but come with a cost - especially when the compiler
won't be able to optimize those.

There is an attempt be the SERIES package to
have such computations optimizable in portable
code.

Unfortunately using LOOP will improve the runtime
of the code, too.

So you have a list

a) then you remove things
b) do computation with the remaining objects
c) and then you reduce the result

this means that you traverse/create a list possibly several times.

(loop with result = initial-result-value
for a in list ; map over the list
when (keep-p a) ; filter out
do (setf result (reduction-op (op a) result)
finally (return result))

for example

(loop with result = ()
for a in list
when (evenp a)
do (setf result (* (sqrt a) result))
finally (return result))

for + the LOOP macro supports SUM

Above LOOP will traverse the list only once
and will not generate intermediate lists.

You can create a nice function or macro out of that code
to keep the LOOPs from invading your code.
A macro would ensure 'inlining', but you can also
create a function and inline it - using
the right declarations and using a Lisp compiler that
would do the inlining (there are several).


Once your algorithms are lighter, you can also
use one of the more advanced compilers to produce
optimization hints (CMUCL, SBCL, SCL, LispWorks, ...).
From there you can see where added declaration would
help to create type-specific code.

--
http://lispm.dyndns.org/

budden

unread,
Dec 11, 2008, 2:07:37 PM12/11/08
to
Hi!
Here is some variant. http://paste.lisp.org/display/71993
It removes most cons production. I tested it on two implementations:
Lispworks Personal 4.4.6 for windows and
SBCL 1.0.23.31. It is up to 2 times faster on Lispworks and up to 5
times faster on SBCL. I sometimes also coerce everything to double-
float to avoid unnecessary type conversions. It is ok on Lispworks,
but my function gives other result than yours in SBCL. Didn't
investigate it in-depth.
But my function is still sub-optimal and allocates lots of temporary
storage. To optimize it further, I think typed arrays are suitable. In
this case it will be more like C or FORTRAN with Lisp syntax. Sorry,
got no more time for this now, maybe find some tomorrow.
My function uses iterate I patched to be more convinient,
http://sourceforge.net/projects/iteratekeywords/

Iterate is better than loop, e.g. it has "multiplying" clause which
loop is missing.

----------------
$4/hour lisp freelancer

budden

unread,
Dec 12, 2008, 4:19:10 AM12/12/08
to
Hi!
I tried to improve function further with vectors, but it was a
failure. Sorry, got no more time for this.
------------
$4/hour cl freelancer

denes.cs...@gmail.com

unread,
Dec 12, 2008, 9:26:04 AM12/12/08
to
Thanks for both of you for your comments.

The code I showed in the original post is the 6th version. I tried to
go as far as I can with the list-container approach with mapping
functions etc. to make up for the speed handicap of constant consing.

Unfortunately at work I can only use CLISP. Installing SBCL in itself
would be a big help with this problem. Anyway I tried to turn the
function into a quasi C-style low level code, with arrays and
iteration only.

(defun lagrange-fn6 (coords)
(declaim (optimize (speed 3) (safety 0)))
(let* ((dimension (length coords))
(coordset (make-array (list dimension 2)
:element-type 'float
:initial-contents coords))
(permanent (make-array dimension :element-type 'float))
(combinations (make-array dimension :element-type 'float)))
(do* ((c 0 (1+ c)))
((= c dimension) c)
(setf (aref permanent c)
(do* ((xn 0 (1+ xn))
(acc 1))
((= xn dimension)
(/ (aref coordset c 1) acc))
(if (/= c xn)
(setf acc (* acc
(- (aref coordset c 0)
(aref coordset xn 0))))))))
(lambda (x)
(block lag-fn
(flet ((update-acc (i dif)
(if (zerop dif)
(return-from lag-fn (aref coordset i 1))
(setf (aref combinations i) dif))))
(let ((proto-denominator
(do* ((i 0 (1+ i))
(xn (aref coordset i 0) (aref coordset i 0))
(dif (- x xn) (- x xn))
(acc (update-acc i dif)
(* acc (update-acc i dif))))
((= i (- dimension 1))
acc))))
(do ((i 0 (1+ i))
(acc 0))
((= i dimension)
acc)
(incf acc (* (/ proto-denominator
(aref combinations i))
(aref permanent i))))))))))

It's basically the same algorithm, but the multiple loop variables of
DO allowed some tightening in the computation path. This version is
almost 2 times faster on my home config with SBCL, but the code is not
too lispful. I don't know if further declarations would do any good -
I don't know them very well since CLISP doesn't seem to react to
declaration.

As for ITERATE-KEYWORDS, I couldn't make it work yet, but the speedup
ratio Budden mentioned looks much better than mine. I'll look into it.

thanks
d

0 new messages