Scalar, vector, and matrix views


Is there some reason we can’t replace our basic views with
something like the following. I know we had some writable vector
views, but can those be put into a different class?

  • Bob

// wrap (scalar OR std::vector OR Eigen vector OR Eigen row
// vector) to present sequence of scalars
struct ScalSeqView { };

struct ScalSeqView {
const T& x_;
explicit SeqView(const T& x) : x_(x) { }
const T& operator[](int n) const { return x_; }

struct ScalSeqView<std::vector > {
const vector& xs_;
explicit SeqView(const vector& xs) : xs_(xs) { }
const T& operator[](int n) const { return xs_[n]; }

struct ScalSeqView<Matrix<T, 1, -1> > {
const Matrix<T, 1, -1>& xs_;
explicit SeqView(const Matrix<T, 1, -1>& xs) : xs_(xs) { }
const T& operator[](int n) const { return xs_[n]; }

struct ScalSeqView<Matrix<T, -1, 1> > {
const Matrix<T, -1, 1>& xs_;
explicit SeqView(const Matrix<T, -1, 1>& xs) : xs_(xs) { }
const T& operator[](int n) const { return xs_[n]; }

// wrap Eigen vector or std::vector of Eigen vector
// to produce Eigen vector; also works for Eigen row vectors
// and Eigen matrices
template <typename S, R, C>
struct MatrixSeqView { };

template <typename T, R, C>
struct MatrixSeqView<Matrix<T, R, C> > {
const Matrix<T, R, C>& x_;
explicit MatrixSeqView(const Matrix<T, R, C>& x) { x_(x); }
const Matrix<T, R, C>& operator[](int n) const { return x_; }

template <typename T, R, C>
MatrixSeqView<vector<Matrix<T, R, C> >, R, C>> {
const vector<Matrix<T, R, C> >& xs_;
MatrixSeqView(const vector<Matrix<T, R, C> >& xs) { xs_(xs); }
const Matrix<T, -1, -1>& operator[](int n) const { return xs_[n]; }


This seems like it drops the “broadcasting” behavior we talked about, where no matter what index is passed in the original T was returned, right?

Given a matrix, is it supposed to iterate through it scalar by scalar?


This was a very difficult place to have you start looking
at Stan as it’s poorly written and named code that’s tied up
in some of our most complicated template metaprogramming. Maybe we
can back off and give you something simpler to work on first?
But let me answer.

There are three different wrappers we need:

(a) return scalars [VectorView]
(b) return Eigen vectors
© return Eigen matrices [VectorViewMvt]

I didn’t write most of this code, so I’m not exactly sure how
everything’s plumbed in. But we need the above eventually to
"vectorize" scalars, vectors, and matrices.

I’m not sure what you mean by “given a matrix”. Which of
(a–c)? I’m not sure we ever use VectorView of Eigen matrices
right now, though I think it’s implemented for that. I’m
not sure how we’re doing (b) now, but it may just be ad-hoc
code in some of the distributions with vector parameters (like
multivariate normal, which also has a matrix parameter).

  • Bob


Gotcha. I am happy to work on whatever you guys find useful and you think would make a good starting point. As a straw man proposal, how about I add tests for the currently used and correct behavior in VectorViewMvt and scalar_type_pre so we can close #93 with it, and leave the refactoring and bug in separate issues?

And then do you have any suggestions for next issues? I took a look and mostly have no idea what I’m looking at, but a couple ideas stood out as maybe easy enough:


That’d be great. One thing we want to do is expose
all of our transforms. The problem we’re having with
a lot of these issues is that we’ve already done the
simple ones. A great project would be to expose the
transforms, but that has tentacles all over the code and
we want to combine it with refactoring them and vectorizing
them. And I think Thel may already be working on it, but
I’m not sure. Thel, are you listening?

  • Bob


Would be a good separate thread to come up with projects that are do-able and priorities without diving into all the innards of Stan’s template programming.


@sakrejda, please don’t hijack threads. And to answer your question / non-question: yes. Please feel free to start a new thread.


@thel, there was a question directed at you.


I don’t know who originally wrote it. I’m guessing it’s missing doc. You’re down in the guts of Stan where we haven’t actually touched a lot of it in a long time. We really need to go back and doc everything and test it. Clearly, it’s hard to move forward without knowing what the current stuff does.

I’m in favor of renaming things as Bob wrote in the original post.


I’m watching the thread and just wanted to chime in that nothing I’m doing is likely to obstruct work on this stuff. I’d be happy to get into rewriting some of it as it goes, please tag me in or add me as a watcher on github.


Rob Trangucci wrote off list:

There is a class we added about 2 years ago called container_view that mirrors the VectorView. It was intended to be used in a multi-variate version of OperandsAndPartials. A multi-variate version of OperandsAndPartials will be very useful for writing a version of multi_normal_cholesky that has analytical gradients.

I’d like to leave the container_view in the math library, but it should probably be refactored to mirror the changes you’ll make to the VectorView class.


I didn’t see container_view used anywhere, just defined. Is it used somewhere outside of stan-dev/math?

VectorViewMvt is the multivariate version of VectorView. We’ll be consolidating that and VectorView into these three classes (and maybe a fourth for row_vector_seq_view if necessary):

  • scalar_seq_view: real, real[], vector, row_vector
    const real& operator[](int)

  • vector_seq_view: vector, vector[]
    const vector& operator[](int)

  • matrix_seq_view: matrix, matrix[]
    const matrix& operator[](int)

What Rob’s talking about is here:


I think I may need some help understanding the motivation of Rob’s container_view. What puzzles me is that I don’t see it used anywhere on the develop branch. Did this get put in on spec for some future work?

If we need it, I think we can improve the design to use specializations for the various container types with reference member variables; that way we won’t even need a Map. I don’t see what the matrix version of this gives you over just using a Map directly, though.

Also, is it intended to be mutable? It will be the way it’s designed, since that Map (in the matrix version) will hold a reference to the original memory.

I want to define a separate set of mutable sequence views. We can then replace the builders used in OperandsAndPartials with those.


The intention was for container_view to offer a matrix/array of matrices generalization one of the 3 behaviors in VectorView: allowing a mutable scalar view of an array of doubles. container_view provides a way to return a Vector, Matrix, or a std::vector of vectors/matrices views of arrays of doubles. This would allow us to build a version of operands_and_partials where the gradients could have matrix structure, intended for use something like multi_normal_cholesky. We haven’t built/modified multivariate operands_and_partials yet, but there is a clear need for something like this to simplify the coding of multivariate densities.

VectorViewMvt only does array broadcasting for vectors and matrices because it wasn’t intended to be integrated into a multivariate version of operands_and_partials. Perhaps container_view isn’t the right name for the class, and should be changed to whatever you call the mutable array version of VectorView.


That’s a different application than T_seq_view, the point of which is to allow a T or a container of T to be viewed as a sequence of T.

What I don’t understand is why you need another class. Can’t you just use the Map directly? Or are you thinking OperandsAndPartials will hold pointers to these map holders which will point into their raw memory?


That’s fine; we can rename container_view whatever you’ll be renaming the mutable version of VectorView on lines 11 - 24:

We could do this. It just isn’t how OperandsAndPartials is currently structured. Will you be rewriting to not rely on VectorView? If so, then we can generalize OperandsAndPartials for multivariate operands and bring the container_view functionality into OperandsAndPartials.


Yes, completely getting rid of vector_view. The OperandsAndPartials class has two unusual use cases for it, so we’ll do those last. My goal is to tease apart the usages into separate classes rather than have these things brimming with template parameters and indirection in attempt to perform multiple functions. So that’s not going to be in the first pull request, I don’t think.


Great! So maybe it would make sense to try to generalize OperandsAndPartials after vector_view. We can rename container_view to something more informative for use in the general version. The general OperandsAndPartials will really simplify writing gradients with respect to matrix-variate parameters.


I’m taking a look at OperandsAndPartials generally (as part of looking at how it uses VectorView) and I have two questions:

  1. Are we okay with removing the guards around the throw_if_accessed behavior in VectorView? This would mean when we end up using OperandsAndPartials for constants we don’t have the protection against the programmer accidentally accessing its elements.
  2. Above we talk a bit about generalizing OperandsAndPartials to the multivariate case. Looking at the implementation of normal_lpdf vs multi_normal_lpdf they use a similar number of lines and I don’t think I have enough underlying math to totally understand how a multivariate version would improve the code in the multinormal case. @rtrangucci would you want to sit down and go over what you’re thinking for OperandsAndPartials and maybe explain a bit how it’ll be used? I think that’d help me a bunch in this area, though I could also just leave this part of it for later or someone else.



Yes, absolutely. Bob convinced me a while ago that we don’t have to really defend against all programmers. It was originally set up that way because it’s easy to instantiate things that run over memory without realizing it. There were a couple bugs that we caught with that sort of behavior, but I think we’re good now. Especially if we’re able to doc it properly.

For the simple case, we would be looking at vectorizing all the multivariate distributions by using a view, not necessarily with operands and partials. We would let autodiff work its magic. The next level up would be to build up the operands and partials using double matrices, but that’s a lot more effort.


We have this currently with the VectorViewMvt. I think we should do the multivariate version of OperandsAndPartials, and I do have time to meet with @seantalts this week. Let’s meet on Thursday after the meeting.