Overloading signature: pass matrix or array of matrix

Developers
1

The motivation is to allow transition matrices in Stan’s hmm functions which vary over time. So instead of passing a single matrix, the user would pass an array of matrices. We would then have two signatures. See Allowing changing transition matrix for hmm suites · Issue #2678 · stan-dev/math · GitHub.

Under the hood, this means rewriting the code using an std::vector<Eigen::Matrix> type for Gamma. I can then write a wrapper which takes in Gamma as an Eigen::Matrix and converts it to an std::vector<Eigen::Matrix> by duplicating the original matrix. This is of course memory intensive. Is there a way to do this without actually duplicating the matrix, i.e. have an std::vector of pointers, with pointers potentially pointing to the same matrix in the case where the users only passes a single matrix? I’m not sure how to do this with std but it seems Eigen has some features supporting this.

2

First, have you checked out how moveHmm does this? They allow covariates for HMMs (with applications to animal movement models with covariates, like direction of nearest water supply).

Is the issue that you want to use a single signature with potentially varying transition matrices even if they’re all the same? The most efficient way to do that is to define your own interface (a so-called “view”) and then implement it for one or a sequence of covariance matrices. Eigen does this for you with Constant, but you can do it yourself:

template <typename T>
class matrix_view {
  const Eigen::Matrix<T, -1, -1> x_;
  matrix_view(const Eigen::Matrix<T, -1, -1> x) : x_(x) { }
  Eigen::Matrix<T, -1, -1> operator[](int n) {
    return x_;
  }
};

then you can use a std::vector<Eigen::Matrix<T, -1, -1>>, and a matrix_view<T> interchangeably (you just have to template out that argument).

A less elegant solution is to instead use a standard vector of pointers,

std::vector<const Eigen::Matrix<T, -1, -1>*> vs(n);

Then if each reference can be to the same matrix and you’ve only allocated n - 1 pointers beyond what’s strictly necessary.

Here’s a simple example doing this with vectors of vectors (so I don’t have to link Eigen to do it).

#include <cmath>
#include <iostream>
#include <vector>

int main() {
  using std::vector;
  vector<vector<double>> x(10);
  for (int i = 0; i < 10; ++i) {
    auto z = std::vector<double>{ 3.9, 2.7 * i, std::pow(2, i) };
    x[i] = z;
  }

  vector<vector<double>*> y(10);
  for (int i = 0; i < 10; ++i)
    y[i] = &x[i];

  std::cout << "y[2][1] = " << (*y[2])[1] << std::endl;
}