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

Question

我想将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

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

Answer 1

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

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

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

你是如何安装犰狳的？

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

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

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