We’re going to do a number of different things to optimise the code developed in the previous article, but before we do that, we need to do some baseline performance measurements, so that we know that our “improvements” really are improving things. The performance results we see will also help to guide us a little in the optimisations that we need to do.

Benchmarking is a fundamental activity when working with numerical code, but it can be quite difficult to do it correctly in a lazy language like Haskell. It’s very easy to get into a situation where a benchmarking loop (perhaps to run a computation 1000 times to collect accurate timing information) performs the computation on the first time round the loop, then reuses the computed result for all the other 999 times round.

We can bypass these worries quite effectively using Bryan O’Sullivan’s Criterion library. As well as providing a clean framework for ensuring that pure calculations really do get rerun for timing purposes, Criterion gives us a nice API for setting up benchmarks, running them and collecting results.

Here’s an example of the most basic use of Criterion:

```
module Main where
import Criterion.Main
import Data.Complex
import Data.Vector
import qualified DFT
import qualified FFT
import qualified Numeric.FFT.Vector.Invertible as FFTW
main :: IO ()
main = do
let v = generate 256 (\i -> realToFrac $ sin (2*pi*fromIntegral i / 256)) ::
Vector (Complex Double)
fftwPlan = FFTW.plan FFTW.dft 256
defaultMain [ bench "DFT" $ nf DFT.dft v
, bench "FFT" $ nf FFT.fft v
, bench "FFTW" $ nf (FFTW.execute fftwPlan) v ]
```

Here, we’re comparing the performance for three FFT implementations on a single input vector length. The `defaultMain`

function from Criterion runs a set of benchmarks one after another after collecting some information about the system clock for the machine you’re using. The `bench`

function is a simple way of defining a benchmark for a pure function – you give a name for the benchmark and a function to call. The tricky bit here is the use of the `nf`

function. This takes two arguments, a function applied to all but its last argument, and the last argument, and it ensures that the result of applying the function to the argument is evaluated to normal form (which basically means it gets evaluated as far as possible, so for our FFT algorithms, we perform the whole of the algorithm to get a final result). Note that for the FFTW library, we perform the “plan creation” step outside the benchmark so that we only measure the time taken by the FFT computation itself. Criterion also has facilities for benchmarking monadic code, but we won’t be using those here.

If we compile and run this code, we get output that looks something like this (the numbers and probably also some of the messages will be different on another machine):

```
[seneca:benchmark-demo] $ ./benchmark-demo
warming up
estimating clock resolution...
mean is 1.170747 us (640001 iterations)
found 96871 outliers among 639999 samples (15.1%)
4 (6.3e-4%) low severe
96867 (15.1%) high severe
estimating cost of a clock call...
mean is 31.18633 ns (12 iterations)
found 2 outliers among 12 samples (16.7%)
2 (16.7%) high mild
benchmarking DFT
mean: 21.87155 ms, lb 21.85268 ms, ub 21.90904 ms, ci 0.950
std dev: 131.1261 us, lb 76.22067 us, ub 212.3907 us, ci 0.950
benchmarking FFT
mean: 878.7273 us, lb 876.9702 us, ub 884.2682 us, ci 0.950
std dev: 14.47592 us, lb 5.361621 us, ub 32.03430 us, ci 0.950
found 12 outliers among 100 samples (12.0%)
5 (5.0%) high mild
7 (7.0%) high severe
variance introduced by outliers: 9.424%
variance is slightly inflated by outliers
benchmarking FFTW
mean: 6.682681 us, lb 6.582545 us, ub 6.768512 us, ci 0.950
std dev: 470.4364 ns, lb 421.4242 ns, ub 508.7335 ns, ci 0.950
variance introduced by outliers: 64.627%
variance is severely inflated by outliers
```

Before running any benchmarks, Criterion measures the timing resolution of the system clock, printing some statistical information. It then runs enough repetitions of each of the benchmarks to give a good estimate of the mean time consumed by the computation. The `defaultMain`

function causes some statistics to be printed for each of the benchmarks.

What we really want to do is to collect timing information for each of the three algorithms for a range of input vector lengths. Here’s how we do this:

```
module Main where
import Prelude hiding (enumFromTo, length, sum)
import Control.Monad
import Control.Monad.IO.Class
import Criterion
import Criterion.Config
import Criterion.Monad
import Criterion.Environment
import Data.Complex
import Data.Vector hiding (forM, forM_, map, (++))
import qualified Data.Vector.Unboxed as VU
import qualified DFT
import qualified FFT
import qualified Numeric.FFT.Vector.Invertible as FFTW
import System.IO
tstvec :: Int -> Vector (Complex Double)
tstvec sz = generate sz (\i -> let ii = fromIntegral i
in sin (2*pi*ii/1024) + sin (2*pi*ii/511))
doit :: Environment
-- ^ Criterion timing environment.
-> Int
-- ^ Problem size to benchmark.
-> Criterion (Int, Double, Double, Double)
-- ^ (Problem size, DFT, my FFT, FFTW) result.
doit env sz = do
let v = tstvec sz
dfts <- runBenchmark env $ nf DFT.dft v
myffts <- runBenchmark env $ nf FFT.fft v
let fftwPlan = FFTW.plan FFTW.dft sz
fftws <- runBenchmark env $ nf (FFTW.execute fftwPlan) v
let mean xs = VU.sum xs / fromIntegral (VU.length xs)
return (sz, mean dfts, mean myffts, mean fftws)
main :: IO ()
main = do
hSetBuffering stdout LineBuffering
putStrLn "Size DFT FFT FFTW"
withConfig (defaultConfig { cfgVerbosity = ljust Quiet }) $ do
env <- measureEnvironment
forM_ [8..1024] $ \sz -> do
(_, dft, fft, fftw) <- doit env sz
liftIO $ putStrLn $ show sz ++
" " ++ show (1.0E6 * dft) ++
" " ++ show (1.0E6 * fft) ++
" " ++ show (1.0E6 * fftw)
```

