11

在 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 版本创建了矩阵的副本（LDLT） - 所以它更适合缓存。

有没有办法在 Eigen 中进行就地计算？或者也许还有其他方法可以改进它？

编辑：

遵循 Fezvez 的建议并使用 LLT 而不是 LDLT 我得到：

Eigen 0.16 s

488184

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

sizeof(A.llt()) == 656

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

编辑：

我通过从 LDLT 子类化（它的内部矩阵受到保护）删除了分配，并直接填充它。现在 LDLT 的结果是：

Eigen 0.26 s
488209

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

LLT 的子类化也有效，但没有提供如此好的效果。

4

2 回答 2

8

你的基准是不公平的，因为 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
于 2012-11-15T07:46:19.307 回答
2

我按照您关于就地计算的提示进行操作：

使用完全相同的代码，并使用函数A.llt().solveInPlace(b);和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 实现有很多开销......）

于 2012-11-14T13:13:41.760 回答