[資料結構] C 語言實作：稀疏矩陣 (Sparse Matrix)，使用陣列 (Array)

前言

在矩陣零值所占比率夠高時，稀疏矩陣在空間上會比傳統矩陣來得節省。一般來說，稀疏矩陣有三種實作方式：

• 陣列 (array)
• 串列的串列 (list of list)
• 映射的映射 (map of map)

本範例程式會展示第一種實作。

稀疏矩陣的抽象資料結構

稀疏矩陣的抽象資料結構如下：

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 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
sub Trans(M)

稀疏矩陣的重點在於內部實作，在外在抽象資料結構上不應和一般矩陣相異，故此處不重覆說明，讀者可看前文的說明。

宣告稀疏矩陣的類別

在內部以陣列儲存元素的前提下，其類別的宣告如下：

class Matrix
    col <- m, m is int and m > 0
    row <- n, n is int and n > 0
    size <- 0
    capacity <- 16
    cols <- a new array with size == capacity
    rows <- a new array with size == capacity
    elements <- a new array with size == capacity

稀疏矩陣內部以三個陣列來儲存元素，其中的 cols 和 rows 分別儲存矩陣元素的 (直) 行和 (橫) 列，而 elements 則儲存矩陣元素實際的值。如下圖：

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

typedef struct matrix matrix_t;

struct matrix {
    size_t col;
    size_t row;
    size_t size;
    size_t capacity;
    size_t *cols;
    size_t *rows;
    double *elements;
};

由於我們在內部會用到動態陣列，故要用 size 和 capacity 分別儲存陣列目前大小和最大容量。

稀疏矩陣的建構函式

以下是 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->size = 0;
    mtx->capacity = 2;

    mtx->cols = calloc(mtx->capacity, sizeof(size_t));
    if (!(mtx->cols)) {
        matrix_delete(mtx);
        mtx = NULL;
        return mtx;
    }

    mtx->rows = calloc(mtx->capacity, sizeof(size_t));
    if (!(mtx->rows)) {
        matrix_delete(mtx);
        mtx = NULL;
        return mtx;
    }

    mtx->elements = calloc(mtx->capacity, sizeof(double));
    if (!(mtx->elements)) {
        matrix_delete(mtx);
        mtx = NULL;
        return mtx;
    }

    return mtx;
}

程式碼看起來長了些，其實只是多建了幾個陣列，程式碼並不複雜。

陣列容量只要是比 2 大的整數即可，一般會用 16 或 32，避免在低矩陣元素時頻搬移陣列。此處故意用 2 是要測試擴展容量是否有問題。

稀疏矩陣的解構函式

以下是 C 語言的稀疏矩陣解構函式：

void matrix_delete(void *self) {
    assert(self);

    size_t *cols = ((matrix_t *) self)->cols;
    if (cols) {
        free(cols);
    }

    size_t *rows = ((matrix_t *) self)->rows;
    if (rows) {
        free(rows);
    }

    double *elements = ((matrix_t *) self)->elements;
    if (elements) {
        free(elements);
    }

    free(self);
}

同樣要先釋放內部陣列後再釋放矩陣物件本身，和先前的差別在於陣列變成三條。

取得稀疏矩陣的大小 (或維度)

以下是取得稀疏矩陣的大小 (或維度) 的虛擬碼：

M is a matrix.

sub Col(M): size
    return M.col

sub Row(M): size
    return M.row

由於在建立矩陣物件時，已儲存其維度，故直接讀取其值即可。以下是等效的 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
    for s (0 to M.size - 1) do
        if cols[s] == i and rows[s] == j then
            return elements[s]
        end if
    end for

    return 0.0

要注意稀疏矩陣取值的觀念和傳統矩陣相異。我們會逐一走訪內部的陣列，若有儲存該位置的值，就將該值回傳；反之，則回傳零。該實作隱藏在函式內部，從外部程式應無法區分兩者的差異。

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

double matrix_at(const matrix_t *self, size_t col, size_t row)
{
    assert(col < matrix_col(self));
    assert(row < matrix_row(self));

    for (size_t i = 0; i < self->size; i++) {
        if (self->cols[i] == col && self->rows[i] == row) {
            return self->elements[i];
        }
    }

    return 0.0;
}

確認兩矩陣相等

以下是確認兩矩陣相等的虛擬碼：

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) - 1) do
        for j (0 to Row(M) - 1) 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 時即相等。讀者可根據自身需求修改此處的定義。

存入稀疏矩陣中任意位置的值

這段虛擬碼分為三部分，較長，請耐心閱讀。存入稀疏矩陣時要考慮 (1) 矩陣內該位置的原始值是否為 0，(2) 存入的值是否為 0；在考量上述條件下，共有四種情境，其虛擬碼如下：

M is a matrix.

