matrices, vectorization, and parallelism

[Matt Hoffman]

Hi all,

Thought I'd revisit this topic. Apologies if it's already been
discussed and I missed it. Looking through (for example)
src/stan/agrad/rev/matrix/multiply.hpp, it seems that we're not
actually calling Eigen to do our matrix multiplies and such. This
obviously precludes our leveraging whatever multicore, SSE, etc.
optimizations that Eigen's got baked in.

My question is: is this an inherent limitation of the current agrad
architecture (which is focused on scalars)? If so, should we consider
revamping the architecture to handle matrices? How much of a colossal
pain would this be? I only ask because the potential payoffs* seem so
large.

Matt

* Potential payoffs that come to mind:
1. We could piggyback on the work that Eigen's done (and will do) to
optimize the hell out of these very important operations.
2. It's fundamentally faster (especially when vectorized instructions
like SSE are available) to work with dense arrays of double values
than with arrays that alternate between packed values and adjoints.
E.g., if we're doing a dot product between two long vectors x and y
the three outputs are dot(x, y), dot(x, y.adj), and dot(x.adj, y). As
it is, all of those dot products has to skip over half of the entries
in each of those arrays, since x[i] is stored right next to x.adj[i].
Those operations should be faster (maybe much faster) if x and y are
stored in separate arrays from x.adj and y.adj.

[Ben Goodrich]

On Sunday, July 7, 2013 9:08:45 PM UTC-4, Matt Hoffman wrote:

My question is: is this an inherent limitation of the current agrad
architecture (which is focused on scalars)? If so, should we consider
revamping the architecture to handle matrices? How much of a colossal
pain would this be? I only ask because the potential payoffs* seem so
large.

I agree that there is no shortage of potential payoffs, and I've always thought that we should be using Eigen as much as possible. My latest thought along the lines you describe would be for Stan to extend the Eigen::Matrix<double,Dynamic,Dynamic> class so that a Stan::Matrix<double,N,K> instance would either include a pointer to a separate Eigen::Matrix<double,N,K> instance that contained the corresponding adjoints or just included a Eigen::Matrix<double,N,K> of adjoints as a member of the Stan::Matrix subclass. But Bob and Marcus would be the ones who would know how much of a pain that would be to integrate into Stan.

Ben

[Matt Hoffman]

Right, I think that's basically what I had in mind—just like we
created var as an extension of double, we could create a type that's
an extension of Eigen::Matrix (mvar?). If we wanted to really go
crazy, we could even replace vars with 1x1 matrices.

> --
> You received this message because you are subscribed to the Google Groups
> "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to stan-dev+u...@googlegroups.com.
> For more options, visit https://groups.google.com/groups/opt_out.
>
>

[Marcus Brubaker]

We could take more advantage of Eigen with the current design if we had a way to allocated large blocks of vari's but it would still be suboptimal.

Changing to something where val and adj were matrices instead of scalars would be a pretty huge change to make but not impossible.  I think the best design would be one where a new class (not inherited from Eigen::Matrix) was created which would store a value and adjoint matrix separately.  In a perfect world this would be able to handle arrays with arbitrary numbers of dimensions and a scalar could become a special 1x1 matrix.

If someone was motivated to do this I would recommend they focus on small test cases first because many of the distributions do some complex things to make them more efficient and it would be easy to get bogged down there early on.

Cheers,

Marcus

[Bob Carpenter]

We have talked about this many times before!

The stumbling block for me is figuring out how
we'd do matrix indexing, which should return an
autodiff var.

On 7/8/13 1:03 AM, Ben Goodrich wrote:
> On Sunday, July 7, 2013 9:08:45 PM UTC-4, Matt Hoffman wrote:
>
> My question is: is this an inherent limitation of the current agrad
> architecture (which is focused on scalars)? If so, should we consider
> revamping the architecture to handle matrices? How much of a colossal
> pain would this be? I only ask because the potential payoffs* seem so
> large.

Yes, they're really big for the matrix multiply part of
the program. It may even be a speedup for derivatives if
we can calculate the partials using matrix ops, too.

> I agree that there is no shortage of potential payoffs, and I've always thought that we should be using Eigen as much as
> possible. My latest thought along the lines you describe would be for Stan to extend the
> Eigen::Matrix<double,Dynamic,Dynamic> class so that a Stan::Matrix<double,N,K> instance would either include a pointer
> to a separate Eigen::Matrix<double,N,K> instance that contained the corresponding adjoints or just included a
> Eigen::Matrix<double,N,K> of adjoints as a member of the Stan::Matrix subclass.