Here, we’re not using Criterion’s `defaultMain`

, but instead we’re setting some configuration options and using the `withConfig`

function to run some Criterion actions. Criterion has a monadic interface (with a monad called, naturally enough, `Criterion`

) that allows us to sequence benchmarking actions ourselves. Before running any benchmarks, we use the `measureEnvironment`

function to measure the system clock resolution. This information is encapsulated in a value of type `Environment`

that we need to pass into the benchmarking functions. We use the usual monadic sequencing operators to loop over the range of input sizes we want to test and, because the `Criterion`

monad is an instance of `MonadIO`

, we can output results from within the `Criterion`

code.

The benchmarks themselves are run by the function `doit`

, which creates a test vector then uses the `runBenchmark`

function to measure the performance of each of the three FFT algorithms in turn. The `runBenchmark`

function takes a timing environment and a `Benchmarkable`

value – pure functions are instances of `Benchmarkable`

, so we can pass a call to `nf`

here, set up appropriately to trigger execution of our FFT functions. Criterion runs each of our functions enough times to get a reasonable idea of the distribution of the time the function takes to complete. The results of the benchmarks are returned as values of type `Sample`

, which is basically just an unboxed vector of double values, from which we can extract statistics for plotting. (For detailed comparison later on, we’ll calculate some uncertainty bounds on these measurements as well, but for a first look, we’re just taking the mean of each sample.)

And here are some results:

This plot shows results (as time for a single transform for a given input vector length) for input vector lengths in the range $[8, 1024]$. Note that the plot is on a log-log scale. The thin lines show $C N \log N$ (black) and $C N^2$ (grey) for a few choices of constant $C$ to give some idea of the scaling of the different results. Let’s look at the DFT results first (the red line). After a slightly faster than $O(N^2)$ rise for small values of $N$, the scaling seems to settle down to something like an $O(N^2)$ increase with input vector length, which is what we expect. Next, the results for FFTW are show in blue. The times here scale more or less as $O(N \log N)$ (and perhaps even a little better than that!) with relatively little variability. It’s definitely the case that some input vector lengths do better and some worse in terms of scaling, but there are no gross discrepancies in performance as $N$ varies.

The results for our mixed-radix Haskell FFT code are shown in green. There are two main things we can draw from this:

The timing results for our FFT code give a sort of “wedge” shape. The bottom edge of the wedge seems to scale as $O(N \log N)$ while the top edge scales as $O(N^2)$, but there are results for different input vector lengths lying all over the space in between these two extremes. The reason for this is fairly obvious: we fall back on the simple $O(N^2)$ DFT algorithm for the prime factors at the “bottom” of our Fourier matrix decomposition. For prime input vector lengths, we just do a normal DFT, which gives us the $O(N^2)$ top edge of our wedge. For cases where we have a “good” prime factorisation, i.e. the last prime factor is small, we get close to $O(N \log N)$ scaling, which gives us the bottom edge of the wedge.

Even in the best cases, our code is something like 100 times slower than FFTW. That sounds bad, I know, but we’re not really comparing like with like – FFTW is probably one of the best optimised pieces of code in existence; our code is totally unoptimised. Of course, I use GHC’s

`-O2`

flag to enable optimisation in the compiler, but we’ve not made any effort to optimise the algorithm itself (apart from the basic Cooley-Tukey divide and conquer approach to the FFT) and, in particular, we’ve not made any effort to look at allocation, unnecessary data copying, inefficiencies due to laziness, and so on.

Most of the rest of this series of articles is going to be about how we can optimise our code. We’re going to take an empirical measurement-driven approach, making modifications to the code based on some obvious first things to try, then on profiling information, and measuring performance at each step along the way to make sure that we’re actually making things better.

What’s a reasonable goal? I don’t yet know. I’d be very happy if we can get within a factor of ten of the performance of FFTW with a pure Haskell implementation.

What should we do first? It seems pretty obvious that those prime factors are something we need to deal with. No matter what other kinds of optimisation we do, as long as we have $O(N^2)$ behaviour in our code, we’re not going to get very far. There are a number of FFT algorithms for prime-length vectors that have $O(N \log N)$ scaling behaviour, and in the next article we’ll implement one of these, after spending a little bit of time thinking about how an efficient prime-length FFT might work.

After that, we’ll take a closer look at our code and do some profiling to get an idea of what we should look at next. Before doing any profiling at all, I can say that there is too much allocation and copying going on, and that we’ll probably end up rewriting things to use a mutable vector to do the Danielson-Lanczos steps in place, but there are probably other things that we can do as well.