vector-based fwd-mode auto-diff

[Bob Carpenter]

Probably the best thing we could do to speed up
Hessian calcs is to use a vector-based forward mode
that computes multiple forward derivatives at once.
Right now, our Hessian calculation works for a K-dimensional
input by running fvar<var> K times, each of which calculates
a gradient of a single partial (a row of the Hessian matrix).

Now:

fvar<T> {
T val_; // value
T d_; // tangent w.r.t. single input
...
}

With vectors:

fvar_vec<T> {
T val_; // value
Matrix<T,Dynamic,1> ds_; // tangent w.r.t. multiple inputs

fvar(T& val, Matrix<T,Dynamic,1>& ds) : val_(val), ds_(ds) { }
}

For instance, exp() will look like this:

template <typename T>
inline
fvar_vec<T>
exp(const fvar_vec<T>& x) {
T val = exp(x.val_);
return fvar_vec<T>(val, val * x.ds_);
}

Presumably the overhead of dealing with vectors will be
less than repeated evals of exp(). In general, the biggest gains
will be with higher dimensionality [fewer allocations per item] and
more complex functions [more expensive computations eliminated], which
is just what we want to speed up.

We are going to have mixed operations (fvar_vec and double), for which
we'd need to specialize the matrix operators like * one more time to
allow clean fvar_vec implementations.

Hessians are now a single evaluation of fvar_vec<var> with a dimension
of ds equal to the number of parameters.

The new C++11 construct-on-collection operators for std::vector<T>
might be more efficient and outweigh the convenience of doing vector
arithmetic in Eigen.

Unfortunately, this would require an entirely new set of fvar_vec
implementations. The good news is that they're all just vectorized
versions of what we already have. And we can use the non-vector versions
for testing.

Of course, we could also unfold the higher-order derivatives in this
situation, but I can't bear to think of the bookkeeping.

- Bob

[Michael Betancourt]

> Unfortunately, this would require an entirely new set of fvar_vec
> implementations. The good news is that they're all just vectorized
> versions of what we already have. And we can use the non-vector versions
> for testing.
>
> Of course, we could also unfold the higher-order derivatives in this
> situation, but I can't bear to think of the bookkeeping.

That's what LaTeX is for! In convincing myself of how to specialize our existing functions
I went through a general derivation that should also hold for vectorized implementations
like this.

Actually, the derivation is pretty clean and it wouldn't be hard to create an operands_and_partials-like
interface where we implement the function value, gradient, Hessian, and third-derivatives and the
specializations can be automated (well, conceptually at least). None of those would have to be
computed twice, although we wouldn't be able to take advantage of common terms unless we wrote
it by hand.