Use Eigen::Map to replace arr functions

Description

As initially discussed in this PR, there is some code duplication in having the same functions implemented for std::vector and Eigen separately. This could be avoided by using a wrapper which compiles to a no-op for Eigen types and returns an Eigen::Map for std::vector types prior to calling the mat definition:

template <typename Derived>
const auto& as_eigen(const Eigen::MatrixBase<Derived>& v) {
  return v;
}

template <typename T>
const auto as_eigen(const std::vector<T>& v) {
  return Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, 1>>(v.data(),
                                                               v.size());
}

This also means that the templating for the mat definitions needs to be changed from Eigen::Matrix<double, R, C> to Eigen::MatrixBase<Derived> so that it works with Eigen::Map inputs.

After some initial testing, the speed gains are minor:

2019-10-28 20:40:26
Running ./log_sum_exp_arr
Run on (8 X 4000 MHz CPU s)
CPU Caches:
  L1 Data 32K (x4)
  L1 Instruction 32K (x4)
  L2 Unified 256K (x4)
  L3 Unified 8192K (x1)
-----------------------------------------------------
Benchmark              Time           CPU Iterations
-----------------------------------------------------
LogSumExp_Old    8175394 ns    8175190 ns         92
LogSumExp_New    5511970 ns    5511907 ns        124
DotSelf_Old      5136511 ns    5136391 ns        132
DotSelf_New      4497041 ns    4496989 ns        157
Sum_Old          4868005 ns    4867939 ns        139
Sum_New          4608184 ns    4608127 ns        149

Benchmarking code here

The main benefits would be in code reduction for functions with existing arr definitions, as well as making it easier to extend other functions to take std::vector inputs (without requiring additional definitions).

Does this seem like a worthwhile change?

Current Version:

v3.0.0

The way I read this, there is a 30% reduction in execution time for LogSum? 15% for DotSelf and 5 for Sum? I would not call that minor. And even if the execution times would be exactly the same and we "only" gain the code reduction I would be all for it.

So if I understand correct, we can remove the /arr folder completely in this case? If that is the case, I am 100% on board with this idea and please let me know if I can help in any way.

This also means that the templating for the mat definitions needs to be changed from Eigen::Matrix<double, R, C> to Eigen::MatrixBase<Derived> so that it works with Eigen::Map inputs.

This idea has been brewing for some time already anyways for many reasons (most recently sparse matrices I think). This is just another example of why we really need to do this. Some latest discussions on this topic are here: https://discourse.mc-stan.org/t/refactoring-math-with-more-general-signatures/10157/6

What would be the best way forward here? A PR (or multiple PRs) to refactor the /mat signatures and then a PR to add what you are proposing here?

It does, but I think we need to wait until we can squash the prim/may/arr distinction. Right now, we do have the distinction in place and the code base is consistent. We shouldn’t be adding (too many) special cases; it becomes a cognitive burden to know when the rules apply and don’t. I think this is a good idea and we should do it, not just here, but in lots of places!

…

On Mon, Oct 28, 2019 at 9:09 AM Andrew Johnson ***@***.***> wrote: Description As initially discussed in this PR <#1373 (comment)>, there is some code duplication in having the same functions implemented for std::vector and Eigen separately. This could be avoided by using a wrapper which compiles to a no-op for Eigen types and returns an Eigen::Map for std::vector types prior to calling the mat definition: template <typename Derived> const auto& as_eigen(const Eigen::MatrixBase<Derived>& v) { return v; } template <typename T> const auto as_eigen(const std::vector<T>& v) { return Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, 1>>(v.data(), v.size()); } This also means that the templating for the mat definitions needs to be changed from Eigen::Matrix<double, R, C> to Eigen::MatrixBase<Derived> so that it works with Eigen::Map inputs. After some initial testing, the speed gains are minor: 2019-10-28 20:40:26 Running ./log_sum_exp_arr Run on (8 X 4000 MHz CPU s) CPU Caches: L1 Data 32K (x4) L1 Instruction 32K (x4) L2 Unified 256K (x4) L3 Unified 8192K (x1) ----------------------------------------------------- Benchmark Time CPU Iterations ----------------------------------------------------- LogSumExp_Old 8175394 ns 8175190 ns 92 LogSumExp_New 5511970 ns 5511907 ns 124 DotSelf_Old 5136511 ns 5136391 ns 132 DotSelf_New 4497041 ns 4496989 ns 157 Sum_Old 4868005 ns 4867939 ns 139 Sum_New 4608184 ns 4608127 ns 149 Benchmarking code here <https://github.com/andrjohns/perf-math/blob/feature/tbb_test/log_sum_exp_arr.cpp> The main benefits would be in code reduction for functions with existing arr definitions, as well as making it easier to extend other functions to take std::vector inputs (without requiring additional definitions). Does this seem like a worthwhile change? Current Version: v3.0.0 — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#1425?email_source=notifications&email_token=AADH6FYOONN3X5GBE6XX7M3QQ3P67A5CNFSM4JFZ47B2YY3PNVWWK3TUL52HS4DFUVEXG43VMWVGG33NNVSW45C7NFSM4HUYHEXQ>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AADH6FZ5GXL64V37BTSFWELQQ3P67ANCNFSM4JFZ47BQ> .

