This post is an (unpaid) advertisement for the following extremely useful resource:
- Petersen, K. B. and M. S. Pedersen. 2008. The Matrix Cookbook. Tehcnical Report, Technical University of Denmark.
It contains 70+ pages of useful relations and derivations involving matrices. What grabbed my eye was the computation of gradients for matrix operations ranging from eigenvalues and determinants to multivariate normal density functions. I had no idea the multivariate normal had such a clean gradient (see section 8).
We’ve been playing around with Hamiltonian (aka Hybrid) Monte Carlo for sampling from the posterior of hierarchical generalized linear models with lots of interactions. HMC speeds up Metropolis sampling by using the gradient of the log probability to drive samples in the direction of higher probability density, which is particularly useful for correlated parameters that mix slowly with standard Gibbs sampling. Matt “III” Hoffman‘s already got it working in Cython/numpy for the discrete-predictor case.
To really get this going, we need to be able to handle the varying intercept/varying slope models used in the Red State/Blue State analysis of Gelman, Shor, Bafumi and Park, which is also explained by Andrew and Jennifer in their regression book. The slopes and intercepts get multivariate normal priors, the covariance matrices of which require hyperpriors. That means log probability functions involving operations like matrix-inverse products or eigenvalues (depending on how you model the covariance prior).
I’ve been following up some earlier suggestions on this blog to automatic differentiation. I’ve pretty much settled on using David Gay‘s elegant little Reverse Automatic Differentiation (RAD) package, which is a very straightforward C++ template-based library for computing gradients. All you need to do is replace the floating point calcs in code with templated versions.
This in turn means we need a templated C++ vector/matrix library with BLAS/LAPACK-like functionality. An exchange on Justin Domke’s blog led me to the Eigen package, which is a beautifully templated implementation of BLAS and at least what I need from LAPACK. Their doc led me to the cookbook.