Function that accepts both Eigen Dense and Sparse Matrices

Question

I'm working on adding Sparse matrix support to an open source math library and would like to not have duplicated functions for both Dense and Sparse matrix types.

The below example shows an add function. A working example with two functions, then two attempts that failed. A godbolt link to the code examples are available below.

I've looked over the Eigen docs on writing functions that take Eigen types but their answers of using Eigen::EigenBase does not work because both MatrixBase and SparseMatrixBase have particular methods available that do not exist in EigenBase

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

We use C++14, any help and your time is very appreciated!!

#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;

} 

Instead of the two add functions I would just like to have one that takes in both sparse and dense matrices, but the attempts below have not worked out.

Template Template type

An obviously poor attempt on my part, but replacing the two add functions above with a template template type causes an ambiguous base class error.

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

Error:

<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);
      |                                                    ^

I believe It's the same diamond inheritance problem from here:

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

Using std::conditional_t

The below attempts to use conditional_t to deduce the correct input type

#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;

} 

This throws the error

<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);

This seems to be because dependent parameters of dependent types can't be deduced like this link goes over.

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

Godbolt Example

The godbolt below has all of the instances above to play with

https://godbolt.org/z/yKEAsn

Is there some way to only have one function instead of two? We have a lot of functions that can support both sparse and dense matrices so it would be nice to avoid the code duplication.

Edit: Possible Answer

@Max Langhof suggested using

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

The auto keyword is a bit dangerous with Eigen

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

But

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

works, though tbh I'm not entirely sure why returning a plain object works in this scenario

Edit Edit

Several people have mentioned the use of the auto keyword. Sadly Eigen does not play well with auto as referenced in the second on C++11 and auto in the link below

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

It's possible to use auto for some cases, though I'd like to see if there is a generic auto'ish way that is complaint for Eigen's template return types

For an example of a segfault with auto you can try replace add with

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

Satisfy

What C++ version are you using? – Poirer

C++14, thanks I'll add that in the post! – Satisfy

Any particular reason you can't just do template <class Mat> auto add(const Mat& A) { return A + A; } (with possibly some SFINAE to constrain it to Eigen matrices)? Or does this somehow run into issues with expression templates? – Poirer

...LoL. Let me try this out! The only thing I worry about is the auto keyword since Eigen says not to use this in their docs. eigen.tuxfamily.org/dox/TopicPitfalls.html – Satisfy

Ah! So if I have the return be Mat::PlainObject that does work! Though I need to check if there are any issues with that since tbh that doesn't make complete sense to me. Let me explore if that causes any weird goofs and if not I would def accept that as an answer! – Satisfy

Answer 1

If you want to pass EigenBase<Derived>, you can extract the underlying type using .derived() (essentially, this just casts to Derived const&):

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

---

More advanced, for this particular example, since you are using A twice, you can express that using the internal evaluator structure:

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;
}

This has the advantage that when passing a product as A_ it won't get evaluated twice when evaluating A+A, but if A_ is something like a Block<...> it will not get copied unnecessarily. However, using internal functionality is not really recommended (the API of that could change at any time).

Expend

Interesting! I'm about to get on an airplane but I'll try this out when I get a minute – Satisfy

This does seem to work! One Q since it looks like you know a lot about Eigen. Is the eigen_return_t<Derived> a good way to set the return type or should it be something else? For instance in my edited case above with return((A+A).eval()).transpose();? – Satisfy

Also so you can see the fruits of your labor here is a link to the sparse matrix design doc I'm working on github.com/SteveBronder/design-docs/blob/spec/sparse-matrices/… – Satisfy

Returning Derived::PlainObject is always a safe option, i.e., it will return by value. In many cases this can cause unnecessary temporaries, but returning custom expression trees is likely overkill (if you want to go down that road, you should study the internals of Eigen, it is unfortunately not very well documented). But returning by value will usually at least give you return-value-optimization (especially, if you have only one path which returns), i.e., no redundant copy. – Expend

🙏🙏 Excellent I've added this answer to the design doc. Wish there was a way to remove the bit of boilerplate about the derived stuff – Satisfy

Actually looking into this it seems like the A.derived causes a mov and call to happen. I'd like to look into how not to have that copy happen. godbolt.org/z/arPlXs – Satisfy

Answer 2

The problem of your compiler is the following:

couldn't deduce template parameter 'Derived'

Passing the required type for Derived should probably work, like follows:

add<double>(v * v)

However I'm not sure because Eigen::Matrix is not the same type as Eigen::MatrixBase as it appears to me.

However, if you restrict the compiler less on the type, it will be able to figure out the type:

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

Edit:

Just saw in the comments that this solution has already been posted and that the Eigen documentation recommends to not use auto. I am not familiar with Eigen, but as it appears to me from skimming over the documentation, it could be that Eigen produces results which represent expressions - e.g. an object representing the matrix addition as an algorithm; not the matrix addition result itself. In this case, if you know that A + A results in type T (which it actually should for operator+ in my opinion) you could write it like follows:

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

In the matrix example, this should force a matrix result to be returned; not the object representing the expression. However, since you have been originally using eigen_result_t, I'm not 100% sure.

Synthesize

Thanks for the answer! I'd like to avoid having to give the compiler the template type on the callee if possible. Your above solution works if auto is replaced with T:PlainObject! Once I investigate if the T::PlainObject is fine then this is a good answer! Wish there was some way to give you and @max the credit as he answered with the same in the comments. – Satisfy

Yeah when Eigen starts talking about auto and its return types my brain turns into dogWithChemistrySetMeme.jpg PlainObject return sort of makes sense but i want to figure out why they dont use that in their code – Satisfy

That's true, so my reasoning about temporaries is pointless. I'll change that to the "abstract expression" case from the Eigen documentation. – Synthesize

That T::PlainObject seems a bit unintuitive to me. Doesn't simply T also work for the return type? – Synthesize

^Exactly why I want to avoid it. It works but something seems obviously unintuitive / silly – Satisfy

@Synthesize T::PlainObject is necessary, e.g., if T is a Product<X,Y> or a Block<X> object (or anything other than a plain object already). – Expend

Answer 3

I haven't understood all of your code and comments. Anyways, it seems that your problem reduces to finding a way to write a function which can accept serveral matrix types.

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

You can also add 2 matrices of different types:

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;
}

Then, add(A,A) gives the same result as add(A). But the add function with 2 arguments makes more sense I think. And it's more versatile as you can sum a sparse matrix with a dense matrix.

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;
}

EDIT

About the edit where you state that auto should not be used with Eigen. This is quite interesting!

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

This produces a segfault. Why? auto does deduce the type well, but the deduced type is not decltype(A), but a reference of that type. Why? I first thought it was because of the parentheses around the return value (read here if interested), but it seems to be due to the return type of the transpose function.

Anyways, it is easy to overcome that problem. As you suggested, you could remove auto:

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

Or, you could use auto but specifying the desired return type:

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();
}

Now, for this particular add function, the first option (omitting auto) is the best solution. However, for the other add function which takes 2 arguments of different types, this is quite a good solution:

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

Coquina

Thanks for the answer! I edited the post above to link to the Eigen docs where they say not to use auto. Your answer works for the simple add function here, I think I need to update the Q to show an example where auto fails – Satisfy

I added a post to show an example of when auto causes a seg fault – Satisfy