My life would be a lot easier if there were a vectorized ^ operator, that is if this worked:
vector[10] x;
vector[10] a;
vector[10] d;
d = x^a;
and if this worked too:
vector[10] x;
vector[1] a;
vector[10] d;
d = x^a;
Is this part of the current vectorization pull requests? I guess strictly speaking that’s broadcasting but I wanted to check if it’s part of a current pull request.
No, we’re doing the unary ones first. But Rayleigh’s
already working on the binary ones.
Firs tof all, what do you expect
a^b
to do for various inputs and what are the type
restrictions on a and b?
I’d think we’d want a and b to be either scalars or
the exact same container type. That’s more rigid than
the probability functions, which only require the same
shape.
Then, do you expect all the other operators to work
the same way like *? So far, they’re only implemented
in a way that preserves types following standard linear
algebra conventions.
Once we know what ^ should do, we could just implement
it directly, or we can wait to see if it’ll be compatible
with the general binary function vectorization scheme.
Partly it depends on how much of a rush you’re in and whether
you’re willing to write the code for a custom one until
the framework’s finished.
Firs tof all, what do you expect
a^b
to do for various inputs and what are the type
restrictions on a and b?
Good point, I’m expecting elementwise exponentiation with broadcasting so:
real x;
real a;
real b = x^a;
real x;
vector[n] a;
vector b = x^a = [x^a[1], x^a[2], ..., x^a[n] ]
vector[n] x;
real a;
vector b = x^a = [x[1]^a, x[2]^a, ..., x[n]^a ]
vector[m] a;
vector[n] x;
vector b = x^a = [x[1]^a[1], x[2]^a[2], ..., x[n]^a[m] ]
vector[m] a;
vector[1] x;
vector b = x^a = [x[1]^a[1], x[1]^a[2], ..., x[1]^a[m] ]
vector[1] a;
vector[n] x;
vector b = x^a = [x[1]^a[1], x[2]^a[1], ..., x[n]^a[1] ]
I think everything else with vector/real types would be an error (no recycling ala R!) and I think the last two should also NOT be part of the spec unless that’s how the rest of the vectorization works.
Yes, I agree. Also think this addresses the issue of whether to prefer 1-d arrays or vectors and the manual could use a discussion of this (i.e.-look ahead a little and pick one!)
I’m expecting that +, -, .*, and ./ would all work in the same way
and ^ should follow whatever that is.
I already implemented this for the few densities I had coded using the Stan language as-is with looping but that did give me some insight into how nice having more complete vectorization would be!
I’ll wait. For the time being I’m still working on what an alternative to expose_stan_functions could look like. I have a build process with cmake that’s pretty reliable and spits out an R package with Stan functions exposed but I’m still writing interface functions half by hand to figure out what RcppExport can handle so that’s a long-winded way of saying I won’t be picking up another project until this one’s done.
Krzysztof
Responses inline.
sakrejda Developer
October 31
Firs tof all, what do you expect
a^b
to do for various inputs and what are the type
restrictions on a and b?Good point, I’m expecting elementwise exponentiation with broadcasting so:
real x;
real a;
real b = x^a;real x;
vector[n] a;
vector b = x^a = [x^a[1], x^a[2], …, x^a[n] ]vector[n] x;
real a;
vector b = x^a = [x[1]^a, x[2]^a, …, x[n]^a ]vector[m] a;
vector[n] x;
vector b = x^a = [x[1]^a[1], x[2]^a[2], …, x[n]^a[m] ]vector[m] a;
vector[1] x;
vector b = x^a = [x[1]^a[1], x[1]^a[2], …, x[1]^a[m] ]vector[1] a;
vector[n] x;
vector b = x^a = [x[1]^a[1], x[2]^a[1], …, x[n]^a[1] ]I think everything else with vector/real types would be an error (no recycling ala R!) and I think the last two should also NOT be part of the spec unless that’s how the rest of the vectorization works.
The only other vectorization now is in the densities, and
that’s slightly different because it sums. We did not
allow the equivalent of the last two cases—those are the
edge case for R’s recycling. Instead, we required that
all container arguments had to be the same size and shape.
Presumably row_vector would work the same way, but we wouldn’t
mix. Otherwise, no way to infer a return type for
vector ^ row_vector : ?
In the density case, it’s not a problem because the result just
sums and we allow this kind of mixing.
Bob_Carpenter:
I’d think we’d want a and b to be either scalars or
the exact same container type. That’s more rigid than
the probability functions, which only require the same
shape.Yes, I agree. Also think this addresses the issue of whether to prefer 1-d arrays or vectors and the manual could use a discussion of this (i.e.-look ahead a little and pick one!)
The question in my mind is whether to allow
real ^ real : real
real ^ real : real
real ^ real : real
for arrays. And do we allow matrices, as we do for unary funs?
matrix ^ matrix : matrix
matrix ^ real : matrix
real ^ matrix : matrix
And then what about arrays of matrices and ever larger container
types. OK if shapes match or if one arg is a scalar?
Bob_Carpenter:
Then, do you expect all the other operators to work
the same way like *? So far, they’re only implemented
in a way that preserves types following standard linear
algebra conventions.I’m expecting that +, -, .*, and ./ would all work in the same way
and ^ should follow whatever that is.
But the confusion here is that + and * are matrix operations. I’m
uncomfortable with:
vector * vector : vector
because we have
row_vector * vector : real
vector * row_vector: matrix
already.
Manual says we have Elementwise Arithmetic Operations .* and ./
so why not add .^ ?
Aki
Those are elementwise but they don’t broadcast reals so that’s more
straightforward than what I’m suggesting i do think we should keep that
notation for elementwise with broadcasting. K
sakrejda Developer
November 3
Those are elementwise but they don’t broadcast reals so that’s more
straightforward than what I’m suggesting i do think we should keep that
notation for elementwise with broadcasting. K
Keep what notation for what?
The .* and ./ are specifically to distinguish them from the
matrix multiplication and matrix division operations.
We don’t have a matrix exponentiation operator. If we did,
it’d be
matrix^int: matrix
to multiply a matrix times itself the specified number of times.