Subtraction

Developers
1

If I do this,

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(StanHeaders)]]
#include <Rcpp.h>
#include <RcppEigen.h>
#include <stan/math/mix/mat/functor/hessian.hpp>
#include <stan/math/fwd/mat/fun/dot_self.hpp>
#include <stan/math.hpp> 

// [[Rcpp::export]]
Eigen::MatrixXd H(Eigen::VectorXd x, Eigen::VectorXd a) { // Hessian by AD using Stan
  double fx;
  Eigen::VectorXd grad_fx;
  Eigen::MatrixXd H;
  stan::math::hessian([&a](auto x) {
                        auto z = x.eval();
                        for (int i = 0; i < x.rows(); i++) z.coeffRef(i) -= a.coeff(i);
                        return stan::math::dot_self(z);
                      }, x, fx, grad_fx, H);
  return H;
}

and call Rcpp::sourceCpp on that file, then I can call H() from R. But if I change the function definition to

// [[Rcpp::export]]
Eigen::MatrixXd H(Eigen::VectorXd x, Eigen::VectorXd a) { // Hessian by AD using Stan
  double fx;
  Eigen::VectorXd grad_fx;
  Eigen::MatrixXd H;
  stan::math::hessian([&a](auto x) {
    return stan::math::dot_self( (x - a).eval() );
  }, x, fx, grad_fx, H);
  return H;
}

it fails with

/home/ben/R/x86_64-pc-linux-gnu-library/3.3/RcppEigen/include/Eigen/src/Core/CwiseBinaryOp.h:107:7: error: static_assert failed due to requirement ‘Eigen::internal::has_ReturnType<ScalarBinaryOpTraits<typename Lhs::Scalar, typename Rhs::Scalar, scalar_difference_op<fvar, double> > >::value’ “YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY”
EIGEN_CHECK_BINARY_COMPATIBILIY(BinaryOp,typename Lhs::Scalar,typename Rhs::Scalar);

Why is subtracting a double vector from a fvar<double> vector illegal while subtracting a double from a fvar<double> is legal?

2

Can you try by first importing Stan/math/mix/mat.hpp?

All bets are off if you include in the order you’ve included things in. There might be problems if you’ve included Eigen before this file too, but maybe including the right header would fix it.

3

I am getting the same error when the top of the file is

// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(StanHeaders)]]
#include <stan/math/mix/mat.hpp>
#include <Rcpp.h>
#include <RcppEigen.h>
#include <stan/math/fwd/mat/fun/dot_self.hpp>
#include <stan/math.hpp>
4

I’ll try to see if I can reproduce later this weekend.

5

I can definitely reproduce the problem under develop. I don’t know if this ever worked, but right now, it doesn’t look like you can subtract a Eigen Vector of fvar<double> from double. This currently doesn’t work (including either stan/math/fwd/mat.hpp or stan/math/mix/mat.hpp):

  Eigen::VectorXd x(5);
  x << 1, 2, 3, 4, 5;
  Eigen::Matrix<stan::math::fvar<double>, -1, 1> fd(5);
  fd << 4, 5, 6, 7, 8;
  std::cout << fd - x << std::endl;

Has this particular function worked in the past for you?

I wonder if we don’t have the Eigen NumTraits specified correctly for forward mode. (We have specializations for both reverse and forward mode.)

6

No, I have never known that to work but I was attempting to test more things in the StanHeaders vignette.

1 Like
7

@bgoodri, want to help build out forward mode?

@Davor_Sluga: if you’re really getting comfortable with forward mode, this is something that should be looked into. I don’t have the slightest clue where to start, though.

8

We implement operator- for all the scalar types in fvar/core. The error message you’re getting is when Eigen tries to mix types, which it won’t allow (at least not without some more traits hints from us).

9

Oh, forgot the fix. You can use add(x, y) or subtract(x, y), but they do inefficient promotion. See the implementations in the math repo in stan/math/prim/mat/fun/.

10

OK, that works. But I’m still not sure why you can just use the subtraction sign when subtracting a double vector from a var vector but have to use the subtract function when subtracting a double vector from a fvar<var> vector.

1 Like
11

I bet we don’t have the traits defined for Eigen properly for forward mode. I couldn’t get it working for fvar<double> and double either.

12

Right. You need to define the traits to allow mixing scalar types to define what the return types are.

1 Like