diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp
index 3405f36a221..d7be8fc9f6c 100644
--- a/stan/math/prim/fun.hpp
+++ b/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>
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 <Eigen/Sparse>
 #include <Eigen/QR>
 #include <Eigen/src/Core/NumTraits.h>
+#include <Eigen/SVD>
 
     namespace Eigen {
 
diff --git a/stan/math/prim/fun/singular_values.hpp b/stan/math/prim/fun/singular_values.hpp
index cbb89938c49..c177c8867d9 100644
--- a/stan/math/prim/fun/singular_values.hpp
+++ b/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)
diff --git a/stan/math/prim/fun/svd_U.hpp b/stan/math/prim/fun/svd_U.hpp
new file mode 100644
index 00000000000..b69d7d9566d
--- /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 <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
diff --git a/stan/math/prim/fun/svd_V.hpp b/stan/math/prim/fun/svd_V.hpp
new file mode 100644
index 00000000000..ef85274acd7
--- /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 <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
diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp
index d76cafc64ed..0d7491522b4 100644
--- a/stan/math/rev/fun.hpp
+++ b/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>
diff --git a/stan/math/rev/fun/singular_values.hpp b/stan/math/rev/fun/singular_values.hpp
new file mode 100644
index 00000000000..ec259a08cd7
--- /dev/null
+++ b/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
diff --git a/stan/math/rev/fun/svd_U.hpp b/stan/math/rev/fun/svd_U.hpp
new file mode 100644
index 00000000000..589c9c069a4
--- /dev/null
+++ b/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
diff --git a/stan/math/rev/fun/svd_V.hpp b/stan/math/rev/fun/svd_V.hpp
new file mode 100644
index 00000000000..cbd267b0f53
--- /dev/null
+++ b/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
diff --git a/test/unit/math/mix/fun/singular_values_test.cpp b/test/unit/math/mix/fun/singular_values_test.cpp
index d6b757035ef..f1e5ec0b98b 100644
--- a/test/unit/math/mix/fun/singular_values_test.cpp
+++ b/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);
 }
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..e2c10b1b89a
--- /dev/null
+++ b/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);
+}
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..18eaa7c44b6
--- /dev/null
+++ b/test/unit/math/mix/fun/svd_V_test.cpp
@@ -0,0 +1,34 @@
+#include <test/unit/math/test_ad.hpp>
+#include <stdexcept>
+
+TEST(MathMixMatFun, svd_V) {
+  auto f = [](const auto& x) { return stan::math::svd_V(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);
+}
diff --git a/test/unit/math/prim/fun/singular_values_test.cpp b/test/unit/math/prim/fun/singular_values_test.cpp
index 229a174cb03..571f3feec27 100644
--- a/test/unit/math/prim/fun/singular_values_test.cpp
+++ b/test/unit/math/prim/fun/singular_values_test.cpp
@@ -1,13 +1,35 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
+#include <test/unit/util.hpp>
+#include <stdexcept>
 
 TEST(MathMatrixPrimMat, singular_values) {
+  using stan::math::matrix_d;
   using stan::math::singular_values;
+  using stan::math::vector_d;
 
-  stan::math::matrix_d m0(0, 0);
-  EXPECT_NO_THROW(singular_values(m0));
+  matrix_d m0(0, 0);
+  EXPECT_THROW(singular_values(m0), std::invalid_argument);
 
-  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/prim/fun/svd_U_test.cpp b/test/unit/math/prim/fun/svd_U_test.cpp
new file mode 100644
index 00000000000..a0bec909f81
--- /dev/null
+++ b/test/unit/math/prim/fun/svd_U_test.cpp
@@ -0,0 +1,42 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <test/unit/util.hpp>
+#include <stdexcept>
+
+TEST(MathMatrixPrimMat, svd_U) {
+  using stan::math::matrix_d;
+  using stan::math::svd_U;
+
+  // 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.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.
+}
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..42dc5e85be4
--- /dev/null
+++ b/test/unit/math/prim/fun/svd_V_test.cpp
@@ -0,0 +1,42 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <test/unit/util.hpp>
+#include <stdexcept>
+
+TEST(MathMatrixPrimMat, svd_V) {
+  using stan::math::matrix_d;
+  using stan::math::svd_V;
+
+  // 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, 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));  // R's SVD returns different signs than Eigen.
+}
diff --git a/test/unit/math/rev/fun/singular_values_test.cpp b/test/unit/math/rev/fun/singular_values_test.cpp
new file mode 100644
index 00000000000..ea53205e884
--- /dev/null
+++ b/test/unit/math/rev/fun/singular_values_test.cpp
@@ -0,0 +1,24 @@
+#include <stan/math/prim.hpp>
+#include <stan/math/rev/fun/singular_values.hpp>
+#include <stan/math/rev/fun/sum.hpp>
+#include <stan/math/rev/fun/log.hpp>
+#include <gtest/gtest.h>
+
+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 = stan::math::singular_values(a_v);
+  auto logdet = stan::math::sum(stan::math::log(w));
+
+  stan::math::set_zero_all_adjoints();
+  logdet.grad();
+  ASSERT_TRUE(a_inv.val().isApprox(a_v.adj()));
+}