Skip to content

use double loop for operator-/+ #2223

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

Merged
merged 25 commits into from
Dec 12, 2020
Merged

use double loop for operator-/+ #2223

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
4c5f2e5
use double loop for operator-/+ to avoid linear access bug
SteveBronder Dec 1, 2020
3172588
adds tests for x + x.transpose
SteveBronder Dec 1, 2020
be3c985
[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…
stan-buildbot Dec 1, 2020
4f8e4bf
fix for elt multiply and divide
SteveBronder Dec 5, 2020
56e03a7
Merge commit '7d93323f57b958950b67d34122649cd8fe5e9acd' into HEAD
yashikno Dec 5, 2020
cb9af42
[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…
stan-buildbot Dec 5, 2020
fa6cc83
Added some tests (Issue #2213)
bbbales2 Dec 8, 2020
a6e05bd
Merge commit 'b49445de22d97f217fd4711a39b8da86157ab868' into HEAD
yashikno Dec 8, 2020
7908194
[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…
stan-buildbot Dec 8, 2020
90fd5a8
adds return_var_matrix_t to be used with reverse mode matrix speciali…
SteveBronder Dec 8, 2020
7dcc621
make block elt_multiply/divide tests
SteveBronder Dec 8, 2020
3352711
[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…
stan-buildbot Dec 8, 2020
6d3c6d1
make block elt_multiply/divide tests
SteveBronder Dec 9, 2020
2bce5b6
update multiply for single to double loop
SteveBronder Dec 9, 2020
d1f80e1
Changed test names (Issue #2245)
bbbales2 Dec 9, 2020
57ae03f
Updated docs on return and promote meta functions (Issue #2245)
bbbales2 Dec 9, 2020
fd1e84e
Fixed some indexing in test_ad_matvar (Issue #2245)
bbbales2 Dec 9, 2020
7e713b2
fix softmax return
SteveBronder Dec 9, 2020
483b9ce
remove plain_type from operator subtract
SteveBronder Dec 9, 2020
ccc6111
fix loop in expect_ad_matvar when testing gradients
SteveBronder Dec 9, 2020
6e46026
[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…
stan-buildbot Dec 9, 2020
5ba84ce
fix return type for colums_dot_self
SteveBronder Dec 9, 2020
92bf783
[Jenkins] auto-formatting by clang-format version 6.0.0-1ubuntu2~16.0…
stan-buildbot Dec 9, 2020
5c53bea
update inverse to return same type for early return
SteveBronder Dec 11, 2020
02eaf83
Merge branch 'bugfix/linear-access-add-subtract' of github.com:stan-d…
SteveBronder Dec 11, 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
31 changes: 18 additions & 13 deletions stan/math/rev/core/operator_addition.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,15 +119,17 @@ template <typename VarMat1, typename VarMat2,
inline auto add(const VarMat1& a, const VarMat2& b) {
check_matching_dims("add", "a", a, "b", b);
using op_ret_type = decltype(a.val() + b.val());
using ret_type = promote_var_matrix_t<op_ret_type, VarMat1, VarMat2>;
using ret_type = return_var_matrix_t<op_ret_type, VarMat1, VarMat2>;
arena_t<VarMat1> arena_a(a);
arena_t<VarMat2> arena_b(b);
arena_t<ret_type> ret(arena_a.val() + arena_b.val());
reverse_pass_callback([ret, arena_a, arena_b]() mutable {
for (Eigen::Index i = 0; i < ret.size(); ++i) {
const auto ref_adj = ret.adj().coeffRef(i);
arena_a.adj().coeffRef(i) += ref_adj;
arena_b.adj().coeffRef(i) += ref_adj;
for (Eigen::Index j = 0; j < ret.cols(); ++j) {
for (Eigen::Index i = 0; i < ret.rows(); ++i) {
const auto ref_adj = ret.adj().coeffRef(i, j);
arena_a.adj().coeffRef(i, j) += ref_adj;
arena_b.adj().coeffRef(i, j) += ref_adj;
}
}
});
return ret_type(ret);
Expand All @@ -151,7 +153,7 @@ inline auto add(const VarMat& a, const Arith& b) {
}
using op_ret_type
= decltype((a.val().array() + as_array_or_scalar(b)).matrix());
using ret_type = promote_var_matrix_t<op_ret_type, VarMat>;
using ret_type = return_var_matrix_t<op_ret_type, VarMat>;
arena_t<VarMat> arena_a(a);
arena_t<ret_type> ret(arena_a.val().array() + as_array_or_scalar(b));
reverse_pass_callback(
Expand Down Expand Up @@ -188,7 +190,7 @@ template <typename Var, typename EigMat,
require_var_vt<std::is_arithmetic, Var>* = nullptr,
require_eigen_vt<std::is_arithmetic, EigMat>* = nullptr>
inline auto add(const Var& a, const EigMat& b) {
using ret_type = promote_scalar_t<var, EigMat>;
using ret_type = return_var_matrix_t<EigMat>;
arena_t<ret_type> ret(a.val() + b.array());
reverse_pass_callback([ret, a]() mutable { a.adj() += ret.adj().sum(); });
return ret_type(ret);
Expand Down Expand Up @@ -224,16 +226,19 @@ template <typename Var, typename VarMat,
require_var_vt<std::is_arithmetic, Var>* = nullptr,
require_rev_matrix_t<VarMat>* = nullptr>
inline auto add(const Var& a, const VarMat& b) {
using ret_type = return_var_matrix_t<VarMat>;
arena_t<VarMat> arena_b(b);
arena_t<VarMat> ret(a.val() + arena_b.val().array());
arena_t<ret_type> ret(a.val() + arena_b.val().array());
reverse_pass_callback([ret, a, arena_b]() mutable {
for (Eigen::Index i = 0; i < arena_b.size(); ++i) {
const auto ret_adj = ret.adj().coeffRef(i);
a.adj() += ret_adj;
arena_b.adj().coeffRef(i) += ret_adj;
for (Eigen::Index j = 0; j < ret.cols(); ++j) {
for (Eigen::Index i = 0; i < ret.rows(); ++i) {
const auto ret_adj = ret.adj().coeffRef(i, j);
a.adj() += ret_adj;
arena_b.adj().coeffRef(i, j) += ret_adj;
}
}
});
return plain_type_t<VarMat>(ret);
return ret_type(ret);
}

/**
Expand Down
50 changes: 29 additions & 21 deletions stan/math/rev/core/operator_subtraction.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -131,15 +131,17 @@ template <typename VarMat1, typename VarMat2,
inline auto subtract(const VarMat1& a, const VarMat2& b) {
check_matching_dims("subtract", "a", a, "b", b);
using op_ret_type = decltype(a.val() - b.val());
using ret_type = promote_var_matrix_t<op_ret_type, VarMat1, VarMat2>;
using ret_type = return_var_matrix_t<op_ret_type, VarMat1, VarMat2>;
arena_t<VarMat1> arena_a = a;
arena_t<VarMat2> arena_b = b;
arena_t<ret_type> ret((arena_a.val() - arena_b.val()));
reverse_pass_callback([ret, arena_a, arena_b]() mutable {
for (Eigen::Index i = 0; i < ret.size(); ++i) {
const auto ret_adj = ret.adj().coeffRef(i);
arena_a.adj().coeffRef(i) += ret_adj;
arena_b.adj().coeffRef(i) -= ret_adj;
for (Eigen::Index j = 0; j < ret.cols(); ++j) {
for (Eigen::Index i = 0; i < ret.rows(); ++i) {
const auto ref_adj = ret.adj().coeffRef(i, j);
arena_a.adj().coeffRef(i, j) += ref_adj;
arena_b.adj().coeffRef(i, j) -= ref_adj;
}
}
});
return ret_type(ret);
Expand All @@ -163,7 +165,7 @@ inline auto subtract(const VarMat& a, const Arith& b) {
}
using op_ret_type = plain_type_t<decltype(
(a.val().array() - as_array_or_scalar(b)).matrix())>;
using ret_type = promote_var_matrix_t<op_ret_type, VarMat>;
using ret_type = return_var_matrix_t<op_ret_type, VarMat>;
arena_t<VarMat> arena_a = a;
arena_t<ret_type> ret(arena_a.val().array() - as_array_or_scalar(b));
reverse_pass_callback(
Expand All @@ -189,7 +191,7 @@ inline auto subtract(const Arith& a, const VarMat& b) {
}
using op_ret_type = plain_type_t<decltype(
(as_array_or_scalar(a) - b.val().array()).matrix())>;
using ret_type = promote_var_matrix_t<op_ret_type, VarMat>;
using ret_type = return_var_matrix_t<op_ret_type, VarMat>;
arena_t<VarMat> arena_b = b;
arena_t<ret_type> ret(as_array_or_scalar(a) - arena_b.val().array());
reverse_pass_callback(
Expand All @@ -210,7 +212,7 @@ template <typename Var, typename EigMat,
require_var_vt<std::is_arithmetic, Var>* = nullptr,
require_eigen_vt<std::is_arithmetic, EigMat>* = nullptr>
inline auto subtract(const Var& a, const EigMat& b) {
using ret_type = promote_scalar_t<var, EigMat>;
using ret_type = return_var_matrix_t<EigMat>;
arena_t<ret_type> ret(a.val() - b.array());
reverse_pass_callback([ret, a]() mutable { a.adj() += ret.adj().sum(); });
return ret_type(ret);
Expand All @@ -229,7 +231,7 @@ template <typename EigMat, typename Var,
require_eigen_vt<std::is_arithmetic, EigMat>* = nullptr,
require_var_vt<std::is_arithmetic, Var>* = nullptr>
inline auto subtract(const EigMat& a, const Var& b) {
using ret_type = promote_scalar_t<var, EigMat>;
using ret_type = return_var_matrix_t<EigMat>;
arena_t<ret_type> ret(a.array() - b.val());
reverse_pass_callback([ret, b]() mutable { b.adj() -= ret.adj().sum(); });
return ret_type(ret);
Expand All @@ -249,16 +251,19 @@ template <typename Var, typename VarMat,
require_var_vt<std::is_arithmetic, Var>* = nullptr,
require_rev_matrix_t<VarMat>* = nullptr>
inline auto subtract(const Var& a, const VarMat& b) {
using ret_type = return_var_matrix_t<VarMat>;
arena_t<VarMat> arena_b(b);
arena_t<VarMat> ret(a.val() - arena_b.val().array());
arena_t<ret_type> ret(a.val() - arena_b.val().array());
reverse_pass_callback([ret, a, arena_b]() mutable {
for (Eigen::Index i = 0; i < ret.size(); ++i) {
auto ret_adj = ret.adj().coeff(i);
a.adj() += ret_adj;
arena_b.adj().coeffRef(i) -= ret_adj;
for (Eigen::Index j = 0; j < ret.cols(); ++j) {
for (Eigen::Index i = 0; i < ret.rows(); ++i) {
auto ret_adj = ret.adj().coeff(i, j);
a.adj() += ret_adj;
arena_b.adj().coeffRef(i, j) -= ret_adj;
}
}
});
return plain_type_t<VarMat>(ret);
return ret_type(ret);
}

/**
Expand All @@ -275,16 +280,19 @@ template <typename Var, typename VarMat,
require_rev_matrix_t<VarMat>* = nullptr,
require_var_vt<std::is_arithmetic, Var>* = nullptr>
inline auto subtract(const VarMat& a, const Var& b) {
using ret_type = return_var_matrix_t<VarMat>;
arena_t<VarMat> arena_a(a);
arena_t<VarMat> ret(arena_a.val().array() - b.val());
arena_t<ret_type> ret(arena_a.val().array() - b.val());
reverse_pass_callback([ret, b, arena_a]() mutable {
for (Eigen::Index i = 0; i < ret.size(); ++i) {
const auto ret_adj = ret.adj().coeff(i);
arena_a.adj().coeffRef(i) += ret_adj;
b.adj() -= ret_adj;
for (Eigen::Index j = 0; j < ret.cols(); ++j) {
for (Eigen::Index i = 0; i < ret.rows(); ++i) {
const auto ret_adj = ret.adj().coeff(i, j);
arena_a.adj().coeffRef(i, j) += ret_adj;
b.adj() -= ret_adj;
}
}
});
return plain_type_t<VarMat>(ret);
return ret_type(ret);
}

template <typename T1, typename T2,
Expand Down
2 changes: 1 addition & 1 deletion stan/math/rev/fun/cholesky_decompose.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -271,7 +271,7 @@ inline auto cholesky_decompose(const EigMat& A) {
template <typename T, require_var_matrix_t<T>* = nullptr>
inline auto cholesky_decompose(const T& A) {
check_symmetric("cholesky_decompose", "A", A.val());
T L = cholesky_decompose(A.val());
plain_type_t<T> L = cholesky_decompose(A.val());
if (A.rows() <= 35) {
reverse_pass_callback(internal::unblocked_cholesky_lambda(L.val(), L, A));
} else {
Expand Down
7 changes: 3 additions & 4 deletions stan/math/rev/fun/columns_dot_product.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -63,10 +63,9 @@ template <typename Mat1, typename Mat2,
require_any_var_matrix_t<Mat1, Mat2>* = nullptr>
inline auto columns_dot_product(const Mat1& v1, const Mat2& v2) {
check_matching_sizes("columns_dot_product", "v1", v1, "v2", v2);

using return_t = promote_var_matrix_t<
decltype((v1.val().array() * v2.val().array()).colwise().sum().matrix()),
Mat1, Mat2>;
using inner_return_t = decltype(
(value_of(v1).array() * value_of(v2).array()).colwise().sum().matrix());
using return_t = return_var_matrix_t<inner_return_t, Mat1, Mat2>;

if (!is_constant<Mat1>::value && !is_constant<Mat2>::value) {
arena_t<promote_scalar_t<var, Mat1>> arena_v1 = v1;
Expand Down
5 changes: 3 additions & 2 deletions stan/math/rev/fun/columns_dot_self.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,9 @@ inline Eigen::Matrix<var, 1, Mat::ColsAtCompileTime> columns_dot_self(
*/
template <typename Mat, require_var_matrix_t<Mat>* = nullptr>
inline auto columns_dot_self(const Mat& x) {
using ret_type = plain_type_t<decltype(x.val().colwise().squaredNorm())>;
var_value<ret_type> res = x.val().colwise().squaredNorm();
using ret_type
= return_var_matrix_t<decltype(x.val().colwise().squaredNorm()), Mat>;
arena_t<ret_type> res = x.val().colwise().squaredNorm();
if (x.size() >= 0) {
reverse_pass_callback([res, x]() mutable {
x.adj() += x.val() * (2 * res.adj()).asDiagonal();
Expand Down
2 changes: 1 addition & 1 deletion stan/math/rev/fun/divide.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -191,7 +191,7 @@ template <typename Mat, typename Scal, require_var_matrix_t<Mat>* = nullptr,
inline auto divide(const Mat& m, const Scal& c) {
double invc = 1.0 / value_of(c);

Mat res = invc * m.val();
plain_type_t<Mat> res = invc * m.val();

reverse_pass_callback([m, c, res, invc]() mutable {
m.adj() += invc * res.adj();
Expand Down
16 changes: 9 additions & 7 deletions stan/math/rev/fun/elt_divide.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,17 +28,19 @@ auto elt_divide(const Mat1& m1, const Mat2& m2) {
check_matching_dims("elt_divide", "m1", m1, "m2", m2);
using inner_ret_type
= decltype((value_of(m1).array() / value_of(m2).array()).matrix());
using ret_type = promote_var_matrix_t<inner_ret_type, Mat1, Mat2>;
using ret_type = return_var_matrix_t<inner_ret_type, Mat1, Mat2>;
if (!is_constant<Mat1>::value && !is_constant<Mat2>::value) {
arena_t<promote_scalar_t<var, Mat1>> arena_m1 = m1;
arena_t<promote_scalar_t<var, Mat2>> arena_m2 = m2;
arena_t<ret_type> ret(arena_m1.val().array() / arena_m2.val().array());
reverse_pass_callback([ret, arena_m1, arena_m2]() mutable {
for (Eigen::Index i = 0; i < arena_m2.size(); ++i) {
const auto ret_adj = ret.adj().coeffRef(i);
arena_m1.adj().coeffRef(i) += ret_adj / arena_m2.val().coeff(i);
arena_m2.adj().coeffRef(i)
-= ret.val().coeff(i) * ret_adj / arena_m2.val().coeff(i);
for (Eigen::Index j = 0; j < arena_m2.cols(); ++j) {
for (Eigen::Index i = 0; i < arena_m2.rows(); ++i) {
const auto ret_div
= ret.adj().coeff(i, j) / arena_m2.val().coeff(i, j);
arena_m1.adj().coeffRef(i, j) += ret_div;
arena_m2.adj().coeffRef(i, j) -= ret.val().coeff(i, j) * ret_div;
}
}
});
return ret_type(ret);
Expand Down Expand Up @@ -76,7 +78,7 @@ auto elt_divide(const Mat1& m1, const Mat2& m2) {
template <typename Scal, typename Mat, require_stan_scalar_t<Scal>* = nullptr,
require_var_matrix_t<Mat>* = nullptr>
auto elt_divide(Scal s, const Mat& m) {
Mat res = value_of(s) / m.val().array();
plain_type_t<Mat> res = value_of(s) / m.val().array();

reverse_pass_callback([m, s, res]() mutable {
m.adj().array() -= res.val().array() * res.adj().array() / m.val().array();
Expand Down
12 changes: 7 additions & 5 deletions stan/math/rev/fun/elt_multiply.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -27,16 +27,18 @@ template <typename Mat1, typename Mat2,
auto elt_multiply(const Mat1& m1, const Mat2& m2) {
check_matching_dims("elt_multiply", "m1", m1, "m2", m2);
using inner_ret_type = decltype(value_of(m1).cwiseProduct(value_of(m2)));
using ret_type = promote_var_matrix_t<inner_ret_type, Mat1, Mat2>;
using ret_type = return_var_matrix_t<inner_ret_type, Mat1, Mat2>;
if (!is_constant<Mat1>::value && !is_constant<Mat2>::value) {
arena_t<promote_scalar_t<var, Mat1>> arena_m1 = m1;
arena_t<promote_scalar_t<var, Mat2>> arena_m2 = m2;
arena_t<ret_type> ret(arena_m1.val().cwiseProduct(arena_m2.val()));
reverse_pass_callback([ret, arena_m1, arena_m2]() mutable {
for (Eigen::Index i = 0; i < arena_m2.size(); ++i) {
const auto ret_adj = ret.adj().coeffRef(i);
arena_m1.adj().coeffRef(i) += arena_m2.val().coeff(i) * ret_adj;
arena_m2.adj().coeffRef(i) += arena_m1.val().coeff(i) * ret_adj;
for (Eigen::Index j = 0; j < arena_m2.cols(); ++j) {
for (Eigen::Index i = 0; i < arena_m2.rows(); ++i) {
const auto ret_adj = ret.adj().coeffRef(i, j);
arena_m1.adj().coeffRef(i, j) += arena_m2.val().coeff(i, j) * ret_adj;
arena_m2.adj().coeffRef(i, j) += arena_m1.val().coeff(i, j) * ret_adj;
}
}
});
return ret_type(ret);
Expand Down
11 changes: 6 additions & 5 deletions stan/math/rev/fun/inverse.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,22 +20,23 @@ namespace math {
* @throw std::invalid_argument if the matrix is not square.
*/
template <typename T, require_rev_matrix_t<T>* = nullptr>
inline plain_type_t<T> inverse(const T& m) {
inline auto inverse(const T& m) {
check_square("inverse", "m", m);

if (m.size() == 0) {
return m;
using ret_type = return_var_matrix_t<T>;
if (unlikely(m.size() == 0)) {
return ret_type(m);
}

arena_t<T> arena_m = m;
arena_t<promote_scalar_t<double, T>> res_val = arena_m.val().inverse();
arena_t<T> res = res_val;
arena_t<ret_type> res = res_val;

reverse_pass_callback([res, res_val, arena_m]() mutable {
arena_m.adj() -= res_val.transpose() * res.adj_op() * res_val.transpose();
});

return res;
return ret_type(res);
}

} // namespace math
Expand Down
6 changes: 3 additions & 3 deletions stan/math/rev/fun/matrix_power.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,8 @@ inline plain_type_t<T> matrix_power(const T& M, const int n) {
for (size_t i = 2; i <= n; ++i) {
arena_powers[i] = arena_powers[1] * arena_powers[i - 1];
}

arena_t<plain_type_t<T>> res = arena_powers[arena_powers.size() - 1];
using ret_type = return_var_matrix_t<T>;
arena_t<ret_type> res = arena_powers[arena_powers.size() - 1];

reverse_pass_callback([arena_M, n, res, arena_powers]() mutable {
const auto& M_val = arena_powers[1];
Expand All @@ -66,7 +66,7 @@ inline plain_type_t<T> matrix_power(const T& M, const int n) {
arena_M.adj() += adj_M + adj_C;
});

return res;
return ret_type(res);
Copy link
Member

Choose a reason for hiding this comment

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

Good catch!

}

} // namespace math
Expand Down
Loading