接受Eigen稠密矩阵和稀疏矩阵的函数

Question

我正在为一个开源数学库添加稀疏矩阵支持，并希望不必为Dense和Sparse矩阵类型分别编写重复的函数。

下面的示例显示了一个add函数。其中包括两个函数的工作示例，以及两次失败的尝试。代码示例的godbolt链接在下面提供。

我查看了Eigen文档中关于编写接受Eigen类型的函数的部分，但是它们使用Eigen::EigenBase的方法不可行，因为MatrixBase和SparseMatrixBase都有特定的方法可用，而这些方法在EigenBase中不存在。

https://eigen.tuxfamily.org/dox/TopicFunctionTakingEigenTypes.html

我们使用C++14，非常感谢您的帮助和时间!!

#include <Eigen/Core>
#include <Eigen/Sparse>
#include <iostream>

// Sparse matrix helper
using triplet_d = Eigen::Triplet<double>;
using sparse_mat_d = Eigen::SparseMatrix<double>;
std::vector<triplet_d> tripletList;

// Returns plain object
template <typename Derived>
using eigen_return_t = typename Derived::PlainObject;

// Below two are the generics that work
template <class Derived>
eigen_return_t<Derived> add(const Eigen::MatrixBase<Derived>& A) {
    return A + A;
}

template <class Derived>
eigen_return_t<Derived> add(const Eigen::SparseMatrixBase<Derived>& A) {
    return A + A;
}

int main()
{
  // Fill up the sparse and dense matrices
  tripletList.reserve(4);
  tripletList.push_back(triplet_d(0, 0, 1));
  tripletList.push_back(triplet_d(0, 1, 2));
  tripletList.push_back(triplet_d(1, 0, 3));
  tripletList.push_back(triplet_d(1, 1, 4));

  sparse_mat_d mat(2, 2);
  mat.setFromTriplets(tripletList.begin(), tripletList.end());

  Eigen::Matrix<double, -1, -1> v(2, 2);
  v << 1, 2, 3, 4;

  // Works fine
  sparse_mat_d output = add(mat * mat);
  std::cout << output;

  // Works fine
  Eigen::Matrix<double, -1, -1> output2 = add(v * v);
  std::cout << output2;

} 

我希望有一个函数可以同时处理稀疏矩阵和密集矩阵的加法，而不是两个不同的函数。但是下面的尝试并没有成功。

模板模板类型

这是我显然很失败的一次尝试，用模板模板类型替换上述的两个add函数会导致一个模糊的基类错误。

template <template <class> class Container, class Derived>
Container<Derived> add(const Container<Derived>& A) {
    return A + A;    
}

错误：

<source>: In function 'int main()':
<source>:35:38: error: no matching function for call to 'add(const Eigen::Product<Eigen::SparseMatrix<double, 0, int>, Eigen::SparseMatrix<double, 0, int>, 2>)'
   35 |   sparse_mat_d output = add(mat * mat);
      |                                      ^
