ginac for symbolic math in C++

[Bob Carpenter]

On 3/20/12 12:10 PM, Ben Goodrich wrote:
> Quoting Bob Carpenter <ca...@alias-i.com>:
>> Makes me wonder what else Ben has up his sleeves ...
>
> Another symbolic possibility (C++) is ginac
>
> http://www.ginac.de/
>
> I didn't suggest it before because it
>
> -- Is GPL3+, so we couldn't distribute it with stan
> -- Has its own symbolic scalar type ("ex")

Did you (Ben) swallow a math software encyclopedia? (Not
that there's anything wrong with that.)

My understanding of GPL and BSD "compatibility" is
derived from:

http://www.gnu.org/licenses/gpl-faq.html#WhatDoesCompatMean

It seems only the "combination" need be licensed as
GPL, whereas the original software can remain BSD.

So I was thinking that we'd distribute the rstan code
with Stan. Stan remains BSD, but rstan itself must be GPL-ed,
as it combines Stan (BSD) and R (GPL) code. I have to check
again whether we can put rstan on Google Code under those
conditions, or if we could somehow BSD-license rstan. It's
not like we're going to distribute R or other R libraries like Rcpp.

> However, iff we could get mixed scalar type operations to work, then it
> would take very little additional effort to have an internal test suite
> (that we didn't distribute) which basically used the existing templated
> stan code but with ex as the scalar type instead of var. Then you can
> differentiate with respect to the primitives, plug in numerical values,
> etc.

> But, at the moment, writing up a small sympy script takes less effort.

Definitely. Thanks again for the pointer and
examples.

- Bob

[Ben Goodrich]

On Tuesday, March 20, 2012 2:55:59 PM UTC-4, Bob Carpenter wrote:

>> Makes me wonder what else Ben has up his sleeves ...
>
> Another symbolic possibility (C++) is ginac
>
> http://www.ginac.de/
>
> I didn't suggest it before because it
>
> -- Is GPL3+, so we couldn't distribute it with stan
> -- Has its own symbolic scalar type ("ex")

Did you (Ben) swallow a math software encyclopedia?  (Not
that there's anything wrong with that.)

No, but I loosely follow sage and they include essentially all GPL2-compatible math software out there. And sage does symbolics with a home-brew python interface to ginac, although they also support sympy as an add-on package.
 

My understanding of GPL and BSD "compatibility" is
derived from:

     http://www.gnu.org/licenses/gpl-faq.html#WhatDoesCompatMean

It seems only the "combination" need be licensed as
GPL, whereas the original software can remain BSD.

IANAL, nor have a swallowed a software licensing encyclopedia, and have been told by actual lawyers to STFU many, many times, but that isn't how I would interpret the FSF. I think that link is analogous to what we call "type promotion" in that if you are working on a GPL2+ project and you want to link to a BSD library that is okay because there is nothing in BSD that is inconsistent with GPL2+. But if you are working on a BSD project (like stan) and want to link to a GPL3+ library (like ginac), then the result would be a derived work of the GPL3+ library and need to be licensed under GPL3+. Of course, we could have a separate GPL3+ project that did this or do it internally but not distribute it. My interpretation is based on the first paragraph of

http://www.gnu.org/licenses/lgpl-java.html
 

So I was thinking that we'd distribute the rstan code
with Stan.  Stan remains BSD, but rstan itself must be GPL-ed,
as it combines Stan (BSD) and R (GPL) code.  I have to check
again whether we can put rstan on Google Code under those
conditions, or if we could somehow BSD-license rstan.  It's
not like we're going to distribute R or other R libraries like Rcpp.

Its hard for me to see how rstan could be anything but GPL2+ due to linking against R and / or Rcpp headers, but I don't know if that means we should separate the repository or host rstan somewhere besides Google Code.

Ben

[Bob Carpenter]

I think we agree. I can't even talk about this
stuff coherently. No wonder legalese has so many
proper names. What Ben's suggesting license-wise is:

* R uses GPL.
This is R's decision and we can't change it.

* Stan uses BSD.
Stan only depends on Eigen and Boost,
which have LGPL and Boost licenses repsectively,
which I'm pretty sure are compatible with our
leaving Stan as BSD.

* rstan uses GPL.
It depends on both Stan and R.

I think this will work for us. The important thing is
that rstan acts as an IP firewall the GPL of R and
the BSD of Stan.

I just checked and Google code allows GPL v3 (I seem to
remember them formerly not allowing GPL). But you only
get a single combobox on the admin page requiring you
to select a single license for the "code".

So let's start a second project for rstan. It's no
big deal and cleaner in many ways. We'll probably have
to do the same for MATLAB (god only knows what kind of
proprietary license hell they impose) and Python
interfaces we promised to implement as part of the grant.

I was hoping that rstan could remain BSD and we could
distribute it that way as long as we didn't distribute R or
other GPL-ed code with the understanding that as soon as you
plugged it together with R, it became GPL. (This is
starting to sound like alchemy.) But on an umpteenth read,
I think Ben's got the right end of the stick and we'll have
to bite the bullet and separate the package distributions.

- Bob

[Marcus Brubaker]

IANAL, but I have been around open source for around 15 years and seen these discussions many many times...

On Tue, Mar 20, 2012 at 3:11 PM, Bob Carpenter <ca...@alias-i.com> wrote:

I think we agree.  I can't even talk about this
stuff coherently.  No wonder legalese has so many
proper names.  What Ben's suggesting license-wise is:

 * R uses GPL.
   This is R's decision and we can't change it.

 * Stan uses BSD.
   Stan only depends on Eigen and Boost,
   which have LGPL and Boost licenses repsectively,
   which I'm pretty sure are compatible with our
   leaving Stan as BSD.

 * rstan uses GPL.
   It depends on both Stan and R.

I think this will work for us.  The important thing is
that rstan acts as an IP firewall the GPL of R and
the BSD of Stan.

So the licensing issues of the GPL are based on the idea of a derived work.  If you took R, changed it (e.g., forked it) then you have a derived work and it would have to be GPL.  If, on the other hand, you create a GUI which calls R but does not link to R, then you don't have a derived work and the GUI can have any license it chooses.  The issue of plugins can be tricky, but if they're dynamically loaded and not compiled in, then in general they're not considered derived works.  However, if you use an (R provided) API that makes it easier to write plugins which is GPL, then it is a derived work and it needs to be GPL as well.  If that GPL plugin calls a library which is BSD, then nothing matters.  In essence, so long as the BSD code doesn't call the R API, then it's fine.  As a result, many projects make plugin APIs LGPL or provide explicit exceptions for plugins and external interfaces.  The Linux kernel is one such example of this.

 

I just checked and Google code allows GPL v3 (I seem to
remember them formerly not allowing GPL).  But you only
get a single combobox on the admin page requiring you
to select a single license for the "code".

So I poked around and found this thread:

https://groups.google.com/d/msg/google-code-hosting/9szR61CySTU/wGdmji0x6EoJ

It seems like the answer is just select GPL for the project and indicate that most of the code is actually BSD and only the R interface is GPL.

 

So let's start a second project for rstan.  It's no
big deal and cleaner in many ways.  We'll probably have
to do the same for MATLAB (god only knows what kind of
proprietary license hell they impose) and Python
interfaces we promised to implement as part of the grant.

I'm not sure I see the problem here.  I've never heard anyone complain (or even talk) about issues with using GPL/BSD/etc code in MEX functions.  In fact, there are a large number of, for instance, GPL licensed MEX interfaces out there.  E.G., see http://www.dii.unisi.it/hybrid/tools/mex/downloads.html  Certainly there is no license issues which will impact anything other than the MEX glue code.
 

I was hoping that rstan could remain BSD and we could
distribute it that way as long as we didn't distribute R or
other GPL-ed code with the understanding that as soon as you
plugged it together with R, it became GPL.  (This is
starting to sound like alchemy.)  But on an umpteenth read,
I think Ben's got the right end of the stick and we'll have
to bite the bullet and separate the package distributions.

I would check with the R project specifically about their external interface.  It may have a different license than the program as a whole.

Cheers,
Marcus

 

[Bob Carpenter]

On 3/20/12 5:05 PM, Marcus Brubaker wrote:
> IANAL, but I have been around open source for around 15 years and seen
> these discussions many many times...

Me too, but they seem to go in circles.

> So I poked around and found this thread:
>
> https://groups.google.com/d/msg/google-code-hosting/9szR61CySTU/wGdmji0x6EoJ
>
> It seems like the answer is just select GPL for the project and indicate
> that most of the code is actually BSD and only the R interface is GPL.

I think that may wind up being awfully confusing
for people looking to download just Stan. They'll
erroneously infer that it's all GPL-ed.

> So the licensing issues of the GPL are based on the idea of a derived
> work. If you took R, changed it (e.g., forked it) then you have a
> derived work and it would have to be GPL. If, on the other hand, you
> create a GUI which calls R but does not link to R, then you don't have a
> derived work and the GUI can have any license it chooses. The issue of
> plugins can be tricky, but if they're dynamically loaded and not
> compiled in, then in general they're not considered derived works.
> However, if you use an (R provided) API that makes it easier to write
> plugins which is GPL, then it is a derived work and it needs to be GPL
> as well. If that GPL plugin calls a library which is BSD, then nothing
> matters. In essence, so long as the BSD code doesn't call the R API,
> then it's fine. As a result, many projects make plugin APIs LGPL or
> provide explicit exceptions for plugins and external interfaces. The
> Linux kernel is one such example of this.

I've heard this spun both ways. But it was definitely
what I was thinking about. So maybe we could've gotten away
with it if we only called the R executables and didn't connect
directly. But we want to connect directly, and really, we might
as well bite the bullet and GPL rstan because no one can use
it other than with the GPL-ed R package anyway.

> So let's start a second project for rstan. It's no
> big deal and cleaner in many ways. We'll probably have
> to do the same for MATLAB (god only knows what kind of
> proprietary license hell they impose) and Python
> interfaces we promised to implement as part of the grant.
>
>
> I'm not sure I see the problem here. I've never heard anyone complain
> (or even talk) about issues with using GPL/BSD/etc code in MEX
> functions. In fact, there are a large number of, for instance, GPL
> licensed MEX interfaces out there. E.G., see
> http://www.dii.unisi.it/hybrid/tools/mex/downloads.html Certainly there
> is no license issues which will impact anything other than the MEX glue
> code.

Everyone -- MEX (MATLAB executable) is for plugins.

Marcus -- do you have a pointer to an actual license?
All I can find with a search on [mex matlab license]
is licenses for MATLAB or for specific plugins people
wrote using MEX.

I was hoping they'd have a license that'd let it work
like an R plugin.

The Wikipedia has the following FUD regarding MATLAB's license:

MATLAB is a proprietary product of MathWorks, so users are subject to
vendor lock-in.[2][17] Although MATLAB Builder can deploy MATLAB
functions as library files which can be used with .NET or Java
application building environment, future development will still be
tied to the MATLAB language.

Each toolbox is purchased separately. If an evaluation licence is
requested, the MathWorks sales department requires detailed
information about the project for which MATLAB is to be evaluated.
Overall the process of acquiring a licence is expensive in terms of
money and time[citation needed]. If at all granted (which happens
often), the evaluation licence is valid for two to four weeks. A
student version of MATLAB is also available.

The licensing of Mathworks software is under investigation for
anti-trust behavior by the EU.[18]

So as long as we don't run afoul of the EU, or the East District
of Texas, we should be OK?

> I was hoping that rstan could remain BSD and we could
> distribute it that way as long as we didn't distribute R or
> other GPL-ed code with the understanding that as soon as you
> plugged it together with R, it became GPL. (This is
> starting to sound like alchemy.) But on an umpteenth read,
> I think Ben's got the right end of the stick and we'll have
> to bite the bullet and separate the package distributions.
>
>
> I would check with the R project specifically about their external
> interface. It may have a different license than the program as a whole.

The whole thing's GPL-ed as far as I know, but I'm sure
Ben will know for sure.

- Bob

[Ben Goodrich]

On Tuesday, March 20, 2012 5:40:53 PM UTC-4, Bob Carpenter wrote:

> I would check with the R project specifically about their external
> interface.  It may have a different license than the program as a whole.

Here is the headers for the API:

http://cran.r-project.org/doc/manuals/R-exts.html#Organization-of-header-files

I think it is technically possible to call dynamically loaded compiled code from R that does not link against the R headers. But the compiled code would not be able to print to R's screen, throw-and-catch errors that get to the user, interrupt long-running jobs, etc. For that reason, the manual says

"The recommended procedure to use these is to include the header file R.h in your C code by #include <R.h>"

whose license is LGPL, as explained

https://svn.r-project.org/R/trunk/doc/COPYRIGHTS

So, rstan couldn't be BSD, although we could make it LGPL. But I was under the impression that we were planning to link against the Rcpp headers, which are GPL2+ only.

Ben

[Bob Carpenter]

Oops -- just sent this to Matt. Here's one for
the list.

I've also been meaning to mention that I've
come up with a design for auto-dif that will
take more space (twice the existing space, max),
but won't require any virtual functions. I'm thinking
it'd be MUCH faster as it'll just execute a very
tight loop. We can talk about it when you're back.

More comments below.

On 3/24/12 2:19 PM, Matt Hoffman wrote:
> I think this would be great to have, although I agree it's absolutely
> not a priority for Stan 1.0. The obvious advantages are

Right. NOT for version 1.

> 1. improved compiler optimization
> 2. possibly less memory overhead
> 3. lots more can be put on the stack instead of the heap
> 4. more inlining
> 5. fewer virtual function calls.
> Furthermore, a user who's truly desperate for maximum performance
> could optimize the generated code. Which is a lot easier than starting
> from scratch.
>
> Downsides: it's potentially a lot of work, and it could be a problem
> to call arbitrary functions in other libraries.

Right. We couldn't re-use exactly what we have, just
its form. And the big problem is algorithms with
branching and/or iteration. Other than that ...

> One intermediate option that might make sense to think about is having
> a more limited code gen facility that could play well with the AD
> facility that's currently in place. We could use that to automatically
> generate, e.g., log-probability vari classes that implement the
> appropriate chain() functions. E.g., we could generate a
> normal_vvv_vari class, a normal_vvd_vari class, etc., each of which
> does all of the work that is currently being done by calling the
> square(), operator/, operator-, operator*, log() functions and their
> corresponding chain() functions. These normal_xxx_vari's would produce
> the same results, but with dramatically fewer calls to the various
> chain() functions.
>
> Just a thought.

Absolutely!

That's what I've been saying from the get go. The best place
for optimization in the current design is in partial evaluation
of the derivatives in the log probability functions and in the
matrix functions.

We could also help with numerical stability by adding more
basic distributions expressed directly in terms of the unconstrained
parameters, like

bernoulli_logit(x|alpha) =def= bernoulli(x|ilogit(alpha))

Everything's pretty stable at this point (other than the builds
on various different compilers), so it's OK with me if anyone wants
to go ahead and start doing that. Just make sure to add all the
tests. Our API doc's all over the place, so I'm not so worried about
that at this point.

I know Ben wanted to do some other multivariate distribution
parameterizations.

I'm in the middle of writing the manual. Which is actually a lot
of fun.

- Bob

[Marcus Brubaker]

On Sat, Mar 24, 2012 at 4:56 PM, Bob Carpenter <ca...@alias-i.com> wrote:

Oops -- just sent this to Matt.  Here's one for
the list.

I've also been meaning to mention that I've
come up with a design for auto-dif that will
take more space (twice the existing space, max),
but won't require any virtual functions.  I'm thinking
it'd be MUCH faster as it'll just execute a very
tight loop.  We can talk about it when you're back.

I'd be curious to hear more about this if you think it could be helped by added hands/eyes on the code.  I realize it's not a target for the 1.0 release, but it's something I might have some time to work on.
 

More comments below.

On 3/24/12 2:19 PM, Matt Hoffman wrote:

I think this would be great to have, although I agree it's absolutely
not a priority for Stan 1.0. The obvious advantages are

Right.  NOT for version 1.

1. improved compiler optimization
2. possibly less memory overhead
3. lots more can be put on the stack instead of the heap
4. more inlining
5. fewer virtual function calls.
Furthermore, a user who's truly desperate for maximum performance
could optimize the generated code. Which is a lot easier than starting
from scratch.

Downsides: it's potentially a lot of work, and it could be a problem
to call arbitrary functions in other libraries.

Right.  We couldn't re-use exactly what we have, just
its form.  And the big problem is algorithms with
branching and/or iteration.  Other than that ...

One question though is how do you easily mix the autodiff code with hand coded derivatives?  It seems like only way to do it is to completely overload grad_log_prob or am I missing something?

 

One intermediate option that might make sense to think about is having
a more limited code gen facility that could play well with the AD
facility that's currently in place. We could use that to automatically
generate, e.g., log-probability vari classes that implement the
appropriate chain() functions. E.g., we could generate a
normal_vvv_vari class, a normal_vvd_vari class, etc., each of which
does all of the work that is currently being done by calling the
square(), operator/, operator-, operator*, log() functions and their
corresponding chain() functions. These normal_xxx_vari's would produce
the same results, but with dramatically fewer calls to the various
chain() functions.

Just a thought.

Absolutely!

That's what I've been saying from the get go.  The best place
for optimization in the current design is in partial evaluation
of the derivatives in the log probability functions and in the
matrix functions.

We could also help with numerical stability by adding more
basic distributions expressed directly in terms of the unconstrained
parameters, like

  bernoulli_logit(x|alpha) =def= bernoulli(x|ilogit(alpha))

So the idea would be to implement a template specialization of, e.g., normal_log, where the returned var points to one of these normal_xxx_vari's?  I might give that a shot for the multi_normal case.

Cheers,
Marcus

[Matt Hoffman]

Yes, at the moment the only way to hand-optimize the gradient
calculation is by overloading grad_log_prob(). Which is a little bit
of a pain if you don't get it exactly right on the first try, which is
why I'm interested in code generation.

Didn't we already have a specialized normal_xxx_vari at some point? Or
did I just (try and possibly fail to) do that on my copy?

In any case, what I'd really like is to have something that
automatically generated vari classes for these distributions, without
someone needing to hand-code several versions of normal, gamma,
bernoulli, bernoulli_logit, chi, chisq, etc., etc., etc. It seems like
it shouldn't be such a huge amount of work to hack something together
to do that. For the univariate distributions, we could probably even
do the code gen using sympy.

Best,
Matt

[Bob Carpenter]

The new auto-dif design is a major change, but compared
to all the other moving parts, not that difficult. The
issue is whether it'll be worth using more space (up to
twice as much in some cases) to cut down on auto-dif time.
Assuming the extra space doesn't kill us, the time should
be faster despite the increased memory usage.
I'll send out some design notes in a followup e-mail.

The rest of this reply is about what we do now.

Our initial design was based on David Gay's RAD system,
which is now part of the Sacado module of the Trilinos
package. I found all these papers pretty hard to follow,
but it might help:

http://www.cs.sandia.gov/~dmgay/ad04_paper.pdf

The Wikipedia page was also helpful.

http://en.wikipedia.org/wiki/Automatic_differentiation

More comments inline below.

On 3/26/12 4:34 PM, Marcus Brubaker wrote:
> On Sat, Mar 24, 2012 at 4:56 PM, Bob Carpenter <ca...@alias-i.com
> <mailto:ca...@alias-i.com>> wrote:
>
> One question though is how do you easily mix the autodiff code with hand
> coded derivatives? It seems like only way to do it is to completely
> overload grad_log_prob or am I missing something?

You need to either

a. write your own gradients using grad_log_prob instead
of grad_log_prob_ad, or

b. write gradients for functions and then use those
as part of auto-dif in grad_log_prob_ad

Stan's set up to let you do (a) and make it fairly easy
to do (b). I was referring to (b), which lets you accelerate
subcomputations by unfolding their auto-difs by hand.

To define a new function with hand-coded partials and
propagation, just follow the examples that are already in
Stan. I go over an example in a bit more detail below.

Pretty much. Except you need eight template specializations
if you want to catch accelerate all the ways in normal_log may
be called.

stan::agrad::var is just a pointer to the implementation,
which needs to be of type stan::agrad::vari. stan::agrad::var
is designed with normal copy syntax.

It gets tricky with multiple arguments as you need to consider all
the possible double/var instantiations. For instance,
here's the code to extend std::pow(double,double) to handle
auto-difs. Custom user functions would be handled the same
way.

namespace stan {

namespace agrad {

namespace { // makes classes anonymous

class pow_vv_vari : public op_vv_vari {
public:
pow_vv_vari(vari* avi, vari* bvi) :
op_vv_vari(std::pow(avi->val_,bvi->val_),avi,bvi) {
}
void chain() {
if (avi_->val_ == 0.0) return;
avi_->adj_ += adj_ * bvi_->val_ * val_ / avi_->val_;
bvi_->adj_ += adj_ * std::log(avi_->val_) * val_;
}
};

class pow_vd_vari : public op_vd_vari {
public:
pow_vd_vari(vari* avi, double b) :
op_vd_vari(std::pow(avi->val_,b),avi,b) {
}
void chain() {
if (avi_->val_ == 0.0) return;
avi_->adj_ += adj_ * bd_ * val_ / avi_->val_;
}
};

class pow_dv_vari : public op_dv_vari {
public:
pow_dv_vari(double a, vari* bvi) :
op_dv_vari(std::pow(a,bvi->val_),a,bvi) {
}
void chain() {
if (adj_ == 0.0) return;
bvi_->adj_ += adj_ * std::log(ad_) * val_;
}
};
}

inline var pow(const var& base, const var& exponent) {
return var(new pow_vv_vari(base.vi_,exponent.vi_));
}

inline var pow(const var& base, const double exponent) {
if (exponent == 0.5)
return sqrt(base);
if (exponent == 1.0)
return base;
if (exponent == 2.0)
return base * base;
return var(new pow_vd_vari(base.vi_,exponent));
}

inline var pow(const double base, const var& exponent) {
return var(new pow_dv_vari(base,exponent.vi_));
}

}
}

It's a lot of code, but it's all pretty simple and organized
consistently by convention. You get
three function definitions that take arguments pow(var,var),
pow(var,double) and pow(double,var). Each of these has
a corresponding implementation class extending vari.

The implementation classes (a) compute the value, (b)
store the arguments by passing to the appropriate parent
class, and (c) implement the chain() method to
pass partial derivatives. Note that in (b), the first
argument is stored by the superclass as avi_ if it is a
stan::agrad::var and as ad_ if it is a double.

It's important to note the adj_ * in all of these. The chain()
method applies the chain rule and adj_ is the partial of the
root expression with respect to the current expression, and chain()
increments the adj_ value of subexpressions (it needs to increment here
rather than set because subexpressions can be reused).

The pow(var,double) adds some extra optimization to deal with
powers for which there are built-in functions or trivial answers.
I should really replace (base * base) with square(base). It's
another unfolding I did because we had so many squares floating
around everywhere. square(double) has just one argument, so
there's less to store.

- Bob

[Bob Carpenter]

On 3/26/12 5:55 PM, Matt Hoffman wrote:
> ...

> In any case, what I'd really like is to have something that
> automatically generated vari classes for these distributions, without
> someone needing to hand-code several versions of normal, gamma,
> bernoulli, bernoulli_logit, chi, chisq, etc., etc., etc. It seems like
> it shouldn't be such a huge amount of work to hack something together
> to do that. For the univariate distributions, we could probably even
> do the code gen using sympy.

It does seem like it'd be possible to auto-generate
some of this code. But then I'm a computer scientist,
so I'm an easy sell on automation.

There are a couple (OK, a few) things to keep in mind.

1. We need 2^N signatures for an N-ary function to
efficiently handle all the ways it might get instantiated.
Some of those you might skip because they'll be efficient
as handled now, or perhaps just not called very often.

2. The biggest savings come from clever algebra in
the gradient calc that can eliminate or combine subexpression
calculations.

Having said that, the densities aren't too bad the way they
are, but there's definitely room for speedup.

3. Vectorization is totally different. If we have
a whole bunch of normal observations, we only need
to compute log(sigma) once, and that's a huge savings.
You don't even need clever algebra here -- just a little
summation.

4. Where we're really inefficient now is in things
like dot-products. That stores a gazillion intermediate
expressions of type vari on the stack.

- Bob

[Matt Hoffman]

Hmm. Maybe a solution to vectorization is to encourage people do
everything with eigen vectors (two words!), a la Matlab/R/python. Not
that I necessarily think we should rewrite Matlab, but there's
certainly something to be said for that sort of compact notation and
the easy optimization it enables.

[Bob Carpenter]

On 3/26/12 6:46 PM, Matt Hoffman wrote:
> Hmm. Maybe a solution to vectorization is to encourage people do
> everything with eigen vectors (two words!), a la Matlab/R/python. Not
> that I necessarily think we should rewrite Matlab, but there's
> certainly something to be said for that sort of compact notation and
> the easy optimization it enables.

In order for that to work, there need to be vectorized
versions of the functions to call.

We could try to do this kind of thing automatically in
the Stan Modeling Language compiler, where if we have
a function

double foo(double);

automatically extend it elementwise to:

matrix foo(matrix);

But we don't really want

vector normal_log(vector,double,double);

to be automatically derived from

double normal_log(double,double,double);

Instead, we want to sum the result to get a
double back.

- Bob

[Matt Hoffman]

Could we do both? That is, is c++ smart enough to realize that if
we're assigning to a double then it should use the version that
returns a double and if we're assigning to a vector it should use the
version that returns a vector? Not that I can really think of a great
use case for the vector function in Stan...

[Bob Carpenter]

On 3/26/12 6:54 PM, Matt Hoffman wrote:
> Could we do both? That is, is c++ smart enough to realize that if
> we're assigning to a double then it should use the version that
> returns a double and if we're assigning to a vector it should use the
> version that returns a vector?

Afraid not. Functions are resolved based only on
their arguments. So you can't overload only on
return type.

-Bob

[Matt Hoffman]

Ah well, no great loss I think.

[Marcus Brubaker]

You could switch to returning values in the parameters.  E.G.,

void normal_log(double&,double,double,double)

instead of

double normal_log(double,double,double)

Which would then allow overloading based on return type.

[Marcus Brubaker]

Doing it automatically seems unlikely to be able to actually get efficient implementations.

Another way to do this would be to encourage people to use multivariate distributions in cases like this.  The only other piece would be to add specializations of, e.g., multi_normal_log to handle passing diagonal or scalar times identity for the covariance to make it as fast as possible.  The same strategy could be used to vectorize multi_normal by defining a different distribution, e.g., vec_multi_normal(matrix,vector,matrix).

Alternatively, we just implement vectorized versions of the distributions.  A trivial (non-optimized) version would be really easy to do which could be swapped out as people implemented better versions.

I would vote for the latter myself.

Cheers,
Marcus

[Bob Carpenter]

As you can see, the who-said-what indentations are
getting lost in the responses. Not that it matters much.

> Hmm. Maybe a solution to vectorization is to encourage people do
> everything with eigen vectors (two words!), a la Matlab/R/python.
> Not that I necessarily think we should rewrite Matlab, but there's
> certainly something to be said for that sort of compact notation and
> the easy optimization it enables.
>
>
> In order for that to work, there need to be vectorized versions of
> the functions to call.
>
> We could try to do this kind of thing automatically in the Stan

> Modeling Language compiler, ...

> Doing it automatically seems unlikely to be able to actually get
> efficient implementations.

Certainly not directly. I think Matt's point was that it might
be a good place to start thinking about optimizations.

> Another way to do this would be to encourage people to use
> multivariate distributions in cases like this.

The basic approach wouldn't be any more efficient.
Less so, in fact, because of wrappers.

> The only other piece
> would be to add specializations of, e.g., multi_normal_log to handle
> passing diagonal or scalar times identity for the covariance to make
> it as fast as possible. The same strategy could be used to vectorize
> multi_normal by defining a different distribution, e.g.,
> vec_multi_normal(matrix,vector,matrix).

This could work. Ben's been itching to be smarter about
vectors. It seems so close given how Eigen handles them.
Easy enough to do in specific cases, but I don't see how to
make the resulting expression templates play nicely with everything
else without a whole lot of work. This is the same problem
we're having trying to extend Eigen to different arg types.

The multivariate with scalar reduces to the next suggestion...

> Alternatively, we just implement vectorized versions of the
> distributions. A trivial (non-optimized) version would be really
> easy to do which could be swapped out as people implemented better
> versions.
>
> I would vote for the latter myself.

That's the route down which we started. So far,
there's only the implementation for

double normal_log(vector ys, double mu, double sigma);

Going back to the last point, we might want a vector
for mu, too, because it's often a regression,

mu[n] = beta * x[n].

- Bob

[Matt Hoffman]

This might just be my Matlab background talking, but it seems to me
like it'd be best to have everything be vectorized every which way.
And the suggestion about using multivariate distributions makes sense
in the context of the normal, but less so for distributions that don't
have such a natural multivariate analog (e.g. gamma, Poisson, etc.).
It's a pain (and somewhat error-prone) to try to code up 4 or 8
versions of every density (16 if we support the generalized
inverse-Gaussian), which is why I suggested that it'd be nice to have
a way of generating them automatically.

It seems to me like there shouldn't necessarily be an obstacle to
getting (relatively) efficient implementations automatically, at least
if we're comparing with what we're doing now. Probably they won't be
100% optimal, but they should at least be better than using autodiff
for intermediate adjoint calculations and throwing a bunch of varis on
the heap. Although that may be less true once Bob's finished the new
autodiff implementation.

Best,
Matt

[Andrew Gelman]

Ultimately it will be great for this to be set up so any model can be done well. But remember goal #1 of getting the 50 or so Bugs examples working fast.
A

[Bob Carpenter]

On 3/27/12 5:39 PM, Matt Hoffman wrote:
> This might just be my Matlab background talking,

It's a bit different in C++ than in R or MATLAB.
There's no cost for iteration/access in the C++
world (or in the compiled Stan world for that matter).

On the other hand, sometimes the vectorization is
a huge savings if you can consolidate operations.
For instance, the kind of thing Ben was talking about
for a Cholesky representation along the lines of
Eigen's LLT class.

> but it seems to me
> like it'd be best to have everything be vectorized every which way.

The problem is that vectorizing can really hurt
you if you do too many allocations/constructions or
indirections as a result. It was the motivating
case for the development of expression templates!

What you don't want to do is encourage people to write
this:

vector x, u, v;
double c;
x = u + c * v;

if you're going to have to allocate an intermediate vector
to represent (c * v). You're better off with the for loop:

for (i = 0; i < x.size(); ++i)
x[i] = u[i] + c * v[i];

With expression templates, as used in Eigen, you can get
code almost this efficient while writing it in vectorized
forms (and in other cases, much more efficient code from
writing in vectorized form). You get an expression template
for scalar multiplication class on the stack pointing
to (c,v) and another one pointing to (u[i],(c,v)). Then
operator[](int) works out to be defined like x[i] above
once all the dust settles.

Now here's the real problem. In order to get the advantages
of both expression templates *and* partial evaluation of
derivatives, we'll have to

(a) rewrite Eigen to support mixed-type binary ops; they're
up for it, but I'm not sure we are

(b) push the auto-dif partial evaluations through to the Eigen's
expression templates

The way stan::agrad::matrix is implemented now, we can only
get the advantages of partial evaluation of derivatives.
Just adding mixed binary types to Eigen only gives us the
advantage of expression templates. We need to do both (a) and
(b) to get both.

The bigger speedup will come from partially evaluating
the auto-dif, because of all the memory manipulation savings.
(And it's not just the malloc time, it's the time to put things
on the stack, run their constructors, and then run their virtual chain()
methods.)

> And the suggestion about using multivariate distributions makes sense
> in the context of the normal, but less so for distributions that don't
> have such a natural multivariate analog (e.g. gamma, Poisson, etc.).
> It's a pain (and somewhat error-prone) to try to code up 4 or 8
> versions of every density (16 if we support the generalized
> inverse-Gaussian), which is why I suggested that it'd be nice to have
> a way of generating them automatically.

I agree. The reason they're not all unfolded
every which way now is the pain and errorfulness of doing
it. I'm still hoping to nail down version 1 before
any of this.

> It seems to me like there shouldn't necessarily be an obstacle to
> getting (relatively) efficient implementations automatically, at least
> if we're comparing with what we're doing now. Probably they won't be
> 100% optimal, but they should at least be better than using autodiff
> for intermediate adjoint calculations and throwing a bunch of varis on
> the heap. Although that may be less true once Bob's finished the new
> autodiff implementation.

I'm not even sure the new auto-dif idea is a good one.
See the last message. It'll depend on the balance of
multivariate and univariate ops. The "new" auto-dif is
actually closer to some of the other versions out
there than what's in Stan. Stan's is based on Matt's
original function-pointer version, which stored pointers
to the daughters the same way the classes do now. I converted
it to an OO rep and built out the library and implemented
all the operators. What's neat is that we really leaned on
the optimizing compiler to do it -- compared to CppAD, another
tool out there, Stan runs a bit faster at -O3 than CppAD at -O3,
but about half as slow at -O0. CppAD is full of case statements,
making it very hard to extend, too. Whereas Stan now has
drop-in auto-dif extensions.

I need to doc what needs to be done to drop in a user-defined
function and/or auto-dif calculation, but it's totally possible
with the existing Stan language.

- Bob

[Bob Carpenter]

On 3/27/12 5:55 PM, Andrew Gelman wrote:
> Ultimately it will be great for this to be set up so any model can be done well. But remember goal #1 of getting the 50 or so Bugs examples working fast.
> A

I think we're pretty much there on the BUGS models
thanks to Jiqiang's efforts at coding them and
Daniel's efforts at organizing and testing them.

Daniel was just talking about setting up timing tests,
and when I left today was grinding through the effective
sample size calcs from Kenny.

Also remember goal #0 of handling the kinds of complex
and large-scale models on which BUGS and JAGS failed to
load or ground to a halt.

- Bob

[Matt Hoffman]

Yup, I think some variant of the no-free-lunch theorem applies here.
We already have a package that's optimized for the easy BUGS models.
(I.e., BUGS.)

[Andrew Gelman]

1. Bugs is not optimized for anything!
2. I don't want Stan to be optimzed for these models, I just want it to beat Bugs in wall time.

[Matt Hoffman]

Maybe it would be more accurate to say that the models are optimized for BUGS.

My point about no free lunches was just that there's no such thing as
a universally optimal algorithm (provably, for optimization).

Also, have we ever compared against BUGS, rather than JAGS? I'm
curious how they compare, speedwise.

[Ben Goodrich]

On Tuesday, March 27, 2012 7:21:21 PM UTC-4, Bob Carpenter wrote:

Now here's the real problem.  In order to get the advantages

of both expression templates *and* partial evaluation of
derivatives, we'll have to

(a) rewrite Eigen to support mixed-type binary ops;  they're
up for it, but I'm not sure we are

Reminder that we almost hacked in (a) in the mixed_type branch. In particular, it should work for Array ops like

    ArrayXd x, u, v;
    // allocate u and v

    double c = ???;

    x = u + c * v;

It just doesn't work (yet) with matrix ops because we / I haven't figured out how to specialize the gebp_traits class to work like the various combinations of complex and double. Also, I don't understand the gebp stuff to begin with.

(b) push the auto-dif partial evaluations through to the Eigen's
expression templates

True, and this would be a good bit of work, but no one is going to start on (b) until (a) happens.

Ben

[Jiqiang Guo]

>
> Also, have we ever compared against BUGS, rather than JAGS? I'm
> curious how they compare, speedwise.

Plummer has done some comparison of JAGS vs BUGS

http://martynplummer.wordpress.com/2010/09/20/how-fast-is-jags/

[Andrew Gelman]

If our graph doesn't look better than that, I'll be pretty annoyed!
A

[Bob Carpenter]

On 3/27/12 9:05 PM, Ben Goodrich wrote:
> On Tuesday, March 27, 2012 7:21:21 PM UTC-4, Bob Carpenter wrote:
>
> Now here's the real problem. In order to get the advantages
>
> of both expression templates *and* partial evaluation of
> derivatives, we'll have to
>
> (a) rewrite Eigen to support mixed-type binary ops; they're
> up for it, but I'm not sure we are
>
> Reminder that we almost hacked in (a) in the mixed_type branch. In particular, it should work for Array ops like

I'm not as optimistic as you are about how close
we got. I'm pretty sure we need to either
(a) follow their lead in std::complex and hand-build
a bunch of cases, or (b) follow their suggestion in
the feature request and generalize the way traits
are specified for operations.

I don't think their current traits can handle the mixed
operations. If they did, then I think they would've done
std::complex that way.

I don't know how far they've pushed std::complex
through things like solvers, but we'd need to do that
for our applications (at least most of them are
not mixed operations, though if we did a/B or B\a it
they could be).

I'm not 100% sure on any of this, though. It's so
intricate that I find it rather daunting.

> ArrayXd x, u, v;
> // allocate u and v
>
> double c = ???;
> x = u + c * v;

I don't get this example. Isn't ArrayXd specialized for double?

> It just doesn't work (yet) with matrix ops because we / I haven't figured out how to specialize the gebp_traits class to
> work like the various combinations of complex and double. Also, I don't understand the gebp stuff to begin with.

I'll repeat that I find it all rather daunting.

I do think the Eigen developers are willing to help out
if we're serious and we get stuck. They were really
helpful when I got stuck the first time trying to provide
operator overrides everywhere.

> (b) push the auto-dif partial evaluations through to the Eigen's
> expression templates
>
> True, and this would be a good bit of work, but no one is going to start on (b) until (a) happens.

Indeed. Getting both of these done would be
the best thing to do. I think there'd be a fairly
decent audience for it, too, completely independently of
Stan.

- Bob

[Bob Carpenter]

As long as we can choose the models, we'll do fine.
Seriously, though, these models all run pretty fast.
If those timings are seconds, the worst is under
4 minutes.

That was JAGS vs. OpenBUGS. 1.5 years ago. Presumably
on a single chain. To 1000 ESS evaluated by CODA.

To emphasize Matt's point about universality, look at the
comparison of GLM vs. non-GLM for JAGS itself. RATS
is almost twice as slow with the GLM module on, whereas
some examples are over twice as fast with it on.

- Bob

[Ben Goodrich]

On Wednesday, March 28, 2012 2:02:28 AM UTC-4, Bob Carpenter wrote:

On 3/27/12 9:05 PM, Ben Goodrich wrote:
> On Tuesday, March 27, 2012 7:21:21 PM UTC-4, Bob Carpenter wrote:
>
>     Now here's the real problem. In order to get the advantages
>
>     of both expression templates *and* partial evaluation of
>     derivatives, we'll have to
>
>     (a) rewrite Eigen to support mixed-type binary ops; they're
>     up for it, but I'm not sure we are
>
> Reminder that we almost hacked in (a) in the mixed_type branch. In particular, it should work for Array ops like

I'm not as optimistic as you are about how close
we got.  I'm pretty sure we need to either
(a) follow their lead in std::complex and hand-build
a bunch of cases, or (b) follow their suggestion in
the feature request and generalize the way traits
are specified for operations.

It is quite possible that there will be another compiler error once we get past the current compiler error, but only one compiler error (repeatedly) occurs at the moment.
 

> ArrayXd x, u, v;

> // allocate u and v
>
> double c = ???;
> x = u + c * v;

I don't get this example.  Isn't ArrayXd specialized for double?

Urgh. I meant to say that v was an Array of vars. But that should work, which is to say that we can fool Eigen into thinking that an Array of vars is close enough to an Array of doubles for it to proceed, but we have not yet fooled Eigen into aligning things in an acceptable way for matrix ops that involve gebp_traits.
 

Indeed.  Getting both of these done would be

the best thing to do.  I think there'd be a fairly
decent audience for it, too, completely independently of
Stan.

It would be better if Eigen implemented genuine mixed type operations. For example, it might be somewhat useful for us to post-multiply a matrix of ints by a vector of double-precision parameters. But I doubt Eigen would be happy with my hacks, which are really just a way to implement mixed type operations between doubles and glorified doubles.

Ben

 

[Bob Carpenter]

I'm not sure we can get the right result type. Even
their approach to complex only implemented one-sided
multiplication, if I'm recalling properly.

>
> Indeed. Getting both of these done would be
>
> the best thing to do. I think there'd be a fairly
> decent audience for it, too, completely independently of
> Stan.
>
> It would be better if Eigen implemented genuine mixed type operations.
> For example, it might be somewhat useful for us to post-multiply a
> matrix of ints by a vector of double-precision parameters. But I doubt
> Eigen would be happy with my hacks, which are really just a way to
> implement mixed type operations between doubles and glorified doubles.

Well, everything has to be a glorified double at some level.
The issue's computing the result type for (double * var) and
(var * double) args.

- Bob