Feature: L1 and L2 norms

stan/math/prim/fun.hpp

@@ -227,6 +227,8 @@
 #include <stan/math/prim/fun/multiply_log.hpp>
 #include <stan/math/prim/fun/multiply_lower_tri_self_transpose.hpp>
 #include <stan/math/prim/fun/norm.hpp>
+#include <stan/math/prim/fun/norm1.hpp>
+#include <stan/math/prim/fun/norm2.hpp>
 #include <stan/math/prim/fun/num_elements.hpp>
 #include <stan/math/prim/fun/offset_multiplier_constrain.hpp>
 #include <stan/math/prim/fun/offset_multiplier_free.hpp>

stan/math/prim/fun/norm1.hpp

@@ -0,0 +1,30 @@
+#ifndef STAN_MATH_PRIM_FUN_NORM1_HPP
+#define STAN_MATH_PRIM_FUN_NORM1_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <cstddef>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns L1 norm of a vector. For vectors that equals the
+ * sum of magnitudes of its individual elements.
+ *
+ * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase)
+ * @param v Vector.
+ * @return L1 norm of v.
+ */
+template <typename Container, require_container_t<Container>* = nullptr,
+          require_not_st_var<Container>* = nullptr>
+inline auto norm1(const Container& x) {
+  return apply_vector_unary<ref_type_t<Container>>::reduce(
+      to_ref(x), [](const auto& v) { return v.template lpNorm<1>(); });
+}
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/prim/fun/norm2.hpp

@@ -0,0 +1,31 @@
+#ifndef STAN_MATH_PRIM_FUN_NORM2_HPP
+#define STAN_MATH_PRIM_FUN_NORM2_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <cstddef>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns L2 norm of a vector. For vectors that equals the square-root of the
+ * sum of squares of the elements.
+ *
+ * @tparam T type of the vector (must be derived from \c Eigen::MatrixBase)
+ * @param v Vector.
+ * @return L2 norm of v.
+ */
+template <typename Container, require_container_t<Container>* = nullptr,
+          require_not_st_var<Container>* = nullptr>
+inline auto norm2(const Container& x) {
+  return apply_vector_unary<ref_type_t<Container>>::reduce(
+      to_ref(x), [](const auto& v) { return v.template lpNorm<2>(); });
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif

stan/math/rev/fun.hpp

@@ -130,6 +130,8 @@
 #include <stan/math/rev/fun/multiply_log.hpp>
 #include <stan/math/rev/fun/multiply_lower_tri_self_transpose.hpp>
 #include <stan/math/rev/fun/norm.hpp>
+#include <stan/math/rev/fun/norm1.hpp>
+#include <stan/math/rev/fun/norm2.hpp>
 #include <stan/math/rev/fun/ordered_constrain.hpp>
 #include <stan/math/rev/fun/owens_t.hpp>
 #include <stan/math/rev/fun/polar.hpp>

stan/math/rev/fun/norm1.hpp

@@ -0,0 +1,60 @@
+#ifndef STAN_MATH_REV_FUN_NORM1_HPP
+#define STAN_MATH_REV_FUN_NORM1_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/core/typedefs.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/sign.hpp>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the L1 norm of a vector of var.
+ *
+ * @tparam T type of the vector (must have one compile-time dimension equal to
+ * 1)
+ * @param[in] v Vector.
+ * @return L1 norm of v.
+ */
+template <typename T, require_eigen_vector_vt<is_var, T>* = nullptr>
+inline var norm1(const T& v) {
+  const auto& v_ref = to_ref(v);
+  arena_t<T> arena_v(v_ref.size());
+  var res(0.);
+  for (size_t i = 0; i < arena_v.size(); ++i) {
+    arena_v.coeffRef(i) = v_ref.coeffRef(i);
+    res += abs(arena_v.coeffRef(i).val());
+  }
+
+  reverse_pass_callback([res, arena_v]() mutable {
+    arena_v.adj().array() += res.adj() * sign(arena_v.val().array());
+  });
+
+  return res;
+}
+
+/**
+ * Returns the L1 norm of a `var_value<Vector>`.
+ *
+ * @tparam A `var_value<>` whose inner type has one compile-time row or column.
+ * @param[in] v Vector.
+ * @return L1 norm of v.
+ */
+template <typename T, require_var_matrix_t<T>* = nullptr>
+inline var norm1(const T& v) {
+  var res = v.val().array().abs().sum();
+
+  reverse_pass_callback([res, v]() mutable {
+    v.adj().array() += res.adj() * sign(v.val().array());
+  });
+
+  return res;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

stan/math/rev/fun/norm2.hpp

@@ -0,0 +1,61 @@
+#ifndef STAN_MATH_REV_FUN_NORM2_HPP
+#define STAN_MATH_REV_FUN_NORM2_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/rev/core/typedefs.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/sign.hpp>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Returns the L2 norm of a vector of var.
+ *
+ * @tparam T type of the vector (must have one compile-time dimension equal to
+ * 1)
+ * @param[in] v Vector.
+ * @return L2 norm of v.
+ */
+template <typename T, require_eigen_vector_vt<is_var, T>* = nullptr>
+inline var norm2(const T& v) {
+  const auto& v_ref = to_ref(v);
+  arena_t<T> arena_v(v_ref.size());
+  value_type_t<T> res_val = 0;
+
+  for (size_t i = 0; i < arena_v.size(); ++i) {
+    arena_v.coeffRef(i) = v_ref.coeffRef(i);
+    res_val += square(arena_v.coeffRef(i).val());
+  }
+  var res(sqrt(res_val));
+  reverse_pass_callback([res, arena_v]() mutable {
+    arena_v.adj().array() += res.adj() * (arena_v.val().array() / res.val());
+  });
+
+  return res;
+}
+
+/**
+ * Returns the L2 norm of a `var_value<Vector>`.
+ *
+ * @tparam A `var_value<>` whose inner type has one compile-time row or column.
+ * @param[in] v Vector.
+ * @return L2 norm of v.
+ */
+template <typename T, require_var_matrix_t<T>* = nullptr>
+inline var norm2(const T& v) {
+  var res = sqrt(v.val().array().square().sum());
+
+  reverse_pass_callback([res, v]() mutable {
+    v.adj().array() += res.adj() * (v.val().array() / res.val());
+  });
+
+  return res;
+}
+
+}  // namespace math
+}  // namespace stan
+#endif

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

