Multi-indexing of matrices (returning a vector)

Summary:

I've put together a quick proof-of-concept for a simple function that could be used for multi-indexing of matrices, but I'm not sure whether it fits in the Stan language or which repo it belongs in.

Description:

Something I've frequently wanted in Stan is the ability to extract specified coefficients of a matrix as a vector, without needing to copy those coefficients to a separate object. For example, something like the following:

matrix[5,4] test_mat = ...;

// Extract coefficients at {0,0}, {1,2}, {1,0}
int ii[3] = {0, 1, 1};
int jj[3] = {0, 2, 0};

// Extract coefficients specified by the ii & jj pair 
vector[3] test_vec = test_mat[ii, jj];

A use-case I have for this is combining vectorised distributions with sparse matrices (due to missingness), so instead of the loop:

for(m in 1:M) {
  y[m] ~ ordered_logistic(ystar[ii[m], jj[m]], thresholds[m]);
}

I could use:

  y ~ ordered_logistic(ystar[ii, jj], thresholds);

Example C++

I've put together an example function that could do this on godbolt: https://godbolt.org/z/fqbj1f

This example function also has the benefit of returning an Eigen expression for the indexing, which would fit nicely in the current ecosystem.

Flaws

A key issue here is that a user might expect ystar[ii, jj] to return a matrix, rather than a vector.

So alternatively, this could be a to_vector overload with the signature:

to_vector(matrix, int[], int[])

In which case this should be added through the math repo, rather than stan.

Current Version:

v2.26.0

Start with the linked Godbolt C++ proof of concept and review the proposed matrix indexing and to_vector overload. Determine whether the API belongs in Stan or the math repository, then define the expected return shape and usage; done means an agreed design, repository scope, implementation plan, and tests.