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

Question

在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的子类也有效,但没有提供如此大的效果.

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

Answer 2

我按照你关于就地计算的提示：

使用完全相同的代码，并使用函数A.llt().solveInPlace(b);and A.ldlt().solveInPlace(b);（将 b 替换为 x），我得到了

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 实现有很多开销......）