diff --git a/.gitignore b/.gitignore
index c36591f03d8..63af199a251 100644
--- a/.gitignore
+++ b/.gitignore
@@ -106,3 +106,4 @@ lib/tbb/**
 
 # gdb debug history
 .gdb_history
+.positai
diff --git a/stan/math/fwd/fun/log_softmax.hpp b/stan/math/fwd/fun/log_softmax.hpp
index a9685d7ca00..8b96e9a231c 100644
--- a/stan/math/fwd/fun/log_softmax.hpp
+++ b/stan/math/fwd/fun/log_softmax.hpp
@@ -4,9 +4,9 @@
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/meta.hpp>
-#include <stan/math/fwd/fun/softmax.hpp>
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/log_softmax.hpp>
+#include <stan/math/prim/fun/softmax_log_sum_exp.hpp>
 #include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/functor/apply_vector_unary.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,29 +29,24 @@ 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));
-  if (x_ref.size() == 0) {
-    return Eigen::Matrix<fvar<T>, Rows, Cols>{};
+template <typename Mat, require_eigen_vt<is_fvar, Mat>* = nullptr>
+inline plain_type_t<Mat> log_softmax(Mat&& x) {
+  if (x.size() == 0) {
+    return {};
   }
-  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());
-  result.val() = s.array().log().matrix();
-  result.d() = (d_in.array() - dot_sd).matrix();
-  return result;
+  decltype(auto) x_ref = to_ref(std::forward<Mat>(x));
+  const auto x_val = value_of(x_ref).eval();
+  const auto d_in = x_ref.d().eval();
+  const auto softmax_lse = internal::softmax_log_sum_exp(x_val);
+  return to_fvar(x_val.array() - softmax_lse.second,
+                 d_in.array() - softmax_lse.first.cwiseProduct(d_in).sum());
 }
 
 }  // namespace math
diff --git a/stan/math/fwd/fun/log_sum_exp.hpp b/stan/math/fwd/fun/log_sum_exp.hpp
index 926838d143c..1e5f01a2259 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_log_sum_exp.hpp>
 #include <cmath>
 #include <vector>
 
@@ -56,11 +57,9 @@ 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 softmax_lse = internal::softmax_log_sum_exp(vals);
+        return fvar<T_fvar_inner>(softmax_lse.second,
+                                  v.d().cwiseProduct(softmax_lse.first).sum());
       });
 }
 
diff --git a/stan/math/fwd/fun/softmax.hpp b/stan/math/fwd/fun/softmax.hpp
index 97d019b6c23..34e35adbc80 100644
--- a/stan/math/fwd/fun/softmax.hpp
+++ b/stan/math/fwd/fun/softmax.hpp
@@ -29,27 +29,21 @@ 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));
-  if (x_ref.size() == 0) {
-    return Eigen::Matrix<fvar<T>, Rows, Cols>{};
+template <typename Mat, require_eigen_vt<is_fvar, Mat>* = nullptr>
+inline plain_type_t<Mat> softmax(Mat&& x) {
+  if (x.size() == 0) {
+    return {};
   }
-  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());
-  result.val() = s;
-  result.d() = (s.array() * (d_in.array() - dot_sd)).matrix();
-  return result;
+  decltype(auto) x_ref = to_ref(std::forward<Mat>(x));
+  const auto s = softmax(value_of(x_ref)).eval();
+  const auto d_in = x_ref.d().eval();
+
+  return to_fvar(s, s.array() * (d_in.array() - s.cwiseProduct(d_in).sum()));
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp
index 8ef4204c81b..4d3c3d017d3 100644
--- a/stan/math/prim/fun.hpp
+++ b/stan/math/prim/fun.hpp
@@ -296,6 +296,7 @@
 #include <stan/math/prim/fun/size_mvt.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
 #include <stan/math/prim/fun/softmax.hpp>
+#include <stan/math/prim/fun/softmax_log_sum_exp.hpp>
 #include <stan/math/prim/fun/sort_asc.hpp>
 #include <stan/math/prim/fun/sort_desc.hpp>
 #include <stan/math/prim/fun/sort_indices.hpp>
diff --git a/stan/math/prim/fun/log_softmax.hpp b/stan/math/prim/fun/log_softmax.hpp
index d6c9017edb6..c3760135df7 100644
--- a/stan/math/prim/fun/log_softmax.hpp
+++ b/stan/math/prim/fun/log_softmax.hpp
@@ -12,8 +12,9 @@ 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)
@@ -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
+ * @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 log softmax of the input, preserving the container structure; an
- *   empty result if any input vector is empty
+ *   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..6cac82507be 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) {
@@ -58,7 +58,7 @@ inline auto softmax(Container&& x) {
               if (v.size() == 0) {
                 return v;
               }
-              const auto theta = (v.array() - v.maxCoeff()).exp();
+              const auto theta = (v.array() - v.maxCoeff()).exp().eval();
               return (theta / theta.sum()).matrix();
             });
       },
diff --git a/stan/math/prim/fun/softmax_log_sum_exp.hpp b/stan/math/prim/fun/softmax_log_sum_exp.hpp
new file mode 100644
index 00000000000..d4ac4f094d9
--- /dev/null
+++ b/stan/math/prim/fun/softmax_log_sum_exp.hpp
@@ -0,0 +1,69 @@
+#ifndef STAN_MATH_PRIM_FUN_SOFTMAX_LOG_SUM_EXP_HPP
+#define STAN_MATH_PRIM_FUN_SOFTMAX_LOG_SUM_EXP_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/constants.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/log_sum_exp.hpp>
+#include <stan/math/prim/fun/softmax.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <cmath>
+#include <utility>
+
+namespace stan {
+namespace math {
+namespace internal {
+
+/**
+ * Return the softmax and the log sum of exponentials of the specified
+ * vector or matrix, using a single evaluation of the exponential per
+ * element. For a matrix, both are taken over all elements.
+ *
+ * With \f$m = \max(x)\f$ and \f$\theta = \exp(x - m)\f$,
+ *
+ * \f$\mbox{softmax}(x) = \theta / \sum_{n=1}^N \theta_n\f$ and
+ * \f$\mbox{log\_sum\_exp}(x) = m + \log \sum_{n=1}^N \theta_n\f$.
+ *
+ * The results match `softmax(x)` and `log_sum_exp(x)`, including for
+ * non-finite input: if \f$m\f$ is not finite, the log sum of exponentials
+ * is \f$m\f$ and the softmax is NaN.
+ *
+ * @tparam T type of Eigen vector or matrix
+ * @param x Eigen vector or matrix
+ * @return pair of the softmax of `x` and the log sum of exponentials of
+ *   `x`. For empty input, returns an empty softmax and negative infinity
+ */
+template <typename T, require_eigen_vt<std::is_arithmetic, T>* = nullptr>
+inline std::pair<plain_type_t<T>, double> softmax_log_sum_exp(const T& x) {
+  if (x.size() == 0) {
+    return {x, NEGATIVE_INFTY};
+  }
+  const auto& x_ref = to_ref(x);
+  const double max = x_ref.maxCoeff();
+  const auto theta = (x_ref.array() - max).exp().eval();
+  const double theta_sum = theta.sum();
+  const double lse = std::isfinite(max) ? max + std::log(theta_sum) : max;
+  return {plain_type_t<T>((theta / theta_sum).matrix()), lse};
+}
+
+/**
+ * Return the softmax and the log sum of exponentials of the specified
+ * vector or matrix of autodiff values. This happens with nested autodiff,
+ * e.g. the values of an `fvar<var>` vector, and falls back to the
+ * autodiff specializations of `softmax()` and `log_sum_exp()`.
+ *
+ * @tparam T type of Eigen vector or matrix with autodiff scalar
+ * @param x vector or matrix
+ * @return pair of the softmax of `x` and the log sum of exponentials of `x`
+ */
+template <typename T, require_eigen_vt<is_autodiff, T>* = nullptr>
+inline std::pair<plain_type_t<T>, value_type_t<T>> softmax_log_sum_exp(
+    const T& x) {
+  return {plain_type_t<T>(softmax(x)), log_sum_exp(x)};
+}
+
+}  // namespace internal
+}  // namespace math
+}  // namespace stan
+
+#endif
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..fc7f770c5d1 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_log_sum_exp.hpp>
 #include <cmath>
 #include <vector>
 
