Multivariate Function with Known Gradients - RStan

Developers
1

Hi,

As in the issue Custom Function with Known Gradients, I got a function to work which maps double -> double. Eventually I want to use a function that maps vector -> vector, but for now I have difficulties to implement a function that maps vector -> double and is callable from a Stan file.

The function is declared as

namespace tests {

    double vecToDouble(const std::vector<double>& x)
    {
        double res = 0;
        for (int i=0; i<x.size(); ++i) res += x[i] * x[i];
        return res;
    }

    std::vector<double> grad_vecToDouble(const std::vector<double>& x)
    {
        std::vector<double> res;
        for (int i=0; i<x.size(); ++i) res.push_back(2 * x[i]);
        return res;
    }
}

and for the Stan interface I currently have

namespace stan {
    namespace math {

        double vecToDouble(const std::vector<double>& x) {
            return tests::vecToDouble(x);
        }

        var vecToDouble(const Eigen::Matrix<var, Eigen::Dynamic, 1>& x) {
            std::vector<var> b(x.data(), x.data() + x.rows() * x.cols());

            std::vector<double> a = value_of(b);
            double fa  = tests::vecToDouble(a);
            std::vector<double> grad_fa = tests::grad_vecToDouble(a);

            return precomputed_gradients(fa, b, grad_fa);
        }

    }
}

It works in pure C++ as in

int main() {

    int N = 3;

    Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1> x(N);
    std::vector<double> xx;
    for (int i=0; i<N; ++i)
    {
        stan::math::var xi = (double) i;
        x(i, 0)  = xi;
        xx.push_back((double) i);
    }

    printf("function value: %.2f\n", stan::math::vecToDouble(xx));

    stan::math::var y = stan::math::vecToDouble(x);
    printf("vecToDouble(x) = %.2f\n", y.val());

    std::vector<double> g;
    std::vector<stan::math::var> in(x.data(), x.data() + x.rows() * x.cols());
    y.grad(in, g);
    for (int i=0; i<g.size(); ++i)
    {
        printf("g[%d] = %.2f\n", i, g[i]);
    }

    return 0;
}

However, once I want to call it from Stan it leads to this error

…/Library/R/3.4/library/RcppEigen/include/Eigen/src/Core/functors/AssignmentFunctors.h:24:104: error: assigning to ‘double’ from incompatible type 'const stan::math::var’
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void assignCoeff(DstScalar& a, const SrcScalar& b) const { a = b; }

and a couple of errors follow after that which seem to be related to this error.

As a side note, I first declared the Stan interface as described in https://github.com/stan-dev/math/wiki/Adding-a-new-function-with-known-gradients like:

namespace stan {
    namespace math {

        double vecToDouble(const std::vector<double>& x);

        var vecToDouble(const std::vector<var>& x);

    }
}

But this lead to an error saying that the function vecToDouble wasn’t declared. I could see in the generated Cpp file that the variable passed into the function is declared as Eigen::Matrix<stan::math::var, Eigen::Dynamic, 1>, so I changed the types and the function can apparently be found.

Any help will be greatly appreciated!

Regards,
Jonas

Operating System: OS X 10.12
Interface Version: RStan 2.16.2
Compiler/Toolkit: clang-900.0.38

2

The solution to this issue is quite simple as I got know in the thread Barriers to New Developers. Basically, one should not expose a vector-valued function based on std::vector, but always based on Eigen::Matrix. Additionally, there always has to be an overloaded function definition for both double and var.

namespace stan {
    namespace math {

        inline double
        vecToDoubleTwo(const Eigen::Matrix<double, Eigen::Dynamic, 1>& m) {

            double val = 0;
            for (int i=0; i<m.size(); ++i) val += m(i) * m(i);

            return val;
        }


        inline var
        vecToDoubleTwo(const Eigen::Matrix<var, Eigen::Dynamic, 1>& m) {
            using Eigen::Matrix;

            Matrix<double, Eigen::Dynamic, 1> m_d(m.rows());
            for (int i = 0; i < m.size(); ++i)
                m_d(i) = m(i).val();

            double val = stan::math::vecToDouble(m_d);

            vari** varis
            = ChainableStack::memalloc_.alloc_array<vari*>(m.size());
            for (int i = 0; i < m.size(); ++i)
                varis[i] = m(i).vi_;

            double* gradients
            = ChainableStack::memalloc_.alloc_array<double>(m.size());
            for (int i = 0; i < m.size(); ++i)
                gradients[i] = 9.9 * m_d(i);

            return var(new precomputed_gradients_vari(val, m.size(),
                                                      varis, gradients));
        }

    }
}
3

The alternative to providing overloads for var and double is to write the function templated in a way that it “does the right” thing automagically depending on the types passed into the function. For these types of functions, the stan-math library has invented the operands_and_partials thing. It is quite heavy in the use of template stuff, but I figured it out by looking at how it is used in the stan-math library. For example the definition of the normal lpdf density can serve as an example. The code gets a lot cleaner when using this.

2 Likes
4

I’ll have a look at that. Thanks!