为什么 boosts 矩阵乘法比我的慢？答案

【问题标题】：Why is boosts matrix multiplication slower than mine?为什么 boosts 矩阵乘法比我的慢？
【发布时间】：2012-06-22 01:28:08
【问题描述】：

我已经用boost::numeric::ublas::matrix 实现了一个矩阵乘法（参见my full, working boost code）

Result result = read ();

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

另一种使用标准算法（见full standard code）：

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 秒，而提升实施获得了超过 4 分钟！

edit: 用

编译后

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 比我的实现慢这么多？

【问题讨论】：

重新发明轮子是一个好主意的唯一时间是你可以制造一个更好的轮子......
Boost.uBLAS 是一个标准的接口，而不是一个强大的实现，所以不要期望它很快，除非你是使用例如LAPACK 后端。
Boost uBLAS 有一些可选的调试检查会减慢速度。请参阅此常见问题解答boost.org/doc/libs/1_49_0/libs/numeric/ublas/doc/index.htm，并检查预处理器宏 BOOST_UBLAS_NDEBUG 和 NDEBUG
虽然读取几个 2k×2k 矩阵不需要 4 分钟。
@Mysticial 棘手的部分是，大多数重新发明轮子的人通常都相信不管它是否真的更好，否则他们可能一开始就不会这样做。 :-D

标签： c++ performance boost ublas boost-ublas

【解决方案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 都比您的实现要快一些。仅测量不使用 I/O 的计算时间、删除不必要的复制以及仔细选择 block_prod（我使用 64）的块大小可以使差异更加深刻。

另请参阅 uBLAS FAQ 和 Effective uBlas and general code optimization。

【讨论】：

您能使用提供的版本 OP 运行相同的测试吗？
@mfontanini：当然，我更新了答案。请注意，我使用了较小的 (1000x1000) 随机矩阵，因此所有时间都更小。

【解决方案2】：

我相信，您的编译器优化不够。 uBLAS 代码大量使用模板，而模板需要大量使用优化。我在发布模式下通过 MS VC 7.1 编译器为 1000x1000 矩阵运行了您的代码，它给了我

10.064s 用于 uBLAS

7.851s 用于矢量

差异仍然存在，但绝不是压倒性的。 uBLAS 的核心概念是惰性评估，因此prod(A, B) 仅在需要时评估结果，例如prod(A, B)(10,100) 将立即执行，因为实际上只会计算一个元素。因此实际上~~没有可以优化的整个矩阵乘法的专用算法~~（见下文）。但是你可以帮助图书馆一点，声明

matrix<int, column_major> B;

将运行时间减少到4.426s，这会在一只手被束缚的情况下击败你的功能。该声明使矩阵相乘时对内存的访问更加有序，从而优化了缓存的使用。

附：阅读 uBLAS 文档到最后 ;)，您应该已经发现实际上有一个专用函数可以一次将 whole 矩阵相乘。 2 个功能 - axpy_prod 和 opb_prod。所以

opb_prod(A, B, C, true);

即使在未优化的 row_major B 矩阵上，也可以在 8.091 秒内执行，并且与您的向量算法相当

附言还有更多优化：

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

在4.4s 中执行，无论 B 是 column_ 还是 row_ 专业。考虑一下描述：“函数 block_prod 是为 large dense 矩阵设计的。”为特定任务选择特定工具！

【讨论】：

正如我已经评论过的，在我的机器/编译器组合 (VS 9) 上，完全优化，当仅计时计算（无 IO）时，op 的 boost 版本实际上比向量版本运行得更快。从反汇编来看，我猜向量版本可以通过 gcc 更好地内联/流线化，并使用 for 循环展开等。另一方面，vector < vector > 需要多个分配（可能进行优化？），其中 boost 可以使用一个整个矩阵。
在我看来，这两次都很大，在我的机器矢量版本上，1000x1000 随机矩阵只需要 1.3 秒。你在什么机器上测试？
@vitaut，这是一台 Pentium M 1600 笔记本 :)
+1 根据我的测量结果，我现在可以得出结论：使用iterators 代替整数 + 元素访问在 VC9 上快 3 倍，但在 g++(4.6.2) 上却不快。无论是否使用迭代器，应用于 boost 矩阵的 ikj 算法都比 vector+iterators 慢（MSVC 为 6x，gcc 为 10x）。 block_prod 与 vector+iterators 与两个编译器一样快。时间安排是粗略的，不包括 IO 和分配。

【解决方案3】：

我创建了一个小网站Matrix-Matrix Product Experiments with uBLAS。这是关于将矩阵矩阵产品的新实现集成到 uBLAS 中。如果您已经拥有 boost 库，则它仅包含额外的 4 个文件。所以它几乎是独立的。

如果其他人可以在不同的机器上运行简单的基准测试，我会很感兴趣。

【讨论】：

以上链接失效