<source>:20:20: note: candidate: 'template<template<class> class Container, class Derived> Container<Derived> add(const Container<Derived>&)'
   20 | Container<Derived> add(const Container<Derived>& A) {
      |                    ^~~
<source>:20:20: note:   template argument deduction/substitution failed:
<source>:35:38: note:   'const Container<Derived>' is an ambiguous base class of 'const Eigen::Product<Eigen::SparseMatrix<double, 0, int>, Eigen::SparseMatrix<double, 0, int>, 2>'
   35 |   sparse_mat_d output = add(mat * mat);
      |                                      ^
<source>:40:52: error: no matching function for call to 'add(const Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0>)'
   40 |   Eigen::Matrix<double, -1, -1> output2 = add(v * v);
      |                                                    ^
<source>:20:20: note: candidate: 'template<template<class> class Container, class Derived> Container<Derived> add(const Container<Derived>&)'
   20 | Container<Derived> add(const Container<Derived>& A) {
      |                    ^~~
<source>:20:20: note:   template argument deduction/substitution failed:
<source>:40:52: note:   'const Container<Derived>' is an ambiguous base class of 'const Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0>'
   40 |   Eigen::Matrix<double, -1, -1> output2 = add(v * v);
      |                                                    ^

我认为这是来自这里的相同的钻石继承问题：

https://www.fluentcpp.com/2017/05/19/crtp-helper/

使用 std::conditional_t

下面的示例尝试使用 conditional_t 推断正确的输入类型。

#include <Eigen/Core>
#include <Eigen/Sparse>
#include <iostream>

// Sparse matrix helper
using triplet_d = Eigen::Triplet<double>;
using sparse_mat_d = Eigen::SparseMatrix<double>;
std::vector<triplet_d> tripletList;


// Returns plain object
template <typename Derived>
using eigen_return_t = typename Derived::PlainObject;

// Check it Object inherits from DenseBase
template<typename Derived>
using is_dense_matrix_expression = std::is_base_of<Eigen::DenseBase<std::decay_t<Derived>>, std::decay_t<Derived>>;

// Check it Object inherits from EigenBase
template<typename Derived>
using is_eigen_expression = std::is_base_of<Eigen::EigenBase<std::decay_t<Derived>>, std::decay_t<Derived>>;

// Alias to deduce if input should be Dense or Sparse matrix
template <typename Derived>
using eigen_matrix = typename std::conditional_t<is_dense_matrix_expression<Derived>::value,
 typename Eigen::MatrixBase<Derived>, typename Eigen::SparseMatrixBase<Derived>>;

template <typename Derived>
eigen_return_t<Derived> add(const eigen_matrix<Derived>& A) {
    return A + A;
}

int main()
{
  tripletList.reserve(4);

  tripletList.push_back(triplet_d(0, 0, 1));
  tripletList.push_back(triplet_d(0, 1, 2));
  tripletList.push_back(triplet_d(1, 0, 3));
  tripletList.push_back(triplet_d(1, 1, 4));

  sparse_mat_d mat(2, 2);
  mat.setFromTriplets(tripletList.begin(), tripletList.end());
  sparse_mat_d output = add(mat * mat);

  std::cout << output;
  Eigen::Matrix<double, -1, -1> v(2, 2);
  v << 1, 2, 3, 4;
  Eigen::Matrix<double, -1, -1> output2 = add(v * v);
  std::cout << output2;

} 

这会抛出错误

<source>: In function 'int main()':
<source>:94:38: error: no matching function for call to 'add(const Eigen::Product<Eigen::SparseMatrix<double, 0, int>, Eigen::SparseMatrix<double, 0, int>, 2>)'
   94 |   sparse_mat_d output = add(mat * mat);
      |                                      ^
<source>:79:25: note: candidate: 'template<class Derived> eigen_return_t<Derived> add(eigen_matrix<Derived>&)'
   79 | eigen_return_t<Derived> add(const eigen_matrix<Derived>& A) {
      |                         ^~~
<source>:79:25: note:   template argument deduction/substitution failed:
<source>:94:38: note:   couldn't deduce template parameter 'Derived'
   94 |   sparse_mat_d output = add(mat * mat);
      |                                      ^
<source>:99:52: error: no matching function for call to 'add(const Eigen::Product<Eigen::Matrix<double, -1, -1>, Eigen::Matrix<double, -1, -1>, 0>)'
   99 |   Eigen::Matrix<double, -1, -1> output2 = add(v * v);
      |                                                    ^
<source>:79:25: note: candidate: 'template<class Derived> eigen_return_t<Derived> add(eigen_matrix<Derived>&)'
   79 | eigen_return_t<Derived> add(const eigen_matrix<Derived>& A) {
      |                         ^~~
<source>:79:25: note:   template argument deduction/substitution failed:
<source>:99:52: note:   couldn't deduce template parameter 'Derived'
   99 |   Eigen::Matrix<double, -1, -1> output2 = add(v * v);

这似乎是因为依赖类型的相关参数不能像链接中所述那样被推导出来。

https://deque.blog/2017/10/12/why-template-parameters-of-dependent-type-names-cannot-be-deduced-and-what-to-do-about-it/

Godbolt示例

下面的Godbolt示例包含了以上所有实例，可以进行操作

https://godbolt.org/z/yKEAsn

有没有仅使用一个函数而不是两个的方法？我们有很多既支持稀疏矩阵又支持密集矩阵的函数，因此避免代码重复将是很好的选择。

编辑：可能的答案

@Max Langhof建议使用

template <class Mat>
auto add(const Mat& A) {
 return A + A; 
}

auto关键字在Eigen中有一定的危险性。

https://eigen.tuxfamily.org/dox/TopicPitfalls.html

但是

template <class Mat> 
typename Mat::PlainObject add(const Mat& A) { 
    return A + A; 
}

虽然说实话我不完全确定为什么在这种情况下返回一个普通对象是可行的

编辑编辑

有几个人提到了使用auto关键字。但遗憾的是，Eigen与auto不兼容，可以参考下面链接中关于C++11和auto的第二部分。

https://eigen.tuxfamily.org/dox/TopicPitfalls.html

有些情况下可以使用auto，但我想看看是否存在一种通用的“auto”方式，适用于Eigen的模板返回类型。

如果您想要一个关于使用auto导致段错误的示例，您可以尝试将add替换为

template <typename T1>
auto add(const T1& A) 
{
    return ((A+A).eval()).transpose();
}

- Steve Bronder

---

你使用的是哪个版本的C++？ - Max Langhof

C++14，谢谢，我会在帖子中添加！ - Steve Bronder

你为什么不能直接这样做呢：template <class Mat> auto add(const Mat& A) { return A + A; }（可能需要一些SFINAE来限制它只能用于Eigen矩阵），有什么特别的原因吗？或者这样会遇到表达式模板的问题吗？ - Max Langhof

...哈哈，让我试试这个！唯一让我担心的是 auto 关键字，因为 Eigen 在他们的文档中说不要使用它。https://eigen.tuxfamily.org/dox/TopicPitfalls.html - Steve Bronder

啊！所以如果我将返回值设为Mat :: PlainObject，那就可以工作了！虽然我需要检查一下是否存在任何问题，因为说实话这对我来说并不完全合理。让我探索一下是否会导致任何奇怪的错误，如果没有，我肯定会接受它作为答案！ - Steve Bronder

Answer 1

如果您想通过 EigenBase<Derived>，则可以使用 .derived() 提取底层类型（实际上，这只是将其转换为 Derived const&）：

template <class Derived>
eigen_return_t<Derived> add(const Eigen::EigenBase<Derived>& A_) {
    Derived const& A = A_.derived();
    return A + A;
}

---

对于这个特定的例子，更高级的做法是，由于您两次使用了A，您可以使用内部评估器结构来表示：

template <class Derived>
eigen_return_t<Derived> add2(const Eigen::EigenBase<Derived>& A_) {
    // A is used twice:
    typedef typename Eigen::internal::nested_eval<Derived,2>::type NestedA;
    NestedA A (A_.derived());
    return A + A;
}

这样做的好处是当将一个产品作为A_传递时，在计算A+A时不会被计算两次，但如果A_是像Block<...>这样的东西，它将不会不必要地被复制。然而，使用internal功能并不是真正推荐的（其API随时可能更改）。

- chtz

---

有趣！我马上要登机了，但我一有时间就会试试看。 - Steve Bronder

这似乎可行！有一个问题，因为您似乎对Eigen很了解。eigen_return_t <Derived>是设置返回类型的好方法还是应该是其他东西？例如，在我的编辑案例中，使用return((A+A).eval()).transpose();？ - Steve Bronder

此外，为了让您看到自己的努力成果，这里提供了一个链接，指向我正在工作的稀疏矩阵设计文档 https://github.com/SteveBronder/design-docs/blob/spec/sparse-matrices/designs/0004-sparse-matrices.md - Steve Bronder

返回 Derived::PlainObject 总是一个安全的选择，即它将按值返回。在许多情况下，这可能会导致不必要的临时变量，但返回自定义表达式树很可能过度设计（如果您想走这条路，应该研究Eigen的内部结构，但遗憾的是文档不是很好）。但按值返回通常至少会给您返回值优化（特别是如果只有一条路径返回），即没有冗余副本。 - chtz

很好，我已将此答案添加到设计文档中。但希望有一种方法可以删除关于派生内容的样板文件。 - Steve Bronder

实际上，仔细研究后发现 A.derived 会导致 mov 和 call 的发生。我想研究一下如何避免这种复制的发生。https://godbolt.org/z/arPlXs - Steve Bronder

Answer 2

你的编译器出现了以下问题：

无法推断模板参数“Derived”

传递所需的类型给Derived可能会起作用，如下所示：

传递所需的类型给Derived可能会起作用，如下所示：

add<double>(v * v)

然而，我不确定，因为对我来说Eigen::Matrix与Eigen::MatrixBase不是相同的类型。

但是，如果你对类型限制较少，编译器将能够找出类型：

最初的回答：

template <typename T>
auto add(const T& A) {
    return A + A;
}

编辑： 在评论中看到已经有人发表了这个解决方案，并且Eigen文档建议不要使用auto。我对Eigen不熟悉，但从我粗略浏览文档的结果来看，可能是Eigen产生表示表达式的结果 - 例如，将矩阵加法表示为算法的对象；而不是矩阵加法本身的结果。在这种情况下，如果你知道A + A的结果是类型T（在我看来对于operator+应该如此），那么可以按以下方式编写它：

template <typename T>
T add(const T& A) {
    return A + A;
}

在矩阵示例中，这应该强制返回一个矩阵结果；而不是表示表达式的对象。但是，由于您最初使用的是eigen_result_t，我不能百分之百确定。

最初的回答: 在这个矩阵的例子中，这将会强制返回一个矩阵结果，而不是表示表达式的对象。然而，由于你一开始使用的是 eigen_result_t，我不完全确定。

- j00hi

---

谢谢你的回答！如果可能的话，我想避免在被调用者上给编译器模板类型。如果将auto替换为T:PlainObject，你上面的解决方案就可以工作了！一旦我调查确认T::PlainObject是可行的，这就是一个好答案！希望有办法让你和@max获得荣誉，因为他在评论中也给出了相同的答案。 - Steve Bronder

是的，当Eigen开始谈论auto及其返回类型时，我的大脑就会变成“dogWithChemistrySetMeme.jpg”。PlainObject返回类型有点合理，但我想弄清楚为什么他们不在他们的代码中使用它。 - Steve Bronder

没错，所以我关于临时变量的推理是无意义的。我会将其更改为Eigen文档中的“抽象表达式”情况。 - j00hi

“T::PlainObject” 对我来说有点不直观。仅使用 “T” 作为返回类型不也可以吗？ - j00hi

正是因为这个原因，我想避免它。虽然它能工作，但显然有些不直观/愚蠢。 - Steve Bronder

@j00hi T::PlainObject 是必须的，例如当 T 为 Product<X,Y> 或 Block<X> 对象（或任何不是普通对象的内容）时。 - chtz

Answer 3

我没有理解你所有的代码和注释。不管怎样，似乎你的问题归结为找到一种编写可以接受多个矩阵类型的函数的方法。

template <typename T>
auto add(const T& A)
{
    return 2*A;
}

您也可以添加两个不同类型的矩阵：

template <typename T1, typename T2>
auto add(const T1& A, const T2& B) -> decltype(A+B) // decltype can be omitted since c++14
{
    return A + B;
}

然后，add(A,A) 的结果与 add(A) 相同。但我认为带有 2 个参数的 add 函数更有意义。而且它更加灵活，因为您可以将稀疏矩阵与密集矩阵相加。

int main()
{
    constexpr size_t size = 10;
    Eigen::SparseMatrix<double> spm_heap(size,size);
    Eigen::MatrixXd m_heap(size,size);
    Eigen::Matrix<double,size,size> m_stack; 

    // fill the matrices

    std::cout << add(spm_heap,m_heap);
    std::cout << add(spm_heap,m_stack);

    return 0;
}

编辑

关于您提到的在使用Eigen时不应该使用auto的编辑，这非常有趣！

template <typename T>
auto add(const T& A) 
{
    return ((A+A).eval()).transpose();
}

这会产生一个“段错误”。为什么？“auto”可以很好地推断出类型，但推断出的类型不是“decltype(A)”，而是该类型的“引用”。为什么？我最初认为这是因为返回值周围有括号（如果感兴趣，请阅读此处），但似乎是由于“transpose”函数的返回类型。

无论如何，很容易克服这个问题。正如你建议的那样，你可以删除“auto”：

template <typename T>
T add(const T& A) 
{
    return ((A+A).eval()).transpose();
}

或者，您可以使用 auto 但要指定所需的返回类型：

template <typename T>
auto add(const T& A) -> typename std::remove_reference<decltype(A)>::type // or simply decltype(A.eval())
{
    return ((A+A).eval()).transpose();
}

现在，对于这个特定的add函数，第一个选项（省略auto）是最好的解决方案。然而，对于另一个接受不同类型的2个参数的add函数，这是相当不错的解决方案：

template <typename T1, typename T2>
auto add(const T1& A, const T2& B) -> decltype((A+B).eval())
{
    return ((A+B).eval()).transpose();
}

- mfnx

---

谢谢你的回答！我编辑了上面的帖子，链接到了Eigen文档，在那里他们说不要使用auto。你的答案对于这里的简单加法函数是有效的，但我认为我需要更新问题，以展示一个auto失败的例子。 - Steve Bronder

我添加了一篇文章，展示了auto导致段错误的一个例子。 - Steve Bronder