dot_product and bernoulli_logit

[Matt Hoffman]

Hi folks,

I went ahead and implemented a vectorized version of the dot_product()
function for arrays (it's in stan/agrad/matrix.hpp). In an informal
test, it seems to be around 25x faster than the old version, so if
you're doing anything that can be expressed as a dot product you
should probably be using the vectorized version.

Also, I was unsatisfied with the domain errors that I was getting when
trying to test moderately large (i.e., O(100)) values for predictors
in logistic regressions, so I implemented bernoulli_logit_log(), which
is equivalent to bernoulli_log(inv_logit(theta)) except that it can
handle large values of theta (e.g., abs(theta) > 20) correctly.

So, when implementing a logistic regression in Stan, I recommend
replacing code that looks like this:

temp <- 0;
for (d in 1:D)
temp <- temp + beta[d] * x[n, d];
y[n] ~ bernoulli(inv_logit(temp));

with code that looks like this:

y[n] ~ bernoulli_logit(dot_product(beta, x[n]));

It's safer and faster.

Matt

[Ben Goodrich]

Thanks, Matt. If a logit model has a constant, is it better to include a one in the x vector or to add the constant outside the dot product?

y[i] ~ bernoulli_log(dot_product(beta, x_with_one[i,]);                     # option 1
y[i] ~ bernoulli_log(constant + dot_product(beta, x_without_one[i,]) # option 2

In other words, is the dot_product() mechanism so much faster that it is worthwhile to repeatedly do beta[1] * 1.0 ?

Also, Bob, what is currently the best way to read in a matrix from R, as opposed to particular columns of it? Do you just dump the matrix as a long vector and then in the .stan file do

data {
  matrix(rows,cols) X;
}

Ben

[Matt Hoffman]

I think it shouldn't make much of a difference; either using an
explicit constant or adding a column of 1s in the design matrix should
be about the same. I tried using an explicit intercept term instead of
a column of ones in a 24-dimensional (plus intercept) logistic
regression:

y[i] ~ bernoulli_logit(alpha + dot_product(beta, x_without_one[i]))

and the result was maybe around 5% slower (there's a nontrivial
uncertainty interval around that number). So I wouldn't sweat it too
much.

[Bob Carpenter]

On 4/6/12 11:16 AM, Matt Hoffman wrote:
> Hi folks,
>
> I went ahead and implemented a vectorized version of the dot_product()
> function for arrays (it's in stan/agrad/matrix.hpp). In an informal
> test, it seems to be around 25x faster than the old version, so if
> you're doing anything that can be expressed as a dot product you
> should probably be using the vectorized version.

Awesome, as I said yesterday. We've known for ages that
unfolding the auto-difs in the vector ops would be a huge
win. You won't get this much gain for very short vectors,
but otherwise, the speedup should be the same once you've
"paid for" the allocation of the dot-product vari object.

We should also do this for the dot_product implementations for
Eigen vectors. They'll be able to use the same vari class
the way Matt implemented it.

Then, those dot product vari can be used as the values in
matrix multiplication, which will be speeded up by an equivalent
amount.

> Also, I was unsatisfied with the domain errors that I was getting when
> trying to test moderately large (i.e., O(100)) values for predictors
> in logistic regressions, so I implemented bernoulli_logit_log(), which
> is equivalent to bernoulli_log(inv_logit(theta)) except that it can
> handle large values of theta (e.g., abs(theta)> 20) correctly.

Nice. We'd been talking about doing just that
while you (Matt) were gone.

http://lingpipe-blog.com/2012/02/16/howprevent-overflow-underflow-logistic-regression/

I'm still not exactly sure
on how to best handle the final arithmetic ops,
so there are a couple of questions inlined.

We need to deal with extreme values of the predictor
and (I'm pretty sure) keep sensible gradient
information. Matt's computed log bernoulli(n|theta) with:

T_prob ntheta = (2*n-1) * theta;
// Handle extreme values gracefully using Taylor approximations.
const static double cutoff = 20.0;
if (ntheta > cutoff)
return -exp(-ntheta);
else if (ntheta < -cutoff)
return ntheta;
else
return -log(1 + exp(-ntheta));

I am handling a similar case for the simplex transform
using a new function I added, log1p_exp():

log1p_exp(x) = log(1 + exp(x)) = log_sum_exp(0,x)

Here's how it's defined:

inline double log1p_exp(const double& a) {
// like log_sum_exp below with b=0.0
if (a > 0.0)
return a + log1p(exp(-a));
return log1p(exp(a));
}

Looking at Matt's code, I'm thinking this'll be better
as:

inline double log1p_exp(const double& a) {
// like log_sum_exp below with b=0.0
if (a > 20.0)
return a;
if (a > 0.0)
return a + log1p(exp(-a));
return log1p(exp(a));
}

What I'm wondering about is the wisdom of adding something
like:

if (a < -20) return 0;

That'll be OK if we continue to handle the derivative the
usual way.

> So, when implementing a logistic regression in Stan, I recommend
> replacing code that looks like this:
>
> temp<- 0;
> for (d in 1:D)
> temp<- temp + beta[d] * x[n, d];
> y[n] ~ bernoulli(inv_logit(temp));
>
> with code that looks like this:
>
> y[n] ~ bernoulli_logit(dot_product(beta, x[n]));

I was really thinking something like this would've
originally been written:

vector(D) beta;
matrix(N,D) x;
... dot_product(beta,x[n]) ...

> It's safer and faster.

Both very good things.

[Bob Carpenter]

On 4/6/12 12:50 PM, Ben Goodrich wrote:
> Thanks, Matt. If a logit model has a constant, is it better to include a one in the x vector or to add the constant
> outside the dot product?
>
> y[i] ~ bernoulli_log(dot_product(beta, x_with_one[i,]); # option 1
> y[i] ~ bernoulli_log(constant + dot_product(beta, x_without_one[i,]) # option 2
>
> In other words, is the dot_product() mechanism so much faster that it is worthwhile to repeatedly do beta[1] * 1.0 ?

Since it's only one entry, there'll only be a small difference
one way or the other.

You'd really just have "alpha" instead of "beta[1] * 1.0", right?
Then it's going to be close to a wash, I would think -- you'd have
to profile to test.

(Even if you did beta[1] * 1.0, I put in a test in the operator*
implementation for constant 1.0 multiplications and they just get
ignored.)

> Also, Bob, what is currently the best way to read in a matrix from R, as opposed to particular columns of it? Do you
> just dump the matrix as a long vector and then in the .stan file do
>
> data {
> matrix(rows,cols) X;
> }

If you can make the doc, I fully documented
the dump format that stan uses for reads.
It basically a subset of the R produces with dump.
The reason I say subset, as I mention in the doc,
is that if you dump out a 2x3 matrix in R,
it writes the dimensions as "2:3", whereas if you
dump out a 2x4 matrix, it writes it as "c(2,4)".
I didn't write Stan's dump reader generally enough
to deal with "2:3". At least not yet.

- Bob

[Bob Carpenter]

Ha ha! I should've read Matt's post before responding
to Ben's.

Keep in mind that these numbers will vary across platforms,
compilers, optimization levels (you are using O3, right?).

Also, if we enable thread-local auto-dif stacks, the extra
allocation involved in the "alpha + " form will be relatively
more expensive and the advantages of vectorization greater.

- Bob

[Matt Hoffman]

Hmm. Maybe we should have only one implementation of log1p_exp(),
since the difference between

log(1/(1+exp(-a))) = -log(1 + exp(-a))

and

log(1 + exp(a))

is only a couple of minus signs. Then again, if we expect
bernoulli_logit_log() to be called enough then the couple-of-percent
performance gain from rolling those minus signs into the adjoint
calculation is maybe worth it.

On the subject of implementation details, I believe that using

if (a < -20)
return exp(a);

will produce more accurate results than

if (a < -20)
return 0;

since log(1 + exp(-20)) is in fact (very slightly) greater than 0.
(This is probably not important very often, but sometimes there's a
big difference between 0 and 2e-9.) The approximation comes from
taking a 1st-order Taylor approximation of log() around 1.

Matt

[Matt Hoffman]

Yup, O3.

I think the virtual function call overhead introduced by an autodiffed
addition trumps the cost of an additional (probably well pipelined)
multiplication by 1. But again, it's a pretty trivial cost.

Matt

[Bob Carpenter]

On 4/6/12 2:10 PM, Matt Hoffman wrote:
> Hmm. Maybe we should have only one implementation of log1p_exp(),
> since the difference between
>
> log(1/(1+exp(-a))) = -log(1 + exp(-a))
>
> and
>
> log(1 + exp(a))
>
> is only a couple of minus signs.

That's what I was thinking. But the other difference
is that you're pre-unfolding the Taylor expansion here.
In general, we want to keep going with these unfoldings
to the extent that the extra code difficulty (affecting
comprehensibility and hence maintainability) is worth it.

For logistic regression, I'm thinking every little bit will help.
Of course, those - will get dwarfed if the dimensionality of
the predictors is high, as the other e-mail thread's been
pointing out.

The outer - is free because it's in a log_prob -= x
context, and it has to be a plus or a minus.

We could also add a log(1+exp(-a)) function :-)
It's going to require some discipline to choose the right
boundary here for unfolding and building in new expressions.

Ideally, we could do some of this with expression templates
that would convert log(1 + exp(a)) to log1p_exp(a) behind
the scenes.

> Then again, if we expect
> bernoulli_logit_log() to be called enough then the couple-of-percent
> performance gain from rolling those minus signs into the adjoint
> calculation is maybe worth it.
>
> On the subject of implementation details, I believe that using
>
> if (a< -20)
> return exp(a);
>
> will produce more accurate results than
>
> if (a< -20)
> return 0;
>
> since log(1 + exp(-20)) is in fact (very slightly) greater than 0.
> (This is probably not important very often, but sometimes there's a
> big difference between 0 and 2e-9.) The approximation comes from
> taking a 1st-order Taylor approximation of log() around 1.

My bad. I just put in the wrong boundary in a
cut-and-paste braino. The bound's more like -800 or so:

if (a < -800) return 0;

Below -800, double-precision exp underflows to 0
and the whole thing reduces to log(1) = 0.

We just need to make sure the gradient doesn't
also disappear!

- Bob

P.S. I knew -O3 is the way Matt rolls (quoting him from
a way-back e-mail message). I just wanted to make
sure everyone else was. Daniel was surprised how
much of a difference it made in running our entire
set of BUGS models for the model regression (not the stats
sense) test. It can easily be a factor of 10 in some
of the cases we're looking at. We leaned very heavily on
compiler optimizations to keep the code cleaner.

[Marcus Brubaker]

Hey Matt,

Is there any reason to not go ahead and implement this for Eigen::Matrix type vectors as well as std::vectors?  If you haven't already done this, let me know and I'll get started.

Cheers,
Marcus

[Matt Hoffman]

Nope, I'm just not that familiar with the Eigen types, so I didn't do
it. But it should be straightforward, I think.

Matt

[Daniel Lee]

We might be able to write a templated version to handle both std::vector and Eigen::Matrix. Although, looking at the different versions of dot_product() we have, it may not be as straightforward as I envisioned.

Marcus, I won't hold you back. Go ahead and implement it. I'll give templating some more thought.




Daniel

[Ben Goodrich]

Not that we shouldn't have a templated version or else versions for each thing separately, but stan often uses a std::vector, an Eigen::vector, or an Eigen::Array interchangeably and arbitrarily. And we haven't gotten around to standardizing on matrix_v, vector_v, etc. So, whenever that happens we should probably standardize on some conventions for when to use std::vector, an Eigen::vector, or an Eigen::Array. Personally, I would say that a function be defined on an Eigen::vector or an Eigen::matrix iff the function is doing matrix algebra such as dot products.

Also, it would be good if there were an optimized implementation for the dot product of a vector with itself.

Ben

[Bob Carpenter]

Definitely Eigen Matrix vectors are a good idea.
There are a bunch of cases.

dot_product(vector,vector)
dot_product(row_vector,vector)
dot_product(vector,row_vector)
dot_product(row_vector,row_vector)

These, along with std::vector, might all be templatable.
They can definitely use the same vari as Matt used.

Then it'd be good to push this into multiply() for
all the matrix types. Again, there's a whole lot of them.
Matrix-matrix multiplication would have one dot-product-like
value in each cell. In fact, the same vari could
be used again.

- Bob

On 4/8/12 9:53 PM, Marcus Brubaker wrote:
> Hey Matt,
>
> Is there any reason to not go ahead and implement this for Eigen::Matrix type vectors as well as std::vectors? If you
> haven't already done this, let me know and I'll get started.
>
> Cheers,
> Marcus
>

[Bob Carpenter]

Ben brings up a very good point.

On 4/9/12 12:21 AM, Ben Goodrich wrote:
> ... we should probably standardize on some conventions

> for when to use std::vector, an Eigen::vector, or an Eigen::Array. Personally, I would say that
> a function be defined on
> an Eigen::vector or an Eigen::matrix iff the function is doing matrix algebra such as dot products.

I agree. That was my original intention.

Does anyone object? It'll mean moving
Matt's unfolded dot product operation to Eigen
vectors.

The operations like sum() and mean() aren't really
matrix ops per se, but they do take whole ranges of
values. There's no way to do a sum of a 4-D array
without allowing things like sum() to apply to
array structures. (We don't have that function now,
but general summation could be added through template
metaprogramming the same way we do general gets.)

It's very easy to add features, but it'll be hard to
take them away once Stan gets into the wild. That's
one reason I've been trying to put the brakes on new
functionality (the other being that we don't have all
the old functionality cleaned up yet, though I believe
everything's at least working at this point).

But then I see the advantage of having
extra functions for arrays.
The other way to go would be to let them be
assigned or even promoted to each other, which could
also be done, though not my preference.

> Not that we shouldn't have a templated version or else versions for each thing separately, but stan often uses a
> std::vector, an Eigen::vector, or an Eigen::Array interchangeably and arbitrarily.

It's rather like Eigen::Matrix vs. Eigen::Array
that way.

Given that arrays and vectors/matrices are only
interchangeable for some operations (like array access
and sum()) but not for others (most of them, including
assignment), I'd say Stan's data types are more "arbitrary"
than "interchangeable".

> And we haven't gotten around to
> standardizing on matrix_v, vector_v, etc.

This doesn't affect the BUGS-like Stan language, but
I do think it's worth keeing our API code clean, too.

You can't do that for templated code, unfortunately.
That was the problem with the other cleaning scheme
we tried (vestigest of which may still be left) -- type
inference wouldn't run through it and we don't want to
make users overdeclare types on function calls.

Without templates, we should use these typedefs, though.
But where do we put them?

We also need to get rid of all the namespace and
file-level includes that are lurking about.

- Bob

[Marcus Brubaker]

Way ahead of you.  I've written templated versions of dot_product and pushed it though to a templated version of multiply.  It all uses the same vari that Matt implemented (with some added functions).  It's probably not the most efficient way to do it but should work fine.  Right now I'm testing things to make sure I haven't broken anything and then I'll do some bench marks, clean it up and commit it.  Hopefully by tomorrow afternoon.

Cheers,
Marcus

[Marcus Brubaker]

Ok, so I've tested the changes and, combined with some changes to multi_normal, have reduced the computation for my multivariate time-series model by more than a factor of 2.  Everything passes the unit tests but before I push this I wanted to check about two changes that I made in the process:

1. There was a specialized multiply() function which returned a scalar when the parameters were appropriately sized vectors instead of returning a 1x1 matrix.  In order to template multiply(), this had to be removed.  As far as I could see, only multi_student_t and some unit tests relied on this functionality and I've replaced it in multi_student_t with a call to dot_product and updated the unit tests.  Does anyone object to this?  If we really want to be able to treat a 1x1 matrix as a scalar, then perhaps we should include an implicit type conversion?

2. I added solve(A,b) and trisolve(A,b) functions which compute A^-1 b.  I'm not sold on the naming (and strictly speaking I only needed trisolve in the end for the changes to multi_normal), linsolve() instead of solve() would be more consistent with Matlab but I'm not sure that's significant really.

Cheers,
Marcus

[Bob Carpenter]

Thanks for running this by us. Comments inline.

On 4/9/12 1:57 PM, Marcus Brubaker wrote:
> Ok, so I've tested the changes and, combined with some changes to
> multi_normal, have reduced the computation for my multivariate
> time-series model by more than a factor of 2. Everything passes the
> unit tests but before I push this I wanted to check about two changes
> that I made in the process:
>
> 1. There was a specialized multiply() function which returned a scalar
> when the parameters were appropriately sized vectors instead of
> returning a 1x1 matrix. In order to template multiply(), this had to be
> removed. As far as I could see, only multi_student_t and some unit
> tests relied on this functionality and I've replaced it in
> multi_student_t with a call to dot_product and updated the unit tests.
>
> Does anyone object to this?

Yes, I do! If I understand you correctly, your
proposed change would fundamentally break our type
inference process.

I believe you're talking about the multiply() function
with signature:

real multiply(row_vector, vector);

If you change that to return a matrix, then we can't do:

real x;
vector(K) y;
x <- y' * y;

because (y' * y) is now a 1 x 1 matrix.

(Yes, I could change the type promotions so that a 1 x 1
matrix can be demoted to a scalar [or a scalar can be promoted
to a 1x1 matrix], a 1 x N matrix demoted to a row vector,
and an M x 1 matrix demoted to a vector, but I'd really rather
not as it makes the entire system much much more complex when
trying to figure out which function to apply and what its
return type should be.)

> If we really want to be able to treat a 1x1
> matrix as a scalar, then perhaps we should include an implicit type
> conversion?

I've designed this the other way around. If you really
want your vector products to be a matrix, declare the
vectors as matrices. That is, you can already do it
all with matrices if you declare them right:

matrix(1,1) x;
matrix(K,1) y;
x <- y' * y;

The only drawback is that matrix(K,1) doesn't support any
vector-specific operations like dot_product because it's not
a vector.

> 2. I added solve(A,b) and trisolve(A,b) functions which compute A^-1 b.

That's great.

Why are there two?

I thought we'd have for matrix m, vector v,

m \ v == inv(m) * v

and for row vector u and matrix m,

u / m == u * inv(m)

That'd be two functions

vector divide_left(matrix m, vector v);

and

row_vector divide_right(row_vector u, matrix m);

I'm happy to add the syntax linking up the division sign in
the parser/compiler.

> I'm not sold on the naming (and strictly speaking I only needed trisolve
> in the end for the changes to multi_normal), linsolve() instead of
> solve() would be more consistent with Matlab but I'm not sure that's
> significant really.

What's linsolve()?

- Bob

[Marcus Brubaker]

On Mon, Apr 9, 2012 at 2:20 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Thanks for running this by us.  Comments inline.

On 4/9/12 1:57 PM, Marcus Brubaker wrote:

Ok, so I've tested the changes and, combined with some changes to
multi_normal, have reduced the computation for my multivariate
time-series model by more than a factor of 2.  Everything passes the
unit tests but before I push this I wanted to check about two changes
that I made in the process:

1. There was a specialized multiply() function which returned a scalar
when the parameters were appropriately sized vectors instead of
returning a 1x1 matrix.  In order to template multiply(), this had to be
removed.  As far as I could see, only multi_student_t and some unit
tests relied on this functionality and I've replaced it in
multi_student_t with a call to dot_product and updated the unit tests.

>

Does anyone object to this?

Yes, I do!  If I understand you correctly, your
proposed change would fundamentally break our type
inference process.

Ah, I figured I was missing something.
 

I believe you're talking about the multiply() function
with signature:

  real multiply(row_vector, vector);

If you change that to return a matrix, then we can't do:

  real x;
  vector(K) y;
  x <- y' * y;

because (y' * y) is now a 1 x 1 matrix.

(Yes, I could change the type promotions so that a 1 x 1
matrix can be demoted to a scalar [or a scalar can be promoted
to a 1x1 matrix], a 1 x N matrix demoted to a row vector,
and an M x 1 matrix demoted to a vector, but I'd really rather
not as it makes the entire system much much more complex when
trying to figure out which function to apply and what its
return type should be.)

OK, I'll put that function back and see if I can still make the templating work for multiply.  If not I'll default back to having several multiply() functions.
 

If we really want to be able to treat a 1x1
matrix as a scalar, then perhaps we should include an implicit type
conversion?

I've designed this the other way around.  If you really
want your vector products to be a matrix, declare the
vectors as matrices.  That is, you can already do it
all with matrices if you declare them right:

  matrix(1,1) x;
  matrix(K,1) y;
  x <- y' * y;

The only drawback is that matrix(K,1) doesn't support any
vector-specific operations like dot_product because it's not
a vector.

Actually, the easy way to template dot_product() allows appropriately sized matrices to be passed in C++.  I haven't pushed this down into the modelling language though.
 

2. I added solve(A,b) and trisolve(A,b) functions which compute A^-1 b.

That's great.

Why are there two?

solve() is for a general square matrix.  trisolve() is for the case that A is triangular.
 

I thought we'd have for matrix m, vector v,

  m \ v  ==  inv(m) * v

and for row vector u and matrix m,

  u / m  ==  u * inv(m)

That'd be two functions

  vector divide_left(matrix m, vector v);

and

  row_vector divide_right(row_vector u, matrix m);

I'm happy to add the syntax linking up the division sign in
the parser/compiler.

Right, I believe there had been some email about doing that before.  I'll rework the code into that form (and add a divide_right) and leave changes to the parser to you, as I haven't taken the time to really understand what's going on in there.
 

I'm not sold on the naming (and strictly speaking I only needed trisolve
in the end for the changes to multi_normal), linsolve() instead of
solve() would be more consistent with Matlab but I'm not sure that's
significant really.

What's linsolve()?

Matlab has two basic ways to solve a linear system.  The "backslash" operators \ (and /) which correspond to functions mldivide and mrdivide for "matrix (left|right) divide".  These operators will do all kinds of magic, including handling non-square matrices in reasonable ways and trying to detect structure in square matrices (symmetry, positive definite, diagonal, triangular, etc) to speed up the solver.  linsolve on the other hand takes a third parameter which describes the properties of the matrix and the function, therefore, uses the specific form of solver implied by those properties.

Cheers,
Marcus
 

- Bob

[Bob Carpenter]

Comments inline.

On 4/9/12 2:53 PM, Marcus Brubaker wrote:
> On Mon, Apr 9, 2012 at 2:20 PM, Bob Carpenter <ca...@alias-i.com

> <mailto:ca...@alias-i.com>> wrote:

...

> I believe you're talking about the multiply() function
> with signature:
>
> real multiply(row_vector, vector);

...

> OK, I'll put that function back and see if I can still make the
> templating work for multiply. If not I'll default back to having
> several multiply() functions.

I might be able to help with the templating if you have
a specific case to look at. It's all just a matter of making
sure there's a specific best match for any instantiation that can
be worked out from the argument types alone.

> If we really want to be able to treat a 1x1
> matrix as a scalar, then perhaps we should include an implicit type
> conversion?
>
>
> I've designed this the other way around. If you really
> want your vector products to be a matrix, declare the
> vectors as matrices. That is, you can already do it
> all with matrices if you declare them right:
>
> matrix(1,1) x;
> matrix(K,1) y;
> x <- y' * y;
>
> The only drawback is that matrix(K,1) doesn't support any
> vector-specific operations like dot_product because it's not
> a vector.
>
>
>
> Actually, the easy way to template dot_product() allows appropriately
> sized matrices to be passed in C++.

I meant the dot() operation in Eigen, which is only defined
on the Vector specializations of the Matrix template class.
And it's defined to return a scalar, not a 1 x 1 matrix.

> I haven't pushed this down into the
> modelling language though.

No problem -- the API and language are two separate
concerns (although there's a dependency from the language
to the API). I'm thinking of the modeling language as
encapsulating the agrad::var vs. double distinction for
the users, while leaving them with control over
the shape types.

The variable typing and function type inference is all
explained in the current doc (which needs a proofreading
pass for clarity and typos, I'm sure).

Thanks for all the info on MATLAB. That's really helpful.
I'm psyched about the solve functions.

>
>
> 2. I added solve(A,b) and trisolve(A,b) functions which compute
> A^-1 b.
>
>
> That's great.
>
> Why are there two?
>
>
> solve() is for a general square matrix. trisolve() is for the case that
> A is triangular.

Aha. Do you implement trisolve by taking a lower or upper
triangular view of the matrix in Eigen? If so, it should be
well defined for any matrix. On the other hand, if you'd rather
have it throw a domain error if it gets a matrix which is not
triangular, that would make sense, too. I don't know if anyone
else has an opinion.

I'm guessing we'll eventually add a triangular matrix data
type, in which case it should just take that as an argument.

>
> I thought we'd have for matrix m, vector v,
>
> m \ v == inv(m) * v
>
> and for row vector u and matrix m,
>
> u / m == u * inv(m)
>
> That'd be two functions
>
> vector divide_left(matrix m, vector v);
>
> and
>
> row_vector divide_right(row_vector u, matrix m);
>
> I'm happy to add the syntax linking up the division sign in
> the parser/compiler.
>
>
> Right, I believe there had been some email about doing that before.
> I'll rework the code into that form (and add a divide_right) and leave
> changes to the parser to you, as I haven't taken the time to really
> understand what's going on in there.

Perfect. The Boost Spirit Qi parser is quite the piece
of template lazy evaluation magic combined with some serious
metaprogramming.

The basic idea's that the parser parses down to an abstract
syntax tree, which is then used to generate. Although the
grammar looks reasonably clean, the details of semantic actions
and type propagation are monstrously complicated. At least as
complex as Eigen's matrix arithmetic evaluation policies and traits
and linked metaprograms.

Thanks again for being so patient with this crazy, sprawling
system.

- Bob