More efficient sampling on the Stiefel manifold

[sethaxen]

The common way to draw points from the Hausdorff measure on the Stiefel manifold is using the SVD decomposition. However, one can also draw points using the unique QR decomposition, which is much faster. Surprisingly, Chikuse's book devotes only a few lines to this approach, and not by name.

What's unfortunate is that the method computes n × k random numbers but in the end returns a matrix with (in the real case) only nk - k(k+1)/2 degrees of freedom. That is, k(k+1)/2 of those random numbers are wasted. Moreover, the method requires computing the QR decomposition for each random draw, which is O(n^2 * k).

Stewart noted that in the unique QR decomposition, the Q and R matrices are independent and used this to work out a technique for expressing the Q matrix as a product of Householder matrices, which can be compactly represented using a strict triangle of an n × k matrix:

G. W. Stewart
The Efficient Generation of Random Orthogonal Matrices with an Application to Condition Estimators. SIAM J. Numer. Anal., 17(3), 403–409.
https://epubs.siam.org/doi/abs/10.1137/0717034

The method only generates k excess random numbers and can be expressed in 3 stages:

1. fill a triangular n × k matrix T with IID standard normal elements.
2. bijectively map T ↦ (Q, D), where Q is drawn from the Hausdorff measure on St(n, k), and D is a diagonal matrix with positive entries
3. discard D

The Q matrix should if possible be represented using one of the special compact types in LinearAlgebra, probably LinearAlgebra.QRPackedQ, so that the Q matrix is only explicitly built if necessary.
Every orthogonal matrix can be expressed as the product of a sequence of Householder transformations and a diagonal matrix, so in general we cannot use the QR compact representations to store Q. We can either devise our own, or, more cleanly, explicitly multiply Q * Diagonal(sign.(R_diag)) to generate the dense matrix.

This gives us the framework for e.g. sampling semi-orthogonal matrices in a PPL with HMC. If Bijectors adds supports for non-bijective maps (which I understand is planned), then we can add the bijective map in (2) and the non-bijective part in (3).