Possible performance issue: cov_matrix_constrain works in row-order on a column-major matrix

Description

This is not a bug, but it might be responsible for a fair number of cache misses if the matrix is large enough. The function in question is stan/math/prim/fun/cov_matrix_constrain.hpp, although the pattern may repeat. This function takes the internal unconstrained representation of a covariance matrix and returns an actual covariance matrix.

The problem is that the key loop looks like this:

  Matrix<value_type_t<T>, Dynamic, Dynamic> L(K, K);
  // stuff
  const Eigen::Ref<const plain_type_t<T>>& x_ref = x;
  int i = 0;
  for (Eigen::Index m = 0; m < K; ++m) {
    L.row(m).head(m) = x_ref.segment(i, m);
    i += m;
    L.coeffRef(m, m) = exp(x_ref.coeff(i++));
    L.row(m).tail(K - m - 1).setZero();
  }

The problem is that the type Matrix<value_type_t<T>, Dynamic, Dynamic> is a matrix stored in column major order (this is the default option in Eigen and the ordering for which they "guarantee" that their algorithms work, so we don't want to change this). This means that when the matrix is flattened to a 1D array internally, it stores the whole of the first column, then the whole of the second column, etc etc.

Hence it is cache efficient to access the matrix by going down the columns and can be quite inefficient to access it by going along the rows. Unfortunately, that loop goes along the rows.

I don't have the time or capacity to check this right now (sorry), but for big models that use covariance matrices, there's probably a good efficiency gain to get from this.

This whole thing is a bit of a mess and, I would guess, at least partly related to an early implementation of the reverse mode AD for the cholesky that followed Giles and used row-major operations (this code is currently only used for small matrices). It is also likely a thing because C/C++ arrays are stored in row-major order. (Here Eigen is following the Fortran standard.) Basically matrices are always a mess.

This means that instead of doing the ldlt decomp as A = L * D * L' (L is lower triangular) the internal parameterisation uses A = U' * D * U (U is upper triangular). This is in all of the positive definite constraints. (search the code for .matrixU() and anything that uses math/stan/math/prim/fun/factor_U.hpp.

This is likely to be a bit of an ordeal to change, so definitely do some speed tests firsts (ie do not ever ever ever trust me). But lower triangular matrices definitely make more sense when working in column major order.

Start by reading stan/math/prim/fun/cov_matrix_constrain.hpp and measuring the current behavior with representative large covariance matrices. Search for .matrixU() and uses of math/stan/math/prim/fun/factor_U.hpp to assess whether the pattern repeats. Done means speed tests establish whether column-major access provides a meaningful gain and clarify the feasible scope without changing behavior.