都知道计算机执行矩阵相乘运算很麻烦,使用传统数学算法的时间复杂度是O(n^3)。
这里不讨论数学上的优化算法。在计算机计算时有个cache-miss问题,这里讨论一下。
C++里的矩阵,其实是二维数组。在存储的时候是按行存的,cache在读取的时候,也是按行取的(下面的代码可以证明)。如果按照正常算法执行矩阵相乘,依次计算新矩阵每个位置的结果,用第一个矩阵的行去乘第二个矩阵的列,然后累加求和,那第二个矩阵每次运算都跨行了,这涉及到3层循环,第1层表示行,第2层表示列,第3层表示累加。那第3层循环里,每次都会有cache-miss。
如果能设计一个算法,让第3层循环的不要跨行,而是一直在执行某一行的计算,该行的所有元素都取完了再去执行下一行,那就可以避免cache-miss了。
这里给出3组方法,第1种是原始的方法。第2和第3种方法调整了循环顺序(reordering)。
第二个方法是重排序V1版本,第三种方法是重排序的V2版本。
2个版本的区别是:最内层循环的读取矩阵的顺序,V1是逐行读取,V2是逐列读取。
#include <iostream>
#include <vector>
#include <chrono>
using namespace std;
typedef vector<vector<int>> Matrix;
Matrix matrixMultiplyReorderedV2(const Matrix& A, const Matrix& B) {
int rowsA = A.size();
int rowsB = B.size();
int colsB = B[0].size();
Matrix C(rowsA, vector<int>(colsB, 0)); // 初始化结果矩阵
for (int k = 0; k < rowsB; k++) {
for (int j = 0; j < colsB; j++) {
for (int i = 0; i < rowsA; i++) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
return C;
}
Matrix matrixMultiplyReorderedV1(const Matrix& A, const Matrix& B) {
int rowsA = A.size();
int rowsB = B.size();
int colsB = B[0].size();
Matrix C(rowsA, vector<int>(colsB, 0)); // 初始化结果矩阵
for (int k = 0; k < rowsB; k++) {
for (int i = 0; i < rowsA; i++) {
for (int j = 0; j < colsB; j++) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
return C;
}
Matrix matrixMultiplyOriginal(const Matrix& A, const Matrix& B) {
int rowsA = A.size();
int rowsB = B.size();
int colsB = B[0].size();
Matrix C(rowsA, vector<int>(colsB, 0)); // 初始化结果矩阵
for (int i = 0; i < rowsA; i++) {
for (int j = 0; j < colsB; j++) {
for (int k = 0; k < rowsB; k++) {
C[i][j] += A[i][k] * B[k][j];
}
}
}
return C;
}
int main() {
int sz = 1000;
int rowsA = sz, rowsB = sz, colsB = sz; // 矩阵大小
Matrix A(rowsA, vector<int>(rowsB, 1));
Matrix B(rowsB, vector<int>(colsB, 1));
auto start = chrono::high_resolution_clock::now();
Matrix C = matrixMultiplyOriginal(A, B);
auto end = chrono::high_resolution_clock::now();
chrono::duration<double> duration = end - start;
cout << "Original version time: " << duration.count() << " seconds" << endl;
start = chrono::high_resolution_clock::now();
C = matrixMultiplyReorderedV1(A, B);
end = chrono::high_resolution_clock::now();
duration = end - start;
cout << "matrixMultiplyReorderedV1 version time: " << duration.count() << " seconds" << endl;
start = chrono::high_resolution_clock::now();
C = matrixMultiplyReorderedV2(A, B);
end = chrono::high_resolution_clock::now();
duration = end - start;
cout << "matrixMultiplyReorderedV1 version time: " << duration.count() << " seconds" << endl;
return 0;
}
运行结果是:
Original version time: 19.3017 seconds
matrixMultiplyReorderedV1 version time: 11.8494 seconds
matrixMultiplyReorderedV1 version time: 22.6975 seconds
可以看到,V1版本明显缩短了运算时间,减少了cache-miss。V2版本没效果,反而比原始版本的还差一点。这可以说明:cache是按行读取内存的矩阵的,而不是列。
另外,使用多线程可以缩短计算时间,我在另一篇文章(多线程实现矩阵相乘_C++)里有详细说明。