something git-related is out of whack
Ben Goodrich
Can everyone post the result of
git branch -a -v
which will show the tip commit for each branch? For example, mine is
master!stan> git branch -a -v
* master a0099e9 [ahead 33] Revert "attempt to enable mixed type operations"
mixed_type 9d4c274 Merge branch 'master' of https://code.google.com/p/stan into mixed_type
remotes/origin/HEAD -> origin/master
remotes/origin/master 9193958 make a pseudo-code stan.R file
Below is the massive diffstat for the 33 extra commits I have.
Ben
master!stan> git diff --stat origin/master > /tmp/diffstat.txt
master!stan> cat /tmp/diffstat.txt
doc/manuals/gm.tex | 87 +-
doc/notes/TO-DO.txt | 78 +-
make/mac | 4 +-
make/win | 10 +-
rstan/DESCRIPTION | 13 +
rstan/NAMESPACE | 3 +
rstan/R/stan.R | 32 +
rstan/R/zzz.R | 6 +
rstan/Read-and-delete-me | 10 +
rstan/man/rcpp_hello_world.Rd | 17 +
rstan/man/rstan-package.Rd | 41 +
rstan/rstan/DESCRIPTION | 13 -
rstan/rstan/NAMESPACE | 3 -
rstan/rstan/R/stan.R | 32 -
rstan/rstan/R/zzz.R | 6 -
rstan/rstan/Read-and-delete-me | 10 -
rstan/rstan/man/rcpp_hello_world.Rd | 17 -
rstan/rstan/man/rstan-package.Rd | 41 -
rstan/rstan/src-i386/Makevars | 27 -
rstan/rstan/src-i386/Makevars.win | 4 -
rstan/rstan/src-i386/hello.cpp | 9 -
rstan/rstan/src-i386/rcpp_hello_world.cpp | 11 -
rstan/rstan/src-i386/rcpp_hello_world.h | 19 -
rstan/rstan/src-i386/rcpp_module.cpp | 82 --
rstan/rstan/src-x86_64/Makevars | 27 -
rstan/rstan/src-x86_64/Makevars.win | 4 -
rstan/rstan/src-x86_64/hello.cpp | 9 -
rstan/rstan/src-x86_64/rcpp_hello_world.cpp | 11 -
rstan/rstan/src-x86_64/rcpp_hello_world.h | 19 -
rstan/rstan/src-x86_64/rcpp_module.cpp | 82 --
rstan/rstan/src/Makevars | 27 -
rstan/rstan/src/Makevars.win | 4 -
rstan/rstan/src/hello.cpp | 9 -
rstan/rstan/src/rcpp_hello_world.cpp | 11 -
rstan/rstan/src/rcpp_hello_world.h | 19 -
rstan/rstan/src/rcpp_module.cpp | 93 --
rstan/src-i386/Makevars | 27 +
rstan/src-i386/Makevars.win | 4 +
rstan/src-i386/hello.cpp | 9 +
rstan/src-i386/rcpp_hello_world.cpp | 11 +
rstan/src-i386/rcpp_hello_world.h | 19 +
rstan/src-i386/rcpp_module.cpp | 82 ++
rstan/src-x86_64/Makevars | 27 +
rstan/src-x86_64/Makevars.win | 4 +
rstan/src-x86_64/hello.cpp | 9 +
rstan/src-x86_64/rcpp_hello_world.cpp | 11 +
rstan/src-x86_64/rcpp_hello_world.h | 19 +
rstan/src-x86_64/rcpp_module.cpp | 82 ++
rstan/src/Makevars | 27 +
rstan/src/Makevars.win | 4 +
rstan/src/hello.cpp | 9 +
rstan/src/rcpp_hello_world.cpp | 11 +
rstan/src/rcpp_hello_world.h | 19 +
rstan/src/rcpp_module.cpp | 93 ++
.../bugs_examples/vol3/hepatitis/hepatitisME.stan | 1 -
src/models/temp.stan | 28 +
src/stan/agrad/matrix.hpp | 359 ++++-
src/stan/gm/command.hpp | 5 +-
src/stan/gm/function_signatures.hpp | 31 +
src/stan/gm/generator.hpp | 5 +-
src/stan/gm/parser.hpp | 224 +++-
src/stan/io/writer.hpp | 4 +-
src/stan/math/error_handling.hpp | 436 ++++--
src/stan/math/matrix.hpp | 232 +++-
src/stan/math/matrix_error_handling.hpp | 451 +++++-
src/stan/math/special_functions.hpp | 4 +
src/stan/mcmc/nuts.hpp | 2 +-
.../multivariate/continuous/dirichlet.hpp | 46 +-
.../multivariate/continuous/inv_wishart.hpp | 61 +-
.../multivariate/continuous/lkj_corr.hpp | 91 +-
.../multivariate/continuous/lkj_cov.hpp | 117 ++-
.../multivariate/continuous/multi_normal.hpp | 212 ++--
.../multivariate/continuous/multi_student_t.hpp | 52 +-
.../multivariate/continuous/wishart.hpp | 41 +-
.../multivariate/discrete/categorical.hpp | 41 +-
.../multivariate/discrete/multinomial.hpp | 37 +-
.../distributions/univariate/continuous/beta.hpp | 37 +-
.../distributions/univariate/continuous/cauchy.hpp | 48 +-
.../univariate/continuous/chi_square.hpp | 37 +-
.../univariate/continuous/double_exponential.hpp | 49 +-
.../univariate/continuous/exponential.hpp | 51 +-
.../distributions/univariate/continuous/gamma.hpp | 84 +-
.../univariate/continuous/inv_chi_square.hpp | 40 +-
.../univariate/continuous/inv_gamma.hpp | 50 +-
.../univariate/continuous/logistic.hpp | 44 +-
.../univariate/continuous/lognormal.hpp | 66 +-
.../distributions/univariate/continuous/normal.hpp | 113 ++-
.../distributions/univariate/continuous/pareto.hpp | 47 +-
.../continuous/scaled_inv_chi_square.hpp | 46 +-
.../univariate/continuous/student_t.hpp | 56 +-
.../univariate/continuous/uniform.hpp | 39 +-
.../univariate/continuous/weibull.hpp | 80 +-
.../univariate/discrete/bernoulli.hpp | 38 +-
.../univariate/discrete/beta_binomial.hpp | 38 +-
.../distributions/univariate/discrete/binomial.hpp | 43 +-
.../univariate/discrete/hypergeometric.hpp | 40 +-
.../univariate/discrete/neg_binomial.hpp | 41 +-
.../distributions/univariate/discrete/poisson.hpp | 35 +-
src/stan/stan.dox | 14 +-
src/test/agrad/matrix_test.cpp | 1537 +++++++++++++-------
src/test/math/matrix_test.cpp | 279 ++++
src/test/prob/transform_test.cpp | 10 +-
102 files changed, 5017 insertions(+), 1782 deletions(-)
Marcus Brubaker
mbrubake@rincewind stan]$ git branch -a -v
* master f280c67 [ahead 2] Typo fix in initializing the sampler.
remotes/origin/HEAD -> origin/master
remotes/origin/master 60e8c4e Fixes so that new version of multi_normal_log passes unit tests.
Cheers,
Marcus
Bob Carpenter
-------------------------------------------
~/stan>git branch -a -v
* master 60e8c4e Fixes so that new version of
multi_normal_log passes unit tests.
remotes/origin/HEAD -> origin/master
remotes/origin/master 60e8c4e Fixes so that new version of
multi_normal_log passes unit tests.
~/stan>git status
# On branch master
nothing to commit (working directory clean)
-----------------------------
Could this be related to the work in R we did on
a branch on my machine at work? I don't know that
there's much in rstan/rstan other than boilerplate
generated by rcpp.
I wrote a bunch of extra functions for the distributions
to get rid of the C++-11 features. Can you see a diff
for any of these files?
- Bob
Ben Goodrich
git checkout master
git pull origin master
to pull what you all had done into my local master branch. I thought that was the same as
git checkout master
git pull
but apparently the former does not update my laptop with the state of the Google server.
So, now when I do
master!stan> git pull --dry-run
From https://code.google.com/p/stan
9193958..60e8c4e master -> origin/master
master!stan> git pull
From https://code.google.com/p/stan
9193958..60e8c4e master -> origin/master
Already up-to-date.
master!stan> git branch -a -v
mixed_type 9d4c274 Merge branch 'master' of https://code.google.com/p/stan into mixed_type
remotes/origin/HEAD -> origin/master
master!stan>
It looks approximately the same as what Markus and Bob show, and the last three entries of my git log are now
master!stan> git log --format=oneline -n 3
a0099e98fb93754c81f73672bdb80f0295c89bae Revert "attempt to enable mixed type operations"
3d22db52e5dab25c6b038f02b17f652c7237b573 attempt to enable mixed type operations
60e8c4e76eff76878a7abcfbf0dcbd31784e6546 Fixes so that new version of multi_normal_log passes unit tests.
The reason I am now ahead by 2 commits is that I accidentally committed a potentially invasive patch to my local master branch instead of leaving it in my mixed_type branch, so I then reverted it in my master branch but not in my mixed_type branch.
So, carry on.
Ben
Bob Carpenter
stan/lib/Eigen/src/Core/CwiseBinaryOp.h
you'll find the macro
EIGEN_CHECK_BINARY_COMPATIBILITY
which is what's throwing our error:
YOU_MIXED_DIFFERENT_NUMERIC_TYPES_...
The doc above says you can't mix different types.
But they do have a special case for complex inside
the first EIGEN_STATIC_ASSERT.
When they say "vectorization", I think they're
referring to their fixed size optimizations
and things like GPU optimizations (some of which
only work on 32-bit float, I think).
You can see the template specializations for the product
type in the following file starting at line 119:
Core/Product.h
For instance, gemm is vectorized floating point operations.
I'm pretty sure this determines the expression template
return type for a product and is responsible for figuring
out whether to actually multiply or use something like
gemm vectorization. Note that Large=2 and Small=3 are just
enum values separated from 1, and all Dynamic dimensions
are Large given the categorization template metaprogram
up at line 62 of Prouct.h.
This all plays out in
Core/products/GeneralMatrixMatrix.h
where they actually define the access operators that eventually
resolve the expression templates.
Every time I dive into this stuff I get scared and back off.
- Bob
Marcus Brubaker
If you look at line 109 of
stan/lib/Eigen/src/Core/CwiseBinaryOp.h
you'll find the macro
EIGEN_CHECK_BINARY_COMPATIBILITY
which is what's throwing our error:
YOU_MIXED_DIFFERENT_NUMERIC_TYPES_...
The doc above says you can't mix different types.
But they do have a special case for complex inside
the first EIGEN_STATIC_ASSERT.
When they say "vectorization", I think they're
referring to their fixed size optimizations
and things like GPU optimizations (some of which
only work on 32-bit float, I think).
I believe vectorization probably more aptly refers to SIMD instructions (think SSE or, if your memory is that good, MMX) which work on both 32bit and 64bit platforms with both floats and doubles. However, they are only relevant for basic numeric instructions and loosing vectorization for the double OP vari cases doesn't really matter at this point.
You can see the template specializations for the product
type in the following file starting at line 119:
Core/Product.h
For instance, gemm is vectorized floating point operations.
I'm pretty sure this determines the expression template
return type for a product and is responsible for figuring
out whether to actually multiply or use something like
gemm vectorization. Note that Large=2 and Small=3 are just
enum values separated from 1, and all Dynamic dimensions
are Large given the categorization template metaprogram
up at line 62 of Prouct.h.
This all plays out in
Core/products/GeneralMatrixMatrix.h
where they actually define the access operators that eventually
resolve the expression templates.
Every time I dive into this stuff I get scared and back off.
It seems like the limiting factor here is really that we're defining a var scalar type and constructing an Eigen::Matrix<var> when Eigen really doesn't support non-basic types (except for complex) for entries.
One way around this is might be to define a matrix variant of the vari class. E.G., something like:
template<int Rows, int Cols>
class vari_matrix : public chainable {
const Eigen::Matrix<double,Rows,Cols> val_;
...
}
This would allow for fast (and vectorized) implementations of things like multiply(matrix_d,matrix_v). However, this seems likely to be a rather large design change and would end up being much messier than if Eigen just allowed (non-vectorized) mixed-type operations.
Cheers,
Marcus
Ben Goodrich
I believe vectorization probably more aptly refers to SIMD instructions (think SSE or, if your memory is that good, MMX) which work on both 32bit and 64bit platforms with both floats and doubles. However, they are only relevant for basic numeric instructions and loosing vectorization for the double OP vari cases doesn't really matter at this point.
I had speculated that disabling vectorization via one of the options listed at
http://eigen.tuxfamily.org/dox/TopicPreprocessorDirectives.html
would avoid the compiler errors at least, but no luck. The traits.madd() stuff in GeneralBlockPanelKernel.h still gets called and still yields
error: non-const lvalue reference to type 'AccPacket' (aka 'stan::agrad::var') cannot bind to a value of unrelated type
'RhsPacket' (aka 'double')
traits.madd(A0,B0,C0,T0);
On to plan Q.
Ben
Bob Carpenter
On 3/12/12 12:33 PM, Marcus Brubaker wrote:
> On Mon, Mar 12, 2012 at 2:05 AM, Bob Carpenter <ca...@alias-i.com
> <mailto:ca...@alias-i.com>> wrote:
> I believe vectorization probably more aptly refers to SIMD instructions
> (think SSE or, if your memory is that good, MMX) which work on both
> 32bit and 64bit platforms with both floats and doubles.
That makes sense.
> However, they
> are only relevant for basic numeric instructions and loosing
> vectorization for the double OP vari cases doesn't really matter at this
> point.
There's also all kinds of parallelism directives tied in on
a platform-specific basis at this low level.
The code's pretty dense to work through.
>
> You can see the template specializations for the product
> type in the following file starting at line 119:
>
> Core/Product.h
>
> For instance, gemm is vectorized floating point operations.
> I'm pretty sure this determines the expression template
> return type for a product and is responsible for figuring
> out whether to actually multiply or use something like
> gemm vectorization. Note that Large=2 and Small=3 are just
> enum values separated from 1, and all Dynamic dimensions
> are Large given the categorization template metaprogram
> up at line 62 of Prouct.h.
>
> This all plays out in
>
> Core/products/__GeneralMatrixMatrix.h
>
> where they actually define the access operators that eventually
> resolve the expression templates.
>
> Every time I dive into this stuff I get scared and back off.
>
>
> It seems like the limiting factor here is really that we're defining a
> var scalar type and constructing an Eigen::Matrix<var> when Eigen really
> doesn't support non-basic types (except for complex) for entries.
Just to be clear. Eigen's fine with our var type, just
not with mixing var and double types.
> One way around this is might be to define a matrix variant of the vari
> class. E.G., something like:
>
> template<int Rows, int Cols>
> class vari_matrix : public chainable {
> const Eigen::Matrix<double,Rows,Cols> val_;
> ...
> }
This is the first thing I tried to design, but no luck.
There are several problems with this specific suggestion.
One, we almost certainly want Rows and Cols to be runtime,
not compile time. Two, we need a similarly shaped adjoint.
So let's assume we really have this:
class vari_matrix : public chainable {
Eigen::matrix<double,Dynamic,Dynamic> val_;
Eigen::matrix<double,Dynamic,Dynamic> adj_;
}
Now let's suppose we define the operators to work as we
want, so that we have:
var& vari_matrix::operator()(int m, int n);
and wrap a var_matrix around it:
class var_matrix {
vari_matrix* vmi_;
}
So when you look at this, you realize maybe we could just template
it all:
template <typename Val>
class vari {
Val val_;
Val adj_;
}
template <typename Val>
class var {
vari<Val>* vi_;
}
We have to construct that var using a new vari given val_(m,n)
and adj_(m,n), which is going to be a real drag. But this should
really be a reference, so that we could do
vari_matrix x;
x(m,n) = y;
I never figured out how to do that without a whole lot of work
on a new kind of var that acts like a view back into the matrix
and remains very inefficient.
Then how do we string together expressions like:
(x - y) * Sigma * (x - y)'
and make sure that they get evaluated with the proper expression
templates AND propagate the right derivatives? We'd have to
build up our own expressions, I think.
So what I did is just went with Matrix<var,Dynamic,Dynamic> for
matrices, wrote our own matrix ops, and just punted on expressions.
I wish to continue to emphasize that the savings from partially
evaluating the derivatives in the var-based ops will dwarf any
savings from getting the expression templates working within Eigen.
That only solves one of our problems, and even if we do it, makes
evaluating the derivatives actually harder, as we have to start
specializing Eigen's template expressions for Scalar = var.
The really slow part is the promotion of double to var, but
that can all pretty much be eliminated by just implementing our
own basic matrix ops. The only mixed ops are the very simple
ones.
At that point, we'd basically be leaning on Eigen for the linear
algebra and the basic data types, but not for the basic matrix
operations other than get/set.
> This would allow for fast (and vectorized) implementations of things
> like multiply(matrix_d,matrix_v). However, this seems likely to be a
> rather large design change and would end up being much messier than if
> Eigen just allowed (non-vectorized) mixed-type operations.
Right. But how do we get the derivatives? I think we still have
to write all of the propagation ourselves as extensions of vari_matrix
as we did for other ops by extending vari.
- Bob
Ben Goodrich
The proper git way is to create a branch on the Google server called mixed_type that we can all look at, contribute to, and possibly get to work.
[Bob, if you currently have a local branch called mixed_type, you can just throw it away with
git branch -d mixed_type
it may give you a message saying to pass -D instead]
Now, anyone should be able to work on it by doing
git checkout remotes/origin/mixed_type
# read long message regarding detached heads
git checkout -b mixed_type
At this point, your local mixed_type branch should be in the same state as the mixed_type branch on the Google server. All should be well as long as we do not accidentally commit to the master branch experimental things that are supposed to go into the mixed_type branch (unless and until we get it working).
Our test file is named problem-eigen-agrad.cpp and is just
#include <stdio.h>
#include <vector>
#include <stan/agrad/matrix.hpp>
#include <stan/agrad/agrad.hpp>
#include <Eigen/Dense>
int main() {
using namespace Eigen;
using namespace stan::agrad;
Matrix<var,Dynamic,Dynamic> a(1,2);
a << 1.0, 2.0;
Matrix<double,Dynamic,Dynamic> b(2,1);
b << -1.0, -10.0;
std::vector<var> x;
x.push_back(a(0,0));
x.push_back(a(0,1));
Matrix<var,Dynamic,Dynamic> c = a * b;
std::cout << "c = " << std::endl << c << std::endl;
return 0;
}
At the moment,
goodrich@CYBERPOWERPC:/tmp$ clang++ -I/opt/stan/src -I/opt/stan/lib problem-eigen-agrad.cpp -o problem-eigen-agrad.o
yields this compiler error:
In file included from problem-eigen-agrad.cpp:3:
In file included from /opt/stan/src/stan/agrad/matrix.hpp:7:
In file included from /opt/stan/lib/Eigen/Dense:1:
In file included from /opt/stan/lib/Eigen/Core:327:
/opt/stan/lib/Eigen/src/Core/products/GeneralBlockPanelKernel.h:561:33: error: cannot initialize a parameter of type 'const RhsScalar *' (aka 'const stan::agrad::var *') with an rvalue of
type 'const double *'
traits.unpackRhs(depth*nr,&blockB[j2*strideB+offsetB*nr],unpackedB);
^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
So, I think we need a var constructor somewhere that takes a pointer to a double. This "progress" was achieved by specializing the gebp_traits class and trying to promote double to var internally with
template<bool _ConjLhs, bool _ConjRhs>
class gebp_traits< double, stan::agrad::var, _ConjLhs, _ConjRhs >
{
public:
typedef stan::agrad::var RhsScalar;
typedef stan::agrad::var LhsScalar;
typedef typename scalar_product_traits<LhsScalar, RhsScalar>::ReturnType ResScalar;
// copy and paste logic from double * double case
And at this point Bob explains why it is harder than it looks ...
Ben
Marcus Brubaker
On 3/12/12 12:33 PM, Marcus Brubaker wrote:
There's also all kinds of parallelism directives tied in on
However, they
are only relevant for basic numeric instructions and loosing
vectorization for the double OP vari cases doesn't really matter at this
point.
a platform-specific basis at this low level.
The code's pretty dense to work through.
True, but I suspect that parallelism at that level won't help much until we're working with REALLY big matrices.
It seems like the limiting factor here is really that we're defining a
var scalar type and constructing an Eigen::Matrix<var> when Eigen really
doesn't support non-basic types (except for complex) for entries.
Just to be clear. Eigen's fine with our var type, just
not with mixing var and double types.
That's more or less what I was meaning, but yes.
This is the first thing I tried to design, but no luck.
One way around this is might be to define a matrix variant of the vari
class. E.G., something like:
template<int Rows, int Cols>
class vari_matrix : public chainable {
const Eigen::Matrix<double,Rows,Cols> val_;
...
}
There are several problems with this specific suggestion.
One, we almost certainly want Rows and Cols to be runtime,
not compile time.
If Eigen::Matrix was used as the base then Rows and Cols could be Dynamic right?
Ah, a good point. I hadn't considered the indexing issue. In essence you would need the matrix to contain references to vari-like structures.
Then how do we string together expressions like:
(x - y) * Sigma * (x - y)'
and make sure that they get evaluated with the proper expression
templates AND propagate the right derivatives? We'd have to
build up our own expressions, I think.
That's true. The one advantage that would have is allowing handling things like inverse(A) better by exploiting that dinv(A) = -inv(A) * dA * inv(A) as opposed to pushing them through an LU factorization and similarly for determinant(A)/log_determinant(A). On the whole, probably not worth the trouble.
I wish to continue to emphasize that the savings from partially
evaluating the derivatives in the var-based ops will dwarf any
savings from getting the expression templates working within Eigen.
That only solves one of our problems, and even if we do it, makes
evaluating the derivatives actually harder, as we have to start
specializing Eigen's template expressions for Scalar = var.
The really slow part is the promotion of double to var, but
that can all pretty much be eliminated by just implementing our
own basic matrix ops. The only mixed ops are the very simple
ones.
I think there are two independent motivating factors here, at least from my perspective. One is performance (e.g., working to avoid unneeded promotion to vars) and the other is simplicity of code (e.g., allowing us to simply write A<var> OP B<double> as opposed to reimplementing every op explicitly). That getting mixed types working would allow the templated expressions to potentially improve performance is really a bonus.
Right. But how do we get the derivatives? I think we still have
This would allow for fast (and vectorized) implementations of things
like multiply(matrix_d,matrix_v). However, this seems likely to be a
rather large design change and would end up being much messier than if
Eigen just allowed (non-vectorized) mixed-type operations.
to write all of the propagation ourselves as extensions of vari_matrix
as we did for other ops by extending vari.
Cheers,
Marcus
Marcus Brubaker
Unfortunately I'm out of time for the next few days to dig into this but I'll take a look at the branch once I get a few space cycles.
Cheers,
Marcus
Ben Goodrich
stan/src/lib/Eigen/src/Core/products/GeneralBlockPanelKernel.h
where the gebp_traits class is defined and specialized for various combinations of multiplication involving real and complex numbers. Bob thinks our situation is similar to the complex * real -> complex situation in that var * double -> var. Without specializing the gebp_traits class, we ran into compiler errors where it tried to assign a var to a double, which is impossible.
Almost all the changes I have made so far have been to the top of stan/src/stan/agrad/matrix.hpp .
Ben
Bob Carpenter
On 3/12/12 3:12 PM, Marcus Brubaker wrote:
> On Mon, Mar 12, 2012 at 1:47 PM, Bob Carpenter <ca...@alias-i.com
> <mailto:ca...@alias-i.com>> wrote:
>
>
>
> On 3/12/12 12:33 PM, Marcus Brubaker wrote:
>
> One way around this is might be to define a matrix variant of
> the vari
> class. E.G., something like:
>
> template<int Rows, int Cols>
> class vari_matrix : public chainable {
> const Eigen::Matrix<double,Rows,__Cols> val_;
> ...
> }
>
>
> This is the first thing I tried to design, but no luck.
>
> There are several problems with this specific suggestion.
> One, we almost certainly want Rows and Cols to be runtime,
> not compile time.
>
>
> If Eigen::Matrix was used as the base then Rows and Cols could be
> Dynamic right?
No, Eigen::Matrix has half a dozen or so template parameters,
the first of which is the type, and the second two of which are
rows and columns, which can either be fixed values or Dynamic.
matrix_d : Matrix<double,Dynamic,Dynamic>
vector_d : Matrix<double,Dynamic,1>
row_vector_d : Matrix<double,1,Dynamic>
When you provide 1 as a fixed column or row size, Eigen uses a
template specialization to implement all the vector operations
like dot().
We could get away with a templated Val. But it makes
matrix indexing by row/col very expensive because we'd need
a vari<double> that referenced a position in the matrix.
We could implement the vari<double> view into vari<Matrix>
as an extension of vari with a no-op chain(), but it's very
heavy in terms of alloc, virtual no-op calls during auto-dif,
and access during use.
>
> Then how do we string together expressions like:
>
> (x - y) * Sigma * (x - y)'
>
> and make sure that they get evaluated with the proper expression
> templates AND propagate the right derivatives? We'd have to
> build up our own expressions, I think.
>
>
>
> That's true. The one advantage that would have is allowing handling
> things like inverse(A) better by exploiting that dinv(A) = -inv(A) * dA
> * inv(A) as opposed to pushing them through an LU factorization and
> similarly for determinant(A)/log_determinant(A). On the whole, probably
> not worth the trouble.
We can implement our own expressions.
And it's precisiely the kind of handling of derivatives
that this is going to make very hard.
> I wish to continue to emphasize that the savings from partially
> evaluating the derivatives in the var-based ops will dwarf any
> savings from getting the expression templates working within Eigen.
> That only solves one of our problems, and even if we do it, makes
> evaluating the derivatives actually harder, as we have to start
> specializing Eigen's template expressions for Scalar = var.
>
> The really slow part is the promotion of double to var, but
> that can all pretty much be eliminated by just implementing our
> own basic matrix ops. The only mixed ops are the very simple
> ones.
>
>
> I think there are two independent motivating factors here, at least from
> my perspective. One is performance (e.g., working to avoid unneeded
> promotion to vars) and the other is simplicity of code (e.g., allowing
> us to simply write A<var> OP B<double> as opposed to reimplementing
> every op explicitly). That getting mixed types working would allow the
> templated expressions to potentially improve performance is really a bonus.
I completely agree that the syntax would be nice.
We can get rid of the promotion through better implementations
of multiply(), add(), etc.
The problem is that even if we implement every OP explicitly,
it's not enough. We have to implement all of the expression
templates created by those ops. When you multiply matrices in
Eigen, or multiply a scalar by a matrix, or add two matrices,
you get an expression template, not a matrix. And they're all
different expression templates depending on the op. It's only
when you copy or assign that you get everything written out to
a matrix. So let's say I implement:
Matrix<var,Dynamic,Dynamic>
operator*(Matrix<var,Dynamic,Dynamic>,
Matrix<double,Dynamic,Dynamic>);
Then I can write:
Matrix<var,Dynamic,Dynamic> A(3,3);
Matrix<double,Dynamic,Dynamic> B(3,3);
Matrix<var,Dynamic,Dynamic> C(3,3);
...
C = A * B; // OK
C = A.transpose() * B; // FAIL
C = (A * A) * B; // FAIL
The problem is that A.transpose() is not of type
Matrix<var,Dynamic,Dynamic>, nor is (A * A). They
each have their own expression template type that's
only related to the base type through the curiously
recursive template pattern. So you need to keep
implementing for all of these expression teplates.
This is what's going to be required to change Eigen to
support mixed operations internally.
I'm forcing a write to Matrix after every op in Stan's
matrix lib. So you could use Stan's matrix transpose() to
do
C = multiply(transpose(A),B);
or
C = multiply(multiply(A,B),C);
The cost for auto-dif is going to dwarf the savings from
expression templates anyway, which is why I did it this way.
In Stan's own language, you get the expressions back, but
they're implemented with the explicit-copy Stan implementations.
It's not impossible to do both right, but just getting the mixed
ops working is only the beginning. If they're done like complex,
we'll need way more derivative-specific stuff before they're
more efficient than the non-expression-templated way I have it now
with derivatives pushed through.
>
> This would allow for fast (and vectorized) implementations of things
> like multiply(matrix_d,matrix_v). However, this seems likely to
> be a
> rather large design change and would end up being much messier
> than if
> Eigen just allowed (non-vectorized) mixed-type operations.
>
>
> Right. But how do we get the derivatives? I think we still have
> to write all of the propagation ourselves as extensions of vari_matrix
> as we did for other ops by extending vari.
>
>
>
> Yeah, that would likely be necessary. You might be able to avoid some
> duplication by treating scalars as 1x1 matrices but it would certainly
> be a big change and not a priority right now if, as you've alluded to
> elsewhere, the goal is to get something mostly fully featured out there
> for people to start using.
But Eigen doesn't implement scalars as 1x1 matrices. I don't
know if it specializes Matrix<T,1,1> for added efficiency in the
really small case. Probably.
- Bob
Ben Goodrich
implementing for all of these expression templates.
http://eigen.tuxfamily.org/dox/TopicLazyEvaluation.html
Like
C = A.transpose().eval() * B; // FAIL?
C = (A * A).eval() * B; // FAIL?
BenBob Carpenter
Bob Carpenter
anything.
Using eval() should work, assuming
internal::eval<Derived>::type
evaluates to the right kind of plain-old matrix.
A matrix constructor or static cast would also work,
but wouldn't have the clean ref passing that eval()
seems to have in case of a no-op.
This is essentially what all of our plain-old functions
in matrix are doing by assigning to an Eigen::Matrix
return type.
- Bob
On 3/12/12 11:10 PM, Ben Goodrich wrote:
Ben Goodrich
Using eval() should work, assuminginternal::eval<Derived>::type
evaluates to the right kind of plain-old matrix.
(And to work around whatever goes wrong next).
Ben