将矩阵的行乘以向量（低级优化）？答案

【问题标题】：Multiply rows of the matrix by a vector (low-level optimization)?将矩阵的行乘以向量（低级优化）？
【发布时间】：2017-09-23 16:15:30
【问题描述】：

我正在优化一个函数，我想摆脱缓慢的 for 循环。我正在寻找一种更快的方法来将矩阵的每一行乘以一个向量。

我不是在寻找“经典”乘法。

例如。我有一个包含 1024 列和 20 行的矩阵和一个长度为 1024 的向量。因此，我想要矩阵 1024 x 20，每行乘以向量。

我现在正在做的事情是在 for 循环中遍历矩阵行并使用 mkl v?Mul 执行当前矩阵行和向量的逐个元素乘法。有什么想法可以改进吗？

问题是 Multiply rows of matrix by vector? 的副本，但对于具有可能低级优化和 MKL 的 C++，而不是 R

【问题讨论】：

我假设你的意思是 1024 行和 20 列？ 20 是固定的（还是在编译时已知并保证是 4 的倍数）？你的矩阵存储的是rowmajor还是columnmajor？
@chtz 1024 列 - 这些是特征。在另一种情况下，52 列。两者都是固定的，是 4 的倍数。20 是批量大小。我选择了它，但它不能很大。而且有很多这样的乘法迭代。
我可能误解了你想要做什么。如果您想对每一行进行逐个元素的乘法运算，那么您的矩阵的列数应该与向量的元素数一样多，不是吗？（如果不是你想要的，请写一些伪代码，或者你现在使用的代码）
@chtz 对不起我的复制粘贴错误。我已经纠正了这个问题。如果我们有 1024 列矩阵，则向量大小应为 1024。
这是内存带宽限制，所以我认为除了使用明显的解决方案和启用矢量化进行编译（即更多的努力是过早的优化）之外，您没有太多收获。

标签： c++ optimization 64-bit intel-mkl avx2

【解决方案1】：

使用Eigen matrix library，您所做的实际上是乘以对角矩阵。如果您有一个任意多行和 20 列的矩阵，您可以编写以下代码（不值得为此编写函数）：

void multRows(Eigen::Matrix<double, Eigen::Dynamic, 20>& mat,
              const Eigen::Matrix<double,20,1>& vect)
{
    mat = mat * vect.asDiagonal();
}

如果编译器启用，Eigen 会生成 AVX2 代码。您可能想要试验在您的用例中存储mat 行主或列主是否更有效。

附录（由于已编辑的问题）：如果您有（很多）超过 20 列，则应该一起使用动态大小的矩阵：

void multRows(Eigen::MatrixXd& mat, const Eigen::VectorXd& vect)
{
    mat = mat * vect.asDiagonal();
}

【讨论】：

在数学上你是对的，但是：1）对角矩阵的额外内存分配是昂贵的（所有这些零）2）矩阵乘法比元素乘法更昂贵（因为缓存未命中）。所以我几乎可以肯定这会比我目前的方法慢
Eigen 不会分配一个完整的矩阵来存储对角线（事实上，该乘法根本不会发生任何分配）。这基本上针对两个元素级标量产品循环进行了优化，尽可能利用 SIMD。
你是对的，这种方法很快。同时 mat.array().rowwise() *= vect.array();相同或稍快。并且比 for loop + mkl v?Mul 快 15%
理论上，mat.array().rowwise() *= vect.array(); 行应该在（相当好）优化编译器上生成完全相同的代码（这意味着您应该使用更易读的变体）。但是对于某些表达式，您无疑会在使用不同的编译器（或不同版本的 Eigen）时体验到不同的结果。

【解决方案2】：

大多数最新的处理器都支持AVX 技术。它提供了一个包含 4 个双精度数（256 位寄存器）的向量。因此，这种优化的解决方案可能是使用 AVX。为此，我使用x86intrin.h 库实现了它，它是GCC 编译器的一部分。我还使用OpenMP 使解决方案成为多线程的。

//gcc -Wall  -fopenmp -O2 -march=native -o "MatrixVectorMultiplication" "MatrixVectorMultiplication.c" 
//gcc 7.2, Skylake Corei7-6700 HQ
//The performance improvement is significant (5232 Cycle in my machine) but MKL is not available to test
#include <stdio.h>
#include <x86intrin.h>
double A[20][1024] __attribute__(( aligned(32))) = {{1.0, 2.0, 3.0, 3.5, 1.0, 2.0, 3.0, 3.5}, {4.0, 5.0, 6.0, 6.5,4.0, 5.0, 6.0, 6.5},{7.0, 8.0, 9.0, 9.5, 4.0, 5.0, 6.0, 6.5 }};//The 32 is for 256-bit registers of AVX
double B[1024]  __attribute__(( aligned(32))) = {2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0 }; //the vector
double C[20][1024] __attribute__(( aligned(32)));//the results are stored here

int main()
{
    int i,j;
    __m256d vec_C1, vec_C2, vec_C3, vec_C4;

    //begin_rdtsc
    //get the start time here
    #pragma omp parallel for
    for(i=0; i<20;i++){
        for(j=0; j<1024; j+=16){

            vec_C1 = _mm256_mul_pd(_mm256_load_pd(&A[i][j]), _mm256_load_pd(&B[j]));
            _mm256_store_pd(&C[i][j], vec_C1);

            vec_C2 = _mm256_mul_pd(_mm256_load_pd(&A[i][j+4]), _mm256_load_pd(&B[j+4]));
            _mm256_store_pd(&C[i][j+4], vec_C2);

            vec_C3 = _mm256_mul_pd(_mm256_load_pd(&A[i][j+8]), _mm256_load_pd(&B[j+8]));
            _mm256_store_pd(&C[i][j+8], vec_C3);

            vec_C4 = _mm256_mul_pd(_mm256_load_pd(&A[i][j+12]), _mm256_load_pd(&B[j+12]));
            _mm256_store_pd(&C[i][j+12], vec_C4);

        }
    }
    //end_rdtsc
    //calculate the elapsead time

    //print the results
    for(i=0; i<20;i++){
        for(j=0; j<1024; j++){
            //printf(" %lf", C[i][j]);
        }
        //printf("\n");
    }

    return 0;
}

【讨论】：