@rok-cesnovar

Ha. @rok-cesnovar, I think our messages crossed. And you're right -- we can remove the arr implementation completely and just use the mat implementation. We still need to flatten everything, but that's a different issue and we don't need to hold up anything for that.

@SteveBronder

Great! I'll go ahead then.

I think Rok's workflow was a good idea:

Generalise mat signatures
Replace arr functions

After some testing, I've realised that the signatures need to be more general than MatrixBase<Derived>, since they need to be able to take std::vectors as inputs. At the moment, I'm thinking of something like this, using log_sum_exp as an example:

template <typename T, typename = require_eigen_st<std::is_arithmetic,T>>
double log_sum_exp(const T& x) {
...
}

Then, once the we're ready to collapse the arr functions, the template can be updated to:

template <typename T, typename = require_vector_like_st<std::is_arithmetic,T>>
double log_sum_exp(const T& x) {
...
}

PS @SteveBronder your require generics are an absolute lifesaver here!

This looks cool! Questions:

How does this compare to scalar_seq_view in https://github.com/stan-dev/math/blob/develop/stan/math/prim/scal/prob/normal_lpdf.hpp#L62 ? (I guess that one can absorb scalars so maybe I answered myself -- but they seem very similar?)
Can this mechanism absorb row_vectors as well?
Can this be used to write functions that take a std::vector of any of a std::vector, an Eigen vector or an Eigen row vector?

No pressure. Just asking cause these seem related.

