C語言導航員～系統開發的最強資源～

C語言完整解析矩陣乘法｜從基礎到最佳化與應用

函式的使用方式
目次

1. 前言

C語言是廣泛用於系統開發與嵌入式程式設計的語言，矩陣計算也是其中重要的一環。本文將說明使用 C 語言實作矩陣乘法的方法，從基礎到最佳化與加速，全面覆蓋。 本文的目標讀者：
  • 正在學習 C 語言的初學者
  • 想了解矩陣計算實作方法的人
  • 對矩陣乘法的最佳化與加速有興趣的人
閱讀本文後，您將掌握在 C 語言中實作矩陣乘法的基本方法以及提升計算效率的技巧。首先，讓我們從矩陣乘法的基礎開始學習。

2. 矩陣乘積的基礎

2.1 什麼是矩陣？以易懂方式說明矩陣乘積的基本

矩陣（Matrix）是數值以縱橫排列的結構所組成的資料集合。一般而言，以具有行與列的二維陣列來表示。 例如，有以下 2×3 的矩陣：
A =
[ 1  2  3 ]
[ 4  5  6 ]
矩陣乘積是求兩個矩陣相乘的結果的運算，必須滿足以下條件：
  • 左側矩陣 A 的列數 與 右側矩陣 B 的行數 相等。
例如，將 2×3 矩陣與 3×2 矩陣相乘，會得到 2×2 的矩陣。
A =
[ 1  2  3 ]
[ 4  5  6 ]

B =
[ 7  8 ]
[ 9 10 ]
[11 12 ]
此時的矩陣乘積 C = A × B 計算方式如下。
C =
[ 58  64 ]
[139 154 ]

2.2 矩陣乘積的計算方法｜以具體例子與公式理解

矩陣乘積的計算是將每一行的元素與每一列的元素相乘並加總以求得。 一般而言，若將矩陣 A（尺寸 m × n）與矩陣 B（尺寸 n × p）相乘，得到的矩陣 C 尺寸為 m × p。 矩陣乘積的一般計算式：
C[i][j] = Σ(A[i][k] * B[k][j])  (k=0 至 n-1)
此處，
  • C[i][j] 為矩陣 C 的第 i 行第 j 列的元素
  • A[i][k] 為矩陣 A 的第 i 行第 k 列的元素
  • B[k][j] 為矩陣 B 的第 k 行第 j 列的元素

2.3 矩陣乘積的性質（結合律、分配律等）

矩陣乘積具有以下性質。
  1. 結合律（Associativity）
   (A × B) × C = A × (B × C)
但矩陣的乘法不可交換。
  1. 分配律（Distributive Law）
   A × (B + C) = A × B + A × C
  1. 單位矩陣（Identity Matrix） 單位矩陣 I 是所有對角元素為 1、其他元素為 0 的矩陣，對任意矩陣 A 成立以下關係。
   A × I = I × A = A
透過了解矩陣乘積的性質，可提升計算的最佳化與應用的廣度。

3. C語言實作矩陣乘積的方法｜附範例程式碼

3.1 C語言以二維陣列表示矩陣的方法

使用固定大小的二維陣列

以 固定大小的二維陣列 表示矩陣的方法簡單易懂，適合初學者。
#include <stdio.h>

#define ROWS 2
#define COLS 3

int main() {
    int matrix[ROWS][COLS] = {
        {1, 2, 3},
        {4, 5, 6}
    };

    // 顯示矩陣
    for (int i = 0; i < ROWS; i++) {
        for (int j = 0; j < COLS; j++) {
            printf("%d ", matrix[i][j]);
        }
        printf("n");
    }
    return 0;
}

動態記憶體配置

若想彈性處理矩陣大小，malloc() 用於動態配置記憶體。
#include <stdio.h>
#include <stdlib.h>

int main() {
    int rows = 2, cols = 3;
    int **matrix;

    // 記憶體配置
    matrix = (int **)malloc(rows * sizeof(int *));
    for (int i = 0; i < rows; i++) {
        matrix[i] = (int *)malloc(cols * sizeof(int));
    }

    // 輸入資料
    int value = 1;
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            matrix[i][j] = value++;
        }
    }

    // 顯示矩陣
    for (int i = 0; i < rows; i++) {
        for (int j = 0; j < cols; j++) {
            printf("%d ", matrix[i][j]);
        }
        printf("n");
    }

    // 釋放記憶體
    for (int i = 0; i < rows; i++) {
        free(matrix[i]);
    }
    free(matrix);

    return 0;
}

3.2 矩陣乘積的 C 語言實作｜基本程式碼說明

此處介紹 計算兩個矩陣乘積的基本程式。
#include <stdio.h>

#define N 3  // 矩陣的大小

