稀疏矩陣稀疏矩陣是一種特殊型別的矩陣，即矩陣中包括較多的零元素。對於稀疏矩陣的這種特性，在MATLAB中可以只儲存矩陣中非零元素及非零元素在矩陣中的位置。在用稀疏矩陣進行計算時，通過消去零元素可以減少計算的時間。7.1 稀疏矩陣的儲存方式對一般矩陣而言，MATLAB儲存矩陣內的每一個元素，矩陣中的零元素與其他元素一樣，需要佔用同樣大小的記憶體空間。但對於稀疏矩陣，MATLAB僅儲存稀疏矩陣中的非零元素及其對應的位置，其他空餘位置只是在訪問時以預設的零元素來填充。對於一個含有大量零元素的大型矩陣，採用這種方法可以大大地減少資料佔據的記憶體空間。MATLAB採用3個內部陣列來儲存元素為實數的稀疏矩陣。稀疏矩陣也可用於儲存複數。當稀疏矩陣用於儲存複數資料時，需用第4個內部陣列儲存非零複數的虛部。一個複數非零，是指其實部或虛部至少其中一個不為零。【例2-16】 稀疏矩陣與一般矩陣的記憶體佔用對比。>> M_full = magic(1100); % 建立一個1100´1100 矩陣>> M_full(M_full > 50) = 0; % 將>50的元素設定為0>> M_sparse = sparse(M_full); % 建立稀疏矩陣>> whos Name Size Bytes Class Attributes M_full 1100x1100 9680000 double M_sparse 1100x1100 9608 double sparse 本例中，M_full和 M_sparse兩個變數儲存的實際上是同一個矩陣，但是二者因為採用的儲存形式分別為一般矩陣和稀疏矩陣，所以佔用的記憶體量卻相差了近1000倍。因為MATLAB版本不同，作業系統不同（例如32位和64位），內部儲存格式也有些變化，但總體上來說佔用的記憶體空間比一般矩陣小很多。7.2 稀疏矩陣的建立MATLAB決不會自動地建立一個稀疏矩陣，這需要使用者來決定是否使用稀疏矩陣。在建立一個矩陣前，使用者需要根據此矩陣中是否包含較多的零元素，採用稀疏矩陣技術是否有利，來決定是否採用稀疏矩陣的形式。把矩陣中非零元素的個數除以所有元素的個數，就叫做矩陣的密度，密度越小的矩陣採用稀疏矩陣的格式越有利。要將一般矩陣轉換為稀疏矩陣，可以使用函式sparse，如s=sparse (A)，是指將矩陣A轉換為稀疏矩陣。另外，使用函式full則可把稀疏矩陣轉換為一般矩陣。【例2-17】 一般矩陣與稀疏矩陣的轉換示例。>> A=[0 0 0 1;0 1 0 0;1 2 0 0;0 0 3 0]A = 0 0 0 1 0 1 0 0 1 2 0 0 0 0 3 0>> s=sparse(A)s = (3,1) 1 (2,2) 1 (3,2) 2 (4,3) 3 (1,4) 1>> B=full(s)B = 0 0 0 1 0 1 0 0 1 2 0 0 0 0 3 0從本例的結果中可以看出所有s的非零元素列表及其對應的行列序號。所有非零元素儲存在一列中，反映了資料的內部結構。稀疏矩陣的建立一般有以下幾種方式。1．直接建立稀疏矩陣使用函式sparse，可以用一組非零元素直接建立一個稀疏矩陣。該函式呼叫格式為：S=sparse(i,j,s,m,n)其中i和j都為向量，分別是指矩陣中非零元素的行號與列號，s是一個全部為非零元素向量，元素在矩陣中排列的位置為(i,j)。m為輸出的稀疏矩陣的行數，n為輸出的稀疏矩陣的列數。【例2-18】 稀疏矩陣的建立。>> S=sparse([1 3 2 1 4],[3 1 4 1 4],[1 2 3 45],4,4)S = (1,1) 4 (3,1) 2 (1,3) 1 (2,4) 3 (4,4) 5>> full(S)ans = 4 0 1 0 0 0 0 3 2 0 0 0 0 0 0 5本例中通過sparse函式直接建立了稀疏矩陣S。sparse函式中的前兩個輸入變數[1 3 2 1 4]和[3 1 4 1 4]就是元素在矩陣中排列的位置，第3個輸入變數[1 2 3 4 5]就是稀疏矩陣前面兩個輸入變數中的位置所對應的元素的值，而最後的兩個輸入變數4是指輸出的稀疏矩陣的行數是4，輸出的稀疏矩陣的列數同樣也為4。通過full函式把稀疏矩陣轉換為一般矩陣，這樣就可以清楚地看出sparse函式輸入和輸出之間的關係。需要指出的是：函式sparse還有一個變化形式，可以設定最大數目的非零元素。如有必要，可以在函式sparse中新增第6個輸入引數，設定稀疏矩陣中非零元素的最大數目，這樣以後要在矩陣中新增非零元素，就無需再修改矩陣的結構。具體的使用方法請查閱help文件。2．從對角線元素中建立稀疏矩陣要將一個矩陣的對角線元素儲存在一個稀疏矩陣中，可以使用函式spdiags。其呼叫格式為：S=spdiags(B,d,m,n)函式spdiags用於建立一個尺寸為m行n列的稀疏矩陣S，其非零元素來自矩陣B中的元素且按對角線排列，引數d指定矩陣B中用於生成稀疏矩陣B的對角線位置。矩陣的主對角線可以認為是第0條對角線，每向右上移動一條對角線編號加1，向左下移動一條對角線編號減1，也就是說B中的j列填充向量d元素，j指定對角線。【例2-19】 稀疏矩陣的建立。>> B=[1 2 3;4 5 6;7 8 9;10 11 12]B = 1 2 3 4 5 6 7 8 9 10 11 12>> d=[-3 0 2]d = -3 0 2>> A=spdiags(B,d,7,4) A = (1,1) 2 (4,1) 1 (2,2) 5 (5,2) 4 (1,3) 9 (3,3) 8 (6,3) 7 (2,4) 12 (4,4) 11 (7,4) 10>> full(A)ans = 2 0 9 0 0 5 0 12 0 0 8 0 1 0 0 11 0 4 0 0 0 0 7 0 0 0 0 10本例生成了一個7行4列的稀疏矩陣A。B的第1列元素排列在主對角線以下的第3條對角線上，第2列元素排列在主對角線上，第3列中的非零元素排列在主對角線上方的第2條對角線上。這裡需要注意B中的第三列並沒有全部分佈在第2條對角線上，而是最後兩個元素9和12排列在該對角線上。3．從外部檔案中匯入稀疏矩陣用外部檔案建立的文字檔案，如果該檔案中的資料按3或者4列排列，則可將這個文字檔案載入記憶體，用於建立一個稀疏矩陣。【例2-20】 稀疏矩陣的建立。假如有這樣一個檔案uphill.dat（使用者可以通過記事本開啟、編輯其內容），檔案內含有以下資料：1 1 1.0000000000000001 2 0.5000000000000002 2 0.3333333333333331 3 0.3333333333333332 3 0.2500000000000003 3 0.2000000000000001 4 0.2500000000000002 4 0.2000000000000003 4 0.1666666666666674 4 0.1428571428571434 4 0.000000000000000那麼通過使用load命令可以將此資料檔案載入MATLAB，然後對其進行操作。實際中使用者可以在命令視窗輸入：>> load uphill.dat % 用load命令將資料的文字檔案uphill.dat載入工作空間>> H = spconvert(uphill)H = (1,1) 1.0000 (1,2) 0.5000 (2,2) 0.3333 (1,3) 0.3333 (2,3) 0.2500 (3,3) 0.2000 (1,4) 0.2500 (2,4) 0.2000 (3,4) 0.1667 (4,4) 0.1429> >> full(H)ans = 1.0000 0.5000 0.3333 0.2500 0 0.3333 0.2500 0.2000 0 0 0.2000 0.1667 0 0 0 0.1429本例首先使用load函式匯入了一個3列資料的文字檔案uphill.dat，使用者可以通過在命令列中輸入變數名uphill 可以檢視資料uphill 中的具體內容來驗證資料讀取是否正確，然後呼叫spconvert將uphill 轉換為相應的稀疏矩陣H。通過呼叫full函式可以直觀地看出得到的稀疏矩陣。MATLAB使用load函式來匯入外部資料檔案，具體的用法可以參閱第10章。7.3 稀疏矩陣的運算多數MATLAB自帶的數學函式都可用於處理稀疏矩陣，此時可以將稀疏矩陣當做一般矩陣看待。另外MATLAB也提供有一些專門針對稀疏矩陣的函式。處理稀疏矩陣時，計算的複雜程度與稀疏矩陣中非零元素的數目成正比，也與矩陣的行列尺寸有關。比如稀疏矩陣的乘法、乘方、包含一定次數的線性方程組等，都是比較複雜的運算。用函式處理稀疏矩陣時，計算結果要遵循以下一些原則。（1）MATLAB函式處理一個矩陣時，不管這個矩陣是一般矩陣還是稀疏矩陣，其返回值為一個數值或矩陣。返回值都按一般矩陣方式進行儲存，並不會根據接受的引數是稀疏矩陣，而將結果儲存為稀疏矩陣。（2）函式處理一個數值或向量返回一個矩陣時，如果矩陣為零矩陣、元素全為1的矩陣、隨機矩陣或單位矩陣，這些矩陣全為一般矩陣形式。對於零矩陣，有一種類似稀疏矩陣的儲存方法，因為零矩陣中沒有非零元素，所以不能將一個零矩陣轉換為一個稀疏矩陣，但指令zeros(m,n)和sparse(m,n)是可用的。對於單位矩陣和隨機矩陣，可以使用類似稀疏矩陣的操作指令，即speye和sprand。對於元素全為1的矩陣，則沒有類似的操作指令。（3）以矩陣為引數返回矩陣或向量的一元函式，返回值的儲存型別與引數的儲存型別相同。例如矩陣S的cholesky分解，如果S為一般矩陣，結果也為一般矩陣；如果S為稀疏矩陣，結果也為稀疏矩陣。對於列向處理矩陣的函式，如求各列最大值的函式max，求各列之和的函式sum等，也都返回與引數相同的儲存型別。如果引數是稀疏矩陣，即使返回的矩陣或向量全為非零元素，也用稀疏方式表示。例外情況只有函式sparse和full，因為它們用於一般矩陣和稀疏矩陣之間的轉換。（4）對於有兩個輸入引數為矩陣的情況，如果輸入的兩個矩陣都為稀疏矩陣，則輸出仍為稀疏矩陣；都為一般矩陣，結果也為一般矩陣。如果輸入引數一個為稀疏矩陣，一個為一般矩陣，結果通常為一般矩陣，但在能夠保證矩陣稀疏性不變的運算中，結果則為稀疏矩陣。（5）使用方括號對矩陣進行組合時，如果組合的矩陣中有稀疏矩陣，結果則為稀疏矩陣。（6）子矩陣在右邊的賦值操作，返回值為右邊子矩陣的儲存型別，子矩陣在左邊賦值不改變其儲存型別。【例2-21】 稀疏矩陣的組合。>> A=[1 0 0;0 0 1;1 2 0]A = 1 0 0 0 0 1 1 2 0>> B=sparse(A)B = (1,1) 1 (3,1) 1 (3,2) 2 (2,3) 1>> C=[A(:,1),B(:,2)]C = (1,1) 1 (3,1) 1 (3,2) 2本例將矩陣A的第1列和矩陣B的第2列組成了新的矩陣C，從結果可知，C為稀疏矩陣。【例2-22】 稀疏矩陣子矩陣的賦值。>> A=[1 0 0;0 0 1;1 2 0];>> B=sparse(A);>> C=sparse(cat(1,full(B),A))C = (1,1) 1 (3,1) 1 (4,1) 1 (6,1) 1 (3,2) 2 (6,2) 2 (2,3) 1 (5,3) 1>> i=[1 2 3];>> j=[1 2 3];>> T=C(i,j)T = (1,1) 1 (3,1) 1 (3,2) 2 (2,3) 1>> C(j,i)=full(T) % 將一般矩陣賦值給一稀疏矩陣，仍返回稀疏矩陣C = (1,1) 1 (3,1) 1 (4,1) 1 (6,1) 1 (3,2) 2 (6,2) 2 (2,3) 1 (5,3) 17.4 稀疏矩陣的交換與重新排序稀疏矩陣S的行交換與列交換可以用以下兩種方法表示。（1）對於交換矩陣P，對稀疏矩陣S進行行交換可表示為PS，進行列交換可表示為P

