一阶矩阵乘以四阶矩阵怎么做
Posted
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了一阶矩阵乘以四阶矩阵怎么做相关的知识,希望对你有一定的参考价值。
参考技术A 矩阵相乘,就一个限定,即:m1*n1 X m2*n2(这里m1*n1表示一个二维矩阵)
要求n1=m2
即:
1*4 X 4*4 才可以。本回答被提问者采纳
将矩阵的行乘以向量(低级优化)?
【中文标题】将矩阵的行乘以向量(低级优化)?【英文标题】: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。 这是内存带宽限制,所以我认为除了使用明显的解决方案和启用矢量化进行编译(即更多的努力是过早的优化)之外,您没有太多收获。 【参考方案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;
【讨论】:
以上是关于一阶矩阵乘以四阶矩阵怎么做的主要内容,如果未能解决你的问题,请参考以下文章