Currying and multiple integrals

Question

I am interested in learning an elegant way to use currying in a functional programming language to numerically evaluate multiple integrals. My language of choice is F#.

If I want to integrate f(x,y,z)=8xyz on the region [0,1]x[0,1]x[0,1] I start by writing down a triple integral of the differential form 8xyz dx dy dz. In some sense, this is a function of three ordered arguments: a (float -> float -> float -> float).

I take the first integral and the problem reduces to the double integral of 4xy dx dy on [0,1]x[0,1]. Conceptually, we have curried the function to become a (float -> float -> float).

After the second integral I am left to take the integral of 2x dx, a (float -> float), on the unit interval.

After three integrals I am left with the result, the number 1.0.

Ignoring optimizations of the numeric integration, how could I succinctly execute this? I would like to write something like:

let diffForm = (fun x y z -> 8 * x * y * z)

let result =
    diffForm
    |> Integrate 0.0 1.0
    |> Integrate 0.0 1.0
    |> Integrate 0.0 1.0

Is this doable, if perhaps impractical? I like the idea of how closely this would capture what is going on mathematically.

Answer 1

I'm not entirely sure how you would implement this in a normal way, so this might not fully solve the problem, but here are some ideas.

To do the numerical integration, you'll (I think?) need to call the original function diffForm at various points as specified by the Integrate calls in the pipeline - but you actually need to call it at a product of the ranges - so if I wanted to call it only at the borders, I would still need to call it 2x2x2 times to cover all possible combinations (diffForm 0 0 0, diffForm 0 0 1, diffForm 0 1 0 etc.) and then do some calcualtion on the 8 results you get.

The following sample (at least) shows how to write similar code that calls the specified function with all combinations of the argument values that you specify.

The idea is to use continuations which can be called multiple times (and so when we get a function, we can call it repeatedly at multiple different points).

// Our original function
let diffForm x y z = 8.0 * x * y * z

// At the first step, we just pass the function to a continuation 'k' (once)
let diffFormK k = k diffForm

// This function takes a function that returns function via a continuation
// (like diffFormK) and it fixes the first argument of the function
// to 'lo' and 'hi' and calls its own continuation with both options
let range lo hi func k = 
  // When called for the first time, 'f' will be your 'diffForm'
  // and here we call it twice with 'lo' and 'hi' and pass the 
  // two results (float -> float -> float) to the next in the pipeline
  func (fun f -> k (f lo)) 
  func (fun f -> k (f hi)) 

// At the end, we end up with a function that takes a continuation
// and it calls the continuation with all combinations of results
// (This is where you need to do something tricky to aggregate the results :-))
let integrate result =
  result (printfn "%f")

// Now, we pass our function to 'range' for every argument and
// then pass the result to 'integrate' which just prints all results
let result =
    diffFormK
    |> range 0.0 1.0
    |> range 0.0 1.0
    |> range 0.0 1.0
    |> integrate

This might be pretty confusing (because continuations take a lot of time to get used to), but perhaps you (or someone else here?) can find a way to turn this first attempt into a real numerical integration :-)

Answer 2

I like the idea of how closely this would capture what is going on mathematically.

I'm afraid your premise is false: The pipe operator threads a value through a chain of functions and is closely related to function composition. Integrating over an n-dimensional domain however is analogous to n nested loops, i.e. in your case something like

for x in x_grid_nodes do
    for y in y_grid_nodes do
        for z in z_grid_nodes do
            integral <- integral + ... // details depend on integration scheme

You cannot easily map that to a chain of three independet calls to some Integrate function and thus the composition integrate x1 x2 >> integrate y1 y2 >> integrate z1 z2 is actually not what you do when you integrate f. That is why Tomas' solution—if I understood it correctly (and I am not sure about that...)—essentially evaluates your function on an implicitly defined 3D grid and passes that to the integration function. I suspect that is as close as you can get to your original question.

You did not ask for it, but if you do want to evaluate a n-dimensional integral in practice, look into Monte Carlo integration, which avoids another problem commonly known as the "curse of dimensionality", i.e. that fact that the number of required sample points grows exponentially with n with classic integration schemes.

Update

You can implement iterated integration, but not with a single integrate function, because the type of the function to be integrated is different for each step of the integration (i.e. each step turns an n-ary function to an (n - 1)-ary one):

let f = fun x y z -> 8.0 * x * y * z

// numerically integrate f on [x1, x2]
let trapRule f x1 x2 = (x2 - x1) * (f x1 + f x2) / 2.0 

// uniform step size for simplicity
let h = 0.1

// integrate an unary function f on a given discrete grid
let integrate grid f =
    let mutable integral = 0.0
    for x1, x2 in Seq.zip grid (Seq.skip 1 grid) do
        integral <- integral + trapRule f x1 x2
    integral

// integrate a 3-ary function f with respect to its last argument
let integrate3 lower upper f =
    let grid = seq { lower .. h .. upper }
    fun x y -> integrate grid (f x y)

// integrate a 2-ary function f with respect to its last argument
let integrate2 lower upper f =
    let grid = seq { lower .. h .. upper }
    fun x -> integrate grid (f x)

// integrate an unary function f on [lower, upper]
let integrate1 lower upper f =
    integrate (seq { lower .. h .. upper }) f

With your example function f

f |> integrate3 0.0 1.0 |> integrate2 0.0 1.0 |> integrate1 0.0 1.0

yields 1.0.