  Motivated by applications to Bayesian inference for statistical models with
orthogonal matrix parameters, we present $\textit{polar expansion},$ a general
approach to Monte Carlo simulation from probability distributions on the
Stiefel manifold. To bypass many of the well-established challenges of
simulating from the distribution of a random orthogonal matrix
$\boldsymbol{Q},$ we construct a distribution for an unconstrained random
matrix $\boldsymbol{X}$ such that $\boldsymbol{Q}_X,$ the orthogonal component
of the polar decomposition of $\boldsymbol{X},$ is equal in distribution to
$\boldsymbol{Q}.$ The distribution of $\boldsymbol{X}$ is amenable to Markov
chain Monte Carlo (MCMC) simulation using standard methods, and an
approximation to the distribution of $\boldsymbol{Q}$ can be recovered from a
Markov chain on the unconstrained space. When combined with modern MCMC
software, polar expansion allows for routine and flexible posterior inference
in models with orthogonal matrix parameters. We find that polar expansion with
adaptive Hamiltonian Monte Carlo is an order of magnitude more efficient than
competing MCMC approaches in a benchmark protein interaction network
application. We also propose a new approach to Bayesian functional principal
components analysis which we illustrate in a meteorological time series
application.
