C++ 'slices' or 'vectorview' read-only non-allocating view of a std::vector

Developers
1

Hey all,

Today in the meeting @wds15 was mentioning that he’s working on a reduce implementation that needs to know what subset of the data to operate on, and Breck mentioned the idea of a slice or view, which is common in other languages. Unfortunately in C++ the community seems to have decided that the real abstraction you want here is a pair of iterators, but that’s not that easy to deal with in our situation where we want e.g. repeated random access. However, with Eigen we could use their Map to create a view of an Eigen Vector:
https://eigen.tuxfamily.org/dox/group__TutorialMapClass.html

Basically we’d declare a map pointing to the underlying data in a std::vector or eigen vector and then we could use it just like an Eigen Vector. I think this would satisfy all our use-cases, though we’d have to template or write functions expecting Eigen::DenseBase. Thoughts on this? Anyone else have other ideas?

2

I think gsl::span (or std::span in C++20) may be what you are looking for. See here: https://github.com/isocpp/CppCoreGuidelines/blob/master/docs/gsl-intro.md

1 Like
3

Yes! Glad to see the C++ community has come around on this issue. Is there no way to get it without including Microsoft’s (MIT licensed) GSL implementation (or waiting for a 2020 compiler, hah).

4

There is GSL-lite: https://github.com/martinmoene/gsl-lite

5

Oh nice, that looks relatively easy to include. Any other contenders? @wds15 what do you think about these two options?

6

We (Bob, myself, and I believe @jpritikin and perhaps a handful of others) evaluated a long time ago whether we could switch over to using Eigen::DenseBase. Although the doc makes it seem like it’s a drop-in replacement for Eigen::Matrix, unfortunately it isn’t. If you start implementing, you’ll find that numerical results don’t line up. I believe it was a design decision in Eigen for certain things to behave differently (there should be a post on Eigen’s mailing list from me with a couple interactions), and so in order to do this correctly, we’d have to carry around an implementation of both Eigen::DenseBase and Eigen::Matrix.

Correction… we evaluated Eigen::MatrixBase, not Eigen::DenseBase. At the bottom, I’ve attached an email that I sent out (which I can’t find on the Eigen mailing list); I think I received a comment back at some point.

But… if I’m to guess at the purpose, do we actually need these to be in Eigen data structures for linear algebra? Or are we just looking to access slices? If it’s the later, then we could do whatever would make the implementation work… pointers and iterators do seem to be the natural way to do this.

From: Daniel Lee
Date: March 12, 2015 at 4:44:36 PM EDT
To: eigen@lists.tuxfamily.org
Subject: help using foo(const MatrixBase&) with (u * v)

We’re trying to follow the advice on this page:

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

and running into trouble passing the expression template (u * v)
to a function declared as suggested above.

It seems to work fine for transpose, for Map, for MatrixXd, but
here’s a minimal example of the problem we’re running into with
products:

#include <Eigen/Dense>

template <typename Derived>
double foo(const Eigen::MatrixBase<Derived>& a, int m, int n) {
  return a(m,n);
}

int main() {
  Eigen::MatrixXd u(2,2);
  u << 1, 2, 3, 4;

  Eigen::MatrixXd v(2,3);
  v << 1, 2, 3, 4, 5, 6;

  foo(u*v, 0, 0);

  return 0;
}

~/temp2$ clang++ -I ~/stan/lib/eigen_3.2.2 eigen-test.cpp

~/temp2$ ./a.out
Assertion failed: (this->rows() == 1 && this->cols() == 1), function coeff, file /Users/carp/stan/lib/eigen_3.2.2/Eigen/src/Core/ProductBase.h, line 141.
Abort trap: 6

~/temp2$ clang++ --version
Apple LLVM version 6.0 (clang-600.0.56) (based on LLVM 3.5svn)
Target: x86_64-apple-darwin14.1.0
Thread model: posix

Here’s the relevant part of Eigen/src/Core/ProductBase.h:

    typename Base::CoeffReturnType coeff(Index row, Index col) const
    {
#ifdef EIGEN2_SUPPORT
      return lhs().row(row).cwiseProduct(rhs().col(col).transpose()).sum();
#else
      EIGEN_STATIC_ASSERT_SIZE_1x1(Derived)
      eigen_assert(this->rows() == 1 && this->cols() == 1);   // *** LINE 141 ***
      Matrix<Scalar,1,1> result = *this;
      return result.coeff(row,col);
#endif
    }

