張智星 (Roger Jang) jang@mirlab.org http://mirlab.org/jang 台大資工系 多媒體檢索實驗室 MATLAB 程式設計進階篇 稀疏矩陣 張智星 (Roger Jang) jang@mirlab.org http://mirlab.org/jang 台大資工系 多媒體檢索實驗室

5-1 稀疏矩陣的建立 根據內部資料結構的不同，MATLAB 矩陣分為兩種 : 完全矩陣（Full Matrix） 5-1 稀疏矩陣的建立 根據內部資料結構的不同，MATLAB 矩陣分為兩種 : 完全矩陣（Full Matrix） 每個元素都有自己的儲存空間 每個元素的佔用空間大小依其資料型態而有所不同 double8 bytes, int324 bytes, int162 bytes, int81 bytes 一個 m×n 的完全矩陣，所佔用的記憶體空間是 k×m×n bytes 稀疏矩陣（Sparse Matrix） 大部分的元素都是 0 只須儲存「非零元素的位置」及其「元素值」 優點 節省記憶體儲存空間稀疏矩陣的最重要任務！ 節省某一些不必要的運算

稀疏矩陣範例-1 (I) sparse 指令可將一個完全矩陣轉換成稀疏矩陣 範例5-1：sparse01.m (1,1) 2 (3,2) 4 (2,4) 1 可看出，S 是一個稀疏矩陣，MATLAB 只儲存其各個非零元素的位置（即其二維下標 (1,1)、(3,2)、(2,4)）和元素值（2、4、1） clear all; % 清除所有的變數 A = [2 0 0 0; 0 0 0 1; 0 4 0 0]; % 完全矩陣 S = sparse(A) % 將完全矩陣 A 轉換成稀疏矩陣 S

稀疏矩陣範例-1 (II) 比較矩陣 A 和 S 所佔用的記憶體大小 : >> whos Name Size Bytes Class Attributes A 3x4 96 double S 3x4 88 double sparse 看出稀疏矩陣 S 佔用記憶體的位元組數目（88 Bytes）比 矩陣A（佔用 96 Bytes）還要小。 當矩陣變大、密度變低時，這個差異會越來越大。

稀疏矩陣範例-2(I) S = sparse(i, j, s, m, n); 使用 sparse 指令來直接產生稀疏矩陣 i 是列索引 m 是 s 的列維度 n 是 s 的行維度 i、j、s 都是長度相同的向量 s(k) 的二維下標即是 i(k)及 j(k)

稀疏矩陣範例-2(II) >> S = sparse([1 3 2], [1 2 4], [2 4 1], 3, 4) S = (1,1) 2 (3,2) 4 (2,4) 1 也可以在 sparse 指令加上第六種輸入變數，代表最多可以容納的非零元素個素，使得您可以後續再加入非零元素，而不必改變整個稀疏矩陣的結構

稀疏矩陣範例-3(I) 指令 spdiags，可由對角線元素來建構一個稀疏矩陣 S = spdiags(D, p, m, n)

稀疏矩陣範例-3(II) 範例5-2：sparse02.m S = D = [3 2 9; 2 4 9; 1 1 4]; (1,1) 2 (2,1) 9 (2,2) 4 (3,2) 9 (1,3) 1 (3,3) 1 (4,3) 4 D = [3 2 9; 2 4 9; 1 1 4]; d = [2 0 -1]; S = spdiags(D, d, 4, 3)

稀疏矩陣範例-3(III) 以完全矩陣來顯示，可得 : >> A = full(S) 2 0 1 9 4 0 0 9 1 0 0 4 可看出矩陣 A 的三個非零對角線向量分別是 D 的三個行向量。（在D的第一個行向量中，只用到最後一個元素1。） 提示： D = 3 2 9 2 4 9 1 1 4 d = 2 0 -1

稀疏矩陣範例-4(I) 一般的 load 及 save 指令，也可以處理稀疏矩陣，並儲存於二進制的 MAT 檔案 spconvert指令，可將一個 m×3 的矩陣轉換成稀疏矩陣 第一直行代表列索引 第二直行代表行索引 第三直行則是非零的元素值 檔案 spmat.dat 的內容可顯示如下： >> type spmat.dat 1 1 2 3 2 4 2 4 1 8 3 5

稀疏矩陣範例-4(II) 建構此稀疏矩陣 範例5-3：spconvert01.m S = load spmat.dat (1,1) 2 (1,1) 2 (3,2) 4 (8,3) 5 (2,4) 1 load spmat.dat S = spconvert(spmat)

