higher-order derivatives in prob/transform.hpp

[jpritikin]

I re-generated the C++ code for a simple model (attached) and reduced it to a minimum test case (also attached). I still get compiler errors expanding the templates.

See https://groups.google.com/forum/#!topic/stan-users/AGDRGUoqjBE

Note: The problem is in the transforms.

[rtrangucci]

@jpritikin Your .stan file compiles fine for me in the develop branches of cmdstan and stan. It doesn't look like the .cpp file was generated by the latest cmdstan.

I'm going to close this issue

[jpritikin]

Yes, of course the .stan file compiles fine. The point of this issue is to obtain 1st and 2nd derivatives for use in another application, OpenMx. I want to use stan as a library, much the same way that stan uses Eigen as a library.

Specifically, I wrote the main() function myself. If you comment out stan::agrad::hessian and uncomment stan::agrad::gradient then you will find that the code works fine and obtains the gradient. I wish to obtain both the hessian and the gradient. Can you tell me how to do this?

[jpritikin]

Here is an earlier conversation with Bob Carpenter that can give you some context, https://groups.google.com/forum/#!searchin/stan-users/pritikin/stan-users/S6Ss-49PVd8/cbq-C-DSftcJ

[bob-carpenter]

As far as I can see, this code should compile. It looks like the main problem is this:

/opt/cmdstan/stan/src/stan/prob/transform.hpp:1784:49: error: call of overloaded ‘multiply_lower_tri_self_transpose(Eigen::Matrix<stan::agrad::fvar<stan::agrad::var>, -1, -1, 0, -1, -1>&)’ is ambiguous 
       return multiply_lower_tri_self_transpose(L); 

I think this problem can be solved with some disable_if templates. Or maybe the double-based one should really just be defined for double and not templated.

[bob-carpenter]

P.S. We probably never tested the transforms in src/prob/transform.hpp for higher-order types.

P.P.S. We also need to vectorize the autodiff in transform.hpp --- maybe these need to be tackled at the same time.

[rtrangucci]

@bob-carpenter, agreed. The offending function in transform.hpp is the cov_matrix_constrain function, and while the cov_matrix_constrain function is templated, we haven't included the fvar version of multiply_lower_tri_self_transpose(). The math/multiply_lower_tri_self_transpose() isn't templated, there's a templated version in fwd/matrix/multiply_lower_tri_self_transpose.

[rtrangucci]

Also, @jpritikin #676 was meant to generalize OperandsAndPartials for higher order derivatives and propagate the generalization through to the math functions and distributions. We haven't yet generalized the any other functions for higher order derivatives, although that is certainly on our to-do list.

[jpritikin]

Good to know.

[rtrangucci]

@bob-carpenter @syclik do we want to modify transform.hpp to allow for higher order derivates in 2.6? If not, I think we should move this issue to the milestone 2.6.0++ rather than 2.6.

[syclik]

I'm thinking we take care of it before releasing 2.6. I'd rather have the next tagged version have a fully functioning higher order autodiff instead of a mostly working higher order autodiff.

On Fri, Dec 19, 2014 at 6:23 PM, Rob Trangucci notifications@github.com wrote:

@bob-carpenter https://github.com/bob-carpenter @syclik https://github.com/syclik do we want to modify transform.hpp to allow for higher order derivates in 2.6? If not, I think we should move this issue to the milestone 2.6.0++ rather than 2.6.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/1175#issuecomment-67711964.

[bob-carpenter]

Agreed. I'm not thinking we'll get this release out before New Years. The upside is that it'll be an awesome release!

I'm about halfway through writing up reverse-mode autodiff from a client and developer perspective. We then need to add performance results. We could also tack on higher-order autodiff, but the performance doesn't look so hot comparatively.

I think rewriting so that the value and tangent are separately parameterized would help immensely --- in all the higher order uses we're building redundant expression graph components due to the values being autodiff variables where we only need the tangents to be.

template <typename V, typename T> struct fvar { V val; T d; };

So we'd use fvar<double,rev> for second order.

To get away with fvar<fvar<double,double>, fvar<double,rev> > for third order, we'd need to allow mixed type multiplications, in sketch form:

fvar<promote<T1,T2>::type, promote<T3,T4>::type> operator*(const fvar<T1,T2>&, const fvar<T3,T4>&)

• Bob

On Dec 19, 2014, at 6:26 PM, Daniel Lee notifications@github.com wrote:

I'm thinking we take care of it before releasing 2.6. I'd rather have the next tagged version have a fully functioning higher order autodiff instead of a mostly working higher order autodiff.

On Fri, Dec 19, 2014 at 6:23 PM, Rob Trangucci notifications@github.com wrote:

@bob-carpenter https://github.com/bob-carpenter @syclik https://github.com/syclik do we want to modify transform.hpp to allow for higher order derivates in 2.6? If not, I think we should move this issue to the milestone 2.6.0++ rather than 2.6.

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/1175#issuecomment-67711964.

— Reply to this email directly or view it on GitHub.

[syclik]

