From b5ed86ea3858752cd383dac6092a250494d59a68 Mon Sep 17 00:00:00 2001 From: spinkney Date: Tue, 1 Dec 2020 13:49:45 -0500 Subject: [PATCH 001/102] Added generalized_inverse files Added generalized_inverse files --- stan/math/prim/fun/generalized_inverse.hpp | 72 +++++++++++++++++++ .../math/prim/fun/generalied_inverse_test.cpp | 17 +++++ 2 files changed, 89 insertions(+) create mode 100644 stan/math/prim/fun/generalized_inverse.hpp create mode 100644 test/unit/math/prim/fun/generalied_inverse_test.cpp diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp new file mode 100644 index 00000000000..ff2c26c9ebf --- /dev/null +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -0,0 +1,72 @@ +#ifndef STAN_MATH_PRIM_FUN_GENERALIZED_INVERSE_HPP +#define STAN_MATH_PRIM_FUN_GENERALIZED_INVERSE_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Returns the Moore-Penrose generalized inverse of the specified matrix. + * + * @tparam T type of elements in the matrix + * @tparam n number of rows, can be Eigen::Dynamic + * @tparam m number of columns, can be Eigen::Dynamic + * + * @param m specified matrix + * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has + * size zero). + */ +template * = nullptr> +inline Eigen::Matrix, EigMat::RowsAtCompileTime, + EigMat::ColsAtCompileTime> +generalized_inverse(const EigMat& M) { + if (M.size() == 0) { + return {}; + } + + constexpr int n = EigMat::RowsAtCompileTime; + constexpr int m = EigMat::ColsAtCompileTime; + + if (n == m) { + return M.inverse(); + } + + bool transp(false); + constexpr int mn; + + if (n < m) { + transp = true; + mn = n; + Eigen::Matrix A = tcrossprod(G); + } else { + mn = m; + Eigen::Matrix A = = crossprod(G); + } + + A = add_diag(A, rep_vector(1e-10, mn)); + + Eigen::Matrix L = cholesky_decompose(A); + Eigen::Matrix M = chol2inv(L); + + if (transp){ + return transpose(G) * tcrossprod(L * M); + } else { + return tcrossprod(L * M) * transpose(G); + } +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/prim/fun/generalied_inverse_test.cpp b/test/unit/math/prim/fun/generalied_inverse_test.cpp new file mode 100644 index 00000000000..09d4c510d0f --- /dev/null +++ b/test/unit/math/prim/fun/generalied_inverse_test.cpp @@ -0,0 +1,17 @@ +#include +#include + +TEST(MathMatrixPrim, gen_inverse) { + stan::math::matrix_d m0(0, 0); + stan::math::matrix_d inv = stan::math::generalized_inverse(m0); + EXPECT_EQ(0, inv.rows()); + EXPECT_EQ(0, inv.cols()); +} + +TEST(MathMatrixPrim, inverse_exception) { + using stan::math::generalized_inverse; + + stan::math::matrix_d m1(2, 3); + m1 << 1, 2, 3, 4, 5, 6; + EXPECT_THROW(inverse(m1), std::invalid_argument); +} From 16698b220397e427697275ba22edaabb98ca5a7e Mon Sep 17 00:00:00 2001 From: spinkney Date: Tue, 1 Dec 2020 15:54:13 -0500 Subject: [PATCH 002/102] finished genealized_inverse test --- .../math/prim/fun/generalied_inverse_test.cpp | 17 ---------------- .../prim/fun/generalized_inverse_test.cpp | 20 +++++++++++++++++++ 2 files changed, 20 insertions(+), 17 deletions(-) delete mode 100644 test/unit/math/prim/fun/generalied_inverse_test.cpp create mode 100644 test/unit/math/prim/fun/generalized_inverse_test.cpp diff --git a/test/unit/math/prim/fun/generalied_inverse_test.cpp b/test/unit/math/prim/fun/generalied_inverse_test.cpp deleted file mode 100644 index 09d4c510d0f..00000000000 --- a/test/unit/math/prim/fun/generalied_inverse_test.cpp +++ /dev/null @@ -1,17 +0,0 @@ -#include -#include - -TEST(MathMatrixPrim, gen_inverse) { - stan::math::matrix_d m0(0, 0); - stan::math::matrix_d inv = stan::math::generalized_inverse(m0); - EXPECT_EQ(0, inv.rows()); - EXPECT_EQ(0, inv.cols()); -} - -TEST(MathMatrixPrim, inverse_exception) { - using stan::math::generalized_inverse; - - stan::math::matrix_d m1(2, 3); - m1 << 1, 2, 3, 4, 5, 6; - EXPECT_THROW(inverse(m1), std::invalid_argument); -} diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp new file mode 100644 index 00000000000..a0443047329 --- /dev/null +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -0,0 +1,20 @@ +#include +#include + +TEST(MathMatrixPrim, generalized_inverse) { + stan::math::matrix_d m0(0, 0); + stan::math::matrix_d ginv = stan::math::generalized_inverse(m0); + EXPECT_EQ(0, ginv.rows()); + EXPECT_EQ(0, ginv.cols()); +} + +TEST(MathMatrixPrim, inverse_correct) { + using stan::math::generalized_inverse; + + stan::math::matrix_d m1(3, 2); + m1 << 1, 2, 2, 4, 1, 2; + + stan::math::matrix_d m2(2, 3); + m2 << 1/30, 1/15, 1/30, 1/15, 2/15, 1/15; + EXPECT_MATRIX_FLOAT_EQ(generalized_inverse(m1), m2); +} From 593942245dcb9e41f90fa81ec1e5769af70273b2 Mon Sep 17 00:00:00 2001 From: spinkney Date: Tue, 1 Dec 2020 16:24:49 -0500 Subject: [PATCH 003/102] added generalized_inverse in fun.hpp --- stan/math/prim/fun.hpp | 1 + test/unit/math/prim/fun/generalized_inverse_test.cpp | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 04cc3990bc7..db436f213ad 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -102,6 +102,7 @@ #include #include #include +#include #include #include #include diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index a0443047329..47cff7230b7 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -1,14 +1,14 @@ #include #include -TEST(MathMatrixPrim, generalized_inverse) { +TEST(MathMatrixPrim, Zero) { stan::math::matrix_d m0(0, 0); stan::math::matrix_d ginv = stan::math::generalized_inverse(m0); EXPECT_EQ(0, ginv.rows()); EXPECT_EQ(0, ginv.cols()); } -TEST(MathMatrixPrim, inverse_correct) { +TEST(MathMatrixPrim, Equal) { using stan::math::generalized_inverse; stan::math::matrix_d m1(3, 2); From d3337c2e5fc0906c8967f7ecc03de591d9278d86 Mon Sep 17 00:00:00 2001 From: spinkney Date: Tue, 1 Dec 2020 16:41:06 -0500 Subject: [PATCH 004/102] fixed docs on generalized_inverse --- stan/math/prim/fun/generalized_inverse.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index ff2c26c9ebf..7779791e5b8 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -22,12 +22,11 @@ namespace math { * @tparam n number of rows, can be Eigen::Dynamic * @tparam m number of columns, can be Eigen::Dynamic * - * @param m specified matrix + * @param M specified matrix * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has * size zero). */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& M) { From fdc6c1b439197dca9a0e01ee5934021a0bdfb117 Mon Sep 17 00:00:00 2001 From: Katya Date: Tue, 1 Dec 2020 20:28:21 -0500 Subject: [PATCH 005/102] update --- stan/math/prim/fun/generalized_inverse.hpp | 48 +++++++++---------- .../prim/fun/generalized_inverse_test.cpp | 8 +++- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 7779791e5b8..858def0ca11 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -11,6 +11,8 @@ #include #include #include +#include +#include namespace stan { namespace math { @@ -29,40 +31,34 @@ namespace math { template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> -generalized_inverse(const EigMat& M) { - if (M.size() == 0) { +generalized_inverse(const EigMat& G) { + using value_t = value_type_t; + if (G.size() == 0) { return {}; } - constexpr int n = EigMat::RowsAtCompileTime; - constexpr int m = EigMat::ColsAtCompileTime; + const auto n = G.rows(); + const auto m = G.cols(); - if (n == m) { - return M.inverse(); - } - - bool transp(false); - constexpr int mn; + if (G.rows() == G.cols()) { + return G.inverse(); + } if (n < m) { - transp = true; - mn = n; - Eigen::Matrix A = tcrossprod(G); + Eigen::Matrix A = tcrossprod(G); + A.diagonal().array() += Eigen::Array::Constant(n, 1e-10); + Eigen::Matrix L = cholesky_decompose(A); + Eigen::Matrix L_inv = mdivide_left_tri_low(L); + Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); + return transpose(G) * tcrossprod(L * M); } else { - mn = m; - Eigen::Matrix A = = crossprod(G); + Eigen::Matrix A = crossprod(G); + A.diagonal().array() += Eigen::Array::Constant(m, 1e-10); + Eigen::Matrix L = cholesky_decompose(A); + Eigen::Matrix L_inv = mdivide_left_tri_low(L); + Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); + return tcrossprod(L * M) * transpose(G); } - - A = add_diag(A, rep_vector(1e-10, mn)); - - Eigen::Matrix L = cholesky_decompose(A); - Eigen::Matrix M = chol2inv(L); - - if (transp){ - return transpose(G) * tcrossprod(L * M); - } else { - return tcrossprod(L * M) * transpose(G); - } } } // namespace math diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 47cff7230b7..20cf92a934e 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -1,4 +1,5 @@ #include +#include #include TEST(MathMatrixPrim, Zero) { @@ -15,6 +16,9 @@ TEST(MathMatrixPrim, Equal) { m1 << 1, 2, 2, 4, 1, 2; stan::math::matrix_d m2(2, 3); - m2 << 1/30, 1/15, 1/30, 1/15, 2/15, 1/15; - EXPECT_MATRIX_FLOAT_EQ(generalized_inverse(m1), m2); + m2 << 1 / 30.0, 1 / 15.0, 1 / 30.0, 1 / 15.0, 2 / 15.0, 1 / 15.0; + stan::math::matrix_d m3 = generalized_inverse(m1); + for (int j = 0; j < m3.cols(); j++) + for (int i = 0; i < m3.rows(); i++) + EXPECT_NEAR(m2(i, j), m3(i, j), 1e-5); } From 1c8a51eb551775b06e4ef4d2ac9e0d31fcfc367c Mon Sep 17 00:00:00 2001 From: spinkney Date: Tue, 1 Dec 2020 13:49:45 -0500 Subject: [PATCH 006/102] Added generalized_inverse files Added generalized_inverse files --- stan/math/prim/fun/generalized_inverse.hpp | 72 +++++++++++++++++++ .../math/prim/fun/generalied_inverse_test.cpp | 17 +++++ 2 files changed, 89 insertions(+) create mode 100644 stan/math/prim/fun/generalized_inverse.hpp create mode 100644 test/unit/math/prim/fun/generalied_inverse_test.cpp diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp new file mode 100644 index 00000000000..ff2c26c9ebf --- /dev/null +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -0,0 +1,72 @@ +#ifndef STAN_MATH_PRIM_FUN_GENERALIZED_INVERSE_HPP +#define STAN_MATH_PRIM_FUN_GENERALIZED_INVERSE_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Returns the Moore-Penrose generalized inverse of the specified matrix. + * + * @tparam T type of elements in the matrix + * @tparam n number of rows, can be Eigen::Dynamic + * @tparam m number of columns, can be Eigen::Dynamic + * + * @param m specified matrix + * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has + * size zero). + */ +template * = nullptr> +inline Eigen::Matrix, EigMat::RowsAtCompileTime, + EigMat::ColsAtCompileTime> +generalized_inverse(const EigMat& M) { + if (M.size() == 0) { + return {}; + } + + constexpr int n = EigMat::RowsAtCompileTime; + constexpr int m = EigMat::ColsAtCompileTime; + + if (n == m) { + return M.inverse(); + } + + bool transp(false); + constexpr int mn; + + if (n < m) { + transp = true; + mn = n; + Eigen::Matrix A = tcrossprod(G); + } else { + mn = m; + Eigen::Matrix A = = crossprod(G); + } + + A = add_diag(A, rep_vector(1e-10, mn)); + + Eigen::Matrix L = cholesky_decompose(A); + Eigen::Matrix M = chol2inv(L); + + if (transp){ + return transpose(G) * tcrossprod(L * M); + } else { + return tcrossprod(L * M) * transpose(G); + } +} + +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/prim/fun/generalied_inverse_test.cpp b/test/unit/math/prim/fun/generalied_inverse_test.cpp new file mode 100644 index 00000000000..09d4c510d0f --- /dev/null +++ b/test/unit/math/prim/fun/generalied_inverse_test.cpp @@ -0,0 +1,17 @@ +#include +#include + +TEST(MathMatrixPrim, gen_inverse) { + stan::math::matrix_d m0(0, 0); + stan::math::matrix_d inv = stan::math::generalized_inverse(m0); + EXPECT_EQ(0, inv.rows()); + EXPECT_EQ(0, inv.cols()); +} + +TEST(MathMatrixPrim, inverse_exception) { + using stan::math::generalized_inverse; + + stan::math::matrix_d m1(2, 3); + m1 << 1, 2, 3, 4, 5, 6; + EXPECT_THROW(inverse(m1), std::invalid_argument); +} From 172f7e9beb9d7195de2b9d3868f8b9e3cc5d4f6a Mon Sep 17 00:00:00 2001 From: spinkney Date: Tue, 1 Dec 2020 15:54:13 -0500 Subject: [PATCH 007/102] finished genealized_inverse test --- .../math/prim/fun/generalied_inverse_test.cpp | 17 ---------------- .../prim/fun/generalized_inverse_test.cpp | 20 +++++++++++++++++++ 2 files changed, 20 insertions(+), 17 deletions(-) delete mode 100644 test/unit/math/prim/fun/generalied_inverse_test.cpp create mode 100644 test/unit/math/prim/fun/generalized_inverse_test.cpp diff --git a/test/unit/math/prim/fun/generalied_inverse_test.cpp b/test/unit/math/prim/fun/generalied_inverse_test.cpp deleted file mode 100644 index 09d4c510d0f..00000000000 --- a/test/unit/math/prim/fun/generalied_inverse_test.cpp +++ /dev/null @@ -1,17 +0,0 @@ -#include -#include - -TEST(MathMatrixPrim, gen_inverse) { - stan::math::matrix_d m0(0, 0); - stan::math::matrix_d inv = stan::math::generalized_inverse(m0); - EXPECT_EQ(0, inv.rows()); - EXPECT_EQ(0, inv.cols()); -} - -TEST(MathMatrixPrim, inverse_exception) { - using stan::math::generalized_inverse; - - stan::math::matrix_d m1(2, 3); - m1 << 1, 2, 3, 4, 5, 6; - EXPECT_THROW(inverse(m1), std::invalid_argument); -} diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp new file mode 100644 index 00000000000..a0443047329 --- /dev/null +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -0,0 +1,20 @@ +#include +#include + +TEST(MathMatrixPrim, generalized_inverse) { + stan::math::matrix_d m0(0, 0); + stan::math::matrix_d ginv = stan::math::generalized_inverse(m0); + EXPECT_EQ(0, ginv.rows()); + EXPECT_EQ(0, ginv.cols()); +} + +TEST(MathMatrixPrim, inverse_correct) { + using stan::math::generalized_inverse; + + stan::math::matrix_d m1(3, 2); + m1 << 1, 2, 2, 4, 1, 2; + + stan::math::matrix_d m2(2, 3); + m2 << 1/30, 1/15, 1/30, 1/15, 2/15, 1/15; + EXPECT_MATRIX_FLOAT_EQ(generalized_inverse(m1), m2); +} From 32ec5177a7b192e38b92f4550b35b39326e597de Mon Sep 17 00:00:00 2001 From: spinkney Date: Tue, 1 Dec 2020 16:24:49 -0500 Subject: [PATCH 008/102] added generalized_inverse in fun.hpp --- stan/math/prim/fun.hpp | 1 + test/unit/math/prim/fun/generalized_inverse_test.cpp | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index f4e51d5d75b..ef110d5fa34 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -102,6 +102,7 @@ #include #include #include +#include #include #include #include diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index a0443047329..47cff7230b7 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -1,14 +1,14 @@ #include #include -TEST(MathMatrixPrim, generalized_inverse) { +TEST(MathMatrixPrim, Zero) { stan::math::matrix_d m0(0, 0); stan::math::matrix_d ginv = stan::math::generalized_inverse(m0); EXPECT_EQ(0, ginv.rows()); EXPECT_EQ(0, ginv.cols()); } -TEST(MathMatrixPrim, inverse_correct) { +TEST(MathMatrixPrim, Equal) { using stan::math::generalized_inverse; stan::math::matrix_d m1(3, 2); From 270dde66a2d2252555b9ff7ddaf6d473dc79f69f Mon Sep 17 00:00:00 2001 From: spinkney Date: Tue, 1 Dec 2020 16:41:06 -0500 Subject: [PATCH 009/102] fixed docs on generalized_inverse --- stan/math/prim/fun/generalized_inverse.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index ff2c26c9ebf..7779791e5b8 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -22,12 +22,11 @@ namespace math { * @tparam n number of rows, can be Eigen::Dynamic * @tparam m number of columns, can be Eigen::Dynamic * - * @param m specified matrix + * @param M specified matrix * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has * size zero). */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& M) { From 037ca41324fd071e9e02748b3d357a15c70bcc8b Mon Sep 17 00:00:00 2001 From: Katya Date: Tue, 1 Dec 2020 20:28:21 -0500 Subject: [PATCH 010/102] update --- stan/math/prim/fun/generalized_inverse.hpp | 48 +++++++++---------- .../prim/fun/generalized_inverse_test.cpp | 8 +++- 2 files changed, 28 insertions(+), 28 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 7779791e5b8..858def0ca11 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -11,6 +11,8 @@ #include #include #include +#include +#include namespace stan { namespace math { @@ -29,40 +31,34 @@ namespace math { template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> -generalized_inverse(const EigMat& M) { - if (M.size() == 0) { +generalized_inverse(const EigMat& G) { + using value_t = value_type_t; + if (G.size() == 0) { return {}; } - constexpr int n = EigMat::RowsAtCompileTime; - constexpr int m = EigMat::ColsAtCompileTime; + const auto n = G.rows(); + const auto m = G.cols(); - if (n == m) { - return M.inverse(); - } - - bool transp(false); - constexpr int mn; + if (G.rows() == G.cols()) { + return G.inverse(); + } if (n < m) { - transp = true; - mn = n; - Eigen::Matrix A = tcrossprod(G); + Eigen::Matrix A = tcrossprod(G); + A.diagonal().array() += Eigen::Array::Constant(n, 1e-10); + Eigen::Matrix L = cholesky_decompose(A); + Eigen::Matrix L_inv = mdivide_left_tri_low(L); + Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); + return transpose(G) * tcrossprod(L * M); } else { - mn = m; - Eigen::Matrix A = = crossprod(G); + Eigen::Matrix A = crossprod(G); + A.diagonal().array() += Eigen::Array::Constant(m, 1e-10); + Eigen::Matrix L = cholesky_decompose(A); + Eigen::Matrix L_inv = mdivide_left_tri_low(L); + Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); + return tcrossprod(L * M) * transpose(G); } - - A = add_diag(A, rep_vector(1e-10, mn)); - - Eigen::Matrix L = cholesky_decompose(A); - Eigen::Matrix M = chol2inv(L); - - if (transp){ - return transpose(G) * tcrossprod(L * M); - } else { - return tcrossprod(L * M) * transpose(G); - } } } // namespace math diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 47cff7230b7..20cf92a934e 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -1,4 +1,5 @@ #include +#include #include TEST(MathMatrixPrim, Zero) { @@ -15,6 +16,9 @@ TEST(MathMatrixPrim, Equal) { m1 << 1, 2, 2, 4, 1, 2; stan::math::matrix_d m2(2, 3); - m2 << 1/30, 1/15, 1/30, 1/15, 2/15, 1/15; - EXPECT_MATRIX_FLOAT_EQ(generalized_inverse(m1), m2); + m2 << 1 / 30.0, 1 / 15.0, 1 / 30.0, 1 / 15.0, 2 / 15.0, 1 / 15.0; + stan::math::matrix_d m3 = generalized_inverse(m1); + for (int j = 0; j < m3.cols(); j++) + for (int i = 0; i < m3.rows(); i++) + EXPECT_NEAR(m2(i, j), m3(i, j), 1e-5); } From 21cbcc1035cdf92db77b8fc18518510ce2ad7f0d Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 2 Dec 2020 11:47:33 +0000 Subject: [PATCH 011/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 11 ++++++----- test/unit/math/prim/fun/generalized_inverse_test.cpp | 2 +- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 858def0ca11..267600044fe 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -25,10 +25,11 @@ namespace math { * @tparam m number of columns, can be Eigen::Dynamic * * @param M specified matrix - * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has - * size zero). + * @return Generalized inverse of the matrix (an empty matrix if the specified + * matrix has size zero). */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G) { @@ -41,8 +42,8 @@ generalized_inverse(const EigMat& G) { const auto m = G.cols(); if (G.rows() == G.cols()) { - return G.inverse(); - } + return G.inverse(); + } if (n < m) { Eigen::Matrix A = tcrossprod(G); diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 20cf92a934e..4642a6d52fd 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -17,7 +17,7 @@ TEST(MathMatrixPrim, Equal) { stan::math::matrix_d m2(2, 3); m2 << 1 / 30.0, 1 / 15.0, 1 / 30.0, 1 / 15.0, 2 / 15.0, 1 / 15.0; - stan::math::matrix_d m3 = generalized_inverse(m1); + stan::math::matrix_d m3 = generalized_inverse(m1); for (int j = 0; j < m3.cols(); j++) for (int i = 0; i < m3.rows(); i++) EXPECT_NEAR(m2(i, j), m3(i, j), 1e-5); From aeb78b51c0629ad47b83886a505c994f0be8a743 Mon Sep 17 00:00:00 2001 From: spinkney Date: Wed, 2 Dec 2020 06:59:31 -0500 Subject: [PATCH 012/102] Fixed function doc --- stan/math/prim/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 267600044fe..540a3c95684 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -24,7 +24,7 @@ namespace math { * @tparam n number of rows, can be Eigen::Dynamic * @tparam m number of columns, can be Eigen::Dynamic * - * @param M specified matrix + * @param G specified matrix * @return Generalized inverse of the matrix (an empty matrix if the specified * matrix has size zero). */ From 2192b906ec5001d3e3167fd47d97c1969416b784 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 2 Dec 2020 08:28:59 -0500 Subject: [PATCH 013/102] Update test/unit/math/prim/fun/generalized_inverse_test.cpp MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Tadej Ciglarič --- test/unit/math/prim/fun/generalized_inverse_test.cpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 4642a6d52fd..8e079504f84 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -18,7 +18,5 @@ TEST(MathMatrixPrim, Equal) { stan::math::matrix_d m2(2, 3); m2 << 1 / 30.0, 1 / 15.0, 1 / 30.0, 1 / 15.0, 2 / 15.0, 1 / 15.0; stan::math::matrix_d m3 = generalized_inverse(m1); - for (int j = 0; j < m3.cols(); j++) - for (int i = 0; i < m3.rows(); i++) - EXPECT_NEAR(m2(i, j), m3(i, j), 1e-5); + EXPECT_MATRIX_NEAR(m2, m3, 1e-5); } From bc582000a20c5fbc908fe33a3fb5a84c741cae93 Mon Sep 17 00:00:00 2001 From: spinkney Date: Wed, 2 Dec 2020 08:35:08 -0500 Subject: [PATCH 014/102] Changed -1 to Eigen::Dynamic --- stan/math/prim/fun/generalized_inverse.hpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 540a3c95684..2803b71db1f 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -46,18 +46,18 @@ generalized_inverse(const EigMat& G) { } if (n < m) { - Eigen::Matrix A = tcrossprod(G); - A.diagonal().array() += Eigen::Array::Constant(n, 1e-10); - Eigen::Matrix L = cholesky_decompose(A); - Eigen::Matrix L_inv = mdivide_left_tri_low(L); - Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); + Eigen::Matrix A = tcrossprod(G); + A.diagonal().array() += Eigen::Array::Constant(n, 1e-10); + Eigen::Matrix L = cholesky_decompose(A); + Eigen::Matrix L_inv = mdivide_left_tri_low(L); + Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); return transpose(G) * tcrossprod(L * M); } else { - Eigen::Matrix A = crossprod(G); - A.diagonal().array() += Eigen::Array::Constant(m, 1e-10); - Eigen::Matrix L = cholesky_decompose(A); - Eigen::Matrix L_inv = mdivide_left_tri_low(L); - Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); + Eigen::Matrix A = crossprod(G); + A.diagonal().array() += Eigen::Array::Constant(m, 1e-10); + Eigen::Matrix L = cholesky_decompose(A); + Eigen::Matrix L_inv = mdivide_left_tri_low(L); + Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); return tcrossprod(L * M) * transpose(G); } } From 5ce32766f40505a3d2d9c1582f7d4a1605358db0 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 2 Dec 2020 13:42:20 +0000 Subject: [PATCH 015/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 24 ++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 2803b71db1f..fa4972fc667 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -47,17 +47,25 @@ generalized_inverse(const EigMat& G) { if (n < m) { Eigen::Matrix A = tcrossprod(G); - A.diagonal().array() += Eigen::Array::Constant(n, 1e-10); - Eigen::Matrix L = cholesky_decompose(A); - Eigen::Matrix L_inv = mdivide_left_tri_low(L); - Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); + A.diagonal().array() + += Eigen::Array::Constant(n, 1e-10); + Eigen::Matrix L + = cholesky_decompose(A); + Eigen::Matrix L_inv + = mdivide_left_tri_low(L); + Eigen::Matrix M + = multiply_lower_tri_self_transpose(L_inv); return transpose(G) * tcrossprod(L * M); } else { Eigen::Matrix A = crossprod(G); - A.diagonal().array() += Eigen::Array::Constant(m, 1e-10); - Eigen::Matrix L = cholesky_decompose(A); - Eigen::Matrix L_inv = mdivide_left_tri_low(L); - Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); + A.diagonal().array() + += Eigen::Array::Constant(m, 1e-10); + Eigen::Matrix L + = cholesky_decompose(A); + Eigen::Matrix L_inv + = mdivide_left_tri_low(L); + Eigen::Matrix M + = multiply_lower_tri_self_transpose(L_inv); return tcrossprod(L * M) * transpose(G); } } From a862f36bfec33b9e50d22b476d70dbcdb6c0624c Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 2 Dec 2020 12:12:24 -0500 Subject: [PATCH 016/102] Using chol2inv adds a bit more numerical stability --- stan/math/prim/fun/generalized_inverse.hpp | 19 +++++++------------ .../prim/fun/generalized_inverse_test.cpp | 2 +- 2 files changed, 8 insertions(+), 13 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 2803b71db1f..59de64aeea8 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -5,14 +5,11 @@ #include #include #include -#include -#include #include #include #include #include -#include -#include +#include namespace stan { namespace math { @@ -47,18 +44,16 @@ generalized_inverse(const EigMat& G) { if (n < m) { Eigen::Matrix A = tcrossprod(G); - A.diagonal().array() += Eigen::Array::Constant(n, 1e-10); + A.diagonal().array() += Eigen::Array::Constant(n, 3.712035-7); Eigen::Matrix L = cholesky_decompose(A); - Eigen::Matrix L_inv = mdivide_left_tri_low(L); - Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); - return transpose(G) * tcrossprod(L * M); + Eigen::Matrix M = chol2inv(L); + return transpose(G) * quad_form(A, M); } else { Eigen::Matrix A = crossprod(G); - A.diagonal().array() += Eigen::Array::Constant(m, 1e-10); + A.diagonal().array() += Eigen::Array::Constant(m, 3.712035e-7); Eigen::Matrix L = cholesky_decompose(A); - Eigen::Matrix L_inv = mdivide_left_tri_low(L); - Eigen::Matrix M = multiply_lower_tri_self_transpose(L_inv); - return tcrossprod(L * M) * transpose(G); + Eigen::Matrix M = chol2inv(L); + return quad_form(A, M) * transpose(G); } } diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 8e079504f84..1ab4bbebea1 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -18,5 +18,5 @@ TEST(MathMatrixPrim, Equal) { stan::math::matrix_d m2(2, 3); m2 << 1 / 30.0, 1 / 15.0, 1 / 30.0, 1 / 15.0, 2 / 15.0, 1 / 15.0; stan::math::matrix_d m3 = generalized_inverse(m1); - EXPECT_MATRIX_NEAR(m2, m3, 1e-5); + EXPECT_MATRIX_NEAR(m2, m3, 1e-9); } From fedfd8a94946d4a73427bb3718cfa18cbf329b42 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 2 Dec 2020 17:25:30 +0000 Subject: [PATCH 017/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 59de64aeea8..3e59b247a56 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -44,14 +44,18 @@ generalized_inverse(const EigMat& G) { if (n < m) { Eigen::Matrix A = tcrossprod(G); - A.diagonal().array() += Eigen::Array::Constant(n, 3.712035-7); - Eigen::Matrix L = cholesky_decompose(A); + A.diagonal().array() + += Eigen::Array::Constant(n, 3.712035 - 7); + Eigen::Matrix L + = cholesky_decompose(A); Eigen::Matrix M = chol2inv(L); return transpose(G) * quad_form(A, M); } else { Eigen::Matrix A = crossprod(G); - A.diagonal().array() += Eigen::Array::Constant(m, 3.712035e-7); - Eigen::Matrix L = cholesky_decompose(A); + A.diagonal().array() + += Eigen::Array::Constant(m, 3.712035e-7); + Eigen::Matrix L + = cholesky_decompose(A); Eigen::Matrix M = chol2inv(L); return quad_form(A, M) * transpose(G); } From 0e6a898b443b068e5bd48c528cf58783a53410a3 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 2 Dec 2020 12:27:53 -0500 Subject: [PATCH 018/102] fixed small typo --- stan/math/prim/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 59de64aeea8..7f374c00f3d 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -44,7 +44,7 @@ generalized_inverse(const EigMat& G) { if (n < m) { Eigen::Matrix A = tcrossprod(G); - A.diagonal().array() += Eigen::Array::Constant(n, 3.712035-7); + A.diagonal().array() += Eigen::Array::Constant(n, 3.712035e-7); Eigen::Matrix L = cholesky_decompose(A); Eigen::Matrix M = chol2inv(L); return transpose(G) * quad_form(A, M); From e2fc0f9eef46beb6d609ef50cc7afa3c0f1bf87e Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 2 Dec 2020 17:29:48 +0000 Subject: [PATCH 019/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 4edf550156d..5711289dc69 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -44,8 +44,10 @@ generalized_inverse(const EigMat& G) { if (n < m) { Eigen::Matrix A = tcrossprod(G); - A.diagonal().array() += Eigen::Array::Constant(n, 3.712035e-7); - Eigen::Matrix L = cholesky_decompose(A); + A.diagonal().array() + += Eigen::Array::Constant(n, 3.712035e-7); + Eigen::Matrix L + = cholesky_decompose(A); Eigen::Matrix M = chol2inv(L); return transpose(G) * quad_form(A, M); } else { From 5601b009b5cbc55acc39e7a52b97b2519ef14c3a Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 2 Dec 2020 14:06:54 -0500 Subject: [PATCH 020/102] optimizing the algebra --- stan/math/prim/fun/generalized_inverse.hpp | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 5711289dc69..7e9d9e66568 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -10,6 +10,9 @@ #include #include #include +#include +#include +#include namespace stan { namespace math { @@ -45,19 +48,13 @@ generalized_inverse(const EigMat& G) { if (n < m) { Eigen::Matrix A = tcrossprod(G); A.diagonal().array() - += Eigen::Array::Constant(n, 3.712035e-7); - Eigen::Matrix L - = cholesky_decompose(A); - Eigen::Matrix M = chol2inv(L); - return transpose(G) * quad_form(A, M); + += Eigen::Array::Constant(n, 1e-8); + return transpose(mdivide_left_spd(A, G)); } else { Eigen::Matrix A = crossprod(G); A.diagonal().array() - += Eigen::Array::Constant(m, 3.712035e-7); - Eigen::Matrix L - = cholesky_decompose(A); - Eigen::Matrix M = chol2inv(L); - return quad_form(A, M) * transpose(G); + += Eigen::Array::Constant(m, 1e-8); + return transpose(mdivide_right_spd(G, A)); } } From 9066a94ea756232259346082a7cbdfa29a36c02c Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 2 Dec 2020 19:08:15 +0000 Subject: [PATCH 021/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 7e9d9e66568..3fae6345dca 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -54,7 +54,7 @@ generalized_inverse(const EigMat& G) { Eigen::Matrix A = crossprod(G); A.diagonal().array() += Eigen::Array::Constant(m, 1e-8); - return transpose(mdivide_right_spd(G, A)); + return transpose(mdivide_right_spd(G, A)); } } From 86d67ec691d71c3a473bd93d464af78e311237eb Mon Sep 17 00:00:00 2001 From: spinkney Date: Wed, 2 Dec 2020 16:48:46 -0500 Subject: [PATCH 022/102] Added jitter alias and extra test --- stan/math/prim/fun/generalized_inverse.hpp | 43 +++++++++++++++++-- .../prim/fun/generalized_inverse_test.cpp | 14 +++++- 2 files changed, 53 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 3fae6345dca..6cf274b6988 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -13,6 +13,7 @@ #include #include #include +#include namespace stan { namespace math { @@ -42,18 +43,54 @@ generalized_inverse(const EigMat& G) { const auto m = G.cols(); if (G.rows() == G.cols()) { - return G.inverse(); + return inverse(G); + } + + if (n < m) { + return transpose(mdivide_left_spd(A, tcrossprod(G))); + } else { + return transpose(mdivide_right_spd(G, crossprod(A))); + } +} + +/** + * Returns the Moore-Penrose generalized inverse of the specified matrix. + * + * @tparam T type of elements in the matrix + * @tparam n number of rows, can be Eigen::Dynamic + * @tparam m number of columns, can be Eigen::Dynamic + * + * @param G specified matrix + * @param a diagonal jitter + * @return Generalized inverse of the matrix (an empty matrix if the specified + * matrix has size zero). + */ +template * = nullptr> +inline Eigen::Matrix, EigMat::RowsAtCompileTime, + EigMat::ColsAtCompileTime> +generalized_inverse(const EigMat& G, double a) { + using value_t = value_type_t; + if (G.size() == 0) { + return {}; + } + + const auto n = G.rows(); + const auto m = G.cols(); + + if (G.rows() == G.cols()) { + return inverse(G); } if (n < m) { Eigen::Matrix A = tcrossprod(G); A.diagonal().array() - += Eigen::Array::Constant(n, 1e-8); + += Eigen::Array::Constant(n, a); return transpose(mdivide_left_spd(A, G)); } else { Eigen::Matrix A = crossprod(G); A.diagonal().array() - += Eigen::Array::Constant(m, 1e-8); + += Eigen::Array::Constant(m, a); return transpose(mdivide_right_spd(G, A)); } } diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 1ab4bbebea1..402f2377f1a 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -9,7 +9,7 @@ TEST(MathMatrixPrim, Zero) { EXPECT_EQ(0, ginv.cols()); } -TEST(MathMatrixPrim, Equal) { +TEST(MathMatrixPrim, EqualJitter) { using stan::math::generalized_inverse; stan::math::matrix_d m1(3, 2); @@ -17,6 +17,18 @@ TEST(MathMatrixPrim, Equal) { stan::math::matrix_d m2(2, 3); m2 << 1 / 30.0, 1 / 15.0, 1 / 30.0, 1 / 15.0, 2 / 15.0, 1 / 15.0; + stan::math::matrix_d m3 = generalized_inverse(m1, 1e-8); + EXPECT_MATRIX_NEAR(m2, m3, 1e-9); +} + +TEST(MathMatrixPrim, Equal) { + using stan::math::generalized_inverse; + + stan::math::matrix_d m1(2, 3); + m1 << 1, 2, 3, 4, 5, 6; + + stan::math::matrix_d m2(2, 3); + m2 << -4.0 / 3.0, 13.0 / 12.0, -1.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0, -5.0 / 12.0; stan::math::matrix_d m3 = generalized_inverse(m1); EXPECT_MATRIX_NEAR(m2, m3, 1e-9); } From 664c4fa909d5281909f4dabb5868b9d6e8669771 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Wed, 2 Dec 2020 21:49:48 +0000 Subject: [PATCH 023/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 6cf274b6988..884da0658a8 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -66,7 +66,7 @@ generalized_inverse(const EigMat& G) { * matrix has size zero). */ template * = nullptr> + require_eigen_vt* = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G, double a) { From 0b7732668bb6356926d7f5acfab82bd6ccd22acf Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 2 Dec 2020 16:54:00 -0500 Subject: [PATCH 024/102] typos --- stan/math/prim/fun/generalized_inverse.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 884da0658a8..6d94833fd6a 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -47,9 +47,9 @@ generalized_inverse(const EigMat& G) { } if (n < m) { - return transpose(mdivide_left_spd(A, tcrossprod(G))); + return transpose(mdivide_left_spd(tcrossprod(G), G)); } else { - return transpose(mdivide_right_spd(G, crossprod(A))); + return transpose(mdivide_right_spd(G, crossprod(G))); } } From 2792baaa17af8a62962d578d7415f785a76bcca9 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 2 Dec 2020 16:57:23 -0500 Subject: [PATCH 025/102] test fix --- test/unit/math/prim/fun/generalized_inverse_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 402f2377f1a..f4e2cc09ecc 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -27,7 +27,7 @@ TEST(MathMatrixPrim, Equal) { stan::math::matrix_d m1(2, 3); m1 << 1, 2, 3, 4, 5, 6; - stan::math::matrix_d m2(2, 3); + stan::math::matrix_d m2(3, 2); m2 << -4.0 / 3.0, 13.0 / 12.0, -1.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0, -5.0 / 12.0; stan::math::matrix_d m3 = generalized_inverse(m1); EXPECT_MATRIX_NEAR(m2, m3, 1e-9); From 9788f2aeb4801419cd77980e30962d07e75e5359 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Wed, 2 Dec 2020 17:01:48 -0500 Subject: [PATCH 026/102] fixed test matrix --- test/unit/math/prim/fun/generalized_inverse_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index f4e2cc09ecc..477fde2eea5 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -25,7 +25,7 @@ TEST(MathMatrixPrim, Equal) { using stan::math::generalized_inverse; stan::math::matrix_d m1(2, 3); - m1 << 1, 2, 3, 4, 5, 6; + m1 << 1, 3, 5, 2, 4, 6; stan::math::matrix_d m2(3, 2); m2 << -4.0 / 3.0, 13.0 / 12.0, -1.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0, -5.0 / 12.0; From 836ba0ed8ba1833d928b03418dbf6a4a7f4f6408 Mon Sep 17 00:00:00 2001 From: spinkney Date: Thu, 3 Dec 2020 11:31:41 -0500 Subject: [PATCH 027/102] clean up docs --- stan/math/prim/fun/generalized_inverse.hpp | 32 ++++++++++++++++------ 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 6d94833fd6a..81703724d96 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -4,12 +4,10 @@ #include #include #include -#include #include #include #include #include -#include #include #include #include @@ -21,13 +19,21 @@ namespace math { /** * Returns the Moore-Penrose generalized inverse of the specified matrix. * - * @tparam T type of elements in the matrix - * @tparam n number of rows, can be Eigen::Dynamic - * @tparam m number of columns, can be Eigen::Dynamic + * The method is based on the Cholesky computation of the transform as specified in + * + *
  • Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse Matrices. + * arXiv 0804.4809
