为什么在以下示例中，Eigen比ublas慢5倍？

Question

我使用的Eigen版本使用“true”固定大小的矩阵和向量，采用更好的算法（在uBlas中采用LDLT而不是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进行子类化也可以工作，但效果没有那么好。

- qble

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

- ggael

---

版本可以轻松调整以实现此目的，但是有更好的选择 - solveInPlace。但这完全是关于方程中的向量，而不是矩阵 - 所以根本没有影响。 "还要确保启用了优化编译" - 是的，我检查过：1）我正在运行x64 - 默认情况下启用sse，2）我甚至尝试显式启用sse - 结果相同。 "使用g ++编译" - 问题是关于MSVC 2010的。 - qble

我也试图避免矩阵复制，但结果没有任何收益。 - 你是通过子类化来做到的吗？“另外，增加n会增加速度差异（这里n = 90）：”- 在这个问题的上下文中，这样增加是毫无意义的。当您将N增加到相当高的数字时，解决方案本身的复杂度O（N ^ 3）开始占主导地位。当N = 9时，我正在击中Eigen情况下的某些缓存边界，并且这会影响性能.. - qble

使用MSVC编译器编译需要2.66秒，而使用GCC编译器只需要0.08秒，这意味着问题出在编译器上，而不是Eigen代码本身。通常这种开销是由于内联不良引起的。当您对LDLT进行子类化时获得的巨大加速表明，您的修改有助于MSVC生成更好的代码。如果您可以将初始代码生成的ASM发送给我，我们可以确定确切的问题并进行上游修复。您的修改也可能很有趣。 - ggael

使用MSVC编译器的时间为2.66秒，使用GCC编译器的时间为0.08秒，对于相同的代码来说，这意味着问题出在编译器上，而不是Eigen代码。一些编译器可以比其他编译器更积极地进行优化。然而，在那种特定情况下，库具有避免不必要复制的所有工具，而不依赖于强大的优化器。所需的只是因式分解操作，它允许修改原始矩阵-就是这样，Boost.uBlas正是这样做的，我的修改也完全相同。而Eigen会复制矩阵：“m_matrix = a;”-这就是问题所在，无需检查汇编即可识别。 - qble

很抱歉要坚持，但是这个副本不是性能下降的源头：1）我可以看到GCC没有优化掉这个副本，并且使用gcc没有减速；2）LLT求解器也执行了一次复制，但在LLT情况下您没有观察到任何减速；3）您认为额外的O（n ^ 2）操作可能会解释一个O（n ^ 3）操作的10倍减速，这是无意义的；4）在堆栈上分配的9x9矩阵太小，不能归咎于缓存未命中。 - ggael

1. 一些其他的东西可以被优化，比如更好的堆栈布局等等。
2. LLT使用的数据较少 - 只是矩阵的副本，但LDLT使用了额外的数据，它的大小更大，请查看源代码。当我通过添加额外的数据来修改LLT结构时，它也会变慢。
3. 我并不是在争论复制操作本身是一个问题，我是说内存布局是一个问题，而不是复制本身（当我增加LLT结构的大小时 - 没有任何复制操作，只是原始数据大小的改变）。
4. 我对我的LDLT进行的更改使其更快，主要是关于原地计算，即删除一个分配。

- qble

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

- B. Decoster

---

solveInPlace 只使用向量 inplace，而不是矩阵 inplace。矩阵仍然被分配。 - qble

是的，但我想告诉你的是，我按照你的提示尝试解决就地问题。当然，我很失望地看到它只是替换了b而不是进行一些矩阵就地操作，但我仍然尝试了一下。然后我发现另一个结果，即LLT似乎比LDLT在这个特定问题上更好 =) - B. Decoster

我非常难过看到Eigen的LDLT实现有很多开销 - 我也很难过，我向许多人推荐了Eigen。 - qble

我绕过了堆栈分配，现在正在进行原地LDLT - 现在Eigen比uBlas快2倍，请查看底部问题的编辑。 - qble