本征稀疏矩阵的性能调整
performance tuning on Eigen sparse matrix
我已经使用 Eigen
的 SparseMatrix
实现了一些东西,基本上是这样的,
SparseMatrix W;
...
W.row(i) += X.row(j); // X is another SparseMatrix, both W and X are row major.
...
我通过google-pprof
对代码做了一些perf-profiling,我认为上面的代码有问题,见下图,
图 1
然后 图 2
最后 图 3
看起来 operator+=
带来了很多 memory-copy 东西。
我不太了解SparseMatrix
操作的内部原理,但是有什么推荐的方法可以优化上面的代码吗?
如果 X 的稀疏度是 W 的稀疏度的子集,那么您可以编写自己的函数来就地进行加法:
namespace Eigen {
template<typename Dst, typename Src>
void inplace_sparse_add(Dst &dst, const Src &src)
{
EIGEN_STATIC_ASSERT( ((internal::evaluator<Dst>::Flags&RowMajorBit) == (internal::evaluator<Src>::Flags&RowMajorBit)),
THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH);
using internal::evaluator;
evaluator<Dst> dst_eval(dst);
evaluator<Src> src_eval(src);
assert(dst.rows()==src.rows() && dst.cols()==src.cols());
for (Index j=0; j<src.outerSize(); ++j)
{
typename evaluator<Dst>::InnerIterator dst_it(dst_eval, j);
typename evaluator<Src>::InnerIterator src_it(src_eval, j);
while(src_it)
{
while(dst_it && dst_it.index()!=src_it.index())
++dst_it;
assert(dst_it);
dst_it.valueRef() += src_it.value();
++src_it;
}
}
}
}
这是一个用法示例:
int main()
{
int n = 10;
MatrixXd R = MatrixXd::Random(n,n);
SparseMatrix<double, RowMajor> A = R.sparseView(0.25,1), B = 0.5*R.sparseView(0.65,1);
cout << A.toDense() << "\n\n" << B.toDense() << "\n\n";
inplace_sparse_add(A, B);
cout << A.toDense() << "\n\n";
auto Ai = A.row(2);
inplace_sparse_add(Ai, B.row(2));
cout << A.toDense() << "\n\n";
}
我已经使用 Eigen
的 SparseMatrix
实现了一些东西,基本上是这样的,
SparseMatrix W;
...
W.row(i) += X.row(j); // X is another SparseMatrix, both W and X are row major.
...
我通过google-pprof
对代码做了一些perf-profiling,我认为上面的代码有问题,见下图,
图 1
然后 图 2
最后 图 3
看起来 operator+=
带来了很多 memory-copy 东西。
我不太了解SparseMatrix
操作的内部原理,但是有什么推荐的方法可以优化上面的代码吗?
如果 X 的稀疏度是 W 的稀疏度的子集,那么您可以编写自己的函数来就地进行加法:
namespace Eigen {
template<typename Dst, typename Src>
void inplace_sparse_add(Dst &dst, const Src &src)
{
EIGEN_STATIC_ASSERT( ((internal::evaluator<Dst>::Flags&RowMajorBit) == (internal::evaluator<Src>::Flags&RowMajorBit)),
THE_STORAGE_ORDER_OF_BOTH_SIDES_MUST_MATCH);
using internal::evaluator;
evaluator<Dst> dst_eval(dst);
evaluator<Src> src_eval(src);
assert(dst.rows()==src.rows() && dst.cols()==src.cols());
for (Index j=0; j<src.outerSize(); ++j)
{
typename evaluator<Dst>::InnerIterator dst_it(dst_eval, j);
typename evaluator<Src>::InnerIterator src_it(src_eval, j);
while(src_it)
{
while(dst_it && dst_it.index()!=src_it.index())
++dst_it;
assert(dst_it);
dst_it.valueRef() += src_it.value();
++src_it;
}
}
}
}
这是一个用法示例:
int main()
{
int n = 10;
MatrixXd R = MatrixXd::Random(n,n);
SparseMatrix<double, RowMajor> A = R.sparseView(0.25,1), B = 0.5*R.sparseView(0.65,1);
cout << A.toDense() << "\n\n" << B.toDense() << "\n\n";
inplace_sparse_add(A, B);
cout << A.toDense() << "\n\n";
auto Ai = A.row(2);
inplace_sparse_add(Ai, B.row(2));
cout << A.toDense() << "\n\n";
}