elementwise (vectorized) functions

[syclik]

From @bob-carpenter on May 27, 2015 19:3

All of our unary scalar functions should be vectorized to work elementwise with the same names on any container, including

• double, int
• vector, row-vector, matrix
• arrays of any of these

For binary or larger operations, this should work so that either the dimensions match of the two arguments or one argument is a scalar that is broadcast. For example, pow(matrix,matrix) or
pow(matrix,double) or pow(double,matrix) should all return a matrix.

An even more ambitious goal is to do this kind of vectorization and broadcasting for arguments that aren't simple scalars. So we'd have matrix[] inverse(matrix[]) apply inverse elementwise to the matrices in the array.

The output will always be the same shape as the input and the operations applied elementwise.

Seth Flaxman pointed out we want to be careful not to conflict with matrix exponentials, matrix power, etc. So we'll need a separate matrix matrix_pow(matrix,int) operation.

This issue can be approached in parts. Create a smaller issue and link into here if you just want to implement some of the functions, but leave this one standing until it's complete in the language.

Ideally, we'd have some template programs to do this vectorization given a functor (or perhaps a static class function). Please ping me for design help with this before taking it on.

Copied from original issue: stan-dev/stan#1464

[syclik]

From @bgoodri on May 27, 2015 20:23

According to
http://eigen.tuxfamily.org/dox/group__QuickRefPage.html
"It is also very simple to apply any user defined function foo using DenseBase::unaryExpr together with std::ptr_fun:

mat1.unaryExpr(std::ptr_fun(foo))

That just leaves std::vector containers.

[syclik]

From @bob-carpenter on May 27, 2015 23:36

The Eigen method unaryExpr is exactly the kind of thing
I'm talking about --- we'd want to do the same thing for
std::vector, but make it recursive so it deals with arrays
of dimensionality.

The ptr_fun functor wrapper is deprecated in C++11, but
it won't work with scalar types anyway --- but that's not
the critical bit of Ben's example.

I was thinking of using a static class function specified
by template parameter:

  template <typename F, typename T, int R, int C>
  inline Matrix<T, R, C>
  unary(const Matrix<T, R, C>& x) {
    Matrix<T, R, C> fx(x.rows(), x.cols());
    for (int j = 0; j < x.rows(); ++j)
      for (int i = 0; i < x.cols(); ++i)
        fx(i,j) = F::apply(x(i,j));
    return fx;
  }

In this scheme, there's no functor passed around or constructed---instead
the class F must implement the static function:

  static T apply(const T& x);

It saves on building constructors. It's less good if we need to put
things together using standard functor tools or have an automatic
way to build a functor out of parts (like various function pointers).

Eigen's might be a lot faster if they can somehow construct and return
a matrix all at once, though I think the compilers are getting better.
I'll have to see how they implemented it for pointers and also do a speed
comparison.

Then we need base cases like this:

  template <typename F>
  inline double unary(double x) {
    return F::apply(x);
  }

and a recursive case like:

  template <typename F, typename T>
  inline std::vector<T> unary(const std::vector<T>& x) {
    std::vector<T> fx;
    for (size_t i = 0; i < x.size(); ++i)
      fx[i] = unary<F>(x[i]);
    return fx;
  }

And then we'd need to create appropriate class wrappers, such as:

  class cos_fun {
    template <typename T>
    T cos(const T& x) {
      using std::cos;
      return cos(x);
    }  
  }

And then we define:

  template <typename T>
  T cos(const T& x) {
    return apply<cos_fun>(x);
  }

and we're done We have all the vectorized forms of cos
we could ever want. We just have to be careful that we don't
get ambiguous calls, but I don't think that'll be a problem for
this overload because all of our definitions are going to be more
specific. If not, we'll be in for some enable_if wrangling that
will restrict T to std::vector or matrix types.

And these all play nice and generally with our variants of cos.

• Bob

On May 27, 2015, at 4:23 PM, bgoodri [email protected] wrote:

According to
http://eigen.tuxfamily.org/dox/group__QuickRefPage.html
"It is also very simple to apply any user defined function foo using DenseBase::unaryExpr together with std::ptr_fun:

mat1.unaryExpr(std::ptr_fun(foo))

That just leaves std::vector containers.

—
Reply to this email directly or view it on GitHub.

[bgoodri]

This is duplicative of #202