为什么 Boost 的矩阵乘法比我的慢？

Question

我已经使用boost::numeric::ublas::matrix实现了一个矩阵乘法（请参见我的完整、可工作的boost代码）

Result result = read ();

boost::numeric::ublas::matrix<int> C;
C = boost::numeric::ublas::prod(result.A, result.B);

还有一个使用标准算法的实现（参见完整的标准代码）：

vector< vector<int> > ijkalgorithm(vector< vector<int> > A, 
                                    vector< vector<int> > B) {
    int n = A.size();

    // initialise C with 0s
    vector<int> tmp(n, 0);
    vector< vector<int> > C(n, tmp);

    for (int i = 0; i < n; i++) {
        for (int k = 0; k < n; k++) {
            for (int j = 0; j < n; j++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
    return C;
}

这是我测试速度的方法：

time boostImplementation.out > boostResult.txt
diff boostResult.txt correctResult.txt

time simpleImplementation.out > simpleResult.txt
diff simpleResult.txt correctResult.txt

这两个程序读取一个硬编码的文本文件，其中包含两个 2000 x 2000 的矩阵。 这两个程序都使用了相同的编译标志：

g++ -std=c++98 -Wall -O3 -g $(PROBLEM).cpp -o $(PROBLEM).out -pedantic

我的实现需要15秒，而使用boost库则需要超过4分钟！

编辑：在编译后进行了修改。

g++ -std=c++98 -Wall -pedantic -O3 -D NDEBUG -DBOOST_UBLAS_NDEBUG library-boost.cpp -o library-boost.out

我用ikj算法得到了28.19秒的结果，而使用Boost则需要60.99秒。因此，Boost仍然明显比我的实现慢。

为什么Boost比我的实现慢这么多？

- Martin Thoma

---

只有当你能制造出更好的轮子时，重新发明轮子才是一个好主意... - Mysticial

Boost.uBLAS旨在成为一个标准的接口，而不是一个强大的实现，因此除非您使用例如LAPACK后端，否则不要期望它运行速度很快。 - ildjarn

Boost uBLAS有一些可选的调试检查，会使程序变慢。请参考此FAQ http://www.boost.org/doc/libs/1_49_0/libs/numeric/ublas/doc/index.htm，并检查预处理宏BOOST_UBLAS_NDEBUG和NDEBUG。 - TJD

虽然阅读几个2k-by-2k矩阵不应该需要4分钟。 - Mysticial

大多数重新发明轮子的人通常相信这样做更好，无论它是否确实如此，否则他们可能不会去这么做。:-D - stinky472

Answer 1

如TJD所指出的那样，uBLAS版本的速度较慢可能部分原因是其调试功能。

以下是开启调试功能后 uBLAS 版本所需的时间：

real    0m19.966s
user    0m19.809s
sys     0m0.112s

以下是关闭调试模式后 uBLAS 版本所用的时间 (-DNDEBUG -DBOOST_UBLAS_NDEBUG编译器标志已添加):

real    0m7.061s
user    0m6.936s
sys     0m0.096s

关闭调试后，uBLAS版本快了近3倍。

剩余的性能差异可以通过引用uBLAS FAQ中的以下部分来解释："为什么uBLAS比（atlas-）BLAS慢得多"

ublas的一个重要设计目标是尽可能通用。

这种通用性往往是有代价的。特别地，prod函数模板可以处理不同类型的矩阵，如稀疏或三角形矩阵。幸运的是，uBLAS提供了针对密集矩阵乘法进行优化的替代方案，特别是 axpy_prod和block_prod。以下是比较不同方法的结果：

ijkalgorithm   prod   axpy_prod  block_prod
   1.335       7.061    1.330       1.278

正如您所看到的，axpy_prod和block_prod都比您的实现稍微快一些。仅测量计算时间，去除不必要的复制，并注意选择block_prod的块大小（我使用了64），可以使差异更加明显。

另请参阅uBLAS FAQ和Effective uBlas和一般代码优化。

- vitaut

---

你能否使用OP提供的版本运行相同的测试？ - mfontanini

@mfontanini：好的，我更新了答案。请注意，我使用了较小的（1000x1000）随机矩阵，因此所有时间都更短。 - vitaut

Answer 2

我相信你的编译器还没有进行足够的优化。uBLAS代码大量使用模板，而模板需要进行大量的优化。我在MS VC 7.1编译器中以1000x1000矩阵为例运行了你的代码（在发布模式下），结果如下：

- uBLAS花费10.064秒 - vector花费7.851秒

虽然仍有差异，但并不是很明显。uBLAS的核心概念是惰性求值，因此只有在需要时prod(A, B)才会计算结果，例如prod(A, B)(10,100)将立即执行，因为只有该元素需要计算。 因此，实际上不存在专门针对整个矩阵乘法进行优化的算法（请参见下文）。但是你可以在声明时帮助库一下，例如：

matrix<int, column_major> B;

这将会减少执行时间至4.426秒，轻松胜过您的函数。该声明使得在矩阵相乘时访问内存更加连续，优化了缓存的使用。

P.S. 已经仔细阅读uBLAS文档到最后 ;)，你应该已经发现了实际上有专门用于一次性乘法的函数，分别为axpy_prod和opb_prod。所以

opb_prod(A, B, C, true);

即使在未经过优化的row_major B矩阵上执行，也能在8.091秒内完成，并且与您的向量算法不相上下。

P.S. 还有更多的优化：

C = block_prod<matrix<int>, 1024>(A, B);

无论B是列主还是行主都在4.4秒内执行。 考虑这个描述：“函数block_prod适用于大型密集矩阵。”选择特定的工具来完成特定的任务！

- panda-34

---

正如我之前评论的那样，在我的机器/编译器组合（VS 9）中，完全优化时，op的boost版本实际上比向量版本更快，仅计时计算（无IO）。从反汇编来看，我猜测向量版本可能已经被gcc进行了更好的内联/优化处理，例如循环展开等。另一方面，“vector<vector>”需要进行多次分配（是否可以进行优化？），而boost可以为整个矩阵使用一个分配。 - dyp

在我的机器上，两个时间都看起来很长。对于一个1000x1000的随机矩阵，向量版本只需要1.3秒。你测试的是什么机器？ - vitaut

@vitaut，这是一台Pentium M 1600笔记本电脑 :) - panda-34

从我的测量结果来看，我现在可以得出结论：在VC9上使用iterator而不是整数+元素访问速度快3倍，但在g++(4.6.2)上并没有更快。应用于boost矩阵的ikj算法比使用vector+iterator慢（在MSVC上慢6倍，在gcc上慢10倍），无论是否使用迭代器。block_prod与两个编译器一起使用时与vector+iterator一样快。计时大致进行，不包括IO和分配。 - dyp

Answer 3

我创建了一个小网站Matrix-Matrix Product Experiments with uBLAS，它涉及将矩阵乘积的新实现集成到uBLAS中。如果您已经拥有boost库，则只需要额外的4个文件。因此，它基本上是自包含的。

我很想知道其他人是否可以在不同的机器上运行简单的基准测试。

- Michael Lehn

---

链接已失效。 - Bryan Fok