Instead, we can turn to automatic differentiation (AD, also sometimes known as algorithmic differentiation). A detailed explanation of AD is beyond the scope of this paper, so we refer the reader to Griewank et al. (2008), or Section 8.2 in Nocedal and Wright (2006). Put simply, AD treats a function as a composite of subfunctions, and computes derivatives by repeatedly applying the chain rule. In practical terms, AD involves coding the log posterior using a specialized numerical library that keeps track of the derivatives of these subfunctions. When we compile the function that computes the log posterior, the AD library will “automatically” generate additional functions that return the gradient, the Hessian, and even higher-order derivatives. 4 The remarkable feature of AD is that computing the gradient of a scalar-valued function takes no more than five times as long as computing the log posterior, regardless of the number of parameters (Griewank et al. 2008, p. xii). If the cost of the log posterior grows linearly in N , so will the cost of the gradient. In most statistical software packages, like R , the default storage mode for any matrix is in a “dense” format; each element in the matrix is stored explicitly, regardless of the value. For a model with n variables, this matrix consists of n 2 numbers, each consuming eight bytes of mem- ory at double precision. If we have a dataset in which N = 10000, k = 5 and p is relatively small, the Hessian for this model with 50, 000 + p variables will consume more than 20GB of RAM. Fur- thermore, the computational effort for matrix-vector multiplication is quadratic in the number of columns, and matrix-matrix multiplication is cubic. To the extent that either of these operations is used in the mode-finding or sampling steps, the computational effort will grow much faster than the size of the dataset. Since multiplying a triangular matrix is roughly one-sixth as expensive as multiplying a full matrix, we could gain some efficiency by working with the Cholesky decompo- sition of the Hessian instead. However, the complexity of the Cholesky decomposition algorithm itself will still be cubic in N (Golub and Van Loan 1996, Ch. 1). For our purposes, the source of scalability is in the sparsity of the Hessian. If the vast majority of elements in a matrix are zero, we do not need to store them explicitly. Instead, we need to store only the non-zero values, the row numbers of those values, and the index of the values that be- $^{4}$There are a number of established AD tools available for researchers to use for many different programming envi- ronments. For C++ , we use CppAD (Bell 2013), although ADOL-C (Walther and Griewank 2012) is also popular. We call our C++ functions from R (R Development Core Team 2014) using the Rcpp package (Eddelbuettel and François 2011; Bates and Eddelbuettel 2013). CppAD is also available for Python . Matlab users have access to ADMAT (Coleman and Verma 2000), among other options. 21