[資料結構] 使用 C 語言：稀疏矩陣 (Sparse Matrix)，使用串列的串列 (List of List)

前言

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

• 陣列 (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 Column
    col <- n
    data <- value
    next <- null as pointer to Column

class Row
    row <- n
    nextRow <- null as pointer to Row
    nextCol <- null as pointer to Column

class Matrix
    col <- m
    row <- n
    head <- null as pointer to Row

要注意此處有兩種異質內部節點，分別儲存矩陣的 (橫) 列和 (直) 行，這和先前的串列相異；在實作時也會比一般的串列來得複雜一點。以下是等效的 C 語言宣告：

typedef struct col_node col_node_t;

struct col_node {
    double data;
    size_t col;
    col_node_t *next;
};

typedef struct row_node row_node_t;

struct row_node {
    size_t row;
    col_node_t *nextCol;
    row_node_t *nextRow;
};

typedef struct matrix matrix_t;

struct matrix {
    size_t col;
    size_t row;
    row_node_t *head;
};

稀疏矩陣的建構函式

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

static col_node_t * col_node_new(size_t col, double value)
{
    col_node_t *node = (col_node_t *) malloc(sizeof(col_node_t));
    if (!node) {
        return node;
    }

    node->col = col;
    node->data = value;
    node->next = NULL;

    return node;
}

以下是內部的 (橫) 列節點的建構函式：

static row_node_t * row_node_new(size_t row)
{
    row_node_t *node = (row_node_t *) malloc(sizeof(row_node_t));
    if (!node) {
        return node;
    }

    node->row = row;
    node->nextCol = NULL;
    node->nextRow = NULL;

    return node;
}

以下是內部的 (直) 行節點的建構函式：

matrix_t * matrix_new(size_t col , size_t row)
{
    matrix_t *mtx = (matrix_t *) malloc(sizeof(matrix_t));
    if (!mtx) {
        return mtx;
    }

    mtx->col = col;
    mtx->row = row;
    mtx->head = NULL;

    return mtx;
}

在這個實作中，共有三種異質物件，故建構函式有三個。

稀疏矩陣的解構函式

以下是 C 語言的解構函式：

void matrix_delete(void *self)
{
    if (!self) {
        return;
    }

    row_node_t *p = NULL;
    row_node_t *q = ((matrix_t *) self)->head;
    while (q) {
        col_node_t *s = NULL;
        col_node_t *t = q->nextCol;
        while (t) {
            s = t;
            t = t->next;
            free(s);
        }

        p = q;
        q = q->nextRow;
        free(p);
    }

    free(self);
}

由於每個 (橫) 列節點都會帶著一串 (直) 行節點，故需走訪兩層迴圈，逐一將內部節點釋放掉。最後再將矩陣物件本身釋放掉。

得知矩陣的大小

以下是取得矩陣大小的虛擬碼：

M is a matrix.

sub Col(M):
    return M.col

sub Row(M):
    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
    assert 0 <= i < M.col
    assert 0 <= j < M.row

    pRow <- M.head
    while (pRow != null) do
        if pRow.row == j then
            pCol <- pRow.nextCol
            while (pCol != null) do
                if pCol.col == i then
                    return pCol.data
                end if

                pCol <- pCol.next
            end while
        end if

        pRow <- pRow.nextRow
    end while

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

    row_node_t *p = self->head;
    while (p) {
        if (row == p->row) {
            col_node_t *q = p->nextCol;
            while (q) {
                if (col == q->col) {
                    return q->data;
                }

                q = q->next;
            }
        }

        p = p->nextRow;
    }

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

    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) {
                return false;
            }
        }
    }

    return true;
}

設置矩陣中任意位置的值

設置矩陣元素時，要考慮 (1) 原本元素是否為零，(2) 新設置元素是否為零；根據這些條件，共有四種情境。以下是設置矩陣中任意位置的值的虛擬碼，較長，請耐心閱讀：

M is a matrix.

sub SetAt(M, i, j, value): void
    assert 0 <= i < M.col
    assert 0 <= j < M.row

    pRow <- null
    qRow <- M.head

    if qRow == null then
        if value != 0 then
            M.head <- Row(j)
            M.head.nextCol <- Col(i, value)
        end if

        return
    end if

    while qRow != null do
        if qRow.row == j then
            pCol <- null
            qCol <- qRow.nextCol

            if qCol == null then
                if value == 0.0 then
                    if pRow != null then
                        pRow.nextRow <- qRow.nextRow
                        delete qRow
                    else
                        delete qRow
                        self.head <- null
                    end if
                else        
                    qRow.nextCol <- Col(i, value) 
                end if

                return
            end if

            while qCol != null do
                if qCol.col == i then
                    if value == 0 then
                        if pCol != null then
                            pCol.next <- qCol.next
                            delete qCol
                        else
                            delete qCol
                            q.nextCol <- null
                        end if
                    else
                        qCol.data <- value
                    end if

                    return
                end if

                pCol <- qCol
                qCol <- qCol.next
            end while

            if qCol == null then
                if value == 0 then
                    return
                end if

                qCol.next <- Col(i, value)

                return
            end if
        end if

        pRow <- qRow
        qRow <- qRow.nextRow
    end while

    if qRow == null then
        qRow.nextRow <- Row(j)
        qRow.nextRow.nextCol <- Col(i, value)
    end if

