for loops

[Ben Goodrich]

There are a ton of them in Stan :) . I have long thought that there was a problem with many of them, and that they could be replaced by Eigen's template expressions (unless they involve mixed operations or special cases relating to auto-diff). I pointed this out to Bob many months ago, and he claimed that the C-style for loops were essentially what the compiler would produce for the equivalent Eigen template expression. But, from my (limited) understanding of

http://eigen.tuxfamily.org/dox/TopicInsideEigenExample.html

this is not correct.

They go into exhaustive detail about how Eigen (can) take aligned floating-point numbers and loop over packets of 4 values, rather than values. In particular, to add vectors v and w of length 50 (4 * 12 + 2) and assign them to the vector u, their expression templates ultimately yield code like

  for(int i = 0; i < 4*(size/4); i+=4) u.packet(i)  = v.packet(i) + w.packet(i);
  for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i];

which is almost 4 times faster than

  for(int i = 0; i < size; i++) u[i] = v[i] + w[i];

But in Stan, there are many places where we are "Doing the C Thing" instead of "Doing the Eigen Thing". To take a simple example, starting at line 1007 of src/stan/math/matrix.hpp, we have

    template <typename T>
    inline T sum(const std::vector<T>& xs) {
      if (xs.size() == 0) return 0;
      T sum(xs[0]);
      for (size_t i = 1; i < xs.size(); ++i)
        sum += xs[i];
      return sum;
    }

which, if I understanding the packet-aligning business that Eigen does (and have the Map / Array syntax right) should be

    template <typename T>
    inline T sum(const std::vector<T>& xs) {
      return Eigen::Map< Eigen::Array<T, Dynamic, 1> >(&xs[0], xs.size(), 1).sum();
    }

And there are lots of C-style loops in Stan. But I really don't understand compilers or vectorization in the computer-science sense, so if I am wrong, can someone tell me why?

Ben

[Marcus Brubaker]

On Sat, Oct 13, 2012 at 6:22 PM, Ben Goodrich <goodri...@gmail.com> wrote:

There are a ton of them in Stan :) . I have long thought that there was a problem with many of them, and that they could be replaced by Eigen's template expressions (unless they involve mixed operations or special cases relating to auto-diff). I pointed this out to Bob many months ago, and he claimed that the C-style for loops were essentially what the compiler would produce for the equivalent Eigen template expression. But, from my (limited) understanding of

http://eigen.tuxfamily.org/dox/TopicInsideEigenExample.html

this is not correct.

They go into exhaustive detail about how Eigen (can) take aligned floating-point numbers and loop over packets of 4 values, rather than values. In particular, to add vectors v and w of length 50 (4 * 12 + 2) and assign them to the vector u, their expression templates ultimately yield code like

  for(int i = 0; i < 4*(size/4); i+=4) u.packet(i)  = v.packet(i) + w.packet(i);

  for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i];

which is almost 4 times faster than

  for(int i = 0; i < size; i++) u[i] = v[i] + w[i];

This is a result of explicitly using SIMD operations (implied by the packet(i) calls) which can be ~ 2-4 times faster for an operation like that.  Note however that SIMD instructions only work for basic types, namely float and double.  There is little to be done easily to exploit SIMD instructions for operations on var types.
 

But in Stan, there are many places where we are "Doing the C Thing" instead of "Doing the Eigen Thing". To take a simple example, starting at line 1007 of src/stan/math/matrix.hpp, we have

    template <typename T>
    inline T sum(const std::vector<T>& xs) {
      if (xs.size() == 0) return 0;
      T sum(xs[0]);
      for (size_t i = 1; i < xs.size(); ++i)
        sum += xs[i];
      return sum;
    }

which, if I understanding the packet-aligning business that Eigen does (and have the Map / Array syntax right) should be

    template <typename T>
    inline T sum(const std::vector<T>& xs) {
      return Eigen::Map< Eigen::Array<T, Dynamic, 1> >(&xs[0], xs.size(), 1).sum();
    }

And there are lots of C-style loops in Stan. But I really don't understand compilers or vectorization in the computer-science sense, so if I am wrong, can someone tell me why?

