构建编译时矩阵表达式

Question

我刚开始使用Eigen，并从他们的文档中了解到，最佳性能来自对矩阵表达式的惰性评估。因此，这样的表达式一经评估就非常有效：

Eigen::Matrix<float, 3, 1> a;
a << 0, 1, 2;

Eigen::Matrix<float, 3, 1> b;
b << 3, 4, 5;

Eigen::Matrix<float, 3, 1> c;
c << (a + b).sum(),
     (a - b).sum(),
     a.sum();

std::cout << c << std::endl;

在构造其列数取决于模板参数的矩阵时遇到问题。例如：

template <std::size_t w>
auto buildGradient() {
  Eigen::Matrix<float, 3, w> matrix;
  matrix << /* ? */;
  return matrix;
}

我的初衷是使用递归c ++模板来做到这一点。

template <std::size_t w, typename Functor>
auto buildGradientExpr(Functor functor) {
  if constexpr (w == 0) {
    return;
  } else if constexpr (w == 1) {
    return functor();
  } else {
    return functor(), buildGradientExpr<w - 1, Functor>(functor);
  }
}

但是使用此方法会导致Eigen发出运行时错误，因为该表达式只有一个初始化程序。

template <std::size_t w>
auto buildGradient() {
  Eigen::Matrix<float, 3, w> gradient;
  /* Emits an error about too few coefficients being passed to the initializer. */
  gradient << buildGradientExpr<w>([]() { /* Return 3x1 matrix */ });
  return gradient;
}

这是一个完整的可运行示例。

#include <Eigen/Dense>

#include <iostream>

#include <cstddef>

namespace {

template <std::size_t w, typename Functor>
auto buildGradientExpr(Functor functor) {
  if constexpr (w == 0) {
    return;
  } else if constexpr (w == 1) {
    return functor(w);
  } else {
    return functor(w), buildGradientExpr<w - 1, Functor>(functor);
  }
}

template <std::size_t w, typename Functor>
auto buildGradient(Functor functor) {
  Eigen::Matrix<float, 3, w> gradient;
  gradient << buildGradientExpr<w>(functor);
  return gradient;
}

} // namespace

int main() {

  constexpr std::size_t gradient_width = 10;

  auto gradient_functor = [](std::size_t w) {
    return Eigen::Matrix<float, 3, 1>::Constant(float(w) / gradient_width);
  };

  auto gradient = buildGradient<gradient_width>(gradient_functor);

  std::cout << gradient << std::endl;

  return 0;
}

是否有一种方法可以构造大小取决于模板参数的矩阵，而无需借助for循环？没有什么适合for循环的，这是我在此期间使用的。我只想知道是否有一种使用模板循环初始化矩阵表达式的方法。

Edit：*我更新了示例，因为假定梯度函子返回的是向量而不是标量。但是，出现示例问题。

Answer 1

如果可以使用c ++ 11，那么使用variadic templates怎么办？

例如：

#include <iostream>
#include <eigen3/Eigen/Dense>

template<class... Args>
Eigen::Matrix<double, 3, sizeof...(Args)/3> buildGradient(Args&&... args)
{
    static_assert(sizeof...(Args) % 3 == 0, "Wrong number of elements");
    return Eigen::Matrix<double, 3, sizeof...(Args)/3>{ args... };
}

int main()
{
    std::cout << buildGradient(1.0,2.0,3.0) << std::endl;
}

已扩展Eigen::Matrix类，如here所述，以包括std::initializer_list中的构造函数：

// File "Eigen_plugin.h"
template<class T>
Matrix(std::initializer_list<T> elems) : Base()
{
    Base::_check_template_params();

    int idx = 0;
    for(auto el : elems)
        coeffRef(idx++) = el;
}