Random orthogonal matrices play an important role in probability and statistics, arising in multivariate analysis, directional statistics, and models of physical systems, among other areas. Calculations involving random orthogonal matrices are complicated by their constrained support. Accordingly, we parametrize the Stiefel and Grassmann manifolds, represented as subsets of orthogonal matrices, in terms of Euclidean parameters using the Cayley transform. We derive the necessary Jacobian terms for change of variables formulas. Given a density defined on the Stiefel or Grassmann manifold, these allow us to specify the corresponding density for the Euclidean parameters, and vice versa. As an application, we describe and illustrate through examples a Markov chain Monte Carlo approach to simulating from distributions on the Stiefel and Grassmann manifolds. Finally, we establish an asymptotic independent normal approximation for the distribution of the Euclidean parameters which corresponds to the uniform distribution on the Stiefel manifold. This result contributes to the growing literature on normal approximations to the entries of random orthogonal matrices or transformations thereof.Closely related to the Stiefel manifold is the Grassmann manifold G(k, p), the set of k-dimensional linear subspaces of R p . The Grassmann manifold has dimension d G = (p − k)k. Typically, points in the Grassmann manifold are thought of as equivalence classes of V(k, p), where two orthogonal matrices belong to the same class if they share the same column space or, equivalently, if one matrix can be obtained from the other through right multiplication by an element of O(k). In Section 3.2, we elaborate and expand upon the contributions of Shepard et al. [42] to provide another representation of the Grassmann manifold G(k, p) as the subset V + (k, p) of p × k orthogonal matrices having a symmetric positive definite (SPD) top block. In this article, we focus on orthogonal matrices having fewer columns than rows.Both the Stiefel and Grassmann manifolds can be equipped with a uniform probability measure, also know as an invariant or Haar measure. The uniform distribution P V(k,p) on V(k, p) is characterized by its invariance to left and right multiplication by orthogonal matrices: If Q ∼ P V(k,p) , then U QV dist.

No abstract

Motivated by applications to Bayesian inference for statistical models with orthogonal matrix parameters, we present 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 Q, we construct a distribution for an unconstrained random matrix X such that Q X , the orthogonal component of the polar decomposition of X, is equal in distribution to Q. The distribution of X is amenable to Markov chain Monte Carlo (MCMC) simulation using standard methods, and an approximation to the distribution of 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.

No abstract