Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Improve tests and documentation for quad_form and quad_form_sym #1698

Merged
merged 8 commits into from
Feb 13, 2020
Merged

Improve tests and documentation for quad_form and quad_form_sym #1698

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
ff16115
Split quad_form and quad_form_sym tests
mcol Feb 11, 2020
30b042d
Cleanup test files
mcol Feb 11, 2020
8f3256f
Fix quad_form_sym mix tests
mcol Feb 11, 2020
8b47978
Use require_var_t to simplify code
mcol Feb 11, 2020
5b5183b
Add fwd version of quad_form and mix tests
mcol Feb 11, 2020
065a1cd
Consolidate implementations and add docs
mcol Feb 11, 2020
702acb1
Improve documentation
mcol Feb 12, 2020
23eded0
Add more tests for boundary conditions
mcol Feb 12, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions stan/math/fwd/fun.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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>
Expand Down
66 changes: 66 additions & 0 deletions stan/math/fwd/fun/quad_form.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,66 @@
#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 {

/**
* 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
* @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, 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
*/
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));
}

/**
* 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
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

[question]
Can these be anything other than Eigen::Dynamic?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I thought so. In reality we can't compile if we hardcode sizes because of check_square:

In file included from test/unit/math/prim/fun/quad_form_test.cpp:1:
In file included from ./stan/math/prim.hpp:10:
In file included from ./stan/math/prim/fun.hpp:238:
./stan/math/prim/fun/quad_form.hpp:31:3: error: no matching function for call to
      'check_square'
  check_square("quad_form", "A", A);
  ^~~~~~~~~~~~
test/unit/math/prim/fun/quad_form_test.cpp:22:31: note: in instantiation of
      function template specialization 'stan::math::quad_form<5, 5, 5, 5, double>'
      requested here
  matrix_d res2 = stan::math::quad_form(A, A);
                              ^
./stan/math/prim/err/check_square.hpp:21:13: note: candidate template ignored:
      could not match -1 against 5
inline void check_square(
            ^
1 error generated.

* @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 (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
67 changes: 43 additions & 24 deletions stan/math/fwd/fun/quad_form_sym.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,40 +8,59 @@
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) {
/**
* 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
* @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, 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
*/
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));
bob-carpenter marked this conversation as resolved.
Show resolved Hide resolved
}

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) {
/**
* 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
* @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 (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));
bob-carpenter marked this conversation as resolved.
Show resolved Hide resolved
}
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

Expand Down
33 changes: 31 additions & 2 deletions stan/math/prim/fun/quad_form.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,23 @@ namespace stan {
namespace math {

/**
* Compute 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
* @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, 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
*/
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) {
Expand All @@ -18,6 +33,20 @@ inline Eigen::Matrix<T, CB, CB> quad_form(const Eigen::Matrix<T, RA, CA>& A,
return B.transpose() * 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
* @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 (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) {
Expand Down
31 changes: 31 additions & 0 deletions stan/math/prim/fun/quad_form_sym.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,23 @@
namespace stan {
namespace math {

/**
* 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
* @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, 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
*/
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) {
Expand All @@ -16,6 +33,20 @@ inline Eigen::Matrix<T, CB, CB> quad_form_sym(
return T(0.5) * (ret + ret.transpose());
}

/**
* 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
* @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 (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) {
Expand Down
53 changes: 42 additions & 11 deletions stan/math/rev/fun/quad_form.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,12 +95,29 @@ 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) {
/**
* 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
* @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, 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
*/
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);

Expand All @@ -110,11 +127,25 @@ 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) {
/**
* 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
* @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 (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,
const Eigen::Matrix<Tb, Rb, 1>& B) {
check_square("quad_form", "A", A);
check_multiplicable("quad_form", "A", A, "B", B);

Expand Down
33 changes: 33 additions & 0 deletions stan/math/rev/fun/quad_form_sym.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,24 @@
namespace stan {
namespace math {

/**
* 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
* @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, 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
*/
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(
Expand All @@ -24,6 +42,21 @@ inline Eigen::Matrix<var, Cb, Cb> quad_form_sym(
return baseVari->impl_->C_;
}

/**
* 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
* @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 (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,
Expand Down
Loading