Stan has matrix operator syntax

[Bob Carpenter]

Done. We now have, in groups of descending
precedence:

---------------------------
A' transpose
---------------------------
-A negate
---------------------------
A * B product
A .* B element-wise product
A ./ B element-wise division
---------------------------
A + B sum
A - B difference
---------------------------

This includes all types that are appropriate. The
result is the minimal type that can be inferred from
the type of the arguments.

So you can write multivariate expressions like:

(mu - y)' * Sigma * (mu - y)

and

sigma * Omega * sigma'

I'm going to doc this on an op-by-op basis. Product's
the worst:

(real * real) : real
(real * vector) : vector
(real * row_vector) : row_vector
(real * matrix) : matrix

(vector * real) : vector
(vector * row_vector) : matrix

(row_vector * real) : row_vector
(row_vector * vector) : real
(row_vector * matrix) : row_vector

(matrix * real) : matrix
(matrix * vector) : vector
(matrix * matrix) : matrix

And remember that int values are promoted to real where
necessary.

We'll catch as syntax errors, things like (matrix * row_vector).
Even though the sizes are dynamic and not type checked, we
don't allow this operation even if the matrix is (N x 1).

The only other op that seems reasonable to add is
power:

(real ^ real) : real (pow)

and matrix power:

(matrix ^ int) : matrix (matrix power)

and elementwise powers:

(matrix .^ real) : matrix
(vector .^ real) : vector
(row_vector .^ real) : row_vector

which has higher precedence than product.

No one seems to use an operator for inverse, but we
probably want to discourage it in favor of more
stable alternatives.

- Bob

P.S. Math terminology quiz: Product is to multiplication
as ??? is to division?

[Marcus Brubaker]

On Tue, Mar 6, 2012 at 1:56 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Done.  We now have, in groups of descending
precedence:

Awesome!  That should make my models much more readable. :)
 

We'll catch as syntax errors, things like (matrix * row_vector).
Even though the sizes are dynamic and not type checked, we
don't allow this operation even if the matrix is (N x 1).

Any reason in particular for this?  While it seems like an unusual use case, it is potentially valid.

No one seems to use an operator for inverse, but we
probably want to discourage it in favor of more
stable alternatives.

Matlab uses \ and / to denote the solution of linear systems, i.e., x = A\b.  It helps to discourage people doing things like x = inv(A)*b.

P.S.  Math terminology quiz: Product is to multiplication
     as ??? is to division?

Quotient?

Cheers,
Marcus

[Bob Carpenter]

On 3/6/12 2:32 PM, Marcus Brubaker wrote:
> On Tue, Mar 6, 2012 at 1:56 PM, Bob Carpenter <ca...@alias-i.com
> <mailto:ca...@alias-i.com>> wrote:
>
> Done. We now have, in groups of descending
> precedence:
>
>
> Awesome! That should make my models much more readable. :)
>
> We'll catch as syntax errors, things like (matrix * row_vector).
> Even though the sizes are dynamic and not type checked, we
> don't allow this operation even if the matrix is (N x 1).
>
>
> Any reason in particular for this? While it seems like an unusual use
> case, it is potentially valid.

The main reason is to catch type errors in the programs
at compile time rather than run time.

The Eigen types implementing the basic types are:

vector_d ==> Eigen<double,Dynamic,1>
row_vector_d ==> Eigen<double,1,Dynamic>
matrix_d ==> Eigen<double,Dynamic,Dynamic>

with vector_v, row_vector_v, and matrix_v being the
same with double replaced by stan::agrad::var auto-dif
variables (Stan figures out which to use on a case-by-case
basis at compile time.)

Because Eigen specializes the class template for the
situation where one of the dimensions is 1, these are
really different types and behave differently (though
implement many of the same concepts, including matrix).

Eigen allows you to use a vector like a matrix, but
not vice versa, even if it has vector dimensions.
So if you declare a matrix variable

Eigen<double,Dynamic,Dynamic> v(5,1);

even though you have a (5 x 1) matrix, aka a vector, it
doesn't support vector operations like dot product that
are built into the specialized vector types.

I could allow the same thing as Eigen, and keep the one
way "vectors may be used as matrices, but not vice-versa",
but I thought it'd be cleaner just to fully separate
the types here. I think that makes more sense for the kind
of models we foresee people building.

Having said that, we could support lifting of vectors
to matrices:

to_matrix(vector) : matrix
to_matrix(row_vector) : matrix

which would always work. The other way around could
throw runtime errors if the dimensions are wrong:

