问题 在下面的例子中,为什么Eigen比ublas慢5倍?


在Eigen版本,我使用“真正的”固定大小的矩阵和向量,更好的算法(LDLT与uBlas的LU),它在内部使用SIMD指令。那么,为什么它在下面的例子中比uBlas慢?

我确信,我做错了什么 - Eigen 必须 更快,或至少可比较。

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/progress.hpp>
#include <Eigen/Dense>
#include <iostream>

using namespace boost;
using namespace std;
const int n=9;
const int total=100000;

void test_ublas()
{
    using namespace boost::numeric::ublas;
    cout << "Boost.ublas ";
    double r=1.0;
    {
        boost::progress_timer t;
        for(int j=0;j!=total;++j)
        {
            //symmetric_matrix< double,lower,row_major,bounded_array<double,(1+n)*n/2> > A(n,n);
            matrix<double,row_major,bounded_array<double,n*n> > A(n,n);
            permutation_matrix< unsigned char,bounded_array<unsigned char,n> > P(n);
            bounded_vector<double,n> v;
            for(int i=0;i!=n;++i)
                for(int k=0;k!=n;++k)
                    A(i,k)=0.0;
            for(int i=0;i!=n;++i)
            {
                A(i,i)=1.0+i;
                v[i]=i;
            }
            lu_factorize(A,P);
            lu_substitute(A,P,v);
            r+=inner_prod(v,v);
        }
    }
    cout << r << endl;
}

void test_eigen()
{
    using namespace Eigen;
    cout << "Eigen ";
    double r=1.0;
    {
        boost::progress_timer t;
        for(int j=0;j!=total;++j)
        {
            Matrix<double,n,n> A;
            Matrix<double,n,1> b;
            for(int i=0;i!=n;++i)
            {
                A(i,i)=1.0+i;
                b[i]=i;
            }
            Matrix<double,n,1> x=A.ldlt().solve(b);
            r+=x.dot(x);
        }
    }
    cout << r << endl;
}

int main()
{
    test_ublas();
    test_eigen();
}

输出是:

Boost.ublas 0.50 s

488184
Eigen 2.66 s

488184

(Visual Studio 2010 x64发布)


编辑

对于

const int n=4;
const int total=1000000;

输出是:

Boost.ublas 0.67 s

1.25695e+006
Eigen 0.40 s

5.4e+007

我想,这种行为是由于uBlas版本计算就地分解,而Eigen版本创建了矩阵的COPY(LDLT) - 因此它适合缓存更糟糕。

有没有办法在Eigen中进行实地计算?或者还有其他方法可以改善它?


编辑:

根据Fezvez的建议并使用LLT代替LDLT,我得到:

Eigen 0.16 s

488184

这很好,但它仍然做不必要的矩阵堆栈分配:

sizeof(A.llt()) == 656

我宁愿避免它 - 它应该更快。


编辑:

我已经通过LDLT的子类化(它的内部矩阵受到保护)删除了分配,并直接填充它。现在LDLT的结果是:

Eigen 0.26 s
488209

它有效,但它是解决方法 - 不是一个真正的解决方案......

LLT的子类也有效,但没有提供如此大的效果。


2804
2017-11-14 10:39


起源



答案:


你的基准测试是不公平的,因为ublas版本在内部解决,而Eigen版本可以很容易地调整到这样做:

b=A.ldlt().solve(b);
r+=b.dot(b);

使用g ++ - 4.6 -O2 -DNDEBUG进行编译,得到(在2.3GHz CPU上):

Boost.ublas 0.15 s
488184

Eigen 0.08 s
488184

还要确保在启用优化的情况下编译,如果运行32位系统或(32位编译链),则启用SSE2。

编辑:我也试图避免矩阵复制,但这导致零增益。

另外,增加n,增加速度差(这里n = 90):

Boost.ublas 0.47 s
Eigen 0.13 s

8
2017-11-15 07:46



“可以很容易地调整版本” - 有更好的替代方案 - solveInPlace。但它完全是关于方程中的向量,而不是关于矩阵 - 因此根本没有影响。 “还要确保你在编译时启用了优化” - 是的,我检查了它:1)我正在运行x64 - sse默认启用,2)我试图甚至显式启用sse - 相同的结果。 “用g ++编译” - 问题是关于MSVC 2010。 - qble
“我也试图避免矩阵复制,但这会导致零增益。” - 你究竟是怎么做到的?通过子类化? “此外,增加n,增加速度差异(此处n = 90):” - 在这个问题的背景下,这种增加是毫无意义的。当你将N增加到如此高的数字时,解决方案本身O(N ^ 3)的复杂性开始占主导地位。在N = 9时 - 我在Eigen案例中遇到了一些缓存边界,并且达到了性能...... - qble
对于相同的代码,具有MSVC的2.66s和具有GCC的0.08s意味着问题在于编译器,而不是本征代码。一般来说,这种开销是由于内联不良引起的。在子类化LDLT时获得的巨大加速表明您的修改有助于MSVC生成更好的代码。如果您可以向我发送为初始代码生成的ASM,我们可以确定准确的问题并将其修复到上游。你的修改也可能很有趣 - ggael
“对于相同的代码,MSVC为2.66s,GCC为0.08s,意味着问题在于编译器,而不是本征代码” - 一些编译器可以比其他编译器更积极地进行优化。但是,在该特定情况下,库具有避免不必要的副本的所有工具,而不依赖于强大的优化器。所需要的只是分解操作,允许修改原始矩阵 - 就是这样,Boost.uBlas正是如此,我的修改也做了完全相同的事情。虽然Eigen复制了矩阵:“m_matrix = a;” - 这就是问题,如果没有ASM检查,它是可识别的。 - qble
很抱歉坚持,但没有这个副本不是性能下降的来源:1)我可以看到这个副本没有被GCC优化掉,并且gcc没有减速; 2)LLT求解器也执行复制,但您没有观察到LLT情况下的任何减速; 3)你认为额外的O(n ^ 2)op可以解释O(n ^ 3)操作的10倍减速,这是无意义的; 4)在堆栈上分配的9x9矩阵太小,无法判断缓存未命中。 - ggael