sub SetAt(M, i, j, value): void
    for s (0 to M.size - 1) do
        if cols[s] == i and rows[s] == j then
            if value != 0 then
                elements[s] <- value
                end sub
            else
                DeleteAt(M, i, j)
                return
            end if
        end if
    end for

    if value == 0 then
        return
    end if

    Expand(M)

    cols[M.size] <- i
    rows[M.size] <- j
    elements[M.size] <- value
    M.size += 1

要先逐一走訪矩陣內的非零元素，當位置相符時，若新設值不為零，直接修改即可；反之，則將該元素從非零元素中刪掉。有關 DeleteAt(M, i, j) 部分的程式碼詳見下文。

若目前矩陣在該位置沒有非零元素，而且新設值為零，就不需加入元素；反之，則需加入元素。在加入元素時，需視需求擴展內部陣列；有關 Expand(M) 部分的程式碼詳見下文。

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

static void matrix_delete_at(matrix_t *self, size_t col, size_t row);
static bool matrix_expand(matrix_t *self);

bool matrix_set_at(matrix_t *self, size_t col, size_t row, double value)
{
    assert(col < matrix_col(self));
    assert(row < matrix_row(self));

    for (size_t i = 0; i < self->size; i++) {
        if (self->cols[i] == col && self->rows[i] == row) {
            if (fabs(value) > 0.000001) {
                self->elements[i] = value;
            }
            else {
                matrix_delete_at(self, col, row);
            }

            return true;
        }
    }

    if (fabs(value) < 0.000001) {
        return true;
    }

    if (!matrix_expand(self)) {
        return false;
    }

    self->cols[self->size] = col;
    self->rows[self->size] = row;
    self->elements[self->size] = value;
    self->size++;

    return true;
}

以下是 DeleteAt(M, i, j) 部分的虛擬碼：

M is a matrix.

sub DeleteAt(M, i, j): void
    if M.size == 1 then
        M.size -= 1
        return
    end if

    k <- 0
    match <- false
    afterMatch <- false
    while k < M.size - 1 do
        if oldCols[k] == i and oldRows[k] == j then
            match <- true
        end if

        if matched then
            M.cols[k] <- M.cols[k+1]
            M.rows[k] <- M.rows[k+1]
            M.elements[k] <- M.elements[k+1]
        end if

        k += 1
    end while

    M.size -= 1

當索引值 i、j 的位置符合時，代表就是去除的目標。去除的方式是從該點開始逐一將後一個元素向前搬移 (覆寫) 即可。

此處沒有縮減內部陣列的容量，如果想要縮減容量的話，建議在陣列長度小於容量的二分之一時將陣列容量縮小一半。讀者可自行嘗試實作看看，此處不列出其程式碼。

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

static void matrix_delete_at(matrix_t *self, size_t col, size_t row)
{
    assert(col < matrix_col(self));
    assert(row < matrix_row(self));

    if (self->size <= 1) {
        self->size--;
        return;
    }

    size_t i = 0;
    bool matched = false;
    while (i < self->size - 1) {
        if (self->cols[i] == col && self->rows[i] == row) {
            matched = true;
        }

        if (matched) {
            self->cols[i] = self->cols[i+1];
            self->rows[i] = self->rows[i+1];
            self->elements[i] = self->elements[i+1];
        }

        i++;
    }

    self->size--;
}

以下是 Expand(M) 部分的虛擬碼：

M is a matrix.

sub Expand(M): void
    if M.size < M.capacity then
        return
    end if

    M.capacity <- M.capacity * 2

    oldCols <- M.cols
    newCols <- a new array which size == M.capacity

    oldRows <- M.rows
    newRows <- a new array which size == M.capacity

    oldElements <- M.elements
    newElements <- a new array which size == M.capacity

    i <- 0
    while i < M.size do
        newCols[i] <- oldCols[i]
        newRows[i] <- oldRows[i]
        newElements[i] <- oldElements[i]

        i += 1
    end while

    M.cols <- newCols
    delete oldCols
    M.rows <- newRows
    delete oldRows
    M.elements <- newElements
    delete oldElements

當內部陣列當前長度小於陣列容量時，不需擴展容量，直接結束函式即可。

當需要擴展容量時，建立三個兩倍容量的新陣列後，逐一將元素從舊陣列搬到新陣列即可。此處僅使用一般陣列，沒用到環狀陣列，索引值從 0 開始遞增即可，不需做額外處理。

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

