stan-dev
Cleanup accumulator?
Description
It's an extra cost, but when we use accumulator.add(Vector) method in a Stan model should the sum we do there use log_sum_exp() since we are summing a bunch of logged values? Or could we have a log_accumulator that does that? That accumulator also has a buffer_ it feeds things into and then does a sum after the buffer is full. We could also use log_sum_exp() when we do the sum over those as well.
Expected Output
Hopefully nicer numerics?
Current Version:
v4.4.0
We can't use log_sum_exp here since we're after the sum of the logs, not the log of the sum of the exps.
But it did make me have a look at the accumulator() and sum() code, and I think we should change the std::vector<T> overload to use Eigen::Map<> and then sum, rather than std::accumulate: https://godbolt.org/z/zhe3eMWv7
Especially since accumulator sums its buffer using the std::vector<> method, which doesn't look to get the vectorised benefits that Eigen does
A couple of other accumulator thoughts.
At the moment, Eigen types are summed and then appended to the buffer, but std::vectorstd::vector<container> types:
template <typename S, require_matrix_t<S>* = nullptr>
inline void add(const S& m) {
buf_.push_back(stan::math::sum(m));
}
template <typename S>
inline void add(const std::vector<S>& xs) {
for (size_t i = 0; i < xs.size(); ++i) {
this->add(xs[i]);
}
}
If we change the templates such that the Eigen overload also takes non-nested std:vectors and the std::vector overload is only for nested containers, then this should dramatically improve performance with larger std::vector inputs. It would also reduce memory usage, since the non-nested std::vectors are no longer just being appended to the buffer.
One last thought, would there be any benefit from std::move()-ing the result?:
template <typename S, require_matrix_t<S>* = nullptr>
inline void add(const S& m) {
buf_.push_back(std::move(stan::math::sum(m)));
}
We can't use log_sum_exp here since we're after the sum of the logs, not the log of the sum of the exps.
Ooops yes totally forgot, forget this idea.
But it did make me have a look at the accumulator() and sum() code, and I think we should change the std::vector<T> overload to use Eigen::Map<> and then sum, rather than std::accumulate: https://godbolt.org/z/zhe3eMWv7
Especially since accumulator sums its buffer using the std::vector<> method, which doesn't look to get the vectorised benefits that Eigen does
Yes I like this idea!
If we change the templates such that the Eigen overload also takes non-nested std:vectors and the std::vector overload is only for nested containers, then this should dramatically improve performance with larger std::vector inputs. It would also reduce memory usage, since the non-nested std::vectors are no longer just being appended to the buffer.
Yes I agree that would be a good change!
One last thought, would there be any benefit from std::move()-ing the result?:
Not here, we only need to move if the container we are moving has a good amount of dynamically allocated memory. In the code above sum() returns a scalar so it would not need to be moved.
You can think of moves like it's giving ownership of deleting dynamically allocated memory to a new object from an old object.
struct foo {
double* data_;//
std::size_t n_;
// Copy makes a copy of the memory
foo(foo& x) : data_(copy(x.data_, n)), n_(x.n_) {}
foo(foo&& x) : data_(x.data_), n_(x.n_) {
// Moves takes the pointer and then nulls out the input's pointer
// This means the new object is now in charge of calling delete
x.data_ = nullptr;
}
}
Great! I've just about finished up a PR with these changes, will open that later today