答案:


你的基准测试是不公平的,因为ublas版本在内部解决,而Eigen版本可以很容易地调整到这样做:

b=A.ldlt().solve(b);
r+=b.dot(b);

使用g ++ - 4.6 -O2 -DNDEBUG进行编译,得到(在2.3GHz CPU上):

Boost.ublas 0.15 s
488184

Eigen 0.08 s
488184

还要确保在启用优化的情况下编译,如果运行32位系统或(32位编译链),则启用SSE2。

编辑:我也试图避免矩阵复制,但这导致零增益。

另外,增加n,增加速度差(这里n = 90):

Boost.ublas 0.47 s
Eigen 0.13 s

8
2017-11-15 07:46



“可以很容易地调整版本” - 有更好的替代方案 - solveInPlace。但它完全是关于方程中的向量,而不是关于矩阵 - 因此根本没有影响。 “还要确保你在编译时启用了优化” - 是的,我检查了它:1)我正在运行x64 - sse默认启用,2)我试图甚至显式启用sse - 相同的结果。 “用g ++编译” - 问题是关于MSVC 2010。 - qble
“我也试图避免矩阵复制,但这会导致零增益。” - 你究竟是怎么做到的?通过子类化? “此外,增加n,增加速度差异(此处n = 90):” - 在这个问题的背景下,这种增加是毫无意义的。当你将N增加到如此高的数字时,解决方案本身O(N ^ 3)的复杂性开始占主导地位。在N = 9时 - 我在Eigen案例中遇到了一些缓存边界,并且达到了性能...... - qble
对于相同的代码,具有MSVC的2.66s和具有GCC的0.08s意味着问题在于编译器,而不是本征代码。一般来说,这种开销是由于内联不良引起的。在子类化LDLT时获得的巨大加速表明您的修改有助于MSVC生成更好的代码。如果您可以向我发送为初始代码生成的ASM,我们可以确定准确的问题并将其修复到上游。你的修改也可能很有趣 - ggael
“对于相同的代码,MSVC为2.66s,GCC为0.08s,意味着问题在于编译器,而不是本征代码” - 一些编译器可以比其他编译器更积极地进行优化。但是,在该特定情况下,库具有避免不必要的副本的所有工具,而不依赖于强大的优化器。所需要的只是分解操作,允许修改原始矩阵 - 就是这样,Boost.uBlas正是如此,我的修改也做了完全相同的事情。虽然Eigen复制了矩阵:“m_matrix = a;” - 这就是问题,如果没有ASM检查,它是可识别的。 - qble
很抱歉坚持,但没有这个副本不是性能下降的来源:1)我可以看到这个副本没有被GCC优化掉,并且gcc没有减速; 2)LLT求解器也执行复制,但您没有观察到LLT情况下的任何减速; 3)你认为额外的O(n ^ 2)op可以解释O(n ^ 3)操作的10倍减速,这是无意义的; 4)在堆栈上分配的9x9矩阵太小,无法判断缓存未命中。 - ggael


我按照你的提示进行了计算:

使用完全相同的代码,并使用功能 A.llt().solveInPlace(b); 和 A.ldlt().solveInPlace(b); (用x代替b),我明白了

There were  100000  Eigen standard LDLT linear solvers applied in  12.658  seconds 
There were  100000  Eigen standard LLT linear solvers applied in  4.652  seconds 

There were  100000  Eigen in place LDLT linear solvers applied in  12.7  seconds 
There were  100000  Eigen in place LLT linear solvers applied in  4.396  seconds 

对于这类问题,LLT求解器可能比LDLT求解器更合适吗? (我可以看到你正在处理大小为9的矩阵)

(顺便说一句,我回答了你之前关于9号线性求解器的问题,我很遗憾地发现Eigen的LDLT实现有很多开销......)


2
2017-11-14 13:13



solveInPlace只使用vector inplace,而不是matrix inplace。矩阵仍然分配。 - qble
是的,但我想告诉你的是,我遵循了你的解决方案。当然,我很失望地读到它只是替换了b而不是做一些矩阵内置的东西,但我仍然尝试过。我偶然发现了另一个结果,那就是LLT似乎比LDLT对这个特殊问题更好=) - Fezvez
“看到Eigen对LDLT的实施有很多开销,我感到非常难过” - 我也很遗憾地看到,我向许多人推荐了Eigen。 - qble
我解决堆栈分配,现在做就地LDLT - 现在比uBlas快2倍,检查问题底部的编辑。 - qble