@@ -68,11 +69,11 @@ template <typename T, require_eigen_st<is_var, T>* = nullptr,
           require_not_var_matrix_t<T>* = nullptr>
 inline var log_sum_exp(T&& v) {
   arena_t<T> arena_v = std::forward<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();
+  const auto softmax_lse = internal::softmax_log_sum_exp(arena_v.val());
+  auto arena_softmax = to_arena(softmax_lse.first);
+  var res = softmax_lse.second;
+  reverse_pass_callback([arena_v, arena_softmax, res]() mutable {
+    arena_v.adj() += res.adj() * arena_softmax;
   });
 
   return res;
@@ -86,9 +87,12 @@ 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();
-  });
+  const auto softmax_lse = internal::softmax_log_sum_exp(x.val());
+  auto arena_softmax = to_arena(softmax_lse.first);
+  return make_callback_vari(softmax_lse.second,
+                            [x, arena_softmax](const auto& res) mutable {
+                              x.adj() += res.adj() * arena_softmax;
+                            });
 }
 
 /**
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..dd713ed1641
--- /dev/null
+++ b/test/unit/math/fwd/fun/log_softmax_test.cpp
@@ -0,0 +1,79 @@
+#include <stan/math/fwd.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+
+#include <cmath>
+
+TEST(AgradFwd, log_softmax_large_dynamic_range) {
+  using stan::math::fvar;
+  using stan::math::log_softmax;
+  using stan::math::vector_fd;
+
+  Eigen::VectorXd val(3);
+  val << 0, -100, -1000;
+
+  Eigen::VectorXd d(3);
+  d << 1, 2, 3;
+
+  vector_fd x = stan::math::to_fvar(val, d);
+
+  auto y = log_softmax(x);
+
+  EXPECT_MATRIX_EQ(y.val(), val);
+  EXPECT_MATRIX_EQ(y.d(), d.array() - (1.0 + 2.0 * std::exp(-100.0)));
+}
+
+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);
+
+  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()));
+}
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..3da59876486
--- /dev/null
+++ b/test/unit/math/fwd/fun/log_sum_exp_test.cpp
@@ -0,0 +1,32 @@
+#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);
+  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()));
+
+  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/prim/fun/softmax_log_sum_exp_test.cpp b/test/unit/math/prim/fun/softmax_log_sum_exp_test.cpp
new file mode 100644
index 00000000000..4a17014a3a6
--- /dev/null
+++ b/test/unit/math/prim/fun/softmax_log_sum_exp_test.cpp
@@ -0,0 +1,59 @@
+#include <stan/math/prim.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+#include <cmath>
+#include <limits>
+
+TEST(MathFunctions, softmax_log_sum_exp_matches_separate_calls) {
+  using stan::math::log_sum_exp;
+  using stan::math::softmax;
+  using stan::math::internal::softmax_log_sum_exp;
+
+  Eigen::VectorXd v(4);
+  v << 629.7901581243797, 31.52411463, 608.19720553, 120.94829574;
+  auto res = softmax_log_sum_exp(v);
+  EXPECT_MATRIX_EQ(softmax(v), res.first);
+  EXPECT_FLOAT_EQ(log_sum_exp(v), res.second);
+
+  Eigen::RowVectorXd rv(3);
+  rv << -1.0, 1.0, 10.0;
+  auto res_rv = softmax_log_sum_exp(rv);
+  EXPECT_MATRIX_EQ(softmax(rv), res_rv.first);
+  EXPECT_FLOAT_EQ(log_sum_exp(rv), res_rv.second);
+
+  Eigen::MatrixXd m(2, 2);
+  m << 0.0, -100.0, -1000.0, 3.0;
+  auto res_m = softmax_log_sum_exp(m);
+  EXPECT_MATRIX_EQ(softmax(m), res_m.first);
+  EXPECT_FLOAT_EQ(log_sum_exp(m), res_m.second);
+}
+
+TEST(MathFunctions, softmax_log_sum_exp_empty) {
+  using stan::math::internal::softmax_log_sum_exp;
+
+  Eigen::VectorXd v(0);
+  auto res = softmax_log_sum_exp(v);
+  EXPECT_EQ(0, res.first.size());
+  EXPECT_EQ(stan::math::NEGATIVE_INFTY, res.second);
+}
+
+TEST(MathFunctions, softmax_log_sum_exp_non_finite) {
+  using stan::math::internal::softmax_log_sum_exp;
+  constexpr double inf = std::numeric_limits<double>::infinity();
+
+  Eigen::VectorXd v(3);
+  v << 1.0, -inf, 2.0;
+  auto res = softmax_log_sum_exp(v);
+  EXPECT_FLOAT_EQ(2.0 + std::log(std::exp(-1.0) + 1.0), res.second);
+  EXPECT_EQ(0.0, res.first(1));
+
+  v << 1.0, inf, 2.0;
+  res = softmax_log_sum_exp(v);
+  EXPECT_EQ(inf, res.second);
+  EXPECT_TRUE(res.first.array().isNaN().all());
+
+  v << -inf, -inf, -inf;
+  res = softmax_log_sum_exp(v);
+  EXPECT_EQ(-inf, res.second);
+  EXPECT_TRUE(res.first.array().isNaN().all());
+}
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);
+  }
+}