Newton-Raphson Made Easy with Scala

Introduction

Some years ago I received a brochure for a very expensive software that proudly announced that it could perform financial calculus with the Newton-Raphson method. Wow, sounds good, but was it really something so special? Well, Isaac Newton and Joseph Raphson invented a rather simple but brilliant idea for handling certain mathematical problems. On the other hand, its implementation is software is rather straightforward. In fact, with functional programming in a language like Scala it can be remarkably easy to achieve.

The Newton-Raphson method explained

Ever since primary school, we were told that an equation is something that can be resolved by analysis. For instance, if x2 - x = 1, we know that x must be either 1 or -1. Or that if ln(x) - 1 = 0, x equals the number e.  But occasionally it can be difficult or even impossible to resolve an equation analytically. Some of these cases are related to normal distributions, the Black-Scholes formula or yield curves, all of which are very important in financial calculus.

An alternative to analytical solution is to make some kind of trial and error approximation. This is exactly what Newton-Raphson is all about and it is based on a simple principle. Let’s say we want to find a value of some variable x that fulfils the condition y = f(x) = 0. We start with a trial number which we call the seed value. If it overshoots, so that f(x_k) > 0, and the function is an upward-sloping curve at x_k, we know that we should try with a smaller x_k+1. Similar statements can be made for the cases when f(x_k) < 0 or the curve Is downward-sloping curve at x_k. The formal definition of the Newton-Raphson method is as follows:

x_k+1 = x_k – f(x_k) / f’(x_k)

y = f(x_k) + f’(x_k) * (x_k+1 - x_k)

If y is zero, we have found our value of x. Otherwise we repeat the algorithm with x_k+1 as the new trial value.

A concise implementation in Scala

The algorithm above can easily be represented in a Scala program. To cover a more general scope, we’ll make use the possibility to pass functions as variables to other functions. To start with, we define an equation y = f(x) as a Scala method func that takes a double precision value as argument and returns the value of f at x:

def func= (Double) => { … }

We also need to define a method deriv, which also takes a double precision value as argument but returns the value of the first order derivative f’ at x:

def deriv = (Double) => { … }

We can now create a concise iterative method that takes care of the rest:

def newtonRaphson(f: Double => Double, d: Double => Double, xk: Double): Double = {
  val x = xk - f(xk)/d(xk)
  val y = f(x) + d(x) * (x - xk)
  if (y == 0) xk // Resolved
  else newtonRphson(f, d, x) // Not resolved, run next iteration
}

To run the algorithm, we just have to call the method with a fairly good guess on the value of x for the initial trial:

val seedValue = 3.0 // Important: double precision number format
val approximation = newton(func, deriv, seedValue)

It’s noteworthy that our method will work for any reasonable definition of func and deriv.

Multivariable too, with a little help from Breeze

We might use the Newton-Raphson method for a system of n functions of n variables as well, by expressing the conditions in matrix form:

X_k+1 = X_k – the inverse of J_k * F(X_k)

Y = F(X_k ) + J_k * (X_k+1 – X_k)

where X, Y and F(X) are vectors, for example (₁x, ₂x, … , _nx)., and J_k is the Jacobian matrix of first order derivatives, evaluated at X_k.

This can be defined in a Scala program by including the object DenseMatrix from the ScalaNLP Breeze library, which allows for some convenient matrix algebra. To start with, we define the array of functions and the Jacobian matrix as:

def func = (m: DenseMatrix[Double]) => DenseMatrix(f1(m11, m12, …), f2(), … )

def jacob = (m:  DenseMatrix[Double]) => DenseMatrix(m1(m11, m12, …), f2(), … )

The respective sizes of these should of course be n * 1 for func and n * n for jacob. The next step is to define the iterative method:

def newtonRaphson (f: DenseMatrix[Double] => DenseMatrix[Double], j: DenseMatrix[Double] => DenseMatrix[Double], xk: DenseMatrix[Double]): DenseMatrix[Double] = {
  val x = xk - inv(j(xk)) * f(xk)
  val y = f(x) + j(x) * (x - xk)
  if (y) == DenseMatrix.zeros[Double](y.rows, y.cols)) x //Resolved, we return the array of x
  else newtonRaphson(f, j, x) // Not resolved, run next iteration
}

As before, we pass functions of our choice to the Newton-Raphson method:

val seedValues = DenseMatrix(3.0, … )  // Important: double precision number format
val approximation = newtonRaphson (func, jacon, seedValues)

More information

A code example based on this article is available at Github.