Introduction

Main principles of development

The Julia Klara package provides an interoperable generic engine for a breadth of Markov Chain Monte Carlo (MCMC) methods.

The idea that there exists no unique optimal MCMC methodology for all purposes has been a development cornerstone. Along these lines, interest is in providing a wealth of Monte Carlo strategies and options, letting the user decide which algorithm suits their use case. Such “agnostic” approach to coding permeates Klara from top to bottom, offering a variety of methods and detailed configuration, ultimately leading to rich functionality. Klara’s wide range of functionality makes itself useful in applications and as a test bed for comparative methodological research. It also gives the flexibility to connect Klara with various other packages, exploiting different levels of ongoing developments in them.

In fact, interoperability has been another central principle of development. A high-level API enables the user to implement their application effectively, while it facilitates connectivity to Julia packages. Whenever deemed necessary, minor wrappers are provided to integrate Klara with other packages such as ReverseDiffSource and ForwardDiff.

The high-level API sits atop of a low-level one. The latter aims at providing an extensible codebase for developers interested in adding new functionality and offers an alternative interface for users who prefer more hands-on access to the underlying routines.

Speed of execution has been another motivation behind the low-level API. Passing from the higher to the lower level interface allows to exploit Julia’s meta-programming capabilities by generating code dynamically and by substituting dictionaries by vectors internally. Memory footprint and garbage collection have been kept to a minimum without compromising ease of use thanks to the duality of higher and lower level APIs.

Features

A summary of Klara’s main features follows:

• Graph-based model specification. Representing the model as a graph widens the scope of accommodated models and enables exploiting graph algorithms from the Graphs package.
• Diverse options for defining model parameters. Parameters can be defined on the basis of a log-target or they can be introduced in a Bayesian fashion via their log-likelihood and log-prior. Parameter targets, likelihoods and priors can be specified via functions or distributions. Klara’s integration with the Distributions package facilitates parameter definition via distributions.
• Job-centric simulations. The concept of MCMC simulation has been separated from that of model specification. Job types indicate the context in which a model is simulated. For example, a BasicMCJob instance determines how to sample an MCMC chain for a model with a single parameter, whereas a GibbsJob provide Gibbs sampling for more complex models involving several parameters.
• Customized job flow. Job control flow comes in two flavors, as it can be set to ordinary loop-based flow or it can be managed by Julia’s tasks (coroutines). Job management with tasks allows MCMC simulations to be suspended and resumed in a flexible manner.
• Wide range of Monte Carlo samplers. A range of MCMC samplers is available, including accept-reject and slice sampling, Metropolis-Hastings algorithm, No-U-Turn (NUTS) sampling, and geometric MCMC schemes, such as Riemann manifold Langevin and Hamiltonian Monte Carlo. Adaptive samplers and empirical tuning are included in Klara as a means to faster convergence. It is noted that most of these samplers need to be ported from the older version of Klara, which is work in progress.
• MCMC summary statistics and convergence diagnostics. Main routines for computing the effective sampling size and integrated autocorrelation time have been coded, while there is a roadmap to provide more convergence diagnostics tools (note to user; this functionality will also be ported soon from the older version of Klara).
• States and chains. Proposed Monte Carlo samples are organized systematically with the help of a state and chain type system. This way, values can be passed around and stored without re-allocating memory. At the same time, the state/chain type system offers scope for extending the current functionality if it is required to store less usual components.
• Detailed configuration of output storage in memory or in file. The chain resulting from a Monte Carlo simulation can be saved in memory or can be written directly to a file stream. Detailed output configuration is possible, allowing to select which elements to save and which to omit from the final output.
• Automatic differentiation for MCMC sampling. Some Monte Carlo methods require the gradient or higher order derivatives of the log-target. If these derivatives are not user-inputted explicitly, Klara can optionally compute them using reverse or forward mode automatic differentiation. For this purpose, Klara uses ReverseDiffSource and ForwardDiff under the hood.

Preliminary exposition of graph models

Klara’s graph model is presented concisely in the current section as a smooth introduction to subsequent elaborate chapters. A single model type, named GenericModel, serves as the sole entry point for defining any graph model in Klara. GenericModel has been inspired by and operates on par with GenericGraph, a versatile graph type of the Graphs package.

GenericModel can be conceptualized as a graph whose nodes represent the underlying model’s variables and its edges specify the dependencies between these variables. As it becomes obvious, the main front end of GenericModel consists of its vertices and edges fields, which are of type Vector{Variable} and Vector{Dependence} respectively. Without going into details, each vertex is defined as constant, data, transformation or parameter, all being Variable subtypes. A single non-abstract Dependence type suffices to describe variable dependencies.

Typical probabilistic graphical models, such as directed acyclic graphs (DAGs) and factor graphs, are permitted in GenericModel. Similarly to GenericGraph, the is_directed field of GenericModel dictates whether the model is directed or not. Practically, Gibbs sampling is not affected by the distinction between directed and non-directed graphs. However, effort has been made to define GenericModel generically in order to provide scope for future developments if the need arises to distinguish between DAGs and factor graphs in programming practice.

The graph-oriented definition of GenericModel finds its main utility in Klara‘s optional declarative model specification. In other words, it is possible to delegate responsibility of variable ordering to GenericModel via topological sorting of the graph. Furthermore, the statistical model is easier to disseminate by visualizing GenericModel as a graph.

Topological sorting and graph visualization are achieved via “outsourcing”. In particular, converting GenericModel to its corresponding GenericGraph allows to harness sorting routines in the Graphs package. Moreover, GenericModel is convertible to DOT format, thus making it possible to use the DOT graph description language for model visualization.