VectorXreal map

[kristianmeyerr]

Hi,

I have data stored in a std::vector, and would like to use the data with the autodiff::VectorXreal type. I wonder if there is a better way than copying it like this:

std::vector<double> x{1, 2, 3}; 
VectorXreal xr(3);
for (unsigned int i = 0; i < x.size(); ++i) {
       xr(i) = x[i];
}

When using Eigen i typically use a map:

Eigen::Map<Eigen::Matrix<double,1,3> > x_map(x.data());

Can we do somethings similar with the VectorXreal class ?

Looking forward.

[allanleal]

Hi Kristian, please try and let me know if it does not work. If it does not work, please provide a minimal reproducible example.

[kristianmeyerr]

Hi Allan,

The above code snippet works, the problem is that it copies the content from x into xr - I wonder if there is a way to construct a VectorXreal using a pointer? Like what you typically do with Eigen::Map?

Best regards Kristian

[jeanchristopheruel]

Here is a minimal reproducible example of the desirable behavior

auto f = [](Eigen::Matrix<autodiff::real, 1, 8> const& x)
{
    return (x * x.transpose()).sum();
};
Eigen::Matrix<autodiff::real, 2, 4> xx({
        {1, 2, 3, 4},
        {5, 6, 7, 8}
});

Eigen::Map<const Eigen::Matrix<autodiff::real, 1, 8>> x(xx.data(), xx.size());
autodiff::real u;  

Eigen::VectorXd gx = gradient(f, wrt(x), at(x), u); 

and the compile error:

_deps/autodiff-src/autodiff/forward/utils/derivative.hpp:96:24: error: static assertion failed: It is not possible to compute higher-order derivatives with order greater than that of the autodiff number (e.g., using wrt(x, x, y, z) will fail if the autodiff numbers in use have order below 4).
   96 |     static_assert(size <= N, "It is not possible to compute higher-order derivatives with order greater than that of the autodiff number (e.g., using wrt(x, x, y, z) will fail if the autodiff numbers in use have order below 4).");
      |                   ~~~~~^~~~

[ipcamit]

I was also wondering the same. Instead of having element wise copy, such eigen map will be lot more convenient

[jeanchristopheruel]

Seems like this is caused by an issue with the Template parameters counting here. I guess we should find a way of subtracting the additional template parameter added by Eigen::Map. Furthermore, we should find a cleaner method than counting all template Args as if they were all related to autodiff number.

Would also be convenient to expose Eigen::Ref in autodiff::gradient and autodiff::jacobian interfaces as it would allow to use a temporary matrix like a block.

[allanleal]

Hi @jeanchristopheruel , thanks for providing a reproducible example and making suggestions above. I've created the PR above that will finally provide support for the use of Eigen::Map together with Eigen::Vector/Matrix of either autodiff::real or autodiff::dual.

Note that for your example to work, you'll need to remove the constness of your map for x, otherwise it is not possible to seed these numbers for the sake of automatic differentiation:

Replace:

Eigen::Map<const Eigen::Matrix<autodiff::real, 1, 8>> x(xx.data(), xx.size());

by

Eigen::Map<Eigen::Matrix<autodiff::real, 1, 8>> x(xx.data(), xx.size());

The following example works now:

#include <iostream>
#include <autodiff/forward/real/real.hpp>
#include <autodiff/forward/real/eigen.hpp>

using namespace autodiff;
using std::cout, std::endl;

int main(int argc, char const *argv[])
{
    auto f = [](Eigen::Matrix<autodiff::real, 1, 8> const& x)
    {
        return (x * x.transpose()).sum();
    };

    Eigen::Matrix<autodiff::real, 2, 4> xx({
            {1, 2, 3, 4},
            {5, 6, 7, 8}
    });

    Eigen::Map<Eigen::Matrix<autodiff::real, 1, 8>> x(xx.data(), xx.size());
    autodiff::real u;

    Eigen::VectorXd gx = gradient(f, wrt(x), at(x), u);

    cout << gx << endl;

    return 0;
}

// Output:
//  2
// 10
//  4
// 12
//  6
// 14
//  8
// 16

[jeanchristopheruel]

This is a big deal, thanks @allanleal ! 👌

[ipcamit]

@allanleal Just one more clarification, can I use Map on double * ? I mean what is the most efficient way for going from double * memory block to VectorXreal? Only way I could think of was

double x[] = {1., 2., 3.};
Eigen::Map<VectorXd> x_vec(x, 3);
auto x_real = static_cast<VectorXreal>(x_vec); 

But it involves making copy of the data, which I want to avoid.

[allanleal]

Hi @ipcamit , I suppose you'll want to compute derivatives of a function with respect to your x data above, currently available in double type. You'll need to convert these double values to an autodiff type number (try real). In addition, you'll want to implement the function using this number type. Otherwise, there is no introspection within your function to compute derivatives automatically. These autodiff numbers make use of chain rules of derivatives during mathematical operations (e.g., +, -, *, /, sin, cos, etc.) to compute not only the actual value of functions (and expressions) but also their derivatives.

[ipcamit]

My actual problem is to port an existing library to have gradients (currently it uses numerical gradients). For that I would like to keep the interface for calculations same. Hence the idea was to accept double *, exactly as it does right now, and convert them to vectorXreal in the library and pass on to new internal function of type real. So in that case x in above example will be an input that I dont have much control over. So that is why I was wondering, that what is the fastest and most efficient way to convert a double * array to real array. It seems like there is no option other then copying the data.

[allanleal]

@ipcamit , there is no way out of this; you'll need to create real objects and copy the double data to them. I would bet, though, that this memory copy will be the least of your concern in terms of performance.