From ff161152bdcee4e637d075d787618f3869a94eb2 Mon Sep 17 00:00:00 2001
From: Marco Colombo <m.colombo@ed.ac.uk>
Date: Tue, 11 Feb 2020 08:39:29 +0000
Subject: [PATCH 1/8] Split quad_form and quad_form_sym tests

---
 .../unit/math/prim/fun/quad_form_sym_test.cpp | 74 +++++++++++++++++
 test/unit/math/prim/fun/quad_form_test.cpp    | 79 +------------------
 2 files changed, 76 insertions(+), 77 deletions(-)
 create mode 100644 test/unit/math/prim/fun/quad_form_sym_test.cpp

diff --git a/test/unit/math/prim/fun/quad_form_sym_test.cpp b/test/unit/math/prim/fun/quad_form_sym_test.cpp
new file mode 100644
index 00000000000..528dbba97b7
--- /dev/null
+++ b/test/unit/math/prim/fun/quad_form_sym_test.cpp
@@ -0,0 +1,74 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+
+TEST(MathMatrixPrim, quad_form_sym_mat) {
+  using stan::math::matrix_d;
+
+  matrix_d ad(4, 4);
+  matrix_d bd(4, 2);
+
+  bd << 100, 10, 0, 1, -3, -3, 5, 2;
+  ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
+      1.0, 112.0;
+
+  // double-double
+  matrix_d resd = stan::math::quad_form_sym(ad, bd);
+  EXPECT_FLOAT_EQ(25433, resd(0, 0));
+  EXPECT_FLOAT_EQ(3396, resd(0, 1));
+  EXPECT_FLOAT_EQ(3396, resd(1, 0));
+  EXPECT_FLOAT_EQ(725, resd(1, 1));
+}
+
+TEST(MathMatrixPrim, quad_form_sym_vec) {
+  using stan::math::matrix_d;
+  using stan::math::vector_d;
+
+  matrix_d ad(4, 4);
+  vector_d bd(4);
+  double res;
+
+  bd << 100, 0, -3, 5;
+  ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
+      1.0, 112.0;
+
+  // double-double
+  res = stan::math::quad_form_sym(ad, bd);
+  EXPECT_FLOAT_EQ(25433, res);
+}
+
+TEST(MathMatrixPrim, quad_form_sym_symmetry) {
+  using stan::math::matrix_d;
+  using stan::math::quad_form_sym;
+
+  matrix_d ad(4, 4);
+  matrix_d bd(4, 2);
+
+  bd << 100, 10, 0, 1, -3, -3, 5, 2;
+  ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
+      1.0, 112.0;
+
+  // double-double
+  matrix_d resd = quad_form_sym(ad, bd);
+  EXPECT_EQ(resd(1, 0), resd(0, 1));
+
+  bd.resize(4, 3);
+  bd << 100, 10, 11, 0, 1, 12, -3, -3, 34, 5, 2, 44;
+  resd = quad_form_sym(ad, bd);
+  EXPECT_EQ(resd(1, 0), resd(0, 1));
+  EXPECT_EQ(resd(2, 0), resd(0, 2));
+  EXPECT_EQ(resd(2, 1), resd(1, 2));
+}
+
+TEST(MathMatrixPrim, quad_form_sym_asymmetric) {
+  using stan::math::matrix_d;
+
+  matrix_d ad(4, 4);
+  matrix_d bd(4, 2);
+
+  bd << 100, 10, 0, 1, -3, -3, 5, 2;
+  ad << 2.0, 3.0, 4.0, 5.0, 6.0, 10.0, 2.0, 2.0, 7.0, 2.0, 7.0, 1.0, 8.0, 2.0,
+      1.0, 112.0;
+
+  // double-double
+  EXPECT_THROW(stan::math::quad_form_sym(ad, bd), std::domain_error);
+}
diff --git a/test/unit/math/prim/fun/quad_form_test.cpp b/test/unit/math/prim/fun/quad_form_test.cpp
index 3af7adeace0..09f82e036bf 100644
--- a/test/unit/math/prim/fun/quad_form_test.cpp
+++ b/test/unit/math/prim/fun/quad_form_test.cpp
@@ -1,7 +1,7 @@
 #include <stan/math/prim.hpp>
 #include <gtest/gtest.h>
 