You're not wrong, most of those loops will be slower than Eigens code.  Sometimes compilers can be clever and use SIMD instructions automatically, but that is hard and doesn't always work depending on how the loop is written.  That said, I expect that the speedup of replacing every loop like that with an appropriate call to Eigen would be relatively small for most models because the overhead cost is going to be dominated by computations with vars which would be unaffected.

Cheers,
Marcus

[Ben Goodrich]

On Saturday, October 13, 2012 6:34:59 PM UTC-4, Marcus Brubaker wrote:

This is a result of explicitly using SIMD operations (implied by the packet(i) calls) which can be ~ 2-4 times faster for an operation like that.  Note however that SIMD instructions only work for basic types, namely float and double.  There is little to be done easily to exploit SIMD instructions for operations on var types.

But it is correct that in a templated function with template <typename T> is invoked when T is var, then the compiler will make a C-style loop that is no slower than what a person would do by hand?

Also, I think there may be some room for improvement in for loops where T is var but the function is only using the .val() portion. Like in normal.hpp now we have

      DoubleVectorView<true,T_scale> inv_sigma(length(sigma));
      DoubleVectorView<include_summand<propto,T_scale>::value,T_scale> log_sigma(length(sigma));
      for (size_t i = 0; i < length(sigma); i++) {
        inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
        if (include_summand<propto,T_scale>::value)
          log_sigma[i] = log(value_of(sigma_vec[i]));
      }

which could be

      Eigen::Map< Eigen::ArrayXd > sigma_Arr(&sigma[0], length(sigma), 1)::cast<double>();
      Eigen::ArrayXd inv_sigma = sigma_Arr.inverse();
      Eigen::ArrayXd log_sigma;
      if (include_summand<propto,T_scale>::value) log_sigma = sigma_Arr.log()


I am still trying to wrap my head around the magic of VectorViews, but they seem to be a bit like Eigen::Maps in the sense that they just point to the leading value of something else.

Ben

[Daniel Lee]

On Oct 13, 2012, at 3:56 PM, Ben Goodrich <goodri...@gmail.com> wrote:

On Saturday, October 13, 2012 6:34:59 PM UTC-4, Marcus Brubaker wrote:

This is a result of explicitly using SIMD operations (implied by the packet(i) calls) which can be ~ 2-4 times faster for an operation like that.  Note however that SIMD instructions only work for basic types, namely float and double.  There is little to be done easily to exploit SIMD instructions for operations on var types.

But it is correct that in a templated function with template <typename T> is invoked when T is var, then the compiler will make a C-style loop that is no slower than what a person would do by hand?

Also, I think there may be some room for improvement in for loops where T is var but the function is only using the .val() portion. Like in normal.hpp now we have

      DoubleVectorView<true,T_scale> inv_sigma(length(sigma));
      DoubleVectorView<include_summand<propto,T_scale>::value,T_scale> log_sigma(length(sigma));
      for (size_t i = 0; i < length(sigma); i++) {
        inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
        if (include_summand<propto,T_scale>::value)
          log_sigma[i] = log(value_of(sigma_vec[i]));
      }

which could be

      Eigen::Map< Eigen::ArrayXd > sigma_Arr(&sigma[0], length(sigma), 1)::cast<double>();
      Eigen::ArrayXd inv_sigma = sigma_Arr.inverse();
      Eigen::ArrayXd log_sigma;
      if (include_summand<propto,T_scale>::value) log_sigma = sigma_Arr.log()


Can you make that compile? Looking at it from my phone and just thinking about it makes me think that the code written would not work in general. See if toucan get the tests to compile after the change. 

I am still trying to wrap my head around the magic of VectorViews, but they seem to be a bit like Eigen::Maps in the sense that they just point to the leading value of something else.

Oh yes, VectorViews are a special type of magic that Matt and Bob figured out. They make loops syntactically valid over all the types that we allow for vectorization: T, vector<T>, matrix<T,dynamic,1>, matrix<T,1,dynamic>, where T can either be double or var. 

Using VectorViews, we can walk over the longest vector once and maintain a single function that takes in any combination of arguments listed. 

Ben

--
 
 

[Marcus Brubaker]

On Sat, Oct 13, 2012 at 6:56 PM, Ben Goodrich <goodri...@gmail.com> wrote:

