Add orthogonal and semi-orthogonal matrix types

Ouverte
  • Quelque chose ne va pas ?
#3,279 3 commentaires 0 réactions 0 personnes assignées Voir sur GitHub

Les mainteneurs répondent en général sous 1 jour

Personne n'a encore pris cette issue.

Évaluation

Difficulté
5/5
Temps estimé
Plus d'une semaine
Accessibilité débutants
35/100
Type d'issue
Fonctionnalité
Clarté
Plutôt claire
Activité
À l'abandon
Stack technique
cpp
Domaine
backend

Piste de recherche

Commencez par lire les discussions liées du forum Stan, la thèse de Michael Jauch et l’implémentation de stan_semiorthogonal_transforms. Comparez les types proposés orthogonal_matrix[K] et semi_orthogonal_matrix[K, N] avec la paramétrisation par réflecteurs et ses règles de comptage des paramètres. Le travail est considéré comme terminé lorsque le projet dispose d’une implémentation convenue couvrant les types demandés ainsi que des opérations efficaces de multiplication, de déterminant et de résolution de systèmes linéaires.

Rédigé par le modèle d'indexation à partir du texte de l'issue.

Description

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

Langage dominant
C++
Étoiles
845
Forks
225
Merge moyen
4 j 21 h
PR mergées (30 j)
13

Préparer son environnement

Par où commencer

  1. Lisez l'issue en entier, puis le guide de contribution du projet.
  2. Signalez en commentaire que vous la prenez — cela évite que deux personnes fassent le même travail.
  3. Forkez le dépôt et travaillez sur une branche.
  4. Ouvrez une pull request qui référence le numéro de l'issue.

Autres issues de stan-dev/math

Toutes les issues de stan-dev/math

Issues similaires

Plus d'issues C++

Recevez les nouvelles issues par e-mail

Un résumé court des issues GitHub adaptées aux débutants.