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

Question

我想将两个矩阵A和B相乘，并想比较三种不同的方法。其中之一是简单地对B的列进行迭代，然后将它们乘以矩阵A，第二个是使用犰狳的函数each_col()，并应用一个lambda，第三个简单地是乘法[ C0]。结果代码如下所示：

A * B

现在，对于

#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

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

Answer 1

这三个函数确实在做相同的事情，结果应该是相同的，除了精度差异无关紧要，因为您将结果与

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

进行比较。

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

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

您是如何安装犰狳的？

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

Armadillo可以在没有blas实现的情况下工作，在该实现中，它将使用自己的（慢速）实现进行矩阵乘法。在没有与blas链接的情况下，我还没有尝试使用它自己的实现。

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