Boost's uBLAS matrix library does not work with autodiff (example attached)

[Daniel Lee]

In an earlier email, Bob outlined the problem:

The issue with lifespan is that Boost's ublas library for matrices defines static class constants everywhere by wrapping the scalar type. And so does the rosenbrock functions.  This won't work with our variable lifespans. We'd need traits more like those in Eigen, which seem to be modeled on the function-based approach of std::numeric_limits.  I was sort of surprised to see Boost using static constants in ublas --- it's usually bad ju-ju in C++ because of initialization issues.  See:

  /boost/numeric/ublas/matrix.hpp:

  ----------------------------
  ...
  static const value_type one_;

The result of uBLAS's use of static constants is that we get a segmentation fault. The attached test, as written, does not seg fault. By accessing any element in the identity matrix (lines 9-12), it will seg fault.

No error messages, no nothing. It all happens before this code executes. It took us a while to figure out how and where the seg fault was coming from, but we found it.

In short, Stan's autodiff library will not work with Boost's uBLAS matrix library.



Daniel

P.S. This means Boost numeric's ODE solvers will not work with Stan. There are other issues to getting those ODE solvers to work.

[Ben Goodrich]

What exactly is the code to compile and link this thing?

[Daniel Lee]

Just save it into src/test/ublas_test.cpp.

Then:
make test/ublas

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

[Ben Goodrich]

Works for me

clang++ -I src -I lib/eigen_3.2.0 -I lib/boost_1.54.0 -Wall -DBOOST_RESULT_OF_USE_TR1 -DBOOST_NO_DECLTYPE -DBOOST_DISABLE_ASSERTS -Wno-unused-function -Wno-uninitialized -Wno-c++11-long-long    -O0 lib/gtest_1.6.0/src/gtest_main.cc test/ublas.o -DGTEST_USE_OWN_TR1_TUPLE -DGTEST_HAS_PTHREAD=0 -Wno-c++11-long-long -I lib/gtest_1.6.0/include -I lib/gtest_1.6.0 -o test/ublas test/libgtest.a -Lbin -lstan -Lbin -lstanc
test/ublas --gtest_output="xml:test/ublas.xml"
Running main() from gtest_main.cc
[==========] Running 3 tests from 1 test case.
[----------] Global test environment set-up.
[----------] 3 tests from ublas
[ RUN      ] ublas.identity_matrix
[       OK ] ublas.identity_matrix (0 ms)
[ RUN      ] ublas.matrix_uninitialized
[       OK ] ublas.matrix_uninitialized (0 ms)
[ RUN      ] ublas.matrix_initialized
[       OK ] ublas.matrix_initialized (0 ms)
[----------] 3 tests from ublas (0 ms total)

[----------] Global test environment tear-down
[==========] 3 tests from 1 test case ran. (0 ms total)
[  PASSED  ] 3 tests.

[Daniel Lee]

Yes, now comment out line 9. 

[Ben Goodrich]

OK, got it now. Looks like unallocated memory. From gdb,

Program received signal SIGSEGV, Segmentation fault.
0x00000000004072d4 in stan::memory::stack_alloc::move_to_next_block (this=0x669a60 <stan::agrad::memalloc_>, len=24) at src/stan/memory/stack_alloc.hpp:104
104               size_t newsize = sizes_.back() * 2;
(gdb) bt
#0  0x00000000004072d4 in stan::memory::stack_alloc::move_to_next_block (this=0x669a60 <stan::agrad::memalloc_>, len=24) at src/stan/memory/stack_alloc.hpp:104
#1  0x00000000004075d8 in stan::memory::stack_alloc::alloc (this=0x669a60 <stan::agrad::memalloc_>, len=24) at src/stan/memory/stack_alloc.hpp:172
#2  0x00000000004078db in stan::agrad::chainable::operator new (nbytes=24) at src/stan/agrad/rev/chainable.hpp:64
#3  0x0000000000407bf6 in stan::agrad::var::var (this=0x6698c0 <boost::numeric::ublas::identity_matrix<stan::agrad::var, std::allocator<stan::agrad::var> >::one_>,
    n=1) at src/stan/agrad/rev/var.hpp:122
#4  0x00000000004065e8 in __static_initialization_and_destruction_0 (__initialize_p=1, __priority=65535) at lib/boost_1.54.0/boost/numeric/ublas/matrix.hpp:2974
#5  0x0000000000406688 in _GLOBAL__sub_I__ZN26ublas_identity_matrix_Test10test_info_E () at src/test/ublas_test.cpp:39
#6  0x000000000044490d in __libc_csu_init ()
#7  0x00007ffff408f925 in __libc_start_main (main=0x404b1d <main(int, char**)>, argc=1, ubp_av=0x7fffffffe248, init=0x4448b0 <__libc_csu_init>, fini=<optimized out>,
    rtld_fini=<optimized out>, stack_end=0x7fffffffe238) at libc-start.c:219
#8  0x0000000000404a53 in _start ()

Ben

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Bob Carpenter]

Exactly. The problem is with static allocation.

I think we could fix this by replacing the static constant
one_ with a call to a static function one(). That's
how a traits-based approach would do it, and that will work
for us as long as nobody calls the function out of turn.
I assume that's what they were talking about in the ublas
package in their comment saying they should use traits
in the vector header file.

It doesn't make sense conceptually in our reverse-mode
auto-diff to have a free-floating value of type
stan::agrad::var. The use pattern is:

-- create var instances for inputs to functions
-- call the function on those instances
-- calculate a gradient in a backward pass
-- free memory

And to boot, it's not thread safe (unless you make
the variable stack thread local, and then there are
performance penalties).

The alternative is to somehow manage multiple memory arenas,
with references from var to the relevant arena. The problem
with that is that the objects are bigger and there's
a dynamic indirection that's not cheap to follow in
the guts of our loop.

Another alternative would be to have long-lived var
instances. They wouldn't need any auto-diff because
they wouldn't depend on any var values. The problem is
that I can't figure out how to do that without
making the value and derivative calls virtual, which
I would also rather not do for speed reasons.

Finally, there's the same mixed-type issue as with Eigen.
They wrap the constants for the algorithm up in
member variable instances of the scalar type, which means
more auto-diff variables that should really be double for us.

What this means is that we don't want the one_ value
to even be a variable --- all we need is a double if we
have double by var matrix multiplication. And seriously,
we don't even need that --- it's probably expression templated
out somewhere because it's a constant!

So it seemed like the easiest thing to do to just
rewrite their solver calling Eigen instead of ublas
and using double instances instead of var types for
constants.

- Bob

> /boost/numeric/ublas/matrix.__hp__p:

> ----------------------------
> ...
> static const value_type one_;
>
>
>
>
> The result of uBLAS's use of static constants is that we get a segmentation fault. The attached
> test, as written, does not seg fault. By accessing any element in the identity matrix (lines 9-12),
> it will seg fault.
>
> No error messages, no nothing. It all happens before this code executes. It took us a while to
> figure out how and where the seg fault was coming from, but we found it.
>
> In short, Stan's autodiff library will not work with Boost's uBLAS matrix library.
>
>
>
> Daniel
>
> P.S. This means Boost numeric's ODE solvers will not work with Stan. There are other issues to
> getting those ODE solvers to work.
>
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing
> list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to

> stan-dev+u...@__googlegroups.com.
> For more options, visit https://groups.google.com/__groups/opt_out
> <https://groups.google.com/groups/opt_out>.

>
>
> --
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to
> stan-dev+u...@googlegroups.com.

> For more options, visit https://groups.google.com/groups/opt_out <https://groups.google.com/groups/opt_out>.