+ * + * @tparam EigMat type of the matrix (must be derived from \c Eigen::MatrixBase) * * @param G specified matrix * @return Generalized inverse of the matrix (an empty matrix if the specified * matrix has size zero). + * @note Because the method inverts a SPD matrix internally that interal matrix may result + in small numerical issues that result in a non-SPD error. There are two + * generalized_inverse functions, one that uses one input matrix (this one) + * and another that works with an input matrix and a small jitter to the diagonal of the internal SPD + * matrix. */ template * = nullptr> @@ -56,14 +62,22 @@ generalized_inverse(const EigMat& G) { /** * Returns the Moore-Penrose generalized inverse of the specified matrix. * - * @tparam T type of elements in the matrix - * @tparam n number of rows, can be Eigen::Dynamic - * @tparam m number of columns, can be Eigen::Dynamic + * The method is based on the Cholesky computation of the transform as specified in + * + *
  • Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse Matrices. + * arXiv 0804.4809
+ * + * @tparam EigMat type of the matrix (must be derived from \c Eigen::MatrixBase) * * @param G specified matrix - * @param a diagonal jitter + * @param a real constant * @return Generalized inverse of the matrix (an empty matrix if the specified * matrix has size zero). + * @note Because the method inverts a SPD matrix internally that interal matrix may result + in small numerical issues that result in a non-SPD error. There are two + * generalized_inverse functions, one that uses one input matrix + * and another (this one) that works with an input matrix and a small jitter to the diagonal of the internal SPD + * matrix. */ template * = nullptr> From 6718fe2cb11b4c69d70af603e9624cca4ee51c72 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 3 Dec 2020 16:38:25 +0000 Subject: [PATCH 028/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 33 ++++++++++++++-------- 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 81703724d96..b24cbbb1b72 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -19,9 +19,11 @@ namespace math { /** * Returns the Moore-Penrose generalized inverse of the specified matrix. * - * The method is based on the Cholesky computation of the transform as specified in + * The method is based on the Cholesky computation of the transform as specified + in * - *
  • Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse Matrices. + *
    • Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse + Matrices. * arXiv 0804.4809
    * * @tparam EigMat type of the matrix (must be derived from \c Eigen::MatrixBase) @@ -29,10 +31,13 @@ namespace math { * @param G specified matrix * @return Generalized inverse of the matrix (an empty matrix if the specified * matrix has size zero). - * @note Because the method inverts a SPD matrix internally that interal matrix may result - in small numerical issues that result in a non-SPD error. There are two - * generalized_inverse functions, one that uses one input matrix (this one) - * and another that works with an input matrix and a small jitter to the diagonal of the internal SPD + * @note Because the method inverts a SPD matrix internally that interal matrix + may result in small numerical issues that result in a non-SPD error. There are + two + * generalized_inverse functions, one that uses one input matrix + (this one) + * and another that works with an input matrix and a small jitter to the + diagonal of the internal SPD * matrix. */ template
  • Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse Matrices. + *
    • Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse + Matrices. * arXiv 0804.4809
    * * @tparam EigMat type of the matrix (must be derived from \c Eigen::MatrixBase) @@ -73,10 +80,12 @@ generalized_inverse(const EigMat& G) { * @param a real constant * @return Generalized inverse of the matrix (an empty matrix if the specified * matrix has size zero). - * @note Because the method inverts a SPD matrix internally that interal matrix may result - in small numerical issues that result in a non-SPD error. There are two - * generalized_inverse functions, one that uses one input matrix - * and another (this one) that works with an input matrix and a small jitter to the diagonal of the internal SPD + * @note Because the method inverts a SPD matrix internally that interal matrix + may result in small numerical issues that result in a non-SPD error. There are + two + * generalized_inverse functions, one that uses one input matrix + * and another (this one) that works with an input matrix and a small jitter to + the diagonal of the internal SPD * matrix. */ template Date: Thu, 3 Dec 2020 12:52:46 -0500 Subject: [PATCH 029/102] adding mix test --- .../math/mix/fun/generalized_inverse_test.cpp | 61 +++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 test/unit/math/mix/fun/generalized_inverse_test.cpp diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp new file mode 100644 index 00000000000..c2af9bfc3fb --- /dev/null +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -0,0 +1,61 @@ +#include +#include + +TEST(mathMixMatFun, generalized_inverse) { + auto f = [](const auto& x) { return stan::math::generalized_inverse(x); }; + + Eigen::MatrixXd t(0, 0); + stan::test::expect_ad(f, t); + stan::test::expect_ad_matvar(f, t); + + Eigen::MatrixXd u(1, 1); + u << 2; + stan::test::expect_ad(f, u); + stan::test::expect_ad_matvar(f, u); + + Eigen::MatrixXd v(2, 3); + v << 1, 3, 5, 2, 4, 6; + stan::test::expect_ad(f, v); + stan::test::expect_ad_matvar(f, v); + v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; + stan::test::expect_ad(f, v); + stan::test::expect_ad_matvar(f, v); + + // issues around zero require looser tolerances for hessians + stan::test::ad_tolerances tols; + tols.hessian_hessian_ = 0.01; + tols.hessian_fvar_hessian_ = 0.01; + + Eigen::MatrixXd w(3, 4); + w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; + stan::test::expect_ad(tols, f, w); + stan::test::expect_ad_matvar(tols, f, w); + + // even lower tolerance, again for cases around zero + stan::test::ad_tolerances tols2; + tols2.hessian_hessian_ = 3.0; + tols2.hessian_fvar_hessian_ = 3.0; + + Eigen::MatrixXd x(4, 4); + x << 2, 3, 4, 5, 9, -1, 2, 2, 4, 3, 7, -1, 0, 1, 19, 112; + stan::test::expect_ad(tols2, f, x); + stan::test::expect_ad_matvar(tols2, f, x); + + Eigen::MatrixXd z(2, 2); + z << 1, 2, 5, std::numeric_limits::quiet_NaN(); + EXPECT_NO_THROW(stan::math::generalized_inverse(z)); + + // autodiff throws, so following fails (throw behavior must match to pass) + // stan::test::expect_ad(f, z); + + Eigen::MatrixXd a(2, 2); + a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); + stan::test::expect_ad(f, a); + stan::test::expect_ad_matvar(f, a); + + // singular matrix, should use the + // alias to input small amount of jitter on the diagonal + Eigen::MatrixXd m(3, 2); + m << 1, 2, 2, 4, 1, 2;; + EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); +} \ No newline at end of file From 192f4e2670db8acd58b9a8344bee471f3aff7d59 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 3 Dec 2020 17:53:37 +0000 Subject: [PATCH 030/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index c2af9bfc3fb..5b9da351495 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -53,9 +53,10 @@ TEST(mathMixMatFun, generalized_inverse) { stan::test::expect_ad(f, a); stan::test::expect_ad_matvar(f, a); - // singular matrix, should use the + // singular matrix, should use the // alias to input small amount of jitter on the diagonal Eigen::MatrixXd m(3, 2); - m << 1, 2, 2, 4, 1, 2;; + m << 1, 2, 2, 4, 1, 2; + ; EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); } \ No newline at end of file From 248d33f139bf1ae022666de3b9a61a2979d70546 Mon Sep 17 00:00:00 2001 From: spinkney Date: Thu, 3 Dec 2020 13:47:59 -0500 Subject: [PATCH 031/102] typo --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index c2af9bfc3fb..6b578ebfd54 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -56,6 +56,6 @@ TEST(mathMixMatFun, generalized_inverse) { // singular matrix, should use the // alias to input small amount of jitter on the diagonal Eigen::MatrixXd m(3, 2); - m << 1, 2, 2, 4, 1, 2;; + m << 1, 2, 2, 4, 1, 2; EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); -} \ No newline at end of file +} From 26d018e879aa3fe74fce87f06f2f25de29a17022 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 3 Dec 2020 18:56:54 +0000 Subject: [PATCH 032/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 779acfa4527..96a51873edf 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -56,6 +56,6 @@ TEST(mathMixMatFun, generalized_inverse) { // singular matrix, should use the // alias to input small amount of jitter on the diagonal Eigen::MatrixXd m(3, 2); - m << 1, 2, 2, 4, 1, 2; + m << 1, 2, 2, 4, 1, 2; EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); } From a1024b6772f24262cf8949a38b155d7cdd5f7a8f Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Thu, 3 Dec 2020 17:24:40 -0500 Subject: [PATCH 033/102] Merge remote-tracking branch 'origin/feature/ginv' into feature/ginv --- stan/math/prim/fun/generalized_inverse.hpp | 9 ++-- .../math/mix/fun/generalized_inverse_test.cpp | 42 +++++++++---------- 2 files changed, 26 insertions(+), 25 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index b24cbbb1b72..3c59831a6b6 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -41,7 +41,7 @@ namespace math { * matrix. */ template * = nullptr> + require_eigen_t* = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G) { @@ -88,11 +88,12 @@ generalized_inverse(const EigMat& G) { the diagonal of the internal SPD * matrix. */ -template * = nullptr> +template * = nullptr, + require_stan_scalar_t* = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> -generalized_inverse(const EigMat& G, double a) { +generalized_inverse(const EigMat& G, const Scal& a) { using value_t = value_type_t; if (G.size() == 0) { return {}; diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 96a51873edf..778709456ae 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -1,37 +1,37 @@ #include -#include +#include TEST(mathMixMatFun, generalized_inverse) { auto f = [](const auto& x) { return stan::math::generalized_inverse(x); }; Eigen::MatrixXd t(0, 0); stan::test::expect_ad(f, t); - stan::test::expect_ad_matvar(f, t); +// stan::test::expect_ad_matvar(f, t); Eigen::MatrixXd u(1, 1); u << 2; stan::test::expect_ad(f, u); - stan::test::expect_ad_matvar(f, u); - +// stan::test::expect_ad_matvar(f, u); +// Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; stan::test::expect_ad(f, v); - stan::test::expect_ad_matvar(f, v); +// stan::test::expect_ad_matvar(f, v); v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; stan::test::expect_ad(f, v); - stan::test::expect_ad_matvar(f, v); +// stan::test::expect_ad_matvar(f, v); - // issues around zero require looser tolerances for hessians - stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 0.01; - tols.hessian_fvar_hessian_ = 0.01; +// issues around zero require looser tolerances for hessians +// stan::test::ad_tolerances tols; +// tols.hessian_hessian_ = 0.01; +// tols.hessian_fvar_hessian_ = 0.01; - Eigen::MatrixXd w(3, 4); - w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; - stan::test::expect_ad(tols, f, w); - stan::test::expect_ad_matvar(tols, f, w); +// Eigen::MatrixXd w(3, 4); +// w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; +// stan::test::expect_ad(tols, f, w); +// stan::test::expect_ad_matvar(tols, f, w); - // even lower tolerance, again for cases around zero +// even lower tolerance, again for cases around zero stan::test::ad_tolerances tols2; tols2.hessian_hessian_ = 3.0; tols2.hessian_fvar_hessian_ = 3.0; @@ -39,22 +39,22 @@ TEST(mathMixMatFun, generalized_inverse) { Eigen::MatrixXd x(4, 4); x << 2, 3, 4, 5, 9, -1, 2, 2, 4, 3, 7, -1, 0, 1, 19, 112; stan::test::expect_ad(tols2, f, x); - stan::test::expect_ad_matvar(tols2, f, x); +// stan::test::expect_ad_matvar(tols2, f, x); Eigen::MatrixXd z(2, 2); z << 1, 2, 5, std::numeric_limits::quiet_NaN(); EXPECT_NO_THROW(stan::math::generalized_inverse(z)); - // autodiff throws, so following fails (throw behavior must match to pass) - // stan::test::expect_ad(f, z); +// autodiff throws, so following fails (throw behavior must match to pass) +// stan::test::expect_ad(f, z); Eigen::MatrixXd a(2, 2); a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); stan::test::expect_ad(f, a); - stan::test::expect_ad_matvar(f, a); +// stan::test::expect_ad_matvar(f, a); - // singular matrix, should use the - // alias to input small amount of jitter on the diagonal +// singular matrix, should use the +// alias to input small amount of jitter on the diagonal Eigen::MatrixXd m(3, 2); m << 1, 2, 2, 4, 1, 2; EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); From 5a21f25d0d661a69bbf04c5cf27b38879690a220 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 3 Dec 2020 22:31:29 +0000 Subject: [PATCH 034/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 8 ++-- .../math/mix/fun/generalized_inverse_test.cpp | 40 +++++++++---------- 2 files changed, 23 insertions(+), 25 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 3c59831a6b6..add6ccd3c61 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -40,8 +40,7 @@ namespace math { diagonal of the internal SPD * matrix. */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G) { @@ -88,9 +87,8 @@ generalized_inverse(const EigMat& G) { the diagonal of the internal SPD * matrix. */ -template * = nullptr, - require_stan_scalar_t* = nullptr> +template * = nullptr, + require_stan_scalar_t* = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G, const Scal& a) { diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 778709456ae..e9982726d08 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -6,32 +6,32 @@ TEST(mathMixMatFun, generalized_inverse) { Eigen::MatrixXd t(0, 0); stan::test::expect_ad(f, t); -// stan::test::expect_ad_matvar(f, t); + // stan::test::expect_ad_matvar(f, t); Eigen::MatrixXd u(1, 1); u << 2; stan::test::expect_ad(f, u); -// stan::test::expect_ad_matvar(f, u); -// + // stan::test::expect_ad_matvar(f, u); + // Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; stan::test::expect_ad(f, v); -// stan::test::expect_ad_matvar(f, v); + // stan::test::expect_ad_matvar(f, v); v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; stan::test::expect_ad(f, v); -// stan::test::expect_ad_matvar(f, v); + // stan::test::expect_ad_matvar(f, v); -// issues around zero require looser tolerances for hessians -// stan::test::ad_tolerances tols; -// tols.hessian_hessian_ = 0.01; -// tols.hessian_fvar_hessian_ = 0.01; + // issues around zero require looser tolerances for hessians + // stan::test::ad_tolerances tols; + // tols.hessian_hessian_ = 0.01; + // tols.hessian_fvar_hessian_ = 0.01; -// Eigen::MatrixXd w(3, 4); -// w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; -// stan::test::expect_ad(tols, f, w); -// stan::test::expect_ad_matvar(tols, f, w); + // Eigen::MatrixXd w(3, 4); + // w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; + // stan::test::expect_ad(tols, f, w); + // stan::test::expect_ad_matvar(tols, f, w); -// even lower tolerance, again for cases around zero + // even lower tolerance, again for cases around zero stan::test::ad_tolerances tols2; tols2.hessian_hessian_ = 3.0; tols2.hessian_fvar_hessian_ = 3.0; @@ -39,22 +39,22 @@ TEST(mathMixMatFun, generalized_inverse) { Eigen::MatrixXd x(4, 4); x << 2, 3, 4, 5, 9, -1, 2, 2, 4, 3, 7, -1, 0, 1, 19, 112; stan::test::expect_ad(tols2, f, x); -// stan::test::expect_ad_matvar(tols2, f, x); + // stan::test::expect_ad_matvar(tols2, f, x); Eigen::MatrixXd z(2, 2); z << 1, 2, 5, std::numeric_limits::quiet_NaN(); EXPECT_NO_THROW(stan::math::generalized_inverse(z)); -// autodiff throws, so following fails (throw behavior must match to pass) -// stan::test::expect_ad(f, z); + // autodiff throws, so following fails (throw behavior must match to pass) + // stan::test::expect_ad(f, z); Eigen::MatrixXd a(2, 2); a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); stan::test::expect_ad(f, a); -// stan::test::expect_ad_matvar(f, a); + // stan::test::expect_ad_matvar(f, a); -// singular matrix, should use the -// alias to input small amount of jitter on the diagonal + // singular matrix, should use the + // alias to input small amount of jitter on the diagonal Eigen::MatrixXd m(3, 2); m << 1, 2, 2, 4, 1, 2; EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); From ec2be9984e83c82f08a2e7fb3b60f8eb74ba457d Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Thu, 3 Dec 2020 18:11:01 -0500 Subject: [PATCH 035/102] updated mix test --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index e9982726d08..1eb6bc09e91 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -22,13 +22,13 @@ TEST(mathMixMatFun, generalized_inverse) { // stan::test::expect_ad_matvar(f, v); // issues around zero require looser tolerances for hessians - // stan::test::ad_tolerances tols; - // tols.hessian_hessian_ = 0.01; - // tols.hessian_fvar_hessian_ = 0.01; + stan::test::ad_tolerances tols; + tols.hessian_hessian_ = 0.1; + tols.hessian_fvar_hessian_ = 0.1; - // Eigen::MatrixXd w(3, 4); - // w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; - // stan::test::expect_ad(tols, f, w); + Eigen::MatrixXd w(3, 4); + w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; + stan::test::expect_ad(tols, f, w); // stan::test::expect_ad_matvar(tols, f, w); // even lower tolerance, again for cases around zero @@ -59,3 +59,4 @@ TEST(mathMixMatFun, generalized_inverse) { m << 1, 2, 2, 4, 1, 2; EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); } + From 61a193d5c622781782360da3706228138e108356 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 3 Dec 2020 23:12:03 +0000 Subject: [PATCH 036/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 1eb6bc09e91..d98c2127b18 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -59,4 +59,3 @@ TEST(mathMixMatFun, generalized_inverse) { m << 1, 2, 2, 4, 1, 2; EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); } - From 99cd3e7ec7ae6f7c3e57a577806076dc9f226d8a Mon Sep 17 00:00:00 2001 From: spinkney Date: Fri, 4 Dec 2020 04:42:06 -0500 Subject: [PATCH 037/102] simplified diagonal add --- stan/math/prim/fun/generalized_inverse.hpp | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index add6ccd3c61..3f0ebd554cd 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -106,13 +106,11 @@ generalized_inverse(const EigMat& G, const Scal& a) { if (n < m) { Eigen::Matrix A = tcrossprod(G); - A.diagonal().array() - += Eigen::Array::Constant(n, a); + A.diagonal().array() += a; return transpose(mdivide_left_spd(A, G)); } else { Eigen::Matrix A = crossprod(G); - A.diagonal().array() - += Eigen::Array::Constant(m, a); + A.diagonal().array() += a; return transpose(mdivide_right_spd(G, A)); } } From c2120e01aac9534fc07298615600a70347520c50 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 05:42:18 -0500 Subject: [PATCH 038/102] using object reference --- stan/math/prim/fun/generalized_inverse.hpp | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 3f0ebd554cd..aadb4b9cb72 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -56,10 +56,12 @@ generalized_inverse(const EigMat& G) { return inverse(G); } + const auto& G_ref = to_ref(G); + if (n < m) { - return transpose(mdivide_left_spd(tcrossprod(G), G)); + return transpose(mdivide_left_spd(tcrossprod(G_ref), G_ref)); } else { - return transpose(mdivide_right_spd(G, crossprod(G))); + return transpose(mdivide_right_spd(G_ref, crossprod(G_ref))); } } @@ -104,14 +106,16 @@ generalized_inverse(const EigMat& G, const Scal& a) { return inverse(G); } + const auto& G_ref = to_ref(G); + if (n < m) { - Eigen::Matrix A = tcrossprod(G); + Eigen::Matrix A = tcrossprod(G_ref); A.diagonal().array() += a; - return transpose(mdivide_left_spd(A, G)); + return transpose(mdivide_left_spd(A, G_ref)); } else { - Eigen::Matrix A = crossprod(G); + Eigen::Matrix A = crossprod(G_ref); A.diagonal().array() += a; - return transpose(mdivide_right_spd(G, A)); + return transpose(mdivide_right_spd(G_ref, A)); } } From 1b62d5ea43cc6db261aabaab78a8f03dc609a349 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 6 Dec 2020 10:43:10 +0000 Subject: [PATCH 039/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index aadb4b9cb72..330e60f005a 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -109,7 +109,8 @@ generalized_inverse(const EigMat& G, const Scal& a) { const auto& G_ref = to_ref(G); if (n < m) { - Eigen::Matrix A = tcrossprod(G_ref); + Eigen::Matrix A + = tcrossprod(G_ref); A.diagonal().array() += a; return transpose(mdivide_left_spd(A, G_ref)); } else { From 90b95d85f00c5702887a65cf034dcdb5a39d40cb Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 06:30:00 -0500 Subject: [PATCH 040/102] add reverse mode differentiation --- stan/math/rev/fun/generalized_inverse.hpp | 124 ++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 stan/math/rev/fun/generalized_inverse.hpp diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp new file mode 100644 index 00000000000..c0ba88bdd01 --- /dev/null +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -0,0 +1,124 @@ +#ifndef STAN_MATH_REV_FUN_AENERALIZED_INVERSE_HPP +#define STAN_MATH_REV_FUN_AENERALIZED_INVERSE_HPP + +auto f = [](const auto& a) { + return stan::math::generalized_inverse(a); +}; +Eigen::MatrixXd A = Eigen::MatrixXd::Random(10, 10); +test_ad(f, A); +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Reverse mode differentiation algorithm reference: + * + * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate. + * Author(s): G. H. Golub and V. Pereyra. + * Source: SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 + * + * Equation 4.12 in the paper + */ +template * = nullptr> +inline auto generalized_inverse(const EigMat& A) { + using value_t = value_type_t; + if (A.size() == 0) { + return {}; + } + + const auto n = A.rows(); + const auto m = A.cols(); + + if (A.rows() == A.cols()) { + return inverse(A); + } + if (n < m) { + arena_t> A_arena(A); + arena_t inv_A(transpose(mdivide_left_spd(tcrossprod(A_arena.val()), A_arena.val()))); + reverse_pass_callback([A_arena, inv_A]() mutable { + A_arena.adj() += -inv_A * A.adj() * inv_A + + tcrossprod(inv_A) * A.adj().transpose() * (1 - A * inv_A.transpose()) + + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + }); + return ret; + } else { + arena_t> A_arena(A); + arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), crossprod(A_arena.val()))); + reverse_pass_callback([A_arena, inv_A]() mutable { + A_arena.adj() += -inv_A * A.adj() * inv_A + + tcrossprod(inv_A) * A.adj().transpose() * (1 - A * inv_A.tranpose()) + + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + }); + return ret; + } +} + +/** + * Reverse mode differentiation algorithm reference: + * + * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate. + * Author(s): G. H. Golub and V. Pereyra. + * Source: SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 + * + * Equation 4.12 in the paper + */ +template * = nullptr, + require_stan_scalar_t* = nullptr> +inline Eigen::Matrix, EigMat::RowsAtCompileTime, + EigMat::ColsAtCompileTime> +generalized_inverse(const EigMat& A, const Scal& a) { + using value_t = value_type_t; + if (A.size() == 0) { + return {}; + } + + const auto n = A.rows(); + const auto m = A.cols(); + + if (A.rows() == A.cols()) { + return inverse(A); + } + + if (n < m) { + arena_t> A_arena(A); + arena_t A_spd = tcrossprod(A_arena.val()); + A_spd.diagonal().array() += a; + + arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); + reverse_pass_callback([A_arena, inv_A]() mutable { + A_arena.adj() += -inv_A * A.adj() * inv_A + + tcrossprod(inv_A) * A.adj().transpose() * (1 - A * inv_A.transpose()) + + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + }); + return ret; + } else { + arena_t> A_arena(A); + arena_t A_spd = crossprod(A_arena.val()); + A_spd.diagonal().array() += a; + + arena_t> A_arena(A); + arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), A_spd)); + reverse_pass_callback([A_arena, inv_A]() mutable { + A_arena.adj() += -inv_A * A.adj() * inv_A + + tcrossprod(inv_A) * A.adj().transpose() * (1 - A * inv_A.tranpose()) + + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + }); + return ret; + } +} + +} // namespace math +} // namespace stan + +#endif \ No newline at end of file From 53544c0e35346ea0393d0c192c76837a04c374a4 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 6 Dec 2020 11:30:53 +0000 Subject: [PATCH 041/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 51 +++++++++++++---------- 1 file changed, 29 insertions(+), 22 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index c0ba88bdd01..af90f6b73c6 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -1,9 +1,7 @@ #ifndef STAN_MATH_REV_FUN_AENERALIZED_INVERSE_HPP #define STAN_MATH_REV_FUN_AENERALIZED_INVERSE_HPP -auto f = [](const auto& a) { - return stan::math::generalized_inverse(a); -}; +auto f = [](const auto& a) { return stan::math::generalized_inverse(a); }; Eigen::MatrixXd A = Eigen::MatrixXd::Random(10, 10); test_ad(f, A); #include @@ -24,9 +22,9 @@ namespace math { /** * Reverse mode differentiation algorithm reference: * - * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate. - * Author(s): G. H. Golub and V. Pereyra. - * Source: SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 + * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems + * Whose Variables Separate. Author(s): G. H. Golub and V. Pereyra. Source: SIAM + * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 * * Equation 4.12 in the paper */ @@ -45,20 +43,25 @@ inline auto generalized_inverse(const EigMat& A) { } if (n < m) { arena_t> A_arena(A); - arena_t inv_A(transpose(mdivide_left_spd(tcrossprod(A_arena.val()), A_arena.val()))); + arena_t inv_A( + transpose(mdivide_left_spd(tcrossprod(A_arena.val()), A_arena.val()))); reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A + - tcrossprod(inv_A) * A.adj().transpose() * (1 - A * inv_A.transpose()) + - (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + A_arena.adj() + += -inv_A * A.adj() * inv_A + + tcrossprod(inv_A) * A.adj().transpose() + * (1 - A * inv_A.transpose()) + + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } else { arena_t> A_arena(A); arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), crossprod(A_arena.val()))); reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A + - tcrossprod(inv_A) * A.adj().transpose() * (1 - A * inv_A.tranpose()) + - (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + A_arena.adj() + += -inv_A * A.adj() * inv_A + + tcrossprod(inv_A) * A.adj().transpose() + * (1 - A * inv_A.tranpose()) + + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } @@ -67,9 +70,9 @@ inline auto generalized_inverse(const EigMat& A) { /** * Reverse mode differentiation algorithm reference: * - * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems Whose Variables Separate. - * Author(s): G. H. Golub and V. Pereyra. - * Source: SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 + * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems + * Whose Variables Separate. Author(s): G. H. Golub and V. Pereyra. Source: SIAM + * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 * * Equation 4.12 in the paper */ @@ -97,9 +100,11 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A + - tcrossprod(inv_A) * A.adj().transpose() * (1 - A * inv_A.transpose()) + - (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + A_arena.adj() + += -inv_A * A.adj() * inv_A + + tcrossprod(inv_A) * A.adj().transpose() + * (1 - A * inv_A.transpose()) + + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } else { @@ -110,9 +115,11 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t> A_arena(A); arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), A_spd)); reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A + - tcrossprod(inv_A) * A.adj().transpose() * (1 - A * inv_A.tranpose()) + - (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + A_arena.adj() + += -inv_A * A.adj() * inv_A + + tcrossprod(inv_A) * A.adj().transpose() + * (1 - A * inv_A.tranpose()) + + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } From 30fdbc2fd5b20afbc3c983b77c3e5ba0515d5cc3 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 06:33:44 -0500 Subject: [PATCH 042/102] fixed typo --- stan/math/rev/fun/generalized_inverse.hpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index af90f6b73c6..e09f7be1c76 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -1,9 +1,6 @@ -#ifndef STAN_MATH_REV_FUN_AENERALIZED_INVERSE_HPP -#define STAN_MATH_REV_FUN_AENERALIZED_INVERSE_HPP +#ifndef STAN_MATH_REV_FUN_GENERALIZED_INVERSE_HPP +#define STAN_MATH_REV_FUN_GENERALIZED_INVERSE_HPP -auto f = [](const auto& a) { return stan::math::generalized_inverse(a); }; -Eigen::MatrixXd A = Eigen::MatrixXd::Random(10, 10); -test_ad(f, A); #include #include #include From d4f400f943c1aa9719414273acd9024fcb8559bd Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 06:39:11 -0500 Subject: [PATCH 043/102] add whitespace --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index e09f7be1c76..d63d630b3e9 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -125,4 +125,4 @@ generalized_inverse(const EigMat& A, const Scal& a) { } // namespace math } // namespace stan -#endif \ No newline at end of file +#endif From cdebea9391c5dfae3e8c4b99c31cd1d30451e1e8 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 06:48:31 -0500 Subject: [PATCH 044/102] make lines <= 80 char --- stan/math/rev/fun/generalized_inverse.hpp | 36 ++++++++++------------- 1 file changed, 16 insertions(+), 20 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index d63d630b3e9..c23ac77f5a9 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -43,22 +43,20 @@ inline auto generalized_inverse(const EigMat& A) { arena_t inv_A( transpose(mdivide_left_spd(tcrossprod(A_arena.val()), A_arena.val()))); reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() - += -inv_A * A.adj() * inv_A - + tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()) - + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + A_arena.adj() += -inv_A * A.adj() * inv_A; + A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() + * (1 - A * inv_A.transpose()); + A_arena.adj() += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } else { arena_t> A_arena(A); arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), crossprod(A_arena.val()))); reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() - += -inv_A * A.adj() * inv_A - + tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.tranpose()) - + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + A_arena.adj() += -inv_A * A.adj() * inv_A; + A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() + * (1 - A * inv_A.transpose()); + A_arena.adj() += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } @@ -97,11 +95,10 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() - += -inv_A * A.adj() * inv_A - + tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()) - + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + A_arena.adj() += -inv_A * A.adj() * inv_A; + A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() + * (1 - A * inv_A.transpose()); + A_arena.adj() += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } else { @@ -112,11 +109,10 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t> A_arena(A); arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), A_spd)); reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() - += -inv_A * A.adj() * inv_A - + tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.tranpose()) - + (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + A_arena.adj() += -inv_A * A.adj() * inv_A; + A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() + * (1 - A * inv_A.transpose()); + A_arena.adj() += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } From e9f25e8e13fd26c4be027df7af15eace84f49d2d Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 6 Dec 2020 11:49:14 +0000 Subject: [PATCH 045/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index c23ac77f5a9..a9dd7be144b 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -45,8 +45,9 @@ inline auto generalized_inverse(const EigMat& A) { reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()); - A_arena.adj() += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + * (1 - A * inv_A.transpose()); + A_arena.adj() + += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } else { @@ -55,8 +56,9 @@ inline auto generalized_inverse(const EigMat& A) { reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()); - A_arena.adj() += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + * (1 - A * inv_A.transpose()); + A_arena.adj() + += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } @@ -97,8 +99,9 @@ generalized_inverse(const EigMat& A, const Scal& a) { reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()); - A_arena.adj() += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + * (1 - A * inv_A.transpose()); + A_arena.adj() + += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } else { @@ -111,8 +114,9 @@ generalized_inverse(const EigMat& A, const Scal& a) { reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()); - A_arena.adj() += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + * (1 - A * inv_A.transpose()); + A_arena.adj() + += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); }); return ret; } From 2467a876eafcd38b4de857350c5378ab351e9afc Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 06:56:27 -0500 Subject: [PATCH 046/102] one more try at <= 80 chars --- stan/math/rev/fun/generalized_inverse.hpp | 28 +++++++++++++---------- 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index a9dd7be144b..5bb9f47abe4 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -42,23 +42,25 @@ inline auto generalized_inverse(const EigMat& A) { arena_t> A_arena(A); arena_t inv_A( transpose(mdivide_left_spd(tcrossprod(A_arena.val()), A_arena.val()))); + arena_t aP = (1 - A * inv_A.transpose()); + arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()); - A_arena.adj() - += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + * aP; + A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); }); return ret; } else { arena_t> A_arena(A); arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), crossprod(A_arena.val()))); + arena_t aP = (1 - A * inv_A.transpose()); + arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()); - A_arena.adj() - += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + * aP; + A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); }); return ret; } @@ -96,12 +98,13 @@ generalized_inverse(const EigMat& A, const Scal& a) { A_spd.diagonal().array() += a; arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); + arena_t aP = (1 - A * inv_A.transpose()); + arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()); - A_arena.adj() - += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + * aP; + A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); }); return ret; } else { @@ -111,12 +114,13 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t> A_arena(A); arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), A_spd)); + arena_t aP = (1 - A * inv_A.transpose()); + arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * (1 - A * inv_A.transpose()); - A_arena.adj() - += (1 - inv_A * A) * A.adj().transpose() * crossprod(inv_A()); + * aP; + A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); }); return ret; } From ee7f8ce261c64f75234855e41ed07353403546d5 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 6 Dec 2020 11:57:14 +0000 Subject: [PATCH 047/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 5bb9f47abe4..c71cc9fa138 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -46,8 +46,7 @@ inline auto generalized_inverse(const EigMat& A) { arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * aP; + A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); }); return ret; @@ -58,8 +57,7 @@ inline auto generalized_inverse(const EigMat& A) { arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * aP; + A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); }); return ret; @@ -102,8 +100,7 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * aP; + A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); }); return ret; @@ -118,8 +115,7 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() - * aP; + A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); }); return ret; From 351e0b08cb7091b2d5d99007c18360c34a00cf03 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 07:05:57 -0500 Subject: [PATCH 048/102] again --- stan/math/rev/fun/generalized_inverse.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index c71cc9fa138..cf8edf14da1 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -40,8 +40,8 @@ inline auto generalized_inverse(const EigMat& A) { } if (n < m) { arena_t> A_arena(A); - arena_t inv_A( - transpose(mdivide_left_spd(tcrossprod(A_arena.val()), A_arena.val()))); + arena_t A_spd = tcrossprod(A_arena.val()); + arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); arena_t aP = (1 - A * inv_A.transpose()); arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { @@ -52,8 +52,9 @@ inline auto generalized_inverse(const EigMat& A) { return ret; } else { arena_t> A_arena(A); - arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), crossprod(A_arena.val()))); - arena_t aP = (1 - A * inv_A.transpose()); + arena_t A_spd = crossprod(A_arena.val()); + arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), A_spd))); + arena_t aP = (1 - A * inv_A.transpose()); arena_t Pa = (1 - inv_A * A); reverse_pass_callback([A_arena, inv_A]() mutable { A_arena.adj() += -inv_A * A.adj() * inv_A; @@ -94,7 +95,7 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t> A_arena(A); arena_t A_spd = tcrossprod(A_arena.val()); A_spd.diagonal().array() += a; - + arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); arena_t aP = (1 - A * inv_A.transpose()); arena_t Pa = (1 - inv_A * A); From 9e6e692cb1f69058e1bf293800d9fa7660f3bc7f Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 6 Dec 2020 12:06:44 +0000 Subject: [PATCH 049/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index cf8edf14da1..50759bca595 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -95,7 +95,7 @@ generalized_inverse(const EigMat& A, const Scal& a) { arena_t> A_arena(A); arena_t A_spd = tcrossprod(A_arena.val()); A_spd.diagonal().array() += a; - + arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); arena_t aP = (1 - A * inv_A.transpose()); arena_t Pa = (1 - inv_A * A); From 113efb3f59e4f3c8045271950f6659cbd129ad97 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 07:18:06 -0500 Subject: [PATCH 050/102] keep name of input the same --- stan/math/rev/fun/generalized_inverse.hpp | 97 +++++++++++------------ 1 file changed, 48 insertions(+), 49 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index cf8edf14da1..4cf715987f5 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -26,40 +26,40 @@ namespace math { * Equation 4.12 in the paper */ template * = nullptr> -inline auto generalized_inverse(const EigMat& A) { +inline auto generalized_inverse(const EigMat& G) { using value_t = value_type_t; - if (A.size() == 0) { + if (G.size() == 0) { return {}; } - const auto n = A.rows(); - const auto m = A.cols(); + const auto n = G.rows(); + const auto m = G.cols(); - if (A.rows() == A.cols()) { - return inverse(A); + if (A.rows() == G.cols()) { + return inverse(G); } if (n < m) { - arena_t> A_arena(A); - arena_t A_spd = tcrossprod(A_arena.val()); - arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); - arena_t aP = (1 - A * inv_A.transpose()); - arena_t Pa = (1 - inv_A * A); - reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; - A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); + arena_t> G_arena(G); + arena_t A_spd = tcrossprod(G_arena.val()); + arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); + arena_t aP = (1 - G * inv_G.transpose()); + arena_t Pa = (1 - inv_G * G); + reverse_pass_callback([G_arena, inv_G]() mutable { + G_arena.adj() += -inv_G * G.adj() * inv_G; + G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; + G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); }); return ret; } else { - arena_t> A_arena(A); - arena_t A_spd = crossprod(A_arena.val()); - arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), A_spd))); - arena_t aP = (1 - A * inv_A.transpose()); - arena_t Pa = (1 - inv_A * A); - reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; - A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); + arena_t> G_arena(G); + arena_t A_spd = crossprod(G_arena.val()); + arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); + arena_t aP = (1 - G * inv_G.transpose()); + arena_t Pa = (1 - inv_G * G); + reverse_pass_callback([G_arena, inv_G]() mutable { + G_arena.adj() += -inv_G * G.adj() * inv_G; + G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; + G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); }); return ret; } @@ -78,46 +78,45 @@ template * = nullptr, require_stan_scalar_t* = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> -generalized_inverse(const EigMat& A, const Scal& a) { +generalized_inverse(const EigMat& G, const Scal& a) { using value_t = value_type_t; - if (A.size() == 0) { + if (G.size() == 0) { return {}; } - const auto n = A.rows(); - const auto m = A.cols(); + const auto n = G.rows(); + const auto m = G.cols(); - if (A.rows() == A.cols()) { - return inverse(A); + if (G.rows() == G.cols()) { + return inverse(G); } if (n < m) { - arena_t> A_arena(A); - arena_t A_spd = tcrossprod(A_arena.val()); + arena_t> G_arena(G); + arena_t A_spd = tcrossprod(G_arena.val()); A_spd.diagonal().array() += a; - arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); - arena_t aP = (1 - A * inv_A.transpose()); - arena_t Pa = (1 - inv_A * A); - reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; - A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); + arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); + arena_t aP = (1 - G * inv_G.transpose()); + arena_t Pa = (1 - inv_G * G); + reverse_pass_callback([G_arena, inv_G]() mutable { + G_arena.adj() += -inv_G * G.adj() * inv_G; + G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; + G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); }); return ret; } else { - arena_t> A_arena(A); - arena_t A_spd = crossprod(A_arena.val()); + arena_t> G_arena(G); + arena_t A_spd = crossprod(G_arena.val()); A_spd.diagonal().array() += a; - arena_t> A_arena(A); - arena_t inv_A(transpose(mdivide_right_spd(A_arena.val(), A_spd)); - arena_t aP = (1 - A * inv_A.transpose()); - arena_t Pa = (1 - inv_A * A); - reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; - A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); + arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd)); + arena_t aP = (1 - G * inv_G.transpose()); + arena_t Pa = (1 - inv_G * G); + reverse_pass_callback([G_arena, inv_G]() mutable { + G_arena.adj() += -inv_G * G.adj() * inv_G; + G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; + G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); }); return ret; } From 29dabbde7f5557c74e36d6bd7a07693d402fe9dd Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 6 Dec 2020 12:19:53 +0000 Subject: [PATCH 051/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index ef38c096839..46d126dae1f 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -96,11 +96,11 @@ generalized_inverse(const EigMat& G, const Scal& a) { arena_t A_spd = tcrossprod(G_arena.val()); A_spd.diagonal().array() += a; <<<<<<< HEAD - + arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); arena_t aP = (1 - G * inv_G.transpose()); arena_t Pa = (1 - inv_G * G); - reverse_pass_callback([G_arena, inv_G]() mutable { + reverse_pass_callback([G_arena, inv_G]() mutable { G_arena.adj() += -inv_G * G.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); From bf572d6e869e9c94dd0c60db4066505abc0e7220 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 07:21:23 -0500 Subject: [PATCH 052/102] fixed merge conflict --- stan/math/rev/fun/generalized_inverse.hpp | 11 ----------- 1 file changed, 11 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 46d126dae1f..674822ccea4 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -95,7 +95,6 @@ generalized_inverse(const EigMat& G, const Scal& a) { arena_t> G_arena(G); arena_t A_spd = tcrossprod(G_arena.val()); A_spd.diagonal().array() += a; -<<<<<<< HEAD arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); arena_t aP = (1 - G * inv_G.transpose()); @@ -104,16 +103,6 @@ generalized_inverse(const EigMat& G, const Scal& a) { G_arena.adj() += -inv_G * G.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); -======= - - arena_t inv_A(transpose(mdivide_left_spd(A_spd, A_arena.val()))); - arena_t aP = (1 - A * inv_A.transpose()); - arena_t Pa = (1 - inv_A * A); - reverse_pass_callback([A_arena, inv_A]() mutable { - A_arena.adj() += -inv_A * A.adj() * inv_A; - A_arena.adj() += tcrossprod(inv_A) * A.adj().transpose() * aP; - A_arena.adj() += Pa * A.adj().transpose() * crossprod(inv_A()); ->>>>>>> 9e6e692cb1f69058e1bf293800d9fa7660f3bc7f }); return ret; } else { From e1ba5c26551b98acc7ed717f9d06594990826655 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 10:26:55 -0500 Subject: [PATCH 053/102] fixed header --- stan/math/rev/fun/generalized_inverse.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 674822ccea4..87e636088c7 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -5,10 +5,10 @@ #include #include #include -#include +#include #include #include -#include +#include #include #include #include From 454a50f6bf2be5f9f241adb18c9beb9caa9d699b Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 10:41:09 -0500 Subject: [PATCH 054/102] more header fix --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 87e636088c7..dd413db9c63 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -7,7 +7,7 @@ #include #include #include -#include +#include #include #include #include From 4af6787eb079a0f6ecfc77c067611f2135c0d921 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 10:56:15 -0500 Subject: [PATCH 055/102] :( --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index dd413db9c63..cc130142f7a 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -10,7 +10,7 @@ #include #include #include -#include +#include #include namespace stan { From 9f713350adf463b5bf0a504b1ea3cad1ba05f71e Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 11:14:14 -0500 Subject: [PATCH 056/102] cleanup --- stan/math/rev/fun/generalized_inverse.hpp | 56 +++++++++++------------ 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index cc130142f7a..8fb2ca9458f 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -40,28 +40,28 @@ inline auto generalized_inverse(const EigMat& G) { } if (n < m) { arena_t> G_arena(G); - arena_t A_spd = tcrossprod(G_arena.val()); + auto A_spd = tcrossprod(G_arena.val()); arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); - arena_t aP = (1 - G * inv_G.transpose()); - arena_t Pa = (1 - inv_G * G); + auto aP = (1 - G_arena.val() * inv_G.transpose()); + auto Pa = (1 - inv_G * G_arena.val()); reverse_pass_callback([G_arena, inv_G]() mutable { - G_arena.adj() += -inv_G * G.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; - G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); + G_arena.adj() += -inv_G * G_arena.adj() * inv_G; + G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; + G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); }); - return ret; + return inv_G; } else { arena_t> G_arena(G); - arena_t A_spd = crossprod(G_arena.val()); + auto A_spd = crossprod(G_arena.val()); arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); - arena_t aP = (1 - G * inv_G.transpose()); - arena_t Pa = (1 - inv_G * G); + auto aP = (1 - G_arena.val() * inv_G.transpose()); + auto Pa = (1 - inv_G * G_arena.val()); reverse_pass_callback([G_arena, inv_G]() mutable { - G_arena.adj() += -inv_G * G.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; - G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); + G_arena.adj() += -inv_G * G_arena.adj() * inv_G; + G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; + G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); }); - return ret; + return inv_G; } } @@ -93,32 +93,32 @@ generalized_inverse(const EigMat& G, const Scal& a) { if (n < m) { arena_t> G_arena(G); - arena_t A_spd = tcrossprod(G_arena.val()); + auto A_spd = tcrossprod(G_arena.val()); A_spd.diagonal().array() += a; arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); - arena_t aP = (1 - G * inv_G.transpose()); - arena_t Pa = (1 - inv_G * G); + auto aP = (1 - G_arena.val() * inv_G.transpose()); + auto Pa = (1 - inv_G * G_arena.val()); reverse_pass_callback([G_arena, inv_G]() mutable { - G_arena.adj() += -inv_G * G.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; - G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); + G_arena.adj() += -inv_G * G_arena.adj() * inv_G; + G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; + G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G); }); - return ret; + return inv_G; } else { arena_t> G_arena(G); - arena_t A_spd = crossprod(G_arena.val()); + auto A_spd = crossprod(G_arena.val()); A_spd.diagonal().array() += a; arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd)); - arena_t aP = (1 - G * inv_G.transpose()); - arena_t Pa = (1 - inv_G * G); + auto aP = (1 - G * inv_G.transpose()); + auto Pa = (1 - inv_G * G_arena.val()); reverse_pass_callback([G_arena, inv_G]() mutable { - G_arena.adj() += -inv_G * G.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G.adj().transpose() * aP; - G_arena.adj() += Pa * G.adj().transpose() * crossprod(inv_G()); + G_arena.adj() += -inv_G * G_arena.adj() * inv_G; + G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; + G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); }); - return ret; + return inv_G; } } From 63784101cd277e1a9ded1f4597bac050cf2198f7 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 11:27:01 -0500 Subject: [PATCH 057/102] more cleanup --- stan/math/rev/fun/generalized_inverse.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 8fb2ca9458f..d5e2bf1181f 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -29,13 +29,13 @@ template * = nullptr> inline auto generalized_inverse(const EigMat& G) { using value_t = value_type_t; if (G.size() == 0) { - return {}; + return G; } const auto n = G.rows(); const auto m = G.cols(); - if (A.rows() == G.cols()) { + if (G.rows() == G.cols()) { return inverse(G); } if (n < m) { @@ -81,7 +81,7 @@ inline Eigen::Matrix, EigMat::RowsAtCompileTime, generalized_inverse(const EigMat& G, const Scal& a) { using value_t = value_type_t; if (G.size() == 0) { - return {}; + return G; } const auto n = G.rows(); @@ -110,7 +110,7 @@ generalized_inverse(const EigMat& G, const Scal& a) { auto A_spd = crossprod(G_arena.val()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd)); + arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); auto aP = (1 - G * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); reverse_pass_callback([G_arena, inv_G]() mutable { From 5866d8889d9dd524b207729da8014de24de90351 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 6 Dec 2020 16:28:03 +0000 Subject: [PATCH 058/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index d5e2bf1181f..ce546b0133e 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -113,7 +113,7 @@ generalized_inverse(const EigMat& G, const Scal& a) { arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); auto aP = (1 - G * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G]() mutable { + reverse_pass_callback([G_arena, inv_G]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); From 27c31ed1718bf80f7c24f17f5d8dfe2d73af95e7 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 12:54:52 -0500 Subject: [PATCH 059/102] fix reverse_pass_callback --- stan/math/rev/fun/generalized_inverse.hpp | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index ce546b0133e..6cab20fc936 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -44,7 +44,7 @@ inline auto generalized_inverse(const EigMat& G) { arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); auto aP = (1 - G_arena.val() * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G]() mutable { + reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); @@ -56,7 +56,7 @@ inline auto generalized_inverse(const EigMat& G) { arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); auto aP = (1 - G_arena.val() * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G]() mutable { + reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); @@ -99,7 +99,7 @@ generalized_inverse(const EigMat& G, const Scal& a) { arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); auto aP = (1 - G_arena.val() * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G]() mutable { + reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G); @@ -113,7 +113,11 @@ generalized_inverse(const EigMat& G, const Scal& a) { arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); auto aP = (1 - G * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); +<<<<<<< Updated upstream reverse_pass_callback([G_arena, inv_G]() mutable { +======= + reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { +>>>>>>> Stashed changes G_arena.adj() += -inv_G * G_arena.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); From e2145be175b76cf5dfb7ba0d93340215f345ed20 Mon Sep 17 00:00:00 2001 From: spinkney Date: Sun, 6 Dec 2020 12:55:49 -0500 Subject: [PATCH 060/102] fix merge --- stan/math/rev/fun/generalized_inverse.hpp | 4 ---- 1 file changed, 4 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 6cab20fc936..f5994cc0116 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -113,11 +113,7 @@ generalized_inverse(const EigMat& G, const Scal& a) { arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); auto aP = (1 - G * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); -<<<<<<< Updated upstream - reverse_pass_callback([G_arena, inv_G]() mutable { -======= reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { ->>>>>>> Stashed changes G_arena.adj() += -inv_G * G_arena.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); From cc21e248eaac4d4c585133171b213cc798ec931b Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 6 Dec 2020 17:57:02 +0000 Subject: [PATCH 061/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index f5994cc0116..c8c7d43598a 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -113,7 +113,7 @@ generalized_inverse(const EigMat& G, const Scal& a) { arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); auto aP = (1 - G * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { + reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); From e9365fa5989365ab0142c91dea1847391dbba6b7 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Mon, 7 Dec 2020 10:06:12 -0500 Subject: [PATCH 062/102] Added tests for jittered function --- stan/math/prim/fun/generalized_inverse.hpp | 24 +++---- stan/math/rev/fun/generalized_inverse.hpp | 26 +++---- .../math/mix/fun/generalized_inverse_test.cpp | 68 ++++++++++--------- 3 files changed, 58 insertions(+), 60 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 330e60f005a..4d30dbb9803 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -4,6 +4,7 @@ #include #include #include +#include #include #include #include @@ -40,21 +41,21 @@ namespace math { diagonal of the internal SPD * matrix. */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G) { using value_t = value_type_t; - if (G.size() == 0) { + + if (G.size() == 0) return {}; - } const auto n = G.rows(); const auto m = G.cols(); - if (G.rows() == G.cols()) { + if (G.rows() == G.cols()) return inverse(G); - } const auto& G_ref = to_ref(G); @@ -89,22 +90,21 @@ generalized_inverse(const EigMat& G) { the diagonal of the internal SPD * matrix. */ -template * = nullptr, - require_stan_scalar_t* = nullptr> +template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> -generalized_inverse(const EigMat& G, const Scal& a) { +generalized_inverse(const EigMat& G, const double a) { using value_t = value_type_t; - if (G.size() == 0) { + + if (G.size() == 0) return {}; - } const auto n = G.rows(); const auto m = G.cols(); - if (G.rows() == G.cols()) { + if (G.rows() == G.cols()) return inverse(G); - } const auto& G_ref = to_ref(G); diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index c8c7d43598a..23ebdcda50d 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -25,19 +25,19 @@ namespace math { * * Equation 4.12 in the paper */ -template * = nullptr> +template * = nullptr> inline auto generalized_inverse(const EigMat& G) { using value_t = value_type_t; - if (G.size() == 0) { - return G; - } + if (G.size() == 0) + return G; + const auto n = G.rows(); const auto m = G.cols(); - if (G.rows() == G.cols()) { + if (G.rows() == G.cols()) return inverse(G); - } + if (n < m) { arena_t> G_arena(G); auto A_spd = tcrossprod(G_arena.val()); @@ -74,22 +74,18 @@ inline auto generalized_inverse(const EigMat& G) { * * Equation 4.12 in the paper */ -template * = nullptr, - require_stan_scalar_t* = nullptr> -inline Eigen::Matrix, EigMat::RowsAtCompileTime, - EigMat::ColsAtCompileTime> -generalized_inverse(const EigMat& G, const Scal& a) { +template * = nullptr> +inline auto generalized_inverse(const EigMat& G, const double a) { using value_t = value_type_t; - if (G.size() == 0) { + + if (G.size() == 0) return G; - } const auto n = G.rows(); const auto m = G.cols(); - if (G.rows() == G.cols()) { + if (G.rows() == G.cols()) return inverse(G); - } if (n < m) { arena_t> G_arena(G); diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index d98c2127b18..41d877ec458 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -1,61 +1,63 @@ #include #include - -TEST(mathMixMatFun, generalized_inverse) { - auto f = [](const auto& x) { return stan::math::generalized_inverse(x); }; - +#include + +TEST(mathMixMatFun, ad_tests) { + using stan::test::expect_ad; + + auto f = [](const auto& G) { return stan::math::generalized_inverse(G); }; + auto fjit = [](const auto& G) { return stan::math::generalized_inverse(G, 0.0); }; + auto fjit2 = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-4); }; + Eigen::MatrixXd t(0, 0); - stan::test::expect_ad(f, t); - // stan::test::expect_ad_matvar(f, t); + expect_ad(f, t); + expect_ad(fjit, t); Eigen::MatrixXd u(1, 1); u << 2; - stan::test::expect_ad(f, u); - // stan::test::expect_ad_matvar(f, u); - // + expect_ad(f, u); + expect_ad(fjit, u); + Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; - stan::test::expect_ad(f, v); - // stan::test::expect_ad_matvar(f, v); + expect_ad(f, v); + expect_ad(fjit, v); + v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; - stan::test::expect_ad(f, v); - // stan::test::expect_ad_matvar(f, v); + expect_ad(f, v); + expect_ad(fjit, v); - // issues around zero require looser tolerances for hessians + // issues around zero require looser tolerances for hessians stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 0.1; - tols.hessian_fvar_hessian_ = 0.1; + tols.hessian_hessian_ = 2.0; + tols.hessian_fvar_hessian_ = 2.0; Eigen::MatrixXd w(3, 4); w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; - stan::test::expect_ad(tols, f, w); - // stan::test::expect_ad_matvar(tols, f, w); - - // even lower tolerance, again for cases around zero - stan::test::ad_tolerances tols2; - tols2.hessian_hessian_ = 3.0; - tols2.hessian_fvar_hessian_ = 3.0; - - Eigen::MatrixXd x(4, 4); - x << 2, 3, 4, 5, 9, -1, 2, 2, 4, 3, 7, -1, 0, 1, 19, 112; - stan::test::expect_ad(tols2, f, x); - // stan::test::expect_ad_matvar(tols2, f, x); + expect_ad(tols, f, w); + expect_ad(tols, fjit, w); Eigen::MatrixXd z(2, 2); z << 1, 2, 5, std::numeric_limits::quiet_NaN(); EXPECT_NO_THROW(stan::math::generalized_inverse(z)); - + EXPECT_NO_THROW(stan::math::generalized_inverse(z, 0.0)); + // autodiff throws, so following fails (throw behavior must match to pass) - // stan::test::expect_ad(f, z); Eigen::MatrixXd a(2, 2); a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); - stan::test::expect_ad(f, a); - // stan::test::expect_ad_matvar(f, a); - + expect_ad(f, a); + expect_ad(fjit, a); + // singular matrix, should use the // alias to input small amount of jitter on the diagonal Eigen::MatrixXd m(3, 2); m << 1, 2, 2, 4, 1, 2; EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); + + // should work with jittered version + stan::test::ad_tolerances tols3; + tols3.hessian_hessian_ = 0.01; + tols3.hessian_fvar_hessian_ = 0.01; + expect_ad(tols3, fjit2, m); } From 6adb8ac2cfa1f9e203903e78770d5e538db2b8a1 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 7 Dec 2020 15:10:36 +0000 Subject: [PATCH 063/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 10 ++++------ stan/math/rev/fun/generalized_inverse.hpp | 14 ++++++------- .../math/mix/fun/generalized_inverse_test.cpp | 20 ++++++++++--------- 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 4d30dbb9803..3305c3b1978 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -41,14 +41,13 @@ namespace math { diagonal of the internal SPD * matrix. */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G) { using value_t = value_type_t; - if (G.size() == 0) + if (G.size() == 0) return {}; const auto n = G.rows(); @@ -90,14 +89,13 @@ generalized_inverse(const EigMat& G) { the diagonal of the internal SPD * matrix. */ -template * = nullptr> +template * = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G, const double a) { using value_t = value_type_t; - if (G.size() == 0) + if (G.size() == 0) return {}; const auto n = G.rows(); diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 23ebdcda50d..477a7cdbf3c 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -29,15 +29,15 @@ template * = nullptr> inline auto generalized_inverse(const EigMat& G) { using value_t = value_type_t; - if (G.size() == 0) + if (G.size() == 0) return G; - + const auto n = G.rows(); const auto m = G.cols(); - if (G.rows() == G.cols()) + if (G.rows() == G.cols()) return inverse(G); - + if (n < m) { arena_t> G_arena(G); auto A_spd = tcrossprod(G_arena.val()); @@ -77,14 +77,14 @@ inline auto generalized_inverse(const EigMat& G) { template * = nullptr> inline auto generalized_inverse(const EigMat& G, const double a) { using value_t = value_type_t; - - if (G.size() == 0) + + if (G.size() == 0) return G; const auto n = G.rows(); const auto m = G.cols(); - if (G.rows() == G.cols()) + if (G.rows() == G.cols()) return inverse(G); if (n < m) { diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 41d877ec458..a76e5cbe99f 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -4,11 +4,13 @@ TEST(mathMixMatFun, ad_tests) { using stan::test::expect_ad; - + auto f = [](const auto& G) { return stan::math::generalized_inverse(G); }; - auto fjit = [](const auto& G) { return stan::math::generalized_inverse(G, 0.0); }; - auto fjit2 = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-4); }; - + auto fjit + = [](const auto& G) { return stan::math::generalized_inverse(G, 0.0); }; + auto fjit2 + = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-4); }; + Eigen::MatrixXd t(0, 0); expect_ad(f, t); expect_ad(fjit, t); @@ -17,17 +19,17 @@ TEST(mathMixMatFun, ad_tests) { u << 2; expect_ad(f, u); expect_ad(fjit, u); - + Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); expect_ad(fjit, v); - + v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; expect_ad(f, v); expect_ad(fjit, v); - // issues around zero require looser tolerances for hessians + // issues around zero require looser tolerances for hessians stan::test::ad_tolerances tols; tols.hessian_hessian_ = 2.0; tols.hessian_fvar_hessian_ = 2.0; @@ -41,14 +43,14 @@ TEST(mathMixMatFun, ad_tests) { z << 1, 2, 5, std::numeric_limits::quiet_NaN(); EXPECT_NO_THROW(stan::math::generalized_inverse(z)); EXPECT_NO_THROW(stan::math::generalized_inverse(z, 0.0)); - + // autodiff throws, so following fails (throw behavior must match to pass) Eigen::MatrixXd a(2, 2); a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); expect_ad(f, a); expect_ad(fjit, a); - + // singular matrix, should use the // alias to input small amount of jitter on the diagonal Eigen::MatrixXd m(3, 2); From b1d7c4f2096d0c89dd6cfe911702bb8b3f6bea16 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Sun, 13 Dec 2020 07:16:44 -0500 Subject: [PATCH 064/102] Update stan/math/rev/fun/generalized_inverse.hpp change G to referenced object G_arena Co-authored-by: Steve Bronder --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 477a7cdbf3c..92aeac180bc 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -107,7 +107,7 @@ inline auto generalized_inverse(const EigMat& G, const double a) { A_spd.diagonal().array() += a; arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); - auto aP = (1 - G * inv_G.transpose()); + auto aP = (1 - G_arena * inv_G.transpose()); auto Pa = (1 - inv_G * G_arena.val()); reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; From f1efc6c35935152edc960b385f832c12825dead5 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Sun, 13 Dec 2020 07:18:03 -0500 Subject: [PATCH 065/102] Update stan/math/prim/fun/generalized_inverse.hpp Doc update Co-authored-by: Steve Bronder --- stan/math/prim/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 3305c3b1978..41b2ece7c74 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -21,7 +21,7 @@ namespace math { * Returns the Moore-Penrose generalized inverse of the specified matrix. * * The method is based on the Cholesky computation of the transform as specified - in + * in * *
    • Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse Matrices. From d2ef3ce424d7b5cb8a02bfcefc4a867a3095e3ce Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 06:12:27 -0500 Subject: [PATCH 066/102] fixed derivative --- stan/math/rev/fun/generalized_inverse.hpp | 95 ++++++++++++++++------- 1 file changed, 65 insertions(+), 30 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 92aeac180bc..485db227845 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -16,7 +16,7 @@ namespace stan { namespace math { -/** +/* * Reverse mode differentiation algorithm reference: * * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems @@ -24,8 +24,17 @@ namespace math { * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 * * Equation 4.12 in the paper + * + * See also + * http://mathoverflow.net/questions/25778/analytical-formula-for-numerical-derivative-of-the-matrix-pseudo-inverse + * + * Implementation is based on + * Title: Uncertainties Python Package + * Author: Eric O. LEBIGOT + * Date: 2020 + * Availability: https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py */ -template * = nullptr> +template * = nullptr> inline auto generalized_inverse(const EigMat& G) { using value_t = value_type_t; @@ -39,33 +48,43 @@ inline auto generalized_inverse(const EigMat& G) { return inverse(G); if (n < m) { - arena_t> G_arena(G); + arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val()); arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); - auto aP = (1 - G_arena.val() * inv_G.transpose()); - auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { + + auto PG = to_arena( -G_arena * inv_G ); + PG.diagonal().array() += 1.0; + + auto GP = to_arena( -inv_G * G_arena ); + GP.diagonal().array() += 1.0; + + reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; - G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); + G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * PG; + G_arena.adj() += GP * G_arena.adj().transpose() * crossprod(inv_G); }); return inv_G; } else { - arena_t> G_arena(G); + arena_t G_arena(G); auto A_spd = crossprod(G_arena.val()); arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); - auto aP = (1 - G_arena.val() * inv_G.transpose()); - auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { + + auto PG = to_arena( -G_arena * inv_G ); + PG.diagonal().array() += 1.0; + + auto GP = to_arena( -inv_G * G_arena ); + GP.diagonal().array() += 1.0; + + reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; - G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); + G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * PG; + G_arena.adj() += GP * G_arena.adj().transpose() * crossprod(inv_G); }); return inv_G; } } -/** +/* * Reverse mode differentiation algorithm reference: * * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems @@ -73,8 +92,16 @@ inline auto generalized_inverse(const EigMat& G) { * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 * * Equation 4.12 in the paper + * See also + * http://mathoverflow.net/questions/25778/analytical-formula-for-numerical-derivative-of-the-matrix-pseudo-inverse + * + * Implementation is based on + * Title: Uncertainties Python Package + * Author: Eric O. LEBIGOT + * Date: 2020 + * Availability: https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py */ -template * = nullptr> +template * = nullptr> inline auto generalized_inverse(const EigMat& G, const double a) { using value_t = value_type_t; @@ -88,31 +115,39 @@ inline auto generalized_inverse(const EigMat& G, const double a) { return inverse(G); if (n < m) { - arena_t> G_arena(G); + arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); - auto aP = (1 - G_arena.val() * inv_G.transpose()); - auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { + + auto PG = to_arena( -G_arena * inv_G ); + PG.diagonal().array() += 1.0; + + auto GP = to_arena( -inv_G * G_arena ); + GP.diagonal().array() += 1.0; + + reverse_pass_callback([G_arena, inv_G, PG, GP]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; - G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G); + G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * PG; + G_arena.adj() += GP * G_arena.adj().transpose() * crossprod(inv_G); }); return inv_G; } else { - arena_t> G_arena(G); + arena_t G_arena(G); auto A_spd = crossprod(G_arena.val()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); - auto aP = (1 - G_arena * inv_G.transpose()); - auto Pa = (1 - inv_G * G_arena.val()); - reverse_pass_callback([G_arena, inv_G, aP, Pa]() mutable { + + auto PG = to_arena( -G_arena * inv_G ); + PG.diagonal().array() += 1.0; + + auto GP = to_arena( -inv_G * G_arena ); + GP.diagonal().array() += 1.0; + + reverse_pass_callback([G_arena, inv_G, PG, GP]() mutable { G_arena.adj() += -inv_G * G_arena.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * aP; - G_arena.adj() += Pa * G_arena.adj().transpose() * crossprod(inv_G()); + G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * PG; + G_arena.adj() += GP * G_arena.adj().transpose() * crossprod(inv_G); }); return inv_G; } From 2f32b8d775705d41cd35d2228269c14506482f3d Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 06:18:43 -0500 Subject: [PATCH 067/102] fixed doc and include --- stan/math/prim/fun/generalized_inverse.hpp | 4 ++-- stan/math/rev/fun.hpp | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 41b2ece7c74..653eaf41f35 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -27,7 +27,7 @@ namespace math { Matrices. * arXiv 0804.4809
    * - * @tparam EigMat type of the matrix (must be derived from \c Eigen::MatrixBase) + * @tparam EigMat type of the matrix (must be derived from `Eigen::MatrixBase`) * * @param G specified matrix * @return Generalized inverse of the matrix (an empty matrix if the specified @@ -75,7 +75,7 @@ generalized_inverse(const EigMat& G) { Matrices. * arXiv 0804.4809
* - * @tparam EigMat type of the matrix (must be derived from \c Eigen::MatrixBase) + * @tparam EigMat type of the matrix (must be derived from `Eigen::MatrixBase`) * * @param G specified matrix * @param a real constant diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index d53279c4588..e7893d15e6b 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -59,6 +59,7 @@ #include #include #include +#include #include #include #include From ee29e910d4b1a3c9c66a93df33c5963309e93db0 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 14 Dec 2020 11:18:46 +0000 Subject: [PATCH 068/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 28 ++++++++++++----------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 485db227845..07eaabbf41e 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -32,7 +32,8 @@ namespace math { * Title: Uncertainties Python Package * Author: Eric O. LEBIGOT * Date: 2020 - * Availability: https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py + * Availability: + * https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py */ template * = nullptr> inline auto generalized_inverse(const EigMat& G) { @@ -51,11 +52,11 @@ inline auto generalized_inverse(const EigMat& G) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val()); arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); - - auto PG = to_arena( -G_arena * inv_G ); + + auto PG = to_arena(-G_arena * inv_G); PG.diagonal().array() += 1.0; - auto GP = to_arena( -inv_G * G_arena ); + auto GP = to_arena(-inv_G * G_arena); GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { @@ -68,11 +69,11 @@ inline auto generalized_inverse(const EigMat& G) { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val()); arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); - - auto PG = to_arena( -G_arena * inv_G ); + + auto PG = to_arena(-G_arena * inv_G); PG.diagonal().array() += 1.0; - auto GP = to_arena( -inv_G * G_arena ); + auto GP = to_arena(-inv_G * G_arena); GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { @@ -99,7 +100,8 @@ inline auto generalized_inverse(const EigMat& G) { * Title: Uncertainties Python Package * Author: Eric O. LEBIGOT * Date: 2020 - * Availability: https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py + * Availability: + * https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py */ template * = nullptr> inline auto generalized_inverse(const EigMat& G, const double a) { @@ -120,10 +122,10 @@ inline auto generalized_inverse(const EigMat& G, const double a) { A_spd.diagonal().array() += a; arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); - auto PG = to_arena( -G_arena * inv_G ); + auto PG = to_arena(-G_arena * inv_G); PG.diagonal().array() += 1.0; - auto GP = to_arena( -inv_G * G_arena ); + auto GP = to_arena(-inv_G * G_arena); GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, inv_G, PG, GP]() mutable { @@ -137,11 +139,11 @@ inline auto generalized_inverse(const EigMat& G, const double a) { auto A_spd = crossprod(G_arena.val()); A_spd.diagonal().array() += a; arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); - - auto PG = to_arena( -G_arena * inv_G ); + + auto PG = to_arena(-G_arena * inv_G); PG.diagonal().array() += 1.0; - auto GP = to_arena( -inv_G * G_arena ); + auto GP = to_arena(-inv_G * G_arena); GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, inv_G, PG, GP]() mutable { From 88a82627da7d38a133e5d6d69c2371bb9e62fb47 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 06:25:10 -0500 Subject: [PATCH 069/102] changed rev template from EigMat to VarMat --- stan/math/rev/fun/generalized_inverse.hpp | 24 +++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 07eaabbf41e..c51da8731d8 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -36,8 +36,8 @@ namespace math { * https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py */ template * = nullptr> -inline auto generalized_inverse(const EigMat& G) { - using value_t = value_type_t; +inline auto generalized_inverse(const VarMat& G) { + using value_t = value_type_t; if (G.size() == 0) return G; @@ -49,9 +49,9 @@ inline auto generalized_inverse(const EigMat& G) { return inverse(G); if (n < m) { - arena_t G_arena(G); + arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val()); - arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); + arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); auto PG = to_arena(-G_arena * inv_G); PG.diagonal().array() += 1.0; @@ -66,9 +66,9 @@ inline auto generalized_inverse(const EigMat& G) { }); return inv_G; } else { - arena_t G_arena(G); + arena_t G_arena(G); auto A_spd = crossprod(G_arena.val()); - arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); + arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); auto PG = to_arena(-G_arena * inv_G); PG.diagonal().array() += 1.0; @@ -104,8 +104,8 @@ inline auto generalized_inverse(const EigMat& G) { * https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py */ template * = nullptr> -inline auto generalized_inverse(const EigMat& G, const double a) { - using value_t = value_type_t; +inline auto generalized_inverse(const VarMat& G, const double a) { + using value_t = value_type_t; if (G.size() == 0) return G; @@ -117,10 +117,10 @@ inline auto generalized_inverse(const EigMat& G, const double a) { return inverse(G); if (n < m) { - arena_t G_arena(G); + arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); + arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); auto PG = to_arena(-G_arena * inv_G); PG.diagonal().array() += 1.0; @@ -135,10 +135,10 @@ inline auto generalized_inverse(const EigMat& G, const double a) { }); return inv_G; } else { - arena_t G_arena(G); + arena_t G_arena(G); auto A_spd = crossprod(G_arena.val()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); + arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); auto PG = to_arena(-G_arena * inv_G); PG.diagonal().array() += 1.0; From 8d4e39cb0787a09f4dd8b2a194f1fa9027647a1f Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:10:03 -0500 Subject: [PATCH 070/102] incorporating Steve's rev compile fixes --- stan/math/prim/fun/generalized_inverse.hpp | 6 +- stan/math/rev/fun/generalized_inverse.hpp | 71 +++++++++---------- .../math/mix/fun/generalized_inverse_test.cpp | 6 +- 3 files changed, 42 insertions(+), 41 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 653eaf41f35..6558f298427 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -41,7 +41,8 @@ namespace math { diagonal of the internal SPD * matrix. */ -template * = nullptr> +template * = nullptr, + require_not_vt_var* = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G) { @@ -89,7 +90,8 @@ generalized_inverse(const EigMat& G) { the diagonal of the internal SPD * matrix. */ -template * = nullptr> +template * = nullptr, + require_not_vt_var* = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G, const double a) { diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index c51da8731d8..087775249e0 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -38,50 +38,48 @@ namespace math { template * = nullptr> inline auto generalized_inverse(const VarMat& G) { using value_t = value_type_t; + using ret_type = promote_var_matrix_t; if (G.size() == 0) - return G; - - const auto n = G.rows(); - const auto m = G.cols(); + return ret_type(G); if (G.rows() == G.cols()) return inverse(G); - if (n < m) { + if (G.rows() < G.cols()) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val()); - arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); + arena_t inv_G(transpose(mdivide_left_spd(A_spd.val_op(), G_arena.val_op()))); - auto PG = to_arena(-G_arena * inv_G); + auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; - auto GP = to_arena(-inv_G * G_arena); + auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { - G_arena.adj() += -inv_G * G_arena.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * PG; - G_arena.adj() += GP * G_arena.adj().transpose() * crossprod(inv_G); + G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); + G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); + G_arena.adj() += GP.val_op() * G_arena.adj().transpose() * crossprod(inv_G.val_op()); }); - return inv_G; + return ret_type(inv_G); } else { arena_t G_arena(G); - auto A_spd = crossprod(G_arena.val()); - arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); + auto A_spd = crossprod(G_arena.val_op()); + arena_t inv_G(transpose(mdivide_right_spd(G_arena.val_op(), A_spd.val_op()))); - auto PG = to_arena(-G_arena * inv_G); + auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; - auto GP = to_arena(-inv_G * G_arena); + auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { - G_arena.adj() += -inv_G * G_arena.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * PG; - G_arena.adj() += GP * G_arena.adj().transpose() * crossprod(inv_G); + G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); + G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); + G_arena.adj() += GP.val_op() * G_arena.adj().transpose() * crossprod(inv_G.val_op()); }); - return inv_G; + return ret_type(inv_G); } } @@ -110,48 +108,45 @@ inline auto generalized_inverse(const VarMat& G, const double a) { if (G.size() == 0) return G; - const auto n = G.rows(); - const auto m = G.cols(); - if (G.rows() == G.cols()) return inverse(G); - if (n < m) { + if (G.rows() < G.cols()) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val()); A_spd.diagonal().array() += a; arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); - auto PG = to_arena(-G_arena * inv_G); + auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; - auto GP = to_arena(-inv_G * G_arena); + auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); GP.diagonal().array() += 1.0; - reverse_pass_callback([G_arena, inv_G, PG, GP]() mutable { - G_arena.adj() += -inv_G * G_arena.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * PG; - G_arena.adj() += GP * G_arena.adj().transpose() * crossprod(inv_G); + reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { + G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); + G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); + G_arena.adj() += GP.val_op() * G_arena.adj().transpose() * crossprod(inv_G.val_op()); }); - return inv_G; + return ret_type(inv_G); } else { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val()); A_spd.diagonal().array() += a; arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); - auto PG = to_arena(-G_arena * inv_G); + auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; - auto GP = to_arena(-inv_G * G_arena); + auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); GP.diagonal().array() += 1.0; - reverse_pass_callback([G_arena, inv_G, PG, GP]() mutable { - G_arena.adj() += -inv_G * G_arena.adj() * inv_G; - G_arena.adj() += tcrossprod(inv_G) * G_arena.adj().transpose() * PG; - G_arena.adj() += GP * G_arena.adj().transpose() * crossprod(inv_G); + reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { + G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); + G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); + G_arena.adj() += GP.val_op() * G_arena.adj().transpose() * crossprod(inv_G.val_op()); }); - return inv_G; + return ret_type(inv_G); } } diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index a76e5cbe99f..e27b5ef2dd3 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -14,20 +14,22 @@ TEST(mathMixMatFun, ad_tests) { Eigen::MatrixXd t(0, 0); expect_ad(f, t); expect_ad(fjit, t); + expect_ad_matvar(f, t); Eigen::MatrixXd u(1, 1); u << 2; expect_ad(f, u); expect_ad(fjit, u); + expect_ad_matvar(f, u); Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); - expect_ad(fjit, v); v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; expect_ad(f, v); expect_ad(fjit, v); + expect_ad_matvar(f, v); // issues around zero require looser tolerances for hessians stan::test::ad_tolerances tols; @@ -38,6 +40,7 @@ TEST(mathMixMatFun, ad_tests) { w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; expect_ad(tols, f, w); expect_ad(tols, fjit, w); + expect_ad_matvar(f, w); Eigen::MatrixXd z(2, 2); z << 1, 2, 5, std::numeric_limits::quiet_NaN(); @@ -50,6 +53,7 @@ TEST(mathMixMatFun, ad_tests) { a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); expect_ad(f, a); expect_ad(fjit, a); + expect_ad_matvar(f, a); // singular matrix, should use the // alias to input small amount of jitter on the diagonal From 709f1df42b0bf125e56222fa2bcbe25c64e2a119 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:13:01 -0500 Subject: [PATCH 071/102] cleanup --- stan/math/rev/fun/generalized_inverse.hpp | 1 + test/unit/math/mix/fun/generalized_inverse_test.cpp | 1 + 2 files changed, 2 insertions(+) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 087775249e0..ceaf664b723 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -104,6 +104,7 @@ inline auto generalized_inverse(const VarMat& G) { template * = nullptr> inline auto generalized_inverse(const VarMat& G, const double a) { using value_t = value_type_t; + using ret_type = promote_var_matrix_t; if (G.size() == 0) return G; diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index e27b5ef2dd3..ab7c7ddd5de 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -4,6 +4,7 @@ TEST(mathMixMatFun, ad_tests) { using stan::test::expect_ad; + using stan::test::expect_ad_matvar; auto f = [](const auto& G) { return stan::math::generalized_inverse(G); }; auto fjit From 5f8f09b63cf6b6b3b4d96c6f6067ac0bfd170e54 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 14 Dec 2020 12:17:02 +0000 Subject: [PATCH 072/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 2 +- stan/math/rev/fun/generalized_inverse.hpp | 30 ++++++++++++++-------- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 6558f298427..f47cab0b760 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -42,7 +42,7 @@ namespace math { * matrix. */ template * = nullptr, - require_not_vt_var* = nullptr> + require_not_vt_var* = nullptr> inline Eigen::Matrix, EigMat::RowsAtCompileTime, EigMat::ColsAtCompileTime> generalized_inverse(const EigMat& G) { diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index ceaf664b723..3708e93f47c 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -49,7 +49,8 @@ inline auto generalized_inverse(const VarMat& G) { if (G.rows() < G.cols()) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val()); - arena_t inv_G(transpose(mdivide_left_spd(A_spd.val_op(), G_arena.val_op()))); + arena_t inv_G( + transpose(mdivide_left_spd(A_spd.val_op(), G_arena.val_op()))); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -59,14 +60,17 @@ inline auto generalized_inverse(const VarMat& G) { reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); - G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj().transpose() * crossprod(inv_G.val_op()); + G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() + * PG.val_op(); + G_arena.adj() += GP.val_op() * G_arena.adj().transpose() + * crossprod(inv_G.val_op()); }); return ret_type(inv_G); } else { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val_op()); - arena_t inv_G(transpose(mdivide_right_spd(G_arena.val_op(), A_spd.val_op()))); + arena_t inv_G( + transpose(mdivide_right_spd(G_arena.val_op(), A_spd.val_op()))); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -76,8 +80,10 @@ inline auto generalized_inverse(const VarMat& G) { reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); - G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj().transpose() * crossprod(inv_G.val_op()); + G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() + * PG.val_op(); + G_arena.adj() += GP.val_op() * G_arena.adj().transpose() + * crossprod(inv_G.val_op()); }); return ret_type(inv_G); } @@ -126,8 +132,10 @@ inline auto generalized_inverse(const VarMat& G, const double a) { reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); - G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj().transpose() * crossprod(inv_G.val_op()); + G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() + * PG.val_op(); + G_arena.adj() += GP.val_op() * G_arena.adj().transpose() + * crossprod(inv_G.val_op()); }); return ret_type(inv_G); } else { @@ -144,8 +152,10 @@ inline auto generalized_inverse(const VarMat& G, const double a) { reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); - G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj().transpose() * crossprod(inv_G.val_op()); + G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() + * PG.val_op(); + G_arena.adj() += GP.val_op() * G_arena.adj().transpose() + * crossprod(inv_G.val_op()); }); return ret_type(inv_G); } From 483d527772cc9ede980457c03e8c75e795c08b69 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:27:15 -0500 Subject: [PATCH 073/102] val to val_op for bug in Eigen's CWiseUnaryView --- stan/math/rev/fun/generalized_inverse.hpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 3708e93f47c..ac751fcf5e2 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -48,7 +48,7 @@ inline auto generalized_inverse(const VarMat& G) { if (G.rows() < G.cols()) { arena_t G_arena(G); - auto A_spd = tcrossprod(G_arena.val()); + auto A_spd = tcrossprod(G_arena.val_op()); arena_t inv_G( transpose(mdivide_left_spd(A_spd.val_op(), G_arena.val_op()))); @@ -120,9 +120,9 @@ inline auto generalized_inverse(const VarMat& G, const double a) { if (G.rows() < G.cols()) { arena_t G_arena(G); - auto A_spd = tcrossprod(G_arena.val()); + auto A_spd = tcrossprod(G_arena.val_op()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val()))); + arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val_op()))); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -140,9 +140,9 @@ inline auto generalized_inverse(const VarMat& G, const double a) { return ret_type(inv_G); } else { arena_t G_arena(G); - auto A_spd = crossprod(G_arena.val()); + auto A_spd = crossprod(G_arena.val_op()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_right_spd(G_arena.val(), A_spd))); + arena_t inv_G(transpose(mdivide_right_spd(G_arena.val_op(), A_spd))); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; From 2c349dfc89cac1844f9855f278724637330d4165 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 14 Dec 2020 12:28:13 +0000 Subject: [PATCH 074/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index ac751fcf5e2..6cb1e891906 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -142,7 +142,8 @@ inline auto generalized_inverse(const VarMat& G, const double a) { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val_op()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_right_spd(G_arena.val_op(), A_spd))); + arena_t inv_G( + transpose(mdivide_right_spd(G_arena.val_op(), A_spd))); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; From 873c1f359c3906bd6e6445b1cda704f22f1100fa Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:31:34 -0500 Subject: [PATCH 075/102] more _op changes --- stan/math/rev/fun/generalized_inverse.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 6cb1e891906..24069391ee0 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -62,7 +62,7 @@ inline auto generalized_inverse(const VarMat& G) { G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj().transpose() + G_arena.adj() += GP.val_op() * G_arena.adj_op().transpose() * crossprod(inv_G.val_op()); }); return ret_type(inv_G); @@ -82,7 +82,7 @@ inline auto generalized_inverse(const VarMat& G) { G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj().transpose() + G_arena.adj() += GP.val_op() * G_arena.adj_op().transpose() * crossprod(inv_G.val_op()); }); return ret_type(inv_G); @@ -134,7 +134,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj().transpose() + G_arena.adj() += GP.val_op() * G_arena.adj_op().transpose() * crossprod(inv_G.val_op()); }); return ret_type(inv_G); @@ -155,7 +155,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj().transpose() + G_arena.adj() += GP.val_op() * G_arena.adj_op().transpose() * crossprod(inv_G.val_op()); }); return ret_type(inv_G); From e3f14c201a9fab311696acbf051d9ce706aaf71d Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:35:30 -0500 Subject: [PATCH 076/102] flip rows/cols inline --- stan/math/prim/fun/generalized_inverse.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index f47cab0b760..c9d0810ae44 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -43,8 +43,8 @@ namespace math { */ template * = nullptr, require_not_vt_var* = nullptr> -inline Eigen::Matrix, EigMat::RowsAtCompileTime, - EigMat::ColsAtCompileTime> +inline Eigen::Matrix, EigMat::ColsAtCompileTime, + EigMat::RowsAtCompileTime> generalized_inverse(const EigMat& G) { using value_t = value_type_t; From f56b9c586ff0c38054914edd29f29818f628e564 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:42:42 -0500 Subject: [PATCH 077/102] small fix remove n, m using .rows()/.cols() --- stan/math/prim/fun/generalized_inverse.hpp | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index c9d0810ae44..23a50a07442 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -51,15 +51,12 @@ generalized_inverse(const EigMat& G) { if (G.size() == 0) return {}; - const auto n = G.rows(); - const auto m = G.cols(); - if (G.rows() == G.cols()) return inverse(G); const auto& G_ref = to_ref(G); - if (n < m) { + if (G.rows() < G.cols()) { return transpose(mdivide_left_spd(tcrossprod(G_ref), G_ref)); } else { return transpose(mdivide_right_spd(G_ref, crossprod(G_ref))); @@ -100,15 +97,12 @@ generalized_inverse(const EigMat& G, const double a) { if (G.size() == 0) return {}; - const auto n = G.rows(); - const auto m = G.cols(); - if (G.rows() == G.cols()) return inverse(G); const auto& G_ref = to_ref(G); - if (n < m) { + if (G.rows() < G.cols()) { Eigen::Matrix A = tcrossprod(G_ref); A.diagonal().array() += a; From e768431ccbcbe157c2ac06d0f09a9858638251e6 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:45:51 -0500 Subject: [PATCH 078/102] small fix --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 24069391ee0..902de90596a 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -113,7 +113,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { using ret_type = promote_var_matrix_t; if (G.size() == 0) - return G; + return ret_type(G); if (G.rows() == G.cols()) return inverse(G); From 2b3ab6f92a8f4bb128ec7ee369fbdc8ce09e927e Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:47:36 -0500 Subject: [PATCH 079/102] more fixes --- stan/math/rev/fun/generalized_inverse.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 902de90596a..736445d06e2 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -44,7 +44,7 @@ inline auto generalized_inverse(const VarMat& G) { return ret_type(G); if (G.rows() == G.cols()) - return inverse(G); + return ret_type(inverse(G)); if (G.rows() < G.cols()) { arena_t G_arena(G); @@ -116,7 +116,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { return ret_type(G); if (G.rows() == G.cols()) - return inverse(G); + return ret_type(inverse(G)); if (G.rows() < G.cols()) { arena_t G_arena(G); From 76b5536a862b4b7ec3ed1a8f1caa7048ba62286e Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 07:59:05 -0500 Subject: [PATCH 080/102] more fixes --- stan/math/rev/fun/generalized_inverse.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 736445d06e2..bbee8f68951 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -50,7 +50,7 @@ inline auto generalized_inverse(const VarMat& G) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val_op()); arena_t inv_G( - transpose(mdivide_left_spd(A_spd.val_op(), G_arena.val_op()))); + mdivide_left_spd(A_spd.val_op(), G_arena.val_op())).transpose(); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -70,7 +70,7 @@ inline auto generalized_inverse(const VarMat& G) { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val_op()); arena_t inv_G( - transpose(mdivide_right_spd(G_arena.val_op(), A_spd.val_op()))); + mdivide_right_spd(G_arena.val_op(), A_spd.val_op())).transpose(); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -122,7 +122,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val_op()); A_spd.diagonal().array() += a; - arena_t inv_G(transpose(mdivide_left_spd(A_spd, G_arena.val_op()))); + arena_t inv_G(mdivide_left_spd(A_spd, G_arena.val_op())).transpose(); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -143,7 +143,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { auto A_spd = crossprod(G_arena.val_op()); A_spd.diagonal().array() += a; arena_t inv_G( - transpose(mdivide_right_spd(G_arena.val_op(), A_spd))); + mdivide_right_spd(G_arena.val_op(), A_spd)).transpose(); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; From 2e3e69cdd42887327cf315b1d59e91595afb4852 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 14 Dec 2020 13:00:02 +0000 Subject: [PATCH 081/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index bbee8f68951..ce01a0b1447 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -49,8 +49,8 @@ inline auto generalized_inverse(const VarMat& G) { if (G.rows() < G.cols()) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val_op()); - arena_t inv_G( - mdivide_left_spd(A_spd.val_op(), G_arena.val_op())).transpose(); + arena_t inv_G(mdivide_left_spd(A_spd.val_op(), G_arena.val_op())) + .transpose(); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -69,8 +69,8 @@ inline auto generalized_inverse(const VarMat& G) { } else { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val_op()); - arena_t inv_G( - mdivide_right_spd(G_arena.val_op(), A_spd.val_op())).transpose(); + arena_t inv_G(mdivide_right_spd(G_arena.val_op(), A_spd.val_op())) + .transpose(); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -122,7 +122,8 @@ inline auto generalized_inverse(const VarMat& G, const double a) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val_op()); A_spd.diagonal().array() += a; - arena_t inv_G(mdivide_left_spd(A_spd, G_arena.val_op())).transpose(); + arena_t inv_G(mdivide_left_spd(A_spd, G_arena.val_op())) + .transpose(); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -142,8 +143,8 @@ inline auto generalized_inverse(const VarMat& G, const double a) { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val_op()); A_spd.diagonal().array() += a; - arena_t inv_G( - mdivide_right_spd(G_arena.val_op(), A_spd)).transpose(); + arena_t inv_G(mdivide_right_spd(G_arena.val_op(), A_spd)) + .transpose(); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; From a118eaad74ad516b84cb6a2e13067b0e14f9b8d7 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 08:02:10 -0500 Subject: [PATCH 082/102] more cleanup --- stan/math/rev/fun/generalized_inverse.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index bbee8f68951..9a49062efd0 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -50,7 +50,7 @@ inline auto generalized_inverse(const VarMat& G) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val_op()); arena_t inv_G( - mdivide_left_spd(A_spd.val_op(), G_arena.val_op())).transpose(); + mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -70,7 +70,7 @@ inline auto generalized_inverse(const VarMat& G) { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val_op()); arena_t inv_G( - mdivide_right_spd(G_arena.val_op(), A_spd.val_op())).transpose(); + mdivide_right_spd(G_arena.val_op(), A_spd.val_op()).transpose()); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -122,7 +122,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val_op()); A_spd.diagonal().array() += a; - arena_t inv_G(mdivide_left_spd(A_spd, G_arena.val_op())).transpose(); + arena_t inv_G(mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -143,7 +143,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { auto A_spd = crossprod(G_arena.val_op()); A_spd.diagonal().array() += a; arena_t inv_G( - mdivide_right_spd(G_arena.val_op(), A_spd)).transpose(); + mdivide_right_spd(G_arena.val_op(), A_spd.val_op()).transpose()); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; From 9a9be26445f55f60982e548bebef9db1b96f17ce Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 14 Dec 2020 13:04:42 +0000 Subject: [PATCH 083/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 570e7c023e5..472403b8c47 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -132,7 +132,8 @@ inline auto generalized_inverse(const VarMat& G, const double a) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val_op()); A_spd.diagonal().array() += a; - arena_t inv_G(mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); + arena_t inv_G( + mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; From a37536dbd31ed986043cad077a4608f9d0f8cd0d Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 08:05:51 -0500 Subject: [PATCH 084/102] fixed merge conflict --- stan/math/rev/fun/generalized_inverse.hpp | 10 ---------- 1 file changed, 10 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 570e7c023e5..9a49062efd0 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -49,13 +49,8 @@ inline auto generalized_inverse(const VarMat& G) { if (G.rows() < G.cols()) { arena_t G_arena(G); auto A_spd = tcrossprod(G_arena.val_op()); -<<<<<<< HEAD arena_t inv_G( mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); -======= - arena_t inv_G(mdivide_left_spd(A_spd.val_op(), G_arena.val_op())) - .transpose(); ->>>>>>> 2e3e69cdd42887327cf315b1d59e91595afb4852 auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -74,13 +69,8 @@ inline auto generalized_inverse(const VarMat& G) { } else { arena_t G_arena(G); auto A_spd = crossprod(G_arena.val_op()); -<<<<<<< HEAD arena_t inv_G( mdivide_right_spd(G_arena.val_op(), A_spd.val_op()).transpose()); -======= - arena_t inv_G(mdivide_right_spd(G_arena.val_op(), A_spd.val_op())) - .transpose(); ->>>>>>> 2e3e69cdd42887327cf315b1d59e91595afb4852 auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; From 6b136cdaeab128a9923c619290d95195b13306c5 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 08:15:04 -0500 Subject: [PATCH 085/102] comment out equal tests --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index ab7c7ddd5de..72b3d1b3b8e 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -12,6 +12,7 @@ TEST(mathMixMatFun, ad_tests) { auto fjit2 = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-4); }; +/* Eigen::MatrixXd t(0, 0); expect_ad(f, t); expect_ad(fjit, t); @@ -22,7 +23,7 @@ TEST(mathMixMatFun, ad_tests) { expect_ad(f, u); expect_ad(fjit, u); expect_ad_matvar(f, u); - +*/ Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); From 27b688a1aabeccec58ab7396e3bc7aca1d7801f7 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 14 Dec 2020 13:15:56 +0000 Subject: [PATCH 086/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- .../math/mix/fun/generalized_inverse_test.cpp | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 72b3d1b3b8e..15e520ea311 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -12,18 +12,18 @@ TEST(mathMixMatFun, ad_tests) { auto fjit2 = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-4); }; -/* - Eigen::MatrixXd t(0, 0); - expect_ad(f, t); - expect_ad(fjit, t); - expect_ad_matvar(f, t); + /* + Eigen::MatrixXd t(0, 0); + expect_ad(f, t); + expect_ad(fjit, t); + expect_ad_matvar(f, t); - Eigen::MatrixXd u(1, 1); - u << 2; - expect_ad(f, u); - expect_ad(fjit, u); - expect_ad_matvar(f, u); -*/ + Eigen::MatrixXd u(1, 1); + u << 2; + expect_ad(f, u); + expect_ad(fjit, u); + expect_ad_matvar(f, u); + */ Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); From 29d536ef29fd036bfc1e9eff916d358bade19295 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 08:16:58 -0500 Subject: [PATCH 087/102] test --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 72b3d1b3b8e..e352d875e66 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -23,11 +23,11 @@ TEST(mathMixMatFun, ad_tests) { expect_ad(f, u); expect_ad(fjit, u); expect_ad_matvar(f, u); -*/ + Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); - +*/ v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; expect_ad(f, v); expect_ad(fjit, v); From e7af4fcdefdb852c24cd0395e2c235f31e4119f2 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 08:17:53 -0500 Subject: [PATCH 088/102] fix --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index e352d875e66..b130c8ed726 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -12,11 +12,11 @@ TEST(mathMixMatFun, ad_tests) { auto fjit2 = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-4); }; -/* - Eigen::MatrixXd t(0, 0); - expect_ad(f, t); - expect_ad(fjit, t); - expect_ad_matvar(f, t); + /* + Eigen::MatrixXd t(0, 0); + expect_ad(f, t); + expect_ad(fjit, t); + expect_ad_matvar(f, t); Eigen::MatrixXd u(1, 1); u << 2; From b3534f5968b1ffa627daad8bdbf58b878f9b3df7 Mon Sep 17 00:00:00 2001 From: spinkney Date: Mon, 14 Dec 2020 08:22:22 -0500 Subject: [PATCH 089/102] Update generalized_inverse_test.cpp --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index b130c8ed726..53f8c3fe8b8 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -23,11 +23,11 @@ TEST(mathMixMatFun, ad_tests) { expect_ad(f, u); expect_ad(fjit, u); expect_ad_matvar(f, u); - +*/ Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); -*/ +/* v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; expect_ad(f, v); expect_ad(fjit, v); @@ -68,4 +68,5 @@ TEST(mathMixMatFun, ad_tests) { tols3.hessian_hessian_ = 0.01; tols3.hessian_fvar_hessian_ = 0.01; expect_ad(tols3, fjit2, m); + */ } From 9c8f0b99917b4cdf573c4572ccaf4422437ddb1b Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Mon, 14 Dec 2020 13:23:06 +0000 Subject: [PATCH 090/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- .../math/mix/fun/generalized_inverse_test.cpp | 70 +++++++++---------- 1 file changed, 35 insertions(+), 35 deletions(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 53f8c3fe8b8..058cd7027c0 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -27,46 +27,46 @@ TEST(mathMixMatFun, ad_tests) { Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); -/* - v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; - expect_ad(f, v); - expect_ad(fjit, v); - expect_ad_matvar(f, v); + /* + v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; + expect_ad(f, v); + expect_ad(fjit, v); + expect_ad_matvar(f, v); - // issues around zero require looser tolerances for hessians - stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 2.0; - tols.hessian_fvar_hessian_ = 2.0; + // issues around zero require looser tolerances for hessians + stan::test::ad_tolerances tols; + tols.hessian_hessian_ = 2.0; + tols.hessian_fvar_hessian_ = 2.0; - Eigen::MatrixXd w(3, 4); - w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; - expect_ad(tols, f, w); - expect_ad(tols, fjit, w); - expect_ad_matvar(f, w); + Eigen::MatrixXd w(3, 4); + w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; + expect_ad(tols, f, w); + expect_ad(tols, fjit, w); + expect_ad_matvar(f, w); - Eigen::MatrixXd z(2, 2); - z << 1, 2, 5, std::numeric_limits::quiet_NaN(); - EXPECT_NO_THROW(stan::math::generalized_inverse(z)); - EXPECT_NO_THROW(stan::math::generalized_inverse(z, 0.0)); + Eigen::MatrixXd z(2, 2); + z << 1, 2, 5, std::numeric_limits::quiet_NaN(); + EXPECT_NO_THROW(stan::math::generalized_inverse(z)); + EXPECT_NO_THROW(stan::math::generalized_inverse(z, 0.0)); - // autodiff throws, so following fails (throw behavior must match to pass) + // autodiff throws, so following fails (throw behavior must match to pass) - Eigen::MatrixXd a(2, 2); - a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); - expect_ad(f, a); - expect_ad(fjit, a); - expect_ad_matvar(f, a); + Eigen::MatrixXd a(2, 2); + a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); + expect_ad(f, a); + expect_ad(fjit, a); + expect_ad_matvar(f, a); - // singular matrix, should use the - // alias to input small amount of jitter on the diagonal - Eigen::MatrixXd m(3, 2); - m << 1, 2, 2, 4, 1, 2; - EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); + // singular matrix, should use the + // alias to input small amount of jitter on the diagonal + Eigen::MatrixXd m(3, 2); + m << 1, 2, 2, 4, 1, 2; + EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); - // should work with jittered version - stan::test::ad_tolerances tols3; - tols3.hessian_hessian_ = 0.01; - tols3.hessian_fvar_hessian_ = 0.01; - expect_ad(tols3, fjit2, m); - */ + // should work with jittered version + stan::test::ad_tolerances tols3; + tols3.hessian_hessian_ = 0.01; + tols3.hessian_fvar_hessian_ = 0.01; + expect_ad(tols3, fjit2, m); + */ } From 438bbf3fa0d85742057761cc87d0740a0a86ea56 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Mon, 14 Dec 2020 19:32:02 -0500 Subject: [PATCH 091/102] Finished tests! --- stan/math/rev/fun/generalized_inverse.hpp | 50 +++++------ .../math/mix/fun/generalized_inverse_test.cpp | 85 ++++++++++--------- 2 files changed, 68 insertions(+), 67 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 8a249c862f2..75299e16056 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -52,19 +52,19 @@ inline auto generalized_inverse(const VarMat& G) { arena_t inv_G( mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); + arena_t res = inv_G; + auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); GP.diagonal().array() += 1.0; - reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { - G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); - G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() - * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj_op().transpose() - * crossprod(inv_G.val_op()); - }); + reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { + G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() * inv_G.val_op().transpose(); + G_arena.adj() += PG.val_op() * res.adj_op().transpose() * tcrossprod(inv_G.val_op()); + G_arena.adj() += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); + }); return ret_type(inv_G); } else { arena_t G_arena(G); @@ -72,18 +72,18 @@ inline auto generalized_inverse(const VarMat& G) { arena_t inv_G( mdivide_right_spd(G_arena.val_op(), A_spd.val_op()).transpose()); + arena_t res = inv_G; + auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); GP.diagonal().array() += 1.0; - reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { - G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); - G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() - * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj_op().transpose() - * crossprod(inv_G.val_op()); + reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { + G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() * inv_G.val_op().transpose(); + G_arena.adj() += PG.val_op() * res.adj_op().transpose() * tcrossprod(inv_G.val_op()); + G_arena.adj() += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); }); return ret_type(inv_G); } @@ -125,18 +125,18 @@ inline auto generalized_inverse(const VarMat& G, const double a) { arena_t inv_G( mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); + arena_t res = inv_G; + auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); GP.diagonal().array() += 1.0; - reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { - G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); - G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() - * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj_op().transpose() - * crossprod(inv_G.val_op()); + reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { + G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() * inv_G.val_op().transpose(); + G_arena.adj() += PG.val_op() * res.adj_op().transpose() * tcrossprod(inv_G.val_op()); + G_arena.adj() += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); }); return ret_type(inv_G); } else { @@ -145,6 +145,8 @@ inline auto generalized_inverse(const VarMat& G, const double a) { A_spd.diagonal().array() += a; arena_t inv_G( mdivide_right_spd(G_arena.val_op(), A_spd.val_op()).transpose()); + + arena_t res = inv_G; auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -152,12 +154,10 @@ inline auto generalized_inverse(const VarMat& G, const double a) { auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); GP.diagonal().array() += 1.0; - reverse_pass_callback([G_arena, inv_G, GP, PG]() mutable { - G_arena.adj() += -inv_G.val_op() * G_arena.adj_op() * inv_G.val_op(); - G_arena.adj() += tcrossprod(inv_G.val_op()) * G_arena.adj_op().transpose() - * PG.val_op(); - G_arena.adj() += GP.val_op() * G_arena.adj_op().transpose() - * crossprod(inv_G.val_op()); + reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { + G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() * inv_G.val_op().transpose(); + G_arena.adj() += PG.val_op() * res.adj_op().transpose() * tcrossprod(inv_G.val_op()); + G_arena.adj() += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); }); return ret_type(inv_G); } diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 058cd7027c0..fdb2c871c03 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -10,63 +10,64 @@ TEST(mathMixMatFun, ad_tests) { auto fjit = [](const auto& G) { return stan::math::generalized_inverse(G, 0.0); }; auto fjit2 - = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-4); }; + = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-6); }; - /* - Eigen::MatrixXd t(0, 0); - expect_ad(f, t); - expect_ad(fjit, t); - expect_ad_matvar(f, t); + Eigen::MatrixXd t(0, 0); + expect_ad(f, t); + expect_ad(fjit, t); + expect_ad_matvar(f, t); Eigen::MatrixXd u(1, 1); u << 2; expect_ad(f, u); expect_ad(fjit, u); expect_ad_matvar(f, u); -*/ + Eigen::MatrixXd v(2, 3); v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); - /* - v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; - expect_ad(f, v); - expect_ad(fjit, v); - expect_ad_matvar(f, v); + expect_ad_matvar(f, v); + + v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; + expect_ad(f, v); + expect_ad(fjit, v); + expect_ad_matvar(f, v); + + // issues around zero require looser tolerances for hessians + stan::test::ad_tolerances tols; + tols.hessian_hessian_ = 2.0; + tols.hessian_fvar_hessian_ = 2.0; - // issues around zero require looser tolerances for hessians - stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 2.0; - tols.hessian_fvar_hessian_ = 2.0; + Eigen::MatrixXd w(3, 4); + w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; + expect_ad(tols, f, w); + expect_ad(tols, fjit, w); + expect_ad_matvar(f, w); - Eigen::MatrixXd w(3, 4); - w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; - expect_ad(tols, f, w); - expect_ad(tols, fjit, w); - expect_ad_matvar(f, w); + Eigen::MatrixXd z(2, 2); + z << 1, 2, 5, std::numeric_limits::quiet_NaN(); + EXPECT_NO_THROW(stan::math::generalized_inverse(z)); + EXPECT_NO_THROW(stan::math::generalized_inverse(z, 0.0)); - Eigen::MatrixXd z(2, 2); - z << 1, 2, 5, std::numeric_limits::quiet_NaN(); - EXPECT_NO_THROW(stan::math::generalized_inverse(z)); - EXPECT_NO_THROW(stan::math::generalized_inverse(z, 0.0)); + // autodiff throws, so following fails (throw behavior must match to pass) - // autodiff throws, so following fails (throw behavior must match to pass) + Eigen::MatrixXd a(2, 2); + a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); + expect_ad(f, a); + expect_ad(fjit, a); + expect_ad_matvar(f, a); - Eigen::MatrixXd a(2, 2); - a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); - expect_ad(f, a); - expect_ad(fjit, a); - expect_ad_matvar(f, a); + // singular matrix, should use the + // alias to input small amount of jitter on the diagonal - // singular matrix, should use the - // alias to input small amount of jitter on the diagonal - Eigen::MatrixXd m(3, 2); - m << 1, 2, 2, 4, 1, 2; - EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); + Eigen::MatrixXd m(3, 2); + m << 1, 2, 2, 4, 1, 2; + EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); - // should work with jittered version - stan::test::ad_tolerances tols3; - tols3.hessian_hessian_ = 0.01; - tols3.hessian_fvar_hessian_ = 0.01; - expect_ad(tols3, fjit2, m); - */ + // should work with jittered version + stan::test::ad_tolerances tols3; + tols3.hessian_hessian_ = 0.01; + tols3.hessian_fvar_hessian_ = 0.01; + expect_ad_matvar(tols3, fjit2, m); + } From ec80ff1cbe8e156552a1787463c00638a1cb8893 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 15 Dec 2020 00:32:56 +0000 Subject: [PATCH 092/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 42 ++++++++++++------- .../math/mix/fun/generalized_inverse_test.cpp | 3 +- 2 files changed, 28 insertions(+), 17 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 75299e16056..6a8dd697e3e 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -53,7 +53,7 @@ inline auto generalized_inverse(const VarMat& G) { mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); arena_t res = inv_G; - + auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); PG.diagonal().array() += 1.0; @@ -61,10 +61,13 @@ inline auto generalized_inverse(const VarMat& G) { GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { - G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() * inv_G.val_op().transpose(); - G_arena.adj() += PG.val_op() * res.adj_op().transpose() * tcrossprod(inv_G.val_op()); - G_arena.adj() += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); - }); + G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() + * inv_G.val_op().transpose(); + G_arena.adj() += PG.val_op() * res.adj_op().transpose() + * tcrossprod(inv_G.val_op()); + G_arena.adj() + += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); + }); return ret_type(inv_G); } else { arena_t G_arena(G); @@ -81,9 +84,12 @@ inline auto generalized_inverse(const VarMat& G) { GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { - G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() * inv_G.val_op().transpose(); - G_arena.adj() += PG.val_op() * res.adj_op().transpose() * tcrossprod(inv_G.val_op()); - G_arena.adj() += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); + G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() + * inv_G.val_op().transpose(); + G_arena.adj() += PG.val_op() * res.adj_op().transpose() + * tcrossprod(inv_G.val_op()); + G_arena.adj() + += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); }); return ret_type(inv_G); } @@ -134,9 +140,12 @@ inline auto generalized_inverse(const VarMat& G, const double a) { GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { - G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() * inv_G.val_op().transpose(); - G_arena.adj() += PG.val_op() * res.adj_op().transpose() * tcrossprod(inv_G.val_op()); - G_arena.adj() += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); + G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() + * inv_G.val_op().transpose(); + G_arena.adj() += PG.val_op() * res.adj_op().transpose() + * tcrossprod(inv_G.val_op()); + G_arena.adj() + += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); }); return ret_type(inv_G); } else { @@ -145,7 +154,7 @@ inline auto generalized_inverse(const VarMat& G, const double a) { A_spd.diagonal().array() += a; arena_t inv_G( mdivide_right_spd(G_arena.val_op(), A_spd.val_op()).transpose()); - + arena_t res = inv_G; auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); @@ -155,9 +164,12 @@ inline auto generalized_inverse(const VarMat& G, const double a) { GP.diagonal().array() += 1.0; reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { - G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() * inv_G.val_op().transpose(); - G_arena.adj() += PG.val_op() * res.adj_op().transpose() * tcrossprod(inv_G.val_op()); - G_arena.adj() += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); + G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() + * inv_G.val_op().transpose(); + G_arena.adj() += PG.val_op() * res.adj_op().transpose() + * tcrossprod(inv_G.val_op()); + G_arena.adj() + += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); }); return ret_type(inv_G); } diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index fdb2c871c03..f7ceea1dcef 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -27,7 +27,7 @@ TEST(mathMixMatFun, ad_tests) { v << 1, 3, 5, 2, 4, 6; expect_ad(f, v); expect_ad_matvar(f, v); - + v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; expect_ad(f, v); expect_ad(fjit, v); @@ -69,5 +69,4 @@ TEST(mathMixMatFun, ad_tests) { tols3.hessian_hessian_ = 0.01; tols3.hessian_fvar_hessian_ = 0.01; expect_ad_matvar(tols3, fjit2, m); - } From d1c1d2bf5d86e787e2057728e2b1532d2a8ac3e8 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Mon, 14 Dec 2020 19:38:46 -0500 Subject: [PATCH 093/102] Update doc for reverse specialization --- stan/math/rev/fun/generalized_inverse.hpp | 27 ++++++++++++++++------- 1 file changed, 19 insertions(+), 8 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 75299e16056..5eb6ba09e9b 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -17,11 +17,16 @@ namespace stan { namespace math { /* - * Reverse mode differentiation algorithm reference: + * Reverse mode specialization of calculating the generalized inverse of a matrix. * - * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems - * Whose Variables Separate. Author(s): G. H. Golub and V. Pereyra. Source: SIAM - * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 + * @param G specified matrix + * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has + * size zero). + * + * @note Reverse mode differentiation algorithm reference: + * + *
  • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems + * Whose Variables Separate. SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432
* * Equation 4.12 in the paper * @@ -90,13 +95,19 @@ inline auto generalized_inverse(const VarMat& G) { } /* - * Reverse mode differentiation algorithm reference: + * Reverse mode specialization of calculating the generalized inverse of a matrix. + * + * @param G specified matrix + * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has + * size zero). * - * The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems - * Whose Variables Separate. Author(s): G. H. Golub and V. Pereyra. Source: SIAM - * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432 + * @note Reverse mode differentiation algorithm reference: + * + *
  • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems + * Whose Variables Separate. SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432
* * Equation 4.12 in the paper + * * See also * http://mathoverflow.net/questions/25778/analytical-formula-for-numerical-derivative-of-the-matrix-pseudo-inverse * From b25eb0a63aef5a5dd22023ecec9f834c416e61c6 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 15 Dec 2020 00:39:35 +0000 Subject: [PATCH 094/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 26 ++++++++++++++--------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index b1b5bad0716..c595dea2d3f 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -17,16 +17,19 @@ namespace stan { namespace math { /* - * Reverse mode specialization of calculating the generalized inverse of a matrix. + * Reverse mode specialization of calculating the generalized inverse of a + * matrix. * * @param G specified matrix - * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has - * size zero). + * @return Generalized inverse of the matrix (an empty matrix if the specified + * matrix has size zero). * * @note Reverse mode differentiation algorithm reference: * - *
  • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems - * Whose Variables Separate. SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432
+ *
  • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses + * and Nonlinear Least Squares Problems Whose Variables Separate. SIAM + * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. + * 413-432
* * Equation 4.12 in the paper * @@ -101,16 +104,19 @@ inline auto generalized_inverse(const VarMat& G) { } /* - * Reverse mode specialization of calculating the generalized inverse of a matrix. + * Reverse mode specialization of calculating the generalized inverse of a + * matrix. * * @param G specified matrix - * @return Generalized inverse of the matrix (an empty matrix if the specified matrix has - * size zero). + * @return Generalized inverse of the matrix (an empty matrix if the specified + * matrix has size zero). * * @note Reverse mode differentiation algorithm reference: * - *
  • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses and Nonlinear Least Squares Problems - * Whose Variables Separate. SIAM Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. 413-432
+ *
  • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses + * and Nonlinear Least Squares Problems Whose Variables Separate. SIAM + * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. + * 413-432
* * Equation 4.12 in the paper * From 361d5836347d8d0b398c415fc45dbe0119dc4102 Mon Sep 17 00:00:00 2001 From: spinkney Date: Thu, 17 Dec 2020 07:31:49 -0500 Subject: [PATCH 095/102] updates using ldlt lots of changes. Removed the jitter as it's no longer needed. Updated everything to use ldlt --- stan/math/prim/fun/generalized_inverse.hpp | 75 +------- stan/math/rev/fun/generalized_inverse.hpp | 160 +++--------------- .../math/mix/fun/generalized_inverse_test.cpp | 36 ++-- .../prim/fun/generalized_inverse_test.cpp | 30 ++-- 4 files changed, 65 insertions(+), 236 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 23a50a07442..8eba7815ca4 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -5,13 +5,6 @@ #include #include #include -#include -#include -#include -#include -#include -#include -#include #include namespace stan { @@ -32,14 +25,6 @@ namespace math { * @param G specified matrix * @return Generalized inverse of the matrix (an empty matrix if the specified * matrix has size zero). - * @note Because the method inverts a SPD matrix internally that interal matrix - may result in small numerical issues that result in a non-SPD error. There are - two - * generalized_inverse functions, one that uses one input matrix - (this one) - * and another that works with an input matrix and a small jitter to the - diagonal of the internal SPD - * matrix. */ template * = nullptr, require_not_vt_var* = nullptr> @@ -57,60 +42,14 @@ generalized_inverse(const EigMat& G) { const auto& G_ref = to_ref(G); if (G.rows() < G.cols()) { - return transpose(mdivide_left_spd(tcrossprod(G_ref), G_ref)); + return (G_ref * G_ref.transpose()) + .ldlt() + .solve(G_ref) + .transpose(); } else { - return transpose(mdivide_right_spd(G_ref, crossprod(G_ref))); - } -} - -/** - * Returns the Moore-Penrose generalized inverse of the specified matrix. - * - * The method is based on the Cholesky computation of the transform as specified - in - * - *
  • Courrieu, Pierre. 2008. Fast Computation of Moore-Penrose Inverse - Matrices. - * arXiv 0804.4809
- * - * @tparam EigMat type of the matrix (must be derived from `Eigen::MatrixBase`) - * - * @param G specified matrix - * @param a real constant - * @return Generalized inverse of the matrix (an empty matrix if the specified - * matrix has size zero). - * @note Because the method inverts a SPD matrix internally that interal matrix - may result in small numerical issues that result in a non-SPD error. There are - two - * generalized_inverse functions, one that uses one input matrix - * and another (this one) that works with an input matrix and a small jitter to - the diagonal of the internal SPD - * matrix. - */ -template * = nullptr, - require_not_vt_var* = nullptr> -inline Eigen::Matrix, EigMat::RowsAtCompileTime, - EigMat::ColsAtCompileTime> -generalized_inverse(const EigMat& G, const double a) { - using value_t = value_type_t; - - if (G.size() == 0) - return {}; - - if (G.rows() == G.cols()) - return inverse(G); - - const auto& G_ref = to_ref(G); - - if (G.rows() < G.cols()) { - Eigen::Matrix A - = tcrossprod(G_ref); - A.diagonal().array() += a; - return transpose(mdivide_left_spd(A, G_ref)); - } else { - Eigen::Matrix A = crossprod(G_ref); - A.diagonal().array() += a; - return transpose(mdivide_right_spd(G_ref, A)); + return (G_ref.transpose() * G_ref) + .ldlt() + .solve(G_ref.transpose()); } } diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index c595dea2d3f..54e2bfb2d20 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -4,104 +4,38 @@ #include #include #include -#include -#include -#include -#include -#include -#include -#include #include namespace stan { namespace math { +namespace internal { /* * Reverse mode specialization of calculating the generalized inverse of a * matrix. - * - * @param G specified matrix - * @return Generalized inverse of the matrix (an empty matrix if the specified - * matrix has size zero). - * - * @note Reverse mode differentiation algorithm reference: - * *
  • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses * and Nonlinear Least Squares Problems Whose Variables Separate. SIAM * Journal on Numerical Analysis, Vol. 10, No. 2 (Apr., 1973), pp. * 413-432
- * - * Equation 4.12 in the paper - * - * See also - * http://mathoverflow.net/questions/25778/analytical-formula-for-numerical-derivative-of-the-matrix-pseudo-inverse - * - * Implementation is based on - * Title: Uncertainties Python Package - * Author: Eric O. LEBIGOT - * Date: 2020 - * Availability: - * https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py */ -template * = nullptr> -inline auto generalized_inverse(const VarMat& G) { - using value_t = value_type_t; - using ret_type = promote_var_matrix_t; - - if (G.size() == 0) - return ret_type(G); - - if (G.rows() == G.cols()) - return ret_type(inverse(G)); - - if (G.rows() < G.cols()) { - arena_t G_arena(G); - auto A_spd = tcrossprod(G_arena.val_op()); - arena_t inv_G( - mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); - - arena_t res = inv_G; - - auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); - PG.diagonal().array() += 1.0; - - auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); - GP.diagonal().array() += 1.0; - - reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { - G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() - * inv_G.val_op().transpose(); - G_arena.adj() += PG.val_op() * res.adj_op().transpose() - * tcrossprod(inv_G.val_op()); +template +inline auto generalized_inverse_lambda(T1& G_arena, T2& inv_G) { +return [G_arena, inv_G]() mutable { G_arena.adj() - += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); - }); - return ret_type(inv_G); - } else { - arena_t G_arena(G); - auto A_spd = crossprod(G_arena.val_op()); - arena_t inv_G( - mdivide_right_spd(G_arena.val_op(), A_spd.val_op()).transpose()); - - arena_t res = inv_G; - - auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); - PG.diagonal().array() += 1.0; - - auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); - GP.diagonal().array() += 1.0; - - reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { - G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() - * inv_G.val_op().transpose(); - G_arena.adj() += PG.val_op() * res.adj_op().transpose() - * tcrossprod(inv_G.val_op()); - G_arena.adj() - += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); - }); - return ret_type(inv_G); - } + += -(inv_G.val_op().transpose() * inv_G.adj_op() + * inv_G.val_op().transpose()) + + (-G_arena.val_op() * inv_G.val_op() + + Eigen::MatrixXd::Identity(G_arena.rows(), inv_G.cols())) + * inv_G.adj_op().transpose() * inv_G.val_op() + * inv_G.val_op().transpose() + + inv_G.val_op().transpose() * inv_G.val_op() + * inv_G.adj_op().transpose() + * (-inv_G.val_op() * G_arena.val_op() + + Eigen::MatrixXd::Identity(inv_G.rows(), + G_arena.cols())); + }; } +} // namespace internal /* * Reverse mode specialization of calculating the generalized inverse of a @@ -123,15 +57,9 @@ inline auto generalized_inverse(const VarMat& G) { * See also * http://mathoverflow.net/questions/25778/analytical-formula-for-numerical-derivative-of-the-matrix-pseudo-inverse * - * Implementation is based on - * Title: Uncertainties Python Package - * Author: Eric O. LEBIGOT - * Date: 2020 - * Availability: - * https://github.com/lebigot/uncertainties/blob/master/uncertainties/unumpy/core.py */ template * = nullptr> -inline auto generalized_inverse(const VarMat& G, const double a) { +inline auto generalized_inverse(const VarMat& G) { using value_t = value_type_t; using ret_type = promote_var_matrix_t; @@ -143,56 +71,22 @@ inline auto generalized_inverse(const VarMat& G, const double a) { if (G.rows() < G.cols()) { arena_t G_arena(G); - auto A_spd = tcrossprod(G_arena.val_op()); - A_spd.diagonal().array() += a; - arena_t inv_G( - mdivide_left_spd(A_spd.val_op(), G_arena.val_op()).transpose()); - - arena_t res = inv_G; - - auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); - PG.diagonal().array() += 1.0; - - auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); - GP.diagonal().array() += 1.0; - - reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { - G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() - * inv_G.val_op().transpose(); - G_arena.adj() += PG.val_op() * res.adj_op().transpose() - * tcrossprod(inv_G.val_op()); - G_arena.adj() - += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); - }); + arena_t inv_G((G_arena.val_op() * G_arena.val_op().transpose()) + .ldlt() + .solve(G_arena.val_op()) + .transpose()); + reverse_pass_callback(internal::generalized_inverse_lambda(G_arena, inv_G)); return ret_type(inv_G); } else { arena_t G_arena(G); - auto A_spd = crossprod(G_arena.val_op()); - A_spd.diagonal().array() += a; - arena_t inv_G( - mdivide_right_spd(G_arena.val_op(), A_spd.val_op()).transpose()); - - arena_t res = inv_G; - - auto PG = to_arena(-G_arena.val_op() * inv_G.val_op()); - PG.diagonal().array() += 1.0; - - auto GP = to_arena(-inv_G.val_op() * G_arena.val_op()); - GP.diagonal().array() += 1.0; - - reverse_pass_callback([G_arena, res, inv_G, GP, PG]() mutable { - G_arena.adj() -= inv_G.val_op().transpose() * res.adj_op() - * inv_G.val_op().transpose(); - G_arena.adj() += PG.val_op() * res.adj_op().transpose() - * tcrossprod(inv_G.val_op()); - G_arena.adj() - += crossprod(inv_G.val_op()) * res.adj_op().transpose() * GP.val_op(); - }); + arena_t inv_G((G_arena.val_op().transpose() * G_arena.val_op()) + .ldlt() + .solve(G_arena.val_op().transpose())); + reverse_pass_callback(internal::generalized_inverse_lambda(G_arena, inv_G)); return ret_type(inv_G); } } } // namespace math } // namespace stan - #endif diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index f7ceea1dcef..746826e3bf0 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -7,20 +7,14 @@ TEST(mathMixMatFun, ad_tests) { using stan::test::expect_ad_matvar; auto f = [](const auto& G) { return stan::math::generalized_inverse(G); }; - auto fjit - = [](const auto& G) { return stan::math::generalized_inverse(G, 0.0); }; - auto fjit2 - = [](const auto& G) { return stan::math::generalized_inverse(G, 1e-6); }; Eigen::MatrixXd t(0, 0); expect_ad(f, t); - expect_ad(fjit, t); expect_ad_matvar(f, t); Eigen::MatrixXd u(1, 1); u << 2; expect_ad(f, u); - expect_ad(fjit, u); expect_ad_matvar(f, u); Eigen::MatrixXd v(2, 3); @@ -30,9 +24,13 @@ TEST(mathMixMatFun, ad_tests) { v << 1.9, 1.3, 2.5, 0.4, 1.7, 0.1; expect_ad(f, v); - expect_ad(fjit, v); expect_ad_matvar(f, v); + Eigen::MatrixXd s(2, 4); + s << 3.4, 2, 5, 1.2, 2, 1, 3.2, 3.1; + expect_ad(f, s); + expect_ad_matvar(f, s); + // issues around zero require looser tolerances for hessians stan::test::ad_tolerances tols; tols.hessian_hessian_ = 2.0; @@ -41,32 +39,26 @@ TEST(mathMixMatFun, ad_tests) { Eigen::MatrixXd w(3, 4); w << 2, 3, 5, 7, 11, 13, 17, 19, 23, 25, 27, 29; expect_ad(tols, f, w); - expect_ad(tols, fjit, w); expect_ad_matvar(f, w); Eigen::MatrixXd z(2, 2); z << 1, 2, 5, std::numeric_limits::quiet_NaN(); EXPECT_NO_THROW(stan::math::generalized_inverse(z)); - EXPECT_NO_THROW(stan::math::generalized_inverse(z, 0.0)); // autodiff throws, so following fails (throw behavior must match to pass) Eigen::MatrixXd a(2, 2); a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); expect_ad(f, a); - expect_ad(fjit, a); expect_ad_matvar(f, a); - // singular matrix, should use the - // alias to input small amount of jitter on the diagonal - - Eigen::MatrixXd m(3, 2); - m << 1, 2, 2, 4, 1, 2; - EXPECT_THROW(stan::math::generalized_inverse(m), std::domain_error); - - // should work with jittered version - stan::test::ad_tolerances tols3; - tols3.hessian_hessian_ = 0.01; - tols3.hessian_fvar_hessian_ = 0.01; - expect_ad_matvar(tols3, fjit2, m); + /* + finite diffs blows up but the ad grad is good + is there a test for this? + Eigen::MatrixXd m(2, 3); + m << 1, 2, 1, 2, 4, 2; + expect_ad(f, m); + expect_ad_matvar(f, m); + EXPECT_NO_THROW(stan::math::generalized_inverse(m)); + */ } diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 477fde2eea5..2f393a180c8 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -9,26 +9,30 @@ TEST(MathMatrixPrim, Zero) { EXPECT_EQ(0, ginv.cols()); } -TEST(MathMatrixPrim, EqualJitter) { +TEST(MathMatrixPrim, Equal1) { using stan::math::generalized_inverse; - stan::math::matrix_d m1(3, 2); - m1 << 1, 2, 2, 4, 1, 2; + stan::math::matrix_d m1(2, 3); + m1 << 1, 3, 5, 2, 4, 6; + + stan::math::matrix_d m2(2, 2); + m2 << 1, 0, 0, 1; - stan::math::matrix_d m2(2, 3); - m2 << 1 / 30.0, 1 / 15.0, 1 / 30.0, 1 / 15.0, 2 / 15.0, 1 / 15.0; - stan::math::matrix_d m3 = generalized_inverse(m1, 1e-8); + stan::math::matrix_d m3 = m1 * generalized_inverse(m1); EXPECT_MATRIX_NEAR(m2, m3, 1e-9); } -TEST(MathMatrixPrim, Equal) { +TEST(MathMatrixPrim, Equal2) { using stan::math::generalized_inverse; - stan::math::matrix_d m1(2, 3); - m1 << 1, 3, 5, 2, 4, 6; + stan::math::matrix_d m1(3, 2); + m1 << 1, 2, 2, 4, 1, 2; - stan::math::matrix_d m2(3, 2); - m2 << -4.0 / 3.0, 13.0 / 12.0, -1.0 / 3.0, 1.0 / 3.0, 2.0 / 3.0, -5.0 / 12.0; - stan::math::matrix_d m3 = generalized_inverse(m1); + stan::math::matrix_d m2(3, 3); + m2 << 1.0 / 6.0, 1.0 / 3.0, 1.0 / 6.0, + 1.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0, + 1.0 / 6.0, 1.0 / 3.0, 1.0 / 6.0; + + stan::math::matrix_d m3 = m1 * generalized_inverse(m1); EXPECT_MATRIX_NEAR(m2, m3, 1e-9); -} +} \ No newline at end of file From 18fb94fd51868157a42650e6bf7e5b0c9225e26e Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 17 Dec 2020 12:32:54 +0000 Subject: [PATCH 096/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/generalized_inverse.hpp | 9 +---- stan/math/rev/fun/generalized_inverse.hpp | 37 +++++++++---------- .../math/mix/fun/generalized_inverse_test.cpp | 20 +++++----- .../prim/fun/generalized_inverse_test.cpp | 7 ++-- 4 files changed, 33 insertions(+), 40 deletions(-) diff --git a/stan/math/prim/fun/generalized_inverse.hpp b/stan/math/prim/fun/generalized_inverse.hpp index 8eba7815ca4..0f6e1bb42bb 100644 --- a/stan/math/prim/fun/generalized_inverse.hpp +++ b/stan/math/prim/fun/generalized_inverse.hpp @@ -42,14 +42,9 @@ generalized_inverse(const EigMat& G) { const auto& G_ref = to_ref(G); if (G.rows() < G.cols()) { - return (G_ref * G_ref.transpose()) - .ldlt() - .solve(G_ref) - .transpose(); + return (G_ref * G_ref.transpose()).ldlt().solve(G_ref).transpose(); } else { - return (G_ref.transpose() * G_ref) - .ldlt() - .solve(G_ref.transpose()); + return (G_ref.transpose() * G_ref).ldlt().solve(G_ref.transpose()); } } diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 54e2bfb2d20..92b7d9a7179 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -20,20 +20,19 @@ namespace internal { */ template inline auto generalized_inverse_lambda(T1& G_arena, T2& inv_G) { -return [G_arena, inv_G]() mutable { - G_arena.adj() - += -(inv_G.val_op().transpose() * inv_G.adj_op() - * inv_G.val_op().transpose()) - + (-G_arena.val_op() * inv_G.val_op() - + Eigen::MatrixXd::Identity(G_arena.rows(), inv_G.cols())) - * inv_G.adj_op().transpose() * inv_G.val_op() - * inv_G.val_op().transpose() - + inv_G.val_op().transpose() * inv_G.val_op() - * inv_G.adj_op().transpose() - * (-inv_G.val_op() * G_arena.val_op() - + Eigen::MatrixXd::Identity(inv_G.rows(), - G_arena.cols())); - }; + return [G_arena, inv_G]() mutable { + G_arena.adj() + += -(inv_G.val_op().transpose() * inv_G.adj_op() + * inv_G.val_op().transpose()) + + (-G_arena.val_op() * inv_G.val_op() + + Eigen::MatrixXd::Identity(G_arena.rows(), inv_G.cols())) + * inv_G.adj_op().transpose() * inv_G.val_op() + * inv_G.val_op().transpose() + + inv_G.val_op().transpose() * inv_G.val_op() + * inv_G.adj_op().transpose() + * (-inv_G.val_op() * G_arena.val_op() + + Eigen::MatrixXd::Identity(inv_G.rows(), G_arena.cols())); + }; } } // namespace internal @@ -72,16 +71,16 @@ inline auto generalized_inverse(const VarMat& G) { if (G.rows() < G.cols()) { arena_t G_arena(G); arena_t inv_G((G_arena.val_op() * G_arena.val_op().transpose()) - .ldlt() - .solve(G_arena.val_op()) - .transpose()); + .ldlt() + .solve(G_arena.val_op()) + .transpose()); reverse_pass_callback(internal::generalized_inverse_lambda(G_arena, inv_G)); return ret_type(inv_G); } else { arena_t G_arena(G); arena_t inv_G((G_arena.val_op().transpose() * G_arena.val_op()) - .ldlt() - .solve(G_arena.val_op().transpose())); + .ldlt() + .solve(G_arena.val_op().transpose())); reverse_pass_callback(internal::generalized_inverse_lambda(G_arena, inv_G)); return ret_type(inv_G); } diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 746826e3bf0..05be0df7c6d 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -27,7 +27,7 @@ TEST(mathMixMatFun, ad_tests) { expect_ad_matvar(f, v); Eigen::MatrixXd s(2, 4); - s << 3.4, 2, 5, 1.2, 2, 1, 3.2, 3.1; + s << 3.4, 2, 5, 1.2, 2, 1, 3.2, 3.1; expect_ad(f, s); expect_ad_matvar(f, s); @@ -52,13 +52,13 @@ TEST(mathMixMatFun, ad_tests) { expect_ad(f, a); expect_ad_matvar(f, a); - /* - finite diffs blows up but the ad grad is good - is there a test for this? - Eigen::MatrixXd m(2, 3); - m << 1, 2, 1, 2, 4, 2; - expect_ad(f, m); - expect_ad_matvar(f, m); - EXPECT_NO_THROW(stan::math::generalized_inverse(m)); - */ + /* + finite diffs blows up but the ad grad is good + is there a test for this? + Eigen::MatrixXd m(2, 3); + m << 1, 2, 1, 2, 4, 2; + expect_ad(f, m); + expect_ad_matvar(f, m); + EXPECT_NO_THROW(stan::math::generalized_inverse(m)); + */ } diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 2f393a180c8..4fcb804a46f 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -29,10 +29,9 @@ TEST(MathMatrixPrim, Equal2) { m1 << 1, 2, 2, 4, 1, 2; stan::math::matrix_d m2(3, 3); - m2 << 1.0 / 6.0, 1.0 / 3.0, 1.0 / 6.0, - 1.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0, - 1.0 / 6.0, 1.0 / 3.0, 1.0 / 6.0; + m2 << 1.0 / 6.0, 1.0 / 3.0, 1.0 / 6.0, 1.0 / 3.0, 2.0 / 3.0, 1.0 / 3.0, + 1.0 / 6.0, 1.0 / 3.0, 1.0 / 6.0; stan::math::matrix_d m3 = m1 * generalized_inverse(m1); EXPECT_MATRIX_NEAR(m2, m3, 1e-9); -} \ No newline at end of file +} \ No newline at end of file From 7fe4f9cf510862619a766d95c3656cdfc754526d Mon Sep 17 00:00:00 2001 From: spinkney Date: Thu, 17 Dec 2020 07:34:38 -0500 Subject: [PATCH 097/102] newline for test --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 05be0df7c6d..28654affbd6 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -62,3 +62,4 @@ TEST(mathMixMatFun, ad_tests) { EXPECT_NO_THROW(stan::math::generalized_inverse(m)); */ } + From 866e9ea3423d5d3929d4eca06b80c8387181bb03 Mon Sep 17 00:00:00 2001 From: spinkney Date: Thu, 17 Dec 2020 07:35:28 -0500 Subject: [PATCH 098/102] new line for test --- test/unit/math/prim/fun/generalized_inverse_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 4fcb804a46f..3f69d575a5a 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -34,4 +34,4 @@ TEST(MathMatrixPrim, Equal2) { stan::math::matrix_d m3 = m1 * generalized_inverse(m1); EXPECT_MATRIX_NEAR(m2, m3, 1e-9); -} \ No newline at end of file +} From 18178077ced9751e5d88fa4e6b46bc4aaed56cc5 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Thu, 17 Dec 2020 12:36:10 +0000 Subject: [PATCH 099/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- test/unit/math/mix/fun/generalized_inverse_test.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 28654affbd6..05be0df7c6d 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -62,4 +62,3 @@ TEST(mathMixMatFun, ad_tests) { EXPECT_NO_THROW(stan::math::generalized_inverse(m)); */ } - From c80f21212ce5f112021cc8b16bf69a28d71c7ef0 Mon Sep 17 00:00:00 2001 From: spinkney Date: Fri, 18 Dec 2020 06:48:50 -0500 Subject: [PATCH 100/102] Removed unnecessary test --- stan/math/rev/fun/generalized_inverse.hpp | 4 +++- test/unit/math/mix/fun/generalized_inverse_test.cpp | 10 ---------- 2 files changed, 3 insertions(+), 11 deletions(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index 92b7d9a7179..e10525384d4 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -44,7 +44,9 @@ inline auto generalized_inverse_lambda(T1& G_arena, T2& inv_G) { * @return Generalized inverse of the matrix (an empty matrix if the specified * matrix has size zero). * - * @note Reverse mode differentiation algorithm reference: + * @note For the derivatives of this function to exist the matrix must be + * of constant rank. + * Reverse mode differentiation algorithm reference: * *
  • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses * and Nonlinear Least Squares Problems Whose Variables Separate. SIAM diff --git a/test/unit/math/mix/fun/generalized_inverse_test.cpp b/test/unit/math/mix/fun/generalized_inverse_test.cpp index 05be0df7c6d..f079f1366eb 100644 --- a/test/unit/math/mix/fun/generalized_inverse_test.cpp +++ b/test/unit/math/mix/fun/generalized_inverse_test.cpp @@ -51,14 +51,4 @@ TEST(mathMixMatFun, ad_tests) { a << 1.9, 0.3, 0.3, std::numeric_limits::infinity(); expect_ad(f, a); expect_ad_matvar(f, a); - - /* - finite diffs blows up but the ad grad is good - is there a test for this? - Eigen::MatrixXd m(2, 3); - m << 1, 2, 1, 2, 4, 2; - expect_ad(f, m); - expect_ad_matvar(f, m); - EXPECT_NO_THROW(stan::math::generalized_inverse(m)); - */ } From a8b1287163ad874a3e2468f21f65db8fa7434fa6 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Fri, 18 Dec 2020 11:55:50 +0000 Subject: [PATCH 101/102] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/generalized_inverse.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/stan/math/rev/fun/generalized_inverse.hpp b/stan/math/rev/fun/generalized_inverse.hpp index e10525384d4..9ec363dd5d5 100644 --- a/stan/math/rev/fun/generalized_inverse.hpp +++ b/stan/math/rev/fun/generalized_inverse.hpp @@ -45,7 +45,7 @@ inline auto generalized_inverse_lambda(T1& G_arena, T2& inv_G) { * matrix has size zero). * * @note For the derivatives of this function to exist the matrix must be - * of constant rank. + * of constant rank. * Reverse mode differentiation algorithm reference: * *
    • Golub, G.H. and Pereyra, V. The Differentiation of Pseudo-Inverses From 20742a87ee767ca1456b959d32397b458f05045b Mon Sep 17 00:00:00 2001 From: spinkney Date: Fri, 18 Dec 2020 12:16:26 -0500 Subject: [PATCH 102/102] prim test for singular --- test/unit/math/prim/fun/generalized_inverse_test.cpp | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/test/unit/math/prim/fun/generalized_inverse_test.cpp b/test/unit/math/prim/fun/generalized_inverse_test.cpp index 3f69d575a5a..5f3bd6a7170 100644 --- a/test/unit/math/prim/fun/generalized_inverse_test.cpp +++ b/test/unit/math/prim/fun/generalized_inverse_test.cpp @@ -9,6 +9,16 @@ TEST(MathMatrixPrim, Zero) { EXPECT_EQ(0, ginv.cols()); } +TEST(MathMatrixPrim, Singular) { + using stan::math::generalized_inverse; + + stan::math::matrix_d m1(2, 3); + m1 << 1, 2, 1, 2, 4, 2; + + stan::math::matrix_d m2 = m1 * generalized_inverse(m1) * m1; + EXPECT_MATRIX_NEAR(m1, m2, 1e-9); +} + TEST(MathMatrixPrim, Equal1) { using stan::math::generalized_inverse;