// 計算矩陣乘積的函式
void matrixMultiplication(int A[N][N], int B[N][N], int C[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            C[i][j] = 0; // 初始化
            for (int k = 0; k < N; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}

// 顯示矩陣的函式
void printMatrix(int matrix[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%d ", matrix[i][j]);
        }
        printf("n");
    }
}

int main() {
    int A[N][N] = {
        {1, 2, 3},
        {4, 5, 6},
        {7, 8, 9}
    };

    int B[N][N] = {
        {9, 8, 7},
        {6, 5, 4},
        {3, 2, 1}
    };

    int C[N][N];  // 結果矩陣

    // 計算矩陣乘積
    matrixMultiplication(A, B, C);

    // 顯示結果
    printf("矩陣 A:n");
    printMatrix(A);
    printf("矩陣 B:n");
    printMatrix(B);
    printf("矩陣乘積 C:n");
    printMatrix(C);

    return 0;
}
在此程式碼中，matrixMultiplication() 函式使用三層迴圈，計算每個元素的乘積與加總。

4. 矩陣乘積的最佳化方法

4.1 透過迴圈最佳化加速矩陣乘積

在基本的矩陣乘積程式碼中，使用三重迴圈進行計算。然而，透過調整迴圈的順序，可以提升快取效率，並加快處理速度。

記憶體存取與快取效率

CPU 的快取在記憶體存取具備高度局部性時能有效運作。在矩陣乘積的計算中，以下的迴圈順序變更是有效的。
#include <stdio.h>

#define N 3  // 矩陣的大小

// 變更迴圈順序以改善快取效率
void optimizedMatrixMultiplication(int A[N][N], int B[N][N], int C[N][N]) {
    for (int i = 0; i < N; i++) {
        for (int k = 0; k < N; k++) {  // 將 k 放在外層
            for (int j = 0; j < N; j++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}

int main() {
    int A[N][N] = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
    int B[N][N] = {{9, 8, 7}, {6, 5, 4}, {3, 2, 1}};
    int C[N][N] = {0};

    optimizedMatrixMultiplication(A, B, C);

    // 顯示結果
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%d ", C[i][j]);
        }
        printf("n");
    }

    return 0;
}

4.2 使用 Strassen 法最佳化矩陣乘積

Strassen 法（Strassen Algorithm）是一種能將矩陣乘積的計算量從 O(N³) → O(N^{2.81}) 縮減的演算法。

Strassen 法的 C 語言實作

#include <stdio.h>
#include <stdlib.h>

// 針對 2×2 矩陣的 Strassen 法
void strassen(int A[2][2], int B[2][2], int C[2][2]) {
    int M1 = (A[0][0] + A[1][1]) * (B[0][0] + B[1][1]);
    int M2 = (A[1][0] + A[1][1]) * B[0][0];
    int M3 = A[0][0] * (B[0][1] - B[1][1]);
    int M4 = A[1][1] * (B[1][0] - B[0][0]);
    int M5 = (A[0][0] + A[0][1]) * B[1][1];
    int M6 = (A[1][0] - A[0][0]) * (B[0][0] + B[0][1]);
    int M7 = (A[0][1] - A[1][1]) * (B[1][0] + B[1][1]);

    C[0][0] = M1 + M4 - M5 + M7;
    C[0][1] = M3 + M5;
    C[1][0] = M2 + M4;
    C[1][1] = M1 - M2 + M3 + M6;
}

int main() {
    int A[2][2] = {{1, 2}, {3, 4}};
    int B[2][2] = {{5, 6}, {7, 8}};
    int C[2][2];

    strassen(A, B, C);

    // 顯示結果
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            printf("%d ", C[i][j]);
        }
        printf("n");
    }

    return 0;
}

4.3 使用 OpenMP 進行平行處理

矩陣乘積的計算因為各元素可以獨立計算，所以可以透過平行處理加速。在 C 語言中，使用 OpenMP 可以輕鬆實作平行處理。
#include <stdio.h>
#include <omp.h>

#define N 3  // 矩陣大小

// 使用 OpenMP 的平行矩陣乘積
void parallelMatrixMultiplication(int A[N][N], int B[N][N], int C[N][N]) {
    #pragma omp parallel for collapse(2)
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            C[i][j] = 0;
            for (int k = 0; k < N; k++) {
                C[i][j] += A[i][k] * B[k][j];
            }
        }
    }
}

int main() {
    int A[N][N] = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
    int B[N][N] = {{9, 8, 7}, {6, 5, 4}, {3, 2, 1}};
    int C[N][N] = {0};

    parallelMatrixMultiplication(A, B, C);

    // 顯示結果
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%d ", C[i][j]);
        }
        printf("n");
    }

    return 0;
}

5. C語言的矩陣乘法中常見的錯誤與處理方法

5.1 防止記憶體洩漏的方法

什麼是記憶體洩漏？