On Saturday, October 13, 2012 6:34:59 PM UTC-4, Marcus Brubaker wrote:

This is a result of explicitly using SIMD operations (implied by the packet(i) calls) which can be ~ 2-4 times faster for an operation like that.  Note however that SIMD instructions only work for basic types, namely float and double.  There is little to be done easily to exploit SIMD instructions for operations on var types.

But it is correct that in a templated function with template <typename T> is invoked when T is var, then the compiler will make a C-style loop that is no slower than what a person would do by hand?

Yes, that should be correct overall.  There may be some clever things that Eigen does with prefetching and blocking to improve cache performance even with T var, but the main speedup comes from the use of SIMD instructions.
 

Also, I think there may be some room for improvement in for loops where T is var but the function is only using the .val() portion. Like in normal.hpp now we have

      DoubleVectorView<true,T_scale> inv_sigma(length(sigma));
      DoubleVectorView<include_summand<propto,T_scale>::value,T_scale> log_sigma(length(sigma));
      for (size_t i = 0; i < length(sigma); i++) {
        inv_sigma[i] = 1.0 / value_of(sigma_vec[i]);
        if (include_summand<propto,T_scale>::value)
          log_sigma[i] = log(value_of(sigma_vec[i]));
      }

which could be

      Eigen::Map< Eigen::ArrayXd > sigma_Arr(&sigma[0], length(sigma), 1)::cast<double>();
      Eigen::ArrayXd inv_sigma = sigma_Arr.inverse();
      Eigen::ArrayXd log_sigma;
      if (include_summand<propto,T_scale>::value) log_sigma = sigma_Arr.log()


I am still trying to wrap my head around the magic of VectorViews, but they seem to be a bit like Eigen::Maps in the sense that they just point to the leading value of something else.

I'm not really familiar with the VectorView stuff so I can't comment on that.

Cheers,
Marcus

[Ben Goodrich]

On Saturday, October 13, 2012 7:11:38 PM UTC-4, Daniel Lee wrote:

which could be

      Eigen::Map< Eigen::ArrayXd > sigma_Arr(&sigma[0], length(sigma), 1)::cast<double>();
      Eigen::ArrayXd inv_sigma = sigma_Arr.inverse();
      Eigen::ArrayXd log_sigma;
      if (include_summand<propto,T_scale>::value) log_sigma = sigma_Arr.log()


Can you make that compile? Looking at it from my phone and just thinking about it makes me think that the code written would not work in general. See if toucan get the tests to compile after the change. 

I left off the last semicolon :). But I will try for real.

Ben

[Bob Carpenter]

On 10/13/12 6:34 PM, Marcus Brubaker wrote:
> On Sat, Oct 13, 2012 at 6:22 PM, Ben Goodrich <goodri...@gmail.com <mailto:goodri...@gmail.com>> wrote:
>
> ...I pointed this out to Bob many months ago, and he claimed that the C-style for loops were essentially

> what the compiler would produce for the equivalent Eigen template expression. But, from my (limited) understanding of
>
> http://eigen.tuxfamily.org/dox/TopicInsideEigenExample.html
>
> this is not correct.
>
> They go into exhaustive detail about how Eigen (can) take aligned floating-point numbers and loop over packets of 4
> values, rather than values. In particular, to add vectors v and w of length 50 (4 * 12 + 2) and assign them to the
> vector u, their expression templates ultimately yield code like
>
> for(int i = 0; i < 4*(size/4); i+=4) u.packet(i) = v.packet(i) + w.packet(i);
>
> for(int i = 4*(size/4); i < size; i++) u[i] = v[i] + w[i];

Very cool. I retract my guess.

The above is for floats (32 bytes); doubles only do
two at a time instead of four. And it looks like
we'd need 16-byte alignment (128 bits). It looks like
Eigen does that under the hood for us.

That page was really useful in understanding other aspects
of how Eigen works. Thanks for pointing it out.

> But in Stan, there are many places where we are "Doing the C Thing" instead of "Doing the Eigen Thing".

It'd be great to fix this.

I only did it this way because I found it simpler
and didn't understand that they had deeper vectorization
inside.

