使用三种不同方法的矩阵乘法会给出不同的结果，具体取决于值的数量

Posted 2023-02-17

【中文标题】使用三种不同方法的矩阵乘法会给出不同的结果，具体取决于值的数量

【英文标题】：Matrix multiplication using three different methods gives different results, depending on amount of values

【发布时间】：2020-05-04 14:36:16

【问题描述】：

我想将A 和B 两个矩阵相乘，并想比较三种不同的方法。其中之一是简单地遍历B 的列并将它们乘以矩阵A，第二个是使用来自犰狳的函数each_col()，并应用一个lambda，第三个是简单的乘法@ 987654326@。生成的代码如下所示：

#include <complex>
#include <iostream>

#include <chrono>

#include <armadillo>

constexpr int num_values = 2048;
constexpr int num_rows = 128;
constexpr int num_cols = num_values / num_rows;

constexpr int bench_rounds = 100;

void test_multiply_loop(const arma::mat &in_mat,
                        const arma::mat &init_mat,
                        arma::mat &out_mat) 
    for(size_t i = 0; i < in_mat.n_cols; ++i) 
        out_mat.col(i) = init_mat * in_mat.col(i);
    


void test_multiply_matrix(const arma::mat &in_mat,
                          const arma::mat &init_mat,
                          arma::mat &out_mat) 
    out_mat = init_mat * in_mat;


void test_multiply_lambda(const arma::mat &in_mat,
                          const arma::mat &init_mat,
                          arma::mat &out_mat) 
    out_mat = in_mat;
    out_mat.each_col([init_mat](arma::colvec &a) 
        a = init_mat * a;
    );



int main()

    std::cout << "Hello World" << "\n";
    //Create matrix
    arma::colvec test_vec = arma::linspace(1, num_values, num_values);
    arma::mat init_mat = arma::reshape(test_vec, num_rows, num_cols);
    arma::mat out_mat_loop = arma::zeros(num_rows, num_cols),
            out_mat_lambda = arma::zeros(num_rows, num_cols),
            out_mat_matrix = arma::zeros(num_rows, num_cols);
    arma::mat test_mat = arma::eye(num_rows, num_rows);
    for(size_t i = 0; i < num_rows; ++i)
        for(size_t j = 0; j < num_rows; ++j)
            test_mat(i, j) *= (i + 1);

    auto t1 = std::chrono::high_resolution_clock::now();
    for(size_t i = 0; i < bench_rounds; ++i)
        test_multiply_loop(init_mat, test_mat, out_mat_loop);
    auto t2 = std::chrono::high_resolution_clock::now();
    auto t3 = std::chrono::high_resolution_clock::now();
    for(size_t i = 0; i < bench_rounds; ++i)
        test_multiply_lambda(init_mat, test_mat, out_mat_lambda);
    auto t4 = std::chrono::high_resolution_clock::now();
    auto t5 = std::chrono::high_resolution_clock::now();
    for(size_t i = 0; i < bench_rounds; ++i)
        test_multiply_matrix(init_mat, test_mat, out_mat_matrix);
    auto t6 = std::chrono::high_resolution_clock::now();
    std::cout << "Multiplication by loop:\t\t" << std::chrono::duration_cast<std::chrono::microseconds>( t2 - t1 ).count() << '\n';
    std::cout << "Multiplication by lambda:\t" << std::chrono::duration_cast<std::chrono::microseconds>( t4 - t3 ).count() << '\n';
    std::cout << "Multiplication by internal:\t" << std::chrono::duration_cast<std::chrono::microseconds>( t6 - t5 ).count() << '\n';
    std::cout << "Loop and matrix are equal:\t" << arma::approx_equal(out_mat_loop, out_mat_matrix, "reldiff", 0.1) << '\n';
    std::cout << "Loop and lambda are equal:\t" << arma::approx_equal(out_mat_loop, out_mat_lambda, "reldiff", 0.1) << '\n';
    std::cout << "Matrix and lambda are equal:\t" << arma::approx_equal(out_mat_matrix, out_mat_lambda, "reldiff", 0.1) << '\n';
    return 0;

现在，对于num_rows = 128，我的输出是

Multiplication by loop:         124525
Multiplication by lambda:       46690
Multiplication by internal:     1270
Loop and matrix are equal:      0
Loop and lambda are equal:      0
Matrix and lambda are equal:    0

但对于num_rows = 64，我的输出是

Multiplication by loop:         32305
Multiplication by lambda:       6517
Multiplication by internal:     56344
Loop and matrix are equal:      1
Loop and lambda are equal:      1
Matrix and lambda are equal:    1

为什么增加列数时输出会如此不同？为什么功能的时间变化如此之大？

【问题讨论】：

@idclev463035818：我刚刚重新检查了，但我仍然看不到我的错误？

代码正确。这只是令人困惑，因为与函数参数的顺序相比，您更改了矩阵乘法的另一个。但这在三个函数中是一致的，因此它们应该返回相同的输出。

@darcamo：这正是问题所在，如果num_rows 高于某个值，则对于矩阵大小的更改，我不会得到相同的结果。我也不明白为什么时间变化如此之大

对不起，误读了代码

如果你只做A * B，时间总是会更好（或相同），因为犰狳（一些 blas 实现）可以以某种“智能方式”执行乘法。许多人已经工作了几十年来优化这个操作（包括使用 CPU 的一些特定功能 - 以 mkl 库为例）。如果你做一个循环，它会起作用，但它可能会阻止使用优化的代码。这也是for_each可能比原始循环更好的原因。

【参考方案1】：

这三个函数确实在做同样的事情，结果应该是相同的，除了精度差异，因为您将结果与arma::approx_equal进行比较，所以这无关紧要。

在我的机器中，对于您提到的尺寸和我尝试过的其他更高值，输出都是正确的。我无法重现该问题。

作为参考，我尝试使用 armadillo 9.870.2 并与 openblas 和 lapack 链接。

你是如何安装犰狳的？

犰狳的大部分功能都使用了 blas 和 lapack。对于矩阵乘法，它使用一些 blas 实现。 blas 有几种实现方式，例如 openblas、mkl 甚至 cublas（用于在 gpu 中运行）等。

犰狳可以在没有 blas 实现的情况下工作，它会使用自己的（较慢的）实现来进行矩阵乘法。我没有尝试使用它自己的实现而不与 blas 链接。

可能相关的另一点是，根据 blas 实现，矩阵乘法可能会使用多个线程，但通常仅适用于大型矩阵，因为对小型矩阵使用多个线程会损害性能。也就是说，用于执行乘法运算的代码路径可能会根据矩阵大小而有所不同（但如果两个代码路径都没有产生相同的答案，这当然会是一个错误）。