Crushing ISAAC

(UPDATE 2020/06/30: the good people at tweag.io have since published a Nix shell environment that appears to make testing arbitrary PRNGs much less of a pain. I recommend you check it out!)

I recently needed a good cryptographically-secure and seedable pseudorandom number generator for Javascript. This didn’t turn out to be as trivial a procedure as I figured it’d be: most Javascript CSPRNGs I found didn’t appear to be manually seedable, instead automatically seeding themselves behind-the-scenes using a trusted high-quality entropy source like /dev/urandom (as per the WebCryptoAPI spec).

But! I want to use my own entropy source. So I could either implement my own cryptographically-secure PRNG, which I would obviously need to test rigorously, or I could find an existing implementation that had already been vetted widely by use. I settled on the ISAAC PRNG/stream cipher, both because it seemed like a reasonable choice of PRNG (it is used in things like coreutils’ shred, and there are no known practical attacks against it), and also because there was a Javascript implementation floating around on the internet. But the interesting thing was that the implementation did not really seem to be vetted widely – it hadn’t been updated in five years or so, and both the Github repository and the npm package seem to show very little activity around it in general.

(Update: Stefano, the author of the fork I used, later emailed me to point out that the original version of this ISAAC code is located here. His fork, on the other hand, is node-friendly.)

I was going to be using this thing in an application that made some nontrivial demands in terms of security. So while the PRNG itself seemed to fit the bill, I wanted to be assured that the implementation satisfied at least some basic criteria for pseudorandomness. I assumed that it probably works – I just wanted to have some reasonable confidence in making that assumption.

There are a number of statistical suites for testing PRNGs out there: I am aware of at least the NIST statistical test suite, the DieHard suite, and TestU01, being most familiar with the latter (this is not saying much). The TestU01 suite is implemented in C; you provide it with a PRNG represented as a function pointer that returns 32-bit integers, and it will run the desired battery of tests (SmallCrush, Crush, or BigCrush – each of increasing size and scope) for you. These more or less consist of frequentist tests involving the chi-square distribution, with more specialised test statistics appearing on occasion. Results are provided in terms of p-values on these statistics.

I found this page that gave some instructions for using TestU01, but the ISAAC implementation I wanted to test is of course written in Javascript. So I knew there was some FFI hackery to be done in order to get the two codebases to play nicely together. I also discovered that Jan de Mooij did some work on testing JS’s basic ‘Math.random’ generator with TestU01 using Emscripten, an LLVM-to-JS compiler, so these two resources seemed a useful place to start.

After several hours of bashing my way through emcc compilation and linker errors careful and methodical programming, I managed to get everything to work. Since the documentation for these tools can be kind of sparse, I figured I’d outline the steps to run the tests here, and hopefully save somebody else a little bit of time and effort in the future. But that said, this is inevitably a somewhat tricky procedure – when trying to reproduce my steps, I ran into a few strange errors that required additional manual fiddling. The following should be about right, but may require some tweaking on your end.

Emscripten and TestU01

Install and activate Emscripten in the current terminal. To go from the official git repo:

$ git clone https://github.com/juj/emsdk.git
$ cd emsdk
$ ./emsdk install latest
$ ./emsdk activate latest
$ source ./emsdk_env.sh

Emscripten provides a couple of wrappers over existing tools; notably, there’s emconfigure and emmake for specialising make builds for compilation via emcc, the Emscripten compiler itself.

In some other directory, grab the TestU01 suite from the Université de Montréal website and extract it:

$ wget http://simul.iro.umontreal.ca/testu01/TestU01.zip
$ unzip -q TestU01.zip

This is some oldish, gnarly academic C code that uses a very wonky autoconf/automake-based build system. There is probably a better way to do it, but the easiest way to get this thing built without too much grief is to build it twice – once as per normal, specifying the appropriate base directory, and once again to specialise it for Emscripten’s use:

$ basedir=`pwd`
$ cd TestU01-1.2.3
$ ./configure --prefix="$basedir"
$ make -j 2
$ make -j 2 install

If all goes well you’ll see ‘bin’, ‘include’, ‘lib’, and ‘share’ directories pop up in ‘basedir’. Repeat the analogous steps for emscripten; note that you’ll get some “no input files” errors here when the configure script checks dynamic linker characteristics, but these are inconsequential:

$ emconfigure ./configure --prefix="$basedir"
$ emmake make -j 2
$ emmake make -j 2 install

Similarly, you’ll notice some warnings re: dynamic linking when building. We’ll handle these later. In the meantime, you can return to your ‘basedir’ to continue working:

$ cd $basedir

A Test Shim for ISAAC

Check out the raw ISAAC code. It’s structured in a sort-of-object-y way; the state of the PRNG is held in a bunch of opaque internal variables, and the whole thing is initialised by calling the ‘isaac’ function and binding the result as a variable. One then iterates the PRNG by calling either the ‘rand’ or ‘random’ property of that variable for a random integer or double, respectively.