> To take a
> simple example, starting at line 1007 of src/stan/math/matrix.hpp, we have
>
> |
> template <typename T>
> inline T sum(const std::vector<T>& xs) {
> if (xs.size() == 0) return 0;
> T sum(xs[0]);
> for (size_t i = 1; i < xs.size(); ++i)
> sum += xs[i];
> return sum;
> }
> |
>
> which, if I understanding the packet-aligning business that Eigen does (and have the Map / Array syntax right) should be
>
> |
> template <typename T>
> inline T sum(const std::vector<T>& xs) {
> return Eigen::Map< Eigen::Array<T, Dynamic, 1> >(&xs[0], xs.size(), 1).sum();
> }
> |
>
> And there are lots of C-style loops in Stan. But I really don't understand compilers or vectorization in the
> computer-science sense, so if I am wrong, can someone tell me why?

Nobody does! It's a rather empirical art from with a basis
partly in compilers and partly in CPU/bus/memory architecture.
The main point is that the CPU can do many operations in
parallel and it can prefetch data IF it can organize the data
into chunks properly and it can predict where it's going.

This book was really useful and I found it fun and it's
free online:

http://www.agner.org/optimize/optimizing_cpp.pdf

For example, in the Eigen-style code, you wind up constructing
new objects on the stack. Their constructors aren't free,
but they're relatively pretty cheap and should be dwarfed by
the vectorization win for larger vectors. Also, using two
loops can be worse than one or better -- sometimes the
compiler will itself do vectorization using the kind of stride
plus remainder business you saw in Eigen.

Partly it depends on context.

It also hugely changes over time as the relative speeds of
CPU, cache, memory and the availability of new instruction
sets, register organizations, etc. emerges.

I wasn't worrying too much about these issues because of
the following point from Marcus:

> You're not wrong, most of those loops will be slower than Eigens code. Sometimes compilers can be clever and use SIMD
> instructions automatically, but that is hard and doesn't always work depending on how the loop is written. That said, I
> expect that the speedup of replacing every loop like that with an appropriate call to Eigen would be relatively small
> for most models because the overhead cost is going to be dominated by computations with vars which would be unaffected.

This is why the unfolding of the auto-dif in the
matrix multiplies and in the probability function
code is such a big deal for speed.

The best way to get more speed in the other ops is to add
more analytic derivatives implemented in single chain()
functions. And then pulling out a double-based
matrix or vector to send to Eigen. At that point, all
the Eigen expression templates and vectorization should
help.

- Bob

[Marcus Brubaker]

On Sat, Oct 13, 2012 at 8:54 PM, Bob Carpenter <ca...@alias-i.com> wrote:

The best way to get more speed in the other ops is to add
more analytic derivatives implemented in single chain()
functions.  And then pulling out a double-based
matrix or vector to send to Eigen.  At that point, all
the Eigen expression templates and vectorization should
help.

One bit of low-hanging fruit for improving speed along these lines is to grep for places where we call x.sum() and replace it with sum(x).  sum() has a specialized vari, but .sum() calls Eigen's code which doesn't use that specialization.  There are a few of these hanging around the code base which could have a measurable impact on performance for bigger models or those with more data.

Cheers,
Marcus

[Ben Goodrich]

 I'll do that since it looks like I am the worst offender

goodrich@CYBERPOWERPC:/opt/stan$ git grep -n "\.sum" src/stan | cat

src/stan/math/matrix.hpp:1024:      return v.sum();

