Boardcast on division operator

[Dawen Liang]

Hi all,

I have a question about division operator. In my model, one of the parameters is a K by T non-negative matrix S, represented as: 

vector<lower=0>[K] S[T];

Each column of S is gamma distributed:

S[t] ~ gamma(alpha, alpha / S[t-1]), where t in 1:T

which is basically time-dependent by setting the expectation as the previous column. Here alpha is of type real<lower=0> and is provided. Since S[t] is of type vector, 

S[t] ~ gamma(alpha, alpha / S[t-1]);
or
S[t]~ gamma(alpha, alpha ./ S[t-1]);

is not a valid expression. From the document, it seems that broadcast is supported for "+" and "-" operators. Therefore it would be great of "/" is also broadcastable. Now my workaround is to have a transformed parameter alpha_vec:

transformed parameters {
  vector[K] alpha_vec;
  for (k in 1:K) {
    alpha_vec[k] <- alpha;
  }
}

which essentially does the broadcast myself by creating a vector with each element as alpha. But as my model is fairly complicated, I hope to optimize my code as best as I can (My model is still running, hopefully will be finished by tomorrow). Any comments about dealing with this will be greatly appreciated! 

Best,

Dawen

[Bob Carpenter]

First, I'd highly recommend just running a model for a few
iterations, then for a few more, and making sure it's
at least not crashing or getting stuck before going on
to days-long runs.

Second, you're right that we don't have an elementwise
division operator ./ defined for a real x and vector y.
It should behave like this:

(x ./ y)[n] = x ./ y[n]

Is that what MATLAB does?

It makes sense to me -- what do others think?

It's certainly easy enough to implement.

We don't need the broadcast multiplication x .* y
because this is just the same thing as regular matrix
multiplication by a scalar, x * y.

We also don't have matrix divison by a scalar. Should
that be for matrix (or vector) m and real x,
m / x or m ./ x?

- Bob

On 3/16/13 8:44 PM, Dawen Liang wrote:
> Hi all,
>
> I have a question about division operator. In my model, one of the parameters is a K by T non-negative matrix S,
> represented as:
>
> |
> vector<lower=0>[K] S[T];
> |
>
> Each column of S is gamma distributed:
>
> S[t] ~ gamma(alpha, alpha / S[t-1]), where t in 1:T
>
> which is basically time-dependent by setting the expectation as the previous column. Here alpha is of type real<lower=0>
> and is provided. Since S[t] is of type vector,
>
> |
> S[t] ~ gamma(alpha, alpha / S[t-1]);
> or

> S[t]~ gamma(alpha, alpha./ S[t-1]);

> |
>
> is not a valid expression. From the document, it seems that broadcast is supported for "+" and "-" operators. Therefore
> it would be great of "/" is also broadcastable. Now my workaround is to have a transformed parameter alpha_vec:
>
> |
> transformed parameters {
> vector[K] alpha_vec;
> for (k in 1:K) {
> alpha_vec[k] <- alpha;
> }
> }
> |
>
> which essentially does the broadcast myself by creating a vector with each element as alpha. But as my model is fairly
> complicated, I hope to optimize my code as best as I can (My model is still running, hopefully will be finished by
> tomorrow). Any comments about dealing with this will be greatly appreciated!
>
> Best,
> Dawen
>

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

[Dawen Liang]

Thanks for the advice. The sampling finished last night, but as I expected the effective sample size/r_hat is far from satisfactory, and it created a stan object occupying 8GB of my RAM (I ran 4 chains with 1000 iterations each). Don't know what to do yet...

As for the MATLAB counterpart, I think y ./ x will do the thing, where y is a number and x is a vector. MATLAB also has an explicit broadcast function bsxfun. But in this case, I think y / x, which is interpreted as broadcast, would be better, because "." means "element-wise" to many people, which implies both arguments are vectors. Anyway, it's just my thought. 

I also noticed that in FAQ about speeding up the code, there mentions that instead of doing (where y is a N*K matrix):

for (k in 1:K) { 

    for (i in 1:N) { 

        y[i,k] ~ normal(0, 1); 

    } 

}

A much faster way would be: 

y ~ normal(0,1);

However, in my model, we have a matrix<lower=0>[F, K] D and I wrote:

log(D) ~ normal(0, 1); 

or

D ~ lognormal(0, 1);