-TEST(MathMatrixPrimMat, quad_form_mat) {
+TEST(MathMatrixPrim, quad_form_mat) {
   using stan::math::matrix_d;
   using stan::math::quad_form;
 
@@ -20,26 +20,7 @@ TEST(MathMatrixPrimMat, quad_form_mat) {
   EXPECT_FLOAT_EQ(725, resd(1, 1));
 }
 
-TEST(MathMatrixPrimMat, quad_form_sym_mat) {
-  using stan::math::matrix_d;
-  using stan::math::quad_form_sym;
-
-  matrix_d ad(4, 4);
-  matrix_d bd(4, 2);
-
-  bd << 100, 10, 0, 1, -3, -3, 5, 2;
-  ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
-      1.0, 112.0;
-
-  // double-double
-  matrix_d resd = quad_form_sym(ad, bd);
-  EXPECT_FLOAT_EQ(25433, resd(0, 0));
-  EXPECT_FLOAT_EQ(3396, resd(0, 1));
-  EXPECT_FLOAT_EQ(3396, resd(1, 0));
-  EXPECT_FLOAT_EQ(725, resd(1, 1));
-}
-
-TEST(MathMatrixPrimMat, quad_form_vec) {
+TEST(MathMatrixPrim, quad_form_vec) {
   using stan::math::matrix_d;
   using stan::math::quad_form;
   using stan::math::vector_d;
@@ -56,59 +37,3 @@ TEST(MathMatrixPrimMat, quad_form_vec) {
   res = quad_form(ad, bd);
   EXPECT_FLOAT_EQ(26033, res);
 }
-
-TEST(MathMatrixPrimMat, quad_form_sym_vec) {
-  using stan::math::matrix_d;
-  using stan::math::quad_form_sym;
-  using stan::math::vector_d;
-
-  matrix_d ad(4, 4);
-  vector_d bd(4);
-  double res;
-
-  bd << 100, 0, -3, 5;
-  ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
-      1.0, 112.0;
-
-  // double-double
-  res = quad_form_sym(ad, bd);
-  EXPECT_FLOAT_EQ(25433, res);
-}
-
-TEST(MathMatrixPrimMat, quad_form_sym_symmetry) {
-  using stan::math::matrix_d;
-  using stan::math::quad_form_sym;
-
-  matrix_d ad(4, 4);
-  matrix_d bd(4, 2);
-
-  bd << 100, 10, 0, 1, -3, -3, 5, 2;
-  ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
-      1.0, 112.0;
-
-  // double-double
-  matrix_d resd = quad_form_sym(ad, bd);
-  EXPECT_EQ(resd(1, 0), resd(0, 1));
-
-  bd.resize(4, 3);
-  bd << 100, 10, 11, 0, 1, 12, -3, -3, 34, 5, 2, 44;
-  resd = quad_form_sym(ad, bd);
-  EXPECT_EQ(resd(1, 0), resd(0, 1));
-  EXPECT_EQ(resd(2, 0), resd(0, 2));
-  EXPECT_EQ(resd(2, 1), resd(1, 2));
-}
-
-TEST(MathMatrixPrimMat, quad_form_sym_asymmetric) {
-  using stan::math::matrix_d;
-  using stan::math::quad_form_sym;
-
-  matrix_d ad(4, 4);
-  matrix_d bd(4, 2);
-
-  bd << 100, 10, 0, 1, -3, -3, 5, 2;
-  ad << 2.0, 3.0, 4.0, 5.0, 6.0, 10.0, 2.0, 2.0, 7.0, 2.0, 7.0, 1.0, 8.0, 2.0,
-      1.0, 112.0;
-
-  // double-double
-  EXPECT_THROW(quad_form_sym(ad, bd), std::domain_error);
-}

From 30b042d3542e1ffe36ee19ef544ccee347ba0857 Mon Sep 17 00:00:00 2001
From: Marco Colombo <m.colombo@ed.ac.uk>
Date: Tue, 11 Feb 2020 15:52:08 +0000
Subject: [PATCH 2/8] Cleanup test files

---
 .../math/prim/fun/quad_form_diag_test.cpp     | 33 +++++++++----------
 test/unit/math/prim/fun/quad_form_test.cpp    |  6 ++--
 2 files changed, 17 insertions(+), 22 deletions(-)

diff --git a/test/unit/math/prim/fun/quad_form_diag_test.cpp b/test/unit/math/prim/fun/quad_form_diag_test.cpp
index 63e981745f7..d7fd747a4c2 100644
--- a/test/unit/math/prim/fun/quad_form_diag_test.cpp
+++ b/test/unit/math/prim/fun/quad_form_diag_test.cpp
@@ -2,52 +2,49 @@
 #include <test/unit/math/prim/fun/expect_matrix_eq.hpp>
 #include <gtest/gtest.h>
 
-using Eigen::Dynamic;
-using Eigen::Matrix;
 using stan::math::quad_form_diag;
 
-TEST(MathMatrixPrimMat, quadFormDiag) {
-  Matrix<double, Dynamic, Dynamic> m(1, 1);
+TEST(MathMatrixPrim, quadFormDiag) {
+  Eigen::MatrixXd m(1, 1);
   m << 3;
 
-  Matrix<double, Dynamic, 1> v(1);
+  Eigen::VectorXd v(1);
   v << 9;
 
-  Matrix<double, Dynamic, Dynamic> v_m(1, 1);
-  v_m << 9;
+  Eigen::MatrixXd v_m = v.asDiagonal();
 
   expect_matrix_eq(v_m * m * v_m, quad_form_diag(m, v));
 }
-TEST(MathMatrixPrimMat, quadFormDiag2) {
-  Matrix<double, Dynamic, Dynamic> m(3, 3);
+
+TEST(MathMatrixPrim, quadFormDiag2) {
+  Eigen::MatrixXd m(3, 3);
   m << 1, 2, 3, 4, 5, 6, 7, 8, 9;
 
-  Matrix<double, Dynamic, 1> v(3);
+  Eigen::VectorXd v(3);
   v << 1, 2, 3;
 
-  Matrix<double, Dynamic, Dynamic> v_m(3, 3);
-  v_m << 1, 0, 0, 0, 2, 0, 0, 0, 3;
+  Eigen::MatrixXd v_m = v.asDiagonal();
 
   expect_matrix_eq(v_m * m * v_m, quad_form_diag(m, v));
 
-  Matrix<double, 1, Dynamic> rv(3);
+  Eigen::RowVectorXd rv(3);
   rv << 1, 2, 3;
   expect_matrix_eq(v_m * m * v_m, quad_form_diag(m, rv));
 }
 
-TEST(MathMatrixPrimMat, quadFormDiagException) {
-  Matrix<double, Dynamic, Dynamic> m(2, 2);
+TEST(MathMatrixPrim, quadFormDiagException) {
+  Eigen::MatrixXd m(2, 2);
   m << 2, 3, 4, 5;
   EXPECT_THROW(quad_form_diag(m, m), std::invalid_argument);
 
-  Matrix<double, Dynamic, 1> v(3);
+  Eigen::VectorXd v(3);
   v << 1, 2, 3;
   EXPECT_THROW(quad_form_diag(m, v), std::invalid_argument);
 
-  Matrix<double, Dynamic, Dynamic> m2(3, 2);
+  Eigen::MatrixXd m2(3, 2);
   m2 << 2, 3, 4, 5, 6, 7;
 
-  Matrix<double, Dynamic, 1> v2(2);
+  Eigen::VectorXd v2(2);
   v2 << 1, 2;
 
   EXPECT_THROW(quad_form_diag(m2, v), std::invalid_argument);
diff --git a/test/unit/math/prim/fun/quad_form_test.cpp b/test/unit/math/prim/fun/quad_form_test.cpp
index 09f82e036bf..92d0d53d94d 100644
--- a/test/unit/math/prim/fun/quad_form_test.cpp
+++ b/test/unit/math/prim/fun/quad_form_test.cpp
@@ -3,7 +3,6 @@
 
 TEST(MathMatrixPrim, quad_form_mat) {
   using stan::math::matrix_d;
-  using stan::math::quad_form;
 
   matrix_d ad(4, 4);
   matrix_d bd(4, 2);
@@ -13,7 +12,7 @@ TEST(MathMatrixPrim, quad_form_mat) {
       1.0, 112.0;
 
   // double-double
-  matrix_d resd = quad_form(ad, bd);
+  matrix_d resd = stan::math::quad_form(ad, bd);
   EXPECT_FLOAT_EQ(26033, resd(0, 0));
   EXPECT_FLOAT_EQ(3456, resd(0, 1));
   EXPECT_FLOAT_EQ(3396, resd(1, 0));
@@ -22,7 +21,6 @@ TEST(MathMatrixPrim, quad_form_mat) {
 
 TEST(MathMatrixPrim, quad_form_vec) {
   using stan::math::matrix_d;
-  using stan::math::quad_form;
   using stan::math::vector_d;
 
   matrix_d ad(4, 4);
@@ -34,6 +32,6 @@ TEST(MathMatrixPrim, quad_form_vec) {
       1.0, 112.0;
 
   // double-double
-  res = quad_form(ad, bd);
+  res = stan::math::quad_form(ad, bd);
   EXPECT_FLOAT_EQ(26033, res);
 }

From 8f3256f7a915a1e453ea5fb9f4b841b3ba89132c Mon Sep 17 00:00:00 2001
From: Marco Colombo <m.colombo@ed.ac.uk>
Date: Tue, 11 Feb 2020 15:58:12 +0000
Subject: [PATCH 3/8] Fix quad_form_sym mix tests

---
 test/unit/math/mix/fun/quad_form_sym_test.cpp | 35 +++++++++++++------
 1 file changed, 24 insertions(+), 11 deletions(-)

diff --git a/test/unit/math/mix/fun/quad_form_sym_test.cpp b/test/unit/math/mix/fun/quad_form_sym_test.cpp
index 03481980e03..3725880cb9d 100644
--- a/test/unit/math/mix/fun/quad_form_sym_test.cpp
+++ b/test/unit/math/mix/fun/quad_form_sym_test.cpp
@@ -1,8 +1,11 @@
 #include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/ad_tolerances.hpp>
 
 TEST(MathMixMatFun, quadFormSym) {
   auto f = [](const auto& x, const auto& y) {
-    return stan::math::quad_form_sym(x, y);
+    // symmetrize the input matrix
+    auto x_sym = ((x + x.transpose()) * 0.5).eval();
+    return stan::math::quad_form_sym(x_sym, y);
   };
 
   Eigen::MatrixXd a00;
@@ -13,11 +16,11 @@ TEST(MathMixMatFun, quadFormSym) {
   Eigen::MatrixXd a22(2, 2);
   a22 << 1, 2, 3, 4;
   Eigen::MatrixXd b22(2, 2);
-  a22 << -3, -2, -10, 112;
-  Eigen::MatrixXd a23(2, 3);
-  a23 << 1, 2, 3, 4, 5, 6;
-  Eigen::MatrixXd a42(4, 2);
-  a42 << 100, 10, 0, 1, -3, -3, 5, 2;
+  b22 << -3, -2, -10, 112;
+  Eigen::MatrixXd b23(2, 3);
+  b23 << 1, 2, 3, 4, 5, 6;
+  Eigen::MatrixXd b42(4, 2);
+  b42 << 100, 10, 0, 1, -3, -3, 5, 2;
   Eigen::MatrixXd a44(4, 4);
   a44 << 2, 3, 4, 5, 6, 10, 2, 2, 7, 2, 7, 1, 8, 2, 1, 112;
 
@@ -29,20 +32,30 @@ TEST(MathMixMatFun, quadFormSym) {
   Eigen::VectorXd v4(4);
   v4 << 100, 0, -3, 5;
 
+  stan::test::ad_tolerances tols;
+  tols.hessian_hessian_ = 2e-1;
+  tols.hessian_fvar_hessian_ = 2e-1;
+
   stan::test::expect_ad(f, a00, a00);
   stan::test::expect_ad(f, a11, b11);
-  stan::test::expect_ad(f, a22, b22);
-  stan::test::expect_ad(f, a44, a42);
+  stan::test::expect_ad(tols, f, a22, b22);
+  stan::test::expect_ad(f, a22, b23);
+  stan::test::expect_ad(tols, f, a44, b42);
 
   stan::test::expect_ad(f, a00, v0);
   stan::test::expect_ad(f, a11, v1);
   stan::test::expect_ad(f, a22, v2);
-  stan::test::expect_ad(f, a44, v4);
+  stan::test::expect_ad(tols, f, a44, v4);
+
+  // asymmetric case should throw
+
+  auto g = [](const auto& x, const auto& y) {
+    return stan::math::quad_form_sym(x, y);
+  };
 
-  // asymmetric case should thorw
   Eigen::MatrixXd u(4, 4);
   u << 2, 3, 4, 5, 6, 10, 2, 2, 7, 2, 7, 1, 8, 2, 1, 112;
   Eigen::MatrixXd v(4, 2);
   v << 100, 10, 0, 1, -3, -3, 5, 2;
-  stan::test::expect_ad(f, u, v);
+  stan::test::expect_ad(g, u, v);
 }

From 8b4797809d5b03e0126c1110023bd62dc43e1cee Mon Sep 17 00:00:00 2001
From: Marco Colombo <m.colombo@ed.ac.uk>
Date: Tue, 11 Feb 2020 16:00:50 +0000
Subject: [PATCH 4/8] Use require_var_t to simplify code

---
 stan/math/rev/fun/quad_form.hpp | 19 ++++++++-----------
 1 file changed, 8 insertions(+), 11 deletions(-)

diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp
index d26edabc2bb..eae0d8a9772 100644
--- a/stan/math/rev/fun/quad_form.hpp
+++ b/stan/math/rev/fun/quad_form.hpp
@@ -95,12 +95,10 @@ class quad_form_vari : public vari {
 };
 }  // namespace internal
 
-template <typename Ta, int Ra, int Ca, typename Tb, int Rb, int Cb>
-inline typename std::enable_if<std::is_same<Ta, var>::value
-                                   || std::is_same<Tb, var>::value,
-                               Eigen::Matrix<var, Cb, Cb> >::type
-quad_form(const Eigen::Matrix<Ta, Ra, Ca>& A,
-          const Eigen::Matrix<Tb, Rb, Cb>& B) {
+template <typename Ta, int Ra, int Ca, typename Tb, int Rb, int Cb,
+          require_any_var_t<Ta, Tb>...>
+inline Eigen::Matrix<var, Cb, Cb> quad_form(
+    const Eigen::Matrix<Ta, Ra, Ca>& A, const Eigen::Matrix<Tb, Rb, Cb>& B) {
   check_square("quad_form", "A", A);
   check_multiplicable("quad_form", "A", A, "B", B);
 
@@ -110,11 +108,10 @@ quad_form(const Eigen::Matrix<Ta, Ra, Ca>& A,
   return baseVari->impl_->C_;
 }
 
-template <typename Ta, int Ra, int Ca, typename Tb, int Rb>
-inline typename std::enable_if<
-    std::is_same<Ta, var>::value || std::is_same<Tb, var>::value, var>::type
-quad_form(const Eigen::Matrix<Ta, Ra, Ca>& A,
-          const Eigen::Matrix<Tb, Rb, 1>& B) {
+template <typename Ta, int Ra, int Ca, typename Tb, int Rb,
+          require_any_var_t<Ta, Tb>...>
+inline var quad_form(const Eigen::Matrix<Ta, Ra, Ca>& A,
+                     const Eigen::Matrix<Tb, Rb, 1>& B) {
   check_square("quad_form", "A", A);
   check_multiplicable("quad_form", "A", A, "B", B);
 

From 5b5183b96c31447437d91856d6273453b4fc6f04 Mon Sep 17 00:00:00 2001
From: Marco Colombo <m.colombo@ed.ac.uk>
Date: Tue, 11 Feb 2020 16:34:08 +0000
Subject: [PATCH 5/8] Add fwd version of quad_form and mix tests

---
 stan/math/fwd/fun.hpp                     |  1 +
 stan/math/fwd/fun/quad_form.hpp           | 64 +++++++++++++++++++++++
 test/unit/math/mix/fun/quad_form_test.cpp | 48 +++++++++++++++++
 3 files changed, 113 insertions(+)
 create mode 100644 stan/math/fwd/fun/quad_form.hpp
 create mode 100644 test/unit/math/mix/fun/quad_form_test.cpp

diff --git a/stan/math/fwd/fun.hpp b/stan/math/fwd/fun.hpp
index f53bb580d55..3b98105fec5 100644
--- a/stan/math/fwd/fun.hpp
+++ b/stan/math/fwd/fun.hpp
@@ -92,6 +92,7 @@
 #include <stan/math/fwd/fun/primitive_value.hpp>
 #include <stan/math/fwd/fun/qr_Q.hpp>
 #include <stan/math/fwd/fun/qr_R.hpp>
+#include <stan/math/fwd/fun/quad_form.hpp>
 #include <stan/math/fwd/fun/quad_form_sym.hpp>
 #include <stan/math/fwd/fun/rising_factorial.hpp>
 #include <stan/math/fwd/fun/round.hpp>
diff --git a/stan/math/fwd/fun/quad_form.hpp b/stan/math/fwd/fun/quad_form.hpp
new file mode 100644
index 00000000000..2e8a4c4af75
--- /dev/null
+++ b/stan/math/fwd/fun/quad_form.hpp
@@ -0,0 +1,64 @@
+#ifndef STAN_MATH_FWD_FUN_QUAD_FORM_HPP
+#define STAN_MATH_FWD_FUN_QUAD_FORM_HPP
+
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/fwd/fun/dot_product.hpp>
+#include <stan/math/fwd/fun/multiply.hpp>
+
+namespace stan {
+namespace math {
+
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam Ta type of elements in the square matrix
+ * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
+ * @tparam Ca number of columns in the square matrix, can be Eigen::Dynamic
+ * @tparam Tb type of elements in the second matrix
+ * @tparam Rb number of rows in the second matrix, can be Eigen::Dynamic
+ * @tparam Cb number of columns in the second matrix, can be Eigen::Dynamic
+ *
+ * @param A square matrix
+ * @param B second matrix
+ * @return The quadratic form B^T A B, which is a symmetric matrix of size Cb
+ * (although symmetry is not guaranteed due to numerical precision).
+ * @throws std::invalid_argument if A is not square, or if A cannot be
+ * multiplied by B
+ */
+template <typename Ta, int Ra, int Ca, typename Tb, int Rb, int Cb,
+          require_any_fvar_t<Ta, Tb>...>
+inline Eigen::Matrix<return_type_t<Ta, Tb>, Cb, Cb> quad_form(
+    const Eigen::Matrix<Ta, Ra, Ca>& A, const Eigen::Matrix<Tb, Rb, Cb>& B) {
+  check_square("quad_form", "A", A);
+  check_multiplicable("quad_form", "A", A, "B", B);
+  return multiply(transpose(B), multiply(A, B));
+}
+
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam Ta type of elements in the square matrix
+ * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
+ * @tparam Ca number of columns in the square matrix, can be Eigen::Dynamic
+ * @tparam Tb type of elements in the vector
+ * @tparam Rb number of rows in the vector, can be Eigen::Dynamic
+ *
+ * @param A square matrix
+ * @param B vector
+ * @return The quadratic form B^T A B, which is a scalar.
+ * @throws std::invalid_argument if A is not square, or if A cannot be
+ * multiplied by B
+ */
+template <typename Ta, int Ra, int Ca, typename Tb, int Rb,
+          require_any_fvar_t<Ta, Tb>...>
+inline return_type_t<Ta, Tb> quad_form(const Eigen::Matrix<Ta, Ra, Ca>& A,
+                                       const Eigen::Matrix<Tb, Rb, 1>& B) {
+  check_square("quad_form", "A", A);
+  check_multiplicable("quad_form", "A", A, "B", B);
+  return dot_product(B, multiply(A, B));
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/test/unit/math/mix/fun/quad_form_test.cpp b/test/unit/math/mix/fun/quad_form_test.cpp
new file mode 100644
index 00000000000..02d78b553c8
--- /dev/null
+++ b/test/unit/math/mix/fun/quad_form_test.cpp
@@ -0,0 +1,48 @@
+#include <test/unit/math/test_ad.hpp>
+#include <test/unit/math/ad_tolerances.hpp>
+
+TEST(MathMixMatFun, quadForm) {
+  using stan::test::relative_tolerance;
+  auto f = [](const auto& x, const auto& y) {
+    return stan::math::quad_form(x, y);
+  };
+
+  Eigen::MatrixXd a00;
+  Eigen::MatrixXd a11(1, 1);
+  a11 << 1;
+  Eigen::MatrixXd b11(1, 1);
+  b11 << -2;
+  Eigen::MatrixXd a22(2, 2);
+  a22 << 1, 2, 3, 4;
+  Eigen::MatrixXd b22(2, 2);
+  b22 << -3, -2, -10, 112;
+  Eigen::MatrixXd b23(2, 3);
+  b23 << 1, 2, 3, 4, 5, 6;
+  Eigen::MatrixXd b42(4, 2);
+  b42 << 100, 10, 0, 1, -3, -3, 5, 2;
+  Eigen::MatrixXd a44(4, 4);
+  a44 << 2, 3, 4, 5, 6, 10, 2, 2, 7, 2, 7, 1, 8, 2, 1, 112;
+
+  Eigen::VectorXd v0(0);
+  Eigen::VectorXd v1(1);
+  v1 << 42;
+  Eigen::VectorXd v2(2);
+  v2 << -3, 13;
+  Eigen::VectorXd v4(4);
+  v4 << 100, 0, -3, 5;
+
+  stan::test::ad_tolerances tols;
+  tols.hessian_hessian_ = 2e-1;
+  tols.hessian_fvar_hessian_ = 2e-1;
+
+  stan::test::expect_ad(f, a00, a00);
+  stan::test::expect_ad(f, a11, b11);
+  stan::test::expect_ad(tols, f, a22, b22);
+  stan::test::expect_ad(f, a22, b23);
+  stan::test::expect_ad(tols, f, a44, b42);
+
+  stan::test::expect_ad(f, a00, v0);
+  stan::test::expect_ad(f, a11, v1);
+  stan::test::expect_ad(f, a22, v2);
+  stan::test::expect_ad(tols, f, a44, v4);
+}

From 065a1cde11c3795c7eebaf5cde26be55f063f473 Mon Sep 17 00:00:00 2001
From: Marco Colombo <m.colombo@ed.ac.uk>
Date: Tue, 11 Feb 2020 16:46:37 +0000
Subject: [PATCH 6/8] Consolidate implementations and add docs

---
 stan/math/fwd/fun/quad_form_sym.hpp  | 66 ++++++++++++++++++----------
 stan/math/prim/fun/quad_form.hpp     | 31 ++++++++++++-
 stan/math/prim/fun/quad_form_sym.hpp | 30 +++++++++++++
 stan/math/rev/fun/quad_form.hpp      | 32 ++++++++++++++
 stan/math/rev/fun/quad_form_sym.hpp  | 32 ++++++++++++++
 5 files changed, 165 insertions(+), 26 deletions(-)

diff --git a/stan/math/fwd/fun/quad_form_sym.hpp b/stan/math/fwd/fun/quad_form_sym.hpp
index c5f069b2e56..536b9cc20c6 100644
--- a/stan/math/fwd/fun/quad_form_sym.hpp
+++ b/stan/math/fwd/fun/quad_form_sym.hpp
@@ -8,40 +8,58 @@
 namespace stan {
 namespace math {
 
-template <int RA, int CA, int RB, int CB, typename T>
-inline Eigen::Matrix<fvar<T>, CB, CB> quad_form_sym(
-    const Eigen::Matrix<fvar<T>, RA, CA>& A,
-    const Eigen::Matrix<double, RB, CB>& B) {
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam TA type of elements in the symmetric matrix
+ * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam CA number of columns in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam TB type of elements in the second matrix
+ * @tparam RB number of rows in the second matrix, can be Eigen::Dynamic
+ * @tparam CB number of columns in the second matrix, can be Eigen::Dynamic
+ *
+ * @param A symmetric matrix
+ * @param B second matrix
+ * @return The quadratic form B^T A B, which is guaranteed to be a symmetric
+ * matrix of size CB.
+ * @throws std::invalid_argument if A is not symmetric, or if A cannot be
+ * multiplied by B
+ */
+template <typename TA, int RA, int CA, typename TB, int RB, int CB,
+          require_any_fvar_t<TA, TB>...>
+inline Eigen::Matrix<return_type_t<TA, TB>, CB, CB> quad_form_sym(
+    const Eigen::Matrix<TA, RA, CA>& A, const Eigen::Matrix<TB, RB, CB>& B) {
+  using T = return_type_t<TA, TB>;
   check_multiplicable("quad_form_sym", "A", A, "B", B);
   check_symmetric("quad_form_sym", "A", A);
-  Eigen::Matrix<fvar<T>, CB, CB> ret(multiply(transpose(B), multiply(A, B)));
+  Eigen::Matrix<T, CB, CB> ret(multiply(transpose(B), multiply(A, B)));
   return T(0.5) * (ret + transpose(ret));
 }
 
-template <int RA, int CA, int RB, typename T>
-inline fvar<T> quad_form_sym(const Eigen::Matrix<fvar<T>, RA, CA>& A,
-                             const Eigen::Matrix<double, RB, 1>& B) {
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam TA type of elements in the symmetric matrix
+ * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam CA number of columns in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam TB type of elements in the vector
+ * @tparam RB number of rows in the vector, can be Eigen::Dynamic
+ *
+ * @param A symmetric matrix
+ * @param B vector
+ * @return The quadratic form B^T A B, which is a scalar.
+ * @throws std::invalid_argument if A is not symmetric, or if A cannot be
+ * multiplied by B
+ */
+template <typename TA, int RA, int CA, typename TB, int RB,
+          require_any_fvar_t<TA, TB>...>
+inline return_type_t<TA, TB> quad_form_sym(const Eigen::Matrix<TA, RA, CA>& A,
+                                           const Eigen::Matrix<TB, RB, 1>& B) {
   check_multiplicable("quad_form_sym", "A", A, "B", B);
   check_symmetric("quad_form_sym", "A", A);
   return dot_product(B, multiply(A, B));
 }
-template <int RA, int CA, int RB, int CB, typename T>
-inline Eigen::Matrix<fvar<T>, CB, CB> quad_form_sym(
-    const Eigen::Matrix<double, RA, CA>& A,
-    const Eigen::Matrix<fvar<T>, RB, CB>& B) {
-  check_multiplicable("quad_form_sym", "A", A, "B", B);
-  check_symmetric("quad_form_sym", "A", A);
-  Eigen::Matrix<fvar<T>, CB, CB> ret(multiply(transpose(B), multiply(A, B)));
-  return T(0.5) * (ret + transpose(ret));
-}
 
-template <int RA, int CA, int RB, typename T>
-inline fvar<T> quad_form_sym(const Eigen::Matrix<double, RA, CA>& A,
-                             const Eigen::Matrix<fvar<T>, RB, 1>& B) {
-  check_multiplicable("quad_form_sym", "A", A, "B", B);
-  check_symmetric("quad_form_sym", "A", A);
-  return dot_product(B, multiply(A, B));
-}
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/fun/quad_form.hpp b/stan/math/prim/fun/quad_form.hpp
index 1f35a99f2c4..d15075451b4 100644
--- a/stan/math/prim/fun/quad_form.hpp
+++ b/stan/math/prim/fun/quad_form.hpp
@@ -8,8 +8,21 @@ namespace stan {
 namespace math {
 
 /**
- * Compute B^T A B
- **/
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam RA number of rows in the square matrix, can be Eigen::Dynamic
+ * @tparam CA number of columns in the square matrix, can be Eigen::Dynamic
+ * @tparam RB number of rows in the second matrix, can be Eigen::Dynamic
+ * @tparam CB number of columns in the second matrix, can be Eigen::Dynamic
+ * @tparam T type of elements
+ *
+ * @param A square matrix
+ * @param B second matrix
+ * @return The quadratic form B^T A B, which is a symmetric matrix of size CB
+ * (although symmetry is not guaranteed due to numerical precision).
+ * @throws std::invalid_argument if A is not square, or if A cannot be
+ * multiplied by B
+ */
 template <int RA, int CA, int RB, int CB, typename T>
 inline Eigen::Matrix<T, CB, CB> quad_form(const Eigen::Matrix<T, RA, CA>& A,
                                           const Eigen::Matrix<T, RB, CB>& B) {
@@ -18,6 +31,20 @@ inline Eigen::Matrix<T, CB, CB> quad_form(const Eigen::Matrix<T, RA, CA>& A,
   return B.transpose() * A * B;
 }
 
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam RA number of rows in the square matrix, can be Eigen::Dynamic
+ * @tparam CA number of columns in the square matrix, can be Eigen::Dynamic
+ * @tparam RB number of rows in the vector, can be Eigen::Dynamic
+ * @tparam T type of elements
+ *
+ * @param A square matrix
+ * @param B vector
+ * @return The quadratic form B^T A B, which is a scalar.
+ * @throws std::invalid_argument if A is not square, or if A cannot be
+ * multiplied by B
+ */
 template <int RA, int CA, int RB, typename T>
 inline T quad_form(const Eigen::Matrix<T, RA, CA>& A,
                    const Eigen::Matrix<T, RB, 1>& B) {
diff --git a/stan/math/prim/fun/quad_form_sym.hpp b/stan/math/prim/fun/quad_form_sym.hpp
index 391ebe3ef4c..4e022365001 100644
--- a/stan/math/prim/fun/quad_form_sym.hpp
+++ b/stan/math/prim/fun/quad_form_sym.hpp
@@ -7,6 +7,22 @@
 namespace stan {
 namespace math {
 
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam CA number of columns in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam RB number of rows in the second matrix, can be Eigen::Dynamic
+ * @tparam CB number of columns in the second matrix, can be Eigen::Dynamic
+ * @tparam TA type of elements
+ *
+ * @param A symmetric matrix
+ * @param B second matrix
+ * @return The quadratic form B^T A B, which is guaranteed to be a symmetric
+ * matrix of size CB.
+ * @throws std::invalid_argument if A is not symmetric, or if A cannot be
+ * multiplied by B
+ */
 template <int RA, int CA, int RB, int CB, typename T>
 inline Eigen::Matrix<T, CB, CB> quad_form_sym(
     const Eigen::Matrix<T, RA, CA>& A, const Eigen::Matrix<T, RB, CB>& B) {
@@ -16,6 +32,20 @@ inline Eigen::Matrix<T, CB, CB> quad_form_sym(
   return T(0.5) * (ret + ret.transpose());
 }
 
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam CA number of columns in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam RB number of rows in the vector, can be Eigen::Dynamic
+ * @tparam T type of elements
+ *
+ * @param A symmetric matrix
+ * @param B vector
+ * @return The quadratic form B^T A B, which is a scalar.
+ * @throws std::invalid_argument if A is not symmetric, or if A cannot be
+ * multiplied by B
+ */
 template <int RA, int CA, int RB, typename T>
 inline T quad_form_sym(const Eigen::Matrix<T, RA, CA>& A,
                        const Eigen::Matrix<T, RB, 1>& B) {
diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp
index eae0d8a9772..78241ac63a2 100644
--- a/stan/math/rev/fun/quad_form.hpp
+++ b/stan/math/rev/fun/quad_form.hpp
@@ -95,6 +95,23 @@ class quad_form_vari : public vari {
 };
 }  // namespace internal
 
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam Ta type of elements in the square matrix
+ * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
+ * @tparam Ca number of columns in the square matrix, can be Eigen::Dynamic
+ * @tparam Tb type of elements in the second matrix
+ * @tparam Rb number of rows in the second matrix, can be Eigen::Dynamic
+ * @tparam Cb number of columns in the second matrix, can be Eigen::Dynamic
+ *
+ * @param A square matrix
+ * @param B second matrix
+ * @return The quadratic form B^T A B, which is a symmetric matrix of size Cb
+ * (although symmetry is not guaranteed due to numerical precision).
+ * @throws std::invalid_argument if A is not square, or if A cannot be
+ * multiplied by B
+ */
 template <typename Ta, int Ra, int Ca, typename Tb, int Rb, int Cb,
           require_any_var_t<Ta, Tb>...>
 inline Eigen::Matrix<var, Cb, Cb> quad_form(
@@ -108,6 +125,21 @@ inline Eigen::Matrix<var, Cb, Cb> quad_form(
   return baseVari->impl_->C_;
 }
 
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam Ta type of elements in the square matrix
+ * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
+ * @tparam Ca number of columns in the square matrix, can be Eigen::Dynamic
+ * @tparam Tb type of elements in the vector
+ * @tparam Rb number of rows in the vector, can be Eigen::Dynamic
+ *
+ * @param A square matrix
+ * @param B vector
+ * @return The quadratic form B^T A B, which is a scalar.
+ * @throws std::invalid_argument if A is not square, or if A cannot be
+ * multiplied by B
+ */
 template <typename Ta, int Ra, int Ca, typename Tb, int Rb,
           require_any_var_t<Ta, Tb>...>
 inline var quad_form(const Eigen::Matrix<Ta, Ra, Ca>& A,
diff --git a/stan/math/rev/fun/quad_form_sym.hpp b/stan/math/rev/fun/quad_form_sym.hpp
index c81501efcf9..1a42500bc94 100644
--- a/stan/math/rev/fun/quad_form_sym.hpp
+++ b/stan/math/rev/fun/quad_form_sym.hpp
@@ -11,6 +11,23 @@
 namespace stan {
 namespace math {
 
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam Ta type of elements in the symmetric matrix
+ * @tparam Ra number of rows in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam Ca number of columns in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam Tb type of elements in the second matrix
+ * @tparam Rb number of rows in the second matrix, can be Eigen::Dynamic
+ * @tparam Cb number of columns in the second matrix, can be Eigen::Dynamic
+ *
+ * @param A symmetric matrix
+ * @param B second matrix
+ * @return The quadratic form B^T A B, which is guaranteed to be a symmetric
+ * matrix of size Cb.
+ * @throws std::invalid_argument if A is not symmetric, or if A cannot be
+ * multiplied by B
+ */
 template <typename Ta, int Ra, int Ca, typename Tb, int Rb, int Cb,
           require_any_var_t<Ta, Tb>...>
 inline Eigen::Matrix<var, Cb, Cb> quad_form_sym(
@@ -24,6 +41,21 @@ inline Eigen::Matrix<var, Cb, Cb> quad_form_sym(
   return baseVari->impl_->C_;
 }
 
+/**
+ * Compute the quadratic form B^T A B.
+ *
+ * @tparam Ta type of elements in the symmetric matrix
+ * @tparam Ra number of rows in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam Ca number of columns in the symmetric matrix, can be Eigen::Dynamic
+ * @tparam Tb type of elements in the vector
+ * @tparam Rb number of rows in the vector, can be Eigen::Dynamic
+ *
+ * @param A symmetric matrix
+ * @param B vector
+ * @return The quadratic form B^T A B, which is a scalar.
+ * @throws std::invalid_argument if A is not symmetric, or if A cannot be
+ * multiplied by B
+ */
 template <typename Ta, int Ra, int Ca, typename Tb, int Rb,
           require_any_var_t<Ta, Tb>...>
 inline var quad_form_sym(const Eigen::Matrix<Ta, Ra, Ca>& A,

From 702acb1ad1e428f82a7f69ad48b9c2e71fe6b86b Mon Sep 17 00:00:00 2001
From: Marco Colombo <m.colombo@ed.ac.uk>
Date: Wed, 12 Feb 2020 11:07:07 +0000
Subject: [PATCH 7/8] Improve documentation

---
 stan/math/fwd/fun/quad_form.hpp      | 12 +++++++-----
 stan/math/fwd/fun/quad_form_sym.hpp  | 11 ++++++-----
 stan/math/prim/fun/quad_form.hpp     | 12 +++++++-----
 stan/math/prim/fun/quad_form_sym.hpp | 11 ++++++-----
 stan/math/rev/fun/quad_form.hpp      | 12 +++++++-----
 stan/math/rev/fun/quad_form_sym.hpp  | 11 ++++++-----
 6 files changed, 39 insertions(+), 30 deletions(-)

diff --git a/stan/math/fwd/fun/quad_form.hpp b/stan/math/fwd/fun/quad_form.hpp
index 2e8a4c4af75..f4d4adc50f3 100644
--- a/stan/math/fwd/fun/quad_form.hpp
+++ b/stan/math/fwd/fun/quad_form.hpp
@@ -9,7 +9,10 @@ namespace stan {
 namespace math {
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$.
+ *
+ * Symmetry of the resulting matrix is not guaranteed due to numerical
+ * precision.
  *
  * @tparam Ta type of elements in the square matrix
  * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
@@ -20,8 +23,7 @@ namespace math {
  *
  * @param A square matrix
  * @param B second matrix
- * @return The quadratic form B^T A B, which is a symmetric matrix of size Cb
- * (although symmetry is not guaranteed due to numerical precision).
+ * @return The quadratic form, which is a symmetric matrix of size Cb.
  * @throws std::invalid_argument if A is not square, or if A cannot be
  * multiplied by B
  */
@@ -35,7 +37,7 @@ inline Eigen::Matrix<return_type_t<Ta, Tb>, Cb, Cb> quad_form(
 }
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$.
  *
  * @tparam Ta type of elements in the square matrix
  * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
@@ -45,7 +47,7 @@ inline Eigen::Matrix<return_type_t<Ta, Tb>, Cb, Cb> quad_form(
  *
  * @param A square matrix
  * @param B vector
- * @return The quadratic form B^T A B, which is a scalar.
+ * @return The quadratic form (a scalar).
  * @throws std::invalid_argument if A is not square, or if A cannot be
  * multiplied by B
  */
diff --git a/stan/math/fwd/fun/quad_form_sym.hpp b/stan/math/fwd/fun/quad_form_sym.hpp
index 536b9cc20c6..a77e2ea65e7 100644
--- a/stan/math/fwd/fun/quad_form_sym.hpp
+++ b/stan/math/fwd/fun/quad_form_sym.hpp
@@ -9,7 +9,9 @@ namespace stan {
 namespace math {
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$ of a symmetric matrix.
+ *
+ * Symmetry of the resulting matrix is guaranteed.
  *
  * @tparam TA type of elements in the symmetric matrix
  * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
@@ -20,8 +22,7 @@ namespace math {
  *
  * @param A symmetric matrix
  * @param B second matrix
- * @return The quadratic form B^T A B, which is guaranteed to be a symmetric
- * matrix of size CB.
+ * @return The quadratic form, which is a symmetric matrix of size CB.
  * @throws std::invalid_argument if A is not symmetric, or if A cannot be
  * multiplied by B
  */
@@ -37,7 +38,7 @@ inline Eigen::Matrix<return_type_t<TA, TB>, CB, CB> quad_form_sym(
 }
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$ of a symmetric matrix.
  *
  * @tparam TA type of elements in the symmetric matrix
  * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
@@ -47,7 +48,7 @@ inline Eigen::Matrix<return_type_t<TA, TB>, CB, CB> quad_form_sym(
  *
  * @param A symmetric matrix
  * @param B vector
- * @return The quadratic form B^T A B, which is a scalar.
+ * @return The quadratic form (a scalar).
  * @throws std::invalid_argument if A is not symmetric, or if A cannot be
  * multiplied by B
  */
diff --git a/stan/math/prim/fun/quad_form.hpp b/stan/math/prim/fun/quad_form.hpp
index d15075451b4..34b1c4ba806 100644
--- a/stan/math/prim/fun/quad_form.hpp
+++ b/stan/math/prim/fun/quad_form.hpp
@@ -8,7 +8,10 @@ namespace stan {
 namespace math {
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$.
+ *
+ * Symmetry of the resulting matrix is not guaranteed due to numerical
+ * precision.
  *
  * @tparam RA number of rows in the square matrix, can be Eigen::Dynamic
  * @tparam CA number of columns in the square matrix, can be Eigen::Dynamic
@@ -18,8 +21,7 @@ namespace math {
  *
  * @param A square matrix
  * @param B second matrix
- * @return The quadratic form B^T A B, which is a symmetric matrix of size CB
- * (although symmetry is not guaranteed due to numerical precision).
+ * @return The quadratic form, which is a symmetric matrix of size CB.
  * @throws std::invalid_argument if A is not square, or if A cannot be
  * multiplied by B
  */
@@ -32,7 +34,7 @@ inline Eigen::Matrix<T, CB, CB> quad_form(const Eigen::Matrix<T, RA, CA>& A,
 }
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$.
  *
  * @tparam RA number of rows in the square matrix, can be Eigen::Dynamic
  * @tparam CA number of columns in the square matrix, can be Eigen::Dynamic
@@ -41,7 +43,7 @@ inline Eigen::Matrix<T, CB, CB> quad_form(const Eigen::Matrix<T, RA, CA>& A,
  *
  * @param A square matrix
  * @param B vector
- * @return The quadratic form B^T A B, which is a scalar.
+ * @return The quadratic form (a scalar).
  * @throws std::invalid_argument if A is not square, or if A cannot be
  * multiplied by B
  */
diff --git a/stan/math/prim/fun/quad_form_sym.hpp b/stan/math/prim/fun/quad_form_sym.hpp
index 4e022365001..a67ba593780 100644
--- a/stan/math/prim/fun/quad_form_sym.hpp
+++ b/stan/math/prim/fun/quad_form_sym.hpp
@@ -8,7 +8,9 @@ namespace stan {
 namespace math {
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$ of a symmetric matrix.
+ *
+ * Symmetry of the resulting matrix is guaranteed.
  *
  * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
  * @tparam CA number of columns in the symmetric matrix, can be Eigen::Dynamic
@@ -18,8 +20,7 @@ namespace math {
  *
  * @param A symmetric matrix
  * @param B second matrix
- * @return The quadratic form B^T A B, which is guaranteed to be a symmetric
- * matrix of size CB.
+ * @return The quadratic form, which is a symmetric matrix of size CB.
  * @throws std::invalid_argument if A is not symmetric, or if A cannot be
  * multiplied by B
  */
@@ -33,7 +34,7 @@ inline Eigen::Matrix<T, CB, CB> quad_form_sym(
 }
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$ of a symmetric matrix.
  *
  * @tparam RA number of rows in the symmetric matrix, can be Eigen::Dynamic
  * @tparam CA number of columns in the symmetric matrix, can be Eigen::Dynamic
@@ -42,7 +43,7 @@ inline Eigen::Matrix<T, CB, CB> quad_form_sym(
  *
  * @param A symmetric matrix
  * @param B vector
- * @return The quadratic form B^T A B, which is a scalar.
+ * @return The quadratic form (a scalar).
  * @throws std::invalid_argument if A is not symmetric, or if A cannot be
  * multiplied by B
  */
diff --git a/stan/math/rev/fun/quad_form.hpp b/stan/math/rev/fun/quad_form.hpp
index 78241ac63a2..406f6398180 100644
--- a/stan/math/rev/fun/quad_form.hpp
+++ b/stan/math/rev/fun/quad_form.hpp
@@ -96,7 +96,10 @@ class quad_form_vari : public vari {
 }  // namespace internal
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$.
+ *
+ * Symmetry of the resulting matrix is not guaranteed due to numerical
+ * precision.
  *
  * @tparam Ta type of elements in the square matrix
  * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
@@ -107,8 +110,7 @@ class quad_form_vari : public vari {
  *
  * @param A square matrix
  * @param B second matrix
- * @return The quadratic form B^T A B, which is a symmetric matrix of size Cb
- * (although symmetry is not guaranteed due to numerical precision).
+ * @return The quadratic form, which is a symmetric matrix of size Cb.
  * @throws std::invalid_argument if A is not square, or if A cannot be
  * multiplied by B
  */
@@ -126,7 +128,7 @@ inline Eigen::Matrix<var, Cb, Cb> quad_form(
 }
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$.
  *
  * @tparam Ta type of elements in the square matrix
  * @tparam Ra number of rows in the square matrix, can be Eigen::Dynamic
@@ -136,7 +138,7 @@ inline Eigen::Matrix<var, Cb, Cb> quad_form(
  *
  * @param A square matrix
  * @param B vector
- * @return The quadratic form B^T A B, which is a scalar.
+ * @return The quadratic form (a scalar).
  * @throws std::invalid_argument if A is not square, or if A cannot be
  * multiplied by B
  */
diff --git a/stan/math/rev/fun/quad_form_sym.hpp b/stan/math/rev/fun/quad_form_sym.hpp
index 1a42500bc94..2254b0c323b 100644
--- a/stan/math/rev/fun/quad_form_sym.hpp
+++ b/stan/math/rev/fun/quad_form_sym.hpp
@@ -12,7 +12,9 @@ namespace stan {
 namespace math {
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$ of a symmetric matrix.
+ *
+ * Symmetry of the resulting matrix is guaranteed.
  *
  * @tparam Ta type of elements in the symmetric matrix
  * @tparam Ra number of rows in the symmetric matrix, can be Eigen::Dynamic
@@ -23,8 +25,7 @@ namespace math {
  *
  * @param A symmetric matrix
  * @param B second matrix
- * @return The quadratic form B^T A B, which is guaranteed to be a symmetric
- * matrix of size Cb.
+ * @return The quadratic form, which is a symmetric matrix of size Cb.
  * @throws std::invalid_argument if A is not symmetric, or if A cannot be
  * multiplied by B
  */
@@ -42,7 +43,7 @@ inline Eigen::Matrix<var, Cb, Cb> quad_form_sym(
 }
 
 /**
- * Compute the quadratic form B^T A B.
+ * Return the quadratic form \f$ B^T A B \f$ of a symmetric matrix.
  *
  * @tparam Ta type of elements in the symmetric matrix
  * @tparam Ra number of rows in the symmetric matrix, can be Eigen::Dynamic
@@ -52,7 +53,7 @@ inline Eigen::Matrix<var, Cb, Cb> quad_form_sym(
  *
  * @param A symmetric matrix
  * @param B vector
- * @return The quadratic form B^T A B, which is a scalar.
+ * @return The quadratic form (a scalar).
  * @throws std::invalid_argument if A is not symmetric, or if A cannot be
  * multiplied by B
  */

From 23eded070ca8dacda30e57feb1ec9c475308bdca Mon Sep 17 00:00:00 2001
From: Marco Colombo <m.colombo@ed.ac.uk>
Date: Wed, 12 Feb 2020 12:30:19 +0000
Subject: [PATCH 8/8] Add more tests for boundary conditions

---
 .../unit/math/prim/fun/quad_form_sym_test.cpp | 66 +++++++++++--------
 test/unit/math/prim/fun/quad_form_test.cpp    | 33 +++++++++-
 2 files changed, 68 insertions(+), 31 deletions(-)

diff --git a/test/unit/math/prim/fun/quad_form_sym_test.cpp b/test/unit/math/prim/fun/quad_form_sym_test.cpp
index 528dbba97b7..a1b450e8fe7 100644
--- a/test/unit/math/prim/fun/quad_form_sym_test.cpp
+++ b/test/unit/math/prim/fun/quad_form_sym_test.cpp
@@ -3,6 +3,24 @@
 
 TEST(MathMatrixPrim, quad_form_sym_mat) {
   using stan::math::matrix_d;
+  using stan::math::quad_form_sym;
+
+  matrix_d resd;
+  matrix_d m0;
+  EXPECT_THROW(quad_form_sym(m0, m0), std::invalid_argument);
+
+  matrix_d m1(1, 1);
+  m1 << 2;
+  resd = quad_form_sym(m1, m1);
+  EXPECT_FLOAT_EQ(8, resd(0));
+
+  matrix_d m2(1, 2);
+  m2 << 1, 2;
+  resd = quad_form_sym(m1, m2);
+  EXPECT_FLOAT_EQ(2, resd(0, 0));
+  EXPECT_FLOAT_EQ(4, resd(0, 1));
+  EXPECT_FLOAT_EQ(4, resd(1, 0));
+  EXPECT_FLOAT_EQ(8, resd(1, 1));
 
   matrix_d ad(4, 4);
   matrix_d bd(4, 2);
@@ -11,18 +29,35 @@ TEST(MathMatrixPrim, quad_form_sym_mat) {
   ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
       1.0, 112.0;
 
-  // double-double
-  matrix_d resd = stan::math::quad_form_sym(ad, bd);
+  resd = quad_form_sym(ad, bd);
   EXPECT_FLOAT_EQ(25433, resd(0, 0));
   EXPECT_FLOAT_EQ(3396, resd(0, 1));
   EXPECT_FLOAT_EQ(3396, resd(1, 0));
   EXPECT_FLOAT_EQ(725, resd(1, 1));
+
+  bd.resize(4, 3);
+  bd << 100, 10, 11, 0, 1, 12, -3, -3, 34, 5, 2, 44;
+  resd = quad_form_sym(ad, bd);
+  EXPECT_EQ(resd(1, 0), resd(0, 1));
+  EXPECT_EQ(resd(2, 0), resd(0, 2));
+  EXPECT_EQ(resd(2, 1), resd(1, 2));
 }
 
 TEST(MathMatrixPrim, quad_form_sym_vec) {
   using stan::math::matrix_d;
+  using stan::math::quad_form_sym;
   using stan::math::vector_d;
 
+  matrix_d m0;
+  vector_d v0;
+  EXPECT_THROW(quad_form_sym(m0, v0), std::invalid_argument);
+
+  matrix_d m1(1, 1);
+  m1 << 2;
+  vector_d v1(1);
+  v1 << 2;
+  EXPECT_FLOAT_EQ(8, quad_form_sym(m1, v1));
+
   matrix_d ad(4, 4);
   vector_d bd(4);
   double res;
@@ -31,34 +66,10 @@ TEST(MathMatrixPrim, quad_form_sym_vec) {
   ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
       1.0, 112.0;
 
-  // double-double
-  res = stan::math::quad_form_sym(ad, bd);
+  res = quad_form_sym(ad, bd);
   EXPECT_FLOAT_EQ(25433, res);
 }
 
-TEST(MathMatrixPrim, quad_form_sym_symmetry) {
-  using stan::math::matrix_d;
-  using stan::math::quad_form_sym;
-
-  matrix_d ad(4, 4);
-  matrix_d bd(4, 2);
-
-  bd << 100, 10, 0, 1, -3, -3, 5, 2;
-  ad << 2.0, 3.0, 4.0, 5.0, 3.0, 10.0, 2.0, 2.0, 4.0, 2.0, 7.0, 1.0, 5.0, 2.0,
-      1.0, 112.0;
-
-  // double-double
-  matrix_d resd = quad_form_sym(ad, bd);
-  EXPECT_EQ(resd(1, 0), resd(0, 1));
-
-  bd.resize(4, 3);
-  bd << 100, 10, 11, 0, 1, 12, -3, -3, 34, 5, 2, 44;
-  resd = quad_form_sym(ad, bd);
-  EXPECT_EQ(resd(1, 0), resd(0, 1));
-  EXPECT_EQ(resd(2, 0), resd(0, 2));
-  EXPECT_EQ(resd(2, 1), resd(1, 2));
-}
-
 TEST(MathMatrixPrim, quad_form_sym_asymmetric) {
   using stan::math::matrix_d;
 
@@ -69,6 +80,5 @@ TEST(MathMatrixPrim, quad_form_sym_asymmetric) {
   ad << 2.0, 3.0, 4.0, 5.0, 6.0, 10.0, 2.0, 2.0, 7.0, 2.0, 7.0, 1.0, 8.0, 2.0,
       1.0, 112.0;
 
-  // double-double
   EXPECT_THROW(stan::math::quad_form_sym(ad, bd), std::domain_error);
 }
diff --git a/test/unit/math/prim/fun/quad_form_test.cpp b/test/unit/math/prim/fun/quad_form_test.cpp
index 92d0d53d94d..11f9fc7ea6b 100644
--- a/test/unit/math/prim/fun/quad_form_test.cpp
+++ b/test/unit/math/prim/fun/quad_form_test.cpp
@@ -3,6 +3,24 @@
 
 TEST(MathMatrixPrim, quad_form_mat) {
   using stan::math::matrix_d;
+  using stan::math::quad_form;
+
+  matrix_d resd;
+  matrix_d m0;
+  EXPECT_THROW(quad_form(m0, m0), std::invalid_argument);
+
+  matrix_d m1(1, 1);
+  m1 << 2;
+  resd = quad_form(m1, m1);
+  EXPECT_FLOAT_EQ(8, resd(0));
+
+  matrix_d m2(1, 2);
+  m2 << 1, 2;
+  resd = quad_form(m1, m2);
+  EXPECT_FLOAT_EQ(2, resd(0, 0));
+  EXPECT_FLOAT_EQ(4, resd(0, 1));
+  EXPECT_FLOAT_EQ(4, resd(1, 0));
+  EXPECT_FLOAT_EQ(8, resd(1, 1));
 
   matrix_d ad(4, 4);
   matrix_d bd(4, 2);
@@ -11,8 +29,7 @@ TEST(MathMatrixPrim, quad_form_mat) {
   ad << 2.0, 3.0, 4.0, 5.0, 6.0, 10.0, 2.0, 2.0, 7.0, 2.0, 7.0, 1.0, 8.0, 2.0,
       1.0, 112.0;
 
-  // double-double
-  matrix_d resd = stan::math::quad_form(ad, bd);
+  resd = quad_form(ad, bd);
   EXPECT_FLOAT_EQ(26033, resd(0, 0));
   EXPECT_FLOAT_EQ(3456, resd(0, 1));
   EXPECT_FLOAT_EQ(3396, resd(1, 0));
@@ -21,8 +38,19 @@ TEST(MathMatrixPrim, quad_form_mat) {
 
 TEST(MathMatrixPrim, quad_form_vec) {
   using stan::math::matrix_d;
+  using stan::math::quad_form;
   using stan::math::vector_d;
 
+  matrix_d m0;
+  vector_d v0;
+  EXPECT_THROW(quad_form(m0, v0), std::invalid_argument);
+
+  matrix_d m1(1, 1);
+  m1 << 2;
+  vector_d v1(1);
+  v1 << 2;
+  EXPECT_FLOAT_EQ(8, quad_form(m1, v1));
+
   matrix_d ad(4, 4);
   vector_d bd(4);
   double res;
@@ -31,7 +59,6 @@ TEST(MathMatrixPrim, quad_form_vec) {
   ad << 2.0, 3.0, 4.0, 5.0, 6.0, 10.0, 2.0, 2.0, 7.0, 2.0, 7.0, 1.0, 8.0, 2.0,
       1.0, 112.0;
 
-  // double-double
   res = stan::math::quad_form(ad, bd);
   EXPECT_FLOAT_EQ(26033, res);
 }