在進行動態記憶體配置 (malloc) 後，如果 free 釋放 沒有執行，則不需要的記憶體在程式結束後仍然被保留，會發生記憶體洩漏。這會導致記憶體使用量持續增加，程式的運作可能變得不穩定。

錯誤的程式碼範例

#include <stdio.h>
#include <stdlib.h>

int main() {
    int rows = 3, cols = 3;
    int **matrix;

    // 記憶體配置
    matrix = (int **)malloc(rows * sizeof(int *));
    for (int i = 0; i < rows; i++) {
        matrix[i] = (int *)malloc(cols * sizeof(int));
    }

    return 0;  // 因忘記 free，導致記憶體洩漏
}

處理方法

已配置的記憶體必須一定要釋放。
#include <stdio.h>
#include <stdlib.h>

int main() {
    int rows = 3, cols = 3;
    int **matrix;

    // 記憶體配置
    matrix = (int **)malloc(rows * sizeof(int *));
    for (int i = 0; i < rows; i++) {
        matrix[i] = (int *)malloc(cols * sizeof(int));
    }

    // 記憶體釋放
    for (int i = 0; i < rows; i++) {
        free(matrix[i]);
    }
    free(matrix);

    return 0;
}

5.2 陣列越界參照錯誤的原因與解決方案

什麼是陣列越界參照？

當存取陣列範圍外的記憶體時，會發生未預期的行為，情況嚴重時會出現 段錯誤（Segmentation Fault） 錯誤。

錯誤的程式碼範例

#include <stdio.h>

#define N 3

int main() {
    int matrix[N][N] = {
        {1, 2, 3},
        {4, 5, 6},
        {7, 8, 9}
    };

    // 越界存取 (當 N = 3 時，matrix[3][0] 不存在)
    printf("%dn", matrix[3][0]);  

    return 0;
}

處理方法

#include <stdio.h>

#define N 3

int main() {
    int matrix[N][N] = {
        {1, 2, 3},
        {4, 5, 6},
        {7, 8, 9}
    };

    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            printf("%d ", matrix[i][j]);
        }
        printf("n");
    }

    return 0;
}

5.3 使用除錯工具 GDB 進行錯誤分析

為了定位錯誤發生位置，使用 C 語言的除錯工具 GDB（GNU Debugger） 會很方便。

GDB 的基本使用方法

①在編譯時加入除錯資訊
   gcc -g program.c -o program
②啟動 GDB
   gdb ./program
③執行程式
   run
      1. 確認錯誤發生時的堆疊追蹤
   backtrace

使用 GDB 的實例

(gdb) run
Starting program: ./program

Program received signal SIGSEGV, Segmentation fault.
0x0000555555555142 in main () at program.c:10
10    printf("%dn", matrix[3][0]);  // 此處發生錯誤

6. 行列乘積的應用範例

6.1 圖形學中的變換矩陣

使用矩陣的 2D·3D 變換

在電腦圖形（CG）中，物件的平移（平行移動）、旋轉、縮放（放大·縮小） 等等是使用矩陣來計算的。 2D的變換矩陣 2維座標變換使用以下的 3×3矩陣。
[x']   [ a  b  tx ]   [x]
[y'] = [ c  d  ty ] * [y]
[ 1]   [ 0  0   1 ]   [1]

在 C 語言中套用 2D 變換矩陣

#include <stdio.h>
#include <math.h>

#define PI 3.14159265

// 2D 座標旋轉函式
void rotate2D(float x, float y, float angle, float *x_out, float *y_out) {
    float rad = angle * PI / 180.0;  // 角度轉換為弧度
    float rotationMatrix[2][2] = {
        {cos(rad), -sin(rad)},
        {sin(rad), cos(rad)}
    };

    // 計算矩陣乘積
    *x_out = rotationMatrix[0][0] * x + rotationMatrix[0][1] * y;
    *y_out = rotationMatrix[1][0] * x + rotationMatrix[1][1] * y;
}

int main() {
    float x = 1.0, y = 0.0;
    float angle = 90.0;
    float x_new, y_new;

    rotate2D(x, y, angle, &x_new, &y_new);
    printf("旋轉後的座標: (%.2f, %.2f)n", x_new, y_new);

    return 0;
}

6.2 機器學習中的神經網路權重計算

在神經網路中，對於輸入資料 X，將權重 W 相乘以計算下一層的輸入。
Y = X × W

在 C 語言中簡單的神經網路權重計算

#include <stdio.h>

// 計算矩陣乘積的函式
void matrixMultiply(float X[2][3], float W[3][2], float Y[2][2]) {
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            Y[i][j] = 0;
            for (int k = 0; k < 3; k++) {
                Y[i][j] += X[i][k] * W[k][j];
            }
        }
    }
}

