diff --git a/stan/math/fwd/fun/log_softmax.hpp b/stan/math/fwd/fun/log_softmax.hpp
index a9685d7ca00..8f03dfd5330 100644
--- a/stan/math/fwd/fun/log_softmax.hpp
+++ b/stan/math/fwd/fun/log_softmax.hpp
@@ -14,10 +14,11 @@ namespace stan {
 namespace math {
 
 /**
- * Return the log softmax of each vector in a container of `fvar` values.
+ * Return the log softmax of each vector or matrix in a container of `fvar`
+ * values.
  *
  * @tparam T `std::vector` whose scalar type is `fvar`
- * @param x container of vectors to transform
+ * @param x container of vectors or matrices to transform
  * @return container of log softmax results
  */
 template <typename T, require_std_vector_st<is_fvar, T>* = nullptr>
@@ -28,27 +29,31 @@ inline auto log_softmax(T&& x) {
 }
 
 /**
- * Return the log softmax of the specified vector of `fvar` values.
+ * Return the log softmax of the specified vector or matrix of `fvar` values.
  *
- * @tparam Vec Eigen vector with `fvar` scalar
- * @param x vector to transform
- * @return log softmax of the vector, or an empty result if the input is empty
+ * @tparam Mat Eigen vector or matrix with `fvar` scalar
+ * @param x vector or matrix to transform
+ * @return log softmax of the vector or matrix, or an empty result if the
+ * input is empty
  */
-template <typename Vec, require_eigen_vector_vt<is_fvar, Vec>* = nullptr>
-inline auto log_softmax(Vec&& x) {
-  using vec = std::decay_t<Vec>;
-  constexpr int Rows = vec::RowsAtCompileTime;
-  constexpr int Cols = vec::ColsAtCompileTime;
-  using T = typename value_type_t<vec>::Scalar;
-  decltype(auto) x_ref = to_ref(std::forward<Vec>(x));
+template <typename Mat, require_eigen_vt<is_fvar, Mat>* = nullptr>
+inline auto log_softmax(Mat&& x) {
+  using mat = std::decay_t<Mat>;
+  constexpr int Rows = mat::RowsAtCompileTime;
+  constexpr int Cols = mat::ColsAtCompileTime;
+  using T = typename value_type_t<mat>::Scalar;
+  decltype(auto) x_ref = to_ref(std::forward<Mat>(x));
   if (x_ref.size() == 0) {
     return Eigen::Matrix<fvar<T>, Rows, Cols>{};
   }
-  const auto s = softmax(value_of(x_ref));
+  const auto x_val = value_of(x_ref);
+  const auto lse = log_sum_exp(x_val);
+  const auto s = softmax(x_val);
   const auto d_in = x_ref.d();
-  const auto dot_sd = s.dot(d_in);
-  Eigen::Matrix<fvar<T>, Rows, Cols> result(x_ref.size());
-  result.val() = s.array().log().matrix();
+  const auto dot_sd = (s.array() * d_in.array()).sum();
+
+  Eigen::Matrix<fvar<T>, Rows, Cols> result(x_ref.rows(), x_ref.cols());
+  result.val() = (x_val.array() - lse).matrix();
   result.d() = (d_in.array() - dot_sd).matrix();
   return result;
 }
diff --git a/stan/math/fwd/fun/log_sum_exp.hpp b/stan/math/fwd/fun/log_sum_exp.hpp
index 926838d143c..812be56d9ee 100644
--- a/stan/math/fwd/fun/log_sum_exp.hpp
+++ b/stan/math/fwd/fun/log_sum_exp.hpp
@@ -8,6 +8,7 @@
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
+#include <stan/math/prim/fun/softmax.hpp>
 #include <cmath>
 #include <vector>
 
@@ -56,11 +57,10 @@ inline auto log_sum_exp(T&& x) {
         using T_fvar_inner = typename value_type_t<decltype(v)>::Scalar;
         using mat_type = Eigen::Matrix<T_fvar_inner, -1, -1>;
         mat_type vals = v.val();
-        mat_type exp_vals = vals.array().exp();
 
-        return fvar<T_fvar_inner>(
-            log_sum_exp(vals),
-            v.d().cwiseProduct(exp_vals).sum() / exp_vals.sum());
+        const auto probs = softmax(vals);
+        return fvar<T_fvar_inner>(log_sum_exp(vals),
+                                  v.d().cwiseProduct(probs).sum());
       });
 }
 
diff --git a/stan/math/fwd/fun/softmax.hpp b/stan/math/fwd/fun/softmax.hpp
index 97d019b6c23..a74b25762e7 100644
--- a/stan/math/fwd/fun/softmax.hpp
+++ b/stan/math/fwd/fun/softmax.hpp
@@ -29,24 +29,25 @@ inline auto softmax(T&& x) {
 /**
  * Return the softmax of the specified vector of `fvar` values.
  *
- * @tparam Vec Eigen vector with `fvar` scalar
- * @param x vector to transform
- * @return softmax of the vector, or an empty result if the input is empty
+ * @tparam Mat Eigen vector or matrix with `fvar` scalar
+ * @param x vector or matrix to transform
+ * @return softmax of the vector, matrix, or an empty result if the input is
+ *   empty
  */
-template <typename Vec, require_eigen_vector_vt<is_fvar, Vec>* = nullptr>
-inline auto softmax(Vec&& x) {
-  using vec = std::decay_t<Vec>;
-  constexpr int Rows = vec::RowsAtCompileTime;
-  constexpr int Cols = vec::ColsAtCompileTime;
-  using T = typename value_type_t<vec>::Scalar;
-  decltype(auto) x_ref = to_ref(std::forward<Vec>(x));
+template <typename Mat, require_eigen_vt<is_fvar, Mat>* = nullptr>
+inline auto softmax(Mat&& x) {
+  using mat = std::decay_t<Mat>;
+  constexpr int Rows = mat::RowsAtCompileTime;
+  constexpr int Cols = mat::ColsAtCompileTime;
+  using T = typename value_type_t<mat>::Scalar;
+  decltype(auto) x_ref = to_ref(std::forward<Mat>(x));
   if (x_ref.size() == 0) {
     return Eigen::Matrix<fvar<T>, Rows, Cols>{};
   }
   const auto s = softmax(value_of(x_ref));
   const auto d_in = x_ref.d();
-  const auto dot_sd = s.dot(d_in);
-  Eigen::Matrix<fvar<T>, Rows, Cols> result(x_ref.size());
+  const auto dot_sd = (s.array() * d_in.array()).sum();
+  Eigen::Matrix<fvar<T>, Rows, Cols> result(x_ref.rows(), x_ref.cols());
   result.val() = s;
   result.d() = (s.array() * (d_in.array() - dot_sd)).matrix();
   return result;
diff --git a/stan/math/prim/fun/log_softmax.hpp b/stan/math/prim/fun/log_softmax.hpp
index d6c9017edb6..f1d1cd0fef5 100644
--- a/stan/math/prim/fun/log_softmax.hpp
+++ b/stan/math/prim/fun/log_softmax.hpp
@@ -12,9 +12,10 @@ namespace stan {
 namespace math {
 
 /**
- * Return the natural logarithm of the softmax of the specified
- * vector, or of each vector in a container.
- *
+ * Return the natural logarithm of the softmax of the specified vector or
+ * matrix, or of each vector or matrix in a container. For a matrix, the
+ * log-softmax is taken over all elements.
+ *  *
  * \f$
  * \log \mbox{softmax}(y)
  * \ = \ y - \log \sum_{k=1}^K \exp(y_k)
@@ -35,17 +36,15 @@ namespace math {
  * \right.
  * \f$
  *
- * @tparam Container type of input: an Eigen vector, `std::vector` of doubles,
- *   or nested container whose scalar type is arithmetic
- * @param x vector or container of vectors to transform
- * @return log softmax of the input, preserving the container structure; an
- *   empty result if any input vector is empty
+ * @tparam Container type of input: an Eigen vector, Eigen matrix,
+ *   `std::vector` of vectors or matrices, or nested container whose scalar
+ *   type is arithmetic
+ * @param x vector, matrix, or container to transform.
+ * @return softmax of the input, preserving the container structure; an empty
+ *   result if any input vector or matrix is empty.
  */
 template <typename Container, require_st_arithmetic<Container>* = nullptr,
-          require_container_t<Container>* = nullptr,
-          require_not_t<bool_constant<
-              is_eigen<std::decay_t<Container>>::value
-              && !is_eigen_vector<std::decay_t<Container>>::value>>* = nullptr>
+          require_container_t<Container>* = nullptr>
 inline auto log_softmax(Container&& x) {
   return make_holder(
       [](auto&& a) {
diff --git a/stan/math/prim/fun/softmax.hpp b/stan/math/prim/fun/softmax.hpp
index 63a51876a60..112abac8612 100644
--- a/stan/math/prim/fun/softmax.hpp
+++ b/stan/math/prim/fun/softmax.hpp
@@ -11,7 +11,9 @@ namespace stan {
 namespace math {
 
 /**
- * Return the softmax of the specified vector, or of each vector in a container.
+ * Return the softmax of the specified vector or matrix, or of each
+ * vector or matrix in a container. For a matrix, the softmax is
+ * taken over all elements.
  *
  * \f$
  * \mbox{softmax}(y)
@@ -38,17 +40,15 @@ namespace math {
  * \end{array}
  * \f$
  *
- * @tparam Container type of input: an Eigen vector, `std::vector` of doubles,
- *   or nested container whose scalar type is arithmetic
- * @param x vector or container of vectors to transform
+ * @tparam Container type of input: an Eigen vector, Eigen matrix,
+ *   `std::vector` of vectors or matrices, or nested container whose scalar
+ *   type is arithmetic
+ * @param x vector, matrix, or container to transform.
  * @return softmax of the input, preserving the container structure; an empty
- *   result if any input vector is empty
+ *   result if any input vector or matrix is empty.
  */
 template <typename Container, require_st_arithmetic<Container>* = nullptr,
-          require_container_t<Container>* = nullptr,
-          require_not_t<bool_constant<
-              is_eigen<std::decay_t<Container>>::value
-              && !is_eigen_vector<std::decay_t<Container>>::value>>* = nullptr>
+          require_container_t<Container>* = nullptr>
 inline auto softmax(Container&& x) {
   return make_holder(
       [](auto&& a) {
diff --git a/stan/math/rev/fun/log_softmax.hpp b/stan/math/rev/fun/log_softmax.hpp
index 4c77c09e15f..845d03aa147 100644
--- a/stan/math/rev/fun/log_softmax.hpp
+++ b/stan/math/rev/fun/log_softmax.hpp
@@ -8,15 +8,17 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/fun/log_softmax.hpp>
+#include <stan/math/prim/fun/softmax.hpp>
 #include <stan/math/prim/functor/apply_vector_unary.hpp>
 
 namespace stan {
 namespace math {
 
 /**
- * Return the log softmax of the specified vector or row vector.
+ * Return the log softmax of the specified vector, row vector, or matrix.
  *
- * @tparam T a `var_value` or Eigen vector/row_vector with `var` scalar
+ * @tparam T a `var_value` or Eigen vector, row_vector, or matrix with
+ *   `var` scalar
  * @param x input
  * @return log softmax of the input, or an empty result if the input is empty
  */
@@ -31,8 +33,8 @@ inline auto log_softmax(T&& x) {
   arena_t<return_t> res = log_softmax(x_arena.val());
   reverse_pass_callback([x_arena, res]() mutable {
     const auto& res_adj = to_ref(res.adj());
-    x_arena.adj().array()
-        += res_adj.array() - res_adj.sum() * res.val().array().exp();
+    const auto s = softmax(x_arena.val());
+    x_arena.adj().array() += res_adj.array() - res_adj.sum() * s.array();
   });
   return res;
 }
diff --git a/stan/math/rev/fun/log_sum_exp.hpp b/stan/math/rev/fun/log_sum_exp.hpp
index ca6a1b483fb..be9cc20db9e 100644
--- a/stan/math/rev/fun/log_sum_exp.hpp
+++ b/stan/math/rev/fun/log_sum_exp.hpp
@@ -9,6 +9,7 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/constants.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
+#include <stan/math/prim/fun/softmax.hpp>
 #include <cmath>
 #include <vector>
 
@@ -71,8 +72,7 @@ inline var log_sum_exp(T&& v) {
   auto arena_v_val = to_arena(arena_v.val());
   var res = log_sum_exp(arena_v_val);
   reverse_pass_callback([arena_v, arena_v_val, res]() mutable {
-    arena_v.adj()
-        += res.adj() * (arena_v_val.array().val() - res.val()).exp().matrix();
+    arena_v.adj() += res.adj() * softmax(arena_v_val);
   });
 
   return res;
@@ -87,7 +87,7 @@ inline var log_sum_exp(T&& v) {
 template <typename T, require_var_matrix_t<T>* = nullptr>
 inline var log_sum_exp(const T& x) {
   return make_callback_vari(log_sum_exp(x.val()), [x](const auto& res) mutable {
-    x.adj() += res.adj() * (x.val().array().val() - res.val()).exp().matrix();
+    x.adj() += res.adj() * softmax(x.val());
   });
 }
 
diff --git a/stan/math/rev/fun/softmax.hpp b/stan/math/rev/fun/softmax.hpp
index cd377c9cb1d..d3d46665774 100644
--- a/stan/math/rev/fun/softmax.hpp
+++ b/stan/math/rev/fun/softmax.hpp
@@ -15,9 +15,10 @@ namespace stan {
 namespace math {
 
 /**
- * Return the softmax of the specified vector or row vector.
+ * Return the softmax of the specified vector, row vector, or matrix.
  *
- * @tparam T a `var_value` or Eigen vector/row_vector with `var` scalar
+ * @tparam T a `var_value` or Eigen vector, row_vector, or matrix with
+ *   `var` scalar
  * @param x input
  * @return softmax of the input, or an empty result if the input is empty
  */
@@ -31,17 +32,17 @@ inline auto softmax(T&& x) {
       = return_var_matrix_t<plain_type_t<decltype(x_arena.val())>, T>;
   arena_t<return_t> res = softmax(x_arena.val());
   reverse_pass_callback([x_arena, res]() mutable {
-    x_arena.adj().array()
-        += res.val().array() * (res.adj().array() - res.val().dot(res.adj()));
+    const auto dot_sd = (res.val().array() * res.adj().array()).sum();
+    x_arena.adj().array() += res.val().array() * (res.adj().array() - dot_sd);
   });
   return res;
 }
 
 /**
- * Return the softmax of each vector in an array.
+ * Return the softmax of each vector or matrix in an array.
  *
  * @tparam T `std::vector` whose scalar type is `var`
- * @param x array of vectors to transform
+ * @param x array of vectors or matrices to transform
  * @return array of softmax results
  */
 template <typename T, require_std_vector_st<is_var, T>* = nullptr>
diff --git a/test/unit/math/fwd/fun/log_softmax_test.cpp b/test/unit/math/fwd/fun/log_softmax_test.cpp
new file mode 100644
index 00000000000..9055b896af2
--- /dev/null
+++ b/test/unit/math/fwd/fun/log_softmax_test.cpp
@@ -0,0 +1,112 @@
+#include <stan/math/fwd.hpp>
+#include <gtest/gtest.h>
+
+#include <cmath>
+
+TEST(AgradFwd, log_softmax_large_dynamic_range) {
+  using stan::math::fvar;
+  using stan::math::log_softmax;
+
+  Eigen::Matrix<fvar<double>, Eigen::Dynamic, 1> x(3);
+
+  x(0) = fvar<double>(0.0, 1.0);
+  x(1) = fvar<double>(-100.0, 2.0);
+  x(2) = fvar<double>(-1000.0, 3.0);
+
+  auto y = log_softmax(x);
+
+  /*
+   * log_softmax(x) = x - log_sum_exp(x).
+   *
+   * exp(-1000) underflows to zero in double precision. Therefore:
+   *
+   *   log(softmax(x)[2])
+   *
+   * gives -inf while
+   *
+   *   x[2] - log_sum_exp(x)
+   *
+   * correctly retains the log-probability.
+   */
+
+  EXPECT_TRUE(std::isfinite(y.val()(0)));
+  EXPECT_TRUE(std::isfinite(y.val()(1)));
+  EXPECT_TRUE(std::isfinite(y.val()(2)));
+
+  EXPECT_DOUBLE_EQ(0.0, y.val()(0));
+  EXPECT_DOUBLE_EQ(-100.0, y.val()(1));
+  EXPECT_DOUBLE_EQ(-1000.0, y.val()(2));
+
+  const double p1 = std::exp(-100.0);
+  const double dot_sd = 1.0 + 2.0 * p1;
+
+  EXPECT_NEAR(1.0 - dot_sd, y(0).d(), 1e-12);
+  EXPECT_NEAR(2.0 - dot_sd, y(1).d(), 1e-12);
+  EXPECT_NEAR(3.0 - dot_sd, y(2).d(), 1e-12);
+}
+
+TEST(AgradFwd, log_softmax_many_small_probabilities) {
+  using stan::math::fvar;
+  using stan::math::log_softmax;
+
+  constexpr int n = 501;
+
+  Eigen::Matrix<fvar<double>, Eigen::Dynamic, 1> x(n);
+
+  x(0) = fvar<double>(0.0, 1.0);
+
+  for (int i = 1; i < n; ++i) {
+    x(i) = fvar<double>(-36.0, 2.0);
+  }
+
+  const auto y = log_softmax(x);
+
+  const double lse = std::log1p(500.0 * std::exp(-36.0));
+  EXPECT_NEAR(-lse, y(0).val(), 1e-14);
+  EXPECT_NEAR(-36.0 - lse, y(1).val(), 1e-14);
+  EXPECT_NEAR(-36.0 - lse, y(500).val(), 1e-14);
+
+  for (int i = 0; i < n; ++i) {
+    EXPECT_TRUE(std::isfinite(y(i).val()));
+    EXPECT_TRUE(std::isfinite(y(i).d()));
+  }
+
+  const double p_small = std::exp(-36.0) / (1.0 + 500.0 * std::exp(-36.0));
+  const double p_large = 1.0 / (1.0 + 500.0 * std::exp(-36.0));
+
+  const double tangent = p_large + 500.0 * 2.0 * p_small;
+
+  EXPECT_NEAR(1.0 - tangent, y(0).d(), 1e-12);
+  EXPECT_NEAR(2.0 - tangent, y(1).d(), 1e-12);
+}
+TEST(AgradFwd, log_softmax_positive_infinity_consistent_with_prim) {
+  using stan::math::fvar;
+  using stan::math::log_softmax;
+
+  constexpr double inf = std::numeric_limits<double>::infinity();
+
+  Eigen::Matrix<fvar<double>, Eigen::Dynamic, 1> x(4);
+  x(0) = fvar<double>(1.0, 1.0);
+  x(1) = fvar<double>(2.0, 2.0);
+  x(2) = fvar<double>(3.0, 3.0);
+  x(3) = fvar<double>(inf, 4.0);
+
+  auto y = log_softmax(x);
+
+  // The value is defined as
+  //
+  //   x - log_sum_exp(x).
+  //
+  // Since log_sum_exp([1, 2, 3, +Inf]) = +Inf:
+  //
+  //   1 - Inf = -Inf
+  //   2 - Inf = -Inf
+  //   3 - Inf = -Inf
+  //   Inf - Inf = NaN
+  //
+
+  EXPECT_EQ(-inf, y(0).val());
+  EXPECT_EQ(-inf, y(1).val());
+  EXPECT_EQ(-inf, y(2).val());
+  EXPECT_TRUE(std::isnan(y(3).val()));
+}
\ No newline at end of file
diff --git a/test/unit/math/fwd/fun/log_sum_exp_test.cpp b/test/unit/math/fwd/fun/log_sum_exp_test.cpp
new file mode 100644
index 00000000000..e6286602bb5
--- /dev/null
+++ b/test/unit/math/fwd/fun/log_sum_exp_test.cpp
@@ -0,0 +1,42 @@
+#include <stan/math/fwd.hpp>
+#include <gtest/gtest.h>
+
+#include <cmath>
+
+TEST(AgradFwd, log_sum_exp_derivative_does_not_overflow) {
+  using stan::math::fvar;
+  using stan::math::log_sum_exp;
+
+  Eigen::Matrix<fvar<double>, Eigen::Dynamic, 1> x(3);
+  x(0) = fvar<double>(1000.0, 1.0);
+  x(1) = fvar<double>(1001.0, 2.0);
+  x(2) = fvar<double>(999.0, 3.0);
+
+  auto y = log_sum_exp(x);
+  // log_sum_exp([1000, 1001, 999])
+  //
+  // = 1001 + log(exp(-1) + 1 + exp(-2)).
+  const double expected_value
+      = 1001.0 + std::log(std::exp(-1.0) + 1.0 + std::exp(-2.0));
+
+  EXPECT_NEAR(expected_value, y.val(), 1e-12);
+  EXPECT_TRUE(std::isfinite(y.val()));
+
+  // The derivative is
+  //
+  // softmax(x) dot x.d().
+  //
+  // The old implementation formed exp(x) directly in the derivative,
+  // which gives inf / inf for these values.
+
+  const double denom = std::exp(-1.0) + 1.0 + std::exp(-2.0);
+
+  const double p0 = std::exp(-1.0) / denom;
+  const double p1 = 1.0 / denom;
+  const double p2 = std::exp(-2.0) / denom;
+
+  const double expected_derivative = p0 + 2.0 * p1 + 3.0 * p2;
+
+  EXPECT_TRUE(std::isfinite(y.d()));
+  EXPECT_NEAR(expected_derivative, y.d(), 1e-12);
+}
diff --git a/test/unit/math/rev/fun/log_softmax_test.cpp b/test/unit/math/rev/fun/log_softmax_test.cpp
index 0008e372b4c..e847af5e04a 100644
--- a/test/unit/math/rev/fun/log_softmax_test.cpp
+++ b/test/unit/math/rev/fun/log_softmax_test.cpp
@@ -60,3 +60,61 @@ TEST_F(AgradRev, log_softmax_var_value_row_vector) {
   EXPECT_FLOAT_EQ(1.0 - 3.0 * std::exp(1.0) / denom, x.adj()(1));
   EXPECT_FLOAT_EQ(1.0 - 3.0 * std::exp(2.0) / denom, x.adj()(2));
 }
+
+TEST_F(AgradRev, log_softmax_negative_infinity_has_finite_adjoint) {
+  using stan::math::log_softmax;
+  using stan::math::sum;
+  using stan::math::var_value;
+
+  const double neg_inf = -std::numeric_limits<double>::infinity();
+
+  Eigen::VectorXd v(4);
+  v << neg_inf, 1.0, 2.0, 3.0;
+
+  var_value<Eigen::VectorXd> x(v);
+
+  auto y = log_softmax(x);
+  sum(y).grad();
+
+  const double e1 = std::exp(1.0);
+  const double e2 = std::exp(2.0);
+  const double e3 = std::exp(3.0);
+  const double denom = e1 + e2 + e3;
+
+  const double p0 = 0.0;
+  const double p1 = e1 / denom;
+  const double p2 = e2 / denom;
+  const double p3 = e3 / denom;
+
+  EXPECT_NEAR(1.0, x.adj()(0), 1e-12);
+  EXPECT_NEAR(1.0 - 4.0 * p1, x.adj()(1), 1e-12);
+  EXPECT_NEAR(1.0 - 4.0 * p2, x.adj()(2), 1e-12);
+  EXPECT_NEAR(1.0 - 4.0 * p3, x.adj()(3), 1e-12);
+}
+
+TEST_F(AgradRev, log_softmax_adjoint_uses_stable_softmax) {
+  using stan::math::log_softmax;
+  using stan::math::sum;
+  using stan::math::var_value;
+
+  Eigen::VectorXd v(4);
+
+  v << 629.7901581243797, 31.52411463, 608.19720553, 120.94829574;
+
+  var_value<Eigen::VectorXd> x(v);
+  auto y = log_softmax(x);
+
+  sum(y).grad();
+
+  // Independent max-shifted softmax calculation
+  Eigen::VectorXd p = (v.array() - v.maxCoeff()).exp();
+  p /= p.sum();
+
+  // d/dx_m sum_k log_softmax(x)_k
+  // = 1 - n * softmax(x)_m.
+  for (Eigen::Index i = 0; i < v.size(); ++i) {
+    const double expected = 1.0 - v.size() * p(i);
+
+    EXPECT_NEAR(expected, x.adj()(i), 1e-12);
+  }
+}
diff --git a/test/unit/math/rev/fun/log_sum_exp_test.cpp b/test/unit/math/rev/fun/log_sum_exp_test.cpp
index f020a0d1dd2..c5f40e51e5c 100644
--- a/test/unit/math/rev/fun/log_sum_exp_test.cpp
+++ b/test/unit/math/rev/fun/log_sum_exp_test.cpp
@@ -61,3 +61,51 @@ TEST_F(AgradRev, log_sum_exp_tests_large_values) {
   output5.grad();
   EXPECT_FLOAT_EQ(a8.adj(), 0.0);
 }
+
+TEST_F(AgradRev, log_sum_exp_negative_infinity_has_zero_adjoint) {
+  using stan::math::log_sum_exp;
+  using stan::math::var_value;
+
+  const double neg_inf = -std::numeric_limits<double>::infinity();
+
+  Eigen::VectorXd v(4);
+  v << neg_inf, 1.0, 2.0, 3.0;
+
+  var_value<Eigen::VectorXd> x(v);
+
+  auto y = log_sum_exp(x);
+  y.grad();
+
+  const double e1 = std::exp(1.0);
+  const double e2 = std::exp(2.0);
+  const double e3 = std::exp(3.0);
+  const double denom = e1 + e2 + e3;
+
+  EXPECT_DOUBLE_EQ(0.0, x.adj()(0));
+  EXPECT_NEAR(e1 / denom, x.adj()(1), 1e-12);
+  EXPECT_NEAR(e2 / denom, x.adj()(2), 1e-12);
+  EXPECT_NEAR(e3 / denom, x.adj()(3), 1e-12);
+}
+
+TEST_F(AgradRev, log_sum_exp_adjoint_uses_stable_softmax) {
+  using stan::math::log_sum_exp;
+  using stan::math::var_value;
+
+  Eigen::VectorXd v(4);
+
+  v << 629.7901581243797, 31.52411463, 608.19720553, 120.94829574;
+
+  var_value<Eigen::VectorXd> x(v);
+
+  auto y = log_sum_exp(x);
+  y.grad();
+
+  // Independent max-shifted softmax calculation
+  Eigen::VectorXd p = (v.array() - v.maxCoeff()).exp();
+  p /= p.sum();
+
+  // The gradient of log_sum_exp is softmax.
+  for (Eigen::Index i = 0; i < v.size(); ++i) {
+    EXPECT_NEAR(p(i), x.adj()(i), 1e-12);
+  }
+}