11/22/2006

10.8.2 稀疏矩陣之應用

稀疏矩陣是精簡一般含有零元素甚多的矩陣,僅就非零元素之位置加以陳述。此項作業可用sparse指令達成。其語法及相關指令如下:


S = sparse(A)
S = sparse(i,j,s,m,n,nzmax)
S = sparse(i,j,s,m,n)
S = sparse(i,j,s)
S = sparse(m,n)
TF = issparse(S)
R = spones(S)


sparse 指令很簡單,它可以將一個含有許多零的元素矩陣改寫為位址表示,但其實質仍如原來之矩陣,功能相同,例如:

>> M=eye(5)

M =

1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1

>> B=sparse(M)

B =

(1,1) 1
(2,2) 1
(3,3) 1
(4,4) 1
(5,5) 1

產生之新矩陣B內僅記載具有非零元素之位置及數值。雖然B具上述型式,但其大小仍然與原來之M矩陣相同,而且可以直接與其他矩陣相互運算,有如與M矩陣運算一般。例如上述之稀疏矩陣B,其大小仍然與原先的一樣:

>> size(B)

ans =

5 5

>> B*magic(5)

ans =

17 24 1 8 15
23 5 7 14 16
4 6 13 20 22
10 12 19 21 3
11 18 25 2 9

若要恢後原來的型式,則可使用full指令復原,如:

>> full(B)

ans =

1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1


除由原有之矩陣轉換外,sparse亦可由參數直接產生,其參數包括:

  L = sparse(i,j,s,m,n,nzmax)

其中,i,j為非零元素之方位向量,而s為其對應之值向量,最終可以產生一個m x n大小之稀疏矩陣,亦即:

  L(i(k),j(k)) = s(k)


這個指令主要用於某些相當大的矩陣,若全部依照傳統規矩輸入可能會使電腦炸掉!但利用這種方式直接轉為稀疏矩陣,則即使在小電腦中也可以從容執行。其中,nzmax為非零元素所能處理之最大數,不特別指定時,其值為length(s);而 i, j與 s之向量長度必須相同,s中若有零值項時會被忽略,連帶對應之i與j也被刪除。若i,j指向重複位置時,相同之值將被相加。

由於這個指令之輸入參數很多,實際上可以配find這個指令,找出某參數之非零元素,因為其輸入之三個參數正好為sparse之前三項輸入,例如:

M =

2 0 9
1 0 0
0 4 0

>> [i,j,s]=find(M);
>> S=sparse(i,j,s,max(i),max(j))

S =

(1,1) 2
(2,1) 1
(3,2) 4
(1,3) 9

>> MM=full(S)

MM =

2 0 9
1 0 0
0 4 0


採用稀疏矩陣主要目的在避免在運算過程中,使用超大矩陣以致記憶體空間不足,例如:
M = sparse(1:10,1:10,1) 產生一個對角元素為1之稀疏矩陣,其結果與M=sparse(eye(10))是相同的,問題是後者一定要透過實際的矩陣,已失去稀疏矩陣的原意。只是執行時,電腦將耗去較多時間,只是以時間換空間,應用上也是一個較好的選擇。

稀疏函數spfun


對於稀疏矩陣之非零元素與函數關係之處理,則可使用spfun,其語法如下:

f = spfun(@fun,S)

例如,設S為前述例中之稀疏矩陣,求其對S之非零元素作正弦函數:

>> f=spfun(@sin,S)

f =

(1,1) 0.9093
(2,1) 0.8415
(3,2) -0.7568
(1,3) 0.4121

此處之fun可用函數之握把,其運算僅針對非零元素,時間上也可因而縮短。至於如何辨別S是否為稀疏矩陣,可以使用issparse(S)進行測試。此外另一個函數 spones(S)則可產生一結構與S相同之矩陣,但S中之非零元素部份均取代為1。例如:

>> spones(S)

ans =

(1,1) 1
(2,1) 1
(3,2) 1
(1,3) 1

矩陣之相關邏輯運算亦可在稀疏矩陣中應用。

處理稀疏矩陣之亂數指令sprandn/sprand


稀疏矩陣之均勻分佈及常態分佈可用sprand及sprandn指令。其語法如下:

R = sprandn(S)|sprand(S)
R = sprandn(m,n,density)/sprand(m,n,density)
R = sprandn(m,n,density,rc)|sprand(m,n,density,rc)


亂數指令 sprandn(S)是產生與S結構相同的亂數,其值在0至1之間。sprandn(m,n,density) 則是產生m-x-n個稀疏矩陣之亂數,其分佈以density*m*n 為基準,其值域為0 <= density <= 1)。其中,rc為反覆條件(reciprocal condition)。rc可為以 lr <= min(m,n)之向量。亂數若以對角線作成對稱的分佈,則可使用sprandsym相關指令為之,其格式如下:

R = sprandsym(S)
R = sprandsym(n,density)
R = sprandsym(n,density,rc)
R = sprandsym(n,density,rc,kind)


轉換用指令spconvert



由於稀疏矩陣之格式與實際存檔資料之格式不同,故自檔案中讀入之資料必須經過spconvert指令處理:

S = spconvert(D)


上述過程中,可先在ASCII資料檔中編輯對應之[i,j,v] 或 [i,j,re,im] 欄位,然後轉為稀疏矩陣。資料需為每列三欄,若為四欄則可產生複數之矩陣。例如下面為 ASCII檔案,其名稱為 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

執行下列指令後,即可產生H之稀疏矩陣:

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