This should be pretty easy for forward-mode, but I
still don't see a clear path to doing it with reverse-mode.

Extending Eigen::Matrix would be problematic because
the results of any operation defined for Eigen::Matrix
would be Eigen::Matrix unless we override everything.

Another option would be to generalize var and vari, which
holds the objects as members rather than extending one and
adding another. It should be pretty easy for fvar.

template <typename T>
struct var {
vari<T>* vi_;
}

template <typename T>
struct vari {
const T val_;
T adj_;
}

Instead of directly using vari, chainable is what's used
on auto-diff stack for reverse mode, so there's some flexibility.

> But Bob and Marcus would be the ones who
> would know how much of a pain that would be to integrate into Stan.

What I can't sort out is how to do indexing. We need to
support the following two methods for indexing, and similarly
for other slices like rows and columns.

var& operator()(int,int);
var& operator()(int);

where the returned references do the right thing.

We could maybe define a specialized new var as a view.
But then we have to figure out how to allow that to
support all the usual operations, so I can write something
like:

exp(sigma[1])

where sigma is a vector. The value for sigma[1] can be
some kind of placeholder, but then we need to define exp()
for that type or somehow figure out how to generalize the definitions.

And once we've done all this, unless we specialize all the
expression templates, we're still back to inefficient compound
operations.

I think the way forward would be to just get it working as a demo
outside of the context of Stan.

- Bob

[Bob Carpenter]

Using 1x1 matrices for scalars would mean generalizing the way we handle
the Eigen types. We'd want to use Eigen::Matrix<T,1,1> if we wanted to
not waste a ton of memory and cause an unacceptable amount of indirection
chasing pointers. It could probably all be done with the kind of abstractions
that the Eigen doc discusses for return types and argument type concepts.

- Bob

On 7/8/13 11:43 AM, Matt Hoffman wrote:
> Right, I think that's basically what I had in mind�just like we

[Matt Hoffman]

Yeah, there's some complexity. I was thinking that one way forward
would be to set up two heaps (one for vals and one for adjs), but
there are tricky things about that too. I agree that the best way to
do it is to try to get it working outside of Stan first. I'd
volunteer, but realistically (and sadly) I won't have time to do this
for months.

On Mon, Jul 8, 2013 at 10:27 AM, Bob Carpenter <ca...@alias-i.com> wrote:
> Using 1x1 matrices for scalars would mean generalizing the way we handle
> the Eigen types. We'd want to use Eigen::Matrix<T,1,1> if we wanted to
> not waste a ton of memory and cause an unacceptable amount of indirection
> chasing pointers. It could probably all be done with the kind of
> abstractions
> that the Eigen doc discusses for return types and argument type concepts.
>
> - Bob
>
>
> On 7/8/13 11:43 AM, Matt Hoffman wrote:
>>

>> Right, I think that's basically what I had in mind—just like we

>> created var as an extension of double, we could create a type that's
>> an extension of Eigen::Matrix (mvar?). If we wanted to really go
>> crazy, we could even replace vars with 1x1 matrices.
>
>
>> On Sun, Jul 7, 2013 at 10:03 PM, Ben Goodrich <goodri...@gmail.com>
>> wrote:
>>>
>>> On Sunday, July 7, 2013 9:08:45 PM UTC-4, Matt Hoffman wrote:
>>>>
>>>>
>>>> My question is: is this an inherent limitation of the current agrad
>>>> architecture (which is focused on scalars)? If so, should we consider
>>>> revamping the architecture to handle matrices? How much of a colossal
>>>> pain would this be? I only ask because the potential payoffs* seem so
>>>> large.
>>>
>>>
>>>
>>> I agree that there is no shortage of potential payoffs, and I've always
>>> thought that we should be using Eigen as much as possible. My latest
>>> thought
>>> along the lines you describe would be for Stan to extend the
>>> Eigen::Matrix<double,Dynamic,Dynamic> class so that a
>>> Stan::Matrix<double,N,K> instance would either include a pointer to a
>>> separate Eigen::Matrix<double,N,K> instance that contained the
>>> corresponding
>>> adjoints or just included a Eigen::Matrix<double,N,K> of adjoints as a
>>> member of the Stan::Matrix subclass. But Bob and Marcus would be the ones
>>> who would know how much of a pain that would be to integrate into Stan.
>
>