dot product templating, static vs. dynamic types, and NDEBUG

225 views
Skip to first unread message

Bob Carpenter

unread,
Apr 21, 2012, 3:58:34 PM4/21/12
to stan...@googlegroups.com
I'm writing the dot_self function to get a bit
more familiar with the templating and wanted to bring
some coding issues to everyone's attention.

Here's the current dot_product in stan/math/matrix.hpp:

template<int R1,int C1,int R2, int C2>
inline double dot_product(const Eigen::Matrix<double, R1, C1>& v1,
const Eigen::Matrix<double, R2, C2>& v2) {
if (v1.rows() != 1 && v1.cols() != 1)
throw std::invalid_argument("v1 must be a vector");
if (v2.rows() != 1 && v2.cols() != 1)
throw std::invalid_argument("v2 must be a vector");
if (v1.size() != v2.size())
throw std::invalid_argument("v1.size() must equal v2.size()");
return v1.dot(v2);
}

Issue 1
-------
This only compiles if one of R1 and C1 is 1
and one of R2 and C2 is 1, so the dynamic vector checks
are redundant (they'll always pass if the code compiled).

The reason is that v1.dot(v2) is only defined if both of
v1 and v2 are Eigen vector template specializations of
Matrix. A type Matrix<T,R,C> is a vector template specialization
if and only if R==1 or C==1.

This means the dynamic checks for row/col size are redundant:

if (v1.rows() != 1 && v1.cols() != 1)
throw std::invalid_argument("v1 must be a vector");
if (v2.rows() != 1 && v2.cols() != 1)
throw std::invalid_argument("v2 must be a vector");

So I removed them.

Issue 2
-------
The dynamic size match test will throw a std::invalid_argument
exception. I'm not sure we want NUTS/HMC to recover form
this. Shouldn't we just stop at this point because it indicates
a programming error? I can't think of any place this would happen
where the situation would be corrected with more sampling.


Issue 3
-------
We're compiling without the -NDEBUG flag, meaning runtime
assertions are turned on.

If we were to get rid of the size check

if (v1.size() != v2.size())
throw std::invalid_argument("v1.size() must equal v2.size()");

then the call

return v1.dot(v2);

will cause assertion failure and crash the program.

I think we want to add the -NDEBUG flag back in for speed
and control the error handling ourselves. For instance,
we don't need to check in v.dot(v) for dot_self()!

- Bob

Ben Goodrich

unread,
Apr 21, 2012, 4:09:23 PM4/21/12
to stan...@googlegroups.com
On Saturday, April 21, 2012 3:58:34 PM UTC-4, Bob Carpenter wrote:
Here's the current dot_product in stan/math/matrix.hpp:

     template<int R1,int C1,int R2, int C2>
     inline double dot_product(const Eigen::Matrix<double, R1, C1>& v1,
                               const Eigen::Matrix<double, R2, C2>& v2) {
       if (v1.rows() != 1 && v1.cols() != 1)
         throw std::invalid_argument("v1 must be a vector");
       if (v2.rows() != 1 && v2.cols() != 1)
         throw std::invalid_argument("v2 must be a vector");
       if (v1.size() != v2.size())
         throw std::invalid_argument("v1.size() must equal v2.size()");
       return v1.dot(v2);
     }

Issue 1
-------
This only compiles if one of R1 and C1 is 1
and one of R2 and C2 is 1, so the dynamic vector checks
are redundant (they'll always pass if the code compiled).

Things like this are one of the reasons that I favored defining all applicable functions on our typedefs, which hasn't gotten done yet. But I agree with everything you said.

Ben

Bob Carpenter

unread,
Apr 21, 2012, 4:19:44 PM4/21/12
to stan...@googlegroups.com
Within the Stan programming language, that's
exactly what happens.

Within the API, we've got all sorts of things
going on.

1. Sometimes there's just a fully-templated
definition in stan/math and nothing in agrad. And
there's not really a good way to template out the scalar.
The problem is that type inference won't pass through
complex typedefs.

2. Sometimes (as for dot_product), we have things
fixed for scalar type and then defined in stan/math for
double and stan/agrad for var, but using templates for
Eigen size types.

3. Sometimes we have everything expanded out for
all of our most specific types. In those cases, we
can definitely use the typedefs.

I would argue that (3) is the easiest to document and
to understand. We have more functions to doc, but
no templates winging around. I think (1) is next simplest
and (2) is most complex to doc and understand.

I think it's a wash for testing purposes.

I'm not sure whether we should just take one approach
and stick to it or decide on a per-function basis.

- Bob

Marcus Brubaker

unread,
Apr 22, 2012, 4:29:24 AM4/22/12
to stan...@googlegroups.com
On Sat, Apr 21, 2012 at 8:58 PM, Bob Carpenter <ca...@alias-i.com> wrote:
I'm writing the dot_self function to get a bit
more familiar with the templating and wanted to bring
some coding issues to everyone's attention.

Here's the current dot_product in stan/math/matrix.hpp:

   template<int R1,int C1,int R2, int C2>
   inline double dot_product(const Eigen::Matrix<double, R1, C1>& v1,
                             const Eigen::Matrix<double, R2, C2>& v2) {
     if (v1.rows() != 1 && v1.cols() != 1)
       throw std::invalid_argument("v1 must be a vector");
     if (v2.rows() != 1 && v2.cols() != 1)
       throw std::invalid_argument("v2 must be a vector");
     if (v1.size() != v2.size())
       throw std::invalid_argument("v1.size() must equal v2.size()");
     return v1.dot(v2);
   }


Are you sure this the latest version you're looking at?  The current version of master doesn't use .dot() but instead uses a dot_product specific varI.
 

Issue 1
-------
This only compiles if one of R1 and C1 is 1
and one of R2 and C2 is 1, so the dynamic vector checks
are redundant (they'll always pass if the code compiled).

The reason is that v1.dot(v2) is only defined if both of
v1 and v2 are Eigen vector template specializations of
Matrix.  A type Matrix<T,R,C> is a vector template specialization
if and only if R==1 or C==1.

This means the dynamic checks for row/col size are redundant:

   if (v1.rows() != 1 && v1.cols() != 1)
       throw std::invalid_argument("v1 must be a vector");
     if (v2.rows() != 1 && v2.cols() != 1)
       throw std::invalid_argument("v2 must be a vector");

So I removed them.

This isn't true with the dot_product vari implementation.  Those tests are needed there in that case.
 

Issue 2
-------
The dynamic size match test will throw a std::invalid_argument
exception.   I'm not sure we want NUTS/HMC to recover form
this.  Shouldn't we just stop at this point because it indicates
a programming error?  I can't think of any place this would happen
where the situation would be corrected with more sampling.

This is up to the caller to determine how to handle.  If the caller wants to, e.g., clean things up before aborting then they should be given the opportunity.
 


Issue 3
-------
We're compiling without the -NDEBUG flag, meaning runtime
assertions are turned on.

If we were to get rid of the size check

    if (v1.size() != v2.size())
       throw std::invalid_argument("v1.size() must equal v2.size()");

then the call

    return v1.dot(v2);

will cause assertion failure and crash the program.

I think we want to add the -NDEBUG flag back in for speed
and control the error handling ourselves.  For instance,
we don't need to check in v.dot(v) for dot_self()!

We probably want to leave those tests in, even with NDEBUG, at least for model compilation as it will catch errors in model specification and potentially produce useful error messages rather than crashing in unpredictable ways.

Cheers,
Marcus

Bob Carpenter

unread,
Apr 22, 2012, 5:47:35 PM4/22/12
to stan...@googlegroups.com


On 4/22/12 4:29 AM, Marcus Brubaker wrote:
> On Sat, Apr 21, 2012 at 8:58 PM, Bob Carpenter <ca...@alias-i.com <mailto:ca...@alias-i.com>> wrote:
>
> I'm writing the dot_self function to get a bit
> more familiar with the templating and wanted to bring
> some coding issues to everyone's attention.
>
> Here's the current dot_product in stan/math/matrix.hpp:
>
> template<int R1,int C1,int R2, int C2>
> inline double dot_product(const Eigen::Matrix<double, R1, C1>& v1,
> const Eigen::Matrix<double, R2, C2>& v2) {
> if (v1.rows() != 1 && v1.cols() != 1)
> throw std::invalid_argument("v1 must be a vector");
> if (v2.rows() != 1 && v2.cols() != 1)
> throw std::invalid_argument("v2 must be a vector");
> if (v1.size() != v2.size())
> throw std::invalid_argument("v1.__size() must equal v2.size()");
> return v1.dot(v2);
> }
>
>
>
> Are you sure this the latest version you're looking at? The current version of master doesn't use .dot() but instead
> uses a dot_product specific varI.

Ah. The problem was that the double version (above) didn't
match the vari version. In order to work on Matrix<double,Dynamic,Dynamic>
instances that just happen to be sized as vectors, I'll rewrite
our own dot product instead of using Eigen's.

But then I'm still not sure if we want to allow this operation.
It doesn't follow the way Eigen works, and has to resort to
runtime type checking. (But then even matrix multiplication
has to result to runtime size checking.)

What does everyone else think? That is, should we allow
the signature:

real dot_product(matrix,matrix)

I asked before, but I don't think we resolved the issue.

We can control this through the stan/gm function signatures.
So it's a separate issue as to whether the API function
can be defined generally as above.

> Issue 1
> -------
> This only compiles if one of R1 and C1 is 1
> and one of R2 and C2 is 1, so the dynamic vector checks
> are redundant (they'll always pass if the code compiled).
>
> The reason is that v1.dot(v2) is only defined if both of
> v1 and v2 are Eigen vector template specializations of
> Matrix. A type Matrix<T,R,C> is a vector template specialization
> if and only if R==1 or C==1.
>
> This means the dynamic checks for row/col size are redundant:
>
> if (v1.rows() != 1 && v1.cols() != 1)
> throw std::invalid_argument("v1 must be a vector");
> if (v2.rows() != 1 && v2.cols() != 1)
> throw std::invalid_argument("v2 must be a vector");
>
> So I removed them.
>
>
> This isn't true with the dot_product vari implementation. Those tests are needed there in that case.

Right -- see above. I figured that out as soon
as I removed them and ran the unit tests :-) That's
why there's nothing checked in yet.

> Issue 2
> -------
> The dynamic size match test will throw a std::invalid_argument
> exception. I'm not sure we want NUTS/HMC to recover form
> this. Shouldn't we just stop at this point because it indicates
> a programming error? I can't think of any place this would happen
> where the situation would be corrected with more sampling.
>
>
> This is up to the caller to determine how to handle. If the caller wants to, e.g., clean things up before aborting then
> they should be given the opportunity.

For the API, the caller's just going to get the exception.

What we should really do is add the Boost-type configurable
type checking throughout, but that's a huge pain.

But what I meant was what should the Stan programming
language do --- NUTS or HMC is the caller in that case.

For things like variates out of support, we're throwing
domain error exceptions from inside the densities and
NUTS/HMC treat those the same way as zero probability
cases and reject the parameters.

I'm thinking we want NUTS/HMC to abort when they hit
this kind of problem because they're not going to be
able to sample their way out of it.

> Issue 3
> -------
> We're compiling without the -NDEBUG flag, meaning runtime
> assertions are turned on.
>
> If we were to get rid of the size check
>
> if (v1.size() != v2.size())
> throw std::invalid_argument("v1.__size() must equal v2.size()");
>
> then the call
>
> return v1.dot(v2);
>
> will cause assertion failure and crash the program.
>
> I think we want to add the -NDEBUG flag back in for speed
> and control the error handling ourselves. For instance,
> we don't need to check in v.dot(v) for dot_self()!
>
>
> We probably want to leave those tests in, even with NDEBUG, at least for model compilation as it will catch errors in
> model specification and potentially produce useful error messages rather than crashing in unpredictable ways.

We're already trapping the main use of asserts in Eigen,
which are array indices out of bounds.

We definitely want to check at some point rather than
just seg fault.

- Bob

Ben Goodrich

unread,
Apr 23, 2012, 2:51:00 PM4/23/12
to stan...@googlegroups.com
On Saturday, April 21, 2012 3:58:34 PM UTC-4, Bob Carpenter wrote:
I'm writing the dot_self function to get a bit
more familiar with the templating

Got some problem here under clang ...

[----------] 1 test from MathMatrix
[ RUN      ] MathMatrix.dot_self
src/test/math/matrix_test.cpp:849: Failure
The difference between 13.0 and dot_self(v2) is 4, which exceeds 1E-12, where
13.0 evaluates to 13,
dot_self(v2) evaluates to 17, and
1E-12 evaluates to 9.9999999999999998e-13.
src/test/math/matrix_test.cpp:856: Failure
The difference between 4.0 and dot_self(rv1) is 29, which exceeds 1E-12, where
4.0 evaluates to 4,
dot_self(rv1) evaluates to 33, and
1E-12 evaluates to 9.9999999999999998e-13.
src/test/math/matrix_test.cpp:862: Failure
The difference between 29.0 and dot_self(rv3) is 13, which exceeds 1E-12, where
29.0 evaluates to 29,
dot_self(rv3) evaluates to 42, and
1E-12 evaluates to 9.9999999999999998e-13.
[  FAILED  ] MathMatrix.dot_self (0 ms)
[----------] 1 test from MathMatrix (0 ms total)

[----------] Global test environment tear-down
[==========] 58 tests from 4 test cases ran. (108 ms total)
[  PASSED  ] 57 tests.
[  FAILED  ] 1 test, listed below:
[  FAILED  ] MathMatrix.dot_self

 1 FAILED TEST
master!stan>

Ben

Bob Carpenter

unread,
Apr 23, 2012, 2:57:36 PM4/23/12
to stan...@googlegroups.com
Hmm. It all works here and presumably on Jenkins.

The function is very simple. I just changed
the Eigen [] access to () -- maybe that'll help
somehow.

- Bob

Daniel Lee

unread,
Apr 23, 2012, 3:07:50 PM4/23/12
to stan...@googlegroups.com
Interesting. I just got a warning regarding dot_self_vari:

In file included from src/test/agrad/distributions/multivariate/continuous/dirichlet_test.cpp:6:0:
src/stan/agrad/matrix.hpp: In static member function 'static double stan::agrad::dot_self_vari::var_dot_self(const Eigen::Matrix<stan::agrad::var, R, C>&)':
src/stan/agrad/matrix.hpp:498:7: warning: no return statement in function returning non-void [-Wreturn-type]


Indeed, there is no return statement. Could that be causing some of your issues?


Daniel

Ben Goodrich

unread,
Apr 23, 2012, 3:12:19 PM4/23/12
to stan...@googlegroups.com
On Monday, April 23, 2012 3:07:50 PM UTC-4, Daniel Lee wrote:
Interesting. I just got a warning regarding dot_self_vari:

In file included from src/test/agrad/distributions/multivariate/continuous/dirichlet_test.cpp:6:0:
src/stan/agrad/matrix.hpp: In static member function 'static double stan::agrad::dot_self_vari::var_dot_self(const Eigen::Matrix<stan::agrad::var, R, C>&)':
src/stan/agrad/matrix.hpp:498:7: warning: no return statement in function returning non-void [-Wreturn-type]


Indeed, there is no return statement. Could that be causing some of your issues?

I was just about to post the same thing. The test passes with g++-4.6 and fails with clang++ 3.1 prerelease.

Ben

Bob Carpenter

unread,
Apr 23, 2012, 7:16:44 PM4/23/12
to stan...@googlegroups.com
OK -- I"ll fix this. I was about to add
the unit tests for it anyway.

I can't believe the other compilers didn't
complain.

- Bob

Ben Goodrich

unread,
Apr 26, 2012, 1:50:21 AM4/26/12
to stan...@googlegroups.com
We've still got scary warnings and failing tests here when g++-4.7 or clang 3.1 prerelease, although g++-4.6 passes with the scary warnings:

In file included from src/test/math/matrix_test.cpp:4:0:
src/stan/math/matrix.hpp: In function ‘double stan::math::dot_self(const Eigen::Matrix<double, R, C>&) [with int R = -1; int C = -1]’:
src/stan/math/matrix.hpp:602:14: warning: ‘sum’ may be used uninitialized in this function [-Wmaybe-uninitialized]
src/stan/math/matrix.hpp: In member function ‘virtual void MathMatrix_dot_self_Test::TestBody()’:
src/stan/math/matrix.hpp:601:9: warning: ‘sum’ may be used uninitialized in this function [-Wuninitialized]
src/stan/math/matrix.hpp:599:14: note: ‘sum’ was declared here
src/stan/math/matrix.hpp:601:9: warning: ‘sum’ may be used uninitialized in this function [-Wuninitialized]
src/stan/math/matrix.hpp:599:14: note: ‘sum’ was declared here
src/test/math/matrix_test.cpp:862:3: warning: ‘sum’ may be used uninitialized in this function [-Wmaybe-uninitialized]
src/test/math/matrix_test.cpp:859:3: warning: ‘sum’ may be used uninitialized in this function [-Wmaybe-uninitialized]
src/test/math/matrix_test.cpp:852:3: warning: ‘sum’ may be used uninitialized in this function [-Wmaybe-uninitialized]
src/test/math/matrix_test.cpp:849:3: warning: ‘sum’ may be used uninitialized in this function [-Wmaybe-uninitialized]
make: *** [test/math/matrix] Error 1

goodrich@CYBERPOWERPC:/opt/stan$ tail -n 35 /tmp/g++47.log
src/test/math/matrix_test.cpp:866: Failure
The difference between 4.0 and dot_self(m1) is nan, which exceeds 1E-12, where
4.0 evaluates to 4,
dot_self(m1) evaluates to nan, and
1E-12 evaluates to 9.9999999999999998e-13.
src/test/math/matrix_test.cpp:869: Failure
The difference between 13.0 and dot_self(m2) is 3, which exceeds 1E-12, where
13.0 evaluates to 13,
dot_self(m2) evaluates to 16, and
1E-12 evaluates to 9.9999999999999998e-13.
src/test/math/matrix_test.cpp:872: Failure
The difference between 29.0 and dot_self(m3) is 4, which exceeds 1E-12, where
29.0 evaluates to 29,
dot_self(m3) evaluates to 33, and
1E-12 evaluates to 9.9999999999999998e-13.
src/test/math/matrix_test.cpp:876: Failure
The difference between 13.0 and dot_self(mm2) is 3, which exceeds 1E-12, where
13.0 evaluates to 13,
dot_self(mm2) evaluates to 16, and
1E-12 evaluates to 9.9999999999999998e-13.
src/test/math/matrix_test.cpp:879: Failure
The difference between 29.0 and dot_self(mm3) is 4, which exceeds 1E-12, where
29.0 evaluates to 29,
dot_self(mm3) evaluates to 33, and
1E-12 evaluates to 9.9999999999999998e-13.
[  FAILED  ] MathMatrix.dot_self (1 ms)
[----------] 1 test from MathMatrix (1 ms total)


[----------] Global test environment tear-down
[==========] 58 tests from 4 test cases ran. (74 ms total)

[  PASSED  ] 57 tests.
[  FAILED  ] 1 test, listed below:
[  FAILED  ] MathMatrix.dot_self

 1 FAILED TEST

Ben

Bob Carpenter

unread,
Apr 26, 2012, 12:37:28 PM4/26/12
to stan...@googlegroups.com
Thanks. I'm not sure why I'm not seeing those warnings,
but then I'm running 4.2 on the Mac.

On 4/26/12 1:50 AM, Ben Goodrich wrote:
> We've still got scary warnings and failing tests here when g++-4.7 or
> clang 3.1 prerelease, although g++-4.6 passes with the scary warnings:
>
> In file included from src/test/math/matrix_test.cpp:4:0:
> src/stan/math/matrix.hpp: In function ‘double stan::math::dot_self(const
> Eigen::Matrix<double, R, C>&) [with int R = -1; int C = -1]’:
> src/stan/math/matrix.hpp:602:14: warning: ‘sum’ may be used
> uninitialized in this function [-Wmaybe-uninitialized]
> src/stan/math/matrix.hpp: In member function ‘virtual void
> MathMatrix_dot_self_Test::TestBody()’:
> src/stan/math/matrix.hpp:601:9: warning: ‘sum’ may be used uninitialized

No maybe about it for a non-empty matrix. I
forgot to initialize.

It failed under g++ on my machine (but passed
under clang++).

I fixed it. Again, I'm amazed this works in
clang++ and on whatever Jenkins is running.

- Bob
Reply all
Reply to author
Forward
0 new messages