It gave me a syntax error on both. And just now, I found that it seems that transpose operator (') doesn't work in RStudio, as it kept asking me to enter another ' to balance. If there is any misunderstanding of the usage, please let me know. 

Best,

Dawen

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/groups/opt_out.

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+unsubscribe@googlegroups.com.

[Bob Carpenter]

On 3/17/13 7:58 PM, Dawen Liang wrote:
> Thanks for the advice. The sampling finished last night, but as I expected the effective sample size/r_hat is far from
> satisfactory,

This might be for other reasons.

If you tell us more about your model, we might be
able to suggest something.

> and it created a stan object occupying 8GB of my RAM (I ran 4 chains with 1000 iterations each). Don't
> know what to do yet...

How many parameters did you have such that 4 chains of 1K
iterations took 8GB of RAM? Is that in RStan or using the
new command from the command-line to calculate R-hat and ESS
and whatnot?

For memory, you can stream the data out to file rather than saving
it. The command-line version streams, and I believe
RStan is set up to turn off internal saving and turn on
streaming (someone else will know how to do it if possible).

You might want to try a smaller version of the model or
less data until you're confident it's all working, then turn
it loose on more data.

> As for the MATLAB counterpart, I think y ./ x will do the thing, where y is a number and x is a vector. MATLAB also has
> an explicit broadcast function bsxfun. But in this case, I think y / x, which is interpreted as broadcast, would be
> better, because "." means "element-wise" to many people, which implies both arguments are vectors. Anyway, it's just my
> thought.

The problem in my mind is that we already use / for
matrix division. And this isn't a matrix division
operation, which would be y / x = y * inv(x).

I guess it's going to be confusing either way.

> I also noticed that in FAQ about speeding up the code, there mentions that instead of doing (where y is a N*K matrix):
>
> for (k in 1:K) {
> for (i in 1:N) {
> y[i,k] ~ normal(0, 1);
> }
> }
>
> A much faster way would be:
>
> y ~ normal(0,1);

Sorry about that.

The FAQ got ahead of what we've implemented.
We've only vectorized probability functions for vectors,
row vectors and 1D arrays. We will be doing matrices
and general arrays in the very near future (we just
laid the metaprogramming groundwork in the code last week).

For now, the best you can do is:

> for (i in 1:N) {

> y[i] ~ normal(0, 1);

It's probably not worth flattening the matrix into a
single size (F * K) array at this point, but it's possible
with some index fiddling.

> However, in my model, we have a matrix<lower=0>[F, K] D and I wrote:
>
> log(D) ~ normal(0, 1);
> or
> D ~ lognormal(0, 1);

These two things are not at all the same. We work through
this very example in section 16.2 of the manual.

> It gave me a syntax error on both.

Right. Until we add matrices to our vectorization,
it needs to be:

for (f in 1:F)
D[f] ~ lognormal(0,1);

It's a bit more efficient here (for now, until we vectorize)
to just declare it as an array,

real<lower=0> D[F,K];

or an array of vectors,

vector<lower=0>[K] D[F];

The reason is that we can then save an allocate/copy
step for D[f].

> And just now, I found that it seems that transpose operator (') doesn't work in
> RStudio, as it kept asking me to enter another ' to balance.

Is that in an editor? Or was it in feeding the
model in as a string argument somewhere? (I don't
like the latter practice -- I think it makes more sense
to keep the models in their own files, but I seem
to be in the minority with this opinion.)

> If there is any misunderstanding of the usage, please let
> me know.

No, it's an apostrophe for transpose.

I have no idea what's up with RStudio. You could always put the
model in another file and just call it by name from the stan()
command.

- Bob

[Marcus Brubaker]

On Sun, Mar 17, 2013 at 5:45 PM, Bob Carpenter <ca...@alias-i.com> wrote:

First, I'd highly recommend just running a model for a few
iterations, then for a few more, and making sure it's
at least not crashing or getting stuck before going on
to days-long runs.

Second, you're right that we don't have an elementwise
division operator ./ defined for a real x and vector y.
It should behave like this:

   (x ./ y)[n] = x ./ y[n]

Is that what MATLAB does?

For a scalar x and matrix/vector y x./y is defined as you said.

 

It makes sense to me -- what do others think?

It's certainly easy enough to implement.

We don't need the broadcast multiplication x .* y
because this is just the same thing as regular matrix
multiplication by a scalar, x * y.

MATLAB only broadcasts scalars in it's operations, but both x * y and x.*y are defined and valid in the scalar/matrix case.

 

We also don't have matrix divison by a scalar.  Should
that be for matrix (or vector) m and real x,
m / x or m ./ x?

If we can do x ./ m we should definitely be able to do m ./ x as well.

Cheers,

Marcus 

 

[Marcus Brubaker]

On Sun, Mar 17, 2013 at 9:10 PM, Bob Carpenter <ca...@alias-i.com> wrote:

As for the MATLAB counterpart, I think y ./ x will do the thing, where y is a number and x is a vector. MATLAB also has
an explicit broadcast function bsxfun. But in this case, I think y / x, which is interpreted as broadcast, would be
better, because "." means "element-wise" to many people, which implies both arguments are vectors. Anyway, it's just my
thought.

The problem in my mind is that we already use / for
matrix division.  And this isn't a matrix division
operation, which would be y / x = y * inv(x).

I guess it's going to be confusing either way.

I'm not sure it's confusing.  In the expressions A\B and B/A is matrix division.  If the operation is to be performed elementwise then it's A./B and B./A whether A and B are scalars, matrices or there's broadcasting involved.

Cheers,

Marcus 

[Dawen Liang]

Hi Bob,

Thanks very much for the timely reply -- 

On Sun, Mar 17, 2013 at 9:10 PM, Bob Carpenter <ca...@alias-i.com> wrote:

This might be for other reasons.

If you tell us more about your model, we might be
able to suggest something.

 

I am attaching my model (as a .txt file). Basically it's kind of a non-negative matrix factorization model. I want to factorize the input X (a F by T non-negative matrix) as:

X = D*(S.*Z) + E

where D is a F*K matrix and S/Z are both K*T matrix. Here "*" denotes matrix multiplication and ".*" denotes element-wise multiplication. D is lognormal distributed and S is gamma distributed with each column has the expectation value as the previous column. Z is a binary matrix (Bernoulli distributed) which is sparsely constructed. E is i.i.d. Gaussian noise. Note that the order of dimensions in the model is not necessarily following what I am describing here, for example, all I want to do with matrix X is to obtain the t'th column so in the model I pass X as X', which makes it possible to just write down X[t]. Hope my description and the attached model together makes it clear. Please let me know if there is any question. 

I just did a toy run with very little data and noticed that the sampled Z is not really binary. If I changed the type of Z from real<lower=0, upper=1> to int<lower=0, upper=1>, then I am getting an error about integer parameters are not allowed -- I guess that's a problem. 

How many parameters did you have such that 4 chains of 1K

iterations took 8GB of RAM?  Is that in RStan or using the
new command from the command-line to calculate R-hat and ESS
and whatnot? 

For memory, you can stream the data out to file rather than saving
it.  The command-line version streams, and I believe
RStan is set up to turn off internal saving and turn on
streaming (someone else will know how to do it if possible).

You might want to try a smaller version of the model or
less data until you're confident it's all working, then turn
it loose on more data.

As you can see, there are three matrices as parameter and on real data, K, T, F can be hundreds even thousands. But this streaming thing seems very promising -- could you please tell me more about it or point me to the reference? 

 

However, in my model, we have a matrix<lower=0>[F, K] D and I wrote:

log(D) ~ normal(0, 1);
or
D ~ lognormal(0, 1);

These two things are not at all the same.  We work through
this very example in section 16.2 of the manual.

Oops, I saw the example without reading the contents carefully. My bad. 

 

>  And just now, I found that it seems that transpose operator (') doesn't work in

RStudio, as it kept asking me to enter another ' to balance.

No, it's an apostrophe for transpose.

I have no idea what's up with RStudio.  You could always put the
model in another file and just call it by name from the stan()
command.

Yes, that's what I am doing now. It won't compile with apostrophe if written in a .R file. 

Best,
Dawen

[Ben Goodrich]

On Sunday, March 17, 2013 11:54:32 PM UTC-4, Dawen Liang wrote:

>  And just now, I found that it seems that transpose operator (') doesn't work in

RStudio, as it kept asking me to enter another ' to balance.

No, it's an apostrophe for transpose.

I have no idea what's up with RStudio.  You could always put the
model in another file and just call it by name from the stan()
command.

Yes, that's what I am doing now. It won't compile with apostrophe if written in a .R file. 

Does it work in RStudio if you use double-quotation marks to define the string that gets passed to model_code? Like

model <-
"
...
y' * x
...
"

stan(model_code = model)
 

Ben

[Dawen Liang]

Yes, it works. I am too used to single-quote to remember double-quote... 

--
You received this message because you are subscribed to the Google Groups "stan users mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.