If I change the #ifdef to an #ifndef, everything works as expected.
Are we somehow supposed to be setting the EIGEN2_SUPPORT property
ourselves? If I do that in the compilation, it works as expected, but
I get a deprecation warning:

~/temp2$ clang++ -I ~/stan/lib/eigen_3.2.2 -DEIGEN2_SUPPORT eigen-test.cpp
In file included from eigen-test.cpp:1:
In file included from /Users/carp/stan/lib/eigen_3.2.2/Eigen/Dense:1:
In file included from /Users/carp/stan/lib/eigen_3.2.2/Eigen/Core:373:
/Users/carp/stan/lib/eigen_3.2.2/Eigen/Eigen2Support:20:2: warning: "Eigen2 support is
      deprecated in Eigen 3.2.x and it will be removed in Eigen 3.3. (Define
      EIGEN_NO_EIGEN2_DEPRECATED_WARNING to disable this warning)" [-W#warnings]
#warning "Eigen2 support is deprecated in Eigen 3.2.x and it will be removed in Eigen 3...
 ^
1 warning generated.
~/temp2$ ./a.out
~/temp2$

Daniel

1 Like
7

This example does not assert when I run it.

2 Likes
8

Cool. Maybe it’s been fixed!

9

The Eigen route sounds reasonable, but as Daniel mentions it was evaluated a while ago. Maybe that assessment needs to be done again since the test is now succeeding - probably other places could benefit from going with densebase.

The span concept looks quite nice I have to say. I am not sure what is better right now. Both can work, but maybe span is even more general and useful as a concept?

10

I like the span concept a lot too! It may suffer the same issue, namely that our specializations are often on std::vector or Eigen::Matrix specifically. I would love to try a dumb find-and-replace of Eigen::Matrix parameters with Eigen::DenseBase (or MatrixBase, I forget which) and see what happens.

@wds15 what were you planning with reduce + start, end indices?

11

I wanted to go with something which is very close to the interface of the TBB reduce:

https://www.threadingbuildingblocks.org/docs/help/index.htm#tbb_userguide/parallel_reduce.html

I would start with the first signature of the two they have. That one looks simpler to start with while the second may offer performance benefits. The TBB then evaluates the user function over a range which is represented by start and end index. I would envision to pass a c++ functor into the reduce. The c++ functor would be called with the start and end index. The data to work over would be internal to the functor. At least this is how I think it can be implemented relatively easy and it would still be very flexible.

This is where I would start given my current thinking. My goal is to rather quickly arrive at some working version in order to start benchmark it. I really wanna see the performance jump which I am expecting to see.

Interestingly, the TBB provides a parallel_reduce and a parallel_deterministic_reduce. It is easy to swap one for the other, but this is an option which should probably give users eventually.

12

Very cool, I’m a big fan of leveraging TBB primitives as much as possible.

I guess my question was more about how we were going to pass a pointer and index range to the underlying math functions that we’re mapping or reducing over. Like if people are still using e.g. normal_lpdf, that is written in terms of either std::vector or Eigen::Vector. If I’m thinking about this right, somewhere we have to translate from that efficient range + pointer into something the Math library is equipped to deal with.

Very relatedly, @Bob_Carpenter and I have been talking a lot about the best way to pre-compile the math library (mostly to help with model compilation times), and Bob is a fan of the pointer+len pair, which we should be able to translate into from a pointer and start, end range by adding the offset to the pointer (eep) and using that as the new pointer. I think this would be roughly equivalent efficiency-wise as writing everything in terms of Eigen::Maps (assuming we use them in a good way), but the latter gives us a much richer interface, so I might prefer that by default.

13

I think we’d create a view of a std::vector. I don’t see why we’d need a view of an Eigen::Vector unless we’re doing some kind of slicing, but then Eigen already has those expression templates built-in.

I’m always confused by their docs on writing functions for Eigen.

https://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html

I think there’s some issues with just using auto and also issues with block (the slicing function for matrices/vectors) being incompatible in some situations.

14

Sorry—somehow responded before seeing the more detailed response from @syclik. I think we also ran into trouble with block operations in other places, but there’s a good chance I’m misremembering.

15

Oops, yes, I meant a view that looked like an Eigen vector (from any chunk of memory, but yeah probably the one underlying a std::vector).