11/22/2006

10.8.3 spdiags稀疏對角指令

前文討論的對角元素處理指令已可隨意就斜角方向置入任意元素。實際上在解題時,尤其進行數值分析過程中,所形成之矩陣雖然很大,但百分之八、九十的空間大部份都是為零。對於資料處理方面是一種相當大的資源浪費。為解決是項問題,必須使用稀鬆矩陣(Sparse matrix)處理,以節省時間與空間。spdiags指令則是針對對角元素部份之存取運作,與前面討論之diag之功能可相輔相成。

稀疏對角處埋指令格式如下:


B = spdiags(A)
[B,d] = spdiags(A)
B = spdiags(A,d)
A = spdiags(B,d,A)
A = spdiags(B,d,m,n)

spdiags(A)指令以不同的方式處理三個矩陣,可以作為輸入與輸出。其組成係利用不同之輸入參數完成各項功能。其中矩陣A之大小可為mxn,只要陳述非零元素之位址及數值即可,B則為一輸出項,旨在描述對角線項之陣列,其大小為min(m,n) x p;d則表示矩陣A中非零元素之之對角序列值。

矩陣A中若大小為m x n,則應有m+n-1條對角線。其序列以d參數表示,其範圍為由-m+1至n-1。例如,A之大小為5x6,則應有10條對角線,此10條線分別以d=-4,-3,-2,-1,0,1,2,3,4,5表示。d=0時表示正好在最長的對角線上。其位置可以由下圖表示。


例:

>> e=ones(4,1);A=diag(1:5)+diag(e,1)+diag(e,-1)

A =

1 1 0 0 0
1 2 1 0 0
0 1 3 1 0
0 0 1 4 1
0 0 0 1 5

>> [B,d]=spdiags(A)

B =

1 1 0
1 2 1
1 3 1
1 4 1
0 5 1


d =

-1
0
1

由結果可知,d指示k值所在,而B值則是非零對角線上之值,其起點係以左邊項為基準,不足的部份補零。就本例而言第一行為 k=-1,故由基線往右下延伸,不足部分填零;第三行為 k=1,由左邊基線開始,有部份不足,須先補零。這些值可以利用下列指令取得B中之值列,例如:由矩陣A中取k=1之對角線項值:

>> C=spdiags(A,1)

C =

0
1
1
1
1

若要在現有之矩陣增加一對角線項,則可利用 spdiags(B,d,A) 這個指令函數,其中B為填入之值,必須以行的型式,所給的值若超過矩陣A之範圍,則多餘的項會被截斷。d則指示置於那一條對角線上:

>> AA=spdiags([3 4 5]',-3,A)

AA =

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

>> full(AA)

ans =

1 1 0 0 0
1 2 1 0 0
0 1 3 1 0
3 0 1 4 1
0 4 0 1 5

上述指令之應用是在現有之矩陣填入對角線資料,下面的指令spdiags(B,d,m,n)則是創造一個新的稀疏矩陣,其大小為m-by-n。 注意陣列B與d之長度應相符合:

>> B = repmat((1:5)',[1 6])

B =

1 1 1 1 1 1
2 2 2 2 2 2
3 3 3 3 3 3
4 4 4 4 4 4
5 5 5 5 5 5

>> d = [-4 -2 -1 0 3 4];
>> A = spdiags(B,d,5,5);
>> full(A)

ans =

1 0 0 4 5
1 2 0 0 5
1 2 3 0 0
0 2 3 4 0
1 0 3 4 5


由所得之A矩陣中可知,其大小受m x n之限制,而左下方之排例由B之高項而下,超過者將被截斷。右上方則截斷原需補足的部份,或其尾項由最後項填起。上述之說明實際上僅討論到m=n的情況,實際上當m與n不同時,有些情況因為截斷位置之不同,結果也不一樣。以下面為例:

>> B=reshape(1:15,5,3)

B =

1 6 11
2 7 12
3 8 13
4 9 14
5 10 15

>> d=[-2 0 2]

d =

-2 0 2

>> A=full(spdiags(B,d,5,5))

A =

6 0 13 0 0
0 7 0 14 0
1 0 8 0 15
0 2 0 9 0
0 0 3 0 10

當m=n=5時,其填入之順序與前例相同,左下依右邊基準填入,再行截斷;右上則依底部填起,上部截斷。

>> A=full(spdiags(B,d,5,4))

A =

6 0 13 0
0 7 0 14
1 0 8 0
0 2 0 9
0 0 3 0

當m>n時,其左下依右邊基準填入,再行截斷,情況與前例相同;但右上因為截斷的關係,可依k=0處之截斷位置(此處為9倒數第二行),則依底部填起,此時底部也要截斷一行,再依序往上,截斷超出者。

>> A=full(spdiags(B,d,4,5))

A =

6 0 11 0 0
0 7 0 12 0
3 0 8 0 13
0 4 0 9 0

當m<n時,填入之順序相反過來,其左下改依底部往上填入,但依k=0處之截斷位置(此處為9倒數第二行),則依底部填起,此時底部也要截斷一行,再依序往上,截斷超出者。右上邊則改由上面之基準填入,再依長度進行截斷。此種情況完全與前面之m=n及m>n之情況相反,須特別注意。

在使用spdiags(A),取出其對角線值置入B矩陣時,也有類似情況,當m<n時,其取出之順序也會相反,例如:

當m=n時:


Matrix A Matrix B

6 0 13 0 0 1 6 0
0 7 0 14 0 2 7 0
1 0 8 0 15 == spdiags => 3 8 13
0 2 0 9 0 0 9 14
0 0 3 0 10 0 10 15

當m>n時:

Matrix A Matrix B

6 0 13 0 1 6 0
0 7 0 14 2 7 0
1 0 8 0 == spdiags => 3 8 13
0 2 0 9 0 9 14
0 0 3 0

當m<n時:

Matrix A Matrix B

6 0 11 0 0 0 6 11
0 7 0 12 0 0 7 12
3 0 8 0 13 == spdiags => 3 8 13
0 4 0 9 0 4 9 0