@@ -0,0 +1,22 @@
+#include <test/unit/math/test_ad.hpp>
+#include <vector>
+
+TEST(MathMixMatFun, norm1) {
+  auto f = [](const auto& y) { return stan::math::norm1(y); };
+
+  Eigen::VectorXd x0(0);
+
+  Eigen::VectorXd x1(1);
+  x1 << 2;
+
+  Eigen::VectorXd x2(2);
+  x2 << 2, 3;
+
+  Eigen::VectorXd x3(3);
+  x3 << 2, 3, 4;
+
+  for (const auto& a : std::vector<Eigen::VectorXd>{x0, x1, x2, x3}) {
+    stan::test::expect_ad(f, a);
+    stan::test::expect_ad_matvar(f, a);
+  }
+}

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

@@ -0,0 +1,22 @@
+#include <test/unit/math/test_ad.hpp>
+#include <vector>
+
+TEST(MathMixMatFun, norm2) {
+  auto f = [](const auto& y) { return stan::math::norm2(y); };
+
+  Eigen::VectorXd x0(0);
+
+  Eigen::VectorXd x1(1);
+  x1 << 2;
+
+  Eigen::VectorXd x2(2);
+  x2 << 2, 3;
+
+  Eigen::VectorXd x3(3);
+  x3 << 2, 3, 4;
+
+  for (const auto& a : std::vector<Eigen::VectorXd>{x0, x1, x2, x3}) {
+    stan::test::expect_ad(f, a);
+    stan::test::expect_ad_matvar(f, a);
+  }
+}

test/unit/math/prim/fun/norm1_test.cpp

@@ -0,0 +1,76 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <cmath>
+#include <limits>
+#include <vector>
+
+TEST(MathFunctions, norm1) {
+  std::vector<double> x(3), y(3);
+  x[0] = 2.33;
+  x[1] = 8.88;
+  x[2] = 9.81;
+  y[0] = 2.46;
+  y[1] = 4.45;
+  y[2] = 1.03;
+
+  EXPECT_FLOAT_EQ(21.02, stan::math::norm1(x));
+  EXPECT_FLOAT_EQ(7.94, stan::math::norm1(y));
+}
+
+TEST(MathFunctions, norm1_nan) {
+  std::vector<double> x(3);
+  x[0] = 2.33;
+  x[1] = 8.88;
+  x[2] = 9.81;
+
+  double nan = std::numeric_limits<double>::quiet_NaN();
+  x[2] = nan;
+
+  EXPECT_TRUE(std::isnan(stan::math::norm1(x)));
+
+  x[0] = nan;
+  x[1] = nan;
+  x[2] = nan;
+  EXPECT_TRUE(std::isnan(stan::math::norm1(x)));
+}
+
+TEST(MathMatrixPrimMat, norm1) {
+  using stan::math::norm1;
+
+  Eigen::Matrix<double, Eigen::Dynamic, 1> v1(1);
+  v1 << 2.0;
+  EXPECT_NEAR(2.0, norm1(v1), 1E-12);
+  Eigen::Matrix<double, Eigen::Dynamic, 1> v2(2);
+  v2 << 2.0, 3.0;
+  EXPECT_NEAR(5.0, norm1(v2), 1E-12);
+  Eigen::Matrix<double, Eigen::Dynamic, 1> v3(3);
+  v3 << 2.0, 3.0, 4.0;
+  EXPECT_NEAR(9.0, norm1(v3), 1E-12);
+
+  Eigen::Matrix<double, 1, Eigen::Dynamic> rv1(1);
+  rv1 << 2.0;
+  EXPECT_NEAR(2.0, norm1(rv1), 1E-12);
+  Eigen::Matrix<double, 1, Eigen::Dynamic> rv2(2);
+  rv2 << 2.0, 3.0;
+  EXPECT_NEAR(5.0, norm1(rv2), 1E-12);
+  Eigen::Matrix<double, 1, Eigen::Dynamic> rv3(3);
+  rv3 << 2.0, 3.0, 4.0;
+  EXPECT_NEAR(9.0, norm1(rv3), 1E-12);
+
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> m1(1, 1);
+  m1 << 2.0;
+  EXPECT_NEAR(2.0, norm1(m1), 1E-12);
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> m2(2, 1);
+  m2 << 2.0, 3.0;
+  EXPECT_NEAR(5.0, norm1(m2), 1E-12);
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> m3(3, 1);
+  m3 << 2.0, 3.0, 4.0;
+  EXPECT_NEAR(9.0, norm1(m3), 1E-12);
+
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> mm2(1, 2);
+  mm2 << 2.0, 3.0;
+  EXPECT_NEAR(5.0, norm1(mm2), 1E-12);
+  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> mm3(1, 3);
+  mm3 << 2.0, 3.0, 4.0;
+  EXPECT_NEAR(9.0, norm1(mm3), 1E-12);
+}