[資料結構] 使用 C 語言:建立數學矩陣 (Matrix)

【分享文章】
Facebook Twitter LinkedIn LINE Skype EverNote GMail Yahoo Email

    數學矩陣的抽象資料結構

    矩陣是線性代數中基本的組成單位。現在有許多程式語言或函式庫,像是 MATLAB 或 R 等,都內建矩陣運算的功能。實務上,當我們要用 C 或 Fortran 進行向量或矩陣運算時,我們會使用 BLAS 系列函式庫,因為使用 BLAS 會比手刻來得有效率。因此,本範例程式重點在於學習其原理,而非重造輪子來用。

    以下是矩陣的抽象資料結構:

    註:由於矩陣運算項目較多,本文僅列出一些基本的矩陣操作。

    M, N, K are matrices
    s is a scalar
    
    sub Col(M): size
    sub Row(M): size
    sub At(M, i, j): result
    sub SetAt(M, i, j, value): void
    sub Trans(M): N
    sub Add(M, N): K
    sub Add(s, M): K
    sub Sub(M, N): K
    sub Sub(s, M): K
    sub Sub(M, s): K
    sub Mul(M, N): K
    sub Mul(s, M): K
    sub Div(M, N): K
    sub Div(s, M): K
    sub Div(M, s): K
    sub Prod(M, N): K
    

    以下是和矩陣的狀態相關的行為:

    • Col(M):矩陣的 (直) 行數
    • Row(M):矩陣的 (橫) 列數
    • At(M, i, j):取得矩陣中任意位置的值
    • SetAt(M, i, j, value):存入矩陣中任意位置的值

    以下是矩陣的操作:

    • Trans(M):將矩陣轉置

    根據矩陣和純量間的相對位置,可得以下三種情境:

    • Add(M, N):矩陣和矩陣相加
    • Add(s, M):純量和矩陣相加,純量在前
    • Add(M, s):矩陣和純量相加,純量在後

    由加法和乘法具有交換律,在後兩個情境者是相等的;但在減法和除法則不等。

    以下是矩陣的其他運算:

    • Prod(M, N):兩矩陣相乘後得純量

    宣告數學矩陣類別

    以下是矩陣的類別宣告:

    class Matrix
        row <- m, m > 0
        col <- n, n > 0
        elements <- a new array which size == m * n
    

    在內部儲存矩陣元素時,可用二維陣列或一維陣列來儲存;使用二維陣列時,索引值和內部陣列可以一比一對應;使用一維陣列時,索引值要經過轉換才可對應到內部的陣列。一般來說,使用一維陣列的效率會稍微好一些些,故我們採用一維陣列來儲存矩陣元素。

    以下是等效的 C 語言程式碼:

    typedef struct matrix_t matrix_t;
    
    struct matrix_t {
        size_t col;
        size_t row;
        double *elements;
    };
    

    矩陣的建構函式

    以下是 C 語言的建構函式:

    matrix_t * matrix_new(size_t col, size_t row)
    {
        assert(0 < col);
        assert(0 < row);
    
        matrix_t *mtx = malloc(sizeof(matrix_t));
        if (!mtx)
            return mtx;
    
        mtx->col = col;
        mtx->row = row;
        mtx->elements = calloc(col * row, sizeof(double));
        if (!(mtx->elements)) {
            matrix_delete(mtx);
            mtx = NULL;
            return mtx;
        }
    
        return mtx;
    }
    

    該建構函式會建立維度為 colrow 的矩陣,並將其初始化為 0.0。該建構函式是比較傳統的,我們也可以在建立矩陣物件時一併將其初始化,使用方式如下:

    matrix_t *mtx = matrix_init(3, 2,
            1.0, 2.0, 3.0,
            4.0, 5.0, 6.0
        );
    

    在這個建構函式 matrix_init 中,前兩個參數是矩陣的維度,之後的參數則是實際要填入矩陣的值。這個建構函式的實作如下:

    matrix_t * matrix_init(                     /*  1 */
        size_t col,                             /*  2 */
        size_t row,                             /*  3 */
        double value, ...)                      /*  4 */
    {                                           /*  5 */
        assert(0 < col);                        /*  6 */
        assert(0 < row);                        /*  7 */
    
        matrix_t *mtx = \
            malloc(sizeof(matrix_t));           /*  8 */
        if (!mtx)                               /*  9 */
            return mtx;                         /* 10 */
    
        mtx->col = col;                         /* 11 */
        mtx->row = row;                         /* 12 */
        mtx->elements = \
            calloc(col * row, sizeof(double));  /* 13 */
        if (!(mtx->elements)) {                 /* 14 */
            matrix_delete(mtx);                 /* 15 */
            mtx = NULL;                         /* 16 */
            return mtx;                         /* 17 */
        }                                       /* 18 */
    
        mtx->elements[0] = value;               /* 19 */
    
        va_list args;                           /* 20 */
        va_start(args, value);                  /* 21 */
    
        for (size_t j = 0; j < row; j++) {      /* 22 */
            for (size_t i = 0; i < col; i++) {  /* 23 */
                if (i == 0 && j == 0) {         /* 24 */
                    continue;                   /* 25 */
                }                               /* 26 */
    
                mtx->elements[i+col*j] = \
                    va_arg(args, double);       /* 27 */
            }                                   /* 28 */
        }                                       /* 29 */
    
        va_end(args);                           /* 30 */
    
        return mtx;                             /* 31 */
    }                                           /* 32 */
    

    第 8 行至第 18 行基本上和傳統的矩陣建構函式無異。但此建構函式在第 19 行至第 30 行間利用不定數量參數將此矩陣物件初始化。

    矩陣的解構函式

    以下是 C 語言的矩陣解構函式:

    void matrix_delete(void *self)
    {
        if (!self) {
            return;
        }
    
        double *elements = \
            ((matrix_t *) self)->elements;
        if (elements) {
            free(elements);
        }
        
        free(self);
    }
    

    先釋放內部用來儲存矩陣元素的陣列,再釋放矩陣物件本身即可。

    得知矩陣的大小

    以下是得知矩陣大小的虛擬碼:

    M is a matrix.
    
    sub Row(M): size
        return M.row
    
    sub Col(M): size
        return M.col
    

    由於我們在建立矩陣物件時,已經儲存矩陣的維度,直接讀取其值即可。以下是等效的 C 語言程式碼:

    size_t matrix_col(const matrix_t *self)
    {
        assert(self);
    
        return self->col;
    }
    
    size_t matrix_row(const matrix_t *self)
    {
        assert(self);
    
        return self->row;
    }
    

    取得矩陣中任意位置的值

    以下是取得矩陣中任意位置的值的虛擬碼:

    M is a matrix.
    
    sub At(M, i, j): result
        assert 0 <= i < M.col
        assert 0 <= j < M.row
        
        return M.elements[i + j * M.col]
    

    先確認索引值是合法的,再進行下一步動作。由於本範例程式內部是以一維矩陣儲存元素,要將座標進行相對應的轉換。轉換的方式有 column major 和 row major 兩種。轉換方式如下:

    (m, n) is the matrix dimension
    
    Column major:
    (i, j) => i + m * j
    
    Row major:
    (i, j) => i * n + j
    

    兩種轉換方式都可以使用,重點是要在整個物件內保持一致。本範例程式採用 column major。

    以下是等效的 C 語言實作:

    double matrix_at(const matrix_t *self, size_t col, size_t row)
    {
        assert(col < matrix_col(self));
        assert(row < matrix_row(self));
    
        return self->elements[col + row * matrix_col(self)];
    }
    

    設置矩陣中任意位置的值

    以下是修改矩陣中任意位置的值的虛擬碼:

    M is a matrix.
    
    sub SetAt(M, i, j, value): void
        assert 0 <= i < M.col
        assert 0 <= j < M.row
        
        M.elements[i + j * M.col] <- value
    

    先確認索引值是合法的,接著將座標進行相對應的轉換即可。由於我們先前採用 column major,這裡也保持一致。以下是等效的 C 語言實作:

    void matrix_set_at(matrix_t *self, size_t col, size_t row, double value)
    {
        assert(col < matrix_col(self));
        assert(row < matrix_row(self));
    
        self->elements[col + row * matrix_col(self)] = value;
    }
    

    檢查兩矩陣是否相等

    以下是檢查兩矩陣是否相等的虛擬碼:

    M, N are matrices.
    
    sub IsEqual(M, N): bool
        if Col(M) != Col(N) then
            return false
        end if
        
        if Row(M) != Row(N) then
            return false
        end if
        
        for i (0 to Col(M)) do
            for j (0 to Row(M)) do
                if At(M, i, j) != At(N, i, j) then
                    return false
                end if
            end for
        end for
        
        return true
    

    先確認兩矩陣的維度相等,再逐一確認兩矩陣的每個元素皆相等,即可確認兩矩陣相等。以下是等效的 C 語言實作:

    bool matrix_is_equal(const matrix_t *m, const matrix_t *n)
    {
        if (!(matrix_col(m) == matrix_col(n))) {
            return false;
        }
    
        if (!(matrix_row(m) == matrix_row(n))) {
            return false;
        }
    
        for (size_t i = 0; i < matrix_col(m); i++) {
            for (size_t j = 0; j < matrix_row(m); j++) {
                if (!(fabs(matrix_at(m, i, j) - matrix_at(n, i, j)) < 0.000001)) {
                    return false;
                }
            }
        }
    
        return true;
    }
    

    由於本範例程式以浮點數儲存元素,會有一些誤差;本範例程式在兩者誤差小於 10^-6 時,將兩元素視為相等。讀者可視自身的需求修改本段程式碼的定義。

    矩陣轉置

    以下是矩陣轉置的虛擬碼:

    M, N are matrices.
    
    sub Trans(M): N
        N <- matrix_t(Row(M), Col(M))
        
        for i (0 to Col(M)) do
            for j (0 to Row(M)) do
                SetAt(N, j, i, At(M, i, j))
            end for
        end for
        
        return N
    

    乍看程式碼略為複雜,但只是按照其數學定義來操作,追蹤一下程式碼即可了解。以下是等效的 C 語言實作:

    matrix_t * matrix_trans(const matrix_t *m)
    {
        matrix_t *out = matrix_new(matrix_row(m), matrix_col(m));
        if (!out) {
            return out;
        }
    
        for (size_t i = 0; i < matrix_col(m); i++) {
            for (size_t j = 0; j < matrix_row(m); j++) {
                matrix_set_at(out, j, i, matrix_at(m, i, j));
            }
        }
    
        return out;
    }
    

    兩矩陣相加

    以下是兩矩陣相加的虛擬碼:

    M, N, K are matrices
    
    sub Add(M, N): K
        assert M.col == N.col
        assert M.row == N.row
        
        K <- a new matrix which size == (M.col, M.row)
        for i (0 to M.col - 1) do
            for j (0 to M.row - 1) do
                SetAt(K, i, j, At(M, i, j) + At(M, i, j))
            end for
        end for
        
        return K
    

    要先確認兩矩陣的維度相等,再進行下一步動作。在本範例程式中,我們另建新的矩陣,這樣比較耗記憶體,如果想要省記憶體,可以將其相加後的值回存在第一個矩陣中。

    在實作矩陣加法時,我們利用 C11 後的泛型巨集敘述來達成多型的特性。以下是其宣告:

    #if __STDC_VERSION__ >= 201112L
    #define matrix_add(m, n) \
        _Generic((m), \
            matrix_t *: _Generic((n), \
                matrix_t *: matrix_add_mm, \
                default: matrix_add_ms), \
            default: matrix_add_sm)((m), (n))
    #else
    matrix_t * matrix_add(const matrix_t *m, const matrix_t *n);
    #endif
    

    當編譯器未支援 C11 時,我們將 matrix_add 改回一般的矩陣加法,做為退路 (fallback)。

    以下是實際的 C 語言程式碼:

    #if __STDC_VERSION__ < 201112L
    matrix_t * matrix_add(const matrix_t *m, const matrix_t *n)
    {
        return matrix_add_mm(m, n);
    }
    #endif
    
    matrix_t * matrix_add_mm(const matrix_t *m, const matrix_t *n)
    {
        assert(matrix_col(m) == matrix_col(n));
        assert(matrix_row(m) == matrix_row(n));
    
        matrix_t *out = matrix_new(matrix_col(m), matrix_row(m));
        if (!out) {
            return out;
        }
    
        for (size_t i = 0; i < matrix_col(m); i++) {
            for (size_t j = 0; j < matrix_row(m); j++) {
                matrix_set_at(out, i, j, matrix_at(m, i, j) + matrix_at(n, i, j));
            }
        }
    
        return out;
    }
    

    矩陣和純量相加

    以下是矩陣和純量相加的虛擬碼:

    M, K are matrices
    s is a scalar
    
    sub Add(s, M): K
        K <- a new matrix which size == (M.col, M.row)
        
        for i (0 to M.col - 1) do
            for j (0 to M.row - 1) do
                SetAt(K, i, j, s + At(M, i, j))
            end for
        end for
        
        return K
    

    按照其數學定義操作即可。矩陣加法具有交換性,故只需實作一次。以下是等效的 C 語言程式碼:

    matrix_t * matrix_add_ms(const matrix_t *m, double s)
    {
        matrix_t *out = matrix_new(matrix_col(m), matrix_row(m));
        if (!out) {
            return out;
        }
    
        for (size_t i = 0; i < matrix_col(m); i++) {
            for (size_t j = 0; j < matrix_row(m); j++) {
                matrix_set_at(out, i, j, matrix_at(m, i, j) + s);
            }
        }
    
        return out;
    }
    
    matrix_t * matrix_add_sm(double s, const matrix_t *m)
    {
        return matrix_add_ms(m, s);
    }
    

    矩陣內積 (Product)

    以下是矩陣內積的虛擬碼:

    M, N, K are matrices.
    
    sub Prod(M, N): K
        assert Col(M) == Row(N)
        
        K <- matrix_t(Row(M), Col(N))
        
        for i (0 to Row(M)) do
            for j (0 to Col(N)) do
                temp <- 0
                
                for k (0 to Col(M)) do
                    temp <- temp + At(M, k, i) * At(N, j, k)
                end for
                
                SetAt(K, j, i, temp)
            end for
        end for
        
        return K
    

    內積的數學定義並沒有特別困難,但實作時很容易搞錯座標,需仔細追蹤程式碼。以下是等效的 C 語言程式碼:

    matrix_t * matrix_prod(const matrix_t *m, const matrix_t *n)
    {
        assert(matrix_col(m) == matrix_row(n));
    
        matrix_t* out = \
            matrix_new(matrix_row(m), matrix_col(n));
        if (!out) {
            return out;
        }
    
        double temp;
        for (size_t i = 0; i < matrix_row(m); i++) {
            for (size_t j = 0; j < matrix_col(n); j++) {
                temp = 0.0;
    
                for (size_t k = 0; k < matrix_col(m); k++) {
                    temp += \
                        matrix_at(m, k, i) * matrix_at(n, j, k);
                }
    
                matrix_set_at(out, j, i, temp);
            }
        }
    
        return out;
    }
    

    演算法上的效率

    以下是和矩陣的狀態相關的行為:

    • Col(M)O(1)
    • Row(M)O(1)
    • At(M, i, j)O(1)
    • SetAt(M, i, j, value)O(1)

    以下是矩陣的基本四則運算:

    • Add(M, N)O(mn)
    • Add(s, M)O(mn)

    以下是矩陣操作的效率:

    • Trans(M)O(mn)

    以下是空間複雜度:

    • 空間使用率:O(mn)

    除了少數特殊的演算法可以增進一些些效率外,基本上矩陣的操作和運算的效率和其維度成正比。

    一般矩陣的實作中,矩陣中所有的值都會存在物件中,需占用較大的記憶體。有些矩陣大部分的值都為零,這時候就會用稀疏矩陣 (sparse matrix) 來節省記憶體。我們將於後文介紹稀疏矩陣。

    【分享文章】
    Facebook Twitter LinkedIn LINE Skype EverNote GMail Yahoo Email
    【追蹤網站】
    Facebook Facebook Twitter Plurk
    【支持本站】
    Buy me a coffeeBuy me a coffee