How to refactor code to make it functional style?

Question

Playing with F#, I'm trying to think of code in a more functional way. A large part of my work happens to be numerical in nature so I'm thinking whether this re-education makes sense. Is writing numerical code in a functional way like trying to fit a square peg in a round hole or is it simply a matter of a steep learning curve irrespective of the application?

For example, lets take a snippet which demonstrates the weak law of large numbers:

open System
open System.IO
open System.Windows.Forms
open System.Windows.Forms.DataVisualization
open FSharp.Data
open FSharp.Charting
open FSharp.Core.Operators
open MathNet.Numerics
open MathNet.Numerics.LinearAlgebra
open MathNet.Numerics.Random
open MathNet.Numerics.Distributions
open MathNet.Numerics.Statistics


let T = 1000

let arr1 = Array.init T (fun i -> float i*0.)
for i in 0 .. T-1 do
    arr1.[i] <- [|for j in 1..i do yield Exponential.Sample(0.1)|] |> Statistics.Mean

let arr2 = Array.init T (fun i -> float i*0.)
for i in 0 .. T-1 do
    arr2.[i] <- arr1.[1 .. i] |> Statistics.Mean

arr2 |> Chart.Line |> Chart.Show

Is there a succinct functional way of expressing the above? How much of a functional paradigm can be incorporated into work like this?

Answer 1

These are actually two questions: one about improving the given code and one about functional numerical code in F#. Since the other answers already focus on the specific code, I'll focus on the more general question.

Is it about performance?

In my experience, the suitability of functional programming in numerics depends on the performance requirements. The more important execution time is, the more you may want to compromise on functional style.

If performance is not an issue, functional code tends to work very well. It's succinct and safe, and it's closer to mathematical writing than imperative programming. Sure, some problems map very well to imperative programs, but overall, functional style makes a good default choice.

If performance is somewhat an issue, you may want to compromise on immutability. The main cost of functional code in F# comes from the garbage collector, especially from objects with an intermediate lifetime. Making costly objects mutable and re-using them can make a huge difference in execution speed. If you want to write stuff like hydrodynamics, n-body simulations, or games, in a succinct and safe way, but aren't aiming for pedal-to-the-metal execution speed, a multiparadigm F# style can be a good way to go.

If performance is everything, chances are, you want GPU execution anyway. Or maybe make good use of CPU vector units, multithreading, and so on. While there are attempts to use F# on the GPU, the language just isn't designed for speed at all cost. It's probably better to use a language that's closer to hardware in such situations.

When the problem is a mixture of these, it's often possible to mix the solutions. For example, yesterday I needed to do a per-pixel computation on a set of images, and execution time was important. So I read the images in F# using a .NET library, then uploaded them to GPU along with a GLSL compute shader that transformed the pixels, then downloaded them back into "F# land". The point here is that the management operations aren't efficient; the code is still copying stuff around for no real reason. But it was just one operation that would've eaten all the performance, so it's reasonable to use a high-performance tool for that one operation, while all the rest happens neatly and safely in F#.

Answer 2

I'd first not separate the call to Array.init and setting the initial values. You can use the form that @s952163 use in their answer, or based on your code:

let arr1 = Array.init T (fun i ->
    [|for j in 1..i do yield Exponential.Sample 0.1 |] |> Statistics.Mean
)

Issue with this is that you are allocating intermediate arrays, which is costly - and you anyway discard them right after computing the mean. Alternative:

let arr1 = Array.init T (fun i ->
    Exponential.Samples 0.1 |> Seq.take (i+1) |> Seq.average
)

Now for the second part: You are computing the mean of elements 1..i repeatedly, which becomes an O(n^2) operation. You can solve it in O(n) by using the fact that the sum of elements 1..i is the sum of elements 1..{i-1} plus the i.th element.

let sums, _ =
    arr1
    |> Array.mapFold (fun sumSoFar xi ->
        let s = sumSoFar + xi
        s, s
    ) 0.0
let arr2 = 
    sums
    |> Array.mapi (fun i sumi -> sumi / (float (i + 1)))

Of course you can all write that in a single pipe.

Alternatively, use the library function Array.scan to compute the cumulative sums, which would in this case give you a result of length T+1, from which you'd then drop the first element:

let arr2 =
    Array.sub (Array.scan (+) 0.0 arr1) 1 T
    |> Array.mapi (fun i sumi -> sumi / (float (i + 1)))

or avoid the intermediate arrays:

Seq.scan (+) 0.0 arr1
|> Seq.skip 1
|> Seq.mapi (fun i sumi -> sumi / (float (i + 1)))
|> Seq.toArray