Merge pull request #904 from yizhang-cae/feature/892_rename_matrix_exp_action

Feature/892 rename matrix exp action

stan/math/prim/mat.hpp

@@ -158,6 +158,7 @@
 #include <stan/math/prim/mat/fun/logit.hpp>
 #include <stan/math/prim/mat/fun/make_nu.hpp>
 #include <stan/math/prim/mat/fun/matrix_exp.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp_multiply.hpp>
 #include <stan/math/prim/mat/fun/max.hpp>
 #include <stan/math/prim/mat/fun/mdivide_left.hpp>
 #include <stan/math/prim/mat/fun/mdivide_left_ldlt.hpp>
@@ -203,6 +204,7 @@
 #include <stan/math/prim/mat/fun/rows.hpp>
 #include <stan/math/prim/mat/fun/rows_dot_product.hpp>
 #include <stan/math/prim/mat/fun/rows_dot_self.hpp>
+#include <stan/math/prim/mat/fun/scale_matrix_exp_multiply.hpp>
 #include <stan/math/prim/mat/fun/sd.hpp>
 #include <stan/math/prim/mat/fun/segment.hpp>
 #include <stan/math/prim/mat/fun/simplex_constrain.hpp>

stan/math/rev/mat/fun/matrix_exp_action_handler.hpp → stan/math/prim/mat/fun/matrix_exp_action_handler.hpp

-#ifndef STAN_MATH_REV_MAT_FUN_MATRIX_EXP_ACTION_HANDLER_HPP
-#define STAN_MATH_REV_MAT_FUN_MATRIX_EXP_ACTION_HANDLER_HPP
+#ifndef STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_ACTION_HANDLER_HPP
+#define STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_ACTION_HANDLER_HPP
 
 #include <stan/math/prim/mat/fun/Eigen.hpp>
-#include <stan/math/rev/core.hpp>
 #include <vector>
 