5-2 稀疏矩陣的儲存空間 一個只包含實數的稀疏矩陣，假設其維度為 m×n，含有 nnz 個非零元素，MATLAB 動用了三個內部陣列來儲存此稀疏矩陣的相關資訊： 第一個陣列：以 double 方式儲存了所有的非零元素，其長度為 nnz，使用的空間為大小 8*nnz 位元組（Bytes）。 第二個陣列：以整數方式儲存了每個元素的列索引，其長度為 nnz，使用的空間大小為 4*nnz 位元組（Bytes）。 第三個陣列：以整數方式儲存了每個直行的起始指標，其長度為 n，使用空間大小為 4*n 位元組（Bytes）。 注意：上述儲存方式僅是用於MATLAB第六版！

稀疏矩陣的儲存空間 MATLAB第七版後，整個稀疏矩陣佔用的空間大小為 16*nnz + 8*n + 8，可以驗證如下。 範例5-4：memorySize01.m totalSize = 34032 Name Size Bytes Class Attributes west0479 479x479 34032 double sparse clear all load west0479.mat [m,n]=size(west0479); totalSize=16*nnz(west0479)+8*n+8 whos west0479

稀疏矩陣的儲存空間 Regarding the formula “16*nnz + 8*n + 8” Can be verified via least-square estimate in memorySize02.m Open question: Why is it so?

稀疏矩陣的相關指令 其他相關指令： nnz：可傳回稀疏矩陣的非零元素個數 nonzeros：傳回一個包含所有非零元素的行向量 nzmax：傳回最大的非零元素個數

提示 在一個稀疏矩陣中，將一個非零元素設定成零時，MATLAB 並不會自動釋放記憶體空間，換包話說，nnz 會隨著非零元素的減少而減少，但 nzmax 並不會隨著改變 但是多加一個非零元素時，若 nnz 已大於 nzmax 時，nzmax 會隨之增大（即 MATLAB 會自動配置記憶體）以儲存新加的元素

5-3 稀疏矩陣的觀看與圖示 spy 指令可用於觀看稀疏矩陣的非零元素分佈情況 範例5-5：spy01.m 5-3 稀疏矩陣的觀看與圖示 spy 指令可用於觀看稀疏矩陣的非零元素分佈情況 範例5-5：spy01.m load west0479.mat % 載入二進位制檔案 west0479.mat spy(west0479) % 觀看稀疏矩陣的非零元素分佈情況

稀疏矩陣的觀看與圖示 矩陣 west0479 的維度是 479*479，但是只包含 1887 個非零元素，因此此矩陣的密度只有 1887/(479*479) = 0.0082

稀疏矩陣的觀看與圖示 稀疏矩陣特別適用於表示一個「無向圖」（Undirected Graph）的「鄰近矩陣」（Adjacency Matrix） 若某圖的第 i 和第 j 個節點有直線連接，則其相對應的鄰近矩陣在第 i 列、第 j 行的元素值為 1，其他元素值則為零 為節省空間，僅用上三角或是下三角矩陣來表示

稀疏矩陣的觀看與圖示 對應的鄰近矩陣可表示成： >> A = spconvert([1 2 1; 2 3 1; 2 4 1; 3 4 1; 3 5 1; 5 6 1; 4 6 1]); >> full(A) ans = 0 1 0 0 0 0 0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 0 1

稀疏矩陣的觀看與圖示 假設這 6 個節點的座標如下： 可用 gplot 指令來畫出上述的無向圖： 範例5-6：gplot01.m >> xy = [0 1; 1 2; 1 0; 2 0; 2 2; 3 1]; % 每個列向量是一組座標 可用 gplot 指令來畫出上述的無向圖： 範例5-6：gplot01.m 其中 '-o' 代表以實線（'-'）及圓圈（'o'）來作圖 A = spconvert([1 2 1; 2 3 1; 2 4 1; 3 2 1; 3 4 1; 3 5 1; 4 2 1; 4 3 1; 4 6 1; 5 3 1; 5 6 1; 6 4 1; 6 5 1]); xy = [0 1; 1 2; 1 0; 2 0; 2 2; 3 1]; % 每一個列向量是一組 (x, y) 座標 gplot(A, xy, '-o') % 畫出無向圖（Undirected Graph）

Adjacency Matrix Real-world examples 台北捷運：n>100

