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.
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/
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/
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>
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
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();
}
A matrix that has been compressed to eliminate zero-values is a sparse matrix, whereas a matrix with zero and nonzero values is a dense matrix.
Operations using standard dense-matrix structures and algorithms are slow and inefficient when applied to large sparse matrices as processing and memory are wasted on the zeros. Sparse data is by nature more easily compressed and thus requires significantly less storage.
Storage: When there is the maximum number of zero elements and the minimum number of non-zero elements then we use a sparse array over a simple array as it requires less memory to store the elements. In the sparse array, we only store the non-zero elements.
A sparse matrix is a matrix that is comprised of mostly zero values. Sparse matrices are distinct from matrices with mostly non-zero values, which are referred to as dense matrices. A matrix is sparse if many of its coefficients are zero.
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).
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.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With