add Givens-based orthonormal matrix transform and inverse plus Jacobian

Description

We would like to have an orthonormal matrix type for constrained variables in Stan. That requires defining the transform and its inverse plus log Jacobian determinant in the math library first.

Among the applications are factor models and parameterizations of positive definite matrices in terms of eigenvectors.

This would amount to a data type orthonormal[P, Q] in Stan and three functions in the math library, which have roughly this signature (see, e.g., corr_matrix_free.hpp and corr_matrix_constrain.hpp in directory stan/math/prim/fun. This amounts to two files with a total of three function signatures.

orthonormal_free.hpp

Eigen::Matrix<T, 1, -1> orthonormal_free(const Matrix<T, -1, -1>& y);

orthonormal_constrain.hpp

// w/o Jacobian adjustment
Eigen::Matrix<T, -1, -1> orthonormal_constrain(const Matrix<T, -1, 1>& x, int P, int Q);

// w Jacobian adjustment
Eigen::Matrix<T, -1, -1> orthonormal_constrain(const Matrix<T, -1, 1>& x, int P, int Q, T& lp);

Example

Reference

We can follow the code and transforms laid out here:

Arya A. Pourzanjani, Richard M. Jiang, Brian Mitchell, Paul J. Atzberger, Linda R. Petzold. 2021. Bayesian Inference over the Stiefel Manifold via the Givens Representation. Bayesian Analysis.

Derivative reduction

Algorithm 1 has the constraining transform, where angles are defined by

theta[i, j] = atan2(y[i, j], y[i, i]).

Equation (4.5) lists the Jacobian adjustment, which involves terms that on the log scale are equal to

(j - i - 1) * log(cos(theta[i, j]))

We can exploit the fact that

d/du log cos(u) = -tan(u)

and the definition of theta[i, j] to reduce to

d/d.theta[i, j] (j - i - 1) * log(cos(theta[i, j])) =  (j - i - 1) * y[i, j] / y[i, i].

Auxiliary variable approach

What I don't see how to do immediately is deal with the overparameterization and unit normal "prior" that makes Pourzanjani et al.'s approach feel very much like Marsaglia's approach to unit vectors that we currently employ.

Unconstraining transform

We also need the transform to go from the constrained to the unconstrained spaces (for initialization). I couldn't find it in the paper.

Original GitHub repo

More details can probably be extracted from Pourzanjani's GitHub repo, TfRotationPca.

Relation to unit vector coding

We should be able to code unit vectors as (p, 1) orthonormal matrices. The two approaches with auxiliary variables are closely related, and may be identical.

Other repos

We also need to add the doc for the transform to the reference manual (docs repo) and the actual data type to the language and code generate (stanc3 repo).

Unit and negative unit determinants

@bgoodri and Pourzanjani were discussing the fact that the Givens approach leads to solutions which can have positive or negative unit determinants. Maybe we could constrain to one in the prior, but in higher dimensions, this turns into combinatorial multimodality.

Alternatives

There is more than one way to do this. Pourzanjani's paper considers the geodesic Monte Carlo (GMC) approach of Byrne and Girolami. It's also possible to use Cayley transforms.

Expected Output

Invertible transforms with Jacobians so we can sample uniformly overly over Stiefel manifolds.

Current Version:

v4.1.0

Start with stan/math/prim/fun/corr_matrix_free.hpp and corr_matrix_constrain.hpp to understand the existing transform conventions and signatures. Then read the cited Pourzanjani et al. paper and its TfRotationPca repository, focusing on the Givens transform, inverse, overparameterization, and Jacobian. Done means the math library has invertible orthonormal transforms with Jacobian adjustments suitable for sampling Stiefel manifolds.