static bool matrix_expand(matrix_t *self)
{
    assert(self);

    if (self->size < self->capacity) {
        return true;
    }

    self->capacity <<= 1;

    size_t *cols = malloc(self->capacity * sizeof(size_t));
    if (!cols) {
        return false;
    }

    size_t *rows = malloc(self->capacity * sizeof(size_t));
    if (!rows) {
        free(cols);
        return false;
    }

    double *elements = malloc(self->capacity * sizeof(double));
    if (!elements) {
        free(cols);
        free(rows);
        return false;
    }

    size_t i = 0;
    while (i < self->size) {
        cols[i] = self->cols[i];
        rows[i] = self->rows[i];
        elements[i] = self->elements[i];

        i++;
    }

    free(self->cols);
    self->cols = cols;
    free(self->rows);
    self->rows = rows;
    free(self->elements);
    self->elements = elements;

    return true;
}

稀疏矩陣的轉置

以下是稀疏矩陣轉置的虛擬碼：

M, N are matrices.

sub Trans(M): N
    N <- matrix_t(Row(M), Col(M))

    for i (0 to Col(M) - 1) do
        for j (0 to Row(M) - 1) do
            n <- At(M, i, j)

            if n == 0 then
                continue
            end if

            SetAt(N, j, i, n)
        end for
    end for

    return N

按照其數學定義操作即可，小心索引不要寫反。當新設值為 0 時，我們就跳過設置矩陣元素的動作，這算是一點點小優化。以下是等效的 C 語言實作：

matrix_t * matrix_trans(const matrix_t *m)
{
    matrix_t *out = matrix_new(matrix_row(m), matrix_col(m));
    if (!out) {
        return out;
    }

    double temp;
    for (size_t i = 0; i < matrix_col(m); i++) {
        for (size_t j = 0; j < matrix_row(m); j++) {
            temp = matrix_at(m, i, j);
            if (fabs(temp) > 0.000001) {
                matrix_set_at(out, j, i, temp);
            }
        }
    }

    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 sparse matrix which size == (M.col, M.row)

    for i (0 to M.col) do
        for j (0 to M.row) do
            a <- At(M, i, j)
            b <- At(M, i, j)

            if a + b == 0 then
                continue
            end if

            SetAt(K, i, j, a + b)
        end for
    end for

    return K

先確認兩矩陣的維度相等，再繼續下一步動作。依照其數學定義操作即可，這裡不會太困難。同樣地，當兩數和為零時，跳過該次設置矩陣的動作。

以 C 語言實作時，我們利用泛型巨集敘述來達到多型的效果：

#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

matrix_t * matrix_add_mm(const matrix_t *m, const matrix_t *n);
matrix_t * matrix_add_ms(const matrix_t *m, double s);
matrix_t * matrix_add_sm(double s, const matrix_t *m);

透過這樣的宣告，這裡的 vector_add 巨集就可以依據不同型別呼叫對應的函式。當 C 編譯器不支援這樣的特性時，以一般的兩向量相加做為其退路 (fallback)。

以下是實作的部分：

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

    double a, b;
    for (size_t i = 0; i < matrix_col(m); i++) {
        for (size_t j = 0; j < matrix_row(m); j++) {
            a = matrix_at(m, i, j);
            b = matrix_at(n, i, j);

            if (fabs(a + b) < 0.000001) {
                continue;
            }

            matrix_set_at(out, i, j, a + b);
        }
    }

    return out;
}

矩陣和純量相加

註：可能會造成效能顯著衰退。

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

sub Add(s, M): K
    K <- a new sparse matrix which size == (M.col, M.row)

    for i (0 to M.col) do
        for j (0 to M.row) do
            if s + At(M, i, j) == 0 then
                continue
            end if

            SetAt(K, 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;
    }

    double n;
    for (size_t i = 0; i < matrix_col(m); i++) {
        for (size_t j = 0; j < matrix_row(m); j++) {
            n = matrix_at(m, i, j);

            if (fabs(n + s) < 0.000001) {
                continue;
            }

            matrix_set_at(out, i, j, n + 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(Row(M), Col(N))

    for i (0 to Row(M) - 1) do
        for j (0 to Col(N) - 1) do
            temp <- 0

            for k (0 to Col(M) - 1) do
                temp <- temp + At(M, k, i) + At(N, j, k)
            end for

            if temp == 0 then
                continue
            end if

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

            if (fabs(temp) > 0.000001) {
                matrix_set_at(out, j, i, temp);
            }
        }
    }

    return out;
}

演算法上的效率

依據本文的範例程式，在演算法上的效率如下：

• Col(M)：O(1)
• Row(M)：O(1)
• At(M, i, j)：O(r)
• IsEqual(M, N)：O(pqr)
• SetAt(M, i, j, value)：O(pq)
• Trans(M)：O(pqr)
• Add(M, N)：O(pqr)
• Add(M, s)：O(pqr)
• Prod(M, N)：O(pqr)
• 使用空間：O(n)

稀疏矩陣在隨機存取的效率反而比一般矩陣差，所以這是一個用時間換取空間的實例。當稀疏矩陣的非零元素的比例過多時，稀疏矩陣就喪失了原有的優勢，這是在使用上需注意的點。