Add orthogonal and semi-orthogonal matrix types

[spinkney]

Add orthogonal and semi-orthogonal matrix types. These are described in https://en.wikipedia.org/wiki/Orthogonal_matrix. There are numerous applications and many of them are detailed in Michael Jauch's PhD thesis https://dukespace.lib.duke.edu/items/da7e8aed-efe9-4006-bbd9-ecc3136e0d3b.

An efficient parameterization is given in https://github.com/sethaxen/stan_semiorthogonal_transforms and the main Stan function file https://github.com/sethaxen/stan_semiorthogonal_transforms/blob/main/transforms/semiorthogonal_reflector_transform.stan. These are compared to other paramerterizations in https://discourse.mc-stan.org/t/parameterizing-an-orthonormal-matrix-the-stiefel-manifold/33211/11?u=spinkney.

The proposal is to have two types

parameters {
 orthogonal_matrix[K] x;
 semi_orthogonal_matrix[K, N] y;
}

The orthogonal matrix is a square K x K matrix with orthonormal

The number of parameters is in the above links as

/**
 * Return the number of unconstrained parameters needed for a N x K semiorthogonal matrix Q.
 * 
 * @param N number of rows of Q
 * @param K number of columns of Q (must be <= N)
 * @param special whether Q is special orthogonal (i.e. det(Q)=+1). Must be in (0,1). Only
 *                checked if N==K.
 * @return nparams Number of unconstrained parameters
 */
int semiorthogonal_reflector_num_params(int N, int K, int special){
  return N * K - (K * (K - 1)) %/% 2 - (N==K && special);
}

This contrasts with the "naive" method of generating semi-orthogonal matrices of K x N. The Householder reflector parameterization has choose(K, 2) fewer parameters. In the case of generating a square orthogonal matrix there is an additional reduction in the number of parameters by 1 due to restricting diagonal sign to be positive. We don't need the "special" keyword as we would fix the sign indeterminancy to positive 1, just like we currently do with the QR decomposition in Stan.

As noted in the forum post by Seth that the implementation would be adding a type and then having efficient ways to do multiplication, determinant, and linear solves with these factors:

Finally, reflectors_factors imagines what if Stan added a native semiorthogonal_matrix[K,N] data type that stored the reflector factors and for which certain operations like multiplication, determinant, and linear solve were implemented in terms of those factors. Efficient algorithms exist for each of these and are implemented in Eigen and LAPACK. This is really the only option with these Householder-inspired transforms if one needs to efficiently infer a large orthogonal matrix.

Adding various questions and applications about orthogonal matrices from the forums

• https://discourse.mc-stan.org/t/how-to-constrain-a-matrix-comprising-several-parameters-to-be-an-orthogonal-matrix/11106?u=spinkney
• https://discourse.mc-stan.org/t/efficient-orthogonal-matrix-parameterization/8685/3
• https://discourse.mc-stan.org/t/orthogonal-not-orthonormal-matrix-parameters/21946/20
• https://discourse.mc-stan.org/t/generating-orthogonal-latents/30920/18
• https://discourse.mc-stan.org/t/new-transform-for-orthonormal-matrices-in-stan/2981/3

[WardBrian]

As noted in the forum post by Seth that the implementation would be adding a type and then having efficient ways to do multiplication, determinant, and linear solves with these factors:

This would be an interesting implementation challenge. Currently, all of our constrained types are type-equivalent to their unconstrained types at runtime (e.g. a ordered vector is just an Eigen::Vector in the end). That choice makes it so things like multiplication don't have to be specialized for every single constrained type, but it also makes it very difficult to do specialization.

As @syclik pointed out to me a few years ago, this also means that we're doing redundant checks at the moment in a lot of cases where you declare a parameter with a constraint and then pass it to a function which has to verify that constraint. In theory we know enough to say that check should be turned off, but in practice we've never had the machinery to do so.

Do you think that this kind of ability to treat the parameter distinctly at runtime is critical to the usefulness of this feature?

[spinkney]

Do you think that this kind of ability to treat the parameter distinctly at runtime is critical to the usefulness of this feature?

No, it is not critical. In fact we can just efficiently construct the matrices using LAPACK or Eigen functions.

We could include additional functions like multiply_left_ortho(Q, X) and the right one as well since Q^T = Q^-1 we don't need to have a solve type function, just the multiplication.

[WardBrian]

@bob-carpenter pointed out to me today that we could define our own Eigen expression template which stored the reflector factors and was able to use them as such, but would naturally .eval() to a normal Eigen matrix