src/stan/math/matrix_error_handling.hpp:386:      if (fabs(1.0 - theta.sum()) > CONSTRAINT_TOLERANCE) {

src/stan/math/matrix_error_handling.hpp:390:        T_prob sum = theta.sum();

src/stan/prob/distributions/multivariate/continuous/dirichlet.hpp:52:        lp += lgamma(alpha.sum());

src/stan/prob/distributions/multivariate/continuous/inv_wishart.hpp:87:        lp += nu * S.llt().matrixLLT().diagonal().array().log().sum();

src/stan/prob/distributions/multivariate/continuous/inv_wishart.hpp:90:        lp -= (nu + k + 1.0) * L.diagonal().array().log().sum();

src/stan/prob/distributions/multivariate/continuous/inv_wishart.hpp:94:        lp -= 0.5 * elt_multiply(S, L).array().sum(); // trace(S * W^-1)

src/stan/prob/distributions/multivariate/continuous/lkj_corr.hpp:78:        lp += (eta - 1.0) * 2.0 * L.diagonal().array().log().sum();

src/stan/prob/distributions/multivariate/continuous/multi_normal.hpp:70:        lp -= L.diagonal().array().log().sum();

src/stan/prob/distributions/multivariate/continuous/multi_normal.hpp:158:        lp -= L.diagonal().array().log().sum() * y.rows();

src/stan/prob/distributions/multivariate/continuous/multi_normal.hpp:175:        lp -= 0.5 * columns_dot_self(half).sum();

src/stan/prob/distributions/multivariate/continuous/multi_student_t.hpp:100:        lp -= L.diagonal().array().log().sum();

src/stan/prob/distributions/multivariate/continuous/wishart.hpp:102:        lp -= nu * L_S.diagonal().array().log().sum();

src/stan/prob/distributions/multivariate/continuous/wishart.hpp:109:        lp -= 0.5 * elt_multiply(L_S, W).array().sum(); // trace(S^-1 * W)

src/stan/prob/distributions/multivariate/continuous/wishart.hpp:113:        lp += (nu - k - 1.0) * L_W.diagonal().array().log().sum();

src/stan/prob/transform.hpp:264:      log_prob += (sds.log().sum() + stan::math::LOG_2) * K;

Ben

[Bob Carpenter]

On 10/13/12 9:06 PM, Marcus Brubaker wrote:

> One bit of low-hanging fruit for improving speed along these lines is to grep for places where we call x.sum() and
> replace it with sum(x). sum() has a specialized vari, but .sum() calls Eigen's code which doesn't use that
> specialization. There are a few of these hanging around the code base which could have a measurable impact on
> performance for bigger models or those with more data.

There are a few we should definitely update right away, but many
of them are going to be optimal again when we convert back from
var to double ops. I put it at the top of the C++ API to do.

I think I'm going to add constraints for vectors/matrices
in declarations before moving on to optimizations.

/stan>grep -R ".sum()" src/stan
src/stan/math/matrix.hpp: return v.sum();
src/stan/math/matrix_error_handling.hpp: if (fabs(1.0 - theta.sum()) > CONSTRAINT_TOLERANCE) {
src/stan/math/matrix_error_handling.hpp: T_prob sum = theta.sum();
src/stan/prob/distributions/multivariate/continuous/dirichlet.hpp: lp += lgamma(alpha.sum());
src/stan/prob/distributions/multivariate/continuous/inv_wishart.hpp: lp += nu *
S.llt().matrixLLT().diagonal().array().log().sum();
src/stan/prob/distributions/multivariate/continuous/inv_wishart.hpp: lp -= (nu + k + 1.0) *
L.diagonal().array().log().sum();
src/stan/prob/distributions/multivariate/continuous/inv_wishart.hpp: lp -= 0.5 * elt_multiply(S,

L).array().sum(); // trace(S * W^-1)

src/stan/prob/distributions/multivariate/continuous/lkj_corr.hpp: lp += (eta - 1.0) * 2.0 *
L.diagonal().array().log().sum();
src/stan/prob/distributions/multivariate/continuous/multi_normal.hpp: lp -= L.diagonal().array().log().sum();
src/stan/prob/distributions/multivariate/continuous/multi_normal.hpp: lp -= L.diagonal().array().log().sum() *
y.rows();
src/stan/prob/distributions/multivariate/continuous/multi_normal.hpp: lp -= 0.5 * columns_dot_self(half).sum();
src/stan/prob/distributions/multivariate/continuous/multi_student_t.hpp: lp -= L.diagonal().array().log().sum();
src/stan/prob/distributions/multivariate/continuous/wishart.hpp: lp -= nu * L_S.diagonal().array().log().sum();
src/stan/prob/distributions/multivariate/continuous/wishart.hpp: lp -= 0.5 * elt_multiply(L_S, W).array().sum();
// trace(S^-1 * W)
src/stan/prob/distributions/multivariate/continuous/wishart.hpp: lp += (nu - k - 1.0) *
L_W.diagonal().array().log().sum();
src/stan/prob/transform.hpp: log_prob += (sds.log().sum() + stan::math::LOG_2) * K;

- Bob