How to constrain a matrix comprising several parameters to be an orthogonal matrix

Modeling
8

Orthogonal matrices have been discussed previously. There’s even code for building one with Householder reflections. However I don’t think that’s going to work well, mainly because there’s a discontinuity where the matrix flips from positive-determinant to negative-determinant.

For dealing with coordinate systems you probably want a special orthogonal matrix, i.e. one that has positive determinant. This additional constraint makes the matrix group at least connected so there’s some hope that Stan could sample it.

It’s not too difficult to create three orthonormal vectors in three dimensions:

  1. take two arbitrary vectors
  2. project one of them to the plane orthogonal to the other
  3. normalize both
  4. take their cross product to create the third vector

This tranlates to Stan program

parameters {
  row_vector[3] a;
  row_vector[3] b;
}
transformed parameters {
  matrix[3,3] r;
  {
    row_vector[3] ax = a/sqrt(dot_self(a));
    row_vector[3] bx = b - dot_product(b, ax) * ax;
    bx = bx/sqrt(dot_self(bx));
    r[1] = ax;
    r[2] = bx;
    r[3] = [ax[2]*bx[3] - ax[3]*bx[2],
            ax[3]*bx[1] - ax[1]*bx[3],
            ax[1]*bx[2] - ax[2]*bx[1]];
  }
}
model {
  a ~ std_normal();
  b ~ std_normal();
  target += -5*square(dot_product(a,b));
}

The target+= statement is rather arbitrary. It doesn’t affect the distribution of the matrix but keeps the underlying vectors from being too collinear. The parametrization is similar to the unit vector and this other thread discussing problems with unit vectors may be relevant.

3 Likes