Martix product

starlitxiling Lv3
  2025-02-04 11:19:41 2025-02-04 11:19:41 Created   2025-02-05 21:04:40 2025-02-05 21:04:40 Updated    1.5k Words  7 Mins  

主定理

主定理用于解决形如 的递归式，也就是使用了分而治之的算法(eg: 归并排序，二分查找等)，其中：

  • 是规模为 n 的问题的运行时间
  • 是子问题的数量
  • 是每个子问题规模的缩小倍数
  • 是合并步骤的复杂度指数

通过比较 和 的大小，可以得出以下结论：

  1. 当 时，时间复杂度为
  2. 当 时，时间复杂度为
  3. 当 时，时间复杂度为

矩阵乘法

先看一下最原始的矩阵乘法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
vector<vector<int>> matrixMultiply(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int rowsA = A.size();
    int colsA = A[0].size();
    int rowsB = B.size();
    int colsB = B[0].size();


    vector<vector<int>> 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 < colsA; ++k) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }

    return C;
}

这样的方法，它计算每个元素的时间复杂度达到了,为了优化矩阵乘法，可以使用分治算法将矩阵分成更小的块，通过递归地对矩阵进行分解，并应用分配律来减少计算复杂度，这里最常见的例子是Strassen算法,如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
// 矩阵加法
vector<vector<int>> add(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    vector<vector<int>> C(n, vector<int>(n));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            C[i][j] = A[i][j] + B[i][j];
    return C;
}

// 矩阵减法
vector<vector<int>> subtract(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    vector<vector<int>> C(n, vector<int>(n));
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < n; ++j)
            C[i][j] = A[i][j] - B[i][j];
    return C;
}

// 矩阵乘法
vector<vector<int>> multiply(const vector<vector<int>>& A, const vector<vector<int>>& B) {
    int n = A.size();
    if (n == 1) {
        vector<vector<int>> C(1, vector<int>(1));
        C[0][0] = A[0][0] * B[0][0];
        return C;
    }

    int mid = n / 2;
    vector<vector<int>> A11(mid, vector<int>(mid)), A12(mid, vector<int>(mid)),
                        A21(mid, vector<int>(mid)), A22(mid, vector<int>(mid));
    vector<vector<int>> B11(mid, vector<int>(mid)), B12(mid, vector<int>(mid)),
                        B21(mid, vector<int>(mid)), B22(mid, vector<int>(mid));

    // 分割矩阵
    for (int i = 0; i < mid; ++i)
        for (int j = 0; j < mid; ++j) {
            A11[i][j] = A[i][j];
            A12[i][j] = A[i][j + mid];
            A21[i][j] = A[i + mid][j];
            A22[i][j] = A[i + mid][j + mid];

            B11[i][j] = B[i][j];
            B12[i][j] = B[i][j + mid];
            B21[i][j] = B[i + mid][j];
            B22[i][j] = B[i + mid][j + mid];
        }

    // 递归计算Strassen的七个乘积
    vector<vector<int>> M1 = multiply(add(A11, A22), add(B11, B22));
    vector<vector<int>> M2 = multiply(add(A21, A22), B11);
    vector<vector<int>> M3 = multiply(A11, subtract(B12, B22));
    vector<vector<int>> M4 = multiply(A22, subtract(B21, B11));
    vector<vector<int>> M5 = multiply(add(A11, A12), B22);
    vector<vector<int>> M6 = multiply(subtract(A21, A11), add(B11, B12));
    vector<vector<int>> M7 = multiply(subtract(A12, A22), add(B21, B22));

    // 组合结果
    vector<vector<int>> C(n, vector<int>(n));
    vector<vector<int>> C11 = add(subtract(add(M1, M4), M5), M7);
    vector<vector<int>> C12 = add(M3, M5);
    vector<vector<int>> C21 = add(M2, M4);
    vector<vector<int>> C22 = add(subtract(add(M1, M3), M2), M6);

    // 合并四个子矩阵
    for (int i = 0; i < mid; ++i)
        for (int j = 0; j < mid; ++j) {
            C[i][j] = C11[i][j];
            C[i][j + mid] = C12[i][j];
            C[i + mid][j] = C21[i][j];
            C[i + mid][j + mid] = C22[i][j];
        }

    return C;
}

这里根据递归主定理，有公式
在Strassen算法中，每次分解成7个子问题，每个子问题的规模是原来的,每一层合并操作的复杂度是，所以递推关系是 ，所以根据主定理，当时，复杂度是，可以计算出：

矩阵乘法理论改进

Coppersmith-Winograd算法

该算法由Coppersmith和Winograd在1987年提出，为矩阵乘法问题提供了一个理论上的时间复杂度改进，理论上为。该算法利用张量分解技术来优化矩阵乘法。（暂时不懂这是什么，后面来填）

Alman & Williams 2020年改进

将Coppersmith-Winograd算法的时间复杂度降低到

理论下限

并行化加速（Parallel Computing）

这里介绍在实际使用程序计算时，可以做的一些优化，并非算法层面的。
矩阵乘法是典型的计算密集型任务，在分块计算的过程中，我们可以在多个处理单元上并行计算这些子矩阵的乘积，有效利用多核处理器。也可以将行或列划分到不同的处理单元，每个处理单元计算自己的行或列的结果。然后就是利用GPU并行处理的能力，通过大规模的线程并行计算加速矩阵乘法。
利用CUDA实现并行矩阵乘法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
#include <cuda_runtime.h>

__global__ void matrixMultiplyKernel(float* A, float* B, float* C, int N) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;

    if (row < N && col < N) {
        float value = 0;
        for (int k = 0; k < N; ++k) {
            value += A[row * N + k] * B[k * N + col];
        }
        C[row * N + col] = value;
    }
}

void matrixMultiply(float* A, float* B, float* C, int N) {
    float *d_A, *d_B, *d_C;
    size_t size = N * N * sizeof(float);

    cudaMalloc((void**)&d_A, size);
    cudaMalloc((void**)&d_B, size);
    cudaMalloc((void**)&d_C, size);

    cudaMemcpy(d_A, A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, B, size, cudaMemcpyHostToDevice);

    dim3 threadsPerBlock(16, 16);
    dim3 numBlocks((N + 15) / 16, (N + 15) / 16);

    matrixMultiplyKernel<<<numBlocks, threadsPerBlock>>>(d_A, d_B, d_C, N);

    cudaMemcpy(C, d_C, size, cudaMemcpyDeviceToHost);

    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
}

还可以使用分布式计算来加速计算过程，常见的分布式计算框架有MPI(Message Passing Interface)和Hadoop。
还有多线程加速，利用多核CPU的多个线程并行处理任务，适用于单台机器上拥有多个处理核心的情况。在C++中可以用std::thread或OpenMP实现多线程矩阵乘法。

  • Title: Martix product
  • Author: starlitxiling
  • Created at : 2025-02-04 11:19:41
  • Updated at : 2025-02-05 21:04:40
  • Link: http://starlitxiling.github.io/2025/02/04/Martix-product/
  • License: This work is licensed under CC BY-NC-SA 4.0.
Comments
On this page
Martix product
  1. 主定理
  2. 矩阵乘法
  3. 矩阵乘法理论改进
  4. 并行化加速（Parallel Computing）