建立線性代數所使用的二維矩陣 (Matrix)

    二維矩陣的抽象資料結構

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

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

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

    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);
    }
    

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

    得知二維矩陣的大小

    由於我們在建立矩陣物件時,已經儲存矩陣的維度,直接讀取其值即可。以下是等效的 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;
    }
    

    取得矩陣中任意位置的值

    先確認索引值是合法的,再進行下一步動作。由於本範例程式內部是以一維矩陣儲存元素,要將座標進行相對應的轉換。轉換的方式有 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)];
    }
    

    設置二維矩陣中任意位置的值

    先確認索引值是合法的,接著將座標進行相對應的轉換即可。由於我們先前採用 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;
    }
    

    檢查兩矩陣是否相等

    先確認兩矩陣的維度相等,再逐一確認兩矩陣的每個元素皆相等,即可確認兩矩陣相等。以下是等效的 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 時,將兩元素視為相等。讀者可視自身的需求修改本段程式碼的定義。

    矩陣轉置

    乍看程式碼略為複雜,但只是按照其數學定義來操作,追蹤一下程式碼即可了解。以下是等效的 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;
    }
    

    兩矩陣相加

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

    在實作矩陣加法時,我們利用 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;
    }
    

    矩陣和純量相加

    按照其數學定義操作即可。矩陣加法具有交換性,故只需實作一次。以下是等效的 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)

    內積的數學定義並沒有特別困難,但實作時很容易搞錯座標,需仔細追蹤程式碼。以下是等效的 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 Yahoo
    【追蹤本站】
    Facebook Facebook Twitter Plurk