-double l1norm(const Eigen::MatrixXd& m) {
-  return m.colwise().lpNorm<1>().maxCoeff();
-}
-
 namespace stan {
 namespace math {
 
@@ -25,8 +20,18 @@ class matrix_exp_action_handler {
   static constexpr int p_max = 8;
   static constexpr int m_max = 55;
   static constexpr double tol = 1.1e-16;
-  static const std::vector<double> theta_m_single_precision;
-  static const std::vector<double> theta_m_double_precision;
+
+  // table 3.1 in the reference
+  const std::vector<double> theta_m_single_precision{
+      1.3e-1, 1.0e0, 2.2e0, 3.6e0, 4.9e0, 6.3e0,
+      7.7e0,  9.1e0, 1.1e1, 1.2e1, 1.3e1};
+  const std::vector<double> theta_m_double_precision{
+      2.4e-3, 1.4e-1, 6.4e-1, 1.4e0, 2.4e0, 3.5e0,
+      4.7e0,  6.0e0,  7.2e0,  8.5e0, 9.9e0};
+
+  double l1norm(const Eigen::MatrixXd& m) {
+    return m.colwise().lpNorm<1>().maxCoeff();
+  }
 
  public:
   /* Constructor
@@ -116,15 +121,6 @@ class matrix_exp_action_handler {
   }
 };
 
-// table 3.1 in the reference
-const std::vector<double> matrix_exp_action_handler::theta_m_single_precision{
-    1.3e-1, 1.0e0, 2.2e0, 3.6e0, 4.9e0, 6.3e0,
-    7.7e0,  9.1e0, 1.1e1, 1.2e1, 1.3e1};
-
-const std::vector<double> matrix_exp_action_handler::theta_m_double_precision{
-    2.4e-3, 1.4e-1, 6.4e-1, 1.4e0, 2.4e0, 3.5e0,
-    4.7e0,  6.0e0,  7.2e0,  8.5e0, 9.9e0};
-
 }  // namespace math
 }  // namespace stan
 

stan/math/prim/mat/fun/matrix_exp_multiply.hpp

+#ifndef STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_MULTIPLY_HPP
+#define STAN_MATH_PRIM_MAT_FUN_MATRIX_EXP_MULTIPLY_HPP
+
+#include <stan/math/prim/mat.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp_action_handler.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return product of exp(A) and B, where A is a NxN double matrix,
+ * B is a NxCb double matrix, and t is a double
+ * @tparam Cb Columns matrix B
+ * @param[in] A Matrix
+ * @param[in] B Matrix
+ * @return exponential of A multiplies B
+ */
+template <int Cb>
+inline Eigen::Matrix<double, -1, Cb> matrix_exp_multiply(
+    const Eigen::MatrixXd& A, const Eigen::Matrix<double, -1, Cb>& B) {
+  check_nonzero_size("scale_matrix_exp_multiply", "input matrix", A);
+  check_nonzero_size("scale_matrix_exp_multiply", "input matrix", B);
+  check_multiplicable("scale_matrix_exp_multiply", "A", A, "B", B);
+  check_square("scale_matrix_exp_multiply", "input matrix", A);
+  return matrix_exp_action_handler().action(A, B);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/prim/mat/fun/scale_matrix_exp_multiply.hpp

+#ifndef STAN_MATH_PRIM_MAT_FUN_SCALE_MATRIX_EXP_MULTIPLY_HPP
+#define STAN_MATH_PRIM_MAT_FUN_SCALE_MATRIX_EXP_MULTIPLY_HPP
+
+#include <stan/math/prim/mat.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp_action_handler.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return product of exp(At) and B, where A is a NxN double matrix,
+ * B is a NxCb double matrix, and t is a double
+ * @tparam Cb Columns matrix B
+ * @param[in] A Matrix
+ * @param[in] B Matrix
+ * @param[in] t double
+ * @return exponential of At multiplies B
+ */
+template <int Cb>
+inline Eigen::Matrix<double, -1, Cb> scale_matrix_exp_multiply(
+    const double& t, const Eigen::MatrixXd& A,
+    const Eigen::Matrix<double, -1, Cb>& B) {
+  check_nonzero_size("scale_matrix_exp_multiply", "input matrix", A);
+  check_nonzero_size("scale_matrix_exp_multiply", "input matrix", B);
+  check_multiplicable("scale_matrix_exp_multiply", "A", A, "B", B);
+  check_square("scale_matrix_exp_multiply", "input matrix", A);
+  return matrix_exp_action_handler().action(A, B, t);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/rev/mat.hpp

@@ -30,6 +30,7 @@
 #include <stan/math/rev/mat/fun/log_determinant_spd.hpp>
 #include <stan/math/rev/mat/fun/log_softmax.hpp>
 #include <stan/math/rev/mat/fun/log_sum_exp.hpp>
+#include <stan/math/rev/mat/fun/matrix_exp_multiply.hpp>
 #include <stan/math/rev/mat/fun/mdivide_left.hpp>
 #include <stan/math/rev/mat/fun/mdivide_left_ldlt.hpp>
 #include <stan/math/rev/mat/fun/mdivide_left_spd.hpp>
@@ -39,6 +40,7 @@
 #include <stan/math/rev/mat/fun/quad_form.hpp>
 #include <stan/math/rev/mat/fun/quad_form_sym.hpp>
 #include <stan/math/rev/mat/fun/rows_dot_product.hpp>
+#include <stan/math/rev/mat/fun/scale_matrix_exp_multiply.hpp>
 #include <stan/math/rev/mat/fun/sd.hpp>
 #include <stan/math/rev/mat/fun/softmax.hpp>
 #include <stan/math/rev/mat/fun/squared_distance.hpp>

stan/math/rev/mat/fun/matrix_exp_action.hpp → stan/math/rev/mat/fun/matrix_exp_multiply.hpp

-#ifndef STAN_MATH_REV_MAT_FUN_MATRIX_EXP_ACTION_HPP
-#define STAN_MATH_REV_MAT_FUN_MATRIX_EXP_ACTION_HPP
+#ifndef STAN_MATH_REV_MAT_FUN_MATRIX_EXP_MULTIPLY_HPP
+#define STAN_MATH_REV_MAT_FUN_MATRIX_EXP_MULTIPLY_HPP
 
+#include <stan/math/rev/mat.hpp>
 #include <stan/math/rev/core.hpp>
-#include <stan/math/rev/mat/fun/matrix_exp_action_handler.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp_action_handler.hpp>
 #include <stan/math/prim/mat/fun/Eigen.hpp>
 #include <stan/math/rev/mat/fun/to_var.hpp>
 #include <stan/math/prim/mat/fun/matrix_exp.hpp>
@@ -20,47 +21,48 @@ namespace math {
  * Calculate adjoint of matrix exponential action
  * exp(At)*B when A is data and B is var, used for chaining action.
  *
- * @tparam N integer rows and cols of matrix A
  * @param Ad double array pointer to the data of matrix A
+ * @param n dim of square matrix A
  * @param adjexpAB MatrixXd adjoint of exp(At)*B
  * @param t double time data
  * @return MatrixXd The adjoint of B.
  */
-template <int N>
-Eigen::MatrixXd exp_action_chain_dv(double* Ad, const Eigen::MatrixXd& adjexpAB,
-                                    const double t) {
+inline Eigen::MatrixXd exp_action_chain_dv(double* Ad, const int& n,
+                                           const Eigen::MatrixXd& adjexpAB,
+                                           const double t) {
   using Eigen::Map;
   matrix_exp_action_handler handle;
-  return handle.action(Map<Eigen::MatrixXd>(Ad, N, N).transpose(), adjexpAB, t);
+  return handle.action(Map<Eigen::MatrixXd>(Ad, n, n).transpose(), adjexpAB, t);
 }
 
 /**
  * Calculate adjoint of matrix exponential action
  * exp(At)*B when A is var and B is data, used for chaining action.
  *
- * @tparam N integer rows and cols of matrix A
  * @param Ad double array pointer to the data of matrix A
  * @param Bd double array pointer to the data of matrix B
+ * @param n dim(nb. of rows) of square matrix A
+ * @param m nb. of cols of matrix B
  * @param adjexpAB MatrixXd adjoint of exp(At)*B
  * @param t double time data
  * @return MatrixXd The adjoint of A.
  */
-template <int N, int M>
-Eigen::MatrixXd exp_action_chain_vd(double* Ad, double* Bd,
-                                    const Eigen::MatrixXd& adjexpAB,
-                                    const double t) {
+inline Eigen::MatrixXd exp_action_chain_vd(double* Ad, double* Bd, const int& n,
+                                           const int& m,
+                                           const Eigen::MatrixXd& adjexpAB,
+                                           const double t) {
   using Eigen::Map;
   using Eigen::Matrix;
   using Eigen::MatrixXd;
-  Eigen::MatrixXd adjexpA = Eigen::MatrixXd::Zero(N, N);
-  Eigen::MatrixXd adjA = Eigen::MatrixXd::Zero(N, N);
+  Eigen::MatrixXd adjexpA = Eigen::MatrixXd::Zero(n, n);
+  Eigen::MatrixXd adjA = Eigen::MatrixXd::Zero(n, n);
 
   // TODO(yizhang): a better way such as complex step approximation
   try {
     start_nested();
 
-    adjexpA = adjexpAB * Map<Matrix<double, N, M> >(Bd).transpose();
-    Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> Av(N, N);
+    adjexpA = adjexpAB * Map<MatrixXd>(Bd, n, m).transpose();
+    Eigen::Matrix<stan::math::var, Eigen::Dynamic, Eigen::Dynamic> Av(n, n);
     for (int i = 0; i < Av.size(); ++i) {
       Av(i) = to_var(Ad[i]);
     }
@@ -163,8 +165,8 @@ class matrix_exp_action_vari : public vari {
     for (size_type i = 0; i < adjexpAB.size(); ++i)
       adjexpAB(i) = variRefexpAB_[i]->adj_;
 
-    MatrixXd adjA = exp_action_chain_vd<N, Cb>(Ad_, Bd_, adjexpAB, t_);
-    MatrixXd adjB = exp_action_chain_dv<N>(Ad_, adjexpAB, t_);
+    MatrixXd adjA = exp_action_chain_vd(Ad_, Bd_, n_, B_cols_, adjexpAB, t_);
+    MatrixXd adjB = exp_action_chain_dv(Ad_, n_, adjexpAB, t_);
 
     for (size_type i = 0; i < A_size_; ++i) {
       variRefA_[i]->adj_ += adjA(i);
@@ -261,7 +263,7 @@ class matrix_exp_action_vari<double, N, Tb, Cb> : public vari {
     for (size_type i = 0; i < adjexpAB.size(); ++i)
       adjexpAB(i) = variRefexpAB_[i]->adj_;
 
-    MatrixXd adjB = exp_action_chain_dv<N>(Ad_, adjexpAB, t_);
+    MatrixXd adjB = exp_action_chain_dv(Ad_, n_, adjexpAB, t_);
 
     for (size_type i = 0; i < B_size_; ++i) {
       variRefB_[i]->adj_ += adjB(i);
@@ -352,7 +354,7 @@ class matrix_exp_action_vari<Ta, N, double, Cb> : public vari {
     for (size_type i = 0; i < adjexpAB.size(); ++i)
       adjexpAB(i) = variRefexpAB_[i]->adj_;
 
-    MatrixXd adjA = exp_action_chain_vd<N, Cb>(Ad_, Bd_, adjexpAB, t_);
+    MatrixXd adjA = exp_action_chain_vd(Ad_, Bd_, n_, B_cols_, adjexpAB, t_);
 
     for (size_type i = 0; i < A_size_; ++i) {
       variRefA_[i]->adj_ += adjA(i);
@@ -388,23 +390,23 @@ matrix_exp_action(const Eigen::Matrix<Ta, N, N>& A,
 }
 
 /**
- * Return product of exp(At) and B, where A is a NxN double matrix,
- * B is a NxCb double matrix, and t is a double
- * @tparam N Rows and cols matrix A, also rows of matrix B
+ * Wrapper of matrix_exp_action function for a more literal name
+ * @tparam Ta scalar type matrix A
+ * @tparam Tb scalar type matrix B
  * @tparam Cb Columns matrix B
  * @param[in] A Matrix
  * @param[in] B Matrix
- * @param[in] t double
- * @return exponential of At multiplies B
+ * @return exponential of A multiplies B
  */
-template <int N, int Cb>
-inline Eigen::Matrix<double, N, Cb> matrix_exp_action(
-    const Eigen::Matrix<double, N, N>& A, const Eigen::Matrix<double, N, Cb>& B,
-    const double& t = 1.0) {
-  Eigen::Matrix<double, N, Cb> expAB;
-  matrix_exp_action_handler handle;
-  expAB = handle.action(A, B, t);
-  return expAB;
+template <typename Ta, typename Tb, int Cb>
+inline Eigen::Matrix<typename stan::return_type<Ta, Tb>::type, -1, Cb>
+matrix_exp_multiply(const Eigen::Matrix<Ta, -1, -1>& A,
+                    const Eigen::Matrix<Tb, -1, Cb>& B) {
+  check_nonzero_size("matrix_exp_multiply", "input matrix", A);
+  check_nonzero_size("matrix_exp_multiply", "input matrix", B);
+  check_multiplicable("matrix_exp_multiply", "A", A, "B", B);
+  check_square("matrix_exp_multiply", "input matrix", A);
+  return matrix_exp_action(A, B);
 }
 
 }  // namespace math

stan/math/rev/mat/fun/scale_matrix_exp_multiply.hpp

+#ifndef STAN_MATH_REV_MAT_FUN_SCALE_MATRIX_EXP_MULTIPLY_HPP
+#define STAN_MATH_REV_MAT_FUN_SCALE_MATRIX_EXP_MULTIPLY_HPP
+
+#include <stan/math/prim/mat.hpp>
+#include <stan/math/rev/mat/fun/matrix_exp_multiply.hpp>
+#include <stan/math/prim/mat/fun/Eigen.hpp>
+#include <boost/math/tools/promotion.hpp>
+#include <boost/utility/enable_if.hpp>
+#include <boost/type_traits.hpp>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Wrapper of matrix_exp_action function for a more literal name.
+ * @tparam Ta scalar type matrix A
+ * @tparam Tb scalar type matrix B
+ * @tparam Cb Columns matrix B
+ * @param[in] A Matrix
+ * @param[in] B Matrix
+ * @param[in] t double
+ * @return exponential of At multiplies B
+ */
+template <typename Ta, typename Tb, int Cb>
+inline Eigen::Matrix<typename stan::return_type<Ta, Tb>::type, -1, Cb>
+scale_matrix_exp_multiply(const double& t, const Eigen::Matrix<Ta, -1, -1>& A,
+                          const Eigen::Matrix<Tb, -1, Cb>& B) {
+  check_nonzero_size("scale_matrix_exp_multiply", "input matrix", A);
+  check_nonzero_size("scale_matrix_exp_multiply", "input matrix", B);
+  check_multiplicable("scale_matrix_exp_multiply", "A", A, "B", B);
+  check_square("scale_matrix_exp_multiply", "input matrix", A);
+  return matrix_exp_action(A, B, t);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

test/unit/math/prim/mat/fun/matrix_exp_multiply_test.cpp

+#include <gtest/gtest.h>
+#include <test/unit/math/prim/mat/util.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp_multiply.hpp>
+#include <vector>
+
+template <int N, int M>
+inline void test_matrix_exp_multiply() {
+  using Eigen::Dynamic;
+  using Eigen::Matrix;
+  using stan::math::matrix_exp_multiply;
+
+  std::srand(1999);
+
+  Eigen::MatrixXd A = Eigen::MatrixXd::Random(N, N);
+  Eigen::Matrix<double, -1, M> B = Eigen::MatrixXd::Random(N, M);
+  Eigen::Matrix<double, Dynamic, Dynamic> A0 = A;
+
+  // brute force
+  Eigen::Matrix<double, N, M> expAB
+      = stan::math::multiply(stan::math::matrix_exp(A0), B);
+
+  // matrix_exp_multiply
+  Eigen::Matrix<double, N, M> res = matrix_exp_multiply(A, B);
+  EXPECT_EQ(res.size(), expAB.size());
+  for (int l = 0; l < res.size(); ++l) {
+    EXPECT_FLOAT_EQ(res(l), expAB(l));
+  }
+}
+
+TEST(MathMatrix, matrix_exp_multiply) {
+  test_matrix_exp_multiply<1, 1>();
+  test_matrix_exp_multiply<1, 5>();
+  test_matrix_exp_multiply<5, 1>();
+  test_matrix_exp_multiply<5, 5>();
+  test_matrix_exp_multiply<20, 2>();
+}
+
+TEST(MathMatrix, matrix_exp_multiply_exception) {
+  using stan::math::matrix_exp_multiply;
+  {  // nonzero size
+    Eigen::MatrixXd A(0, 0);
+    Eigen::MatrixXd B = Eigen::MatrixXd::Random(1, 2);
+    EXPECT_THROW(matrix_exp_multiply(A, B), std::invalid_argument);
+    EXPECT_THROW(matrix_exp_multiply(B, A), std::invalid_argument);
+  }
+
+  {  // multiplicable
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(2, 2);
+    Eigen::MatrixXd B = Eigen::MatrixXd::Random(3, 2);
+    EXPECT_THROW(matrix_exp_multiply(A, B), std::invalid_argument);
+  }
+
+  {  // square
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(2, 3);
+    Eigen::MatrixXd B = Eigen::MatrixXd::Random(3, 2);
+    EXPECT_THROW(matrix_exp_multiply(A, B), std::invalid_argument);
+  }
+}

test/unit/math/prim/mat/fun/scale_matrix_exp_multiply_test.cpp

+#include <gtest/gtest.h>
+#include <test/unit/math/prim/mat/util.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp.hpp>
+#include <stan/math/prim/mat/fun/scale_matrix_exp_multiply.hpp>
+#include <vector>
+
+template <int N, int M>
+inline void test_scale_matrix_exp_multiply() {
+  using Eigen::Dynamic;
+  using Eigen::Matrix;
+  using stan::math::scale_matrix_exp_multiply;
+
+  std::srand(1999);
+
+  Eigen::MatrixXd A = Eigen::MatrixXd::Random(N, N);
+  Eigen::Matrix<double, -1, M> B = Eigen::MatrixXd::Random(N, M);
+  Eigen::Matrix<double, Dynamic, Dynamic> A0 = A;
+
+  // brute force
+  Eigen::Matrix<double, N, M> expAB
+      = stan::math::multiply(stan::math::matrix_exp(A0), B);
+
+  // matrix_exp_multiply
+  const double t = 1.0;
+  Eigen::Matrix<double, N, M> res = scale_matrix_exp_multiply(t, A, B);
+  EXPECT_EQ(res.size(), expAB.size());
+  for (int l = 0; l < res.size(); ++l) {
+    EXPECT_FLOAT_EQ(res(l), expAB(l));
+  }
+}
+
+TEST(MathMatrix, scale_matrix_exp_multiply) {
+  test_scale_matrix_exp_multiply<1, 1>();
+  test_scale_matrix_exp_multiply<1, 5>();
+  test_scale_matrix_exp_multiply<5, 1>();
+  test_scale_matrix_exp_multiply<5, 5>();
+  test_scale_matrix_exp_multiply<20, 2>();
+}
+
+TEST(MathMatrix, scale_matrix_exp_multiply_exception) {
+  using stan::math::scale_matrix_exp_multiply;
+  const double t = 1.0;
+  {  // nonzero size
+    Eigen::MatrixXd A(0, 0);
+    Eigen::MatrixXd B = Eigen::MatrixXd::Random(1, 2);
+    EXPECT_THROW(scale_matrix_exp_multiply(t, A, B), std::invalid_argument);
+    EXPECT_THROW(scale_matrix_exp_multiply(t, B, A), std::invalid_argument);
+  }
+
+  {  // multiplicable
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(2, 2);
+    Eigen::MatrixXd B = Eigen::MatrixXd::Random(3, 2);
+    EXPECT_THROW(scale_matrix_exp_multiply(t, A, B), std::invalid_argument);
+  }
+
+  {  // square
+    Eigen::MatrixXd A = Eigen::MatrixXd::Random(2, 3);
+    Eigen::MatrixXd B = Eigen::MatrixXd::Random(3, 2);
+    EXPECT_THROW(scale_matrix_exp_multiply(t, A, B), std::invalid_argument);
+  }
+}

test/unit/math/rev/mat/fun/matrix_exp_action_handler_test.cpp

+#include <gtest/gtest.h>
+// #include <stan/math/prim/scal/err/check_not_nan.hpp>
+#include <test/unit/math/prim/mat/util.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp.hpp>
+#include <stan/math/prim/mat/fun/matrix_exp_action_handler.hpp>
+#include <vector>
+
+TEST(MathMatrix, matrix_exp_action_diag) {
+  using Eigen::MatrixXd;
+  using Eigen::VectorXd;
+  stan::math::matrix_exp_action_handler handler;
+
+  {
+    double t = 0.5;
+    MatrixXd m1(2, 2);
+    VectorXd b(2);
+    m1 << 2, 0, 0, 2;
+    b << 1, 1;
+    auto res = handler.action(m1, b, t);
+    EXPECT_NEAR(res(0), M_E, 1.e-8);
+    EXPECT_NEAR(res(1), M_E, 1.e-8);
+  }
+
+  {
+    double t = 1.0;
+    MatrixXd m1(2, 2);
+    VectorXd b = VectorXd::Random(2);
+    m1 << 1, 0, 0, 2;
+    auto res = handler.action(m1, b, t);
+    EXPECT_NEAR(res(0), b(0) * M_E, 1.e-8);
+    EXPECT_NEAR(res(1), b(1) * M_E * M_E, 1.e-8);
+  }
+
+  {
+    double t = 1.0;
+    std::srand(1299);
+    MatrixXd m1(2, 2);
+    VectorXd b = VectorXd::Random(2);
+    m1 << -4.0, 0, 0, -5.0;
+    auto res = handler.action(m1, b, t);
+    EXPECT_NEAR(res(0), b(0) / (M_E * M_E * M_E * M_E), 1.e-8);
+    EXPECT_NEAR(res(1), b(1) / (M_E * M_E * M_E * M_E * M_E), 1.e-8);
+  }
+
+  {
+    std::srand(999);
+    double t = static_cast<double>((std::rand()) / RAND_MAX);
+    VectorXd b = VectorXd::Random(5);
+    VectorXd d = VectorXd::Random(5);
+    MatrixXd m = d.asDiagonal();
+    auto res = handler.action(m, b, t);
+    for (int i = 0; i < 5; ++i) {
+      EXPECT_NEAR(res(i), b(i) * std::exp(t * d(i)), 1.e-8);
+    }
+  }
+}
+
+TEST(MathMatrix, matrix_exp_action_vector) {
+  using Eigen::MatrixXd;
+  using Eigen::VectorXd;
+  stan::math::matrix_exp_action_handler handler;
+  std::srand(999);
+
+  for (size_t n = 2; n < 10; ++n) {
+    MatrixXd A = MatrixXd::Random(n, n);
+    VectorXd b = VectorXd::Random(n);
+    VectorXd res = handler.action(A, b);
+    MatrixXd expA = stan::math::matrix_exp(A);
+    VectorXd expb = expA * b;
+    for (size_t i = 0; i < n; ++i) {
+      EXPECT_NEAR(res(i), expb(i), 1.e-6);
+    }
+
+    int m1, s1, m2, s2;
+    const double t1 = 9.9, t2 = 1.0;
+    handler.set_approximation_parameter(A, t1, m1, s1);
+    A *= t1;
+    handler.set_approximation_parameter(A, t2, m2, s2);
+    EXPECT_EQ(m1, m2);
+    EXPECT_EQ(s1, s2);
+  }
+}
+
+TEST(MathMatrix, matrix_exp_action_matrix) {
+  using Eigen::MatrixXd;
+  stan::math::matrix_exp_action_handler handler;
+  std::srand(999);
+
+  constexpr int N = 10;
+  constexpr int M = 4;
+  Eigen::Matrix<double, N, N> A = Eigen::Matrix<double, N, N>::Random();
+  Eigen::Matrix<double, N, M> B = Eigen::Matrix<double, N, M>::Random();
+
+  Eigen::Matrix<double, N, M> res = handler.action(A, B);
+  MatrixXd Ad(A);
+  MatrixXd expa = stan::math::matrix_exp(Ad);
+  Eigen::Matrix<double, N, M> expb = expa * B;
+
+  for (int i = 0; i < N; ++i) {
+    for (int j = 0; j < M; ++j) {
+      EXPECT_FLOAT_EQ(res(i, j), expb(i, j));
+    }
+  }
+}
+
+TEST(MathMatrix, matrix_exp_action_matrix_transpose) {
+  using Eigen::MatrixXd;
+  stan::math::matrix_exp_action_handler handler;
+  std::srand(1999);
+
+  constexpr int N = 10;
+  constexpr int M = 4;
+  Eigen::Matrix<double, N, N> A = Eigen::Matrix<double, N, N>::Random();
+  Eigen::Matrix<double, N, M> B = Eigen::Matrix<double, N, M>::Random();
+
+  Eigen::Matrix<double, N, M> res = handler.action(A.transpose(), B);
+  MatrixXd Ad(A);
+  MatrixXd expa = stan::math::matrix_exp(Ad).transpose();
+  Eigen::Matrix<double, N, M> expb = expa * B;
+
+  for (int i = 0; i < N; ++i) {
+    for (int j = 0; j < M; ++j) {
+      EXPECT_NEAR(res(i, j), expb(i, j), 1.e-6);
+    }
+  }
+}

test/unit/math/rev/mat/fun/matrix_exp_action_test.cpp → test/unit/math/rev/mat/fun/matrix_exp_multiply_test.cpp

@@ -3,161 +3,33 @@
 // #include <stan/math/prim/scal/err/check_not_nan.hpp>
 #include <test/unit/math/rev/mat/fun/util.hpp>
 #include <test/unit/math/rev/mat/util.hpp>
-#include <stan/math/rev/mat/fun/matrix_exp_action_handler.hpp>
-#include <stan/math/rev/mat/fun/matrix_exp_action.hpp>
+#include <stan/math/rev/mat/fun/matrix_exp_multiply.hpp>
 #include <stan/math/rev/mat/fun/to_var.hpp>
 #include <vector>
 
-TEST(MathMatrix, matrix_exp_action_diag) {
-  using MatrixType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
-  using VectorType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
-  stan::math::matrix_exp_action_handler handler;
-
-  {
-    double t = 0.5;
-    MatrixType m1(2, 2);
-    VectorType b(2);
-    m1 << 2, 0, 0, 2;
-    b << 1, 1;
-    auto res = handler.action(m1, b, t);
-    EXPECT_NEAR(res(0), M_E, 1.e-8);
-    EXPECT_NEAR(res(1), M_E, 1.e-8);
-  }
-
-  {
-    double t = 1.0;
-    MatrixType m1(2, 2);
-    VectorType b = VectorType::Random(2);
-    m1 << 1, 0, 0, 2;
-    auto res = handler.action(m1, b, t);
-    EXPECT_NEAR(res(0), b(0) * M_E, 1.e-8);
-    EXPECT_NEAR(res(1), b(1) * M_E * M_E, 1.e-8);
-  }
-
-  {
-    double t = 1.0;
-    std::srand(1299);
-    MatrixType m1(2, 2);
-    VectorType b = VectorType::Random(2);
-    m1 << -4.0, 0, 0, -5.0;
-    auto res = handler.action(m1, b, t);
-    EXPECT_NEAR(res(0), b(0) / (M_E * M_E * M_E * M_E), 1.e-8);
-    EXPECT_NEAR(res(1), b(1) / (M_E * M_E * M_E * M_E * M_E), 1.e-8);
-  }
-
-  {
-    std::srand(999);
-    double t = static_cast<double>((std::rand()) / RAND_MAX);
-    VectorType b = VectorType::Random(5);
-    VectorType d = VectorType::Random(5);
-    MatrixType m = d.asDiagonal();
-    auto res = handler.action(m, b, t);
-    for (int i = 0; i < 5; ++i) {
-      EXPECT_NEAR(res(i), b(i) * std::exp(t * d(i)), 1.e-8);
-    }
-  }
-}
-
-TEST(MathMatrix, matrix_exp_action_vector) {
-  using MatrixType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
-  using VectorType = Eigen::Matrix<double, Eigen::Dynamic, 1>;
-  stan::math::matrix_exp_action_handler handler;
-  std::srand(999);
-
-  for (size_t n = 2; n < 10; ++n) {
-    MatrixType A = MatrixType::Random(n, n);
-    VectorType b = VectorType::Random(n);
-    VectorType res = handler.action(A, b);
-    MatrixType expA = stan::math::matrix_exp(A);
-    VectorType expb = expA * b;
-    for (size_t i = 0; i < n; ++i) {
-      EXPECT_NEAR(res(i), expb(i), 1.e-6);
-    }
-
-    int m1, s1, m2, s2;
-    const double t1 = 9.9, t2 = 1.0;
-    handler.set_approximation_parameter(A, t1, m1, s1);
-    A *= t1;
-    handler.set_approximation_parameter(A, t2, m2, s2);
-    EXPECT_EQ(m1, m2);
-    EXPECT_EQ(s1, s2);
-  }
-}
-
-TEST(MathMatrix, matrix_exp_action_matrix) {
-  using MatrixType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
-  stan::math::matrix_exp_action_handler handler;
-  std::srand(999);
-
-  constexpr int N = 10;
-  constexpr int M = 4;
-  Eigen::Matrix<double, N, N> A = Eigen::Matrix<double, N, N>::Random();
-  Eigen::Matrix<double, N, M> B = Eigen::Matrix<double, N, M>::Random();
-
-  Eigen::Matrix<double, N, M> res = handler.action(A, B);
-  MatrixType Ad(A);
-  MatrixType expa = stan::math::matrix_exp(Ad);
-  Eigen::Matrix<double, N, M> expb = expa * B;
-
-  for (int i = 0; i < N; ++i) {
-    for (int j = 0; j < M; ++j) {
-      EXPECT_FLOAT_EQ(res(i, j), expb(i, j));
-    }
-  }
-}
-
-TEST(MathMatrix, matrix_exp_action_matrix_transpose) {
-  using MatrixType = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>;
-  stan::math::matrix_exp_action_handler handler;
-  std::srand(1999);
-
-  constexpr int N = 10;
-  constexpr int M = 4;
-  Eigen::Matrix<double, N, N> A = Eigen::Matrix<double, N, N>::Random();
-  Eigen::Matrix<double, N, M> B = Eigen::Matrix<double, N, M>::Random();
-
-  Eigen::Matrix<double, N, M> res = handler.action(A.transpose(), B);
-  MatrixType Ad(A);
-  MatrixType expa = stan::math::matrix_exp(Ad).transpose();
-  Eigen::Matrix<double, N, M> expb = expa * B;
-
-  for (int i = 0; i < N; ++i) {
-    for (int j = 0; j < M; ++j) {
-      EXPECT_NEAR(res(i, j), expb(i, j), 1.e-6);
-    }
-  }
-}
-
-template <int N, int M>
-inline void test_matrix_exp_action_dv() {
+inline void test_matrix_exp_multiply_dv(int N, int M) {
   using stan::math::value_of;
   using stan::math::var;
-  stan::math::matrix_exp_action_handler handler;
+
   std::srand(1999);
 
-  Eigen::Matrix<var, N, N> Av = Eigen::Matrix<var, N, N>::Random();
-  Eigen::Matrix<var, N, M> Bv = Eigen::Matrix<var, N, M>::Random();
-  std::vector<stan::math::var> Bvec(Bv.data(), Bv.data() + Bv.size());
-  std::vector<stan::math::var> Avec(Av.data(), Av.data() + Av.size());
-  Eigen::Matrix<double, N, N> A = value_of(Av);
+  Eigen::Matrix<var, -1, -1> Av = Eigen::Matrix<var, -1, -1>::Random(N, N);
+  Eigen::Matrix<var, -1, -1> Bv = Eigen::Matrix<var, -1, -1>::Random(N, M);
+  std::vector<stan::math::var> Avec = stan::math::to_array_1d(Av);
+  std::vector<stan::math::var> Bvec = stan::math::to_array_1d(Bv);
+  Eigen::MatrixXd A = value_of(Av);
 
   // brute force
-  Eigen::Matrix<double, N, N> expA = stan::math::matrix_exp(A);
-  Eigen::Matrix<var, N, M> expAB;
-  for (int i = 0; i < N; ++i) {
-    for (int j = 0; j < M; ++j) {
-      expAB(i, j) = 0.0;
-      for (int k = 0; k < N; ++k) {
-        expAB(i, j) += expA(i, k) * Bv(k, j);
-      }
-    }
-  }
+  Eigen::Matrix<var, -1, -1> expAB
+      = stan::math::multiply(stan::math::matrix_exp(A), Bv);
 
-  // matrix_exp_action
-  Eigen::Matrix<var, N, M> res_dv = stan::math::matrix_exp_action(A, Bv);
+  // matrix_exp_multiply
+  Eigen::Matrix<var, -1, -1> res_dv = stan::math::matrix_exp_multiply(A, Bv);
+  EXPECT_EQ(res_dv.size(), expAB.size());
   for (int l = 0; l < res_dv.size(); ++l) {
     EXPECT_FLOAT_EQ(res_dv(l).val(), expAB(l).val());
   }
+
   // compare adjoints
   std::vector<double> g, g0;
   for (int l = 0; l < M; ++l) {
@@ -184,39 +56,32 @@ inline void test_matrix_exp_action_dv() {
   }
 }
 
-TEST(MathMatrix, matrix_exp_action_dv) {
-  test_matrix_exp_action_dv<1, 1>();
-  test_matrix_exp_action_dv<1, 5>();
-  test_matrix_exp_action_dv<5, 1>();
-  test_matrix_exp_action_dv<5, 5>();
+TEST(MathMatrix, matrix_exp_multiply_dv) {
+  test_matrix_exp_multiply_dv(1, 1);
+  test_matrix_exp_multiply_dv(1, 5);
+  test_matrix_exp_multiply_dv(5, 1);
+  test_matrix_exp_multiply_dv(5, 5);
 }
 
-template <int N, int M>
-inline void test_matrix_exp_action_vd() {
+inline void test_matrix_exp_multiply_vd(int N, int M) {
   using stan::math::value_of;
   using stan::math::var;
-  stan::math::matrix_exp_action_handler handler;
+
   std::srand(1999);
 
-  Eigen::Matrix<var, N, N> Av = Eigen::Matrix<var, N, N>::Random();
-  Eigen::Matrix<var, N, M> Bv = Eigen::Matrix<var, N, M>::Random();
-  std::vector<stan::math::var> Bvec(Bv.data(), Bv.data() + Bv.size());
-  std::vector<stan::math::var> Avec(Av.data(), Av.data() + Av.size());
-  Eigen::Matrix<double, N, M> B = value_of(Bv);
+  Eigen::Matrix<var, -1, -1> Av = Eigen::Matrix<var, -1, -1>::Random(N, N);
+  Eigen::Matrix<var, -1, -1> Bv = Eigen::Matrix<var, -1, -1>::Random(N, M);
+  std::vector<stan::math::var> Avec = stan::math::to_array_1d(Av);
+  std::vector<stan::math::var> Bvec = stan::math::to_array_1d(Bv);
+  Eigen::MatrixXd B = value_of(Bv);
 
   // brute force
-  Eigen::Matrix<var, N, N> expA = stan::math::matrix_exp(Av);
-  Eigen::Matrix<var, N, M> expAB;
-  for (int i = 0; i < N; ++i) {
-    for (int j = 0; j < M; ++j) {
-      expAB(i, j) = 0.0;
-      for (int k = 0; k < N; ++k) {
-        expAB(i, j) += expA(i, k) * B(k, j);
-      }
-    }
-  }
-  // matrix_exp_action
-  Eigen::Matrix<var, N, M> res_vd = stan::math::matrix_exp_action(Av, B);
+  Eigen::Matrix<var, -1, -1> expAB
+      = stan::math::multiply(stan::math::matrix_exp(Av), B);
+
+  // matrix_exp_multiply
+  Eigen::Matrix<var, -1, -1> res_vd = stan::math::matrix_exp_multiply(Av, B);
+  EXPECT_EQ(res_vd.size(), expAB.size());
   for (int l = 0; l < res_vd.size(); ++l) {
     EXPECT_FLOAT_EQ(res_vd(l).val(), expAB(l).val());
   }
@@ -247,44 +112,37 @@ inline void test_matrix_exp_action_vd() {
   }
 }
 
-TEST(MathMatrix, matrix_exp_action_vd) {
-  test_matrix_exp_action_vd<1, 1>();
-  test_matrix_exp_action_vd<1, 5>();
-  test_matrix_exp_action_vd<5, 1>();
-  test_matrix_exp_action_vd<5, 5>();
+TEST(MathMatrix, matrix_exp_multiply_vd) {
+  test_matrix_exp_multiply_vd(1, 1);
+  test_matrix_exp_multiply_vd(1, 5);
+  test_matrix_exp_multiply_vd(5, 1);
+  test_matrix_exp_multiply_vd(5, 5);
 }
 
-template <int N, int M>
-inline void test_matrix_exp_action_vv() {
+inline void test_matrix_exp_multiply_vv(int N, int M) {
   using stan::math::value_of;
   using stan::math::var;
-  stan::math::matrix_exp_action_handler handler;
+
   std::srand(2999);
 
-  Eigen::Matrix<var, N, N> Av = Eigen::Matrix<var, N, N>::Random();
-  Eigen::Matrix<var, N, M> Bv = Eigen::Matrix<var, N, M>::Random();
-  std::vector<stan::math::var> Bvec(Bv.data(), Bv.data() + Bv.size());
-  std::vector<stan::math::var> Avec(Av.data(), Av.data() + Av.size());
+  Eigen::Matrix<var, -1, -1> Av = Eigen::Matrix<var, -1, -1>::Random(N, N);
+  Eigen::Matrix<var, -1, -1> Bv = Eigen::Matrix<var, -1, -1>::Random(N, M);
+  std::vector<stan::math::var> Avec = stan::math::to_array_1d(Av);
+  std::vector<stan::math::var> Bvec = stan::math::to_array_1d(Bv);
 
   // brute force
-  Eigen::Matrix<var, N, N> expA = stan::math::matrix_exp(Av);
-  Eigen::Matrix<var, N, M> expAB;
-  for (int i = 0; i < N; ++i) {
-    for (int j = 0; j < M; ++j) {
-      expAB(i, j) = 0.0;
-      for (int k = 0; k < N; ++k) {
-        expAB(i, j) += expA(i, k) * Bv(k, j);
-      }
-    }
-  }
+  Eigen::Matrix<var, -1, -1> expAB
+      = stan::math::multiply(stan::math::matrix_exp(Av), Bv);
 
-  // matrix_exp_action
-  Eigen::Matrix<var, N, M> res_vv = stan::math::matrix_exp_action(Av, Bv);
+  // matrix_exp_multiply
+  Eigen::Matrix<var, -1, -1> res_vv = stan::math::matrix_exp_multiply(Av, Bv);
+  EXPECT_EQ(res_vv.size(), expAB.size());
   for (int l = 0; l < res_vv.size(); ++l) {
     EXPECT_FLOAT_EQ(res_vv(l).val(), expAB(l).val());
   }
-  Avec.insert(Avec.end(), Bvec.begin(), Bvec.end());
+
   // compare adjoints
+  Avec.insert(Avec.end(), Bvec.begin(), Bvec.end());
   std::vector<double> g, g0;
   for (int l = 0; l < M; ++l) {
     for (int k = 0; k < N; ++k) {
@@ -310,19 +168,41 @@ inline void test_matrix_exp_action_vv() {
   }
 }
 
-TEST(MathMatrix, matrix_exp_action_vv) {
-  test_matrix_exp_action_vv<1, 1>();
-  test_matrix_exp_action_vv<1, 5>();
-  test_matrix_exp_action_vv<5, 1>();
-  test_matrix_exp_action_vv<5, 5>();
-  test_matrix_exp_action_vv<10, 2>();
+TEST(MathMatrix, matrix_exp_multiply_vv) {
+  test_matrix_exp_multiply_vv(1, 1);
+  test_matrix_exp_multiply_vv(1, 5);
+  test_matrix_exp_multiply_vv(5, 1);
+  test_matrix_exp_multiply_vv(5, 5);
+  test_matrix_exp_multiply_vv(8, 2);
+}
+
+TEST(MathMatrix, matrix_exp_multiply_exception) {
+  using stan::math::matrix_exp_multiply;
+  using stan::math::var;
+  {  // nonzero size
+    Eigen::Matrix<var, -1, -1> A(0, 0);
+    Eigen::Matrix<var, -1, -1> B = Eigen::Matrix<var, -1, -1>::Random(1, 2);
+    EXPECT_THROW(matrix_exp_multiply(A, B), std::invalid_argument);
+    EXPECT_THROW(matrix_exp_multiply(B, A), std::invalid_argument);
+  }
+
+  {  // multiplicable
+    Eigen::Matrix<var, -1, -1> A = Eigen::Matrix<var, -1, -1>::Random(2, 2);
+    Eigen::Matrix<var, -1, -1> B = Eigen::Matrix<var, -1, -1>::Random(3, 2);
+    EXPECT_THROW(matrix_exp_multiply(A, B), std::invalid_argument);
+  }
+
+  {  // square
+    Eigen::Matrix<var, -1, -1> A = Eigen::Matrix<var, -1, -1>::Random(2, 3);
+    Eigen::Matrix<var, -1, -1> B = Eigen::Matrix<var, -1, -1>::Random(3, 2);
+    EXPECT_THROW(matrix_exp_multiply(A, B), std::invalid_argument);
+  }
 }
 
-// TEST(MathMatrix, matrix_exp_action_clock) {
+// TEST(MathMatrix, matrix_exp_multiply_clock) {
 //   using Eigen::MatrixXd;
 //   using stan::math::var;
 //   using stan::math::value_of;
-//   stan::math::matrix_exp_action_handler handler;
 //   std::srand(2999);
 
 //   for (int i = 1; i < 50; ++i) {
@@ -347,7 +227,6 @@ TEST(MathMatrix, matrix_exp_action_vv) {
 //   using Eigen::MatrixXd;
 //   using stan::math::var;
 //   using stan::math::value_of;
-//   stan::math::matrix_exp_action_handler handler;
 //   std::srand(2999);
 
 //   for (int i = 1; i < 50; ++i) {

test/unit/math/rev/mat/fun/scale_matrix_exp_multiply_test.cpp

+#include <stan/math/rev/mat.hpp>
+#include <gtest/gtest.h>
+// #include <stan/math/prim/scal/err/check_not_nan.hpp>
+#include <test/unit/math/rev/mat/fun/util.hpp>
+#include <test/unit/math/rev/mat/util.hpp>
+#include <stan/math/rev/mat/fun/scale_matrix_exp_multiply.hpp>
+#include <stan/math/rev/mat/fun/to_var.hpp>
+#include <vector>
+
+inline void test_scale_matrix_exp_multiply_dv(int N, int M) {
+  using stan::math::value_of;
+  using stan::math::var;
+
+  std::srand(1999);
+
+  Eigen::Matrix<var, -1, -1> Av = Eigen::Matrix<var, -1, -1>::Random(N, N);
+  Eigen::Matrix<var, -1, -1> Bv = Eigen::Matrix<var, -1, -1>::Random(N, M);
+  std::vector<stan::math::var> Avec = stan::math::to_array_1d(Av);
+  std::vector<stan::math::var> Bvec = stan::math::to_array_1d(Bv);
+  Eigen::MatrixXd A = value_of(Av);
+
+  // brute force
+  Eigen::Matrix<var, -1, -1> expAB
+      = stan::math::multiply(stan::math::matrix_exp(A), Bv);
+
+  // matrix_exp_multiply
+  const double t = 1.0;
+  Eigen::Matrix<var, -1, -1> res_dv
+      = stan::math::scale_matrix_exp_multiply(t, A, Bv);
+  EXPECT_EQ(res_dv.size(), expAB.size());
+  for (int l = 0; l < res_dv.size(); ++l) {
+    EXPECT_FLOAT_EQ(res_dv(l).val(), expAB(l).val());
+  }
+
+  // compare adjoints
+  std::vector<double> g, g0;
+  for (int l = 0; l < M; ++l) {
+    for (int k = 0; k < N; ++k) {
+      stan::math::set_zero_all_adjoints();
+      res_dv(k, l).grad(Bvec, g);
+      stan::math::set_zero_all_adjoints();
+      expAB(k, l).grad(Bvec, g0);
+      for (size_t j = 0; j < g.size(); ++j) {
+        EXPECT_FLOAT_EQ(g[j], g0[j]);
+      }
+    }
+  }
+
+  // test a single function of expA*B
+  var f = sum(res_dv);
+  var f0 = sum(expAB);
+  stan::math::set_zero_all_adjoints();
+  f.grad(Bvec, g);
+  stan::math::set_zero_all_adjoints();
+  f0.grad(Bvec, g0);
+  for (size_t j = 0; j < g.size(); ++j) {
+    EXPECT_FLOAT_EQ(g[j], g0[j]);
+  }
+}
+
+TEST(MathMatrix, scale_matrix_exp_multiply_dv) {
+  test_scale_matrix_exp_multiply_dv(1, 1);
+  test_scale_matrix_exp_multiply_dv(1, 5);
+  test_scale_matrix_exp_multiply_dv(5, 1);
+  test_scale_matrix_exp_multiply_dv(5, 5);
+}
+
+inline void test_scale_matrix_exp_multiply_vd(int N, int M) {
+  using stan::math::value_of;
+  using stan::math::var;
+
+  std::srand(1999);
+
+  Eigen::Matrix<var, -1, -1> Av = Eigen::Matrix<var, -1, -1>::Random(N, N);
+  Eigen::Matrix<var, -1, -1> Bv = Eigen::Matrix<var, -1, -1>::Random(N, M);
+  std::vector<stan::math::var> Avec = stan::math::to_array_1d(Av);
+  std::vector<stan::math::var> Bvec = stan::math::to_array_1d(Bv);
+  Eigen::MatrixXd B = value_of(Bv);
+
+  // brute force
+  Eigen::Matrix<var, -1, -1> expAB
+      = stan::math::multiply(stan::math::matrix_exp(Av), B);
+
+  // matrix_exp_multiply
+  const double t = 1.0;
+  Eigen::Matrix<var, -1, -1> res_vd
+      = stan::math::scale_matrix_exp_multiply(t, Av, B);
+  EXPECT_EQ(res_vd.size(), expAB.size());
+  for (int l = 0; l < res_vd.size(); ++l) {
+    EXPECT_FLOAT_EQ(res_vd(l).val(), expAB(l).val());
+  }
+
+  // compare adjoints
+  std::vector<double> g, g0;
+  for (int l = 0; l < M; ++l) {
+    for (int k = 0; k < N; ++k) {
+      stan::math::set_zero_all_adjoints();
+      res_vd(k, l).grad(Avec, g);
+      stan::math::set_zero_all_adjoints();
+      expAB(k, l).grad(Avec, g0);
+      for (size_t j = 0; j < g.size(); ++j) {
+        EXPECT_FLOAT_EQ(g[j], g0[j]);
+      }
+    }
+  }
+
+  // test a single function of expA*B
+  var f = sum(res_vd);
+  var f0 = sum(expAB);
+  stan::math::set_zero_all_adjoints();
+  f.grad(Avec, g);
+  stan::math::set_zero_all_adjoints();
+  f0.grad(Avec, g0);
+  for (size_t j = 0; j < g.size(); ++j) {
+    EXPECT_FLOAT_EQ(g[j], g0[j]);
+  }
+}
+
+TEST(MathMatrix, scale_matrix_exp_multiply_vd) {
+  test_scale_matrix_exp_multiply_vd(1, 1);
+  test_scale_matrix_exp_multiply_vd(1, 5);
+  test_scale_matrix_exp_multiply_vd(5, 1);
+  test_scale_matrix_exp_multiply_vd(5, 5);
+}
+
+inline void test_scale_matrix_exp_multiply_vv(int N, int M) {
+  using stan::math::value_of;
+  using stan::math::var;
+
+  std::srand(2999);
+
+  Eigen::Matrix<var, -1, -1> Av = Eigen::Matrix<var, -1, -1>::Random(N, N);
+  Eigen::Matrix<var, -1, -1> Bv = Eigen::Matrix<var, -1, -1>::Random(N, M);
+  std::vector<stan::math::var> Avec = stan::math::to_array_1d(Av);
+  std::vector<stan::math::var> Bvec = stan::math::to_array_1d(Bv);
+
+  // brute force
+  Eigen::Matrix<var, -1, -1> expAB
+      = stan::math::multiply(stan::math::matrix_exp(Av), Bv);
+
+  // matrix_exp_multiply
+  const double t = 1.0;
+  Eigen::Matrix<var, -1, -1> res_vv
+      = stan::math::scale_matrix_exp_multiply(t, Av, Bv);
+  EXPECT_EQ(res_vv.size(), expAB.size());
+  for (int l = 0; l < res_vv.size(); ++l) {
+    EXPECT_FLOAT_EQ(res_vv(l).val(), expAB(l).val());
+  }
+
+  // compare adjoints
+  Avec.insert(Avec.end(), Bvec.begin(), Bvec.end());
+  std::vector<double> g, g0;
+  for (int l = 0; l < M; ++l) {
+    for (int k = 0; k < N; ++k) {
+      stan::math::set_zero_all_adjoints();
+      res_vv(k, l).grad(Avec, g);
+      stan::math::set_zero_all_adjoints();
+      expAB(k, l).grad(Avec, g0);
+      for (size_t j = 0; j < g.size(); ++j) {
+        EXPECT_FLOAT_EQ(g[j], g0[j]);
+      }
+    }
+  }
+
+  // test a single function of expA*B
+  var f = sum(res_vv);
+  var f0 = sum(expAB);
+  stan::math::set_zero_all_adjoints();
+  f.grad(Avec, g);
+  stan::math::set_zero_all_adjoints();
+  f0.grad(Avec, g0);
+  for (size_t j = 0; j < g.size(); ++j) {
+    EXPECT_FLOAT_EQ(g[j], g0[j]);
+  }
+}
+
+TEST(MathMatrix, scale_matrix_exp_multiply_vv) {
+  test_scale_matrix_exp_multiply_vv(1, 1);
+  test_scale_matrix_exp_multiply_vv(1, 5);
+  test_scale_matrix_exp_multiply_vv(5, 1);
+  test_scale_matrix_exp_multiply_vv(5, 5);
+  test_scale_matrix_exp_multiply_vv(8, 2);
+}
+
+TEST(MathMatrix, scale_matrix_exp_multiply_exception) {
+  using stan::math::var;
+  const double t = 1.0;
+  {  // nonzero size
+    Eigen::Matrix<var, -1, -1> A(0, 0);
+    Eigen::Matrix<var, -1, -1> B = Eigen::Matrix<var, -1, -1>::Random(1, 2);
+    EXPECT_THROW(scale_matrix_exp_multiply(t, A, B), std::invalid_argument);
+    EXPECT_THROW(scale_matrix_exp_multiply(t, B, A), std::invalid_argument);
+  }
+
+  {  // multiplicable
+    Eigen::Matrix<var, -1, -1> A = Eigen::Matrix<var, -1, -1>::Random(2, 2);
+    Eigen::Matrix<var, -1, -1> B = Eigen::Matrix<var, -1, -1>::Random(3, 2);
+    EXPECT_THROW(scale_matrix_exp_multiply(t, A, B), std::invalid_argument);
+  }
+
+  {  // square
+    Eigen::Matrix<var, -1, -1> A = Eigen::Matrix<var, -1, -1>::Random(2, 3);
+    Eigen::Matrix<var, -1, -1> B = Eigen::Matrix<var, -1, -1>::Random(3, 2);
+    EXPECT_THROW(scale_matrix_exp_multiply(t, A, B), std::invalid_argument);
+  }
+}