@jpritikin, I'm closing the issue. This isn't a bug. We've separated out the header files so that by default you don't have to include forward mode with a typical Stan program.

To fix your example, you need to add two lines to the very top. Your includes should now look like this:

#include <stan/agrad/fwd.hpp>
#include <stan/agrad/fwd/matrix.hpp>
#include <stan/model/model_header.hpp>

By including that, I was able to run your program (I couldn't get your script to work, but I cobbled something together). The output looks like:

> time ./mvn
lp -12.4637
grad:
   0.0131862
-0.000858941
     10.2412
 -0.00213288
     2.08381
hess:
    -33.0292     -0.38292   -0.0263979   0.00104714  4.33294e-05
    -0.38292     -15.1816 -0.000153168  -0.00900987   0.00171788
  -0.0263979 -0.000153168     -44.4889    -0.380546  7.25181e-05
  0.00104714  -0.00900987    -0.380546     -22.5105   0.00426577
 4.33294e-05   0.00171788  7.25181e-05   0.00426577     -30.1676

real    0m0.006s
user    0m0.002s
sys 0m0.002s

[syclik]

@rtrangucci @bob-carpenter: the higher order is working. We just didn't include the forward mode headers into the model header for all models.

[jpritikin]

That's awesome. I added the headers mentioned, but I still get a compile error compiling against stan dc61f6fb4508b3b6e10 (develop). Should I be using a different branch?

In file included from /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/Core:259:0, from /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/Dense:1, from /opt/cmdstan/stan/src/stan/math/matrix/Eigen.hpp:4, from /opt/cmdstan/stan/src/stan/meta/traits.hpp:13, from /opt/cmdstan/stan/src/stan/agrad/fwd/numeric_limits.hpp:5, from /opt/cmdstan/stan/src/stan/agrad/fwd.hpp:6, from mvn.cpp:1: /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/MathFunctions.h: In instantiation of ‘static NewType Eigen::internal::cast_impl<OldType, NewType>::run(const OldType&) [with OldType = stan::agrad::fvarstan::agrad::var; NewType = int]’: /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/MathFunctions.h:328:44: required from ‘NewType Eigen::internal::cast(const OldType&) [with OldType = stan::agrad::fvarstan::agrad::var; NewType = int]’ /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/IO.h:132:97: required from ‘static int Eigen::internal::significant_decimals_default_impl<Scalar, IsInteger>::run() [with Scalar = stan::agrad::fvarstan::agrad::var; bool IsInteger = false]’ /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/IO.h:180:67: required from ‘std::ostream& Eigen::internal::print_matrix(std::ostream&, const Derived&, const Eigen::IOFormat&) [with Derived = Eigen::Matrixstan::agrad::fvar<stan::agrad::var, -1, 1>; std::ostream = std::basic_ostream]’ /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/IO.h:245:69: required from ‘std::ostream& Eigen::operator<<(std::ostream&, const Eigen::DenseBase&) [with Derived = Eigen::Matrixstan::agrad::fvar<stan::agrad::var, -1, 1>; std::ostream = std::basic_ostream]’ mvn.cpp:70:22: required from ‘T mvn_model_namespace::mvn_model::log_prob(std::**debug::vector&, std::debug::vector&, std::ostream*) const [with bool propto* = true; bool jacobian = true; T = stan::agrad::fvarstan::agrad::var; std::ostream = std::basicostream]’ mvn.cpp:107:78: required from ‘T mvn_model_namespace::mvn_model::log_prob(Eigen::Matrix<T, -1, 1>&, std::ostream) const [with bool propto = true; bool jacobian = true; T_ = stan::agrad::fvarstan::agrad::var; std::ostream = std::basic_ostream]’ mvn.cpp:124:53: required from ‘T model_functional::operator()(Eigen::Matrix<T, -1, 1>&) const [with T = stan::agrad::fvarstan::agrad::var; M = mvn_model_namespace::mvn_model]’ /opt/cmdstan/stan/src/stan/agrad/autodiff.hpp:207:39: required from ‘void stan::agrad::hessian(const F&, const Eigen::Matrix<double, -1, 1>&, double&, Eigen::Matrix<double, -1, 1>&, Eigen::Matrix<double, -1, -1>&) [with F = model_functional]’ mvn.cpp:152:58: required from here /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/MathFunctions.h:319:34: error: invalid static_cast from type ‘const stan::agrad::fvarstan::agrad::var’ to type ‘int’ return static_cast(x); ^

[syclik]

Did you also update CmdStan to develop? (It looks like you're behind, but shouldn't be a problem.)

I couldn't get your script to work, so here's what I did: 1) copy your main into CmdStan's src/cmdstan/main.cpp. 2) added the fvar includes to the top of CmdStan's src/cmdstan/main.cpp 3) regenerated the mvn.cpp by using CmdStan: make mvn.cpp (from within CmdStan) 4) added the custom constructor that you had in the original mvn.cpp file 5) generated the executable: make mvn (from within CmdStan)

Then I was able to run the executable.

I'm not sure what your other flags did. Also, it looks like some of your types are messed up. We've never created a 'stan::agrad::fvarstan::agrad::var' type.