當內部元素皆為零值時，直接加入新元素即可。

當內部含非零元素時，逐一走訪 (橫) 列元素。這時會根據新設置是否為零影響後續的動作。當新設值為零時，就刪除多餘的元素；當新設值非零時，則視情形覆蓋原有值或加入新元素。細節請自行閱讀虛擬碼。

當所有內部非零元素皆不符合索引值時，就在矩陣物件尾端加入新的內部節點。

在新增 (或刪除) 內部節點時，要注意內部有兩種節點，不要誤增 (或刪) 節點。以下是等效的 C 語言程式碼：

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

    row_node_t *p = NULL;
    row_node_t *q = self->head;

    if (!q) {
        if (fabs(value) > 0.000001) {
            self->head = row_node_new(row);
            if (!(self->head)) {
                return false;
            }

            self->head->nextCol = col_node_new(col, value);
            if (!(self->head->nextCol)) {
                free(self->head);
                return false;
            }
        }

        return true;
    }

    while (q) {
        if (row == q->row) {
            col_node_t *s = NULL;
            col_node_t *t = q->nextCol;

            if (!t) {
                if (fabs(value) > 0.000001) {
                    q->nextCol = col_node_new(col, value);
                    if (!(q->nextCol)) {
                        free(q->nextCol);
                        return false;
                    }
                }
                else {
                    if (p) {
                        p->nextRow = q->nextRow;
                        free(q);
                    }
                    else {
                        free(q);
                        self->head = NULL;
                    }
                }

                return true;
            }

            while (t) {
                if (col == t->col) {
                    if (fabs(value) > 0.000001) {
                        t->data = value;
                    }
                    else {
                        if (s) {
                            s->next = t->next;
                            free(t);
                        }
                        else {
                            free(t);
                            q->nextCol = NULL;
                        }
                    }

                    return true;
                }

                s = t;
                t = t->next;
            }

            if (!t) {
                if (fabs(value) > 0.000001) {
                    s->next = col_node_new(col, value);
                    if (!(s->next)) {
                        free(s->next);
                        return false;
                    }
                }

                return true;
            }
        }

        p = q;
        q = q->nextRow;
    }

    if (!q) {
        if (fabs(value) > 0.000001) {
            p->nextRow = row_node_new(row);
            if (!(p->nextRow)) {
                return false;
            }

            p->nextRow->nextCol = col_node_new(col, value);
            if (!(p->nextRow->nextCol)) {
                free(p->nextRow);
                return false;
            }
        }

        return true;
    }

    return false;
}

轉置稀疏矩陣

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

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, j, i)

            if n != 0 then
                SetAt(N, j, i, n)
            end if
        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;
    }

    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 Col(M) == Col(N)
    assert Row(M) == Row(N)

    K <- matrix_t(Col(M), Row(M))

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

            if n != 0 then
                SetAt(K, i, j, n)
            end if
        end for
    end for

    return K

確認兩矩陣維度相等後，將兩矩陣的元素逐一相加即可。當新設值為零時，可以跳過設置值的動作，算是一個小優化。

以 C 語言實作時，我們利用 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

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

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

矩陣和純量相加

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

M, N are matrices.

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

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

            if n != 0 then
                SetAt(N, i, j, n)
            end if
        end for
    end for

    return N

按照其數學定義來操作即可。當新設置為零時，可以略去設置矩陣元素的動作。以下是等效的 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);
}

矩陣內積

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

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) - 1) do
        for j (0 to Col(N) - 1) do
            result <- 0

            for k (0 to Col(M) - 1) do
                result += At(M, k, i) * At(N, j, k)
            end for

            if result != 0 then
                SetAt(K, j, i, result)
            end if
        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(rs)
• IsEqual(M, N)：O(pqrs)
• SetAt(M, i, j, value)：O(rs)
• Trans(M)：O(pqrs)
• Add(M, N)：O(pqrs)
• Add(M, s)：O(pqrs)
• Prod(M, N)：O(pqrs)
• 使用空間：O(n)

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