Add orthogonal and semi-orthogonal matrix types

Nobody has claimed this yet.

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

Start by reading the linked Stan forum discussions, Michael Jauch's thesis, and the stan_semiorthogonal_transforms implementation. Compare the proposed orthogonal_matrix[K] and semi_orthogonal_matrix[K, N] types with the reflector parameterization and its parameter-count rules. Done means the project has an agreed implementation covering the requested types and efficient multiplication, determinant, and linear-solve operations.

Written by the indexing model from the issue text.

Description

Contributor guide

First steps

Research direction

Assessment