vectorizing transforms

[Bob Carpenter]

[changed subject, moved discussion to stan-dev]

[Sorry, fat-fingered the response out of HTML mailer again, which
is all I can access from where I'm at. I still can't
detect how I'm doing it, or I could be more careful.]

On June 21, 2013 04:14:19 PM EDT, Bob Carpenter <ca...@alias-i.com> wrote:

>
>
>
>>> 1. Vectorize the implementation of the transforms.
>>> As is, they're just iterating, so a vector of positive
>>> parameters (as in that compound K-distribution formulation)
>>> is going to be dog slow.
>>>
>> By vectorize, do you mean checking to make sure the length of the
>> vectors are all the same or one and then use vector
>> operations to carryout the transform?
>
> I meant the efficient thing for derivatives -- define a new vari
> type (for reverse mode) for the Jacobian determinant. This should
> be easy for the simple lower and upper bound constraints, but more
> involved for ordered vectors, unit-length vectors, and simplexes,
> and for covariance and correlation matrices.
>
> The implementation should look something like that of the distributions,
> as there may be vector inputs, but there's only a single return.
> If you can do more efficient implementations for forward-mode, that'd
> be good, too. In fact, if you use OperandsAndPartials, you should be
> able to do both at the same time the way the distributions are being
> done now.
>
> - Bob

[Peter Li]

I don't think I understand what you're saying completely. So for the new
vari type, is it supposed to just take in the jacobian matrix and
compute the determinant?

Also, should I separate the transforms into their own files (i.e.
refactor like we did for agrad/fvar/math)?


Peter

[Bob Carpenter]

On 6/24/13 11:04 AM, Peter Li wrote:
> I don't think I understand what you're saying completely. So for the new vari type, is it supposed to just take in the
> jacobian matrix and compute the determinant?

However you have to implement it, there should be a single
vari for the log Jacobian determinant that gets added to the
log probability parameter being passed in.

The point is to make the entire transform stage of Stan faster.

For example, suppose we have:

parameters {
real<lower=L> x[N];
...
}

The transform to constrained space is

x[n] |--> exp(x[n]) + L

So the Jacobian is diagonal so the log determinant is
just sum(x). So you won't even need a new type in
this case.

The lower-bound L may be a variable, but even if it is,
dL/dx = 0 because the bound L can't depend on x because
of variable declaration ordering constraints in Stan (can't
use a variable before it's fully declared).

So why don't you start with this case.

The full set of transforms we use and their Jacobian determinants
are defined in the manual.

> Also, should I separate the transforms into their own files (i.e. refactor like we did for agrad/fvar/math)?

That would be great.

- Bob

[Peter Li]

I understand how the jacobian plays into this, but I'm still confused
about defining a vari for it and how to use it. Right now, I have
something like

template <typename T, typename TL>
inline
typename return_type<T,TL>::type>::type
lb_constrain(const T& x, const TL& lb, T& lp) {

if (!(stan::length(x)))
return 0.0;

double res(0.0);

agrad::OperandsAndPartials<T> operands_and_partials(x);
VectorView<const T> x_vec(x);
size_t N = max_size(x);
for (size_t n = 0; n < N; n++) {
const double x_dbl = value_of(x_vec[n]);

lp += x_dbl;

res += exp(x_dbl) + lb;

if (!is_constant_struct<T>::value)
operands_and_partials.d_x1[n] += exp(x_dbl);
}

return operands_and_partials.to_var(res);
}

But I'm not sure how I should be treating res--should it be a vector at
least type T to match x?


Peter

[Bob Carpenter]

Now I see the issue.

The model:

parameters {
real<lower=0> y[10];
}
model {
y ~ normal(0,1);
}

translates the constraints to:

vector<T__> y;
size_t dim_y_0__ = 10;
y.reserve(dim_y_0__);
for (size_t k_0__ = 0; k_0__ < dim_y_0__; ++k_0__) {
y.push_back(in__.scalar_lb_constrain(0,lp__));
}

So in order to do what I want, we need to change the generated
code to something like:

vector<T__> y = in__.scalar_lb_constrain_array(dim_y_0__,0,lp__);

Then scalar_lb_constrain_array() would generate a single vari
instance to add to lp__ for the Jacobian.

So either you're going to have to dive into the generator code
in stan/gm, or wait until I get to it.

- Bob

[Ben Goodrich]

On Monday, June 24, 2013 5:45:34 PM UTC-4, Bob Carpenter wrote:

So either you're going to have to dive into the generator code
in stan/gm, or wait until I get to it.

While we are at it, perhaps now would be a good time to think about the design for vectors of support? As in, hypothetically,

vector<lower=lb,upper=ub>[K] beta;

where lb and ub are K-vectors (perhaps parameterized) of lower and upper bounds.

Ben