int main() {
    float X[2][3] = {
        {1, 2, 3},
        {4, 5, 6}
    };

    float W[3][2] = {
        {0.1, 0.2},
        {0.3, 0.4},
        {0.5, 0.6}
    };

    float Y[2][2];

    matrixMultiply(X, W, Y);

    // 顯示結果
    for (int i = 0; i < 2; i++) {
        for (int j = 0; j < 2; j++) {
            printf("%.2f ", Y[i][j]);
        }
        printf("n");
    }

    return 0;
}

6.3 物理模擬中的狀態轉移

狀態轉移矩陣

例如，若 2D 物體的位置為 (x, y)、速度為 (vx, vy)，則時間 t 後的新狀態 S' 可用以下矩陣乘表示。
S' = T × S
[x']   [ 1  0  t  0 ]   [ x ]
[y'] = [ 0  1  0  t ] * [ y ]
[vx']  [ 0  0  1  0 ]   [vx ]
[vy']  [ 0  0  0  1 ]   [vy ]

C 語言的實作範例

#include <stdio.h>

// 計算狀態轉移
void stateTransition(float T[4][4], float S[4], float S_new[4]) {
    for (int i = 0; i < 4; i++) {
        S_new[i] = 0;
        for (int j = 0; j < 4; j++) {
            S_new[i] += T[i][j] * S[j];
        }
    }
}

int main() {
    float T[4][4] = {
        {1, 0, 1, 0},
        {0, 1, 0, 1},
        {0, 0, 1, 0},
        {0, 0, 0, 1}
    };

    float S[4] = {0, 0, 1, 1};  // 初始位置(0,0) 速度(1,1)
    float S_new[4];

    stateTransition(T, S, S_new);

    printf("新的狀態: %.2f %.2f %.2f %.2fn", S_new[0], S_new[1], S_new[2], S_new[3]);

    return 0;
}

7. FAQ（常見問題）

7.1 為什麼矩陣尺寸變大時計算速度會變慢？

因為計算量以 O(N³) 增加

基本的矩陣乘積計算使用三重迴圈，計算量為 O(N³)。 當矩陣尺寸變為兩倍時，計算次數會變成 8 倍，因此尺寸越大，處理時間會急劇增加。

解決方案

  1. 利用最佳化方法
  • 變更迴圈順序（提升快取效率）
  • 使用 Strassen 法（將計算量降低至 O(N^{2.81})）
  • 利用 OpenMP 等平行化技術
  1. 使用專用的數值計算函式庫
  • BLAS（Basic Linear Algebra Subprograms）
  • Intel MKL（Math Kernel Library）
  • OpenBLAS（開源高速線性代數函式庫）

7.2 Strassen 法在什麼情況下有效？

計算大規模矩陣乘積的情況下

Strassen 法是一種 矩陣尺寸越大，提速效果越顯著 的演算法。 但有以下限制。

優點

  • 計算量變為 O(N^{2.81})，可快速計算大規模矩陣乘積
  • 相較於一般的 O(N³) 方法，特別在 矩陣 N 較大時更有優勢

缺點

  • 在小矩陣上，開銷較大，反而變慢
  • 因遞迴呼叫，可能導致記憶體使用量增加

7.3 計算矩陣乘積時需注意的數值誤差是什麼？

浮點運算導致的誤差累積

在 C 語言中，使用浮點數（float 或 double）計算時，可能會產生 捨入誤差。 在矩陣乘積中，因為會進行大量的乘法與加法，這些誤差容易累積。

減少數值误差的方法

  1. 使用 double 型別
  • float（32 位元）較 double（64 位元）精度更高
  • 例：float 的精度約為 7 位數，double 約為 15 位數
  1. 利用高精度運算函式庫
  • GNU MP（GMP） 與 Intel MKL 可用於高精度運算

7.4 除了 C 語言之外的語言如何實作矩陣乘積？

Python（使用 NumPy）

在 Python 中，使用 NumPy 可以非常簡單地計算矩陣乘積。
import numpy as np

A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
B = np.array([[9, 8, 7], [6, 5, 4], [3, 2, 1]])

C = np.dot(A, B)  # 矩陣乘積

print(C)

Java（使用標準函式庫）

在 Java 中也可以 使用二維陣列實作。
public class MatrixMultiplication {
    public static void main(String[] args) {
        int[][] A = { {1, 2, 3}, {4, 5, 6}, {7, 8, 9} };
        int[][] B = { {9, 8, 7} , {6, 5, 4} , {3, 2, 1} };
        int[][] C = new int[3][3];

        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                for (int k = 0; k < 3; k++) {
                    C[i][j] += A[i][k] * B[k][j];
                }
            }
        }

        for (int i = 0; i < 3; i++) {
            for (int j = 0; j < 3; j++) {
                System.out.print(C[i][j] + " ");
            }
            System.out.println();
        }
    }
}