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

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

    前言

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

    • 陣列 (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
    

    稀疏矩陣內部以三個陣列來儲存元素,其中的 colsrows 分別儲存矩陣元素的 (直) 行和 (橫) 列,而 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;
    };
    

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

    稀疏矩陣的建構函式

    以下是 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
    

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

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

    以下是等效的 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)

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

    【分享本文】
    Facebook Twitter LinkedIn LINE Skype EverNote GMail Yahoo Email
    【追蹤新文章】
    Facebook Twitter Plurk
    臉書討論區