to_vector(matrix) : vector
to_row_vector(matrix) : row_vector

When in doubt, I tend to go with "don't include it" and
"be strict in the syntax".

> No one seems to use an operator for inverse, but we
> probably want to discourage it in favor of more
> stable alternatives.
>
>
> Matlab uses \ and / to denote the solution of linear systems, i.e., x =
> A\b. It helps to discourage people doing things like x = inv(A)*b.

Cool. That'd be easy to add and potentially useful.

Is the operation left associative at the same precedence
as multiplication.

a / b * c == (a / b) * c

a * b / c == (a * b) / c

That is, do we have:

A \ b * C == (A \ b) * C

C * A \ b == (C * A) \ b

Does anyone ever use the rightward version?

b / A == b * inv(A)

and does it have the precednece I'd expect

C * b / A == (C * b) / A

b / A * C == (b / A) * C

This is bringing back my algebra days -- I used left and
right residuals in the string semigroup for models of natural
language syntax (categorial grammars based on some algebraic
work of Lambek in the 1950s). I didn't realize they were
used for matrices.

> P.S. Math terminology quiz: Product is to multiplication
> as ??? is to division?
>
>
> Quotient?

Thanks. I couldn't think of that word for the life of me.

And more seriously, thanks for the technical input.

- Bob

[Marcus Brubaker]

Fair enough but it seems unlikely to cause problems to include it for completeness.  I look at it in terms of expected behavior.  I think most people would reasonably expect that to work if they felt a need to use it.

Having said that, we could support lifting of vectors
to matrices:

  to_matrix(vector)      : matrix
  to_matrix(row_vector)  : matrix

which would always work.  The other way around could
throw runtime errors if the dimensions are wrong:

  to_vector(matrix)      : vector
  to_row_vector(matrix)  : row_vector

When in doubt, I tend to go with "don't include it" and
"be strict in the syntax".

These I think would definitely be useful regardless of the above.  For one it ensures that anything that can be done with matrices can be done with vectors if the dimensions are right.  It also, potentially, could improve efficiency.  For instance, there are certainly cases where some operations are well defined and fast on a full matrix, but other operations need to be specified on rows/columns of the matrix.

   No one seems to use an operator for inverse, but we
   probably want to discourage it in favor of more
   stable alternatives.


Matlab uses \ and / to denote the solution of linear systems, i.e., x =
A\b.  It helps to discourage people doing things like x = inv(A)*b.

Cool.  That'd be easy to add and potentially useful.

Is the operation left associative at the same precedence
as multiplication.

  a / b * c  == (a / b) * c

  a * b / c  == (a * b) / c

That is, do we have:

  A \ b * C  == (A \ b) * C

  C * A \ b  == (C * A) \ b

Yes, that's right.  See also http://www.mathworks.com/help/techdoc/matlab_prog/f0-40063.html#f0-38155
 

Does anyone ever use the rightward version?

  b / A == b * inv(A)

I've definitely used it in Matlab but it's obviously not as common.

and does it have the precednece I'd expect

  C * b / A == (C * b) / A

  b / A * C == (b / A) * C

Correct.

This is bringing back my algebra days -- I used left and
right residuals in the string semigroup for models of natural
language syntax (categorial grammars based on some algebraic
work of Lambek in the 1950s).  I didn't realize they were
used for matrices.

I haven't had to think about operator precedence and order of evaluation in ages.  Definitely has required me to dust off some long unused parts of my brain. :)

Cheers,
Marcus

[Ben Goodrich]

Having parser problems with the new syntax here. In this case, b_0 is a real vector and R is a matrix and I want the scalar rather than a 1x1 matrix.

no matches for function name="transpose"
    arg 0 type=double[1]
no matches for function name="multiply"
    arg 0 type=ill formed[0]
    arg 1 type=double[2]
no matches for function name="multiply"
    arg 0 type=ill formed[0]
    arg 1 type=double[1]

ERROR PARSING
LOCATION:  file=input; line=80, column=29

    temp <- ( b_0' * R * b_0 )[1];

DIAGNOSTIC(S) FROM PARSER:

expression is ill formed
Parser expecting: <eps>

What do I need to be doing instead?

Thanks,
Ben

[Bob Carpenter]

You need to use Eigen matrices or vectors
instead of double n. That is, instead of

double b_0[N];

you want:

vector(N) b_0;

I'm just sitting down to write that section of the
manual! Things specified with [] are array types.
Some multi-variate ops are defined on arrays, but not
many.

- Bob