位元詩人 [資料結構] 使用 C 語言:建立線性代數所使用的二維矩陣 (Matrix)

Facebook Twitter LinkedIn LINE Skype EverNote GMail Yahoo Email

二維矩陣的抽象資料結構

二維矩陣是線性代數中基本的組成單位。現在有許多程式語言或函式庫,像是 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) 來節省記憶體。我們將於後文介紹稀疏矩陣。

關於作者

身為資訊領域碩士,位元詩人 (ByteBard) 認為開發應用程式的目的是為社會帶來價值。如果在這個過程中該軟體能成為永續經營的項目,那就是開發者和使用者雙贏的局面。

位元詩人喜歡用開源技術來解決各式各樣的問題,但必要時對專有技術也不排斥。閒暇之餘,位元詩人將所學寫成文章,放在這個網站上和大家分享。