Also it probably makes sense to just work this out for one or two representative functions first. No need to try to get rid of everything in array in one go (that'll just make the first review painful for whoever does it).

How does this compare to scalar_seq_view in https://github.com/stan-dev/math/blob/develop/stan/math/prim/scal/prob/normal_lpdf.hpp#L62 ? (I guess that one can absorb scalars so maybe I answered myself -- but they seem very similar?)

Pretty much the same conceptually, in that it's taking different vector types and returning a single type. The main difference is that it needs to return an Eigen type specifically, so that we can directly call their functions (e.g. v.squaredNorm())

Can this mechanism absorb row_vectors as well?

That depends on what you mean, the as_eigen wrapper will return all Eigen objects unchanged, and the require_eigen_st<std::is_arithmetic,T> templating will allow all types under the EigenBase umbrella (which includes row-vectors). So it will work in the sense that the templating will be inclusive of row-vectors, but not in the sense that row-vectors would be mapped to column-vectors. But that's an easy enough inclusion to make as well.

Can this be used to write functions that take a std::vector of any of a std::vector, an Eigen vector or an Eigen row vector?

Tricky one. For a given function it would be easy to add another definition explicitly for nested vectors, but for a single function that handles both vectors and nested vectors, I could look into implementing the as_eigen wrapper in a similar way to vector_seq_view (a la eigen_seq_view). That should make it pretty trivial (touch wood) to add nested vectors to current mat functions as well.

Also it probably makes sense to just work this out for one or two representative functions first. No need to try to get rid of everything in array in one go (that'll just make the first review painful for whoever does it).

Good call!

@andrjohns

The main difference is that it needs to return an Eigen type specifically, so that we can directly call their functions (e.g. v.squaredNorm())

@andrjohns, can you explain this in more detail? I want to make sure I understand. Ok first read, I think what you’ve written here is not correct: if the return type is always an Eigen type, then we can’t mix the arr and mat implementation because this will lose the type information that’s necessary to work with Stan. But maybe what you’ve written isn’t what you were envisioning?

Also, starting with one function and seeing it all the way through is definitely the right way to go.

Sure! It's not that the return type of the function (e.g. log_sum_exp) will always be an Eigen type, it's that the return type of wrapper (i.e.as_eigen) will always be an Eigen type.

Using log_sum_exp as an example. The current definition is:

template <int R, int C>
double log_sum_exp(const Eigen::Matrix<double, R, C>& x) {
  const double max = x.maxCoeff();
  if (!std::isfinite(max)) {
    return max;
  }
  return max + std::log((x.array() - max).exp().sum());
}

What I'm proposing is something like:

template <typename T, typename = require_vector_like_st<std::is_arithmetic,T>>
double log_sum_exp(const T& x) {
  auto v = as_eigen(x);
  const double max = v.maxCoeff();
  if (!std::isfinite(max)) {
    return max;
  }
  return max + std::log((v.array() - max).exp().sum());
}

So the return type (double) doesn't change, it's just that std::vectors are mapped to Eigen types within the function.

You don't have to implement as_eigen(). I think as_column_vector_or_scalar() I implemented when working with GLMs will do exactley what you want.

You don't have to implement as_eigen(). I think as_column_vector_or_scalar() I implemented when working with GLMs will do exactley what you want.

I did have a look at that, the only issue was that it was specifically for vectors, whereas I'm hoping to work with any Eigen object (including Map).

I've put together a rough proof-of-concept here. I've implemented a view framework (similar to vector_seq_view) and updated the mat definition to work with both std::vectors and std::vectors containing other vector types

The templating for the nested views is all explicit at the moment, since I'm still figuring out to navigate the specialisations. But this should give a good overview of what I had in mind.

I like the idea. I think in that case existing function can be reworked to work with any Eigen expression. We don't need two functions that do almost the same thing.

If you want I can do that.

I'm not sure the functions are aiming to do the same thing - since Eigen types are meant to be returned unchanged, not as a column vector. I don't think as_column_vector_or_scalar could be applied to a matrix?

No, it can't be applied to a matrix. But that is not a problem, since a matrix can not be stored in a std::vector. But I see it could be a problem if you want to return row vectors for row vector inputs.

Glad to hear the requires are useful! Are they missing anything conceptually that would be useful? Also for some of the above I think you want the _vt (value_type) instead of _st (scalar_type) since value type only looks one level down. Scalar type looks all the way down so vector of vector would also go in that signature. If as_eigen can work with vector of vectors then _st is the right choice.

... idt we want to think about this in the R way of applying to matrix rows and columns.

Once we get closures, we'll probably want to add those map functions. We already have some specialized cross-product functions.

We definitely need specific-case reductions like log-sum-exp. General ones might also be useful. The algorithms std library has accumulators that Steve's been adding to the Stan math lib code. The accumulators are based on initialization and combination. These will be more useful when get more functional programming in Stan. I'm not going to be the first to say the m-word.

the m-word.

Anyone who mentions those things and forces me to google what they are for the 19th time will be reported for CoC violation

When even the Wikipedia's explanation is a mess, I don't hold out much hope that the 19th time will be the charm.

…

On Nov 11, 2019, at 8:20 PM, Steve Bronder ***@***.***> wrote: the m-word. Anyone who mentions those things and forces me to google what they are for the 19th time will be reported for CoC violation — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

Just a random Idea. Do we for some reason need to separate function implementations from function signatures into a separate struct? I think we could get code that is a bit easier to read if we used generic lambdas and template deduction. So instead of this:

struct log_softmax_fun {
  template <typename T>
  static inline auto fun(const T& v) {
  }
};

template <typename T>
inline auto log_softmax(const T& x) {
  return apply_vector_unary<log_softmax_fun, T>::apply(x);
}

Is there any reason not to do this:

template <typename T>
inline auto log_softmax(const T& x) {
  return apply_vector_unary<T>::apply(x, [](auto& v){
    check_nonzero_size("log_softmax", "v", v);
    return v.array() - log_sum_exp(v);
  });
}

Of cource apply_vector_unary would have to modified a bit to accept functor instead of template type with run member function.

The one advantage to the way it is now is that the functions are all static so that instances never get constructed. That may not be very high overhead compared to everything else. In the testing framework, I've just been using lambdas everywhere.

…

On Nov 12, 2019, at 2:36 AM, Tadej Ciglarič ***@***.***> wrote: Just a random Idea. Do we for some reason need to separate function implementations from function signatures into a separate struct? I think we could get code that is a bit easier to read if we used generic lambdas and template deduction. So instead of this: struct log_softmax_fun { template <typename T> static inline auto fun(const T& v) { } }; template <typename T> inline auto log_softmax(const T& x) { return apply_vector_unary<log_softmax_fun, T>::apply(x); } Is there any reason not to do this: template <typename T> inline auto log_softmax(const T& x) { return apply_vector_unary<T>::apply(x, [](auto& v){ check_nonzero_size("log_softmax", "v", v); return v.array() - log_sum_exp(v); }); } Of cource apply_vector_unary would have to modified a bit to accept functor instead of template type with run member function. — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

They are trivially constructible tho' so idt the overhead would be that much. Could run both through godbolt to see what happens

I think there should not be any overhead in optimized code. But it can't hurt to check.

Uhm... this is a great thread and I am really looking forward to this template magic...given the complexity and size of discussion would it make sense to write a brief design doc and file that in the design repo?

(I really never like when people tell me to write design docs - lots of work - but after having done it things gain in clarity)

Feel free to ignore my comment.

^what does and doesn't need a design doc is still kind of vague, but I think it does fit here. Plus with the new site we can take the approved design doc and add it there so contributors have an easier way to get why we do things a certain way.

I'm on board with the design doc, especially since the scope of this has changed quite a bit from my initial post (in a good way!). I'd just like to clarify the intended direction of this before I write up a doc about the wrong thing.

In my mind, this issue is now concerned with extending the existing prim/mat functions that take Eigen column vectors as inputs to also work with Eigen row vectors and std::vectors, without the need for additional specialisations for each type (and combination) of inputs. There are going to be several different implementations, given the broad taxonomy of vector-based functions we have at the moment:

f(vector) -> scalar
- e.g. log_sum_exp
f(vector) -> vector
- e.g. log_softmax
f(vector) -> matrix
- e.g. diag_matrix
f(vector, scalar) -> vector
- e.g. head
f(vector, vector) -> scalar
- e.g. dot_product
f(vector, vector) -> vector
- e.g. elt_multiply

(Let me know if I've missed a type in there)

In terms of a design, I'm not sure whether it would be possible to have a flexible implementation that would handle more than one of these types of functions, or whether we need a separate one for each (e.g. apply_vector_unary_reduce,apply_vector_unary, and apply_vector_scalar_unary; although it seems like the naming itself might also get out of hand here)

I would very strongly prefer that we do *not* extend matrix operations to std::vector. I know there's been a lot of pressure to make Stan more like R in that it will guess what the shape of your vector is, but I very much want to keep Stan strongly typed. For instance, I do *not* want a std::vector dot product, or an Eigen row vector * std::vector to be a scalar. You should be able to be flexible on the argument types in a meta-programming framework---it's a pretty easy generalization. At worst, you can just take the shape as a template parameter. There are many more signatures than listed when you get into matrix operations. real(row-vector, vector) matrix(vector, row-vector) matrix(matrix, vector) vector(real, vector) .... But those are exactly the ones I don't want to generalize to std::vector. There's also the distribution functions, which are already coded to be flexible about this using views.

I would very strongly prefer that we do not extend matrix operations to std::vector

That sounds good to me, it was going to be a bit of a mess determining whether the to treat them as column or row-vectors.

Navigation Menu

Repository navigation

Use Eigen::Map to replace arr functions #1425

Description

Current Version:

38 remaining items

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Use Eigen::Map to replace arr functions #1425

Description

Description

Current Version:

Activity

rok-cesnovar commented on Oct 28, 2019

syclik commented on Oct 28, 2019

syclik commented on Oct 28, 2019

andrjohns commented on Oct 29, 2019

bbbales2 commented on Oct 29, 2019

andrjohns commented on Oct 29, 2019

syclik commented on Oct 29, 2019

andrjohns commented on Oct 29, 2019

t4c1 commented on Oct 29, 2019

andrjohns commented on Oct 29, 2019

t4c1 commented on Oct 29, 2019

t4c1 commented on Oct 29, 2019

andrjohns commented on Oct 29, 2019

t4c1 commented on Oct 29, 2019

SteveBronder commented on Oct 29, 2019

38 remaining items

bob-carpenter commented on Nov 12, 2019

bob-carpenter commented on Nov 12, 2019

SteveBronder commented on Nov 12, 2019

bob-carpenter commented on Nov 12, 2019

t4c1 commented on Nov 12, 2019

bob-carpenter commented on Nov 13, 2019

SteveBronder commented on Nov 13, 2019

t4c1 commented on Nov 13, 2019

wds15 commented on Nov 13, 2019

SteveBronder commented on Nov 13, 2019

andrjohns commented on Nov 27, 2019

bob-carpenter commented on Nov 27, 2019

andrjohns commented on Nov 27, 2019

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions