Fwd: [stan-users] random effects models with known covariance matrix

[Ben Goodrich]

Do we not expose the dot_self() function to the parser? I was looking for it when I wrote this post on stan-users, but it seems to be C++ only.

Ben

On Tuesday, September 11, 2012 11:04:58 AM UTC-4, Marcus wrote:

With this latest version you should be able to write:

v ~ normal(0,sqrt(sigmasq_u));

And get almost exactly the same performance without having to muck
about with the log probability directly.

Or, alternatively, if you really care about getting the most speed
possible leave the lp__ statements as they are and  replace the loop

>  sum_of_squares <- 0.0; // or sum(pow(v, 2)) may / should work
>  for(i in 1:N) {
>    sum_of_squares <- sum_of_squares + pow(v[i], 2);
>
>  }

with

sum_of_squares <- dot_product(v,v);

which will be a fair bit faster than the loop, particularly for large N.

Cheers,
Marcus

On Tue, Sep 11, 2012 at 9:44 AM, Ben Goodrich <goodri...@gmail.com> wrote:
> If you want a little more juice, you can exploit the fact that
>
> u' * A_inv * u = u' * ( LA_inv * LA_inv' ) * u = (u' * LA_inv) * (LA_inv' *
> u) = v' v
>
> where LA_inv is the Cholesky factor of A_inv and v = LA_inv' * u . Thus, the
> lower triangle of LA_inv' is all zeros. Stan is smart enough to ignore
> multiplication and addition involving zeros and also knows not to carry
> through a bunch of constant autodiffs.
>
> So, the blocks of your .stan file would look like
>
> transformed data {
>   matrix[N,N] A_inv; A_inv <- inverse(A);
>   t_LA_inv <- transpose(cholesky(A_inv));
> }
>
> model {
>  vector[N] v;
>  real sum_of_squares;
>  v <- t_LA_inv * u;
>
>  sum_of_squares <- 0.0; // or sum(pow(v, 2)) may / should work
>  for(i in 1:N) {
>    sum_of_squares <- sum_of_squares + pow(v[i], 2);
>
>  }
>
>  lp__ <- lp__ - 0.5 * N * log(sigmasq_u); // non-constant part of
> log(det(sigmasq_u * A)^-0.5)
>
>  lp__ <- lp__ - sum_of_squares / (2 * sigmasq_u); // log of kernel of
> multivariate normal
>
>  // rest of your likelihood and priors
> }
>
>
> Ben
>
>
> On Tuesday, September 11, 2012 4:32:25 AM UTC-4, Frank wrote:
>>
>> The custom made multivariate normal sampler suggested by Ben really did
>> the trick! Thanks allot.
>>
>> Memory is not an issue anymore, even with N = 1000.
>>
>> I think that Stan is now at least not slower than Jags/MCMCglmm ... on a
>> per iteration base. Considering that Stan is supposed to need much less
>> total number of iterations for a given ESS, it may have quiet an edge.
>>
>> Best regards,
>>
>> Frank
>>
>>
>>
>>
>> Am Dienstag, 11. September 2012 02:30:41 UTC+2 schrieb Ben Goodrich:
>>>
>>> On Monday, September 10, 2012 12:46:53 PM UTC-4, Marcus wrote:
>>>>
>>>> One thing you can try is to use multi_normal_cholesky
>>>>
>>>> I.E., add
>>>>
>>>> transformed data {
>>>> matrix[N,N] A_chol;
>>>> A_chol <- cholesky_decompose(A);
>>>> }
>>>>
>>>> then in the model
>>>>
>>>> u ~ multi_normal_cholesky(0,sqrt(sigmasq_u)*A_chol);
>>>>
>>>> Also, you should be able to use the vectorized version of the normal
>>>> distribution. I.E., instead of the for loop over ys write
>>>>
>>>> y ~ normal(beta0 + u, sqrt(sig2_e));
>>>
>>> Another thing to try (especially if RAM is the binding constraint) is to
>>> do the multivariate normal part yourself in the .stan file with
>>>
>>> lp__ <- lp__ + whatever;
>>>
>>> statements. If A_chol is fixed but sigmasq_u is a parameter, then
>>> sqrt(sigmasq_u)*A_chol is a parameter and multi_normal_cholesky() is doing a
>>> giant triangular solve every time the function gets called which also spews
>>> out a lot of autodiff calculations. If you invert A once at the outset (to
>>> get A_inv) and be careful not to divide by sigmasq_u too soon, it will go
>>> faster. Something like
>>>
>>> lp__ <- lp__ - 0.5 * N * log(sigmasq_u); // non-constant part of
>>> log(det(sigmasq_u * A)^-0.5)
>>> lp__ <- lp__ - (u' * A_inv * u) / ( 2 * sigmasq_u); // log of kernel of
>>> multivariate normal
>>>
>>> Ben
>>>
>

[Marcus Brubaker]

We should expose it if we don't but I think it may be slower than
dot_product(v,v) right now, as it doesn't use a special vari like
dot_product.

Cheers,
Marcus

> --
>
>

[Bob Carpenter]

We probably also want this for Cholesky factors:

multiply_self_transpose(L) = L * L'

and this for general multivariate distros:

multiply_wrap_transpose(y,Sigma) = y' * Sigma * y

or even

normal_thingy(mu,y,SigmaInv) = (y - mu)' * SigmaInv * (y - mu)

implemented with special vari. But then we should just
vectorize all of multi_normal.

- Bob

[Marcus Brubaker]

On Tue, Sep 11, 2012 at 3:14 PM, Bob Carpenter <ca...@alias-i.com> wrote:
> We probably also want this for Cholesky factors:
>
> multiply_self_transpose(L) = L * L'
>
> and this for general multivariate distros:
>
> multiply_wrap_transpose(y,Sigma) = y' * Sigma * y
>
> or even
>
> normal_thingy(mu,y,SigmaInv) = (y - mu)' * SigmaInv * (y - mu)
>
> implemented with special vari. But then we should just
> vectorize all of multi_normal.

A special vari to handle a Mahalanobis norm/distance is probably a
good idea, separate from the vectorized multi_normal. It will be
handy in efficiently implementing any elliptical distribution, like
the multivariate student t distribution. That said, I don't think it
will be much faster than the existing multi_normal code for the
non-vectorized case.

Also, I just checked and found that there is a dot_self vari. I'm not
sure who did that, it might have even been me...

Cheers,
Marcus

> --
>
>

[Bob Carpenter]

On 9/11/12 3:27 PM, Marcus Brubaker wrote:
> On Tue, Sep 11, 2012 at 3:14 PM, Bob Carpenter <ca...@alias-i.com> wrote:
>> We probably also want this for Cholesky factors:
>>
>> multiply_self_transpose(L) = L * L'
>>
>> and this for general multivariate distros:
>>
>> multiply_wrap_transpose(y,Sigma) = y' * Sigma * y
>>
>> or even
>>
>> normal_thingy(mu,y,SigmaInv) = (y - mu)' * SigmaInv * (y - mu)
>>
>> implemented with special vari. But then we should just
>> vectorize all of multi_normal.
>
> A special vari to handle a Mahalanobis norm/distance is probably a
> good idea, separate from the vectorized multi_normal. It will be
> handy in efficiently implementing any elliptical distribution, like
> the multivariate student t distribution. That said, I don't think it
> will be much faster than the existing multi_normal code for the
> non-vectorized case.
>
> Also, I just checked and found that there is a dot_self vari. I'm not
> sure who did that, it might have even been me...

I may have put it in there. I don't even care
enough to figure out "git blame" :-)

Given that we had to run another patch, I'll add it
to the function_signatures to allow it to be used
in models in 1.0.1.

- Bob

[Ben Goodrich]

On Tuesday, September 11, 2012 3:14:08 PM UTC-4, Bob Carpenter wrote:

But then we should just vectorize all of multi_normal.

Actually, I apparently did vectorize multi_normal back in May. I meant to redo it in the new style of normal(), but discovered that I had already pushed it accidentally. It has tests, but I don't think it is exposed to the parser.

Ben

[Andrew Gelman]

This comment (below) reminds me . . . once we have "forward simulation" working, we can talk about implementing the automatic check of inferences as in the Cook/Gelman/Rubin paper.

On Sep 11, 2012, at 7:19 AM, Frank wrote:

> I replied to the wrong post. So here is the code. From simulations it seems that the results produced are very sensible.
>

[Bob Carpenter]

Re: the subject

By "self-cleaning oven", Andrew means testing.

On 9/11/12 5:02 PM, Andrew Gelman wrote:
> This comment (below) reminds me . . . once we have "forward simulation" working, we can talk about implementing the automatic check of inferences as in the Cook/Gelman/Rubin paper.

Yes, that was the earlier context of that post.
When we have RNGs implemented in the models, it
should be easy to move missing data from being
specified as parameters to being specified as
generated quantities.

With subroutines, it should even be easy to share
the relevant code, but that's on the to-do queue
after RNGs.

- Bob

[Bob Carpenter]

You added a version with a vectorized y in
multi_normal(y|mu,Sigma) by letting y be a matrix.

The "new" style of vectorization is much more involved
and supports arbitrary subsets of the arguments being
vectorized and/or agrad vars.

Multivariate normals may be a problem, though, in that
we won't be able to easily compute the gradients (or maybe
they're easier than I suspect?).

- Bob

[Ben Goodrich]

On Tuesday, September 11, 2012 5:22:49 PM UTC-4, Bob Carpenter wrote:

Multivariate normals may be a problem, though, in that
we won't be able to easily compute the gradients (or maybe
they're easier than I suspect?).

The gradient wrt mu is not bad. The gradient wrt the primitive parameters underlying Sigma is more difficult, especially considering that multi_normal() wouldn't know how Sigma was constructed.

Ben

[Marcus Brubaker]

We could vectorize multi_normal_cholesky first and write the
vectorized mult_normal in terms of that. The derivatives wrt L should
be more straight forward in the cholesky form. Both cases would need
to rely on the identities

d inv(A) = - inv(A) * dA * inv(A)

and

d log(det(A)) = Tr(inv(A) dA)

for matrix A. Given those both forms should be relatively straight forward.

Also, on a side note about vectorizing things. For large N

for (i = 1:N) {
a += foo(b[i])
}

is *much* slower than

vector[N] tmp;

for (i = 1:N) {
tmp[i] = foo(b[i]);
}

a += sum(tmp);

This is true regardless of the complexity of foo(). The extra
allocation and assignment operations are cheap compared to the N extra
vari instantiations implied by the repeated += operations. This is
probably not true with a is not a variable.

This applies equally well for a model or for a vectorized
implementation of a function.

Cheers,
Marcus

[Ben Goodrich]

On Tuesday, September 11, 2012 7:10:59 PM UTC-4, Marcus Brubaker wrote:

> The gradient wrt mu is not bad. The gradient wrt the primitive parameters
> underlying Sigma is more difficult, especially considering that
> multi_normal() wouldn't know how Sigma was constructed.
 

We could vectorize multi_normal_cholesky first and write the
vectorized mult_normal in terms of that.  The derivatives wrt L should
be more straight forward in the cholesky form.  Both cases would need
to rely on the identities

d inv(A) = - inv(A) * dA * inv(A)

and

d log(det(A)) = Tr(inv(A) dA)

for matrix A.  Given those both forms should be relatively straight forward.

That is all true. I just didn't know how to make dA or dL fall from the derivative tree.

Ben
 

[Ben Goodrich]

I can run make test-unit with g++-4.6 without hitting this error, so it must not be getting tested.

On Wednesday, September 12, 2012 9:50:43 PM UTC-4, Ben wrote:

~>g++ --version

g++ (Ubuntu/Linaro 4.6.3-1ubuntu5) 4.6.3

Copyright (C) 2011 Free Software Foundation, Inc.

This is free software; see the source for copying conditions.  There is NO

warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

On Wed, Sep 12, 2012 at 7:28 PM, Ben Goodrich <goodri...@gmail.com> wrote:

On Wednesday, September 12, 2012 5:32:25 PM UTC-4, Ben wrote:

r<T, std::allocator<_CharT> >&) [with T = double]':
chains.cpp:150:   instantiated from here
../inst/include/stansrc/stan/prob/autocovariance.hpp:63: error:
return-statement with a value, in function returning 'void'
make: *** [chains.o] Error 1
ERROR: compilation failed for package 'rstan'
* removing '/home/domingue/R/x86_64-redhat-linux-gnu-library/2.15/rstan'
Warning in install.packages("rstan", type = "source") :
  installation of package 'rstan' had non-zero exit status

BTW, can you tell us the output of

g++ --version

?

Thanks,
Ben

[Jiqiang Guo]

But it seems fine that we can write as follows although the `return` in testf is redundant. 

void void_fun() {

}

void testf() {

  return void_fun();  

} 

This happens in src/stan/prob/autocovariance.hpp:63, which seems to be the issue that the log reports. 

And it is a bit weird that the log of installing RStan says his system is redhat, but the g++ --version says it is ubuntu. 

Jiqiang

--
 
 

[Ben Goodrich]

He later said that it was g++-4.4.x on RedHat. But, test-unit still works for me. Anyway, the end of autocovariance.hpp looks messed up.

    template <typename T>
    void autocovariance(const std::vector<T>& y,
                        std::vector<T>& acov,
                        Eigen::FFT<T>& fft) {

      stan::prob::autocorrelation(y, acov, fft);

      T var = stan::math::variance(y) * (y.size()-1) / y.size();
      for (size_t i = 0; i < y.size(); i++) {
        acov[i] *= var;
      }
    }

    template <typename T>
    void autocovariance(const std::vector<T>& y,
                        std::vector<T>& acov) {
      Eigen::FFT<T> fft;
      return autocovariance(y,acov,fft);
    }

Since the second version doesn't have a Eigen::FFT argument like the first version does, I think it shouldn't have return type void.

Ben

[Bob Carpenter]

n 9/12/12 11:07 PM, Ben Goodrich wrote:
> He later said that it was g++-4.4.x on RedHat. But, test-unit still works for me. Anyway, the end of autocovariance.hpp
> looks messed up.

Indeed it is. Thanks for the heads up.

See below for the fix.

>
> template <typename T>
> void autocovariance(const std::vector<T>& y,
> std::vector<T>& acov,
> Eigen::FFT<T>& fft) {
>
> stan::prob::autocorrelation(y, acov, fft);
>
> T var = stan::math::variance(y) * (y.size()-1) / y.size();
> for (size_t i = 0; i < y.size(); i++) {
> acov[i] *= var;
> }
> }
>
> template <typename T>
> void autocovariance(const std::vector<T>& y,
> std::vector<T>& acov) {
> Eigen::FFT<T> fft;
> return autocovariance(y,acov,fft);
> }
>
> Since the second version doesn't have a Eigen::FFT argument like the first version does, I think it shouldn't have
> return type void.

I just patched it by removing "return".

The return type of void was right, but there shouldn't
have been a return statement. All this function does is
allocate the FFT object and then pass it to three-argument
version that fills in the acov vector passed by reference.
The three-argument version is there to allow you to reuse
the FFT object across multiple calls (it does some internal
allocation that I think can be reused, though I haven't
done speed tests).

- Bob

[Bob Carpenter]

On 9/11/12 7:10 PM, Marcus Brubaker wrote:

> Also, on a side note about vectorizing things. For large N
>
> for (i = 1:N) {
> a += foo(b[i])
> }
>
> is *much* slower than
>
> vector[N] tmp;
>
> for (i = 1:N) {
> tmp[i] = foo(b[i]);
> }
>
> a += sum(tmp);
>

> .. The extra

> allocation and assignment operations are cheap compared to the N extra
> vari instantiations implied by the repeated += operations. This is
> probably not true with a is not a variable.

Right -- the speed's from not allocating vari objects and
having to visit them and call their virtual functions during
derivative propagation.

I just dropped a note to that effect in the manual section
on vectorization.

I've thought about this issue before (though not so far
as to vectorize it with sum()), and I think we can
probably use an expression template for the result
of += to lazily drop a sequence of += into a vector sum.

By concentrating on +=, we wouldn't need to
do tons of overloads and conversions as we would with
a more full template expression approach.

sum_expression& operator+=(sum_expression&, var&);

sum_expression& operator+=(var&, var&);

and then we'd need a copy constructor for var
to create the vectorized sum:

var(sum_expression);

> This applies equally well for a model or for a vectorized
> implementation of a function.

Exactly the same principle.

- Bob

[Marcus Brubaker]

I'm not sure that would work because the return type of operator+= is
generally ignored unless someone does (a += b) += c or d = (a += b) or
something similarly obtuse and strange. In fact, += has to be defined
as a member function:

class var {
...
var &operator+=(const var &b);
...
}

so it can't change the type of the LHS of a += expression. Rather, we
would need a sum vari which can have a variable size and have things
added to it. I think this would break the rules of a vari as it's
value cannot change once its instantiated.

I think the key here is to just write smart C++ code taking this into
account. In terms of all the "lp__ +=" operations in a model, we
could just as easily push the RHS of those statements into a vector
and then call a single sum at the end. (I'm not sure if I've
impemented the sum vari operation to handle a vector<var> type, but
that shouldn't be too hard.

Cheers,
Marcus

[Bob Carpenter]

You can define it outside of the class, but in some sense
it still behaves like a member function. I really find
this list of signatures useful:

http://en.wikipedia.org/wiki/Operators_in_C_and_C%2B%2B

>
> class var {
> ...
> var &operator+=(const var &b);
> ...
> }
>
> so it can't change the type of the LHS of a += expression.

Good point -- we can't have functions that vary in return
type but not in argument types. So we'd have to make
the lp__ accumulator of type sum_expression throughout.

> Rather, we
> would need a sum vari which can have a variable size and have things
> added to it. I think this would break the rules of a vari as it's
> value cannot change once its instantiated.

But that's not true of var objects. I'll have to think
about this some more.

> I think the key here is to just write smart C++ code taking this into
> account. In terms of all the "lp__ +=" operations in a model, we
> could just as easily push the RHS of those statements into a vector
> and then call a single sum at the end. (I'm not sure if I've
> impemented the sum vari operation to handle a vector<var> type, but
> that shouldn't be too hard.

We could certainly generate smarter C++ code from stanc, if that's
what you mean. I was trying to figure out how to do it
on the C++ API side.

But we can't just accumulate because users can access intermediate
results, even for lp__. Maybe we could restrict that and instead
of treating lp__ as a variable users have access to, force them
to use something like an increment_lp(real) function.

- Bob

[Matt Hoffman]

Personally, I'm comfortable just telling people to use vectorized
functions. If it's good enough for R/Matlab/Python/Octave/everything
else it's good enough for me. But that's just because it seems easier
to me, which maybe it's not.

On a related note, I wonder if we could make something like an apply()
function that compiles to something efficient? That could be an
alternative to vectorizing every function we've got, but then again it
might be more work than it's worth.

Matt

> --
>
>

[Marcus Brubaker]

On Thu, Sep 13, 2012 at 4:12 PM, Bob Carpenter <ca...@alias-i.com> wrote:
> But we can't just accumulate because users can access intermediate
> results, even for lp__. Maybe we could restrict that and instead
> of treating lp__ as a variable users have access to, force them
> to use something like an increment_lp(real) function.

Hmmm. I'm trying to think of a reason that a user would want or need
to do something other than increment lp__ and not really coming up
with anything likely. The only thing I can think of is multiplying it
by some value to anneal some or all of the posterior. Setting it
directly (without an increment) seems like a very bad idea as it
destroys the correction due to the parameter transforms.

Cheers,
Marcus

[Bob Carpenter]

On 9/13/12 4:17 PM, Matt Hoffman wrote:
> Personally, I'm comfortable just telling people to use vectorized
> functions. If it's good enough for R/Matlab/Python/Octave/everything
> else it's good enough for me. But that's just because it seems easier
> to me, which maybe it's not.

The only problem is where there's not already an array
over which you want to vectorize.

In the current example, there's a running sum and allocating something
to put it into in order to vectorize is a pain and
makes the models hard to read (just the opposite of
the usual vectorization use case, which can make the
models easier to read).

Now if the loop-index-dependent function foo(n,...) can itself
be vectorized, then that's ideal.

> On a related note, I wonder if we could make something like an apply()
> function that compiles to something efficient? That could be an
> alternative to vectorizing every function we've got, but then again it
> might be more work than it's worth.

We could do the apply, but the question would be how to
implement it so that it vectorized rather than just looping.
The issue isn't just moving from Stan to C++, but moving
the += over var to the vectorized vari impl.

- Bob