Add singular_values, svd_U, and svd_V

stan/math/prim/fun.hpp

@@ -312,6 +312,8 @@
 #include <stan/math/prim/fun/sub_row.hpp>
 #include <stan/math/prim/fun/subtract.hpp>
 #include <stan/math/prim/fun/sum.hpp>
+#include <stan/math/prim/fun/svd_U.hpp>
+#include <stan/math/prim/fun/svd_V.hpp>
 #include <stan/math/prim/fun/symmetrize_from_lower_tri.hpp>
 #include <stan/math/prim/fun/tail.hpp>
 #include <stan/math/prim/fun/tan.hpp>

stan/math/prim/fun/Eigen.hpp

@@ -23,6 +23,7 @@
 #include <Eigen/Sparse>
 #include <Eigen/QR>
 #include <Eigen/src/Core/NumTraits.h>
+#include <Eigen/SVD>
 
     namespace Eigen {
 

stan/math/prim/fun/singular_values.hpp

@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/err/check_nonzero_size.hpp>
 
 namespace stan {
 namespace math {
@@ -17,12 +18,11 @@ namespace math {
  * @param m Specified matrix.
  * @return Singular values of the matrix.
  */
-template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr>
+template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
+          require_not_st_var<EigMat>* = nullptr>
 Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, 1> singular_values(
     const EigMat& m) {
-  if (m.size() == 0) {
-    return {};
-  }
+  check_nonzero_size("singular_values", "m", m);
 
   return Eigen::JacobiSVD<Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic,
                                         Eigen::Dynamic> >(m)

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 <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Given input matrix m, return matrix U where `m = UDV^{T}`
+ *
+ * @tparam EigMat type of the matrix
+ * @param m MxN input matrix
+ * @return Orthogonal matrix U
+ */
+template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
+          require_not_st_var<EigMat>* = nullptr>
+Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic> svd_U(
+    const EigMat& m) {
+  check_nonzero_size("svd_U", "m", m);
+
+  return Eigen::JacobiSVD<Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic,
+                                        Eigen::Dynamic> >(m,
+                                                          Eigen::ComputeThinU)
+      .matrixU();
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

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 <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Given input matrix m, return matrix V where `m = UDV^{T}`
+ *
+ * @tparam EigMat type of the matrix
+ * @param m MxN input matrix
+ * @return Orthogonal matrix V
+ */
+template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
+          require_not_st_var<EigMat>* = nullptr>
+Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic> svd_V(
+    const EigMat& m) {
+  check_nonzero_size("svd_V", "m", m);
+
+  return Eigen::JacobiSVD<Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic,
+                                        Eigen::Dynamic> >(m,
+                                                          Eigen::ComputeThinV)
+      .matrixV();
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/rev/fun.hpp

@@ -144,6 +144,9 @@
 #include <stan/math/rev/fun/sd.hpp>
 #include <stan/math/rev/fun/simplex_constrain.hpp>
 #include <stan/math/rev/fun/sin.hpp>
+#include <stan/math/rev/fun/singular_values.hpp>
+#include <stan/math/rev/fun/svd_U.hpp>
+#include <stan/math/rev/fun/svd_V.hpp>
 #include <stan/math/rev/fun/sinh.hpp>
 #include <stan/math/rev/fun/softmax.hpp>
 #include <stan/math/rev/fun/sqrt.hpp>

stan/math/rev/fun/singular_values.hpp

@@ -0,0 +1,48 @@
+#ifndef STAN_MATH_REV_FUN_SINGULAR_VALUES_HPP
+#define STAN_MATH_REV_FUN_SINGULAR_VALUES_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/err/check_nonzero_size.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return the singular values of the specified matrix.
+ *
+ * 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
+ * @return Singular values of matrix
+ */
+template <typename EigMat, require_rev_matrix_t<EigMat>* = nullptr>
+inline auto singular_values(const EigMat& m) {
+  using ret_type = return_var_matrix_t<Eigen::VectorXd, EigMat>;
+  check_nonzero_size("singular_values", "m", m);
+
+  auto arena_m = to_arena(m);
+
+  Eigen::JacobiSVD<Eigen::MatrixXd> svd(
+      arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV);
+
+  arena_t<ret_type> singular_values = svd.singularValues();
+
+  auto arena_U = to_arena(svd.matrixU());
+  auto arena_V = to_arena(svd.matrixV());
+
+  reverse_pass_callback([arena_m, arena_U, singular_values, arena_V]() mutable {
+    arena_m.adj()
+        += arena_U * singular_values.adj().asDiagonal() * arena_V.transpose();
+  });
+
+  return ret_type(singular_values);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/rev/fun/svd_U.hpp

@@ -0,0 +1,73 @@
+#ifndef STAN_MATH_REV_FUN_SVD_U_HPP
+#define STAN_MATH_REV_FUN_SVD_U_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/err/check_nonzero_size.hpp>
+#include <stan/math/prim/fun/typedefs.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/value_of.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Given input matrix m, return matrix U where `m = UDV^{T}`
+ *
+ * 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
+ * @return Orthogonal matrix U
+ */
+template <typename EigMat, require_rev_matrix_t<EigMat>* = nullptr>
+inline auto svd_U(const EigMat& m) {
+  using ret_type = return_var_matrix_t<Eigen::MatrixXd, EigMat>;
+  check_nonzero_size("svd_U", "m", m);
+
+  const int M = std::min(m.rows(), m.cols());
+  auto arena_m = to_arena(m);
+
+  Eigen::JacobiSVD<Eigen::MatrixXd> svd(
+      arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV);
+
+  auto arena_D = to_arena(svd.singularValues());
+
+  arena_t<Eigen::MatrixXd> arena_Fp(M, M);
+
+  for (int i = 0; i < M; i++) {
+    for (int j = 0; j < M; j++) {
+      if (j == i) {
+        arena_Fp(i, j) = 0.0;
+      } else {
+        arena_Fp(i, j)
+            = 1.0 / (arena_D[j] - arena_D[i]) + 1.0 / (arena_D[i] + arena_D[j]);
+      }
+    }
+  }
+
+  arena_t<ret_type> 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();
+  });
+
+  return ret_type(arena_U);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/rev/fun/svd_V.hpp

@@ -0,0 +1,73 @@
+#ifndef STAN_MATH_REV_FUN_SVD_V_HPP
+#define STAN_MATH_REV_FUN_SVD_V_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/err/check_nonzero_size.hpp>
+#include <stan/math/prim/fun/typedefs.hpp>
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/fun/value_of.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Given input matrix m, return matrix V where `m = UDV^{T}`
+ *
+ * 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
+ * @return Orthogonal matrix V
+ */
+template <typename EigMat, require_rev_matrix_t<EigMat>* = nullptr>
+inline auto svd_V(const EigMat& m) {
+  using ret_type = return_var_matrix_t<Eigen::MatrixXd, EigMat>;
+  check_nonzero_size("svd_V", "m", m);
+
+  const int M = std::min(m.rows(), m.cols());
+  auto arena_m = to_arena(m);
+
+  Eigen::JacobiSVD<Eigen::MatrixXd> svd(
+      arena_m.val(), Eigen::ComputeThinU | Eigen::ComputeThinV);
+
+  auto arena_D = to_arena(svd.singularValues());
+
+  arena_t<Eigen::MatrixXd> arena_Fm(M, M);
+
+  for (int i = 0; i < M; i++) {
+    for (int j = 0; j < M; j++) {
+      if (j == i) {
+        arena_Fm(i, j) = 0.0;
+      } else {
+        arena_Fm(i, j)
+            = 1.0 / (arena_D[j] - arena_D[i]) - 1.0 / (arena_D[i] + arena_D[j]);
+      }
+    }
+  }
+
+  auto arena_U = to_arena(svd.matrixU());
+  arena_t<ret_type> 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
+               * (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);
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

test/unit/math/mix/fun/singular_values_test.cpp

@@ -1,30 +1,32 @@
 #include <test/unit/math/test_ad.hpp>
+#include <stdexcept>
 
 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);
-  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(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);
 }

test/unit/math/mix/fun/svd_U_test.cpp

@@ -0,0 +1,34 @@
+#include <test/unit/math/test_ad.hpp>
+#include <stdexcept>
+
+TEST(MathMixMatFun, svd_U) {
+  auto f = [](const auto& x) { return stan::math::svd_U(x); };
+
+  Eigen::MatrixXd m00(0, 0);
+  EXPECT_THROW(f(m00), std::invalid_argument);
+
+  Eigen::MatrixXd m11(1, 1);
+  m11 << 1.1;
+  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(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(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(f, m32);
+  stan::test::expect_ad_matvar(f, m32);
+
+  Eigen::MatrixXd a22(2, 2);
+  a22 << 1, 2, 3, 4;
+  stan::test::expect_ad(f, a22);
+  stan::test::expect_ad_matvar(f, a22);
+}