matlab 求矢量夾角_MATLAB稀疏矩陣

7 稀疏矩陣

稀疏矩陣是一種特殊類型的矩陣，即矩陣中包括較多的零元素。對于稀疏矩陣的這種特性，在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.000000000000000

1 2 0.500000000000000

2 2 0.333333333333333

1 3 0.333333333333333

2 3 0.250000000000000

3 3 0.200000000000000

1 4 0.250000000000000

2 4 0.200000000000000

3 4 0.166666666666667

4 4 0.142857142857143

4 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) 1

7.4 稀疏矩陣的交換與重新排序

稀疏矩陣S的行交換與列交換可以用以下兩種方法表示。

（1）對于交換矩陣P，對稀疏矩陣S進行行交換可表示為P*S，進行列交換可表示為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 =

1

矩陣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的非零元素分布視圖。