We need to write the actual testing code in C. You can get away with the following, which I’ve adapted from M.E. O’Neill’s code – call it something like ‘test-isaac.c’:

#include <emscripten.h>
#include "TestU01.h"

extern void isaac_init(void);
extern unsigned int isaac_rand(void);

int main()
{
    isaac_init();

    unif01_Gen* gen = unif01_CreateExternGenBits("ISAAC", isaac_rand);

    bbattery_SmallCrush(gen);

    unif01_DeleteExternGenBits(gen);

    return 0;
}

Note the two external functions I’m declaring here: the first mimics calling the ‘isaac’ function in Javascript and binding it to a variable, ‘isaac’, and the second mimics a call to isaac.rand(). The testing code follows the same pattern: isaac_init() initialises the generator state, and isaac_rand produces a value from it. The surrounding code passes isaac_rand in as the generator to use for the SmallCrush battery of tests.

C to LLVM Bitcode

We can compile this to LLVM IR as it is, via Emscripten. But first, recall those dynamic linker warnings from the initial setup step. Emscripten deals with a lot of files, compile targets, etc. based on file extension. We thus need to rename all the .dylib files in the ‘lib’ directory, which are in fact all LLVM bitcode, to have the appropraite .bc extension. From the ‘lib’ directory itself, one can just do the following in bash:

$ for old in *.dylib; do mv $old `basename $old .dylib`.bc; done

When that’s done, we can compile the C code to LLVM with emcc:

$ emcc -O3 -o test-isaac.bc \
    -I/$basedir/include \
    -L/$basedir/lib \
    -ltestu01 -lprobdist -lmylib -lm -I/usr/local/include \
    test-isaac.c

Again, Emscripten decides what its compile target should be via its file extension – thus here, an output with the .bc extension means we’re compiling to LLVM IR.

LLVM to Javascript and WASM

Now, to provide the requisite isaac_init and isaac_rand symbols to the compiled LLVM bitcode, we need to pass the ISAAC library itself. This is another finicky procedure, but there is a method to the madness, and one can find a bit of documentation on it here.

Helpfully, Evan Wallace at Figma wrote an emscripten JS library generation helper that makes this task a little less painful. Install that via npm as follows:

$ npm install -g emscripten-library-generator

To wrap the ISAAC code up in the appropriate format, one needs to make a few small modifications to it. I won’t elaborate on this too much, but one needs to:

• Change the String.prototype.toIntArray function declaration to a simple function toIntArray(string) declaration, and alter its use in the code appropriately,
• Change the var isaac = ... function declaration/execution binding to a simple function isaac() declaration, and,

• Declare the two functions used in our C shim:

  function isaac_init() {
    var isaac_initialised = isaac();
  }
  
  function isaac_rand() {
    return isaac_initialised.rand();
  }
  

Then we can package it up in an emscripten-friendly format as follows:

$ emscripten-library-generator isaac-mod.js > isaac-lib.js

You’ll need to make one last adjustment. In the isaac-lib.js file just generated for us, add the following emscripten ‘postset’ instruction above the ‘isaac’ property:

isaac__postset: 'var isaac_initialised = _isaac();',   // add me
isaac: function () {                                   // leave me alone

This makes sure that the isaac_initialised variable referred to in isaac_rand is accessible.

Whew. Ok, with all that done, we’ll compile our LLVM bytecode to a Javascript and wasm pair. You’ll need to bump up the total memory option in order to run the resulting code – I think I grabbed the amount I used from Jan de Mooij’s example:

$ emcc --js-library isaac-lib.js \
    lib/libtestu01.0.0.1.bc \
    -o test-isaac.js \
    -s TOTAL_MEMORY=536870912 \
    test-isaac.bc

Running SmallCrush

That’s about it. If all has gone well, you should have seen no compilation or linking errors when running emcc, and you should have a couple of ‘test-isaac.js’ and ‘test-isaac.wasm’ files kicking around in your ‘basedir’.

To (finally) run the Small Crush suite, execute ‘test-isaac.js’ with node:

$ node test-isaac.js
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
                 Starting SmallCrush
                 Version: TestU01 1.2.3
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx


***********************************************************
HOST = emscripten, Emscripten

ISAAC


smarsa_BirthdaySpacings test:
-----------------------------------------------
   N =  1,  n = 5000000,  r =  0,    d = 1073741824,    t = 2,    p = 1


      Number of cells = d^t = 1152921504606846976
      Lambda = Poisson mean =      27.1051


----------------------------------------------------
Total expected number = N*Lambda      :      27.11
Total observed number                 :      27
p-value of test                       :    0.53


-----------------------------------------------
CPU time used                    :  00:00:00.00

Generator state:

<more output omitted>

========= Summary results of SmallCrush =========

 Version:          TestU01 1.2.3
 Generator:        ISAAC
 Number of statistics:  15
 Total CPU time:   00:00:00.00

 All tests were passed

Et voilà, your SmallCrush battery output is dumped to stdout. You need only to tweak the C code shim and recompile if you want to run the more intensive Crush or BigCrush suites. Similarly, you can dump generator output to stdout with console.log() if you want to reassure yourself that the generator is running correctly.

Fin

So: the Javascript ISAAC PRNG indeed passes TestU01. Nice! It was satisfying enough to get this hacky sequence of steps to actually run, but it was even better to see the tests actually pass, as I’d both hoped and expected.

I did a few extra things to convince myself that ISAAC was really passing the tests as it seemed to be doing. I ran TestU01 on a cheap little xorshift generator (which failed several tests) and also did some ad-hoc analysis of values that I had ISAAC log to stdout. And, I even looked at the code, and compared it with a reference implementation by eye. At this point, I’m convinced it operates as advertised, and I feel very comfortable using it in my application.

Transforming to CPS

I recently picked up Appel’s classic Compiling with Continuations and have been refreshing my continuation-fu more generally.

Continuation-passing style (CPS) itself is nothing uncommon to the functional programmer; it simply involves writing in a manner such that functions never return, instead passing control over to something else (a continuation) to finish the job. The simplest example is just the identity function, which in CPS looks like this:

id :: a -> (a -> b) -> b
id x k = k x

The first argument is the conventional identity function argument – the second is the continuation. I wrote a little about continuations in the context of the Giry monad, which is a somewhat unfamiliar setting, but one that follows the same principles as anything else.

In this post I just want to summarise a few useful CPS transforms and related techniques in one place.

Manual CPS Transformation

Consider a binary tree type. We’ll keep things simple here:

data Tree a =
    Leaf a
  | Branch a (Tree a) (Tree a)

Calculating the depth of a tree is done very easily:

depth :: Tree a -> Int
depth = loop where
  loop tree = case tree of
    Leaf _       -> 1
    Branch _ l r ->
      let dl = loop l
          dr = loop r
      in  succ (max dl dr)

Note however that this is not a tail-recursive function – that is, it does not end with a call to itself (instead it ends with a call to something like ‘succ . uncurry max’). This isn’t necessarily a big deal – the function is easy to read and write and everything, and certainly has fine performance characteristics in Haskell – but it is less easy to deal with for, say, an optimising compiler that may want to handle evaluation in this or that alternative way (primarily related to memory management).

One can construct a tail-recursive (depth-first) version of ‘depth’ via a manual CPS transformation. The looping function is simply augmented to take an additional continuation argument, like so:

depth :: Tree a -> Int
depth tree = loop tree id where
  loop cons k = case cons of
    Leaf _       -> k 1
    Branch _ l r ->
      loop l $ \dl ->
        loop r $ \dr ->
          k (succ (max dl dr))

Notice now that the ‘loop’ function terminates with a call to itself (or just passes control to a supplied continuation), and is thus tail-recursive.

Due to the presence of the continuation argument, ‘loop’ is a higher-order function. This is fine and dandy in Haskell, but there is a neat technique called defunctionalisation that allows us to avoid the jump to higher-order and makes sure things stay KILO (“keep it lower order”), which can be simpler to deal with more generally.

The idea is just to reify the continuations as abstract syntax, and then evaluate them as one would any embedded language. Note the continuation \dl -> .., for example – the free parameters ‘r’ and ‘k’ occuring in the function body correspond to a tree (the right subtree) and another continuation, respectively. And in \dr -> .. one has the free parameters ‘dl’ and ‘k’ – now the depth of the left subtree, and the other continuation again. We also have ‘id’ used on the initial call to ‘loop’. These can all be reified via the following data type:

data DCont a =
    DContL (Tree a) (DCont a)
  | DContR Int (DCont a)
  | DContId

Note that this is a very simple recursive type – it has a simple list-like pattern of recursion, in which each ‘level’ of a value is either a constructor, carrying both a field of some type and a recursive point, or is the ‘DContId’ constructor, which simply terminates the recursion. The reified continuations are, on a suitable level of abstraction, more or less the sequential operations to be performed in the computation. In other words: by reifying the continuations, we also reify the stack of the computation.

Now ‘depth’ can be rewritten such that its looping function is not higher-order; the cost is that another function is needed, one that lets us evaluate items (again, reified continuations) on the stack:

depth :: Tree a -> Int
depth tree = loop tree DContId where
  loop cons k = case cons of
    Leaf _       -> eval k 1
    Branch _ l r -> loop l (DContL r k)

  eval cons d = case cons of
    DContL r k  -> loop r (DContR d k)
    DContR dl k -> eval k (succ (max dl d))
    DContId     -> d

The resulting function is mutually tail-recursive in terms of both ‘loop’ and ‘eval’, neither of which are higher-order.

One can do a little better in this particular case and reify the stack using an actual Haskell list, which simplifies evaluation somewhat – it just requires that the list elements have a type along the lines of ‘(Tree a, Int)’ rather than something like ‘Either (Tree a) Int’, which is effectively what we get from ‘DCont a’. You can see an example of this in this StackOverflow answer by Chris Taylor.

Mechanical CPS Transformation

“Mechanical CPS transformation” might be translated as simply “compiling with continuations.” Matt Might has quite a few posts on this topic; in particular he has one very nice post on mechanical CPS conversion that summarises various transformations described in Appel, etc.

Matt describes three transformations that I think illustrate the general mechanical CPS business very well (he describes more, but they are more specialised). The first is a “naive” transformation, which is simple, but produces a lot of noisy “administrative redexes” that must be cleaned up in another pass. The second is a higher-order transformation, which makes use of the host language’s facilities for function definition and application – it produces simpler code, but some unnecessary noise still leaks through. The last is a “hybrid” transformation, which makes use of both the naive and higher-order transformations, depending on which is more appropriate.

Let’s take a look at these in Haskell. First let’s get some imports out of the way:

{-# LANGUAGE OverloadedStrings #-}

import Data.Monoid
import Data.Text (Text)
import qualified Data.Text as T
import Data.Unique
import qualified Text.PrettyPrint.Leijen.Text as PP

I’ll also make use of a simple, Racket-like ‘gensym’ function:

gensym :: IO Text
gensym = fmap render newUnique where
  render u =
    let hu = hashUnique u
    in  T.pack ("$v" <> show hu)

We’ll use a bare-bones lambda calculus as our input language. Many examples – Appel’s especially – use significantly more complex languages when illustrating CPS transforms, but I think this distracts from the meat of the topic. Lambda does just fine:

data Expr =
    Lam Text Expr
  | Var Text
  | App Expr Expr

I want to render expressions in my input and output languages in a Lisp-like manner. This is very easy to do using a good pretty-printing library; here I’m using the excellent wl-pprint-text, and will omit the ‘Pretty’ instances in the body of my post. But I’ll link to a gist including them at the bottom.

When performing a mechanical CPS transform, one targets both “atomic” expressions – i.e., variables and lambda abstractions – and “complex” expressions, i.e. function application. The target language is thus a combination of the ‘AExpr’ and ‘CExpr’ types:

data AExpr =
    AVar Text
  | ALam [Text] CExpr

data CExpr =
    CApp AExpr [AExpr]

All the mechanical CPS transformations use variants on two functions going by the cryptic names m and t. m is responsible for converting atomic expressions in the input languages (i.e., variables and lambda abstractions) into atomic expressions in the target language (an atomic CPS expression). t is the actual CPS transformation; it converts an expression in the input language into CPS, invoking a specified continuation (already in the target language) on the result.

Let’s look at the naive transform. Here are m and t, prefixed by ‘n’ to indicate that they are naive. First, m:

nm :: Expr -> IO AExpr
nm expr = case expr of
  Lam var cexpr0 -> do
    k      <- gensym
    cexpr1 <- nt cexpr0 (AVar k)
    return (ALam [var, k] cexpr1)

  Var var -> return (AVar var)

  App {} -> error "non-atomic expression"

(N.b. you almost never want to use ‘error’ in a production implementation of anything. It’s trivial to wrap e.g. ‘MaybeT’ around the appropriate functions to handle the bogus pattern match on ‘App’ totally, but I just want to keep the types super simple here.)

The only noteworthy thing that m does here is in the case of a lambda abstraction: a new abstract continuation is generated, and the body of the abstraction is converted to CPS via t, such that the freshly-generated continuation is called on the result. Remember, m is really just mapping atomic expressions in the input language to atomic expressions in the target language.

Here’s t for the naive transform. Remember, t is responsible for converting expressions to continuation-passing style:

nt :: Expr -> AExpr -> IO CExpr
nt expr cont = case expr of
  Lam {} -> do
    aexpr <- m expr
    return (CApp cont [aexpr])

  Var _ -> do
    aexpr <- m expr
    return (CApp cont [aexpr])

  App f e -> do
    fs <- gensym
    es <- gensym

    let aexpr0 = ALam [es] (CApp (AVar fs) [AVar es, cont])

    cexpr <- nt e aexpr0

    let aexpr1 = ALam [fs] cexpr

    nt f aexpr1

For both kinds of atomic expressions (lambda and variable), the expression is converted to the target language via m, and then the supplied continuation is applied to it. Very simple.

In the case of function application (a “complex”, or non-atomic expression), both the function to be applied, and the argument it is to be applied to, must be converted to CPS. This is done by generating two fresh continuations, transforming the argument, and then transforming the function. The control flow here is always handled by stitching continuations together; notice when transforming the function ‘f’ that the continuation to be applied has already handled its argument.

Next, the higher-order transform. Here are m and t:

hom :: Expr -> IO AExpr
hom expr = case expr of
  Lam var e -> do
    k  <- gensym
    ce <- hot e (\rv -> return (CApp (AVar k) [rv]))
    return (ALam [var, k] ce)

  Var n  -> return (AVar n)

  App {} -> error "non-atomic expression"

hot :: Expr -> (AExpr -> IO CExpr) -> IO CExpr
hot expr k = case expr of
  Lam {} -> do
    aexpr <- m expr
    k aexpr

  Var {} -> do
    aexpr <- m expr
    k aexpr

  App f e -> do
    rv      <- gensym
    xformed <- k (AVar rv)

    let cont = ALam [rv] xformed

        cexpr fs = hot e (\es -> return (CApp fs [es, cont]))

    hot f cexpr

Both of these have the same form as they do in the naive transform – the difference here is simply that the continuation to be applied to a transformed expression is expressed in the host language – i.e., here, Haskell. Thus the transform is “higher-order,” in exactly the same sense that higher-order abstract syntax is higher-order.

The final transformation I’ll illustrate here, the hybrid transform, applies the naive transformation to lambda and variable expressions, and applies the higher-order transformation to function applications. Here t is split up into tc and tk to handle these cases accordingly:

m :: Expr -> IO AExpr
m expr = case expr of
  Lam var cexpr -> do
    k       <- gensym
    xformed <- tc cexpr (AVar k)
    return (ALam [var, k] xformed)

  Var n  -> return (AVar n)

  App {} -> error "non-atomic expression"

tc :: Expr -> AExpr -> IO CExpr
tc expr c = case expr of
  Lam {} -> do
    aexpr <- m expr
    return (CApp c [aexpr])

  Var _ -> do
    aexpr <- m expr
    return (CApp c [aexpr])

  App f e -> do
    let cexpr fs = tk e (\es -> return (CApp fs [es, c]))
    tk f cexpr

tk :: Expr -> (AExpr -> IO CExpr) -> IO CExpr
tk expr k = case expr of
  Lam {} -> do
    aexpr <- m expr
    k aexpr

  Var {} -> do
    aexpr <- m expr
    k aexpr

  App f e -> do
    rv <- gensym
    xformed <- k (AVar rv)

    let cont     = ALam [rv] xformed
        cexpr fs = tk e (\es -> return (CApp fs [es, cont]))

    tk f cexpr

Matt illustrates these transformations on a simple expression: (g a). We can do the same:

test :: Expr
test = App (Var "g") (Var "a")

First, the naive transform. Note all the noisy administrative redexes that come along with it:

> cexpr <- nt test (AVar "halt")
> PP.pretty cexpr
((λ ($v1).
  ((λ ($v2).
    ($v1 $v2 halt)) a)) g)

The higher-order transform does better, containing only one such redex (an eta-expansion). Note that since the supplied continuation must be expressed in terms of a Haskell function, we need to write it in a more HOAS-y style:

> cexpr <- hot test (\ans -> return (CApp (AVar "halt") [ans]))
> PP.pretty cexpr
(g a (λ ($v3).
  (halt $v3)))

Finally the hybrid transform, which, here, is literally perfect. We don’t even need to deal with the minor annoyance of the HOAS-style continuation when calling it:

> cexpr <- tc test (AVar "halt")
> PP.pretty cexpr
(g a halt)

Matt goes on to describe a “partioned CPS transform” that can be used to recover a stack, in (seemingly) much the same manner that the defunctionalised manual CPS transform worked in the previous section. Very neat, but something I’ll have to look at in another post.

Fin

CPS is pretty gnarly. My experience in compiling with continuations is not substantial, but I dig learning it. Appel’s book, in particular, is meaty – expect more posts on the subject here eventually, probably.

‘Til next time! I’ve dumped the code from the latter section into a gist.

Embedded DSLs for Bayesian Modelling and Inference: a Retrospective

Why does my blog often feature its typical motley mix of probability, functional programming, and computer science anyway?

From 2011 through 2017 I slogged through a Ph.D. in statistics, working on it full time in 2012, and part-time in every other year. It was an interesting experience. Although everything worked out for me in the end – I managed to do a lot of good and interesting work in industry while still picking up a Ph.D. on the side – it’s not something I’d necessarily recommend to others. The smart strategy is surely to choose one thing and give it one’s maximum effort; by splitting my time between work and academia, both obviously suffered to some degree.

That said, at the end of the day I was pretty happy with the results on both fronts. On the academic side of things, the main product was a dissertation, Embedded Domain-Specific Languages for Bayesian Modelling and Inference, supporting my thesis: that novel and useful DSLs for solving problems in Bayesian statistics can be embedded in statically-typed, purely functional programming languages.

It helps to remember that in this day and age, one can still typically graduate by, uh, “merely” submitting and defending a dissertation. Publishing in academic venues certainly helps focus one’s work, and is obviously necessary for a career in academia (or, increasingly, industrial research). But it’s optional when it comes to getting your degree, so if it doesn’t help you achieve your goals, you may want to reconsider it, as I did.

The problem with the dissertation-first approach, of course, is that nobody reads your work. To some extent I think I’ve mitigated that; most of the content in my dissertation is merely a fleshed-out version of various ideas I’ve written about on this blog. Here I’ll continue that tradition and write a brief, informal summary of my dissertation and Ph.D. more broadly – what I did, how I approached it, and what my thoughts are on everything after the fact.

The Idea

Following the advice of Olin Shivers (by way of Matt Might), I oriented my work around a concrete thesis, which wound up more or less being that embedding DSLs in a Haskell-like language can be a useful technique for solving statistical problems. This thesis wasn’t born into the world fully-formed, of course – it began as quite a vague (or misguided) thing, but matured naturally over time. Using the tools of programming languages and compilers to do statistics and machine learning is the motivation behind probabilistic programming in general; what I was interested in was exploring the problem in the setting of languages embedded in a purely functional host. Haskell was the obvious choice of host for all of my implementations.

It may sound obvious that putting together a thesis is a good strategy for a Ph.D. But here I’m talking about a thesis in the original (Greek) sense of a proposition, i.e. a falsifiable idea or claim (in contrast to a dissertation, from the Latin disserere, i.e. to examine or to discuss). Having a central idea to orient your work around can be immensely useful in terms of focus. When you read a dissertation with a clear thesis, it’s easy to know what the writer is generally on about – without one it can (increasingly) be tricky.

My thesis is pretty easy to defend in the abstract. A DSL really exposes the structure of one’s problem while also constraining it appropriately, and embedding one in a host language means that one doesn’t have to implement an entire compiler toolchain to support it. I reckoned that simply pointing the artillery of “language engineering” at the statistical domain would lead to some interesting insight on structure, and maybe even produce some useful tools. And it did!

The Contributions

Of course, one needs to do a little more defending than that to satisfy his or her examination committee. Doctoral research is supposed to be substantial and novel. In my experience, reviewers are concerned with your answers to the following questions:

• What, specifically, are your claims?
• Are they novel contributions to your field?
• Have you backed them up sufficiently?

At the end of the day, I claimed the following advances from my work.

• Novel probabilistic interpretations of the Giry monad’s algebraic structure. The Giry monad (Lawvere, 1962; Giry, 1981) is the “canonical” probability monad, in a meaningful sense, and I demonstrated that one can characterise the measure-theoretic notion of image measure by its functorial structure, as well as the notion of product measure by its monoidal structure. Having the former around makes it easy to transform the support of a probability distribution while leaving its density structure invariant, and the latter lets one encode probabilistic independence, enabling things like measure convolution and the like. What’s more, the analogous semantics carry over to other probability monads – for example the well-known sampling monad, or more abstract variants.

• A novel characterisation of the Giry monad as a restricted continuation monad. Ramsey & Pfeffer (2002) discussed an “expectation monad,” and I had independently come up with my own “measure monad” based on continuations. But I showed both reduce to a restricted form of the continuation monad of Wadler (1994) – and that indeed, when the return type of Wadler’s continuation monad is restricted to the reals, it is the Giry monad.

  To be precise it’s actually somewhat more general – it permits integration with respect to any measure, not only a probability measure – but that definition strictly subsumes the Giry monad. I also showed that product measure, via the applicative instance, yields measure convolution and associated operations.

• A novel technique for embedding a statically-typed probabilistic programming language in a purely functional language. The general idea itself is well-known to those who have worked with DSLs in Haskell: one constructs a base functor and wraps it in the free monad. But the reason that technique is appropriate in the probabilistic programming domain is that probabilistic models are fundamentally monadic constructs – merely recall the existence of the Giry monad for proof!

  To construct the requisite base functor, one maps some core set of concrete probability distributions denoted by the Giry monad to a collection of abstract probability distributions represented only by unique names. These constitute the branches of one’s base functor, which is then wrapped in the familiar ‘Free’ machinery that gives one access to the functorial, applicative, and monadic structure that I talked about above. This abstract representation of a probabilistic model allows one to implement other probability monads, such as the well-known sampling monad (Ramsey & Pfeffer, 2002; Park et al., 2008) or the Giry monad, by way of interpreters.

  (N.b. Ścibior et al. (2015) did some very similar work to this, although the monad they used was arguably more operational in its flavour.)

• A novel characterisation of execution traces as cofree comonads. The idea of an “execution trace” is that one runs a probabilistic program (typically generating a sample) and then records how it executed – what randomness was used, the execution path of the program, etc. To do Bayesian inference, one then runs a Markov chain over the space of possible execution traces, calculating statistics about the resulting distribution in trace space (Wingate et al., 2011).

  Remarkably, a cofree comonad over the same abstract probabilistic base functor described above allows us to represent an execution trace at the embedded language level itself. In practical terms, that means one can denote a probabilistic model, and then run a Markov chain over the space of possible ways it could have executed, without leaving GHCi. You can alternatively examine and perturb the way the program executes, stepping through it piece by piece, as I believe was originally a feature in Venture (Mansinghka et al., 2014).

  (N.b. this really blew my mind when I first started toying with it, converting programs into execution traces and then manipulating them as first-class values, defining other probabilistic programs over spaces of execution traces, etc. Meta.)

• A novel technique for statically encoding conditional independence of terms in this kind of embedded probabilistic programming language. If you recall that I previously demonstrated the monoidal (i.e. applicative) structure of the Giry monad encodes the notion of product measure, it will not be too surprising to hear that I used the free applicative functor (Capriotti & Kaposi, 2014) (again, over the same kind of abstract probabilistic base functor) to reify applicative expressions such that they can be identified statically.

• A novel shallowly-embedded language for building custom transition operators for use in Markov chain Monte Carlo. MCMC is the de-facto standard way to perform inference on Bayesian models (although it is not limited to Bayesian models in particular). By wrapping a simple state monad transformer around a probability monad, one can denote Markov transition operators, combine them, and transform them in a few ways that are useful for doing MCMC.

  The framework here was inspired by the old parallel “strategies” idea of Trinder et al. (1998). The idea is that you want to “evaluate” a posterior via MCMC, and want to choose a strategy by which to do so – e.g. Metropolis (Metropolis, 1953), slice sampling (Neal, 2003), Hamiltonian (Neal, 2011), etc. Since Markov transition operators are closed under composition and convex combinations, it is easy to write a little shallowly-embedded combinator language for working with them – effectively building evaluation strategies in a manner familiar to those who’ve worked with Haskell’s parallel library.

  (N.b. although this was the most trivial part of my research, theoretical or implementation-wise, it remains the most useful for day-to-day practical work.)

The Execution

One needs to stitch his or her contributions together in some kind of over-arching narrative that supports the underlying thesis. Mine went something like this:

The Giry monad is appropriate for denoting probabilistic semantics in languages with purely-functional hosts. Its functorial, applicative, and monadic structure denote probability distributions, independence, and marginalisation, respectively, and these are necessary and sufficient for encoding probabilistic models. An embedded language based on the Giry monad is type-safe and composable.

Probabilistic models in an embedded language, semantically denoted in terms of the Giry monad, can be made abstract and interpretation-independent by defining them in terms of a probabilistic base functor and a free monad instead. They can be forward-interpreted using standard free monad recursion schemes in order to compute probabilities (via a measure intepretation) or samples (via a sampling interpretation); the latter interpretation is useful for performing limited forms of Bayesian inference, in particular. These free-encoded models can also be transformed into cofree-encoded models, under which they represent execution traces that can be perturbed arbitrarily by standard comonadic machinery. This representation is amenable to more elaborate forms of Bayesian inference. To accurately denote conditional independence in the embedded language, the free applicative functor can also be used.

One can easily construct a shallowly-embedded language for building custom Markov transitions. Markov chains that use these compound transitions can outperform those that use only “primitive” transitions in certain settings. The shallowly embedded language guarantees that transitions can only be composed in well-defined, type-safe ways that preserve the properties desirable for MCMC. What’s more, one can implement “transition transformers” for implementing still more complex inference techniques, e.g. annealing or tempering, over existing transitions.

Thus: novel and useful domain-specific languages for solving problems in Bayesian statistics can be embedded in statically-typed, purely-functional programming languages.

I used the twenty-minute talk period of my defence to go through this narrative and point out my claims, after which I was grilled on them for an hour or two. The defence was probably the funnest part of my whole Ph.D.

The Product

In the end, I mainly produced a dissertation, a few blog posts, and some code. By my count, the following repos came out of the work:

• deanie: An embedded probabilistic programming language.
  http://github.com/jtobin/deanie

• declarative: DIY Markov Chains.
  http://github.com/jtobin/declarative

• flat-mcmc: Painless general-purpose sampling.
  http://github.com/jtobin/flat-mcmc

• hasty-hamiltonian: Speedy traversal through parameter space.
  http://github.com/jtobin/hasty-hamiltonian

• hnuts: Automatic gradient-based sampling.
  http://github.com/jtobin/hnuts

• lazy-langevin: Gradient-based diffusion.
  http://github.com/jtobin/lazy-langevin

• mcmc-types: Common types for implementing MCMC algorithms.
  https://github.com/jtobin/mcmc-types

• measurable: A shallowly-embedded DSL for basic measure wrangling.
  http://github.com/jtobin/measurable

• mighty-metropolis: The Metropolis sampling algorithm.
  http://github.com/jtobin/mighty-metropolis

• mwc-probability: Sampling function-based probability distributions.
  http://github.com/jtobin/mwc-probability

• sampling: Tools for sampling from collections.
  https://github.com/jtobin/sampling

• speedy-slice: Speedy slice sampling.
  http://github.com/jtobin/speedy-slice

If any of this stuff is or was useful to you, that’s great! I still use the declarative libraries, flat-mcmc, mwc-probability, and sampling pretty regularly. They’re fast and convenient for practical work.

Some of the other stuff, e.g. measurable, is useful for building intuition, but not so much in practice, and deanie, for example, is a work-in-progress that will probably not see much more progress (from me, at least). Continuing from where I left off might be a good idea for someone who wants to explore problems in this kind of setting in the future.

General Thoughts

When I first read about probabilistic (functional) programming in Dan Roy’s 2011 dissertation I was absolutely blown away by the idea. It seemed that, since there was such an obvious connection between the structure of Bayesian models and programming languages (via the underlying semantic graph structure, something that has been exploited to some degree as far back as BUGS), it was only a matter of time until someone was able to really create a tool that would revolutionize the practice of Bayesian statistics.

Now I’m much more skeptical. It’s true that probabilistic programming tends to expose some beautiful structure in statistical models, and that a probabilistic programming language that was easy to use and “just worked” for inference would be a very useful tool. But putting something expressive and usable together that also “just works” for that inference step is very, very difficult. Very difficult indeed.

Almost every probabilistic programming framework of the past ten years, from Church down to my own stuff, has more or less wound up as “thesisware,” or remains the exclusive publication-generating mechanism of a single research group. The exceptions are almost unexceptional in of themselves: JAGS and Stan are probably the most-used such frameworks, certainly in statistics (I will mention the very honourable PyMC here as well), but they innovate little, if at all, over the original BUGS in terms of expressiveness. Similarly it’s very questionable whether the fancy MCMC algo du jour is really any better than some combination of Metropolis-Hastings (even plain Metropolis), Gibbs (or its approximate variant, slice sampling), or nested sampling in anything outside of favourably-engineered examples (I will note that Hamiltonian Monte Carlo could probably be counted in there too, but it can still be quite a pain to use, its variants are probably overrated, and it is comparatively expensive).

Don’t get me wrong. I am a militant Bayesian. Bayesian statistics, i.e., as far as I’m concerned, probability theory, describes the world accurately. And there’s nothing wrong with thesisware, either. Research is research, and this is a very thorny problem area. I hope to see more abandoned, innovative software that moves the ball up the field, or kicks it into another stadium entirely. Not less. The more ingenious implementations and sampling schemes out there, the better.

But more broadly, I often find myself in the camp of Leo Breiman, who in 2001 characterised the two predominant cultures in statistics as those of data modelling and algorithmic modelling respectively, the latter now known as machine learning, of course. The crux of the data modelling argument, which is of course predominant in probabilistic programming research and Bayesian statistics more generally, is that a practitioner, by means of his or her ingenuity, is able to suss out the essence of a problem and distill it into a useful equation or program. Certainly there is something to this: science is a matter of creating hypotheses, testing them against the world, and iterating on that, and the “data modelling” procedure is absolutely scientific in principle. Moreover, with a hat tip to Max Dama, one often wants to impose a lot of structure on a problem, especially if the problem is in a domain where there is a tremendous amount of noise. There are many areas where this approach is just the thing one is looking for.

That said, it seems to me that a lot of the data modelling-flavoured side of probabilistic programming, Bayesian nonparametrics, etc., is to some degree geared more towards being, uh, “research paper friendly” than anything else. These are extremely seductive areas for curious folks who like to play at the intersection of math, statistics, and computer science (raises hand), and one can easily spend a lifetime chasing this or that exquisite theoretical construct into any number of rabbit holes. But at the end of the day, the data modelling culture, per Breiman:

.. has at its heart the belief that a statistician, by imagination and by looking at the data, can invent a reasonably good parametric class of models for a complex mechanism devised by nature.

Certainly the traditional statistics that Breiman wrote about in 2001 was very different from probabilistic programming and similar fields in 2018. But I think there is the same element of hubris in them, and to some extent, a similar dissociation from reality. I have cultivated some of the applied bent of a Breiman, or a Dama, or a Locklin, so perhaps this should not be too surprising.

I feel that the 2012-ish resurgence of neural networks jolted the machine learning community out of a large-scale descent into some rather dubious Bayesian nonparametrics research, which, much as I enjoy that subject area, seemed more geared towards generating fun machine learning summer school lectures and NIPS papers than actually getting much practical work done. I can’t help but feel that probabilistic programming may share a few of those same characteristics. When all is said and done, answering the question is this stuff useful? often feels like a stretch.

So: onward & upward and all that, but my enthusiasm has been tempered somewhat, is all.

Fini

Administrative headaches and the existential questions associated with grad school aside, I had a great time working in this area for a few years, if in my own aloof and eccentric way.

If you ever interacted with this area of my work, I hope you got some utility out of it: research ideas, use of my code, or just some blog post that you thought was interesting during a slow day at the office. If you’re working in the area, or are considering it, I wish you success, whether your goal is to build practical tools, or to publish sexy papers. :-)