Functional programming and root finding
I recently discovered Incanter which looks really promising for statistical computing on the JVM. Incanter is written in Clojure, a lisp like functional programming language for the JVM. We have been using Scala, a hybrid OO/functional programming language for the JVM, in one of our applications but I have yet to find a robust statistics API for Scala. We also use R in the same application. It would be nice to stay within the JVM for statistical procedures rather than communicate between the JVM and an R session. I wanted to investigate Incanter, but first I needed to wrap my head around Clojure. Folding and nesting are common procedures in functional programming languages and in numerical methods.
Many algorithms follow the iterate and accumulate procedure that naturally maps to the folding and nesting paradigms. Newton's Method for polynomial root finding follows this paradigm as do many optimization algorithms that converge on an extrema. As an avid Mathematica enthusiast, I often used NestList to implement and visualize the steps of nesting algorithms, but I haven't found the equivalent built into the libraries of any of the other three languages. So, to dive in to Clojure and Incanter, I decided to implement the root finding in Clojure, Scala, and R for comparison. First, we need a NestList equivalent. Then, we can use the NestList function to implement the root finding algorithm which we'll test out on the trivial polynomial x^2-5 which obviously has its root at ~±2.23606797749979. The root finding algorithm using NestList in Mathematica can be viewed here. In Clojure, the implementation and usage is as follows:
(defn nestlist [fn iv n] (take n (iterate fn iv)))
(defn findroot [f df iv n]
(nestlist #(- %1 (/ (f %1) (df %1))) (double iv) n))user=> (findroot #(- (* %1 %1) 5) #(* 2 %1) 1.0 20)
(1.0 3.0 2.3333333333333335 2.238095238095238 2.2360688956433634 2.236067977499978 2.23606797749979 2.23606797749979 2.23606797749979 2.23606797749979)
The iterate function is exactly what I was looking for. The lazy evaluation of the sequence makes it easy to work with and the definition of the nestlist function is just syntactic sugar on the iterate function. In Scala, the implementation of nestlist uses the foldLeft function and accumulates its results
def nestlist(f:(Double)=>Double, iv: Double, n: Int): List[Double]={
(0 until n).foldLeft(List(f(iv)))
((xs,i) => xs ++ List(f(xs.head)))
}
def findroot(f:(Double)=>Double,
df:(Double)=>Double,iv:Double,n:Int):List[Double]={
nestlist((x)=>x-f(x)/df(x),iv,n)
}
The accumulating list is a bit awkward. I'm sure there are cleaner methods for implementing nestlist in Scala. Usage of the method:
findroot((x)=>x*x-5,(x)=>2*x,1,10)
res1: List[Double] = List(3.0, 2.3333333333333335, 2.238095238095238, 2.2360688956433634, 2.236067977499978, 2.23606797749979, 2.23606797749979, 2.23606797749979, 2.23606797749979, 2.23606797749979, 2.23606797749979)In R, I had to resort to a very non-FP for loop. I looked at the apply family of functions and replicate but couldn't come up with a good algorithm quickly so here is the result.
nestlist=function(f, iv, n) {
acc=as.vector(iv)
for(e in 1:n) acc=append(acc, f(acc[e]))
acc
}
findroot=function(f, df, iv, n)
nestlist(function(x) x-f(x)/df(x), iv, n)And its usage:
> findroot(function(x) x^2-5, function(x) 2*x, 1, 10)
[1] 1.000000 3.000000 2.333333 2.238095 2.236069 2.236068 2.236068 2.236068
[9] 2.236068 2.236068 2.236068
Next steps are to start playing around with Incanter and implement some statistical procedures that utilize the root finding algorithm.