2019-08-12
Function that accepts both Eigen Dense and Sparse Matrices
stackoverflow
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>
return A + A;
}

template <class Derived>
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>
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>
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.

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

@Max Langhof suggested using

``````template <class Mat>
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>
{
return ((A+A).eval()).transpose();
}

``````
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>
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>
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).