稀疏矩陣有趣的例子-1 Bucky 球，此圖包含了 60 個三度空間中的點，每一點和他的三個鄰近點都是等距離，可用 bucky 指令來產生這些點的鄰近矩陣，並用 gplot 來顯示圖形 範例5-7：gplot02.m [A, xy] = bucky; % A 為鄰近矩陣，xy 為座標 gplot(A, xy, '-o'); % 畫出無向圖（Undirected Graph） axis equal % 設定 x 軸和 y 軸的刻度一樣）

畫出一棵電腦圖學中的樹 treeplot指令可畫出一棵電腦圖學中的樹 範例5-8：treePlot01.m 使用 parent 向量來代表這一棵樹，其中 parent(1)=0 則代表第一個節點是此樹的根節點（Root），而 parent(i)=j 代表第 i 個節點的父親是第 j 個節點，例如 parent(5)=4 代表第5個節點的父親是第 4 個節點，依此類推 parent = [0 1 2 2 4 4 4 1 8 8 10 10 11 11 11 11]; treeplot(parent);

稀疏矩陣有趣的例子-2 是由 NASA （美國太空總署）所主導的計畫，其中包含計算流過機翼的氣流所造成的作用力，由於必須進行偏微分方程的數值運算，所以必須對二維空間進行三角化切割，其鄰近矩陣即為一個稀疏矩陣，您可在 MATLAB 下執行 airfoil 指令即可產生相關圖形，相關說明則記載於 airfoil.m，可經由 type airfoil.m 來檢視之。

5-4 稀疏矩陣的運算 MATLAB 針對完全矩陣設計的運算與函式，也都適用於稀疏矩陣，而且輸出也是大部分以稀疏矩陣的方式來表示 5-4 稀疏矩陣的運算 MATLAB 針對完全矩陣設計的運算與函式，也都適用於稀疏矩陣，而且輸出也是大部分以稀疏矩陣的方式來表示 若計算過程包含稀疏及完全矩陣，則計算結果的表示方式就依情況而變，其規則可見 MATLAB 線上輔助說明

稀疏矩陣的運算 稀疏矩陣的運算還包含下列幾種： 排列（Permutation）及重排（Reordering） 因子分解（Factorization） 線性聯立方程式的求解 固有值及奇異值的計算 這方面的運算，牽涉到很多數學理論，有興趣的讀者，可參考 MATLAB 的手冊，或在 MATLAB 指令視窗下輸入「help sparfun」以取得線上輔助說明

5-5 本章指令彙整 指令 功能 B = sparse(A) 將一個完全矩陣 A 轉換成稀疏矩陣 B 5-5 本章指令彙整 指令 功能 B = sparse(A) 將一個完全矩陣 A 轉換成稀疏矩陣 B B = sparse(i, j, s, m , n ) 直接產生一個稀疏矩陣 B，其中 ： i 是列索引 j 是行索引 s 是非零元素形成的向量 m 是 B 的列維度 n 是 B 的行維度 i、j、s 都是長度相同的向量 B = spdiags(D, p, m, n) 由矩陣 D 的對角線元素來建構一個稀疏矩陣 B，其中： D 的每一個直行代表矩陣的對角線向量 p 代表對角線的位置（0 代表主對角線，-1 代表向下位移一單位的次對角線，1 代表向上位移一單位的次對角線，依此類推） m 與 n 則分別代表矩陣的列維度與行維度 full(B) 以完全矩陣來顯示矩陣 B

5-5 本章指令彙整 spconvert(C) 將一個 m×3 的矩陣 C，轉換成稀疏矩陣，其中： 第一直行代表列索引 第二直行代表行索引 5-5 本章指令彙整 spconvert(C) 將一個 m×3 的矩陣 C，轉換成稀疏矩陣，其中： 第一直行代表列索引 第二直行代表行索引 第三直行則是非零的元素值 nnz(C) 傳回稀疏矩陣 C 的非零元素個數 nonzeros(C) 傳回稀疏矩陣 C 的所有非零元素形成的行向量 nzmax(C) 傳回稀疏矩陣 C 的目前可容納非零元素個數的最大值，當nnz > nzmax 時，MATLAB 會動態調增配置記憶體給 nzmax，以儲存新增的非零元素 spy(C) 觀看稀疏矩陣 C 的非零元素分佈情況 gplot(A, xy, '-o'); 畫出無向圖（Undirected Graph） treeplot(nodes) 畫出樹狀圖