From 63237eaf78abf0fa29485e9a247e7f8826eac5e2 Mon Sep 17 00:00:00 2001 From: hyunjimoon Date: Mon, 5 Oct 2020 22:53:10 +0900 Subject: [PATCH 01/32] implement interp1 --- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/interp1.hpp | 45 ++++++++++++++++++++++++ test/unit/math/prim/fun/interp1_test.cpp | 27 ++++++++++++++ 3 files changed, 73 insertions(+) create mode 100644 stan/math/prim/fun/interp1.hpp create mode 100644 test/unit/math/prim/fun/interp1_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index bdc196f6e22..aa49846f070 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -128,6 +128,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/interp1.hpp b/stan/math/prim/fun/interp1.hpp new file mode 100644 index 00000000000..92412a5bb98 --- /dev/null +++ b/stan/math/prim/fun/interp1.hpp @@ -0,0 +1,45 @@ +#ifndef STAN_MATH_PRIM_FUN_INTERP1_HPP +#define STAN_MATH_PRIM_FUN_INTERP1_HPP + +#include +using namespace std; +namespace stan { + namespace math { + /** + * Return the one-dimensional interpolation of the second argument of points corresponding to the first argument. + * @tparam T_true type of the true argument + * @tparam T_false type of the false argument + * @param c Boolean condition value. + * @param y_true Value to return if condition is true. + * @param y_false Value to return if condition is false. + */ + inline Eigen::VectorXd + interp1(const Eigen::VectorXd& xData, const Eigen::VectorXd& yData, const Eigen::VectorXd& xVector) + { + + Eigen::VectorXd yVals; + int size_x = xData.size(); + int size_xData = xData.size(); + + for(int j=0; j < size_x; j++) { + double x = xVector[j]; + int i = 0; + if (x >= xData[size_xData - 2]) // special case: beyond right end + { + i = size_xData - 2; + } else { + while (x > xData[i + 1]) + i++; // find left end of interval for interpolation + } + double xL = xData[i], yL = yData[i], xR = xData[i + 1], yR = yData[i + 1]; // points on either side (unless beyond ends) + if (x < xL) yR = yL; // if beyond ends of array + if (x > xR) yL = yR; + double dydx = (yR - yL) / (xR - xL); // gradient + yVals[j] = yL + dydx * (x - xL); + } + return yVals; // linear interpolation + } + } // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/prim/fun/interp1_test.cpp b/test/unit/math/prim/fun/interp1_test.cpp new file mode 100644 index 00000000000..cb00c3c1343 --- /dev/null +++ b/test/unit/math/prim/fun/interp1_test.cpp @@ -0,0 +1,27 @@ +#include +#include +#include +#include + +// This file goes in test/unit/math/prim/fun/ +// +// To run the test type: +// ./runTest.py test/unit/math/prim/fun/interp1_test.cpp + +TEST(MathPrim, interp1_test) { +using stan::test::expect_near_rel; +Eigen::VectorXd x(5); +x << 0.00, 1.25, 2.50, 3.75, 5.00; +Eigen::VectorXd y(5); +y << 0.0000000000, 0.9489846194, 0.5984721441, -0.5715613187, -0.9589242747; + +Eigen::VectorXd x_test(4); +x_test << 1.0, 1.7, 2.1, 4.0; +Eigen::VectorXd y_ref(4); // Reference values from R +y_ref << 0.7591876955, 0.8228001283, 0.7106361362, -0.6490339099; + +Eigen::VectorXd y_test = stan::math::interp1(x, y, x_test); + +stan::test::expect_near_rel("test Interp1", y_test, y_ref); +cout << y_test; +} From 71fe4b6c38de35d7eec3f2de9920d7f8a6b4f305 Mon Sep 17 00:00:00 2001 From: hyunjimoon Date: Sat, 21 Nov 2020 21:45:31 +0900 Subject: [PATCH 02/32] reverse svd and its test --- stan/math/rev/fun/singular_values.hpp | 70 +++++++++++++++++++ .../unit/math/rev/fun/singularvalues_test.cpp | 25 +++++++ 2 files changed, 95 insertions(+) create mode 100644 stan/math/rev/fun/singular_values.hpp create mode 100644 test/unit/math/rev/fun/singularvalues_test.cpp diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp new file mode 100644 index 00000000000..77a55a6ca7d --- /dev/null +++ b/stan/math/rev/fun/singular_values.hpp @@ -0,0 +1,70 @@ +#ifndef MATH_SINGULAR_VALUES_H +#define MATH_SINGULAR_VALUES_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { + namespace math { + +/** + * Return the eigenvalues of the specified symmetric matrix. + *