S’。（2）對於一個交換向量p，p為一般向量包含1到n個自然數的一個排列。對稀疏矩陣進行行交換，可以表示為S(p,:)。S(:,p)為列交換形式。對於矩陣S的某一列進行行交換，可以表示為S(p,n)，如S(p,1)為對第1列進行交換。【例2-23】 稀疏矩陣S的交換。>> p=[1 3 2 4];>> S=eye(4,4)S = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1>> P=S(p,:)P = 1 0 0 0 0 0 1 0 0 1 0 0 0 0 0 1>> V=S(p,2)V = 0 0 1 0矩陣P的第1行為S的第1行，第2行為S的第3行，等等。即對矩陣S的行，按照向量p指定的順序進行調整。>> S1=speye(4,4)S1 = (1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1>> P1=S1(p,:) % 對於稀疏矩陣行列的交換，返回的形式仍為稀疏矩陣P1 = (1,1) 1 (3,2) 1 (2,3) 1 (4,4) 1對於稀疏矩陣S1進行行列的交換，返回的P1仍為稀疏矩陣。對稀疏矩陣的列重新排序，有時可以使矩陣分解的速度更快，最簡單的矩陣排序是根據矩陣中非零元素的個數進行的，這種方法對於元素極不規則的矩陣很有效，特別適用於非零元素在行或列中數目變化較大的矩陣。MATLAB提供有一個非常簡單的函式colperm，可以實現這種排序方法。此函式的M檔案僅有以下幾行：function p=colperm(S)if ~ismatrix(S) % 判斷輸入變數是否是一個矩陣 error(message(‘MATLAB:colperm:invalidInput’)); % 不滿足條件的話返回錯誤資訊end[~,p] = sort(full(sum(S ~= 0, 1)));程式的第5行，實現了以下4個功能。（1）呼叫S ~= 0判斷矩陣中各元素是否為0，若不為0則返回邏輯值1。（2）函式sum求上一步建立的矩陣各列的和，也即為各列中非零元素的個數。（3）函式full將上一步建立的向量轉換為一般向量的格式。（4）使用函式sort對上一步操作建立的向量元素進行升序排序，函式sort的第2個輸出引數p，即為對矩陣S中各列中非零元素的個數進行重新排序的交換向量。【例2-24】 對下面的矩陣A，先用函式colperm獲取一個交換向量p，然後根據向量p對矩陣A的列，按照非零元素的個數升序排序。>> A=[0 1 2 3;3 2 1 0;0 0 2 0;1 0 0 2]A = 0 1 2 3 3 2 1 0 0 0 2 0 1 0 0 2>> p=colperm(A)p = 1 2 4 3>> B=A(:,p)B = 0 1 3 2 3 2 0 1 0 0 0 2 1 0 2 0結果顯示，矩陣B就是將A的各列按照非零元素的個數升序排序的結果。7.5 稀疏矩陣檢視MATLAB提供有spy函式，用於觀察稀疏矩陣非零元素的分佈檢視。本小節舉例來說明spy函式的用法。【例2-25】 稀疏矩陣檢視示例。本例採用spy函式繪製Buckminster Fuller網格球頂的60×60鄰接矩陣檢視。這個矩陣還可用來表示碳60模型和足球。>>B = bucky;>>spy(B)

得到的結果如圖2-2所示。圖中顯示了稀疏矩陣B的非零元素分佈檢視。