Use Eigen::Map to replace arr functions

[andrjohns]

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

[rok-cesnovar]

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?

[syclik]

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> .

[syclik]

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.

[andrjohns]

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!

[bbbales2]

This looks cool! Questions:

1. 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?)

2. Can this mechanism absorb row_vectors as well?

3. 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).

[andrjohns]

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!

[syclik]

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.

[andrjohns]

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.

[t4c1]

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.

[andrjohns]

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.

[t4c1]

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.

[t4c1]

If you want I can do that.

[andrjohns]

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?

[t4c1]

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.

[SteveBronder]

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.