See eigen_decompose() for more information. + * + * @tparam EigMat type of input matrix + * @param m Input matrix. + * @return Eigenvalues of matrix. + */ + template * = nullptr, require_not_vt_var* = nullptr> + inline auto singular_values(const EigMat& m) { + + using ret_type = promote_scalar_t; + check_nonzero_size("singular_values", "m", m); + const int M = m.rows(); + auto arena_m = to_arena(m); + Eigen::MatrixXd m_val = value_of(arena_m); + Eigen::JacobiSVD svd(m_val, Eigen::ComputeFullU | Eigen::ComputeFullV); + Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> D = Eigen::MatrixXd::Zero(M,M); + D.diagonal() = svd.singularValues(); + arena_t sigularvalues = svd.singularValues(); + auto U = svd.matrixU(); + auto V = svd.matrixU(); + // equation (4) from "Diļ¬€erentiable Programming Tensor Networks" + Eigen::MatrixXd Fp(M,M); + Eigen::MatrixXd Fm(M,M); + + for(int i = 0; i < M; i++) { + for(int j = 0; i < M; j++) { + if(j == i) { + Fp[i, j] = 0.0; + Fm[i, j] = 0.0; + } else { + Fp[i, j] = 1.0 / (D[j,j] - D[i,i]) + 1.0 / (D[i,i] + D[j,j]); + Fm[i, j] = 1.0 / (D[j,j] - D[i,i]) - 1.0 / (D[i,i] + D[j,j]); + } + } + } + reverse_pass_callback( + [arena_m, U, V]() mutable { + arena_m.adj() += 0.5 * U * (Fp * (U.transpose() * U.adj() - t(U.adj()) * U)) * V.transpose() // adjU contributions + +(Eigen::MatrixXd::Identity(M, M) - U * U.transpose()) * U.adj() * D.inverse() * V.transpose() + + U * D.adj() * V.transpose() // adjD contributions + + 0.5 * U * (Fm * (V.transpose() * V.adj() - V.adj().transpose() * V)) * V.transpose() // adjV contributions + + U * D.inverse() * V.adj().transpose() * (Eigen::MatrixXd::Identity(M, M) - V * V.transpose()); + }); + return ret_type(sigularvalues); + } + } // namespace math +} // namespace stan + +#endif \ No newline at end of file diff --git a/test/unit/math/rev/fun/singularvalues_test.cpp b/test/unit/math/rev/fun/singularvalues_test.cpp new file mode 100644 index 00000000000..2981bb3c63f --- /dev/null +++ b/test/unit/math/rev/fun/singularvalues_test.cpp @@ -0,0 +1,25 @@ +#include +#include +#include +#include +#include + +TEST(AgradRev, singularvalues_gradient) { +// logdet(A) can be calculated using singularvalues of matrix A +// the derivative of logdet(A) should be inverse(A) + +Eigen::MatrixXd a(4, 4); +// Random symmetric matrix +a << 1.8904, 0.7204, -0.1599, 1.2028, 0.7204, 7.3394, 2.0895, -0.6151, +-0.1599, 2.0895, 2.4601, -1.7219, 1.2028, -0.6151, -1.7219, 4.5260; +stan::math::matrix_d a_inv = stan::math::inverse(a); + +stan::math::matrix_v a_v(a); +auto w = singular_values(a_v); +auto logdet = stan::math::sum(stan::math::log(w)); + +stan::math::set_zero_all_adjoints(); +logdet.grad(); +std::cout<<"a_inv.val()" < Date: Sat, 21 Nov 2020 21:58:49 +0900 Subject: [PATCH 03/32] clean interp1 and add to function header file --- stan/math/prim/fun.hpp | 1 - stan/math/prim/fun/interp1.hpp | 45 ------------------------ stan/math/rev/fun.hpp | 1 + test/unit/math/prim/fun/interp1_test.cpp | 27 -------------- 4 files changed, 1 insertion(+), 73 deletions(-) delete mode 100644 stan/math/prim/fun/interp1.hpp delete mode 100644 test/unit/math/prim/fun/interp1_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index aa49846f070..bdc196f6e22 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -128,7 +128,6 @@ #include #include #include -#include #include #include #include diff --git a/stan/math/prim/fun/interp1.hpp b/stan/math/prim/fun/interp1.hpp deleted file mode 100644 index 92412a5bb98..00000000000 --- a/stan/math/prim/fun/interp1.hpp +++ /dev/null @@ -1,45 +0,0 @@ -#ifndef STAN_MATH_PRIM_FUN_INTERP1_HPP -#define STAN_MATH_PRIM_FUN_INTERP1_HPP - -#include -using namespace std; -namespace stan { - namespace math { - /** - * Return the one-dimensional interpolation of the second argument of points corresponding to the first argument. - * @tparam T_true type of the true argument - * @tparam T_false type of the false argument - * @param c Boolean condition value. - * @param y_true Value to return if condition is true. - * @param y_false Value to return if condition is false. - */ - inline Eigen::VectorXd - interp1(const Eigen::VectorXd& xData, const Eigen::VectorXd& yData, const Eigen::VectorXd& xVector) - { - - Eigen::VectorXd yVals; - int size_x = xData.size(); - int size_xData = xData.size(); - - for(int j=0; j < size_x; j++) { - double x = xVector[j]; - int i = 0; - if (x >= xData[size_xData - 2]) // special case: beyond right end - { - i = size_xData - 2; - } else { - while (x > xData[i + 1]) - i++; // find left end of interval for interpolation - } - double xL = xData[i], yL = yData[i], xR = xData[i + 1], yR = yData[i + 1]; // points on either side (unless beyond ends) - if (x < xL) yR = yL; // if beyond ends of array - if (x > xR) yL = yR; - double dydx = (yR - yL) / (xR - xL); // gradient - yVals[j] = yL + dydx * (x - xL); - } - return yVals; // linear interpolation - } - } // namespace math -} // namespace stan - -#endif diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index ebc92237dd4..9fa75220888 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -128,6 +128,7 @@ #include #include #include +#include #include #include #include diff --git a/test/unit/math/prim/fun/interp1_test.cpp b/test/unit/math/prim/fun/interp1_test.cpp deleted file mode 100644 index cb00c3c1343..00000000000 --- a/test/unit/math/prim/fun/interp1_test.cpp +++ /dev/null @@ -1,27 +0,0 @@ -#include -#include -#include -#include - -// This file goes in test/unit/math/prim/fun/ -// -// To run the test type: -// ./runTest.py test/unit/math/prim/fun/interp1_test.cpp - -TEST(MathPrim, interp1_test) { -using stan::test::expect_near_rel; -Eigen::VectorXd x(5); -x << 0.00, 1.25, 2.50, 3.75, 5.00; -Eigen::VectorXd y(5); -y << 0.0000000000, 0.9489846194, 0.5984721441, -0.5715613187, -0.9589242747; - -Eigen::VectorXd x_test(4); -x_test << 1.0, 1.7, 2.1, 4.0; -Eigen::VectorXd y_ref(4); // Reference values from R -y_ref << 0.7591876955, 0.8228001283, 0.7106361362, -0.6490339099; - -Eigen::VectorXd y_test = stan::math::interp1(x, y, x_test); - -stan::test::expect_near_rel("test Interp1", y_test, y_ref); -cout << y_test; -} From 77178abdff63197f6e2308d02e8a11810bf5a51a Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 22 Nov 2020 14:52:57 -0500 Subject: [PATCH 04/32] Formatting and such (Issue #2192) --- stan/math/prim/fun/singular_values.hpp | 4 +- stan/math/rev/fun/singular_values.hpp | 93 ++++++++++++++------------ 2 files changed, 53 insertions(+), 44 deletions(-) diff --git a/stan/math/prim/fun/singular_values.hpp b/stan/math/prim/fun/singular_values.hpp index cbb89938c49..85035351e29 100644 --- a/stan/math/prim/fun/singular_values.hpp +++ b/stan/math/prim/fun/singular_values.hpp @@ -17,7 +17,9 @@ namespace math { * @param m Specified matrix. * @return Singular values of the matrix. */ -template * = nullptr> + template * = nullptr, + require_not_vt_var* = nullptr> Eigen::Matrix, Eigen::Dynamic, 1> singular_values( const EigMat& m) { if (m.size() == 0) { diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index 77a55a6ca7d..5c5888051cf 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -15,56 +15,63 @@ #include namespace stan { - namespace math { +namespace math { /** - * Return the eigenvalues of the specified symmetric matrix. - *

See eigen_decompose() for more information. + * Return the singular values of the specified symmetric matrix. + * + * For MxN input matrix, M must be <= N * * @tparam EigMat type of input matrix - * @param m Input matrix. - * @return Eigenvalues of matrix. + * @param m MxN input matrix, M <= N + * @return Singular values of matrix */ - template * = nullptr, require_not_vt_var* = nullptr> - inline auto singular_values(const EigMat& m) { +template * = nullptr, + require_not_vt_var* = nullptr> +inline auto singular_values(const EigMat& m) { + using ret_type = promote_scalar_t; + check_nonzero_size("singular_values", "m", m); + + const int M = m.rows(); + auto arena_m = to_arena(m); + + Eigen::MatrixXd m_val = value_of(arena_m); + Eigen::JacobiSVD svd(m_val, Eigen::ComputeFullU | Eigen::ComputeFullV); + + Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> D = Eigen::MatrixXd::Zero(M,M); + D.diagonal() = svd.singularValues(); + arena_t sigular_values = svd.singularValues(); - using ret_type = promote_scalar_t; - check_nonzero_size("singular_values", "m", m); - const int M = m.rows(); - auto arena_m = to_arena(m); - Eigen::MatrixXd m_val = value_of(arena_m); - Eigen::JacobiSVD svd(m_val, Eigen::ComputeFullU | Eigen::ComputeFullV); - Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> D = Eigen::MatrixXd::Zero(M,M); - D.diagonal() = svd.singularValues(); - arena_t sigularvalues = svd.singularValues(); - auto U = svd.matrixU(); - auto V = svd.matrixU(); - // equation (4) from "Diļ¬€erentiable Programming Tensor Networks" - Eigen::MatrixXd Fp(M,M); - Eigen::MatrixXd Fm(M,M); + auto U = svd.matrixU(); + auto V = svd.matrixU(); - for(int i = 0; i < M; i++) { - for(int j = 0; i < M; j++) { - if(j == i) { - Fp[i, j] = 0.0; - Fm[i, j] = 0.0; - } else { - Fp[i, j] = 1.0 / (D[j,j] - D[i,i]) + 1.0 / (D[i,i] + D[j,j]); - Fm[i, j] = 1.0 / (D[j,j] - D[i,i]) - 1.0 / (D[i,i] + D[j,j]); - } - } - } - reverse_pass_callback( - [arena_m, U, V]() mutable { - arena_m.adj() += 0.5 * U * (Fp * (U.transpose() * U.adj() - t(U.adj()) * U)) * V.transpose() // adjU contributions - +(Eigen::MatrixXd::Identity(M, M) - U * U.transpose()) * U.adj() * D.inverse() * V.transpose() - + U * D.adj() * V.transpose() // adjD contributions - + 0.5 * U * (Fm * (V.transpose() * V.adj() - V.adj().transpose() * V)) * V.transpose() // adjV contributions - + U * D.inverse() * V.adj().transpose() * (Eigen::MatrixXd::Identity(M, M) - V * V.transpose()); - }); - return ret_type(sigularvalues); + // equation (4) from "Differentiable Programming Tensor Networks" + Eigen::MatrixXd Fp(M,M); + Eigen::MatrixXd Fm(M,M); + + for(int i = 0; i < M; i++) { + for(int j = 0; i < M; j++) { + if(j == i) { + Fp[i, j] = 0.0; + Fm[i, j] = 0.0; + } else { + Fp[i, j] = 1.0 / (D[j,j] - D[i,i]) + 1.0 / (D[i,i] + D[j,j]); + Fm[i, j] = 1.0 / (D[j,j] - D[i,i]) - 1.0 / (D[i,i] + D[j,j]); } - } // namespace math + } + } + + reverse_pass_callback([arena_m, U, V]() mutable { + arena_m.adj() += 0.5 * U * (Fp * (U.transpose() * U.adj() - t(U.adj()) * U)) * V.transpose() // adjU contributions + +(Eigen::MatrixXd::Identity(M, M) - U * U.transpose()) * U.adj() * D.inverse() * V.transpose() + + U * D.adj() * V.transpose() // adjD contributions + + 0.5 * U * (Fm * (V.transpose() * V.adj() - V.adj().transpose() * V)) * V.transpose() // adjV contributions + + U * D.inverse() * V.adj().transpose() * (Eigen::MatrixXd::Identity(M, M) - V * V.transpose()); + }); + return ret_type(sigularvalues); +} +} // namespace math } // namespace stan -#endif \ No newline at end of file +#endif From dc43ce0fb2d229d32bd0bfdbe4055657f308d467 Mon Sep 17 00:00:00 2001 From: hyunjimoon Date: Sun, 29 Nov 2020 14:03:50 +0900 Subject: [PATCH 05/32] update reverse mode on deleting D and add to_arena --- stan/math/rev/fun/singular_values.hpp | 56 +++++++++++++-------------- 1 file changed, 27 insertions(+), 29 deletions(-) diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index 5c5888051cf..49c46cd67e8 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -18,58 +18,56 @@ namespace stan { namespace math { /** - * Return the singular values of the specified symmetric matrix. + * Return the singular values of the specified matrix. * - * For MxN input matrix, M must be <= N + * Equation (4) from Differentiable Programming Tensor Networks + * is used for MxN input matrix's adjoint calculation. * * @tparam EigMat type of input matrix - * @param m MxN input matrix, M <= N + * @param m MxN input matrix * @return Singular values of matrix */ template * = nullptr, - require_not_vt_var* = nullptr> + require_eigen_matrix_dynamic_t* = nullptr, + require_vt_var* = nullptr> inline auto singular_values(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("singular_values", "m", m); - const int M = m.rows(); + const int M = to_arena(m.rows()); auto arena_m = to_arena(m); Eigen::MatrixXd m_val = value_of(arena_m); - Eigen::JacobiSVD svd(m_val, Eigen::ComputeFullU | Eigen::ComputeFullV); + Eigen::JacobiSVD svd(m_val, Eigen::ComputeFullU | Eigen::ComputeFullV); - Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> D = Eigen::MatrixXd::Zero(M,M); - D.diagonal() = svd.singularValues(); - arena_t sigular_values = svd.singularValues(); - - auto U = svd.matrixU(); - auto V = svd.matrixU(); - - // equation (4) from "Differentiable Programming Tensor Networks" + auto singular_values = to_arena(svd.singularValues()); + auto U = to_arena(svd.matrixU()); + auto V = to_arena(svd.matrixV()); Eigen::MatrixXd Fp(M,M); Eigen::MatrixXd Fm(M,M); - for(int i = 0; i < M; i++) { - for(int j = 0; i < M; j++) { + for(int j = 0; j < M; j++) { if(j == i) { - Fp[i, j] = 0.0; - Fm[i, j] = 0.0; + Fp(i, j) = 0.0; + Fm(i, j) = 0.0; } else { - Fp[i, j] = 1.0 / (D[j,j] - D[i,i]) + 1.0 / (D[i,i] + D[j,j]); - Fm[i, j] = 1.0 / (D[j,j] - D[i,i]) - 1.0 / (D[i,i] + D[j,j]); + Fp(i, j) = 1.0 / (singular_values[j] - singular_values[i]) + 1.0 / (singular_values[i] + singular_values[j]); + Fm(i, j) = 1.0 / (singular_values[j] - singular_values[i]) - 1.0 / (singular_values[j] + singular_values[i]); } } } + auto arena_Fp = to_arena(Fp); + auto arena_Fm = to_arena(Fm); - reverse_pass_callback([arena_m, U, V]() mutable { - arena_m.adj() += 0.5 * U * (Fp * (U.transpose() * U.adj() - t(U.adj()) * U)) * V.transpose() // adjU contributions - +(Eigen::MatrixXd::Identity(M, M) - U * U.transpose()) * U.adj() * D.inverse() * V.transpose() - + U * D.adj() * V.transpose() // adjD contributions - + 0.5 * U * (Fm * (V.transpose() * V.adj() - V.adj().transpose() * V)) * V.transpose() // adjV contributions - + U * D.inverse() * V.adj().transpose() * (Eigen::MatrixXd::Identity(M, M) - V * V.transpose()); - }); - return ret_type(sigularvalues); + reverse_pass_callback( + [arena_m, M, arena_Fp, arena_Fm, U, singular_values, V]() mutable { + arena_m.adj() += 0.5 * U * (arena_Fp * (U.transpose() * U.adj() - U.adj().transpose() * U)) * V.transpose() // adjU contributions + +(Eigen::MatrixXd::Identity(M, M) - U * U.transpose()) * U.adj() * singular_values.asDiagonal().inverse() * V.transpose() + + U * singular_values.adj().asDiagonal() * V.transpose() // adjD contributions + + 0.5 * U * (arena_Fm * (V.transpose() * V.adj() - V.adj().transpose() * V)) * V.transpose() // adjV contributions + + U * singular_values.asDiagonal().inverse() * V.adj().transpose() * (Eigen::MatrixXd::Identity(M, M) - V * V.transpose()); + }); + return ret_type(singular_values); } } // namespace math } // namespace stan From c596c4d7ba44efa1716c7b4684ae20ce5505988e Mon Sep 17 00:00:00 2001 From: Dashadower Date: Thu, 10 Dec 2020 19:26:28 +0900 Subject: [PATCH 06/32] Finish singular_values.hpp, verify tests pass --- stan/math/rev/fun/singular_values.hpp | 20 ++++++++++--------- .../unit/math/rev/fun/singularvalues_test.cpp | 2 +- 2 files changed, 12 insertions(+), 10 deletions(-) diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index 49c46cd67e8..ee620bc62e7 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -40,7 +40,9 @@ inline auto singular_values(const EigMat& m) { Eigen::MatrixXd m_val = value_of(arena_m); Eigen::JacobiSVD svd(m_val, Eigen::ComputeFullU | Eigen::ComputeFullV); - auto singular_values = to_arena(svd.singularValues()); + vector_v singular_values(to_arena(svd.singularValues())); // switch to vector_v to enable adjoint + + //auto singular_values = to_arena(svd.singularValues()); auto U = to_arena(svd.matrixU()); auto V = to_arena(svd.matrixV()); Eigen::MatrixXd Fp(M,M); @@ -51,8 +53,8 @@ inline auto singular_values(const EigMat& m) { Fp(i, j) = 0.0; Fm(i, j) = 0.0; } else { - Fp(i, j) = 1.0 / (singular_values[j] - singular_values[i]) + 1.0 / (singular_values[i] + singular_values[j]); - Fm(i, j) = 1.0 / (singular_values[j] - singular_values[i]) - 1.0 / (singular_values[j] + singular_values[i]); + Fp(i, j) = 1.0 / (singular_values[j].val() - singular_values[i].val()) + 1.0 / (singular_values[i].val() + singular_values[j].val()); + Fm(i, j) = 1.0 / (singular_values[j].val() - singular_values[i].val()) - 1.0 / (singular_values[j].val() + singular_values[i].val()); } } } @@ -60,12 +62,12 @@ inline auto singular_values(const EigMat& m) { auto arena_Fm = to_arena(Fm); reverse_pass_callback( - [arena_m, M, arena_Fp, arena_Fm, U, singular_values, V]() mutable { - arena_m.adj() += 0.5 * U * (arena_Fp * (U.transpose() * U.adj() - U.adj().transpose() * U)) * V.transpose() // adjU contributions - +(Eigen::MatrixXd::Identity(M, M) - U * U.transpose()) * U.adj() * singular_values.asDiagonal().inverse() * V.transpose() - + U * singular_values.adj().asDiagonal() * V.transpose() // adjD contributions - + 0.5 * U * (arena_Fm * (V.transpose() * V.adj() - V.adj().transpose() * V)) * V.transpose() // adjV contributions - + U * singular_values.asDiagonal().inverse() * V.adj().transpose() * (Eigen::MatrixXd::Identity(M, M) - V * V.transpose()); + [arena_m, U, singular_values, V]() mutable { + arena_m.adj() //+= 0.5 * U * (arena_Fp * (U.transpose() * U.adj() - U.adj().transpose() * U)) * V.transpose() // adjU contributions + //+(Eigen::MatrixXd::Identity(M, M) - U * U.transpose()) * U.adj() * singular_values.asDiagonal().inverse() * V.transpose(); + += U * singular_values.adj().asDiagonal() * V.transpose(); // adjD contributions + //+ 0.5 * U * (arena_Fm * (V.transpose() * V.adj() - V.adj().transpose() * V)) * V.transpose() // adjV contributions + //+ U * singular_values.asDiagonal().inverse() * V.adj().transpose() * (Eigen::MatrixXd::Identity(M, M) - V * V.transpose()); }); return ret_type(singular_values); } diff --git a/test/unit/math/rev/fun/singularvalues_test.cpp b/test/unit/math/rev/fun/singularvalues_test.cpp index 2981bb3c63f..1d63699a043 100644 --- a/test/unit/math/rev/fun/singularvalues_test.cpp +++ b/test/unit/math/rev/fun/singularvalues_test.cpp @@ -15,7 +15,7 @@ a << 1.8904, 0.7204, -0.1599, 1.2028, 0.7204, 7.3394, 2.0895, -0.6151, stan::math::matrix_d a_inv = stan::math::inverse(a); stan::math::matrix_v a_v(a); -auto w = singular_values(a_v); +auto w = stan::math::singular_values(a_v); auto logdet = stan::math::sum(stan::math::log(w)); stan::math::set_zero_all_adjoints(); From 2ba5173623822fd5cb40b12497d9b7fd7345ed2d Mon Sep 17 00:00:00 2001 From: Dashadower Date: Thu, 10 Dec 2020 19:42:15 +0900 Subject: [PATCH 07/32] Remove unnecessary code from singular_values.hpp --- stan/math/rev/fun/singular_values.hpp | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index ee620bc62e7..330b889beda 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -45,29 +45,10 @@ inline auto singular_values(const EigMat& m) { //auto singular_values = to_arena(svd.singularValues()); auto U = to_arena(svd.matrixU()); auto V = to_arena(svd.matrixV()); - Eigen::MatrixXd Fp(M,M); - Eigen::MatrixXd Fm(M,M); - for(int i = 0; i < M; i++) { - for(int j = 0; j < M; j++) { - if(j == i) { - Fp(i, j) = 0.0; - Fm(i, j) = 0.0; - } else { - Fp(i, j) = 1.0 / (singular_values[j].val() - singular_values[i].val()) + 1.0 / (singular_values[i].val() + singular_values[j].val()); - Fm(i, j) = 1.0 / (singular_values[j].val() - singular_values[i].val()) - 1.0 / (singular_values[j].val() + singular_values[i].val()); - } - } - } - auto arena_Fp = to_arena(Fp); - auto arena_Fm = to_arena(Fm); reverse_pass_callback( [arena_m, U, singular_values, V]() mutable { - arena_m.adj() //+= 0.5 * U * (arena_Fp * (U.transpose() * U.adj() - U.adj().transpose() * U)) * V.transpose() // adjU contributions - //+(Eigen::MatrixXd::Identity(M, M) - U * U.transpose()) * U.adj() * singular_values.asDiagonal().inverse() * V.transpose(); - += U * singular_values.adj().asDiagonal() * V.transpose(); // adjD contributions - //+ 0.5 * U * (arena_Fm * (V.transpose() * V.adj() - V.adj().transpose() * V)) * V.transpose() // adjV contributions - //+ U * singular_values.asDiagonal().inverse() * V.adj().transpose() * (Eigen::MatrixXd::Identity(M, M) - V * V.transpose()); + arena_m.adj() += U * singular_values.adj().asDiagonal() * V.transpose(); // adjD contributions U.adj() and V.adj() is zero. }); return ret_type(singular_values); } From 81b442aad41107796b24d939a4d8dff8520f81bf Mon Sep 17 00:00:00 2001 From: Dashadower Date: Thu, 10 Dec 2020 21:25:35 +0900 Subject: [PATCH 08/32] Use arena_t and remove redundant comments in singular_values.hpp --- stan/math/rev/fun/singular_values.hpp | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index 330b889beda..118df09560d 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -20,7 +20,7 @@ namespace math { /** * Return the singular values of the specified matrix. * - * Equation (4) from Differentiable Programming Tensor Networks + * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650) * is used for MxN input matrix's adjoint calculation. * * @tparam EigMat type of input matrix @@ -40,15 +40,14 @@ inline auto singular_values(const EigMat& m) { Eigen::MatrixXd m_val = value_of(arena_m); Eigen::JacobiSVD svd(m_val, Eigen::ComputeFullU | Eigen::ComputeFullV); - vector_v singular_values(to_arena(svd.singularValues())); // switch to vector_v to enable adjoint + arena_t singular_values = svd.singularValues(); - //auto singular_values = to_arena(svd.singularValues()); auto U = to_arena(svd.matrixU()); auto V = to_arena(svd.matrixV()); reverse_pass_callback( [arena_m, U, singular_values, V]() mutable { - arena_m.adj() += U * singular_values.adj().asDiagonal() * V.transpose(); // adjD contributions U.adj() and V.adj() is zero. + arena_m.adj() += U * singular_values.adj().asDiagonal() * V.transpose(); }); return ret_type(singular_values); } From 8a3a324f560cb08f7f848193dccf32e857d714c7 Mon Sep 17 00:00:00 2001 From: Ben Date: Wed, 16 Dec 2020 16:53:48 -0500 Subject: [PATCH 09/32] Couple cleanup things (Issue #2192) --- stan/math/rev/fun/singular_values.hpp | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index 118df09560d..80af912ac58 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -1,17 +1,13 @@ -#ifndef MATH_SINGULAR_VALUES_H -#define MATH_SINGULAR_VALUES_H +#ifndef STAN_MATH_REV_FUN_SINGULAR_VALUES_HPP +#define STAN_MATH_REV_FUN_SINGULAR_VALUES_HPP #include #include #include -#include #include #include #include -#include #include -#include -#include #include namespace stan { @@ -28,7 +24,7 @@ namespace math { * @return Singular values of matrix */ template * = nullptr, + require_eigen_matrix_dynamic_t* = nullptr, require_vt_var* = nullptr> inline auto singular_values(const EigMat& m) { using ret_type = promote_scalar_t; @@ -37,8 +33,8 @@ inline auto singular_values(const EigMat& m) { const int M = to_arena(m.rows()); auto arena_m = to_arena(m); - Eigen::MatrixXd m_val = value_of(arena_m); - Eigen::JacobiSVD svd(m_val, Eigen::ComputeFullU | Eigen::ComputeFullV); + Eigen::JacobiSVD svd(arena_m.val(), + Eigen::ComputeFullU | Eigen::ComputeFullV); arena_t singular_values = svd.singularValues(); @@ -47,10 +43,14 @@ inline auto singular_values(const EigMat& m) { reverse_pass_callback( [arena_m, U, singular_values, V]() mutable { - arena_m.adj() += U * singular_values.adj().asDiagonal() * V.transpose(); + arena_m.adj() += U.block(0, 0, U.rows(), singular_values.size()) * + singular_values.adj().asDiagonal() * + V.block(0, 0, V.rows(), singular_values.size()).transpose(); }); + return ret_type(singular_values); } + } // namespace math } // namespace stan From f674c6cf98e78a5f644a2ddc960c80079ea4f3b4 Mon Sep 17 00:00:00 2001 From: Dashadower Date: Sat, 19 Dec 2020 01:19:28 +0900 Subject: [PATCH 10/32] WIP svd_U function and tests --- stan/math/rev/fun/singular_values_U.hpp | 73 +++++++++++++++++++ .../math/rev/fun/singular_values_U_test.cpp | 47 ++++++++++++ 2 files changed, 120 insertions(+) create mode 100644 stan/math/rev/fun/singular_values_U.hpp create mode 100644 test/unit/math/rev/fun/singular_values_U_test.cpp diff --git a/stan/math/rev/fun/singular_values_U.hpp b/stan/math/rev/fun/singular_values_U.hpp new file mode 100644 index 00000000000..12863ac0a67 --- /dev/null +++ b/stan/math/rev/fun/singular_values_U.hpp @@ -0,0 +1,73 @@ +#ifndef MATH_SINGULAR_VALUES_H +#define MATH_SINGULAR_VALUES_H + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Given input matrix m, return matrix U where \mathbf{m} = \mathbf{UDV^{T}} + * + * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650) + * is used for MxN input matrix's adjoint calculation. + * + * @tparam EigMat type of input matrix + * @param m MxN input matrix + * @return Orthogonal matrix U + */ +template * = nullptr, + require_vt_var* = nullptr> +inline auto singular_values_U(const EigMat& m) { + using ret_type = promote_scalar_t; + check_nonzero_size("singular_values", "m", m); + + const int M = to_arena(m.rows()); + auto arena_m = to_arena(m); + + Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeFullU | Eigen::ComputeFullV); + + auto arena_D = to_arena((Eigen::MatrixXd) svd.singularValues().asDiagonal()); + + Eigen::MatrixXd Fp(M, M); + + for(int i = 0; i < M; i++) { + for(int j = 0; j < M; j++) { + if(j == i) { + Fp(i, j) = 0.0; + } else { + Fp(i, j) = 1.0 / (arena_D(j, j) - arena_D(i, i)) + 1.0 / (arena_D(i, i) + arena_D(j, j)); + } + } + } + + auto arena_Fp = to_arena(Fp); + auto arena_U = to_arena(svd.matrixU()); + auto arena_V = to_arena(svd.matrixV()); + + + reverse_pass_callback( + [arena_m, arena_U, arena_D, arena_V, arena_Fp, M]() mutable { + arena_m.adj() += 0.5 * arena_U.val() * (arena_Fp.array() * (arena_U.val().transpose() * arena_U.adj() - arena_U.adj().transpose() * arena_U.val()).array()).matrix() + * arena_V.transpose() + (Eigen::MatrixXd::Identity(M, M) - arena_U.val() * arena_U.val().transpose()) + * arena_U.adj() * arena_D.inverse() * arena_V.transpose(); + }); + + return ret_type(arena_U); +} +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/rev/fun/singular_values_U_test.cpp b/test/unit/math/rev/fun/singular_values_U_test.cpp new file mode 100644 index 00000000000..6aef7e34ba8 --- /dev/null +++ b/test/unit/math/rev/fun/singular_values_U_test.cpp @@ -0,0 +1,47 @@ +#include +#include +#include +#include +#include +#include +#include +TEST(AgradRev, singularvalues_U_gradient) { + using namespace stan::math; + Eigen::MatrixXd A(4, 4); + // Random symmetric matrix + A << 1.8904, 0.7204, -0.1599, 1.2028, 0.7204, 7.3394, 2.0895, -0.6151, + -0.1599, 2.0895, 2.4601, -1.7219, 1.2028, -0.6151, -1.7219, 4.5260; + + Eigen::JacobiSVD svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV); + + auto U = svd.matrixU(); + Eigen::MatrixXd D = svd.singularValues().asDiagonal(); + auto V = svd.matrixV(); + + Eigen:: MatrixXd Fp(4, 4); + for(int i = 0; i < 4; i++){ + for(int j = 0; j < 4; j++){ + if(i == j){ + Fp(i, j) = 0.0; + } + else{ + Fp(i, j) = 1.0 / (D(j, j) - D(i, i)) + 1.0 / (D(i, i) + D(j, j)); + } + } + } + + Eigen::MatrixXd adj_A_fd = Eigen::MatrixXd::Zero(4, 4); + + Eigen::MatrixXd adj_U_fd = Eigen::MatrixXd::Ones(4, 4); + + adj_A_fd += 0.5 * U * (Fp.cwiseProduct(U.transpose() * adj_U_fd - adj_U_fd.transpose() * U)).matrix() * V.transpose() + + (Eigen::MatrixXd::Identity(4, 4) - U * U.transpose()) * adj_U_fd * D.inverse() * V.transpose(); + + + matrix_v arena_a(A); + auto stan_u = singular_values_U(arena_a); + + EXPECT_MATRIX_FLOAT_EQ(U, stan_u.val()); + + EXPECT_MATRIX_FLOAT_EQ(adj_A_fd, arena_a.adj()) +} \ No newline at end of file From 09e7ecf897da4e38984bf81ffe0f08a2d992eed5 Mon Sep 17 00:00:00 2001 From: Dashadower Date: Sat, 19 Dec 2020 01:26:12 +0900 Subject: [PATCH 11/32] Change to use matrix_v for matrix U --- stan/math/rev/fun/singular_values_U.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/rev/fun/singular_values_U.hpp b/stan/math/rev/fun/singular_values_U.hpp index 12863ac0a67..d0d12ee9e06 100644 --- a/stan/math/rev/fun/singular_values_U.hpp +++ b/stan/math/rev/fun/singular_values_U.hpp @@ -54,7 +54,8 @@ inline auto singular_values_U(const EigMat& m) { } auto arena_Fp = to_arena(Fp); - auto arena_U = to_arena(svd.matrixU()); + //auto arena_U = to_arena(svd.matrixU()); + arena_t arena_U = (svd.matrixU()); auto arena_V = to_arena(svd.matrixV()); From c097d118373a35d67c502536ed8a36ce40050c94 Mon Sep 17 00:00:00 2001 From: Ben Date: Fri, 18 Dec 2020 16:02:09 -0500 Subject: [PATCH 12/32] Merging in develop, adding .val_op() and .adj_op() --- stan/math/rev/fun/singular_values_U.hpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/stan/math/rev/fun/singular_values_U.hpp b/stan/math/rev/fun/singular_values_U.hpp index d0d12ee9e06..b740af3a2c9 100644 --- a/stan/math/rev/fun/singular_values_U.hpp +++ b/stan/math/rev/fun/singular_values_U.hpp @@ -3,10 +3,7 @@ #include #include -#include -#include #include -#include #include #include #include @@ -54,16 +51,17 @@ inline auto singular_values_U(const EigMat& m) { } auto arena_Fp = to_arena(Fp); - //auto arena_U = to_arena(svd.matrixU()); - arena_t arena_U = (svd.matrixU()); + arena_t arena_U = svd.matrixU(); auto arena_V = to_arena(svd.matrixV()); reverse_pass_callback( [arena_m, arena_U, arena_D, arena_V, arena_Fp, M]() mutable { - arena_m.adj() += 0.5 * arena_U.val() * (arena_Fp.array() * (arena_U.val().transpose() * arena_U.adj() - arena_U.adj().transpose() * arena_U.val()).array()).matrix() - * arena_V.transpose() + (Eigen::MatrixXd::Identity(M, M) - arena_U.val() * arena_U.val().transpose()) - * arena_U.adj() * arena_D.inverse() * arena_V.transpose(); + Eigen::MatrixXd UUadjT = arena_U.val_op().transpose() * arena_U.adj_op(); + arena_m.adj() += 0.5 * arena_U.val_op() * + (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()).matrix() + * arena_V.transpose() + + (Eigen::MatrixXd::Identity(M, M) - arena_U.val_op() * arena_U.val_op().transpose()) * arena_U.adj_op() * arena_D.inverse() * arena_V.transpose(); }); return ret_type(arena_U); From 3b33ec8f4744c3e6242f24661b1a76a239807b69 Mon Sep 17 00:00:00 2001 From: Dashadower Date: Sun, 20 Dec 2020 23:17:40 +0900 Subject: [PATCH 13/32] Minor tweaks to svd_U and corresponding test suite --- stan/math/rev/fun/singular_values_U.hpp | 4 +++- test/unit/math/rev/fun/singular_values_U_test.cpp | 12 ++++++------ 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/stan/math/rev/fun/singular_values_U.hpp b/stan/math/rev/fun/singular_values_U.hpp index b740af3a2c9..39d34b0a7b4 100644 --- a/stan/math/rev/fun/singular_values_U.hpp +++ b/stan/math/rev/fun/singular_values_U.hpp @@ -17,6 +17,8 @@ namespace math { /** * Given input matrix m, return matrix U where \mathbf{m} = \mathbf{UDV^{T}} * + * This function internally calls Eigen::JacobiSVD to compute decomposition. + * * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650) * is used for MxN input matrix's adjoint calculation. * @@ -27,7 +29,7 @@ namespace math { template * = nullptr, require_vt_var* = nullptr> -inline auto singular_values_U(const EigMat& m) { +inline auto svd_U(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("singular_values", "m", m); diff --git a/test/unit/math/rev/fun/singular_values_U_test.cpp b/test/unit/math/rev/fun/singular_values_U_test.cpp index 6aef7e34ba8..7f879b367c6 100644 --- a/test/unit/math/rev/fun/singular_values_U_test.cpp +++ b/test/unit/math/rev/fun/singular_values_U_test.cpp @@ -5,7 +5,7 @@ #include #include #include -TEST(AgradRev, singularvalues_U_gradient) { +TEST(AgradRev, svd_U_gradient) { using namespace stan::math; Eigen::MatrixXd A(4, 4); // Random symmetric matrix @@ -34,14 +34,14 @@ TEST(AgradRev, singularvalues_U_gradient) { Eigen::MatrixXd adj_U_fd = Eigen::MatrixXd::Ones(4, 4); - adj_A_fd += 0.5 * U * (Fp.cwiseProduct(U.transpose() * adj_U_fd - adj_U_fd.transpose() * U)).matrix() * V.transpose() + adj_A_fd += 0.5 * U * + (Fp.array() * (U.transpose() * adj_A_fd - adj_A_fd.transpose() * U).array()).matrix() + (Eigen::MatrixXd::Identity(4, 4) - U * U.transpose()) * adj_U_fd * D.inverse() * V.transpose(); + stan::math::set_zero_all_adjoints(); matrix_v arena_a(A); - auto stan_u = singular_values_U(arena_a); - + auto stan_u = svd_U(arena_a); EXPECT_MATRIX_FLOAT_EQ(U, stan_u.val()); - - EXPECT_MATRIX_FLOAT_EQ(adj_A_fd, arena_a.adj()) + EXPECT_MATRIX_FLOAT_EQ(adj_A_fd, arena_a.adj_op()); } \ No newline at end of file From 156aceaaa1b96fb096c5fd7bb7ed5fe288832fcd Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 20 Dec 2020 09:52:45 -0500 Subject: [PATCH 14/32] Renamed `singular_values_U` to `svd_U`, add prim version of `svd_U`, and added mixed mode tests (Issue #2101) --- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/svd_U.hpp | 32 +++++++++++++++++++ stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/singular_values.hpp | 3 -- .../fun/{singular_values_U.hpp => svd_U.hpp} | 6 ++-- test/unit/math/mix/fun/svd_U_test.cpp | 30 +++++++++++++++++ 6 files changed, 67 insertions(+), 6 deletions(-) create mode 100644 stan/math/prim/fun/svd_U.hpp rename stan/math/rev/fun/{singular_values_U.hpp => svd_U.hpp} (94%) create mode 100644 test/unit/math/mix/fun/svd_U_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index f4e51d5d75b..866e6bbba5d 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -305,6 +305,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/svd_U.hpp b/stan/math/prim/fun/svd_U.hpp new file mode 100644 index 00000000000..57054b7596b --- /dev/null +++ b/stan/math/prim/fun/svd_U.hpp @@ -0,0 +1,32 @@ +#ifndef STAN_MATH_PRIM_FUN_SVD_U_HPP +#define STAN_MATH_PRIM_FUN_SVD_U_HPP + +#include +#include + +namespace stan { +namespace math { + +/** + * Given input matrix m, return matrix U where \mathbf{m} = \mathbf{UDV^{T}} + * + * @tparam EigMat type of the matrix + * @param m MxN input matrix + * @return Orthogonal matrix U + */ +template * = nullptr, + require_not_vt_var* = nullptr> +Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> +svd_U(const EigMat& m) { + check_nonzero_size("svd_U", "m", m); + + return Eigen::JacobiSVD, Eigen::Dynamic, + Eigen::Dynamic> >(m, Eigen::ComputeFullU) + .matrixU(); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 9d37cfdee0d..601fb34ffe1 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -133,6 +133,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index 80af912ac58..eb3970e73d5 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -3,9 +3,6 @@ #include #include -#include -#include -#include #include #include #include diff --git a/stan/math/rev/fun/singular_values_U.hpp b/stan/math/rev/fun/svd_U.hpp similarity index 94% rename from stan/math/rev/fun/singular_values_U.hpp rename to stan/math/rev/fun/svd_U.hpp index 39d34b0a7b4..b2780ccc323 100644 --- a/stan/math/rev/fun/singular_values_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -1,5 +1,5 @@ -#ifndef MATH_SINGULAR_VALUES_H -#define MATH_SINGULAR_VALUES_H +#ifndef STAN_MATH_REV_FUN_SVD_U_HPP +#define STAN_MATH_REV_FUN_SVD_U_HPP #include #include @@ -27,7 +27,7 @@ namespace math { * @return Orthogonal matrix U */ template * = nullptr, + require_eigen_matrix_dynamic_t* = nullptr, require_vt_var* = nullptr> inline auto svd_U(const EigMat& m) { using ret_type = promote_scalar_t; diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp new file mode 100644 index 00000000000..dd5c2779e64 --- /dev/null +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -0,0 +1,30 @@ +#include + +TEST(MathMixMatFun, svd_U) { + auto f = [](const auto& x) { return stan::math::svd_U(x); }; + + stan::test::ad_tolerances tols; + //tols.hessian_hessian_ = 1e-2; + //tols.hessian_fvar_hessian_ = 1e-2; + + Eigen::MatrixXd m00(0, 0); + stan::test::expect_ad(f, m00); + + Eigen::MatrixXd m11(1, 1); + m11 << 1.1; + stan::test::expect_ad(tols, f, m11); + + /*Eigen::MatrixXd m22(2, 2); + m22 << 3, -5, 7, 11; + stan::test::expect_ad(tols, f, m22); + + Eigen::MatrixXd m23(2, 3); + m23 << 3, 5, -7, -11, 13, -17; + Eigen::MatrixXd m32 = m23.transpose(); + stan::test::expect_ad(tols, f, m23); + stan::test::expect_ad(tols, f, m32); + + Eigen::MatrixXd a22(2, 2); + a22 << 1, 2, 3, 4; + stan::test::expect_ad(tols, f, a22);*/ +} From 7ba6df0a582aa3204ee23f92c00e7781ca96c23e Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 20 Dec 2020 09:56:20 -0500 Subject: [PATCH 15/32] Updated tests again (Issue #2101) --- test/unit/math/mix/fun/svd_U_test.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index dd5c2779e64..53135ad0a53 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -14,17 +14,17 @@ TEST(MathMixMatFun, svd_U) { m11 << 1.1; stan::test::expect_ad(tols, f, m11); - /*Eigen::MatrixXd m22(2, 2); + Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; stan::test::expect_ad(tols, f, m22); - Eigen::MatrixXd m23(2, 3); + /*Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; Eigen::MatrixXd m32 = m23.transpose(); stan::test::expect_ad(tols, f, m23); - stan::test::expect_ad(tols, f, m32); + stan::test::expect_ad(tols, f, m32);*/ Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; - stan::test::expect_ad(tols, f, a22);*/ + stan::test::expect_ad(tols, f, a22); } From aa43a00bd2a93bb0f067602458e16daf9159fa5b Mon Sep 17 00:00:00 2001 From: Dashadower Date: Mon, 21 Dec 2020 01:38:03 +0900 Subject: [PATCH 16/32] Add svd_V function --- stan/math/prim/fun.hpp | 1 + stan/math/prim/fun/svd_V.hpp | 32 ++++++++ stan/math/rev/fun.hpp | 1 + stan/math/rev/fun/svd_U.hpp | 2 +- stan/math/rev/fun/svd_V.hpp | 74 +++++++++++++++++++ test/unit/math/mix/fun/svd_V_test.cpp | 30 ++++++++ .../math/rev/fun/singular_values_U_test.cpp | 47 ------------ 7 files changed, 139 insertions(+), 48 deletions(-) create mode 100644 stan/math/prim/fun/svd_V.hpp create mode 100644 stan/math/rev/fun/svd_V.hpp create mode 100644 test/unit/math/mix/fun/svd_V_test.cpp delete mode 100644 test/unit/math/rev/fun/singular_values_U_test.cpp diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp index 866e6bbba5d..004e2934737 100644 --- a/stan/math/prim/fun.hpp +++ b/stan/math/prim/fun.hpp @@ -306,6 +306,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/prim/fun/svd_V.hpp b/stan/math/prim/fun/svd_V.hpp new file mode 100644 index 00000000000..d261eb16e26 --- /dev/null +++ b/stan/math/prim/fun/svd_V.hpp @@ -0,0 +1,32 @@ +#ifndef STAN_MATH_PRIM_FUN_SVD_V_HPP +#define STAN_MATH_PRIM_FUN_SVD_V_HPP + +#include +#include + +namespace stan { +namespace math { + +/** + * Given input matrix m, return matrix V where \mathbf{m} = \mathbf{UDV^{T}} + * + * @tparam EigMat type of the matrix + * @param m MxN input matrix + * @return Orthogonal matrix V + */ +template * = nullptr, + require_not_vt_var* = nullptr> +Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> +svd_V(const EigMat& m) { + check_nonzero_size("svd_V", "m", m); + + return Eigen::JacobiSVD, Eigen::Dynamic, + Eigen::Dynamic> >(m, Eigen::ComputeFullV) + .matrixV(); +} + +} // namespace math +} // namespace stan + +#endif diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp index 601fb34ffe1..0aa512881e8 100644 --- a/stan/math/rev/fun.hpp +++ b/stan/math/rev/fun.hpp @@ -134,6 +134,7 @@ #include #include #include +#include #include #include #include diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index b2780ccc323..df4cd6e41bb 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -31,7 +31,7 @@ template * = nullptr> inline auto svd_U(const EigMat& m) { using ret_type = promote_scalar_t; - check_nonzero_size("singular_values", "m", m); + check_nonzero_size("svd_U", "m", m); const int M = to_arena(m.rows()); auto arena_m = to_arena(m); diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp new file mode 100644 index 00000000000..e129c050308 --- /dev/null +++ b/stan/math/rev/fun/svd_V.hpp @@ -0,0 +1,74 @@ +#ifndef STAN_MATH_REV_FUN_SVD_V_HPP +#define STAN_MATH_REV_FUN_SVD_V_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace stan { +namespace math { + +/** + * Given input matrix m, return matrix V where \mathbf{m} = \mathbf{UDV^{T}} + * + * This function internally calls Eigen::JacobiSVD to compute decomposition. + * + * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650) + * is used for MxN input matrix's adjoint calculation. + * + * @tparam EigMat type of input matrix + * @param m MxN input matrix + * @return Orthogonal matrix V + */ +template * = nullptr, + require_vt_var* = nullptr> +inline auto svd_V(const EigMat& m) { + using ret_type = promote_scalar_t; + check_nonzero_size("svd_V", "m", m); + + const int M = to_arena(m.rows()); + auto arena_m = to_arena(m); + + Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeFullU | Eigen::ComputeFullV); + + auto arena_D = to_arena((Eigen::MatrixXd) svd.singularValues().asDiagonal()); + + Eigen::MatrixXd Fm(M, M); + + for(int i = 0; i < M; i++) { + for(int j = 0; j < M; j++) { + if(j == i) { + Fm(i, j) = 0.0; + } else { + Fm(i, j) = 1.0 / (arena_D(j, j) - arena_D(i, i)) - 1.0 / (arena_D(i, i) + arena_D(j, j)); + } + } + } + + auto arena_Fm = to_arena(Fm); + auto arena_U = to_arena(svd.matrixU()); + arena_t arena_V = svd.matrixV(); + + + reverse_pass_callback( + [arena_m, arena_U, arena_D, arena_V, arena_Fm, M]() mutable { + Eigen::MatrixXd VTVadj = arena_V.val_op().transpose() * arena_V.adj_op(); + arena_m.adj() += 0.5 * arena_U.val_op() * + (arena_Fm.array() * (VTVadj - VTVadj.transpose()).array()).matrix() * + arena_V.val_op().transpose() + arena_U.val_op() * arena_D.val_op().inverse() * arena_V.adj_op().transpose() * + (Eigen::MatrixXd::Identity(M, M) - arena_V.val_op() * arena_V.val_op().transpose()); + }); + + return ret_type(arena_V); +} +} // namespace math +} // namespace stan + +#endif diff --git a/test/unit/math/mix/fun/svd_V_test.cpp b/test/unit/math/mix/fun/svd_V_test.cpp new file mode 100644 index 00000000000..90b00d3092e --- /dev/null +++ b/test/unit/math/mix/fun/svd_V_test.cpp @@ -0,0 +1,30 @@ +#include + +TEST(MathMixMatFun, svd_V) { + auto f = [](const auto& x) { return stan::math::svd_V(x); }; + + stan::test::ad_tolerances tols; + //tols.hessian_hessian_ = 1e-2; + //tols.hessian_fvar_hessian_ = 1e-2; + + Eigen::MatrixXd m00(0, 0); + stan::test::expect_ad(f, m00); + + Eigen::MatrixXd m11(1, 1); + m11 << 1.1; + stan::test::expect_ad(tols, f, m11); + + Eigen::MatrixXd m22(2, 2); + m22 << 3, -5, 7, 11; + stan::test::expect_ad(tols, f, m22); + + /*Eigen::MatrixXd m23(2, 3); + m23 << 3, 5, -7, -11, 13, -17; + Eigen::MatrixXd m32 = m23.transpose(); + stan::test::expect_ad(tols, f, m23); + stan::test::expect_ad(tols, f, m32);*/ + + Eigen::MatrixXd a22(2, 2); + a22 << 1, 2, 3, 4; + stan::test::expect_ad(tols, f, a22); +} diff --git a/test/unit/math/rev/fun/singular_values_U_test.cpp b/test/unit/math/rev/fun/singular_values_U_test.cpp deleted file mode 100644 index 7f879b367c6..00000000000 --- a/test/unit/math/rev/fun/singular_values_U_test.cpp +++ /dev/null @@ -1,47 +0,0 @@ -#include -#include -#include -#include -#include -#include -#include -TEST(AgradRev, svd_U_gradient) { - using namespace stan::math; - Eigen::MatrixXd A(4, 4); - // Random symmetric matrix - A << 1.8904, 0.7204, -0.1599, 1.2028, 0.7204, 7.3394, 2.0895, -0.6151, - -0.1599, 2.0895, 2.4601, -1.7219, 1.2028, -0.6151, -1.7219, 4.5260; - - Eigen::JacobiSVD svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV); - - auto U = svd.matrixU(); - Eigen::MatrixXd D = svd.singularValues().asDiagonal(); - auto V = svd.matrixV(); - - Eigen:: MatrixXd Fp(4, 4); - for(int i = 0; i < 4; i++){ - for(int j = 0; j < 4; j++){ - if(i == j){ - Fp(i, j) = 0.0; - } - else{ - Fp(i, j) = 1.0 / (D(j, j) - D(i, i)) + 1.0 / (D(i, i) + D(j, j)); - } - } - } - - Eigen::MatrixXd adj_A_fd = Eigen::MatrixXd::Zero(4, 4); - - Eigen::MatrixXd adj_U_fd = Eigen::MatrixXd::Ones(4, 4); - - adj_A_fd += 0.5 * U * - (Fp.array() * (U.transpose() * adj_A_fd - adj_A_fd.transpose() * U).array()).matrix() - + (Eigen::MatrixXd::Identity(4, 4) - U * U.transpose()) * adj_U_fd * D.inverse() * V.transpose(); - - - stan::math::set_zero_all_adjoints(); - matrix_v arena_a(A); - auto stan_u = svd_U(arena_a); - EXPECT_MATRIX_FLOAT_EQ(U, stan_u.val()); - EXPECT_MATRIX_FLOAT_EQ(adj_A_fd, arena_a.adj_op()); -} \ No newline at end of file From 7df73072571cdea630ab8b06a0a8d1adeea07511 Mon Sep 17 00:00:00 2001 From: Dashadower Date: Mon, 21 Dec 2020 03:33:59 +0900 Subject: [PATCH 17/32] Address non-square inputs for svd_U and svd_V --- stan/math/rev/fun/svd_U.hpp | 25 ++++++++++++++----------- stan/math/rev/fun/svd_V.hpp | 24 +++++++++++++----------- test/unit/math/mix/fun/svd_U_test.cpp | 8 ++++---- test/unit/math/mix/fun/svd_V_test.cpp | 8 ++++---- 4 files changed, 35 insertions(+), 30 deletions(-) diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index df4cd6e41bb..0aa071cdeb4 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -33,12 +33,12 @@ inline auto svd_U(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("svd_U", "m", m); - const int M = to_arena(m.rows()); - auto arena_m = to_arena(m); + const int M = to_arena(std::min(m.rows(), m.cols())); + auto arena_m = to_arena(m); // N by P Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeFullU | Eigen::ComputeFullV); - auto arena_D = to_arena((Eigen::MatrixXd) svd.singularValues().asDiagonal()); + auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M Eigen::MatrixXd Fp(M, M); @@ -47,23 +47,26 @@ inline auto svd_U(const EigMat& m) { if(j == i) { Fp(i, j) = 0.0; } else { - Fp(i, j) = 1.0 / (arena_D(j, j) - arena_D(i, i)) + 1.0 / (arena_D(i, i) + arena_D(j, j)); + Fp(i, j) = 1.0 / (arena_D[j] - arena_D[i]) + 1.0 / (arena_D[i] + arena_D[j]); } } } auto arena_Fp = to_arena(Fp); - arena_t arena_U = svd.matrixU(); - auto arena_V = to_arena(svd.matrixV()); + arena_t arena_U = svd.matrixU(); // N by N + auto arena_V = to_arena(svd.matrixV()); // P by P reverse_pass_callback( [arena_m, arena_U, arena_D, arena_V, arena_Fp, M]() mutable { - Eigen::MatrixXd UUadjT = arena_U.val_op().transpose() * arena_U.adj_op(); - arena_m.adj() += 0.5 * arena_U.val_op() * - (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()).matrix() - * arena_V.transpose() + - (Eigen::MatrixXd::Identity(M, M) - arena_U.val_op() * arena_U.val_op().transpose()) * arena_U.adj_op() * arena_D.inverse() * arena_V.transpose(); + auto U_block = arena_U.block(0, 0, arena_m.rows(), M); // N by M + auto V_block = arena_V.block(0, 0, arena_m.cols(), M); // P by M + Eigen::MatrixXd UUadjT = U_block.val_op().transpose() * U_block.adj_op(); // M by M + arena_m.adj() += 0.5 * U_block.val_op() * + (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()).matrix() + * V_block.transpose() + + (Eigen::MatrixXd::Identity(U_block.rows(), U_block.rows()) - U_block.val_op() * U_block.val_op().transpose()) * + U_block.adj_op() * arena_D.asDiagonal().inverse() * V_block.transpose(); }); return ret_type(arena_U); diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index e129c050308..30ec1b22198 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -33,12 +33,12 @@ inline auto svd_V(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("svd_V", "m", m); - const int M = to_arena(m.rows()); - auto arena_m = to_arena(m); + const int M = to_arena(std::min(m.rows(), m.cols())); + auto arena_m = to_arena(m); // N by P Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeFullU | Eigen::ComputeFullV); - auto arena_D = to_arena((Eigen::MatrixXd) svd.singularValues().asDiagonal()); + auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M Eigen::MatrixXd Fm(M, M); @@ -47,23 +47,25 @@ inline auto svd_V(const EigMat& m) { if(j == i) { Fm(i, j) = 0.0; } else { - Fm(i, j) = 1.0 / (arena_D(j, j) - arena_D(i, i)) - 1.0 / (arena_D(i, i) + arena_D(j, j)); + Fm(i, j) = 1.0 / (arena_D[j] - arena_D[i]) - 1.0 / (arena_D[i] + arena_D[j]); } } } auto arena_Fm = to_arena(Fm); - auto arena_U = to_arena(svd.matrixU()); - arena_t arena_V = svd.matrixV(); + auto arena_U = to_arena(svd.matrixU()); // N by N + arena_t arena_V = svd.matrixV(); // P by P reverse_pass_callback( [arena_m, arena_U, arena_D, arena_V, arena_Fm, M]() mutable { - Eigen::MatrixXd VTVadj = arena_V.val_op().transpose() * arena_V.adj_op(); - arena_m.adj() += 0.5 * arena_U.val_op() * - (arena_Fm.array() * (VTVadj - VTVadj.transpose()).array()).matrix() * - arena_V.val_op().transpose() + arena_U.val_op() * arena_D.val_op().inverse() * arena_V.adj_op().transpose() * - (Eigen::MatrixXd::Identity(M, M) - arena_V.val_op() * arena_V.val_op().transpose()); + auto U_block = arena_U.block(0, 0, arena_m.rows(), M); // N by M + auto V_block = arena_V.block(0, 0, arena_m.cols(), M); // P by M + Eigen::MatrixXd VTVadj = V_block.val_op().transpose() * V_block.adj_op(); // M by M + arena_m.adj() += 0.5 * U_block * + (arena_Fm.array() * (VTVadj - VTVadj.transpose()).array()).matrix() * + V_block.val_op().transpose() + U_block * arena_D.asDiagonal().inverse() * V_block.adj_op().transpose() * + (Eigen::MatrixXd::Identity(arena_m.cols(), arena_m.cols()) - V_block.val_op() * V_block.val_op().transpose()); }); return ret_type(arena_V); diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index 53135ad0a53..ca496797d99 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -4,8 +4,8 @@ TEST(MathMixMatFun, svd_U) { auto f = [](const auto& x) { return stan::math::svd_U(x); }; stan::test::ad_tolerances tols; - //tols.hessian_hessian_ = 1e-2; - //tols.hessian_fvar_hessian_ = 1e-2; + tols.hessian_hessian_ = 1e-2; + tols.hessian_fvar_hessian_ = 1e-2; Eigen::MatrixXd m00(0, 0); stan::test::expect_ad(f, m00); @@ -18,11 +18,11 @@ TEST(MathMixMatFun, svd_U) { m22 << 3, -5, 7, 11; stan::test::expect_ad(tols, f, m22); - /*Eigen::MatrixXd m23(2, 3); + Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; Eigen::MatrixXd m32 = m23.transpose(); stan::test::expect_ad(tols, f, m23); - stan::test::expect_ad(tols, f, m32);*/ + stan::test::expect_ad(tols, f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; diff --git a/test/unit/math/mix/fun/svd_V_test.cpp b/test/unit/math/mix/fun/svd_V_test.cpp index 90b00d3092e..d1bbdc80f3b 100644 --- a/test/unit/math/mix/fun/svd_V_test.cpp +++ b/test/unit/math/mix/fun/svd_V_test.cpp @@ -4,8 +4,8 @@ TEST(MathMixMatFun, svd_V) { auto f = [](const auto& x) { return stan::math::svd_V(x); }; stan::test::ad_tolerances tols; - //tols.hessian_hessian_ = 1e-2; - //tols.hessian_fvar_hessian_ = 1e-2; + tols.hessian_hessian_ = 1e-2; + tols.hessian_fvar_hessian_ = 1e-2; Eigen::MatrixXd m00(0, 0); stan::test::expect_ad(f, m00); @@ -18,11 +18,11 @@ TEST(MathMixMatFun, svd_V) { m22 << 3, -5, 7, 11; stan::test::expect_ad(tols, f, m22); - /*Eigen::MatrixXd m23(2, 3); + Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; Eigen::MatrixXd m32 = m23.transpose(); stan::test::expect_ad(tols, f, m23); - stan::test::expect_ad(tols, f, m32);*/ + stan::test::expect_ad(tols, f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; From e1c7dcca63ddb215557db6df90484677f48e37ae Mon Sep 17 00:00:00 2001 From: Dashadower Date: Mon, 21 Dec 2020 03:54:24 +0900 Subject: [PATCH 18/32] Add prim tests for svd_V, svd_U --- test/unit/math/prim/fun/svd_U_test.cpp | 14 ++++++++++++++ test/unit/math/prim/fun/svd_V_test.cpp | 15 +++++++++++++++ 2 files changed, 29 insertions(+) create mode 100644 test/unit/math/prim/fun/svd_U_test.cpp create mode 100644 test/unit/math/prim/fun/svd_V_test.cpp diff --git a/test/unit/math/prim/fun/svd_U_test.cpp b/test/unit/math/prim/fun/svd_U_test.cpp new file mode 100644 index 00000000000..bc19df5cb20 --- /dev/null +++ b/test/unit/math/prim/fun/svd_U_test.cpp @@ -0,0 +1,14 @@ +#include +#include +#include + +TEST(MathMatrixPrimMat, svd_U) { + using stan::math::svd_U; + using stan::math::matrix_d; + matrix_d m23(2, 3); + m23 << 1, 3, -5, 7, 9, -11; + + matrix_d m23_U(2, 2); + m23_U << -0.3378432, -0.9412024, -0.9412024, 0.3378432; //Generated using R base::svd + EXPECT_MATRIX_FLOAT_EQ(m23_U, svd_U(m23)); +} \ No newline at end of file diff --git a/test/unit/math/prim/fun/svd_V_test.cpp b/test/unit/math/prim/fun/svd_V_test.cpp new file mode 100644 index 00000000000..733c82c2a47 --- /dev/null +++ b/test/unit/math/prim/fun/svd_V_test.cpp @@ -0,0 +1,15 @@ +#include +#include +#include + +TEST(MathMatrixPrimMat, svd_V) { + using stan::math::svd_V; + using stan::math::matrix_d; + matrix_d m23(2, 3); + m23 << 1, 3, -5, 7, 9, -11; + + matrix_d m23_V(3, 3); + m23_V << -0.4117624, 0.8147301, -0.4082483, -0.5638395, + 0.1241705, 0.8164966, 0.7159167, 0.5663891, 0.4082483; //Generated using R base::svd + EXPECT_MATRIX_FLOAT_EQ(m23_V, svd_V(m23)); +} \ No newline at end of file From 7ec29e83b70cf6f11b093edc7c202253f6c5e8cb Mon Sep 17 00:00:00 2001 From: Ben Date: Mon, 21 Dec 2020 22:04:08 -0500 Subject: [PATCH 19/32] Updated svd V test --- test/unit/math/prim/fun/svd_V_test.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/unit/math/prim/fun/svd_V_test.cpp b/test/unit/math/prim/fun/svd_V_test.cpp index 733c82c2a47..ed6010b71b6 100644 --- a/test/unit/math/prim/fun/svd_V_test.cpp +++ b/test/unit/math/prim/fun/svd_V_test.cpp @@ -9,7 +9,8 @@ TEST(MathMatrixPrimMat, svd_V) { m23 << 1, 3, -5, 7, 9, -11; matrix_d m23_V(3, 3); - m23_V << -0.4117624, 0.8147301, -0.4082483, -0.5638395, - 0.1241705, 0.8164966, 0.7159167, 0.5663891, 0.4082483; //Generated using R base::svd + m23_V << -0.41176240532160857, 0.81473005032163681, -0.40824829046386291, + -0.56383954240865775, 0.12417046246885260, 0.81649658092772603, + 0.71591667949570703, 0.56638912538393127, 0.40824829046386307; //Generated using R base::svd EXPECT_MATRIX_FLOAT_EQ(m23_V, svd_V(m23)); -} \ No newline at end of file +} From a7c21ec61c49ba06c482ea70bdd7489304dedb1d Mon Sep 17 00:00:00 2001 From: Dashadower Date: Sun, 27 Dec 2020 12:14:01 +0900 Subject: [PATCH 20/32] Update tests for prim/fun/{svd_U, svd_V}.hpp. Change rev/fun/{svd_U, svd_V} functions to return thinned matrices instead of full --- stan/math/prim/fun/svd_U.hpp | 2 +- stan/math/prim/fun/svd_V.hpp | 2 +- stan/math/rev/fun/svd_U.hpp | 20 +++++++--------- stan/math/rev/fun/svd_V.hpp | 16 ++++++------- test/unit/math/mix/fun/svd_U_test.cpp | 10 +++++--- test/unit/math/prim/fun/svd_U_test.cpp | 29 ++++++++++++++++++++-- test/unit/math/prim/fun/svd_V_test.cpp | 33 ++++++++++++++++++++++---- 7 files changed, 80 insertions(+), 32 deletions(-) diff --git a/stan/math/prim/fun/svd_U.hpp b/stan/math/prim/fun/svd_U.hpp index 57054b7596b..dee604102a8 100644 --- a/stan/math/prim/fun/svd_U.hpp +++ b/stan/math/prim/fun/svd_U.hpp @@ -22,7 +22,7 @@ svd_U(const EigMat& m) { check_nonzero_size("svd_U", "m", m); return Eigen::JacobiSVD, Eigen::Dynamic, - Eigen::Dynamic> >(m, Eigen::ComputeFullU) + Eigen::Dynamic> >(m, Eigen::ComputeThinU) .matrixU(); } diff --git a/stan/math/prim/fun/svd_V.hpp b/stan/math/prim/fun/svd_V.hpp index d261eb16e26..ac07059fce8 100644 --- a/stan/math/prim/fun/svd_V.hpp +++ b/stan/math/prim/fun/svd_V.hpp @@ -22,7 +22,7 @@ svd_V(const EigMat& m) { check_nonzero_size("svd_V", "m", m); return Eigen::JacobiSVD, Eigen::Dynamic, - Eigen::Dynamic> >(m, Eigen::ComputeFullV) + Eigen::Dynamic> >(m, Eigen::ComputeThinV) .matrixV(); } diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index 0aa071cdeb4..0cf4bbea75b 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -36,7 +36,7 @@ inline auto svd_U(const EigMat& m) { const int M = to_arena(std::min(m.rows(), m.cols())); auto arena_m = to_arena(m); // N by P - Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeFullU | Eigen::ComputeFullV); + Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M @@ -53,20 +53,18 @@ inline auto svd_U(const EigMat& m) { } auto arena_Fp = to_arena(Fp); - arena_t arena_U = svd.matrixU(); // N by N - auto arena_V = to_arena(svd.matrixV()); // P by P + arena_t arena_U = svd.matrixU(); // N by M + auto arena_V = to_arena(svd.matrixV()); // P by M reverse_pass_callback( [arena_m, arena_U, arena_D, arena_V, arena_Fp, M]() mutable { - auto U_block = arena_U.block(0, 0, arena_m.rows(), M); // N by M - auto V_block = arena_V.block(0, 0, arena_m.cols(), M); // P by M - Eigen::MatrixXd UUadjT = U_block.val_op().transpose() * U_block.adj_op(); // M by M - arena_m.adj() += 0.5 * U_block.val_op() * - (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()).matrix() - * V_block.transpose() + - (Eigen::MatrixXd::Identity(U_block.rows(), U_block.rows()) - U_block.val_op() * U_block.val_op().transpose()) * - U_block.adj_op() * arena_D.asDiagonal().inverse() * V_block.transpose(); + Eigen::MatrixXd UUadjT = arena_U.val_op().transpose() * arena_U.adj_op(); // M by M + arena_m.adj() += .5 * arena_U.val_op() * + (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()).matrix() + * arena_V.transpose() + + (Eigen::MatrixXd::Identity(arena_m.rows(), arena_m.rows()) - arena_U.val_op() * arena_U.val_op().transpose()) * + arena_U.adj_op() * arena_D.asDiagonal().inverse() * arena_V.transpose(); }); return ret_type(arena_U); diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index 30ec1b22198..a92776d795e 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -36,7 +36,7 @@ inline auto svd_V(const EigMat& m) { const int M = to_arena(std::min(m.rows(), m.cols())); auto arena_m = to_arena(m); // N by P - Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeFullU | Eigen::ComputeFullV); + Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M @@ -53,19 +53,17 @@ inline auto svd_V(const EigMat& m) { } auto arena_Fm = to_arena(Fm); - auto arena_U = to_arena(svd.matrixU()); // N by N - arena_t arena_V = svd.matrixV(); // P by P + auto arena_U = to_arena(svd.matrixU()); // N by M + arena_t arena_V = svd.matrixV(); // P by M reverse_pass_callback( [arena_m, arena_U, arena_D, arena_V, arena_Fm, M]() mutable { - auto U_block = arena_U.block(0, 0, arena_m.rows(), M); // N by M - auto V_block = arena_V.block(0, 0, arena_m.cols(), M); // P by M - Eigen::MatrixXd VTVadj = V_block.val_op().transpose() * V_block.adj_op(); // M by M - arena_m.adj() += 0.5 * U_block * + Eigen::MatrixXd VTVadj = arena_V.val_op().transpose() * arena_V.adj_op(); // M by M + arena_m.adj() += 0.5 * arena_U * (arena_Fm.array() * (VTVadj - VTVadj.transpose()).array()).matrix() * - V_block.val_op().transpose() + U_block * arena_D.asDiagonal().inverse() * V_block.adj_op().transpose() * - (Eigen::MatrixXd::Identity(arena_m.cols(), arena_m.cols()) - V_block.val_op() * V_block.val_op().transpose()); + arena_V.val_op().transpose() + arena_U * arena_D.asDiagonal().inverse() * arena_V.adj_op().transpose() * + (Eigen::MatrixXd::Identity(arena_m.cols(), arena_m.cols()) - arena_V.val_op() * arena_V.val_op().transpose()); }); return ret_type(arena_V); diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index ca496797d99..d8e847dd90e 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -1,5 +1,5 @@ #include - +//#include TEST(MathMixMatFun, svd_U) { auto f = [](const auto& x) { return stan::math::svd_U(x); }; @@ -17,12 +17,16 @@ TEST(MathMixMatFun, svd_U) { Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; stan::test::expect_ad(tols, f, m22); - + Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; Eigen::MatrixXd m32 = m23.transpose(); + /*stan::math::matrix_v test_arr = m32; + std::cout << test_arr.val_op() << std::endl; + std::cout << "----------" << std::endl; + std::cout << f(test_arr) << std:: endl;*/ stan::test::expect_ad(tols, f, m23); - stan::test::expect_ad(tols, f, m32); + //stan::test::expect_ad(tols, f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; diff --git a/test/unit/math/prim/fun/svd_U_test.cpp b/test/unit/math/prim/fun/svd_U_test.cpp index bc19df5cb20..6b3eff6e5a7 100644 --- a/test/unit/math/prim/fun/svd_U_test.cpp +++ b/test/unit/math/prim/fun/svd_U_test.cpp @@ -1,14 +1,39 @@ #include #include #include +#include TEST(MathMatrixPrimMat, svd_U) { using stan::math::svd_U; using stan::math::matrix_d; + + //Values generated using R base::svd + + matrix_d m00(0, 0); + EXPECT_THROW(svd_U(m00), std::invalid_argument); + + matrix_d m11(1, 1); + m11 << 5; + matrix_d m11_U(1,1); + m11_U << 1; + EXPECT_MATRIX_FLOAT_EQ(m11_U, svd_U(m11)); + + matrix_d m22(2, 2); + m22 << 1, 9, -4, 2; + matrix_d m22_U(2, 2); + m22_U << 0.97759158130035961, 0.21051057021124275, 0.21051057021124264, -0.97759158130035972; + EXPECT_MATRIX_FLOAT_EQ(m22_U, svd_U(m22)); + matrix_d m23(2, 3); m23 << 1, 3, -5, 7, 9, -11; - matrix_d m23_U(2, 2); - m23_U << -0.3378432, -0.9412024, -0.9412024, 0.3378432; //Generated using R base::svd + m23_U << -0.33784321032557019, -0.94120240396894062, -0.94120240396894039, 0.33784321032557019; EXPECT_MATRIX_FLOAT_EQ(m23_U, svd_U(m23)); + + matrix_d m32(3, 2); + m32 << 1, 3, -5, 7, 9, -11; + matrix_d m32_U(3, 2); + m32_U << 0.10657276921942949, -0.978015199679778457, 0.51489182605641137, -0.099934023569151917, + -0.85060487438128174, -0.183028577355020677; + EXPECT_MATRIX_FLOAT_EQ(m32_U, svd_U(m32)); // need to solve sign ambiguity } \ No newline at end of file diff --git a/test/unit/math/prim/fun/svd_V_test.cpp b/test/unit/math/prim/fun/svd_V_test.cpp index ed6010b71b6..78fb8174192 100644 --- a/test/unit/math/prim/fun/svd_V_test.cpp +++ b/test/unit/math/prim/fun/svd_V_test.cpp @@ -1,16 +1,39 @@ #include #include #include +#include TEST(MathMatrixPrimMat, svd_V) { using stan::math::svd_V; using stan::math::matrix_d; + + //Values generated using R base::svd + + matrix_d m00(0, 0); + EXPECT_THROW(svd_V(m00), std::invalid_argument); + + matrix_d m11(1, 1); + m11 << 5; + matrix_d m11_V(1,1); + m11_V << 1; + EXPECT_MATRIX_FLOAT_EQ(m11_V, svd_V(m11)); + + matrix_d m22(2, 2); + m22 << 1, 9, -4, 2; + matrix_d m22_V(2, 2); + m22_V << 0.014701114509569043, 0.999891932776825976, 0.999891932776825976, -0.014701114509569043; + EXPECT_MATRIX_FLOAT_EQ(m22_V, svd_V(m22)); + matrix_d m23(2, 3); m23 << 1, 3, -5, 7, 9, -11; - - matrix_d m23_V(3, 3); - m23_V << -0.41176240532160857, 0.81473005032163681, -0.40824829046386291, - -0.56383954240865775, 0.12417046246885260, 0.81649658092772603, - 0.71591667949570703, 0.56638912538393127, 0.40824829046386307; //Generated using R base::svd + matrix_d m23_V(3, 2); + m23_V << -0.41176240532160857, 0.81473005032163681, -0.56383954240865775, + 0.12417046246885260, 0.71591667949570703, 0.56638912538393127; EXPECT_MATRIX_FLOAT_EQ(m23_V, svd_V(m23)); + + matrix_d m32(3, 2); + m32 << 1, 3, -5, 7, 9, -11; + matrix_d m32_V(2, 2); + m32_V << -0.60622380392317887, 0.79529409626685355, -0.79529409626685355, -0.60622380392317887; + EXPECT_MATRIX_FLOAT_EQ(m32_V, svd_V(m32)); // // need to solve sign ambiguity } From e0574574947100f4a22d4a4748de5bcf8d531c33 Mon Sep 17 00:00:00 2001 From: Dashadower Date: Thu, 31 Dec 2020 18:42:32 +0900 Subject: [PATCH 21/32] Update prim tests --- test/unit/math/mix/fun/svd_U_test.cpp | 5 ----- test/unit/math/prim/fun/svd_U_test.cpp | 6 +++--- test/unit/math/prim/fun/svd_V_test.cpp | 4 ++-- 3 files changed, 5 insertions(+), 10 deletions(-) diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index d8e847dd90e..69c6e8a9b73 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -21,12 +21,7 @@ TEST(MathMixMatFun, svd_U) { Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; Eigen::MatrixXd m32 = m23.transpose(); - /*stan::math::matrix_v test_arr = m32; - std::cout << test_arr.val_op() << std::endl; - std::cout << "----------" << std::endl; - std::cout << f(test_arr) << std:: endl;*/ stan::test::expect_ad(tols, f, m23); - //stan::test::expect_ad(tols, f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; diff --git a/test/unit/math/prim/fun/svd_U_test.cpp b/test/unit/math/prim/fun/svd_U_test.cpp index 6b3eff6e5a7..9794e0eeafd 100644 --- a/test/unit/math/prim/fun/svd_U_test.cpp +++ b/test/unit/math/prim/fun/svd_U_test.cpp @@ -33,7 +33,7 @@ TEST(MathMatrixPrimMat, svd_U) { matrix_d m32(3, 2); m32 << 1, 3, -5, 7, 9, -11; matrix_d m32_U(3, 2); - m32_U << 0.10657276921942949, -0.978015199679778457, 0.51489182605641137, -0.099934023569151917, - -0.85060487438128174, -0.183028577355020677; - EXPECT_MATRIX_FLOAT_EQ(m32_U, svd_U(m32)); // need to solve sign ambiguity + m32_U << 0.10657276921942949, 0.978015199679778457, 0.51489182605641137, 0.099934023569151917, + -0.85060487438128174, 0.183028577355020677; + EXPECT_MATRIX_FLOAT_EQ(m32_U, svd_U(m32)); // R's SVD returns different signs than Eigen. } \ No newline at end of file diff --git a/test/unit/math/prim/fun/svd_V_test.cpp b/test/unit/math/prim/fun/svd_V_test.cpp index 78fb8174192..80a8db55399 100644 --- a/test/unit/math/prim/fun/svd_V_test.cpp +++ b/test/unit/math/prim/fun/svd_V_test.cpp @@ -34,6 +34,6 @@ TEST(MathMatrixPrimMat, svd_V) { matrix_d m32(3, 2); m32 << 1, 3, -5, 7, 9, -11; matrix_d m32_V(2, 2); - m32_V << -0.60622380392317887, 0.79529409626685355, -0.79529409626685355, -0.60622380392317887; - EXPECT_MATRIX_FLOAT_EQ(m32_V, svd_V(m32)); // // need to solve sign ambiguity + m32_V << -0.60622380392317887, 0.79529409626685355, 0.79529409626685355, 0.60622380392317887; + EXPECT_MATRIX_FLOAT_EQ(m32_V, svd_V(m32)); // R's SVD returns different signs than Eigen. } From 63defd1772f4bb5dc68dbe3de683d4961267eaeb Mon Sep 17 00:00:00 2001 From: Dashadower Date: Sat, 2 Jan 2021 20:11:04 +0900 Subject: [PATCH 22/32] Update mix mode tests --- test/unit/math/mix/fun/svd_U_test.cpp | 5 ++++- test/unit/math/mix/fun/svd_V_test.cpp | 4 +++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index 69c6e8a9b73..31bcb84b321 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -20,9 +20,12 @@ TEST(MathMixMatFun, svd_U) { Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; - Eigen::MatrixXd m32 = m23.transpose(); stan::test::expect_ad(tols, f, m23); + Eigen::MatrixXd m32(3, 2); + m32 << 1, 3, -5, 7, 9, -11; + stan::test::expect_ad(tols, f, m32); + Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; stan::test::expect_ad(tols, f, a22); diff --git a/test/unit/math/mix/fun/svd_V_test.cpp b/test/unit/math/mix/fun/svd_V_test.cpp index d1bbdc80f3b..eaaaa9f5550 100644 --- a/test/unit/math/mix/fun/svd_V_test.cpp +++ b/test/unit/math/mix/fun/svd_V_test.cpp @@ -20,8 +20,10 @@ TEST(MathMixMatFun, svd_V) { Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; - Eigen::MatrixXd m32 = m23.transpose(); stan::test::expect_ad(tols, f, m23); + + Eigen::MatrixXd m32(3, 2); + m32 << 1, 3, -5, 7, 9, -11; stan::test::expect_ad(tols, f, m32); Eigen::MatrixXd a22(2, 2); From 7873afcef301354aec11ad1781ffd6cda35c2085 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 2 Jan 2021 13:28:55 +0000 Subject: [PATCH 23/32] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/prim/fun/singular_values.hpp | 5 +- stan/math/prim/fun/svd_U.hpp | 14 ++--- stan/math/prim/fun/svd_V.hpp | 14 ++--- stan/math/rev/fun/singular_values.hpp | 27 +++++---- stan/math/rev/fun/svd_U.hpp | 49 +++++++++------- stan/math/rev/fun/svd_V.hpp | 54 ++++++++++-------- test/unit/math/mix/fun/svd_U_test.cpp | 2 +- test/unit/math/prim/fun/svd_U_test.cpp | 57 ++++++++++--------- test/unit/math/prim/fun/svd_V_test.cpp | 53 +++++++++-------- .../unit/math/rev/fun/singularvalues_test.cpp | 28 ++++----- 10 files changed, 159 insertions(+), 144 deletions(-) diff --git a/stan/math/prim/fun/singular_values.hpp b/stan/math/prim/fun/singular_values.hpp index 85035351e29..d8b7abc39d4 100644 --- a/stan/math/prim/fun/singular_values.hpp +++ b/stan/math/prim/fun/singular_values.hpp @@ -17,9 +17,8 @@ namespace math { * @param m Specified matrix. * @return Singular values of the matrix. */ - template * = nullptr, - require_not_vt_var* = nullptr> +template * = nullptr, + require_not_vt_var* = nullptr> Eigen::Matrix, Eigen::Dynamic, 1> singular_values( const EigMat& m) { if (m.size() == 0) { diff --git a/stan/math/prim/fun/svd_U.hpp b/stan/math/prim/fun/svd_U.hpp index dee604102a8..055ac700b88 100644 --- a/stan/math/prim/fun/svd_U.hpp +++ b/stan/math/prim/fun/svd_U.hpp @@ -14,16 +14,16 @@ namespace math { * @param m MxN input matrix * @return Orthogonal matrix U */ -template * = nullptr, - require_not_vt_var* = nullptr> -Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> -svd_U(const EigMat& m) { +template * = nullptr, + require_not_vt_var* = nullptr> +Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> svd_U( + const EigMat& m) { check_nonzero_size("svd_U", "m", m); return Eigen::JacobiSVD, Eigen::Dynamic, - Eigen::Dynamic> >(m, Eigen::ComputeThinU) - .matrixU(); + Eigen::Dynamic> >(m, + Eigen::ComputeThinU) + .matrixU(); } } // namespace math diff --git a/stan/math/prim/fun/svd_V.hpp b/stan/math/prim/fun/svd_V.hpp index ac07059fce8..485a9af14ab 100644 --- a/stan/math/prim/fun/svd_V.hpp +++ b/stan/math/prim/fun/svd_V.hpp @@ -14,16 +14,16 @@ namespace math { * @param m MxN input matrix * @return Orthogonal matrix V */ -template * = nullptr, - require_not_vt_var* = nullptr> -Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> -svd_V(const EigMat& m) { +template * = nullptr, + require_not_vt_var* = nullptr> +Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> svd_V( + const EigMat& m) { check_nonzero_size("svd_V", "m", m); return Eigen::JacobiSVD, Eigen::Dynamic, - Eigen::Dynamic> >(m, Eigen::ComputeThinV) - .matrixV(); + Eigen::Dynamic> >(m, + Eigen::ComputeThinV) + .matrixV(); } } // namespace math diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index eb3970e73d5..b0e5f9e35f3 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -13,37 +13,36 @@ namespace math { /** * Return the singular values of the specified matrix. * - * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650) - * is used for MxN input matrix's adjoint calculation. + * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, + * et al., arXiv:1903.09650) is used for MxN input matrix's adjoint calculation. * * @tparam EigMat type of input matrix * @param m MxN input matrix * @return Singular values of matrix */ -template * = nullptr, - require_vt_var* = nullptr> +template * = nullptr, + require_vt_var* = nullptr> inline auto singular_values(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("singular_values", "m", m); const int M = to_arena(m.rows()); auto arena_m = to_arena(m); - - Eigen::JacobiSVD svd(arena_m.val(), - Eigen::ComputeFullU | Eigen::ComputeFullV); + + Eigen::JacobiSVD svd( + arena_m.val(), Eigen::ComputeFullU | Eigen::ComputeFullV); arena_t singular_values = svd.singularValues(); auto U = to_arena(svd.matrixU()); auto V = to_arena(svd.matrixV()); - reverse_pass_callback( - [arena_m, U, singular_values, V]() mutable { - arena_m.adj() += U.block(0, 0, U.rows(), singular_values.size()) * - singular_values.adj().asDiagonal() * - V.block(0, 0, V.rows(), singular_values.size()).transpose(); - }); + reverse_pass_callback([arena_m, U, singular_values, V]() mutable { + arena_m.adj() + += U.block(0, 0, U.rows(), singular_values.size()) + * singular_values.adj().asDiagonal() + * V.block(0, 0, V.rows(), singular_values.size()).transpose(); + }); return ret_type(singular_values); } diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index 0cf4bbea75b..816e7769a5c 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -19,52 +19,57 @@ namespace math { * * This function internally calls Eigen::JacobiSVD to compute decomposition. * - * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650) - * is used for MxN input matrix's adjoint calculation. + * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, + * et al., arXiv:1903.09650) is used for MxN input matrix's adjoint calculation. * * @tparam EigMat type of input matrix * @param m MxN input matrix * @return Orthogonal matrix U */ -template * = nullptr, - require_vt_var* = nullptr> +template * = nullptr, + require_vt_var* = nullptr> inline auto svd_U(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("svd_U", "m", m); const int M = to_arena(std::min(m.rows(), m.cols())); - auto arena_m = to_arena(m); // N by P - - Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); + auto arena_m = to_arena(m); // N by P - auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M + Eigen::JacobiSVD svd( + arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); + + auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M Eigen::MatrixXd Fp(M, M); - for(int i = 0; i < M; i++) { - for(int j = 0; j < M; j++) { - if(j == i) { + for (int i = 0; i < M; i++) { + for (int j = 0; j < M; j++) { + if (j == i) { Fp(i, j) = 0.0; } else { - Fp(i, j) = 1.0 / (arena_D[j] - arena_D[i]) + 1.0 / (arena_D[i] + arena_D[j]); + Fp(i, j) + = 1.0 / (arena_D[j] - arena_D[i]) + 1.0 / (arena_D[i] + arena_D[j]); } } } auto arena_Fp = to_arena(Fp); - arena_t arena_U = svd.matrixU(); // N by M - auto arena_V = to_arena(svd.matrixV()); // P by M - + arena_t arena_U = svd.matrixU(); // N by M + auto arena_V = to_arena(svd.matrixV()); // P by M reverse_pass_callback( [arena_m, arena_U, arena_D, arena_V, arena_Fp, M]() mutable { - Eigen::MatrixXd UUadjT = arena_U.val_op().transpose() * arena_U.adj_op(); // M by M - arena_m.adj() += .5 * arena_U.val_op() * - (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()).matrix() - * arena_V.transpose() + - (Eigen::MatrixXd::Identity(arena_m.rows(), arena_m.rows()) - arena_U.val_op() * arena_U.val_op().transpose()) * - arena_U.adj_op() * arena_D.asDiagonal().inverse() * arena_V.transpose(); + Eigen::MatrixXd UUadjT + = arena_U.val_op().transpose() * arena_U.adj_op(); // M by M + arena_m.adj() + += .5 * arena_U.val_op() + * (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()) + .matrix() + * arena_V.transpose() + + (Eigen::MatrixXd::Identity(arena_m.rows(), arena_m.rows()) + - arena_U.val_op() * arena_U.val_op().transpose()) + * arena_U.adj_op() * arena_D.asDiagonal().inverse() + * arena_V.transpose(); }); return ret_type(arena_U); diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index a92776d795e..00d68f03e8f 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -19,52 +19,58 @@ namespace math { * * This function internally calls Eigen::JacobiSVD to compute decomposition. * - * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650) - * is used for MxN input matrix's adjoint calculation. + * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, + * et al., arXiv:1903.09650) is used for MxN input matrix's adjoint calculation. * * @tparam EigMat type of input matrix * @param m MxN input matrix * @return Orthogonal matrix V */ -template * = nullptr, - require_vt_var* = nullptr> +template * = nullptr, + require_vt_var* = nullptr> inline auto svd_V(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("svd_V", "m", m); const int M = to_arena(std::min(m.rows(), m.cols())); - auto arena_m = to_arena(m); // N by P - - Eigen::JacobiSVD svd(arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); + auto arena_m = to_arena(m); // N by P - auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M + Eigen::JacobiSVD svd( + arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); + + auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M Eigen::MatrixXd Fm(M, M); - for(int i = 0; i < M; i++) { - for(int j = 0; j < M; j++) { - if(j == i) { + for (int i = 0; i < M; i++) { + for (int j = 0; j < M; j++) { + if (j == i) { Fm(i, j) = 0.0; } else { - Fm(i, j) = 1.0 / (arena_D[j] - arena_D[i]) - 1.0 / (arena_D[i] + arena_D[j]); + Fm(i, j) + = 1.0 / (arena_D[j] - arena_D[i]) - 1.0 / (arena_D[i] + arena_D[j]); } } } auto arena_Fm = to_arena(Fm); - auto arena_U = to_arena(svd.matrixU()); // N by M - arena_t arena_V = svd.matrixV(); // P by M - + auto arena_U = to_arena(svd.matrixU()); // N by M + arena_t arena_V = svd.matrixV(); // P by M - reverse_pass_callback( - [arena_m, arena_U, arena_D, arena_V, arena_Fm, M]() mutable { - Eigen::MatrixXd VTVadj = arena_V.val_op().transpose() * arena_V.adj_op(); // M by M - arena_m.adj() += 0.5 * arena_U * - (arena_Fm.array() * (VTVadj - VTVadj.transpose()).array()).matrix() * - arena_V.val_op().transpose() + arena_U * arena_D.asDiagonal().inverse() * arena_V.adj_op().transpose() * - (Eigen::MatrixXd::Identity(arena_m.cols(), arena_m.cols()) - arena_V.val_op() * arena_V.val_op().transpose()); - }); + reverse_pass_callback([arena_m, arena_U, arena_D, arena_V, arena_Fm, + M]() mutable { + Eigen::MatrixXd VTVadj + = arena_V.val_op().transpose() * arena_V.adj_op(); // M by M + arena_m.adj() + += 0.5 * arena_U + * (arena_Fm.array() * (VTVadj - VTVadj.transpose()).array()) + .matrix() + * arena_V.val_op().transpose() + + arena_U * arena_D.asDiagonal().inverse() + * arena_V.adj_op().transpose() + * (Eigen::MatrixXd::Identity(arena_m.cols(), arena_m.cols()) + - arena_V.val_op() * arena_V.val_op().transpose()); + }); return ret_type(arena_V); } diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index 31bcb84b321..386d04df2ee 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -17,7 +17,7 @@ TEST(MathMixMatFun, svd_U) { Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; stan::test::expect_ad(tols, f, m22); - + Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; stan::test::expect_ad(tols, f, m23); diff --git a/test/unit/math/prim/fun/svd_U_test.cpp b/test/unit/math/prim/fun/svd_U_test.cpp index 9794e0eeafd..1fbd16ed782 100644 --- a/test/unit/math/prim/fun/svd_U_test.cpp +++ b/test/unit/math/prim/fun/svd_U_test.cpp @@ -4,36 +4,39 @@ #include TEST(MathMatrixPrimMat, svd_U) { - using stan::math::svd_U; - using stan::math::matrix_d; + using stan::math::matrix_d; + using stan::math::svd_U; - //Values generated using R base::svd + // Values generated using R base::svd - matrix_d m00(0, 0); - EXPECT_THROW(svd_U(m00), std::invalid_argument); + matrix_d m00(0, 0); + EXPECT_THROW(svd_U(m00), std::invalid_argument); - matrix_d m11(1, 1); - m11 << 5; - matrix_d m11_U(1,1); - m11_U << 1; - EXPECT_MATRIX_FLOAT_EQ(m11_U, svd_U(m11)); + matrix_d m11(1, 1); + m11 << 5; + matrix_d m11_U(1, 1); + m11_U << 1; + EXPECT_MATRIX_FLOAT_EQ(m11_U, svd_U(m11)); - matrix_d m22(2, 2); - m22 << 1, 9, -4, 2; - matrix_d m22_U(2, 2); - m22_U << 0.97759158130035961, 0.21051057021124275, 0.21051057021124264, -0.97759158130035972; - EXPECT_MATRIX_FLOAT_EQ(m22_U, svd_U(m22)); - - matrix_d m23(2, 3); - m23 << 1, 3, -5, 7, 9, -11; - matrix_d m23_U(2, 2); - m23_U << -0.33784321032557019, -0.94120240396894062, -0.94120240396894039, 0.33784321032557019; - EXPECT_MATRIX_FLOAT_EQ(m23_U, svd_U(m23)); + matrix_d m22(2, 2); + m22 << 1, 9, -4, 2; + matrix_d m22_U(2, 2); + m22_U << 0.97759158130035961, 0.21051057021124275, 0.21051057021124264, + -0.97759158130035972; + EXPECT_MATRIX_FLOAT_EQ(m22_U, svd_U(m22)); - matrix_d m32(3, 2); - m32 << 1, 3, -5, 7, 9, -11; - matrix_d m32_U(3, 2); - m32_U << 0.10657276921942949, 0.978015199679778457, 0.51489182605641137, 0.099934023569151917, - -0.85060487438128174, 0.183028577355020677; - EXPECT_MATRIX_FLOAT_EQ(m32_U, svd_U(m32)); // R's SVD returns different signs than Eigen. + matrix_d m23(2, 3); + m23 << 1, 3, -5, 7, 9, -11; + matrix_d m23_U(2, 2); + m23_U << -0.33784321032557019, -0.94120240396894062, -0.94120240396894039, + 0.33784321032557019; + EXPECT_MATRIX_FLOAT_EQ(m23_U, svd_U(m23)); + + matrix_d m32(3, 2); + m32 << 1, 3, -5, 7, 9, -11; + matrix_d m32_U(3, 2); + m32_U << 0.10657276921942949, 0.978015199679778457, 0.51489182605641137, + 0.099934023569151917, -0.85060487438128174, 0.183028577355020677; + EXPECT_MATRIX_FLOAT_EQ( + m32_U, svd_U(m32)); // R's SVD returns different signs than Eigen. } \ No newline at end of file diff --git a/test/unit/math/prim/fun/svd_V_test.cpp b/test/unit/math/prim/fun/svd_V_test.cpp index 80a8db55399..42dc5e85be4 100644 --- a/test/unit/math/prim/fun/svd_V_test.cpp +++ b/test/unit/math/prim/fun/svd_V_test.cpp @@ -4,36 +4,39 @@ #include TEST(MathMatrixPrimMat, svd_V) { - using stan::math::svd_V; - using stan::math::matrix_d; + using stan::math::matrix_d; + using stan::math::svd_V; - //Values generated using R base::svd + // Values generated using R base::svd - matrix_d m00(0, 0); - EXPECT_THROW(svd_V(m00), std::invalid_argument); + matrix_d m00(0, 0); + EXPECT_THROW(svd_V(m00), std::invalid_argument); - matrix_d m11(1, 1); - m11 << 5; - matrix_d m11_V(1,1); - m11_V << 1; - EXPECT_MATRIX_FLOAT_EQ(m11_V, svd_V(m11)); + matrix_d m11(1, 1); + m11 << 5; + matrix_d m11_V(1, 1); + m11_V << 1; + EXPECT_MATRIX_FLOAT_EQ(m11_V, svd_V(m11)); - matrix_d m22(2, 2); - m22 << 1, 9, -4, 2; - matrix_d m22_V(2, 2); - m22_V << 0.014701114509569043, 0.999891932776825976, 0.999891932776825976, -0.014701114509569043; - EXPECT_MATRIX_FLOAT_EQ(m22_V, svd_V(m22)); + matrix_d m22(2, 2); + m22 << 1, 9, -4, 2; + matrix_d m22_V(2, 2); + m22_V << 0.014701114509569043, 0.999891932776825976, 0.999891932776825976, + -0.014701114509569043; + EXPECT_MATRIX_FLOAT_EQ(m22_V, svd_V(m22)); - matrix_d m23(2, 3); - m23 << 1, 3, -5, 7, 9, -11; - matrix_d m23_V(3, 2); - m23_V << -0.41176240532160857, 0.81473005032163681, -0.56383954240865775, + matrix_d m23(2, 3); + m23 << 1, 3, -5, 7, 9, -11; + matrix_d m23_V(3, 2); + m23_V << -0.41176240532160857, 0.81473005032163681, -0.56383954240865775, 0.12417046246885260, 0.71591667949570703, 0.56638912538393127; - EXPECT_MATRIX_FLOAT_EQ(m23_V, svd_V(m23)); + EXPECT_MATRIX_FLOAT_EQ(m23_V, svd_V(m23)); - matrix_d m32(3, 2); - m32 << 1, 3, -5, 7, 9, -11; - matrix_d m32_V(2, 2); - m32_V << -0.60622380392317887, 0.79529409626685355, 0.79529409626685355, 0.60622380392317887; - EXPECT_MATRIX_FLOAT_EQ(m32_V, svd_V(m32)); // R's SVD returns different signs than Eigen. + matrix_d m32(3, 2); + m32 << 1, 3, -5, 7, 9, -11; + matrix_d m32_V(2, 2); + m32_V << -0.60622380392317887, 0.79529409626685355, 0.79529409626685355, + 0.60622380392317887; + EXPECT_MATRIX_FLOAT_EQ( + m32_V, svd_V(m32)); // R's SVD returns different signs than Eigen. } diff --git a/test/unit/math/rev/fun/singularvalues_test.cpp b/test/unit/math/rev/fun/singularvalues_test.cpp index 1d63699a043..d2ec0d105eb 100644 --- a/test/unit/math/rev/fun/singularvalues_test.cpp +++ b/test/unit/math/rev/fun/singularvalues_test.cpp @@ -5,21 +5,21 @@ #include TEST(AgradRev, singularvalues_gradient) { -// logdet(A) can be calculated using singularvalues of matrix A -// the derivative of logdet(A) should be inverse(A) + // logdet(A) can be calculated using singularvalues of matrix A + // the derivative of logdet(A) should be inverse(A) -Eigen::MatrixXd a(4, 4); -// Random symmetric matrix -a << 1.8904, 0.7204, -0.1599, 1.2028, 0.7204, 7.3394, 2.0895, -0.6151, --0.1599, 2.0895, 2.4601, -1.7219, 1.2028, -0.6151, -1.7219, 4.5260; -stan::math::matrix_d a_inv = stan::math::inverse(a); + Eigen::MatrixXd a(4, 4); + // Random symmetric matrix + a << 1.8904, 0.7204, -0.1599, 1.2028, 0.7204, 7.3394, 2.0895, -0.6151, + -0.1599, 2.0895, 2.4601, -1.7219, 1.2028, -0.6151, -1.7219, 4.5260; + stan::math::matrix_d a_inv = stan::math::inverse(a); -stan::math::matrix_v a_v(a); -auto w = stan::math::singular_values(a_v); -auto logdet = stan::math::sum(stan::math::log(w)); + stan::math::matrix_v a_v(a); + auto w = stan::math::singular_values(a_v); + auto logdet = stan::math::sum(stan::math::log(w)); -stan::math::set_zero_all_adjoints(); -logdet.grad(); -std::cout<<"a_inv.val()" < Date: Sun, 3 Jan 2021 01:48:32 +0900 Subject: [PATCH 24/32] Fix lint warnings --- test/unit/math/mix/fun/svd_U_test.cpp | 1 - test/unit/math/prim/fun/svd_U_test.cpp | 2 +- test/unit/math/rev/fun/singularvalues_test.cpp | 2 +- 3 files changed, 2 insertions(+), 3 deletions(-) diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index 386d04df2ee..3a131518224 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -1,5 +1,4 @@ #include -//#include TEST(MathMixMatFun, svd_U) { auto f = [](const auto& x) { return stan::math::svd_U(x); }; diff --git a/test/unit/math/prim/fun/svd_U_test.cpp b/test/unit/math/prim/fun/svd_U_test.cpp index 1fbd16ed782..a0bec909f81 100644 --- a/test/unit/math/prim/fun/svd_U_test.cpp +++ b/test/unit/math/prim/fun/svd_U_test.cpp @@ -39,4 +39,4 @@ TEST(MathMatrixPrimMat, svd_U) { 0.099934023569151917, -0.85060487438128174, 0.183028577355020677; EXPECT_MATRIX_FLOAT_EQ( m32_U, svd_U(m32)); // R's SVD returns different signs than Eigen. -} \ No newline at end of file +} diff --git a/test/unit/math/rev/fun/singularvalues_test.cpp b/test/unit/math/rev/fun/singularvalues_test.cpp index d2ec0d105eb..5e3e54e466e 100644 --- a/test/unit/math/rev/fun/singularvalues_test.cpp +++ b/test/unit/math/rev/fun/singularvalues_test.cpp @@ -22,4 +22,4 @@ TEST(AgradRev, singularvalues_gradient) { logdet.grad(); std::cout << "a_inv.val()" << a_inv.val() << "a_v.adj()" << a_v.adj(); ASSERT_TRUE(a_inv.val().isApprox(a_v.adj())); -} \ No newline at end of file +} From 52889c4d0c9ef2aaf185f5df4ec23c6f7ffd54ef Mon Sep 17 00:00:00 2001 From: Dashadower Date: Sun, 3 Jan 2021 02:04:16 +0900 Subject: [PATCH 25/32] Update doxygen warnings --- stan/math/prim/fun/svd_U.hpp | 2 +- stan/math/prim/fun/svd_V.hpp | 2 +- stan/math/rev/fun/svd_U.hpp | 2 +- stan/math/rev/fun/svd_V.hpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/stan/math/prim/fun/svd_U.hpp b/stan/math/prim/fun/svd_U.hpp index 055ac700b88..5001349f983 100644 --- a/stan/math/prim/fun/svd_U.hpp +++ b/stan/math/prim/fun/svd_U.hpp @@ -8,7 +8,7 @@ namespace stan { namespace math { /** - * Given input matrix m, return matrix U where \mathbf{m} = \mathbf{UDV^{T}} + * Given input matrix m, return matrix U where `m = UDV^{T}` * * @tparam EigMat type of the matrix * @param m MxN input matrix diff --git a/stan/math/prim/fun/svd_V.hpp b/stan/math/prim/fun/svd_V.hpp index 485a9af14ab..9b2611209ac 100644 --- a/stan/math/prim/fun/svd_V.hpp +++ b/stan/math/prim/fun/svd_V.hpp @@ -8,7 +8,7 @@ namespace stan { namespace math { /** - * Given input matrix m, return matrix V where \mathbf{m} = \mathbf{UDV^{T}} + * Given input matrix m, return matrix V where `m = UDV^{T}` * * @tparam EigMat type of the matrix * @param m MxN input matrix diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index 816e7769a5c..8077d15b848 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -15,7 +15,7 @@ namespace stan { namespace math { /** - * Given input matrix m, return matrix U where \mathbf{m} = \mathbf{UDV^{T}} + * Given input matrix m, return matrix U where `m = UDV^{T}` * * This function internally calls Eigen::JacobiSVD to compute decomposition. * diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index 00d68f03e8f..dad1aaa55f0 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -15,7 +15,7 @@ namespace stan { namespace math { /** - * Given input matrix m, return matrix V where \mathbf{m} = \mathbf{UDV^{T}} + * Given input matrix m, return matrix V where \f$\mathbf{m} = \mathbf{UDV^{T}}\f$ * * This function internally calls Eigen::JacobiSVD to compute decomposition. * From c4cb965173744950081fb21f390223752d647043 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sat, 2 Jan 2021 17:11:08 +0000 Subject: [PATCH 26/32] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/svd_V.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index dad1aaa55f0..44b5bc353d5 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -15,7 +15,8 @@ namespace stan { namespace math { /** - * Given input matrix m, return matrix V where \f$\mathbf{m} = \mathbf{UDV^{T}}\f$ + * Given input matrix m, return matrix V where \f$\mathbf{m} = + * \mathbf{UDV^{T}}\f$ * * This function internally calls Eigen::JacobiSVD to compute decomposition. * From 8d41dd46818c9d7e22a9b07cb833fbd5045be4fe Mon Sep 17 00:00:00 2001 From: Dashadower Date: Sun, 3 Jan 2021 02:23:14 +0900 Subject: [PATCH 27/32] Update doxygen docs --- stan/math/rev/fun/svd_V.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index 44b5bc353d5..22b88fcc23c 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -15,8 +15,7 @@ namespace stan { namespace math { /** - * Given input matrix m, return matrix V where \f$\mathbf{m} = - * \mathbf{UDV^{T}}\f$ + * Given input matrix m, return matrix V where `m = UDV^{T}` * * This function internally calls Eigen::JacobiSVD to compute decomposition. * From e8713e4757a6533e52e93a1da22c1c8a76f62e94 Mon Sep 17 00:00:00 2001 From: Dashadower Date: Tue, 5 Jan 2021 23:38:23 +0900 Subject: [PATCH 28/32] Update tests and functions for svd --- stan/math/prim/fun/Eigen.hpp | 1 + stan/math/rev/fun/singular_values.hpp | 17 ++++++------- stan/math/rev/fun/svd_U.hpp | 18 +++++-------- stan/math/rev/fun/svd_V.hpp | 17 +++++-------- .../math/prim/fun/singular_values_test.cpp | 25 +++++++++++++++++-- ...lues_test.cpp => singular_values_test.cpp} | 1 - 6 files changed, 43 insertions(+), 36 deletions(-) rename test/unit/math/rev/fun/{singularvalues_test.cpp => singular_values_test.cpp} (91%) diff --git a/stan/math/prim/fun/Eigen.hpp b/stan/math/prim/fun/Eigen.hpp index b1071b49879..6efb32cbdea 100644 --- a/stan/math/prim/fun/Eigen.hpp +++ b/stan/math/prim/fun/Eigen.hpp @@ -23,6 +23,7 @@ #include #include #include +#include namespace Eigen { diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index b0e5f9e35f3..e03baf3d757 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -13,8 +13,8 @@ namespace math { /** * Return the singular values of the specified matrix. * - * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, - * et al., arXiv:1903.09650) is used for MxN input matrix's adjoint calculation. + * Adjoint update equation comes from Equation (4) in Differentiable Programming + * Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650). * * @tparam EigMat type of input matrix * @param m MxN input matrix @@ -26,22 +26,19 @@ inline auto singular_values(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("singular_values", "m", m); - const int M = to_arena(m.rows()); auto arena_m = to_arena(m); Eigen::JacobiSVD svd( - arena_m.val(), Eigen::ComputeFullU | Eigen::ComputeFullV); + arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); arena_t singular_values = svd.singularValues(); - auto U = to_arena(svd.matrixU()); - auto V = to_arena(svd.matrixV()); + auto arena_U = to_arena(svd.matrixU()); + auto arena_V = to_arena(svd.matrixV()); - reverse_pass_callback([arena_m, U, singular_values, V]() mutable { + reverse_pass_callback([arena_m, arena_U, singular_values, arena_V]() mutable { arena_m.adj() - += U.block(0, 0, U.rows(), singular_values.size()) - * singular_values.adj().asDiagonal() - * V.block(0, 0, V.rows(), singular_values.size()).transpose(); + += arena_U * singular_values.adj().asDiagonal() * arena_V.transpose(); }); return ret_type(singular_values); diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index 8077d15b848..e7234187a9e 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -5,11 +5,8 @@ #include #include #include -#include #include #include -#include -#include namespace stan { namespace math { @@ -17,10 +14,8 @@ namespace math { /** * Given input matrix m, return matrix U where `m = UDV^{T}` * - * This function internally calls Eigen::JacobiSVD to compute decomposition. - * - * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, - * et al., arXiv:1903.09650) is used for MxN input matrix's adjoint calculation. + * Adjoint update equation comes from Equation (4) in Differentiable Programming + * Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650). * * @tparam EigMat type of input matrix * @param m MxN input matrix @@ -32,7 +27,7 @@ inline auto svd_U(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("svd_U", "m", m); - const int M = to_arena(std::min(m.rows(), m.cols())); + const int M = std::min(m.rows(), m.cols()); auto arena_m = to_arena(m); // N by P Eigen::JacobiSVD svd( @@ -40,20 +35,19 @@ inline auto svd_U(const EigMat& m) { auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M - Eigen::MatrixXd Fp(M, M); + arena_t arena_Fp(M, M); for (int i = 0; i < M; i++) { for (int j = 0; j < M; j++) { if (j == i) { - Fp(i, j) = 0.0; + arena_Fp(i, j) = 0.0; } else { - Fp(i, j) + arena_Fp(i, j) = 1.0 / (arena_D[j] - arena_D[i]) + 1.0 / (arena_D[i] + arena_D[j]); } } } - auto arena_Fp = to_arena(Fp); arena_t arena_U = svd.matrixU(); // N by M auto arena_V = to_arena(svd.matrixV()); // P by M diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index 22b88fcc23c..65f8c17d432 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -5,10 +5,8 @@ #include #include #include -#include #include #include -#include #include namespace stan { @@ -17,10 +15,8 @@ namespace math { /** * Given input matrix m, return matrix V where `m = UDV^{T}` * - * This function internally calls Eigen::JacobiSVD to compute decomposition. - * - * Equation (4) from Differentiable Programming Tensor Networks(H. Liao, J. Liu, - * et al., arXiv:1903.09650) is used for MxN input matrix's adjoint calculation. + * Adjoint update equation comes from Equation (4) in Differentiable Programming + * Tensor Networks(H. Liao, J. Liu, et al., arXiv:1903.09650). * * @tparam EigMat type of input matrix * @param m MxN input matrix @@ -32,7 +28,7 @@ inline auto svd_V(const EigMat& m) { using ret_type = promote_scalar_t; check_nonzero_size("svd_V", "m", m); - const int M = to_arena(std::min(m.rows(), m.cols())); + const int M = std::min(m.rows(), m.cols()); auto arena_m = to_arena(m); // N by P Eigen::JacobiSVD svd( @@ -40,20 +36,19 @@ inline auto svd_V(const EigMat& m) { auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M - Eigen::MatrixXd Fm(M, M); + arena_t arena_Fm(M, M); for (int i = 0; i < M; i++) { for (int j = 0; j < M; j++) { if (j == i) { - Fm(i, j) = 0.0; + arena_Fm(i, j) = 0.0; } else { - Fm(i, j) + arena_Fm(i, j) = 1.0 / (arena_D[j] - arena_D[i]) - 1.0 / (arena_D[i] + arena_D[j]); } } } - auto arena_Fm = to_arena(Fm); auto arena_U = to_arena(svd.matrixU()); // N by M arena_t arena_V = svd.matrixV(); // P by M diff --git a/test/unit/math/prim/fun/singular_values_test.cpp b/test/unit/math/prim/fun/singular_values_test.cpp index 229a174cb03..2597e95474a 100644 --- a/test/unit/math/prim/fun/singular_values_test.cpp +++ b/test/unit/math/prim/fun/singular_values_test.cpp @@ -1,13 +1,34 @@ #include #include +#include TEST(MathMatrixPrimMat, singular_values) { using stan::math::singular_values; + using stan::math::matrix_d; + using stan::math::vector_d; - stan::math::matrix_d m0(0, 0); + matrix_d m0(0, 0); EXPECT_NO_THROW(singular_values(m0)); - stan::math::matrix_d m1(1, 1); + matrix_d m1(1, 1); m1 << 1.0; EXPECT_NO_THROW(singular_values(m1)); + + matrix_d m22(2, 2); + m22 << 1, 9, -4, 2; + vector_d m22_D(2); + m22_D << 9.220341788859560, 4.121322275266775; + EXPECT_MATRIX_FLOAT_EQ(m22_D, singular_values(m22)); + + matrix_d m23(2, 3); + m23 << 1, 3, -5, 7, 9, -11; + vector_d m23_D(2); + m23_D << 16.821011215675149, 1.747450051398016; + EXPECT_MATRIX_FLOAT_EQ(m23_D, singular_values(m23)); + + matrix_d m32(3, 2); + m32 << 1, 3, -5, 7, 9, -11; + vector_d m32_D(2); + m32_D << 16.698998232964481, 2.672724829728879; + EXPECT_MATRIX_FLOAT_EQ(m32_D, singular_values(m32)); } diff --git a/test/unit/math/rev/fun/singularvalues_test.cpp b/test/unit/math/rev/fun/singular_values_test.cpp similarity index 91% rename from test/unit/math/rev/fun/singularvalues_test.cpp rename to test/unit/math/rev/fun/singular_values_test.cpp index 5e3e54e466e..ea53205e884 100644 --- a/test/unit/math/rev/fun/singularvalues_test.cpp +++ b/test/unit/math/rev/fun/singular_values_test.cpp @@ -20,6 +20,5 @@ TEST(AgradRev, singularvalues_gradient) { stan::math::set_zero_all_adjoints(); logdet.grad(); - std::cout << "a_inv.val()" << a_inv.val() << "a_v.adj()" << a_v.adj(); ASSERT_TRUE(a_inv.val().isApprox(a_v.adj())); } From f47682b91bd07b453e229ef57024f5a84b8a376d Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Tue, 5 Jan 2021 14:44:36 +0000 Subject: [PATCH 29/32] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/svd_U.hpp | 2 +- test/unit/math/prim/fun/singular_values_test.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index e7234187a9e..a09d5b8c52b 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -35,7 +35,7 @@ inline auto svd_U(const EigMat& m) { auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M - arena_t arena_Fp(M, M); + arena_t arena_Fp(M, M); for (int i = 0; i < M; i++) { for (int j = 0; j < M; j++) { diff --git a/test/unit/math/prim/fun/singular_values_test.cpp b/test/unit/math/prim/fun/singular_values_test.cpp index 2597e95474a..5877b0c6279 100644 --- a/test/unit/math/prim/fun/singular_values_test.cpp +++ b/test/unit/math/prim/fun/singular_values_test.cpp @@ -3,8 +3,8 @@ #include TEST(MathMatrixPrimMat, singular_values) { - using stan::math::singular_values; using stan::math::matrix_d; + using stan::math::singular_values; using stan::math::vector_d; matrix_d m0(0, 0); From c3835c687a5b3b2f0b7adde80039a1a062c12997 Mon Sep 17 00:00:00 2001 From: Dashadower Date: Mon, 11 Jan 2021 00:53:23 +0900 Subject: [PATCH 30/32] Add vari support and relevant tests --- stan/math/prim/fun/singular_values.hpp | 7 +++---- stan/math/prim/fun/svd_U.hpp | 2 +- stan/math/prim/fun/svd_V.hpp | 2 +- stan/math/rev/fun/singular_values.hpp | 5 ++--- stan/math/rev/fun/svd_U.hpp | 5 ++--- stan/math/rev/fun/svd_V.hpp | 5 ++--- test/unit/math/mix/fun/singular_values_test.cpp | 8 +++++++- test/unit/math/mix/fun/svd_U_test.cpp | 9 ++++++++- test/unit/math/mix/fun/svd_V_test.cpp | 8 +++++++- test/unit/math/prim/fun/singular_values_test.cpp | 3 ++- 10 files changed, 35 insertions(+), 19 deletions(-) diff --git a/stan/math/prim/fun/singular_values.hpp b/stan/math/prim/fun/singular_values.hpp index d8b7abc39d4..c177c8867d9 100644 --- a/stan/math/prim/fun/singular_values.hpp +++ b/stan/math/prim/fun/singular_values.hpp @@ -3,6 +3,7 @@ #include #include +#include namespace stan { namespace math { @@ -18,12 +19,10 @@ namespace math { * @return Singular values of the matrix. */ template * = nullptr, - require_not_vt_var* = nullptr> + require_not_st_var* = nullptr> Eigen::Matrix, Eigen::Dynamic, 1> singular_values( const EigMat& m) { - if (m.size() == 0) { - return {}; - } + check_nonzero_size("singular_values", "m", m); return Eigen::JacobiSVD, Eigen::Dynamic, Eigen::Dynamic> >(m) diff --git a/stan/math/prim/fun/svd_U.hpp b/stan/math/prim/fun/svd_U.hpp index 5001349f983..b69d7d9566d 100644 --- a/stan/math/prim/fun/svd_U.hpp +++ b/stan/math/prim/fun/svd_U.hpp @@ -15,7 +15,7 @@ namespace math { * @return Orthogonal matrix U */ template * = nullptr, - require_not_vt_var* = nullptr> + require_not_st_var* = nullptr> Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> svd_U( const EigMat& m) { check_nonzero_size("svd_U", "m", m); diff --git a/stan/math/prim/fun/svd_V.hpp b/stan/math/prim/fun/svd_V.hpp index 9b2611209ac..ef85274acd7 100644 --- a/stan/math/prim/fun/svd_V.hpp +++ b/stan/math/prim/fun/svd_V.hpp @@ -15,7 +15,7 @@ namespace math { * @return Orthogonal matrix V */ template * = nullptr, - require_not_vt_var* = nullptr> + require_not_st_var* = nullptr> Eigen::Matrix, Eigen::Dynamic, Eigen::Dynamic> svd_V( const EigMat& m) { check_nonzero_size("svd_V", "m", m); diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index e03baf3d757..bd3fe7a03b7 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -20,10 +20,9 @@ namespace math { * @param m MxN input matrix * @return Singular values of matrix */ -template * = nullptr, - require_vt_var* = nullptr> +template * = nullptr> inline auto singular_values(const EigMat& m) { - using ret_type = promote_scalar_t; + using ret_type = return_var_matrix_t; check_nonzero_size("singular_values", "m", m); auto arena_m = to_arena(m); diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index a09d5b8c52b..d984c57cfa4 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -21,10 +21,9 @@ namespace math { * @param m MxN input matrix * @return Orthogonal matrix U */ -template * = nullptr, - require_vt_var* = nullptr> +template * = nullptr> inline auto svd_U(const EigMat& m) { - using ret_type = promote_scalar_t; + using ret_type = return_var_matrix_t; check_nonzero_size("svd_U", "m", m); const int M = std::min(m.rows(), m.cols()); diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index 65f8c17d432..3ec9800337a 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -22,10 +22,9 @@ namespace math { * @param m MxN input matrix * @return Orthogonal matrix V */ -template * = nullptr, - require_vt_var* = nullptr> +template * = nullptr> inline auto svd_V(const EigMat& m) { - using ret_type = promote_scalar_t; + using ret_type = return_var_matrix_t; check_nonzero_size("svd_V", "m", m); const int M = std::min(m.rows(), m.cols()); diff --git a/test/unit/math/mix/fun/singular_values_test.cpp b/test/unit/math/mix/fun/singular_values_test.cpp index d6b757035ef..76c88205f71 100644 --- a/test/unit/math/mix/fun/singular_values_test.cpp +++ b/test/unit/math/mix/fun/singular_values_test.cpp @@ -1,4 +1,5 @@ #include +#include TEST(MathMixMatFun, singularValues) { auto f = [](const auto& x) { return stan::math::singular_values(x); }; @@ -8,23 +9,28 @@ TEST(MathMixMatFun, singularValues) { tols.hessian_fvar_hessian_ = 1e-2; Eigen::MatrixXd m00(0, 0); - stan::test::expect_ad(f, m00); + EXPECT_THROW(f(m00), std::invalid_argument); Eigen::MatrixXd m11(1, 1); m11 << 1.1; stan::test::expect_ad(tols, f, m11); + stan::test::expect_ad_matvar(f, m11); Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; stan::test::expect_ad(tols, f, m22); + stan::test::expect_ad_matvar(f, m22); Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; Eigen::MatrixXd m32 = m23.transpose(); stan::test::expect_ad(tols, f, m23); stan::test::expect_ad(tols, f, m32); + stan::test::expect_ad_matvar(f, m23); + stan::test::expect_ad_matvar(f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; stan::test::expect_ad(tols, f, a22); + stan::test::expect_ad_matvar(f, a22); } diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index 3a131518224..fd3a8389861 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -1,4 +1,6 @@ #include +#include + TEST(MathMixMatFun, svd_U) { auto f = [](const auto& x) { return stan::math::svd_U(x); }; @@ -7,25 +9,30 @@ TEST(MathMixMatFun, svd_U) { tols.hessian_fvar_hessian_ = 1e-2; Eigen::MatrixXd m00(0, 0); - stan::test::expect_ad(f, m00); + EXPECT_THROW(f(m00), std::invalid_argument); Eigen::MatrixXd m11(1, 1); m11 << 1.1; stan::test::expect_ad(tols, f, m11); + stan::test::expect_ad_matvar(f, m11); Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; stan::test::expect_ad(tols, f, m22); + stan::test::expect_ad_matvar(f, m22); Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; stan::test::expect_ad(tols, f, m23); + stan::test::expect_ad_matvar(f, m23); Eigen::MatrixXd m32(3, 2); m32 << 1, 3, -5, 7, 9, -11; stan::test::expect_ad(tols, f, m32); + stan::test::expect_ad_matvar(f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; stan::test::expect_ad(tols, f, a22); + stan::test::expect_ad_matvar(f, a22); } diff --git a/test/unit/math/mix/fun/svd_V_test.cpp b/test/unit/math/mix/fun/svd_V_test.cpp index eaaaa9f5550..efe4202c0b3 100644 --- a/test/unit/math/mix/fun/svd_V_test.cpp +++ b/test/unit/math/mix/fun/svd_V_test.cpp @@ -1,4 +1,5 @@ #include +#include TEST(MathMixMatFun, svd_V) { auto f = [](const auto& x) { return stan::math::svd_V(x); }; @@ -8,25 +9,30 @@ TEST(MathMixMatFun, svd_V) { tols.hessian_fvar_hessian_ = 1e-2; Eigen::MatrixXd m00(0, 0); - stan::test::expect_ad(f, m00); + EXPECT_THROW(f(m00), std::invalid_argument); Eigen::MatrixXd m11(1, 1); m11 << 1.1; stan::test::expect_ad(tols, f, m11); + stan::test::expect_ad_matvar(f, m11); Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; stan::test::expect_ad(tols, f, m22); + stan::test::expect_ad_matvar(f, m22); Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; stan::test::expect_ad(tols, f, m23); + stan::test::expect_ad_matvar(f, m23); Eigen::MatrixXd m32(3, 2); m32 << 1, 3, -5, 7, 9, -11; stan::test::expect_ad(tols, f, m32); + stan::test::expect_ad_matvar(f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; stan::test::expect_ad(tols, f, a22); + stan::test::expect_ad_matvar(f, a22); } diff --git a/test/unit/math/prim/fun/singular_values_test.cpp b/test/unit/math/prim/fun/singular_values_test.cpp index 5877b0c6279..571f3feec27 100644 --- a/test/unit/math/prim/fun/singular_values_test.cpp +++ b/test/unit/math/prim/fun/singular_values_test.cpp @@ -1,6 +1,7 @@ #include #include #include +#include TEST(MathMatrixPrimMat, singular_values) { using stan::math::matrix_d; @@ -8,7 +9,7 @@ TEST(MathMatrixPrimMat, singular_values) { using stan::math::vector_d; matrix_d m0(0, 0); - EXPECT_NO_THROW(singular_values(m0)); + EXPECT_THROW(singular_values(m0), std::invalid_argument); matrix_d m1(1, 1); m1 << 1.0; From 02c5a81342d46579dff511c741bdd8a7165d11b0 Mon Sep 17 00:00:00 2001 From: Ben Date: Sun, 10 Jan 2021 11:24:41 -0500 Subject: [PATCH 31/32] Removed extra comments, rearranged includes, and removed custom tolerances in tests (Issue #2101) --- stan/math/rev/fun/singular_values.hpp | 5 ++--- stan/math/rev/fun/svd_U.hpp | 17 +++++++++-------- stan/math/rev/fun/svd_V.hpp | 18 +++++++++--------- .../unit/math/mix/fun/singular_values_test.cpp | 14 +++++--------- test/unit/math/mix/fun/svd_U_test.cpp | 14 +++++--------- test/unit/math/mix/fun/svd_V_test.cpp | 14 +++++--------- 6 files changed, 35 insertions(+), 47 deletions(-) diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp index bd3fe7a03b7..ec259a08cd7 100644 --- a/stan/math/rev/fun/singular_values.hpp +++ b/stan/math/rev/fun/singular_values.hpp @@ -1,11 +1,10 @@ #ifndef STAN_MATH_REV_FUN_SINGULAR_VALUES_HPP #define STAN_MATH_REV_FUN_SINGULAR_VALUES_HPP -#include -#include #include #include -#include +#include +#include namespace stan { namespace math { diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index d984c57cfa4..a85128d9ee5 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -1,12 +1,12 @@ #ifndef STAN_MATH_REV_FUN_SVD_U_HPP #define STAN_MATH_REV_FUN_SVD_U_HPP -#include -#include -#include #include #include #include +#include +#include +#include namespace stan { namespace math { @@ -27,12 +27,12 @@ inline auto svd_U(const EigMat& m) { check_nonzero_size("svd_U", "m", m); const int M = std::min(m.rows(), m.cols()); - auto arena_m = to_arena(m); // N by P + auto arena_m = to_arena(m); Eigen::JacobiSVD svd( arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); - auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M + auto arena_D = to_arena(svd.singularValues()); arena_t arena_Fp(M, M); @@ -47,13 +47,13 @@ inline auto svd_U(const EigMat& m) { } } - arena_t arena_U = svd.matrixU(); // N by M - auto arena_V = to_arena(svd.matrixV()); // P by M + arena_t arena_U = svd.matrixU(); + auto arena_V = to_arena(svd.matrixV()); reverse_pass_callback( [arena_m, arena_U, arena_D, arena_V, arena_Fp, M]() mutable { Eigen::MatrixXd UUadjT - = arena_U.val_op().transpose() * arena_U.adj_op(); // M by M + = arena_U.val_op().transpose() * arena_U.adj_op(); arena_m.adj() += .5 * arena_U.val_op() * (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()) @@ -67,6 +67,7 @@ inline auto svd_U(const EigMat& m) { return ret_type(arena_U); } + } // namespace math } // namespace stan diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index 3ec9800337a..ff41631a541 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -1,13 +1,12 @@ #ifndef STAN_MATH_REV_FUN_SVD_V_HPP #define STAN_MATH_REV_FUN_SVD_V_HPP -#include -#include -#include #include #include #include -#include +#include +#include +#include namespace stan { namespace math { @@ -28,12 +27,12 @@ inline auto svd_V(const EigMat& m) { check_nonzero_size("svd_V", "m", m); const int M = std::min(m.rows(), m.cols()); - auto arena_m = to_arena(m); // N by P + auto arena_m = to_arena(m); Eigen::JacobiSVD svd( arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV); - auto arena_D = to_arena(svd.singularValues()); // size min(N, P) = M + auto arena_D = to_arena(svd.singularValues()); arena_t arena_Fm(M, M); @@ -48,13 +47,13 @@ inline auto svd_V(const EigMat& m) { } } - auto arena_U = to_arena(svd.matrixU()); // N by M - arena_t arena_V = svd.matrixV(); // P by M + auto arena_U = to_arena(svd.matrixU()); + arena_t arena_V = svd.matrixV(); reverse_pass_callback([arena_m, arena_U, arena_D, arena_V, arena_Fm, M]() mutable { Eigen::MatrixXd VTVadj - = arena_V.val_op().transpose() * arena_V.adj_op(); // M by M + = arena_V.val_op().transpose() * arena_V.adj_op(); arena_m.adj() += 0.5 * arena_U * (arena_Fm.array() * (VTVadj - VTVadj.transpose()).array()) @@ -68,6 +67,7 @@ inline auto svd_V(const EigMat& m) { return ret_type(arena_V); } + } // namespace math } // namespace stan diff --git a/test/unit/math/mix/fun/singular_values_test.cpp b/test/unit/math/mix/fun/singular_values_test.cpp index 76c88205f71..f1e5ec0b98b 100644 --- a/test/unit/math/mix/fun/singular_values_test.cpp +++ b/test/unit/math/mix/fun/singular_values_test.cpp @@ -4,33 +4,29 @@ TEST(MathMixMatFun, singularValues) { auto f = [](const auto& x) { return stan::math::singular_values(x); }; - stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; - Eigen::MatrixXd m00(0, 0); EXPECT_THROW(f(m00), std::invalid_argument); Eigen::MatrixXd m11(1, 1); m11 << 1.1; - stan::test::expect_ad(tols, f, m11); + stan::test::expect_ad(f, m11); stan::test::expect_ad_matvar(f, m11); Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; - stan::test::expect_ad(tols, f, m22); + stan::test::expect_ad(f, m22); stan::test::expect_ad_matvar(f, m22); Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; Eigen::MatrixXd m32 = m23.transpose(); - stan::test::expect_ad(tols, f, m23); - stan::test::expect_ad(tols, f, m32); + stan::test::expect_ad(f, m23); + stan::test::expect_ad(f, m32); stan::test::expect_ad_matvar(f, m23); stan::test::expect_ad_matvar(f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; - stan::test::expect_ad(tols, f, a22); + stan::test::expect_ad(f, a22); stan::test::expect_ad_matvar(f, a22); } diff --git a/test/unit/math/mix/fun/svd_U_test.cpp b/test/unit/math/mix/fun/svd_U_test.cpp index fd3a8389861..e2c10b1b89a 100644 --- a/test/unit/math/mix/fun/svd_U_test.cpp +++ b/test/unit/math/mix/fun/svd_U_test.cpp @@ -4,35 +4,31 @@ TEST(MathMixMatFun, svd_U) { auto f = [](const auto& x) { return stan::math::svd_U(x); }; - stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; - Eigen::MatrixXd m00(0, 0); EXPECT_THROW(f(m00), std::invalid_argument); Eigen::MatrixXd m11(1, 1); m11 << 1.1; - stan::test::expect_ad(tols, f, m11); + stan::test::expect_ad(f, m11); stan::test::expect_ad_matvar(f, m11); Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; - stan::test::expect_ad(tols, f, m22); + stan::test::expect_ad(f, m22); stan::test::expect_ad_matvar(f, m22); Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; - stan::test::expect_ad(tols, f, m23); + stan::test::expect_ad(f, m23); stan::test::expect_ad_matvar(f, m23); Eigen::MatrixXd m32(3, 2); m32 << 1, 3, -5, 7, 9, -11; - stan::test::expect_ad(tols, f, m32); + stan::test::expect_ad(f, m32); stan::test::expect_ad_matvar(f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; - stan::test::expect_ad(tols, f, a22); + stan::test::expect_ad(f, a22); stan::test::expect_ad_matvar(f, a22); } diff --git a/test/unit/math/mix/fun/svd_V_test.cpp b/test/unit/math/mix/fun/svd_V_test.cpp index efe4202c0b3..18eaa7c44b6 100644 --- a/test/unit/math/mix/fun/svd_V_test.cpp +++ b/test/unit/math/mix/fun/svd_V_test.cpp @@ -4,35 +4,31 @@ TEST(MathMixMatFun, svd_V) { auto f = [](const auto& x) { return stan::math::svd_V(x); }; - stan::test::ad_tolerances tols; - tols.hessian_hessian_ = 1e-2; - tols.hessian_fvar_hessian_ = 1e-2; - Eigen::MatrixXd m00(0, 0); EXPECT_THROW(f(m00), std::invalid_argument); Eigen::MatrixXd m11(1, 1); m11 << 1.1; - stan::test::expect_ad(tols, f, m11); + stan::test::expect_ad(f, m11); stan::test::expect_ad_matvar(f, m11); Eigen::MatrixXd m22(2, 2); m22 << 3, -5, 7, 11; - stan::test::expect_ad(tols, f, m22); + stan::test::expect_ad(f, m22); stan::test::expect_ad_matvar(f, m22); Eigen::MatrixXd m23(2, 3); m23 << 3, 5, -7, -11, 13, -17; - stan::test::expect_ad(tols, f, m23); + stan::test::expect_ad(f, m23); stan::test::expect_ad_matvar(f, m23); Eigen::MatrixXd m32(3, 2); m32 << 1, 3, -5, 7, 9, -11; - stan::test::expect_ad(tols, f, m32); + stan::test::expect_ad(f, m32); stan::test::expect_ad_matvar(f, m32); Eigen::MatrixXd a22(2, 2); a22 << 1, 2, 3, 4; - stan::test::expect_ad(tols, f, a22); + stan::test::expect_ad(f, a22); stan::test::expect_ad_matvar(f, a22); } From 784284fe3a931be7105834f90fe021c297949ec1 Mon Sep 17 00:00:00 2001 From: Stan Jenkins Date: Sun, 10 Jan 2021 16:26:01 +0000 Subject: [PATCH 32/32] [Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.04.1 (tags/RELEASE_600/final) --- stan/math/rev/fun/svd_U.hpp | 27 +++++++++++++-------------- stan/math/rev/fun/svd_V.hpp | 3 +-- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp index a85128d9ee5..589c9c069a4 100644 --- a/stan/math/rev/fun/svd_U.hpp +++ b/stan/math/rev/fun/svd_U.hpp @@ -50,20 +50,19 @@ inline auto svd_U(const EigMat& m) { arena_t arena_U = svd.matrixU(); auto arena_V = to_arena(svd.matrixV()); - reverse_pass_callback( - [arena_m, arena_U, arena_D, arena_V, arena_Fp, M]() mutable { - Eigen::MatrixXd UUadjT - = arena_U.val_op().transpose() * arena_U.adj_op(); - arena_m.adj() - += .5 * arena_U.val_op() - * (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()) - .matrix() - * arena_V.transpose() - + (Eigen::MatrixXd::Identity(arena_m.rows(), arena_m.rows()) - - arena_U.val_op() * arena_U.val_op().transpose()) - * arena_U.adj_op() * arena_D.asDiagonal().inverse() - * arena_V.transpose(); - }); + reverse_pass_callback([arena_m, arena_U, arena_D, arena_V, arena_Fp, + M]() mutable { + Eigen::MatrixXd UUadjT = arena_U.val_op().transpose() * arena_U.adj_op(); + arena_m.adj() + += .5 * arena_U.val_op() + * (arena_Fp.array() * (UUadjT - UUadjT.transpose()).array()) + .matrix() + * arena_V.transpose() + + (Eigen::MatrixXd::Identity(arena_m.rows(), arena_m.rows()) + - arena_U.val_op() * arena_U.val_op().transpose()) + * arena_U.adj_op() * arena_D.asDiagonal().inverse() + * arena_V.transpose(); + }); return ret_type(arena_U); } diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp index ff41631a541..cbd267b0f53 100644 --- a/stan/math/rev/fun/svd_V.hpp +++ b/stan/math/rev/fun/svd_V.hpp @@ -52,8 +52,7 @@ inline auto svd_V(const EigMat& m) { reverse_pass_callback([arena_m, arena_U, arena_D, arena_V, arena_Fm, M]() mutable { - Eigen::MatrixXd VTVadj - = arena_V.val_op().transpose() * arena_V.adj_op(); + Eigen::MatrixXd VTVadj = arena_V.val_op().transpose() * arena_V.adj_op(); arena_m.adj() += 0.5 * arena_U * (arena_Fm.array() * (VTVadj - VTVadj.transpose()).array())