On Fri, Dec 26, 2014 at 10:03 AM, Joshua Pritikin notifications@github.com wrote:

That's awesome. I added the headers mentioned, but I still get a compile error compiling against stan dc61f6f https://github.com/stan-dev/stan/commit/dc61f6fb4508b3b6e10a7bc6e5d8966b248a3330 (develop). Should I be using a different branch?

In file included from /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/Core:259:0, from /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/Dense:1, from /opt/cmdstan/stan/src/stan/math/matrix/Eigen.hpp:4, from /opt/cmdstan/stan/src/stan/meta/traits.hpp:13, from /opt/cmdstan/stan/src/stan/agrad/fwd/numeric_limits.hpp:5, from /opt/cmdstan/stan/src/stan/agrad/fwd.hpp:6, from mvn.cpp:1: /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/MathFunctions.h: In instantiation of ‘static NewType Eigen::internal::cast_impl::run(const OldType&) [with OldType = stan::agrad::fvarstan::agrad::var; NewType = int]’: /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/MathFunctions.h:328:44: required from ‘NewType Eigen::internal::cast(const OldType&) [with OldType = stan::agrad::fvarstan::agrad::var; NewType = int]’ /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/IO.h:132:97: required from ‘static int Eigen::internal::significant_decimals_default_impl::run() [with Scalar = stan::agrad::fvarstan::agrad::var; bool IsInteger = false]’ /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/IO.h:180:67: required from ‘std::ostream& Eigen::internal::print_matrix(std::ostream&, const Derived&, const Eigen::IOFormat&) [with Derived = Eigen::Matrixstan::agrad::fvar<stan::agrad::var, -1, 1>; std::ostream = std::basic_ostream]’ /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/IO.h:245:69: required from ‘std::ostream& Eigen::operator<<(std::ostream&, const Eigen::DenseBase&) [with Derived = Eigen::Matrixstan::agrad::fvar<stan::agrad::var, -1, 1>; std::ostream = std::basic_ostream]’ mvn.cpp:70:22: required from ‘T mvn_model_namespace::mvn_model::log_prob(std::_debug::vector&, std::_debug::vector&, std::ostream) const [with bool propto* = true; bool jacobian = true; T__ = stan::agrad::fvarstan::agrad::var; std::ostream = std::basicostream]’ mvn.cpp:107:78: required from ‘T mvn_model_namespace::mvn_model::logprob(Eigen::Matrix&, std::ostream*) const [with bool propto = true; bool jacobian = true; T = stan::agrad::fvarstan::agrad::var; std::ostream = std::basic_ostream]’ mvn.cpp:124:53: required from ‘T model_functional::operator()(Eigen::Matrix&) const [with T = stan::agrad::fvarstan::agrad::var; M = mvn_model_namespace::mvn_model]’ /opt/cmdstan/stan/src/stan/agrad/autodiff.hpp:207:39: required from ‘void stan::agrad::hessian(const F&, const Eigen::Matrix&, double&, Eigen::Matrix&, Eigen::Matrix&) [with F = model_functional]’ mvn.cpp:152:58: required from here /opt/cmdstan/stan/lib/eigen_3.2.2/Eigen/src/Core/MathFunctions.h:319:34: error: invalid static_cast from type ‘const stan::agrad::fvarstan::agrad::var’ to type ‘int’ return static_cast(x); ^

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/1175#issuecomment-68144139.

[jpritikin]

It wasn't quite as easy as you suggested, but I did get it working somehow. Thanks for the hints!

[syclik]

Odd. Those were the steps I took.

What did you end up doing? On Dec 26, 2014 4:38 PM, "Joshua Pritikin" notifications@github.com wrote:

It wasn't quite as easy as you suggested, but I did get it working somehow. Thanks for the hints!

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/1175#issuecomment-68159282.

[jpritikin]

On Fri, Dec 26, 2014 at 01:39:52PM -0800, Daniel Lee wrote:

Odd. Those were the steps I took. What did you end up doing?

It wasn't entirely clear where to paste what. Attached is a copy of my main.cpp and mvn.cpp file.

Joshua N. Pritikin Department of Psychology University of Virginia 485 McCormick Rd, Gilmer Hall Room 102 Charlottesville, VA 22904 http://people.virginia.edu/~jnp3bc

[syclik]

Ah. My bad.

Glad you figured it out. On Dec 26, 2014 4:45 PM, "Joshua Pritikin" notifications@github.com wrote:

On Fri, Dec 26, 2014 at 01:39:52PM -0800, Daniel Lee wrote:

Odd. Those were the steps I took. What did you end up doing?

It wasn't entirely clear where to paste what. Attached is a copy of my main.cpp and mvn.cpp file.

Joshua N. Pritikin Department of Psychology University of Virginia 485 McCormick Rd, Gilmer Hall Room 102 Charlottesville, VA 22904 http://people.virginia.edu/~jnp3bc

— Reply to this email directly or view it on GitHub https://github.com/stan-dev/stan/issues/1175#issuecomment-68159530.