c++ - 为什么 Eigen 在下面的例子中比 ublas 慢 5 倍？

Question

在 Eigen 版本中，我使用“真正的”固定大小矩阵和 vector ，更好的算法(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 版本创建矩阵的拷贝 (LDLT) - 因此它更适合缓存。

有什么方法可以在 Eigen 中进行就地计算吗？或者也许还有其他改进方法？

---

编辑:

按照 Fezvez 的建议并使用 LLT 而不是 LDLT，我得到:

Eigen 0.16 s

488184

很好，但它仍然进行不必要的矩阵堆栈分配:

sizeof(A.llt()) == 656

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

---

编辑:

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

Eigen 0.26 s
488209

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

LLT 的子类化也可以，但效果不佳。

Answer 1

您的基准测试不公平，因为 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