Title
スパース行列の最適ピポッティングの諸手法とそれらの
比較研究
Author(s)
喜屋武, 盛基
Citation
琉球大学理工学部紀要. 工学篇 = Bulletin of Science &
Engineering Division, University of the Ryukyus.
Engineering(8): 105-111
Issue Date
1975-01-30
URL
http://hdl.handle.net/20.500.12000/26297
疏球大学理工学部紀要 (工学 籍) 105
スパース行列の最適ピポッティングの
諸手法とそれらの比較研究
喜 屋 武 盛 基 本
On Some Optimal p
i
v
o
t
i
n
g
Algorithms
and
t
h
e
i
r
Comparisons
Surnrnany
The coefficient matrix of a very largesystem of equations isgenerallyvery sparse， i.e.， non-zero elements arevery few， say， lessthan 10percent oftotal number of elements. Thus in computer assisted analysisofsuch a system， sparsity oriental techniques for solution and data storage are very important.
Recently
，
many papers have been published on the sparse matrix techniques. In this paper，
three algorithms of so cal1ed収minimumfil1-in" have been studiedand actually programmed for F ACOM 230・15computer. Comparisions are made from the results of actual computations.
1
.
は じ め に NxNの行列の非零要素の数がN21r.比較して非常 IC.小さしまたNが大きいとき，その行列をスパス行 列 (sparsematrix) と呼ぶ。 このような行列を係数 行列としてもつような大規模な線形システムを電子計 算機で解く場合に重要なことは，行列のもつスパ{ス 性を活用することであり，それにより計算機のメモリ の節約と計算時間の短縮がはかられ，通常の計算方法 では計算機の記憶容量や計算時間等の制約により，物 理的IC.不可能なものでも計算可能にすることができ る。 近年とのスパース性を活用したいわゆる収スパ{ス 技法"の研究が盛んに行なわれ多数の研究報告がなさ れている。(l)~(12) 本文では特iζ非対称なスパ{ス行列を LjU分解す る場合の最適ピボッティング1I頂序決定の三つのアル 法11)，を比較検討した。 スパース技法とは，非零要素だけを計算機の記憶 装置にたくわえて計算を行なう方法であり，従って L/U分解またはガウス消去により消去をかけて行く 過程でもともと零であった要素が非零 Ir.変る (fill.in が発生するという)ことを出来るだけ最小にするよう なピポッティング操作の順序づけが霊要になる。2
.
三手法の理論的背景 2.1 Heishi-Ghausiの deterministicな方法 行列Aが (r-1)の段階の Gauss-J
ordanの消 去(以下GJE) を完了したものとする。 行列A
(
←1)は ) 唱 E -e。 ，
u ， z z、 、 、 E E E B B E E -a a E E E . ， ， ， ， 、 、 B I l l i -E J ' h A F E E B B B R E E Z E -a ' a目 、
n u -噌 ム 司 ム 唱 A A U F E E -E a B a a -E・
E・
-E， 、=
) ' i r ( A ゴリズム， Hsieh とGhausiによる deterministic こ』で Ar= (alj (r-1)) ， i = 1，…， n j= r，…， n な方法12)，同じく Hsieh と Ghausiによる proba- ブ-)レ行列 Brをつぎのように定義する。bilisticな方法13)，および喜屋武，白Jflと尾崎Ir.よ Br= (bij (←1)) るパイパータイトグラフとの対応により導びいた方 ただし 受付:昭和49年 5月9臼 ホ琉球大学理工学部電気工学科 10 : alj(r-l)= 0 btl ={
l
1 : aljい-1)キO (2.2)106 寝屋武:スパ{ス行列の最適ピボッテイングの諸手法とそれらの比較研究 定 理
2-1
Pivot axy(←1)キ O をとり pivoting操作を i行 について行なった際 A(戸 1)から A(r)への i行に ついてのGJEによる変換で A(のに新しく発生する 非零要素 (fill.in) の数は次のようになる。 (r)， (r-1) . (r-1) ( n ;'yl)GlJ = r (~)I~'b
i
y
(2 -3) と 』 で r(x)1 (i = 1. 2.…
，
一対の数である。 n)はつぎのような (r-1) (r-1) a x j 宇O • a i j = 0 ( iキx) したがって A(r-1)から A(r)への G.J.Eによって できる新しい fill-inの総数は次式で示される。 (r) n (r-1) (r-l) (nxll i)ω= l:! rωi b iy j-l(2-4)
上の定理より A (r-l)から A (r)への fill-inの 数を最小にする最適なピポットはC
n劣
?
〕
GJzmin〔
nJ
3
〕
GJ となるようなものである。 (r-l) (2-5) (2 -4)式は ax y キO をピボットとして A (r-l) から A (r)への G.J.Eによって発生する fill-inの総和である。 同式から明らかに (W担一1) (r-1) 個の零が axll宇Oをピポットとしてr段 階 の 消 去を行なった際の行列A
(r)において発生する。 ただし Wyはベクトル Wr=Col(Wr • Wr+l.
…
Wy: Wn)の中の非零要素の数を表わす。 従って.A
(r)にあった非零の正味の変化，すなわち G.J.E によって A(r)11:生じた零と fi11-inの数と の差は次式で与えられる。 (Nx~))GJ
= (n匁
)GJー (Wy-1) (2 -6) G.J.Eの r段 階 の 消 去 に お け る 最 適 な pivotは a Ij(r-1)であり次式のようになる。 (Nii})GJ=min〔
n与
)J
GJ -max CWy-1) したがって， もし (n zla))GJ=(n z27i)GJ= min (n fj))GJ(2-7)
のとき，すなはち minimumfill-inが同数の要素が2個あるとき(tieであるとき). Wy2と Wy1 を比較 して大きい方を最適ピポットとして採用する乙とがだ 当であるζとを示している。
2
-2 H
eishi-Chausiの確率論的な方法 定 理2-2
行列A
(r)の非零要素がランダム i己分布していて (r-l) a x.y
=ドOがr
番目のピボットであると仮定し， Cr-l) そのとき ai.Y 宇 O でr番目の消去の際のi行に できる新しい非零の確率的な数は (r) ^ (r-1) nx'y
i= r(~-)(' = L ((1-k)qx) (2-8) ただし，;
S
T
i
)
は確率的な r8
7
i
)
でk=dE
m+lは行列
A (
r
)
の次数。 関数で例えば.L (2.8)= 2 定 理 2. 3 またLc.)は整数化 となる。 行例 A(r)の中l乙非望書要素がランダム11:分布して (r-l) いて，a
X1 Y キOをr番目のピボットとし，そのと きG.J.E. により発生する fil1-inの最も確率の高い 値は次式で与えられる。〔
A
J
3
〕
GJ=E〔
;(Jfl)〕
(wy-1) (2-9) (r-l) _ ••^
(r-l) 乙乙で E(r (x) r -')は r(x)i の期待値であ る。 定 理 2. 4 行例 A(r)の非零要素がランダムに分布していて Qr. Wr (行および列の非零要素の数)が知られてい るとする。 もし.a xy (r-l)キO.Wy= wmin. qx =qmin であるとすると axy (r-l) が fil1-inを最 小にする最適なピボットである。 もし axy (r-l)=
Oならば minii((qiー1)(Wj-l))が最適なピポッ トとなる。 2-3 喜屋武一白1
1
1
-
尾崎(1)による方法 乙の方法についての理論づけは論文(1)に詳細に記 述しであり，論文 (9)には概略が示しであるので本文 では省略する。3
.
三手法のアルゴリズム 3 -1 Hsieh-Ghausiの deterministicな方法 盟J:Ll行例 A (r-1)のすべての非零要素 axy琉球大学理工学部紀要(工学第) 107 のfi1l.inを計算し，その中で最小の fill.inを与える 要素の集合をKとしつぎへ移れ。 主1L主 IKI=1 なら K= (axy)をr番目の ピボット1<:選ぴ操作4へ移り. 1 K 1
>1
なら操作3 へ移れ。重
.
.
f
L
1
.
Kの各要素の属する列の非零要素の数が最 大のものをピポッ トl乙選ぴ操作4へ移れ。 主.fL._! r=nなら操作完了。 r<nなら r←
r+ 1として A(r)←
A(r-l)の変換をして操作4へ戻れ。 3 -2 Hsieh.Ghausiの確率的な方法 主 主 よ 行列A(r-l)の i行およびj列の非零要素 の数をそれぞれωi.qiとし， ωy=min ω. qx= m m qで axyキOであるなら axyをr番目のピ ボット IL選ぴ操作 3へ移りもし axy=0ならつぎへ 移れ。 主-w2 (qi四 1)(ωi -1)を計算し，その中で最 小値を与える axy宇Oをr番目のピボット1<:選ぴつ ぎへ移れ。 皇-w3 r=nなら操作完了。 r<nなら r←
r +1として A(r)←
A(r-l)の変換をして操作 1へ戻 れ。 3司 3 喜屋武ー白1
1
1
一尾崎の方法 主J1_よ行列 A(r-l)の中で fi1l.in の数が零であ る非零要素 axyがあるかどうかをしらぺあればaxy をr
番目のピボットとして操作5へ移り，なければつ ぎ〈す、む。主主
i
行列 A(r-l)の中でfil1.inの数が最小の 非 零 要 素 の 集 合 をK
l
とし1
K
l
1
=
1ならKl=
( axy)をピボット 1<:選ぴ操作 5へ移り.1
K
l
1
>1 ならばつぎへす』む。盟主.2.
集合K
l
I
乙属する各非零要素をピボットと して消去を行ないその結果得られる行列の中で最小の fil1.inを含む行列を得る非零要素の集合をK
2
として1
K
2
1
=
1ならK
2=
(a町〕をr
番目のピボットと し操作5へ移9
. 1
K2
1
>1ならばつぎへすh
む。 操 作 4 K2 I乙属する各非零要素l乙対してそれぞれ の属する行の要素のもつ fi1l.inの総和を計算し，そ の最大値を与える要素をr番目のピボッ トとしてつぎ へ移れ。盟主
i
r=nならば操作完了。 r<nなら r←
r+1 として A(r)←
A(r-l)の変換を行ない操作1 へ戻れ。4
.
計算結果 上記の三手法をプログラムし非零要素のパターンを 任意にきめた例題14個についてのそれぞれの手法によ る fi1l.inの発生総数を計算した所，表 1の結果を得t
:
。
1 ，-~" _... _'1非零/行列のサイズ￨ 例題番号￨行列のサイズ 1...1--r.J/ I .J:，，~"'::".
.
.
.
'
I
1 ' "' ，._， . . . ( % ) fil1.inの総数 HG-1 HG-2 KSO 一一一一一I 2 3 4 5 6 7 8 9 10 12 13 14 46x46 46x464
6
x
4
6
4
6
x
4
6
46x46 46x4646x46
4
6
x
4
6
4
6
x
4
6
46x46 46x4646x46
46x46
5Ox50 4. 6 18 4. 6 9 4.7 13 4.7 14 5.7 1 5. 7 23 5. 7 26 5. 7 34 5. 7 47 6.1 60 6.5 71 6. 6 96 6. 6 76 5. 3 49 (表 1) 注: HG-1 : Hsieh-Ghausi の deterministicな方法 HG-2: Hshieh-Ghausiの probabiJisticな方法 KSO :喜屋武ー臼JII・尾崎の方法 17 16 9 7 11 11 14 14 1 1 20 20 25 23 32 31 48 47 55 52 68 65 92 88 76 75 47 43108 喜昼武:スパース行列の最適ピポッテイングの諸手法とそれらの比較研究 5. む す び 先ず HsiehGhausi による deterministicな方法と Probabilisticな方法とでは前者が一般に fill.inの発 生が少ないがこれは実際の fill.inの数を計算した前 者に対して，後者は fill.inの最小値を確率的な方法 で得たためである。しかしながら計算法が単純である ため計算時間は後者の万が早いことが期待出来る。 Hsieh.Ghausiの deterministicな方法と喜屋武一白 )11ー尾崎の方法については，両方とも実際 lζ計算した fill-inの最小値をとる点では同じだが“tie"が生じた 場合の最適ピボットの選定基準が前者は proba bilistic な方法をとるのに対して後者は，“tieけの各要素を実 際l乙ピボットとして操作を行なった場合つぎの操作 で最小の fi1l-inを生じるものを選んだ。このように “次の手を読む"ととを可能にしたのは論文 (1)に示 しであるようにパイパータイトグラフに対応させて考 察した結果，“つぎの手"で生じる fill-inの 情 報 を 現段階の情報だけで得られるととがわかったためであ る。 表1から明らかなように，喜屋武一白)11-尾崎の万 法が fi1l圃inの発生の点ですぐれており，十分実用価 値があるものと思はれる。(付録l乙Hsieh-Ghausiの probabi1isticな方法による FORTRAN プログラム をあげる。なお KONEやその他必要なサブルーチン は文献(9)に述べてある(喜屋武一白
)
1
1
-
尾崎の方 法)で用いたものの一部を併用した。) おわりに本研究にあたりプログラムやデータ等の作 成および電算機による計算の実行等IL協力した電気工 学科学生宮里憲幸，玉盛譲の両君に謝意を表する。 参 考 文 献 (1) 喜屋武，白川，尾崎:“線形系解析における最適 pivoting順序とこれに付随するパイパータイト グラブ，電子通信学会回路とシステム理論研資 CT -72-68， (1973-01) (2) H. Y. Hsieh-M. S. Ghausi “On Optimal -pivoting Algorithm in Sparse Matrices，" 付 録HG-l
の手法によるプログラムIEEE， Trans on Circuit Theory， CT・19pp 93・
96(Jan. 1972) (3) H. Y. Hsieh M. S. Ghausi
:
“
A Probabilistic Approach to Optimal-pivoting and Prediction of fi1l-in for Random Sparse Matric巴s"，IEEE， Trans， on Circuit Theory， CT-19 pp 329・336 (July 1972) (4) R. D. Berry:“An Optimal-ordef!ng of Elec -tronic Circuit eguations for a Sparse Matrix Solutions"， IEEE， Trans. on Circuit Theory， CT・18，pp 40-50(Jan 1971) (5) 岸，植竹:“アドミツタンス行列の Sparsity1<: 関する一考察勺電子通信学会論文誌 Vo1.55 -A. No.5 pp 205-212 (1972-05) (6) 翁長，川端，白川，尾崎:“一般佑ガウス消去法 と疎大行列の消去系列最適佑"， 電子通信学会， 回路とシステム理論研資 CT-73-64，(1973-12) (7) 喜屋武，白川，尾崎"一般の非対称線形連立方 程式のスパース技法によるPivotingI!民序決定"， 電子通信学会，回路とシステム理論研資 CT-73-65. (1973-12) (8) 相原，平山:“スパース回路網万程式の最適消去 系列の後退順序づけ勺 電子通信学会，回路とシ ステム理論研資 CT-73-62. (1973-12) (9) 喜屋武，白川，尾崎:“線形連立方程式のスパー シティとこれに付随するパイパータイトグラフH， 昭和48年度電子通信学会全国大会 制喜屋武，臼川，尾崎:“スパース行列のガウス消 去における最適ピボッティングアルゴリズムとそ のフォトランプログラム二 琉球大学理工学部紀 要工学篇 No.7， (1974-3) ω F . G. Gustavson， W. Linger， R.A. Wi1lou -ghly" Symbolic Gereration of an Opimal Crout Algorithm for Spavse Systems of Linear Equations勺 J. ACM， 17， p. 87 (1970)ω
後藤:“スパース処理技法を用いた最短径路問題 の解法"，電子通信学会，回路とシステム研資 CT 73-22 (1973-07) C 判 * * * * 材 料 * * * 材 料 材 料 * * * * * 紳 材 料 * * * * * * * * * * * * * * * * 材 料 材 料 材 料 榊 材 料 材 料 材 料 材 料 材 料 紳 材 料 材 料 榊 SUBROUTINE GHAUSI 1 (NZ， N， KG， K1ST， KCNT)琉球大学理工学部紀要(工学第)
COMMON I1ND (50). I1NZ (250). JJND (50). JJNZ (250). NNUBF (2日)
COMMON I]FL (100). IORDR (50). JORDR (50) DIMENSION KGHAU (50). K1ST (50)
DIMENSION 1L (50). JC (印).K1L (50). KJC (50)
C 紳 *(RECOMPUTE NUBF ACCORDING TO GHAUS1 METHOD)料 *
109
C 榊 *(CALCULA TlON OF NUMBER OF NONZERO AT THE EACH ROW AND COLUMN)帥 *
DO 100 1=1. N 1F (1-1) 1
∞
.
102. 101 102 CONTlNUE 1L (1) =1ND (1) JC (1) =JND (1) GO TO 100 101 1L (1) =1ND (1) -1ND (1-1) JC (1) =JND (1)ー
JND(1-1) 100 CONT1NUEC *料(F1NDINGM1N1MUM NUMBER OF NONZERO 1N ROW)
*
*
*
M1N=50 K1=O DO 200 1=1. N 1F (1L (1)) 201. 201. 202 202 CONTlNUE 1F (M1N -1L (1)) 201. 203. 204 203 Kl=Kl+1 K1L (Kl)=1 GO TO 201 204 M1N = 1L (1) Kl=1 K1L (Kl)=1 201 CONT1NUE 気淘 CONT1NUE 1COUNT=Kl
C 柿 *(F1NDING M1N1MUM NUMBER OF NONZERO 1N COLUMN)柿 $
M1N=50 K1=O DO 12
∞
J=1. N 1F (JC (J)) 1201. 1201. 1202 1202 CONT1NUE 1F (M1N -J C (J)) 1201. 1203. 1204 1203 Kl=K1 +1 KJC(K1)=J GO TO 1201 1204 M1N = JC (J) Kl=1 KJC (K1)=J f110 著書屋武:スパース行列の最適ピポツテイングの諸手法とそれらの比較研究
1201 CONTINUE 12'
∞
CONTINUE JCOUNT=K1C *本*(FINDlNG EXISTENCE OF NONZERO) **取 DO 3
∞
1=1
.
ICOUNT IIL=KIL (1) DO 302 J=
1
.
JCOUNT JJC=KJC (J) IF (IIL-1) 307. 307. 308 ?JJ7 IAl=l GO TO初9 308 IAl=IND (IIL-1)+1 初9 IB=IND (IIL) DO 305 K=IA，l IB IF (INZ (K) -JJC) 305. 310. 305 ?JJ5 CONTINUE ?JJ2 CONTINUE ?JJO CONTINUEC 柿 傘(CALCULATIONOF SEKI AND STORAGE)
輔
.
1=1 IA=1 GO TO 20 1 IA=IND (1-1)+1 20 lB=IND (1) IF (IB-IA) 2
1
.
22. 22 22 CONTINUE DO 90 K=IA. IB J=INZ (K) DO 10 L=IA. lB IF (J-INZ (L)) 10. 2. 10 2 NUBF (L)= (IL(I)ー1)市(JC(J)-1) 10 CONTINUE 90 CONTINUE 21 CONTINUE IF (N-I) 5. 5.，4 4 ， 1=1+1 GO TO 1 5 CONTINUEC 帥 *(FINDlNG MINIMUM SEKI)帥 * NZN=NZ
CALL KONE (NZN. K1ST. KCNT) IF (KCNT-1) 555. 555.550
RETURN 550 DO 14 K=，l KCNT KK=K1ST (K) 琉球大学理工学部紀要(工学篇) CALL RETV (N. KK. 11.JJ) IF (JJ-1) 11. 11. 12 11]A=O GO TO 13 12 ]A=]ND (JJ -1) 13 ]B=]ND (JJ) JDIF=JB-]A KGHAU (K) =JDlF 14 CONTINUE
C 柿 *(FINDlNG MAXMUM NONZERO AT COLUMN IN MIN SEKI)紳 常 MAX=KGHAU (1) MG=1 KCN
