From 73bf8b3c1d692ae916f54f93f39f7fd205a76979 Mon Sep 17 00:00:00 2001 From: Ben Goodrich Date: Thu, 7 Mar 2019 10:43:45 -0500 Subject: [PATCH 1/9] initial implementation of sqrt_spd and inv_sqrt_spd adj_jac_apply() and jacobian() seem to not work together --- stan/math/prim/mat.hpp | 2 + stan/math/prim/mat/fun/inv_sqrt_spd.hpp | 39 ++++++++++++ stan/math/prim/mat/fun/sqrt_spd.hpp | 39 ++++++++++++ stan/math/rev/mat.hpp | 17 ++--- stan/math/rev/mat/fun/sqrt_spd.hpp | 62 +++++++++++++++++++ .../math/prim/mat/fun/inv_sqrt_spd_test.cpp | 15 +++++ test/unit/math/prim/mat/fun/sqrt_spd_test.cpp | 15 +++++ .../math/rev/mat/fun/inv_sqrt_spd_test.cpp | 60 ++++++++++++++++++ test/unit/math/rev/mat/fun/sqrt_spd_test.cpp | 60 ++++++++++++++++++ 9 files changed, 301 insertions(+), 8 deletions(-) create mode 100644 stan/math/prim/mat/fun/inv_sqrt_spd.hpp create mode 100644 stan/math/prim/mat/fun/sqrt_spd.hpp create mode 100644 stan/math/rev/mat/fun/sqrt_spd.hpp create mode 100644 test/unit/math/prim/mat/fun/inv_sqrt_spd_test.cpp create mode 100644 test/unit/math/prim/mat/fun/sqrt_spd_test.cpp create mode 100644 test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp create mode 100644 test/unit/math/rev/mat/fun/sqrt_spd_test.cpp diff --git a/stan/math/prim/mat.hpp b/stan/math/prim/mat.hpp index 0dbdc17d0c8..a5575a86bb5 100644 --- a/stan/math/prim/mat.hpp +++ b/stan/math/prim/mat.hpp @@ -135,6 +135,7 @@ #include #include #include +#include #include #include #include @@ -221,6 +222,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/mat/fun/inv_sqrt_spd.hpp b/stan/math/prim/mat/fun/inv_sqrt_spd.hpp new file mode 100644 index 00000000000..1cc24575f86 --- /dev/null +++ b/stan/math/prim/mat/fun/inv_sqrt_spd.hpp @@ -0,0 +1,39 @@ +#ifndef STAN_MATH_PRIM_MAT_FUN_INV_SQRT_SPD_HPP +#define STAN_MATH_PRIM_MAT_FUN_INV_SQRT_SPD_HPP + +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Returns the inverse symmetric square root of the specified + * symmetric and positive definite matrix, which is defined as + * V * diag_matrix(inv_sqrt(d)) * V' + * where V is a matrix of eigenvectors of the input matrix and + * d is a vector of eigenvalue of the input matrix. If the + * input matrix is not actually positive definite, then the + * output matrix will have NaNs in it. + * + * @param[in] m Specified symmetric positive definite matrix. + * @return Inverse symmetric square root of the matrix. + * @throws If matrix is of size zero or asymmetric + */ + +template +Eigen::Matrix +inv_sqrt_spd(const Eigen::Matrix &m) { + check_nonzero_size("inv_sqrt_spd", "m", m); + check_symmetric("inv_sqrt_spd", "m", m); + + Eigen::SelfAdjointEigenSolver< + Eigen::Matrix> + solver(m); + return solver.operatorInverseSqrt(); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/prim/mat/fun/sqrt_spd.hpp b/stan/math/prim/mat/fun/sqrt_spd.hpp new file mode 100644 index 00000000000..e5cb88b8e75 --- /dev/null +++ b/stan/math/prim/mat/fun/sqrt_spd.hpp @@ -0,0 +1,39 @@ +#ifndef STAN_MATH_PRIM_MAT_FUN_SQRT_SPD_HPP +#define STAN_MATH_PRIM_MAT_FUN_SQRT_SPD_HPP + +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Returns the symmetric square root of the specified symmetric + * and positive definite matrix, which is defined as + * V * diag_matrix(sqrt(d)) * V' + * where V is a matrix of eigenvectors of the input matrix and + * d is a vector of eigenvalue of the input matrix. If the + * input matrix is not actually positive definite, then the + * output matrix will have NaNs in it. + * + * @param[in] m Specified symmetric positive definite matrix. + * @return Symmetric square root of the matrix. + * @throws If matrix is of size zero or asymmetric + */ + +template +Eigen::Matrix +sqrt_spd(const Eigen::Matrix &m) { + check_nonzero_size("sqrt_spd", "m", m); + check_symmetric("sqrt_spd", "m", m); + + Eigen::SelfAdjointEigenSolver< + Eigen::Matrix> + solver(m); + return solver.operatorSqrt(); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/rev/mat.hpp b/stan/math/rev/mat.hpp index 9138a7b9827..c3310f75e94 100644 --- a/stan/math/rev/mat.hpp +++ b/stan/math/rev/mat.hpp @@ -2,16 +2,18 @@ #define STAN_MATH_REV_MAT_HPP #include +#include #include #include -#include #include -#include #include #include +#include +#include +#include #include #include #include @@ -24,8 +26,6 @@ #include #include #include -#include -#include #include #include #include @@ -47,6 +47,7 @@ #include #include #include +#include #include #include #include @@ -61,13 +62,13 @@ #include #include -#include -#include -#include #include +#include +#include +#include #include #include -#include +#include #include #endif diff --git a/stan/math/rev/mat/fun/sqrt_spd.hpp b/stan/math/rev/mat/fun/sqrt_spd.hpp new file mode 100644 index 00000000000..4c0c29e241a --- /dev/null +++ b/stan/math/rev/mat/fun/sqrt_spd.hpp @@ -0,0 +1,62 @@ +#ifndef STAN_MATH_REV_MAT_FUN_SQRT_SPD_HPP +#define STAN_MATH_REV_MAT_FUN_SQRT_SPD_HPP + +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +namespace { +class sqrt_spd_op { + int K_; + double *y_; // Holds the output, which has K_^2 elements + +public: + sqrt_spd_op() : K_(0), y_(nullptr) {} + + template + Eigen::MatrixXd operator()(const std::array &needs_adj, + const Eigen::MatrixXd &m) { + K_ = m.rows(); + Eigen::MatrixXd sqrt_m = sqrt_spd(m); + int KK = K_ * K_; + y_ = ChainableStack::instance().memalloc_.alloc_array(KK); + for (int k = 0; k < KK; ++k) + y_[k] = sqrt_m.coeff(k); + return sqrt_m; + } + + // https://math.stackexchange.com/questions/540361/ + // derivative-or-differential-of-symmetric-square-root-of-a-matrix + + template + std::tuple + multiply_adjoint_jacobian(const std::array &needs_adj, + const Eigen::VectorXd &adj) const { + using Eigen::kroneckerProduct; + using Eigen::MatrixXd; + Eigen::Map sqrt_m(y_, K_, K_); + return std::make_tuple( + (kroneckerProduct(sqrt_m, MatrixXd::Identity(K_, K_)) + + kroneckerProduct(MatrixXd::Identity(K_, K_), sqrt_m)) + // the above is symmetric so can skip the transpose + .ldlt() + .solve(adj) + .eval()); + } +}; +} // namespace + +inline Eigen::Matrix +sqrt_spd(const Eigen::Matrix &m) { + return adj_jac_apply(m); +} + +} // namespace math +} // namespace stan +#endif diff --git a/test/unit/math/prim/mat/fun/inv_sqrt_spd_test.cpp b/test/unit/math/prim/mat/fun/inv_sqrt_spd_test.cpp new file mode 100644 index 00000000000..c4828cd8df6 --- /dev/null +++ b/test/unit/math/prim/mat/fun/inv_sqrt_spd_test.cpp @@ -0,0 +1,15 @@ +#include +#include + +TEST(MathMatrix, inv_sqrt_spd) { + stan::math::matrix_d m0; + stan::math::matrix_d m1(2, 3); + m1 << 1, 2, 3, 4, 5, 6; + stan::math::matrix_d ev_m1(1, 1); + ev_m1 << 4.0; + + using stan::math::inv_sqrt_spd; + EXPECT_THROW(inv_sqrt_spd(m0), std::invalid_argument); + EXPECT_NEAR(0.5, inv_sqrt_spd(ev_m1)(0, 0), 1e-16); + EXPECT_THROW(inv_sqrt_spd(m1), std::invalid_argument); +} diff --git a/test/unit/math/prim/mat/fun/sqrt_spd_test.cpp b/test/unit/math/prim/mat/fun/sqrt_spd_test.cpp new file mode 100644 index 00000000000..013bf96ca81 --- /dev/null +++ b/test/unit/math/prim/mat/fun/sqrt_spd_test.cpp @@ -0,0 +1,15 @@ +#include +#include + +TEST(MathMatrix, sqrt_spd) { + stan::math::matrix_d m0; + stan::math::matrix_d m1(2, 3); + m1 << 1, 2, 3, 4, 5, 6; + stan::math::matrix_d ev_m1(1, 1); + ev_m1 << 4.0; + + using stan::math::sqrt_spd; + EXPECT_THROW(sqrt_spd(m0), std::invalid_argument); + EXPECT_NEAR(2.0, sqrt_spd(ev_m1)(0, 0), 1e-16); + EXPECT_THROW(sqrt_spd(m1), std::invalid_argument); +} diff --git a/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp b/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp new file mode 100644 index 00000000000..a5c70cd424e --- /dev/null +++ b/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp @@ -0,0 +1,60 @@ +#include +#include +#include +#include + +TEST(AgradRevMatrix, check_varis_on_stack) { + using stan::math::inv_sqrt_spd; + using stan::math::matrix_v; + + matrix_v a(3, 3); + a << 1.0, 2.0, 3.0, 2.0, 5.0, 7.9, 3.0, 7.9, 1.08; + test::check_varis_on_stack(inv_sqrt_spd(a)); +} + +struct make_zero { + template + Eigen::Matrix + operator()(const Eigen::Matrix &a) const { + typedef Eigen::Matrix matrix_t; + const int K = static_cast(sqrt(a.rows())); + matrix_t A(K, K); + int pos = 0; + for (int i = 0; i < K; i++) + for (int j = 0; j < K; j++) + A(i, j) = a[pos++]; + + matrix_t inv_sqrt_A = stan::math::inv_sqrt_spd(A); + matrix_t zero = A.inverse() - inv_sqrt_A * inv_sqrt_A; + zero.resize(K * K, 1); + return zero; + } +}; + +TEST(AgradRevMatrix, sqrt_spd) { + using stan::math::sqrt_spd; + using stan::math::matrix_v; + using stan::math::vector_v; + using stan::math::matrix_d; + + int K = 2; + matrix_d A(K, K); + A.setRandom(); + A = A.transpose() * A; + matrix_d J; + Eigen::VectorXd f_x; + A.resize(K * K, 1); + make_zero f; + stan::math::jacobian(f, A, f_x, J); + const double TOL = 1e-14; + int pos = 0; + for (int i = 0; i < K; i++) + for (int j = 0; j < K; j++) + EXPECT_NEAR(f_x(pos++), 0, TOL); + /* + for (int i = 0; i < J.rows(); i++) + for (int j = 0; j < J.cols(); j++) + if (i != j) EXPECT_NEAR(J(i, j), 0, TOL); + */ + EXPECT_TRUE(J.array().any()); +} diff --git a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp new file mode 100644 index 00000000000..51e6f4c3166 --- /dev/null +++ b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp @@ -0,0 +1,60 @@ +#include +#include +#include +#include + +TEST(AgradRevMatrix, check_varis_on_stack) { + using stan::math::sqrt_spd; + using stan::math::matrix_v; + + matrix_v a(3, 3); + a << 1.0, 2.0, 3.0, 2.0, 5.0, 7.9, 3.0, 7.9, 1.08; + test::check_varis_on_stack(sqrt_spd(a)); +} + +/* +struct make_zero { + template + Eigen::Matrix operator()( + const Eigen::Matrix& a) const { + typedef Eigen::Matrix matrix_t; + const int K = static_cast(sqrt(a.rows())); + matrix_t A(K, K); + int pos = 0; + for (int i = 0; i < K; i++) + for (int j = 0; j < K; j++) + A(i, j) = a[pos++]; + + matrix_t sqrt_A = stan::math::sqrt_spd(A); + matrix_t zero = A - sqrt_A * sqrt_A; + zero.resize(K * K, 1); + return zero; + } +}; + +TEST(AgradRevMatrix, sqrt_spd) { + using stan::math::sqrt_spd; + using stan::math::matrix_v; + using stan::math::vector_v; + using stan::math::matrix_d; + + int K = 2; + matrix_d A(K, K); + A.setRandom(); + A = A.transpose() * A; + matrix_d J; + Eigen::VectorXd f_x; + A.resize(K * K, 1); + make_zero f; + stan::math::jacobian(f, A, f_x, J); + const double TOL = 1e-14; + int pos = 0; + for (int i = 0; i < K; i++) + for (int j = 0; j < K; j++) + EXPECT_NEAR(f_x(pos++), 0, TOL); + for (int i = 0; i < J.rows(); i++) + for (int j = 0; j < J.cols(); j++) + if (i != j) EXPECT_NEAR(J(i, j), 0, TOL); + EXPECT_TRUE(J.array().any()); +} +*/ From f85bae9deb71b82402eaed01dbe8aa9698bb5598 Mon Sep 17 00:00:00 2001 From: Ben Bales Date: Thu, 7 Mar 2019 16:46:19 -0800 Subject: [PATCH 2/9] Updated the input and output types of the sqrt_spd adj_jac_apply (Issue #1144) Rearranged some includes in stan/math/rev/mat.hpp to get rid of some type traits errors. Not sure why these were moved around. Might need moved back. --- stan/math/rev/mat.hpp | 6 +++--- stan/math/rev/mat/fun/sqrt_spd.hpp | 22 ++++++++++++-------- test/unit/math/rev/mat/fun/sqrt_spd_test.cpp | 4 ++-- 3 files changed, 18 insertions(+), 14 deletions(-) diff --git a/stan/math/rev/mat.hpp b/stan/math/rev/mat.hpp index c3310f75e94..79243c33f71 100644 --- a/stan/math/rev/mat.hpp +++ b/stan/math/rev/mat.hpp @@ -8,12 +8,10 @@ #include +#include #include #include -#include -#include -#include #include #include #include @@ -26,6 +24,8 @@ #include #include #include +#include +#include #include #include #include diff --git a/stan/math/rev/mat/fun/sqrt_spd.hpp b/stan/math/rev/mat/fun/sqrt_spd.hpp index 4c0c29e241a..fb6c5326f83 100644 --- a/stan/math/rev/mat/fun/sqrt_spd.hpp +++ b/stan/math/rev/mat/fun/sqrt_spd.hpp @@ -35,19 +35,23 @@ class sqrt_spd_op { // derivative-or-differential-of-symmetric-square-root-of-a-matrix template - std::tuple + std::tuple multiply_adjoint_jacobian(const std::array &needs_adj, - const Eigen::VectorXd &adj) const { + const Eigen::MatrixXd &adj) const { using Eigen::kroneckerProduct; using Eigen::MatrixXd; + Eigen::MatrixXd output(K_, K_); + Eigen::Map map_input(adj.data(), adj.size()); + Eigen::Map map_output(output.data(), output.size()); Eigen::Map sqrt_m(y_, K_, K_); - return std::make_tuple( - (kroneckerProduct(sqrt_m, MatrixXd::Identity(K_, K_)) + - kroneckerProduct(MatrixXd::Identity(K_, K_), sqrt_m)) - // the above is symmetric so can skip the transpose - .ldlt() - .solve(adj) - .eval()); + + map_output = (kroneckerProduct(sqrt_m, MatrixXd::Identity(K_, K_)) + + kroneckerProduct(MatrixXd::Identity(K_, K_), sqrt_m)) + // the above is symmetric so can skip the transpose + .ldlt() + .solve(map_input); + + return std::make_tuple(output); } }; } // namespace diff --git a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp index 51e6f4c3166..8f91292cf70 100644 --- a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp +++ b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp @@ -12,7 +12,7 @@ TEST(AgradRevMatrix, check_varis_on_stack) { test::check_varis_on_stack(sqrt_spd(a)); } -/* + struct make_zero { template Eigen::Matrix operator()( @@ -57,4 +57,4 @@ TEST(AgradRevMatrix, sqrt_spd) { if (i != j) EXPECT_NEAR(J(i, j), 0, TOL); EXPECT_TRUE(J.array().any()); } -*/ + From 5480a40a3738e86781fd8dfd9308b033e37ebf40 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 7 Mar 2019 19:53:48 -0500 Subject: [PATCH 3/9] [Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/stable/2017-11-14) --- stan/math/prim/mat/fun/inv_sqrt_spd.hpp | 8 ++--- stan/math/prim/mat/fun/sqrt_spd.hpp | 8 ++--- stan/math/rev/mat/fun/sqrt_spd.hpp | 32 +++++++++---------- .../math/rev/mat/fun/inv_sqrt_spd_test.cpp | 8 ++--- test/unit/math/rev/mat/fun/sqrt_spd_test.cpp | 11 +++---- 5 files changed, 33 insertions(+), 34 deletions(-) diff --git a/stan/math/prim/mat/fun/inv_sqrt_spd.hpp b/stan/math/prim/mat/fun/inv_sqrt_spd.hpp index 1cc24575f86..d1c162a1d2c 100644 --- a/stan/math/prim/mat/fun/inv_sqrt_spd.hpp +++ b/stan/math/prim/mat/fun/inv_sqrt_spd.hpp @@ -23,8 +23,8 @@ namespace math { */ template -Eigen::Matrix -inv_sqrt_spd(const Eigen::Matrix &m) { +Eigen::Matrix inv_sqrt_spd( + const Eigen::Matrix &m) { check_nonzero_size("inv_sqrt_spd", "m", m); check_symmetric("inv_sqrt_spd", "m", m); @@ -34,6 +34,6 @@ inv_sqrt_spd(const Eigen::Matrix &m) { return solver.operatorInverseSqrt(); } -} // namespace math -} // namespace stan +} // namespace math +} // namespace stan #endif diff --git a/stan/math/prim/mat/fun/sqrt_spd.hpp b/stan/math/prim/mat/fun/sqrt_spd.hpp index e5cb88b8e75..8069cf824db 100644 --- a/stan/math/prim/mat/fun/sqrt_spd.hpp +++ b/stan/math/prim/mat/fun/sqrt_spd.hpp @@ -23,8 +23,8 @@ namespace math { */ template -Eigen::Matrix -sqrt_spd(const Eigen::Matrix &m) { +Eigen::Matrix sqrt_spd( + const Eigen::Matrix &m) { check_nonzero_size("sqrt_spd", "m", m); check_symmetric("sqrt_spd", "m", m); @@ -34,6 +34,6 @@ sqrt_spd(const Eigen::Matrix &m) { return solver.operatorSqrt(); } -} // namespace math -} // namespace stan +} // namespace math +} // namespace stan #endif diff --git a/stan/math/rev/mat/fun/sqrt_spd.hpp b/stan/math/rev/mat/fun/sqrt_spd.hpp index fb6c5326f83..cf8d65c8d90 100644 --- a/stan/math/rev/mat/fun/sqrt_spd.hpp +++ b/stan/math/rev/mat/fun/sqrt_spd.hpp @@ -14,9 +14,9 @@ namespace math { namespace { class sqrt_spd_op { int K_; - double *y_; // Holds the output, which has K_^2 elements + double *y_; // Holds the output, which has K_^2 elements -public: + public: sqrt_spd_op() : K_(0), y_(nullptr) {} template @@ -35,32 +35,32 @@ class sqrt_spd_op { // derivative-or-differential-of-symmetric-square-root-of-a-matrix template - std::tuple - multiply_adjoint_jacobian(const std::array &needs_adj, - const Eigen::MatrixXd &adj) const { - using Eigen::kroneckerProduct; + std::tuple multiply_adjoint_jacobian( + const std::array &needs_adj, + const Eigen::MatrixXd &adj) const { using Eigen::MatrixXd; + using Eigen::kroneckerProduct; Eigen::MatrixXd output(K_, K_); Eigen::Map map_input(adj.data(), adj.size()); Eigen::Map map_output(output.data(), output.size()); Eigen::Map sqrt_m(y_, K_, K_); - map_output = (kroneckerProduct(sqrt_m, MatrixXd::Identity(K_, K_)) + - kroneckerProduct(MatrixXd::Identity(K_, K_), sqrt_m)) - // the above is symmetric so can skip the transpose - .ldlt() - .solve(map_input); + map_output = (kroneckerProduct(sqrt_m, MatrixXd::Identity(K_, K_)) + + kroneckerProduct(MatrixXd::Identity(K_, K_), sqrt_m)) + // the above is symmetric so can skip the transpose + .ldlt() + .solve(map_input); return std::make_tuple(output); } }; -} // namespace +} // namespace -inline Eigen::Matrix -sqrt_spd(const Eigen::Matrix &m) { +inline Eigen::Matrix sqrt_spd( + const Eigen::Matrix &m) { return adj_jac_apply(m); } -} // namespace math -} // namespace stan +} // namespace math +} // namespace stan #endif diff --git a/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp b/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp index a5c70cd424e..998b0d8f1a3 100644 --- a/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp +++ b/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp @@ -14,8 +14,8 @@ TEST(AgradRevMatrix, check_varis_on_stack) { struct make_zero { template - Eigen::Matrix - operator()(const Eigen::Matrix &a) const { + Eigen::Matrix operator()( + const Eigen::Matrix &a) const { typedef Eigen::Matrix matrix_t; const int K = static_cast(sqrt(a.rows())); matrix_t A(K, K); @@ -32,10 +32,10 @@ struct make_zero { }; TEST(AgradRevMatrix, sqrt_spd) { - using stan::math::sqrt_spd; + using stan::math::matrix_d; using stan::math::matrix_v; + using stan::math::sqrt_spd; using stan::math::vector_v; - using stan::math::matrix_d; int K = 2; matrix_d A(K, K); diff --git a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp index 8f91292cf70..34e933a33d8 100644 --- a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp +++ b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp @@ -4,15 +4,14 @@ #include TEST(AgradRevMatrix, check_varis_on_stack) { - using stan::math::sqrt_spd; using stan::math::matrix_v; + using stan::math::sqrt_spd; matrix_v a(3, 3); a << 1.0, 2.0, 3.0, 2.0, 5.0, 7.9, 3.0, 7.9, 1.08; test::check_varis_on_stack(sqrt_spd(a)); } - struct make_zero { template Eigen::Matrix operator()( @@ -33,10 +32,10 @@ struct make_zero { }; TEST(AgradRevMatrix, sqrt_spd) { - using stan::math::sqrt_spd; + using stan::math::matrix_d; using stan::math::matrix_v; + using stan::math::sqrt_spd; using stan::math::vector_v; - using stan::math::matrix_d; int K = 2; matrix_d A(K, K); @@ -54,7 +53,7 @@ TEST(AgradRevMatrix, sqrt_spd) { EXPECT_NEAR(f_x(pos++), 0, TOL); for (int i = 0; i < J.rows(); i++) for (int j = 0; j < J.cols(); j++) - if (i != j) EXPECT_NEAR(J(i, j), 0, TOL); + if (i != j) + EXPECT_NEAR(J(i, j), 0, TOL); EXPECT_TRUE(J.array().any()); } - From f89b6a548d64b621953e1a58505fa6f7a5966e6c Mon Sep 17 00:00:00 2001 From: Ben Goodrich Date: Thu, 7 Mar 2019 23:08:57 -0500 Subject: [PATCH 4/9] fix order of headers --- stan/math/rev/mat/fun/sqrt_spd.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/rev/mat/fun/sqrt_spd.hpp b/stan/math/rev/mat/fun/sqrt_spd.hpp index cf8d65c8d90..dfde960f149 100644 --- a/stan/math/rev/mat/fun/sqrt_spd.hpp +++ b/stan/math/rev/mat/fun/sqrt_spd.hpp @@ -1,11 +1,11 @@ #ifndef STAN_MATH_REV_MAT_FUN_SQRT_SPD_HPP #define STAN_MATH_REV_MAT_FUN_SQRT_SPD_HPP -#include #include #include -#include #include +#include +#include #include namespace stan { From ca80605705f5d02db8bf3942419a800bf9b91301 Mon Sep 17 00:00:00 2001 From: Ben Goodrich Date: Thu, 7 Mar 2019 23:09:40 -0500 Subject: [PATCH 5/9] bump size of test up [CI skip] --- test/unit/math/rev/mat/fun/sqrt_spd_test.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp index 34e933a33d8..7f2c25f99f7 100644 --- a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp +++ b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp @@ -14,8 +14,8 @@ TEST(AgradRevMatrix, check_varis_on_stack) { struct make_zero { template - Eigen::Matrix operator()( - const Eigen::Matrix& a) const { + Eigen::Matrix + operator()(const Eigen::Matrix &a) const { typedef Eigen::Matrix matrix_t; const int K = static_cast(sqrt(a.rows())); matrix_t A(K, K); @@ -37,7 +37,7 @@ TEST(AgradRevMatrix, sqrt_spd) { using stan::math::sqrt_spd; using stan::math::vector_v; - int K = 2; + int K = 9; matrix_d A(K, K); A.setRandom(); A = A.transpose() * A; @@ -53,7 +53,6 @@ TEST(AgradRevMatrix, sqrt_spd) { EXPECT_NEAR(f_x(pos++), 0, TOL); for (int i = 0; i < J.rows(); i++) for (int j = 0; j < J.cols(); j++) - if (i != j) - EXPECT_NEAR(J(i, j), 0, TOL); + EXPECT_NEAR(J(i, j), 0, TOL); EXPECT_TRUE(J.array().any()); } From 2fa0a27a84690665242861f25f0c19dec9870098 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 7 Mar 2019 23:10:55 -0500 Subject: [PATCH 6/9] [Jenkins] auto-formatting by clang-format version 6.0.0 (tags/google/stable/2017-11-14) --- test/unit/math/rev/mat/fun/sqrt_spd_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp index 7f2c25f99f7..3b2a46251a8 100644 --- a/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp +++ b/test/unit/math/rev/mat/fun/sqrt_spd_test.cpp @@ -14,8 +14,8 @@ TEST(AgradRevMatrix, check_varis_on_stack) { struct make_zero { template - Eigen::Matrix - operator()(const Eigen::Matrix &a) const { + Eigen::Matrix operator()( + const Eigen::Matrix &a) const { typedef Eigen::Matrix matrix_t; const int K = static_cast(sqrt(a.rows())); matrix_t A(K, K); From 2b1f2fe492102718e356f256a9fcd1b5ca15efe8 Mon Sep 17 00:00:00 2001 From: Ben Goodrich Date: Fri, 8 Mar 2019 01:00:12 -0500 Subject: [PATCH 7/9] add check for NaNs --- test/unit/math/prim/mat/fun/sqrt_spd_test.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/unit/math/prim/mat/fun/sqrt_spd_test.cpp b/test/unit/math/prim/mat/fun/sqrt_spd_test.cpp index 013bf96ca81..6dadc4b1e65 100644 --- a/test/unit/math/prim/mat/fun/sqrt_spd_test.cpp +++ b/test/unit/math/prim/mat/fun/sqrt_spd_test.cpp @@ -12,4 +12,6 @@ TEST(MathMatrix, sqrt_spd) { EXPECT_THROW(sqrt_spd(m0), std::invalid_argument); EXPECT_NEAR(2.0, sqrt_spd(ev_m1)(0, 0), 1e-16); EXPECT_THROW(sqrt_spd(m1), std::invalid_argument); + ev_m1(0, 0) = -1; + EXPECT_TRUE(std::isnan(sqrt_spd(ev_m1)(0, 0))); } From 1751a768db0f811ca8ba1fc67861156106ea18a3 Mon Sep 17 00:00:00 2001 From: Ben Goodrich Date: Sat, 9 Mar 2019 19:46:24 -0500 Subject: [PATCH 8/9] more usage of adj_jac_apply() --- stan/math/rev/mat.hpp | 3 + stan/math/rev/mat/fun/inv_sqrt_spd.hpp | 83 +++++++++++++++++++ stan/math/rev/mat/fun/inverse.hpp | 54 ++++++++++++ stan/math/rev/mat/fun/inverse_spd.hpp | 63 ++++++++++++++ stan/math/rev/mat/fun/sqrt_spd.hpp | 57 ++++++++----- .../math/rev/mat/fun/inv_sqrt_spd_test.cpp | 24 +++--- .../math/rev/mat/fun/inverse_spd_test.cpp | 72 +++++++++++++--- test/unit/math/rev/mat/fun/inverse_test.cpp | 45 ++++++++++ 8 files changed, 353 insertions(+), 48 deletions(-) create mode 100644 stan/math/rev/mat/fun/inv_sqrt_spd.hpp create mode 100644 stan/math/rev/mat/fun/inverse.hpp create mode 100644 stan/math/rev/mat/fun/inverse_spd.hpp diff --git a/stan/math/rev/mat.hpp b/stan/math/rev/mat.hpp index 79243c33f71..f83479fba60 100644 --- a/stan/math/rev/mat.hpp +++ b/stan/math/rev/mat.hpp @@ -24,6 +24,9 @@ #include #include #include +#include +#include +#include #include #include #include diff --git a/stan/math/rev/mat/fun/inv_sqrt_spd.hpp b/stan/math/rev/mat/fun/inv_sqrt_spd.hpp new file mode 100644 index 00000000000..94fcd3b5f74 --- /dev/null +++ b/stan/math/rev/mat/fun/inv_sqrt_spd.hpp @@ -0,0 +1,83 @@ +#ifndef STAN_MATH_REV_MAT_FUN_INV_SQRT_SPD_HPP +#define STAN_MATH_REV_MAT_FUN_INV_SQRT_SPD_HPP + +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +namespace { +class inv_sqrt_spd_op { + int K_; + double *J_; // Holds the Jacobian, which has K_^2 rows & columns + + public: + inv_sqrt_spd_op() : K_(0), J_(nullptr) {} + + template + Eigen::MatrixXd operator()(const std::array &needs_adj, + const Eigen::MatrixXd &m) { + using Eigen::kroneckerProduct; + using Eigen::MatrixXd; + using Eigen::VectorXd; + K_ = m.rows(); + Eigen::SelfAdjointEigenSolver solver(m); + MatrixXd V = solver.eigenvectors(); + VectorXd sqrt_ev = solver.eigenvalues().cwiseSqrt(); + MatrixXd inv_sqrt_m + = V * sqrt_ev.cwiseInverse().asDiagonal() * V.transpose(); + MatrixXd sqrt_m = V * sqrt_ev.asDiagonal() * V.transpose(); + int K2 = K_ * K_; + VectorXd d(K2); + // https://math.stackexchange.com/questions/540361/ + // and then simplifying the inverse of a Kronecker sum + // and multiplying by the Jacobian of the inverse of a SPD matrix + int pos = 0; + for (int i = 0; i < K_; i++) + for (int j = 0; j < K_; j++) + d.coeffRef(pos++) = sqrt_ev.coeff(i) * sqrt_ev.coeff(j) + * (sqrt_ev.coeff(i) + sqrt_ev.coeff(j)); + MatrixXd U = kroneckerProduct(V, V); + MatrixXd J = -U * d.cwiseInverse().asDiagonal() * U.transpose(); + int K4 = K2 * K2; + J_ = ChainableStack::instance().memalloc_.alloc_array(K4); + for (int k = 0; k < K4; ++k) + J_[k] = J.coeff(k); + return inv_sqrt_m; + } + + template + std::tuple multiply_adjoint_jacobian( + const std::array &needs_adj, + const Eigen::MatrixXd &adj) const { + using Eigen::Map; + using Eigen::MatrixXd; + using Eigen::VectorXd; + MatrixXd output(K_, K_); + Map map_output(output.data(), output.size()); + Map map_input(adj.data(), adj.size()); + Map J(J_, adj.size(), adj.size()); + // J is symmetric so do not need to transpose it + map_output = J.selfadjointView() * map_input; + return std::make_tuple(output); + } +}; +} // namespace + +inline Eigen::Matrix inv_sqrt_spd( + const Eigen::Matrix &m) { + check_nonzero_size("inv_sqrt_spd", "m", m); + check_symmetric("inv_sqrt_spd", "m", m); + return adj_jac_apply(m); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/rev/mat/fun/inverse.hpp b/stan/math/rev/mat/fun/inverse.hpp new file mode 100644 index 00000000000..986df3805ca --- /dev/null +++ b/stan/math/rev/mat/fun/inverse.hpp @@ -0,0 +1,54 @@ +#ifndef STAN_MATH_REV_MAT_FUN_INVERSE_HPP +#define STAN_MATH_REV_MAT_FUN_INVERSE_HPP + +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +namespace { +class inverse_op { + int K_; + double *y_; // Holds the inverse, which has K_ rows & columns + + public: + inverse_op() : K_(0), y_(nullptr) {} + + template + Eigen::MatrixXd operator()(const std::array &needs_adj, + const Eigen::MatrixXd &m) { + using Eigen::MatrixXd; + K_ = m.rows(); + MatrixXd inv_m = m.inverse(); + int K2 = K_ * K_; + y_ = ChainableStack::instance().memalloc_.alloc_array(K2); + for (int k = 0; k < K2; ++k) + y_[k] = inv_m.coeff(k); + return inv_m; + } + + template + std::tuple multiply_adjoint_jacobian( + const std::array &needs_adj, + const Eigen::MatrixXd &adj) const { + using Eigen::MatrixXd; + Eigen::Map y(y_, adj.rows(), adj.cols()); + MatrixXd output = -y.transpose() * adj * y.transpose(); + return std::make_tuple(output); + } +}; +} // namespace + +inline Eigen::Matrix inverse( + const Eigen::Matrix &m) { + check_square("inverse", "m", m); + return adj_jac_apply(m); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/rev/mat/fun/inverse_spd.hpp b/stan/math/rev/mat/fun/inverse_spd.hpp new file mode 100644 index 00000000000..f5d67c7ecfc --- /dev/null +++ b/stan/math/rev/mat/fun/inverse_spd.hpp @@ -0,0 +1,63 @@ +#ifndef STAN_MATH_REV_MAT_FUN_INVERSE_SPD_HPP +#define STAN_MATH_REV_MAT_FUN_INVERSE_SPD_HPP + +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +namespace { +class inverse_spd_op { + int K_; + double *y_; // Holds the inverse, which has K_ rows & columns + + public: + inverse_spd_op() : K_(0), y_(nullptr) {} + + template + Eigen::MatrixXd operator()(const std::array &needs_adj, + const Eigen::MatrixXd &m) { + using Eigen::MatrixXd; + K_ = m.rows(); + Eigen::LDLT ldlt(m); + if (ldlt.info() != Eigen::Success) + domain_error("invese_spd", "LDLT factor failed", "", ""); + if (!ldlt.isPositive() || ldlt.vectorD().coeff(K_ - 1) <= 0) + domain_error("invese_spd", "matrix not positive definite", "", ""); + MatrixXd inv_m = ldlt.solve(MatrixXd::Identity(m.rows(), m.cols())); + + int K2 = K_ * K_; + y_ = ChainableStack::instance().memalloc_.alloc_array(K2); + for (int k = 0; k < K2; ++k) + y_[k] = inv_m.coeff(k); + return inv_m; + } + + template + std::tuple multiply_adjoint_jacobian( + const std::array &needs_adj, + const Eigen::MatrixXd &adj) const { + using Eigen::MatrixXd; + Eigen::Map y(y_, K_, K_); + MatrixXd output = y.selfadjointView() * -adj + * y.selfadjointView(); + return std::make_tuple(output); + } +}; +} // namespace + +inline Eigen::Matrix inverse_spd( + const Eigen::Matrix &m) { + check_symmetric("inverse_spd", "m", m); + return adj_jac_apply(m); +} + +} // namespace math +} // namespace stan +#endif diff --git a/stan/math/rev/mat/fun/sqrt_spd.hpp b/stan/math/rev/mat/fun/sqrt_spd.hpp index dfde960f149..9d2abb93081 100644 --- a/stan/math/rev/mat/fun/sqrt_spd.hpp +++ b/stan/math/rev/mat/fun/sqrt_spd.hpp @@ -1,6 +1,8 @@ #ifndef STAN_MATH_REV_MAT_FUN_SQRT_SPD_HPP #define STAN_MATH_REV_MAT_FUN_SQRT_SPD_HPP +#include +#include #include #include #include @@ -14,43 +16,52 @@ namespace math { namespace { class sqrt_spd_op { int K_; - double *y_; // Holds the output, which has K_^2 elements + double *J_; // Holds the Jacobian, which has K_^2 rows & columns public: - sqrt_spd_op() : K_(0), y_(nullptr) {} + sqrt_spd_op() : K_(0), J_(nullptr) {} template Eigen::MatrixXd operator()(const std::array &needs_adj, const Eigen::MatrixXd &m) { + using Eigen::kroneckerProduct; + using Eigen::MatrixXd; + using Eigen::VectorXd; K_ = m.rows(); - Eigen::MatrixXd sqrt_m = sqrt_spd(m); - int KK = K_ * K_; - y_ = ChainableStack::instance().memalloc_.alloc_array(KK); - for (int k = 0; k < KK; ++k) - y_[k] = sqrt_m.coeff(k); + Eigen::SelfAdjointEigenSolver solver(m); + MatrixXd V = solver.eigenvectors(); + VectorXd sqrt_ev = solver.eigenvalues().cwiseSqrt(); + MatrixXd sqrt_m = V * sqrt_ev.asDiagonal() * V.transpose(); + int K2 = K_ * K_; + VectorXd d(K2); + // https://math.stackexchange.com/questions/540361/ + // and then simplifying the inverse of a Kronecker sum + int pos = 0; + for (int i = 0; i < K_; i++) + for (int j = 0; j < K_; j++) + d.coeffRef(pos++) = sqrt_ev.coeff(i) + sqrt_ev.coeff(j); + MatrixXd U = kroneckerProduct(V, V); + MatrixXd J = U * d.cwiseInverse().asDiagonal() * U.transpose(); + int K4 = K2 * K2; + J_ = ChainableStack::instance().memalloc_.alloc_array(K4); + for (int k = 0; k < K4; ++k) + J_[k] = J.coeff(k); return sqrt_m; } - // https://math.stackexchange.com/questions/540361/ - // derivative-or-differential-of-symmetric-square-root-of-a-matrix - template std::tuple multiply_adjoint_jacobian( const std::array &needs_adj, const Eigen::MatrixXd &adj) const { + using Eigen::Map; using Eigen::MatrixXd; - using Eigen::kroneckerProduct; - Eigen::MatrixXd output(K_, K_); - Eigen::Map map_input(adj.data(), adj.size()); - Eigen::Map map_output(output.data(), output.size()); - Eigen::Map sqrt_m(y_, K_, K_); - - map_output = (kroneckerProduct(sqrt_m, MatrixXd::Identity(K_, K_)) - + kroneckerProduct(MatrixXd::Identity(K_, K_), sqrt_m)) - // the above is symmetric so can skip the transpose - .ldlt() - .solve(map_input); - + using Eigen::VectorXd; + MatrixXd output(K_, K_); + Map map_output(output.data(), output.size()); + Map map_input(adj.data(), adj.size()); + Map J(J_, adj.size(), adj.size()); + // J is symmetric so do not need to transpose it + map_output = J.selfadjointView() * map_input; return std::make_tuple(output); } }; @@ -58,6 +69,8 @@ class sqrt_spd_op { inline Eigen::Matrix sqrt_spd( const Eigen::Matrix &m) { + check_nonzero_size("sqrt_spd", "m", m); + check_symmetric("sqrt_spd", "m", m); return adj_jac_apply(m); } diff --git a/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp b/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp index 998b0d8f1a3..740080c7505 100644 --- a/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp +++ b/test/unit/math/rev/mat/fun/inv_sqrt_spd_test.cpp @@ -12,7 +12,7 @@ TEST(AgradRevMatrix, check_varis_on_stack) { test::check_varis_on_stack(inv_sqrt_spd(a)); } -struct make_zero { +struct make_I { template Eigen::Matrix operator()( const Eigen::Matrix &a) const { @@ -25,9 +25,9 @@ struct make_zero { A(i, j) = a[pos++]; matrix_t inv_sqrt_A = stan::math::inv_sqrt_spd(A); - matrix_t zero = A.inverse() - inv_sqrt_A * inv_sqrt_A; - zero.resize(K * K, 1); - return zero; + matrix_t I = inv_sqrt_A * A * inv_sqrt_A; + I.resize(K * K, 1); + return I; } }; @@ -37,24 +37,22 @@ TEST(AgradRevMatrix, sqrt_spd) { using stan::math::sqrt_spd; using stan::math::vector_v; - int K = 2; + int K = 7; matrix_d A(K, K); A.setRandom(); A = A.transpose() * A; matrix_d J; Eigen::VectorXd f_x; A.resize(K * K, 1); - make_zero f; + make_I f; stan::math::jacobian(f, A, f_x, J); - const double TOL = 1e-14; + const double TOL = 1e-13; int pos = 0; for (int i = 0; i < K; i++) for (int j = 0; j < K; j++) - EXPECT_NEAR(f_x(pos++), 0, TOL); - /* - for (int i = 0; i < J.rows(); i++) - for (int j = 0; j < J.cols(); j++) - if (i != j) EXPECT_NEAR(J(i, j), 0, TOL); - */ + EXPECT_NEAR(f_x(pos++), i == j, TOL); + for (int i = 0; i < J.rows(); i++) + for (int j = 0; j < J.cols(); j++) + EXPECT_NEAR(J(i, j), 0, 10 * TOL); EXPECT_TRUE(J.array().any()); } diff --git a/test/unit/math/rev/mat/fun/inverse_spd_test.cpp b/test/unit/math/rev/mat/fun/inverse_spd_test.cpp index bd512df8988..21701d37182 100644 --- a/test/unit/math/rev/mat/fun/inverse_spd_test.cpp +++ b/test/unit/math/rev/mat/fun/inverse_spd_test.cpp @@ -3,6 +3,60 @@ #include #include +TEST(AgradRevMatrix, check_varis_on_stack) { + using stan::math::inverse_spd; + using stan::math::matrix_v; + + matrix_v a(2, 2); + a << 2.0, 3.0, 3.0, 7.0; + + test::check_varis_on_stack(stan::math::inverse_spd(a)); +} + +struct make_I { + template + Eigen::Matrix operator()( + const Eigen::Matrix &a) const { + typedef Eigen::Matrix matrix_t; + const int K = static_cast(sqrt(a.rows())); + matrix_t A(K, K); + int pos = 0; + for (int i = 0; i < K; i++) + for (int j = 0; j < K; j++) + A(i, j) = a[pos++]; + + matrix_t I = A * stan::math::inverse_spd(A); + I.resize(K * K, 1); + return I; + } +}; + +TEST(AgradRevMatrix, inverse) { + using stan::math::matrix_d; + using stan::math::matrix_v; + using stan::math::sqrt_spd; + using stan::math::vector_v; + + int K = 11; + matrix_d A(K, K); + A.setRandom(); + A = A.transpose() * A; + matrix_d J; + Eigen::VectorXd f_x; + A.resize(K * K, 1); + make_I f; + stan::math::jacobian(f, A, f_x, J); + const double TOL = 1e-13; + int pos = 0; + for (int i = 0; i < K; i++) + for (int j = 0; j < K; j++) + EXPECT_NEAR(f_x(pos++), i == j, TOL); + for (int i = 0; i < J.rows(); i++) + for (int j = 0; j < J.cols(); j++) + EXPECT_NEAR(J(i, j), 0, 10 * TOL); + EXPECT_TRUE(J.array().any()); +} + TEST(AgradRevMatrix, inverse_spd_val) { using stan::math::inverse_spd; using stan::math::matrix_v; @@ -26,7 +80,7 @@ TEST(AgradRevMatrix, inverse_spd_val) { a << 1.0, -1.0, -1.0, -1.0; EXPECT_THROW(inverse_spd(a), std::domain_error); } - +/* TEST(AgradRevMatrix, inverse_spd_grad) { using stan::math::inverse_spd; using stan::math::matrix_v; @@ -39,12 +93,13 @@ TEST(AgradRevMatrix, inverse_spd_grad) { AVEC x = createAVEC(ad(0, 0), ad(0, 1), ad(1, 0), ad(1, 1)); matrix_v ad_inv = inverse_spd(ad); + std::cout << "ad_inv = " << std::endl << ad_inv << std::endl; // int k = 0; // int l = 1; - VEC g; + VEC g(4); (0.5 * (ad_inv(k, l) + ad_inv(l, k))).grad(x, g); - + std::cout << "passed" << std::endl; int idx = 0; for (int i = 0; i < 2; ++i) { for (int j = 0; j < 2; ++j) { @@ -85,13 +140,4 @@ TEST(AgradRevMatrix, inverse_spd_inverse_spd_sum) { } } } - -TEST(AgradRevMatrix, check_varis_on_stack) { - using stan::math::inverse_spd; - using stan::math::matrix_v; - - matrix_v a(2, 2); - a << 2.0, 3.0, 3.0, 7.0; - - test::check_varis_on_stack(stan::math::inverse_spd(a)); -} +*/ diff --git a/test/unit/math/rev/mat/fun/inverse_test.cpp b/test/unit/math/rev/mat/fun/inverse_test.cpp index b5b8bfe1247..0349b6020b7 100644 --- a/test/unit/math/rev/mat/fun/inverse_test.cpp +++ b/test/unit/math/rev/mat/fun/inverse_test.cpp @@ -3,6 +3,50 @@ #include #include +struct make_I { + template + Eigen::Matrix operator()( + const Eigen::Matrix &a) const { + typedef Eigen::Matrix matrix_t; + const int K = static_cast(sqrt(a.rows())); + matrix_t A(K, K); + int pos = 0; + for (int i = 0; i < K; i++) + for (int j = 0; j < K; j++) + A(i, j) = a[pos++]; + + matrix_t I = A * stan::math::inverse(A); + I.resize(K * K, 1); + return I; + } +}; + +TEST(AgradRevMatrix, inverse) { + using stan::math::matrix_d; + using stan::math::matrix_v; + using stan::math::sqrt_spd; + using stan::math::vector_v; + + int K = 7; + matrix_d A(K, K); + A.setRandom(); + A = A.transpose() * A; + matrix_d J; + Eigen::VectorXd f_x; + A.resize(K * K, 1); + make_I f; + stan::math::jacobian(f, A, f_x, J); + const double TOL = 1e-13; + int pos = 0; + for (int i = 0; i < K; i++) + for (int j = 0; j < K; j++) + EXPECT_NEAR(f_x(pos++), i == j, TOL); + for (int i = 0; i < J.rows(); i++) + for (int j = 0; j < J.cols(); j++) + EXPECT_NEAR(J(i, j), 0, 100 * TOL); + EXPECT_TRUE(J.array().any()); +} + TEST(AgradRevMatrix, inverse_val) { using stan::math::inverse; using stan::math::matrix_v; @@ -21,6 +65,7 @@ TEST(AgradRevMatrix, inverse_val) { EXPECT_THROW(inverse(matrix_v(2, 3)), std::invalid_argument); } + TEST(AgradRevMatrix, inverse_grad) { using stan::math::inverse; using stan::math::matrix_v; From 1adbee1b20401a3908a1ba81ded824873d169cb3 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 10 Mar 2019 00:47:27 +0000 Subject: [PATCH 9/9] [Jenkins] auto-formatting by clang-format version 5.0.0-3~16.04.1 (tags/RELEASE_500/final) --- stan/math/rev/mat/fun/inv_sqrt_spd.hpp | 2 +- stan/math/rev/mat/fun/sqrt_spd.hpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/rev/mat/fun/inv_sqrt_spd.hpp b/stan/math/rev/mat/fun/inv_sqrt_spd.hpp index 94fcd3b5f74..c3db4cf72cd 100644 --- a/stan/math/rev/mat/fun/inv_sqrt_spd.hpp +++ b/stan/math/rev/mat/fun/inv_sqrt_spd.hpp @@ -24,9 +24,9 @@ class inv_sqrt_spd_op { template Eigen::MatrixXd operator()(const std::array &needs_adj, const Eigen::MatrixXd &m) { - using Eigen::kroneckerProduct; using Eigen::MatrixXd; using Eigen::VectorXd; + using Eigen::kroneckerProduct; K_ = m.rows(); Eigen::SelfAdjointEigenSolver solver(m); MatrixXd V = solver.eigenvectors(); diff --git a/stan/math/rev/mat/fun/sqrt_spd.hpp b/stan/math/rev/mat/fun/sqrt_spd.hpp index 9d2abb93081..3a1fabd037f 100644 --- a/stan/math/rev/mat/fun/sqrt_spd.hpp +++ b/stan/math/rev/mat/fun/sqrt_spd.hpp @@ -24,9 +24,9 @@ class sqrt_spd_op { template Eigen::MatrixXd operator()(const std::array &needs_adj, const Eigen::MatrixXd &m) { - using Eigen::kroneckerProduct; using Eigen::MatrixXd; using Eigen::VectorXd; + using Eigen::kroneckerProduct; K_ = m.rows(); Eigen::SelfAdjointEigenSolver solver(m); MatrixXd V = solver.eigenvectors();