Can stan compute Jacobian as a Eigen::SparseMatrix?

Modeling
1

Apologies if this is already answered in the related dev threads…

Is there a way to (efficiently) compute the sparse Jacobian using stan?

Here’s a toy version of what I’m trying to do. First, here’s it working correctly but dense:

#include <cmath>
#include <stan/math.hpp>
#include <vector>

#include <Eigen/Dense>
#include <Eigen/Sparse>

template <typename T>
Eigen::Matrix<T, Eigen::Dynamic, 1> f( const Eigen::Matrix<T, Eigen::Dynamic, 1>& x)
{
  Eigen::SparseMatrix<T> G(10,10);
  std::vector<Eigen::Triplet<T> > ijv;
  ijv.emplace_back(0,0,x(0));
  ijv.emplace_back(1,1,x(1));
  ijv.emplace_back(2,2,x(2));
  ijv.emplace_back(3,3,x(3));
  for(int d = 4;d<10;d++) ijv.emplace_back(d,d,1);
  G.setFromTriplets(ijv.begin(),ijv.end());
  G.makeCompressed();
  return (Eigen::RowVectorXd::Constant(1,10,1) * G).transpose();
}

int main() 
{
  Eigen::VectorXd x = Eigen::VectorXd::Constant(10,1,1);
  Eigen::MatrixXd J;
  Eigen::VectorXd fx;
  stan::math::jacobian( f<stan::math::var>, x, fx, J);
  std::cout<< Eigen::MatrixXd(J) <<std::endl;
  return 0; 
}

This prints:

1 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0 0 0
0 0 0 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0

If I naively try to switch J to SparseMatrix like this:

Eigen::SparseMatrix<double> J;

Then I get this error:

error: no matching function for call to 'jacobian'
  stan::math::jacobian( f<stan::math::var>, x, fx, J);
  ^~~~~~~~~~~~~~~~~~~~
/Users/ajx/Dropbox/math/stan/math/rev/mat/functor/jacobian.hpp:13:6: note: candidate function not viable: no known conversion from
      'Eigen::SparseMatrix<double>' to 'Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> &' for 4th argument
void jacobian(const F& f, const Eigen::Matrix<double, Eigen::Dynamic, 1>& x,
     ^
1 error generated.

Is there a way to compute df/dx efficiently as a sparse matrix?

2

The fast answer is no. Stan does not support sparse matrices (except for sparse matrix-vector products). It also does not compute sparse Jacobians.

3

Thanks for the very fast answer. Is this kind of computation on the roadmap for sparse stuff that is being discussed on the forum? If not, is there a formal place to request it?

4

I don’t believe that our current autodiff library could be easily extended to building sparse Jacobians so I don’t think there are any immediate plans to do it.