連立 1 次方程式の基礎知識〜 および Gauss の消去法の安定性について 〜
渡部 善隆 y学会や研究集会での発表や、そこでの議論を聞いていると、 (時間が限られていることもあって)ある「知識」を持っていることが前提となって話が進んでいく場合がよくあります。特に数値計算の分野では、「××法」と名のつく方法が山のように登場します。その方法の「意味」がわかっている
なら問題なく話について行けるのですが、何を言っているのかイメージすら掴めないことが往々にしてあります。この場合、そこからの議論に取り残されてしまい、 \ぼやっ〜" と聞いているうちに完全に興味を失ってしまい、「セッションが終ったら何を食べにいこうかなぁ」などと違うことを考えるようになります。それはよろしくありません。
この記事は、科学技術計算の主役の一つである連立 1 次方程式の数値解法に登場する手法および用語の解説です。この分野の用語の特徴は、 \略語だらけ" で、 \よく知らない人の名前が多い" ことです。解説は最低限の意味の説明にとどめ、数学的な証明や具体的なアルゴリズムはすべて参考文献に押しつけています。さらに詳しく知りたい方は、記事の終りに文献紹介および「優れた文献を紹介している文献を紹介」しますので、そちらを御覧下さい。
なお、最後の章で Gauss の消去法の安定性に関する興味ある (ある面で恐ろしい) 話題を紹介します。こちらも「知識」として持っていると、例えば座長をしていて質問に窮したときや、懇親会で知らない人に囲まれたときなどに役立つかも知れません。
この記事によって、多少ともみなさんの \ぼやっ〜" とする時間が減ることになれば幸いです。

目次
1 Introduction

299
299

2.1.14 直交行列 : : : : : : : : : : :

299

1.1

線形作用素 : : : : : : : : : : : : : : 293

2.2

解法の大まかな分類: : : : : : : : :

299

1.2

近似モデル : : : : : : : : : : : : : : 293

2.3

直接法

: : : : : : : : : : : : : : : :

300

1.3

非線形を線形にするからくり

294

2.4

反復法

: : : : : : : : : : : : : : : :300

Newton-Raphson 法と連立 1 次方程

2.5

共役勾配法 : : : : : : : : : : : : : : 300式の関係 : : : : : : : : : : : : : : : 295

2.6

各方法の組合せ

1.4

2

293

2.1.12 M 行列 : : : : : : : : : : : :2.1.13 置換行列 : : : : : : : : : : :: : : :

1.5行列の次元 : : : : : : : : : : : : : : 295

1.6

精度の問題 : : : : : : : : : : : : : : 296

連立 1 次方程式
2.1

3

297

: : : : : : : : : : :300

直接法の概要301

3.1

Cramer の公式 : : : : : : : : : : : :

301

3.2

Hamilton-Cayley の定理に基礎をおく方法

行列の種類 : : : : : : : : : : : : : : 297
3.3: : : : : : : : : : : : : : : :301

三角化に基礎をおく方法 : : : : : : : 301

2.1.1

密行列 : : : : : : : : : : : : 297

2.1.2

疎行列 : : : : : : : : : : : : 297

2.1.3

帯行列 : : : : : : : : : : : : 297

3.3.2

Crout 法 : : : : : : : : : : :303

Doolittle 法: : : : : : : : :303: : : : : : : : :

303

3.3.1

Gauss の消去法

: : : : : : :

302

2.1.4

三重対角行列

297

3.3.3

2.1.5

対称行列 : : : : : : : : : : : 298

3.3.4

Cholesky 法

2.1.6

対角行列 : : : : : : : : : : : 298

3.3.5

修正 Cholesky 法

2.1.7

下三角行列 : : : : : : : : : : 298

3.3.6

Bunch の方法

2.1.8

上三角行列 : : : : : : : : : : 298

3.3.7

越境法 : : : : : : : : : : : : 304

2.1.9

Hessenberg 行列 : : : : : : :298

2.1.10 正定値行列 : : : : : : : : : :298

2.1.11 対角優位行列298: : : : : : : :: : : : : : : :3.4

y 九州大学大型計算機センター・研究開発部

291: : : : : :303: : : : : : : :304

対角化に基礎をおく方法 : : : : : : : 305
3.4.1

Gauss-Jordan の消去法 : : :

305

3.4.2

逆行列の積表現

305

3.4.3

合同変換による方法 : : : : : 305: : : : : : :

3.5

3.6

3.7

4

直交ベクトル法

3.5.2

直交行列法 : : : : : : : : : : 305: : : : : : :: : : : : : : : : : : : : : : :3.6.1

分割計算法 : : : : : : : : : : 306

3.6.2

縁どり法 : : : : : : : : : : : 306

3.6.3

エスカレータ法

3.6.4

階数零化法 : : : : : : : : : : 306

306

三重対角行列専用の方法 : : : : : : : 306
3.7.1

Gauss の消去法, Cholesky 法

306

3.7.2

cyclic reduction 法

306: : : : :疎行列, 帯行列専用の方法 : : : : : : 306

3.9

現在人気の解法: : : : : : : : : : :4.1

反復法の創始者: : : : : : : : : : :

309

4.2

Jacobi 法 : : : : : : : : : : : : : : :

309

4.3

Richardson 法 : : : : : : : : : : : :

310

4.4

Gauss-Seidel 法

310

4.5

SOR 法

: : : : : : : : : : :: : : : : : : : : : : : : : :

310

4.5.1

SLOR 法 : : : : : : : : : : :

310

4.5.2

SSOR 法 : : : : : : : : : : :

311

4.5.3

SBOR 法: : : : : : : : : :

311
311

4.7

ADI 法 : : : : : : : : : : : : : : : :Aitken の  2 法 : : : : : : : : : : :

311

4.8

Chebyshev 加速法 : : : : : : : : : :

311

4.9

モンテカルロ法: : : : : : : : : : :312

4.10 反復改良法 : : : : : : : : : : : : : :

312

5.1

鮮烈なデビューとその後 : : : : : : : 313

5.2

基本的な考え方

5.3

計算の手順 : : : : : : : : : : : : : : 314

5.4

最急降下法 : : : : : : : : : : : : : : 3145.5共役方向法 : : : : : : : : : : : : : : 3145.6共役勾配法 : : : : : : : : : : : : : : 315

5.7

5.8

5.9

: : : : : : : : : : :ICCG 法 : : : : : : : : : : :

316

5.7.2

MICCG 法 : : : : : : : : : :

316

5.8.1

ILUCR 法5.8.2

MILUCR 法 : : : : : : : : :: : : : : : : : : :: : : : : : : : : : : : :5.10 双対共役勾配法: : : : : : : : : : :5.10.1 ILUBCG 法

6.4

ピボット選択 : : : : : : : : : : : : : 323

6.5

ピボット選択の利点

6.6

ピボット選択が必要ない行列

6.7

スケーリング : : : : : : : : : : : : : 324: : : : : : : : :: : : :

323
324

6.7.1

部分ピボット選択の場合 : : : 324

6.7.2

完全ピボット選択の場合 : : : 325

プログラム例 : : : : : : : : : : : : : 325

327

7.1

実行時間の調査: : : : : : : : : : :

327

7.2

行列の値の交換: : : : : : : : : : :

328

7.2.1

親切なコンパイラ

7.2.2

最適化指示行の挿入 : : : : : 328

7.2.3

本当に行と列の値を交換する: : : : : :

328
329

7.3

添字の動き : : : : : : : : : : : : : : 329

7.4

GECP プログラム : : : : : : : : : :

330

7.5

実行時間の比較332: : : : : : : : : : :

333

8.1

消去法の歴史 : : : : : : : : : : : : : 3338.2

後退誤差解析 : : : : : : : : : : : : : 334
334

8.2.1

マシンイプシロン

8.2.2

条件数 : : : : : : : : : : : : 334

8.2.3

誤差評価 : : : : : : : : : : : 335

8.2.4

悪条件の場合の対策 : : : : : 335

8.2.5

Hilbert 行列 : : : : : : : : :: : : : : :

335

8.3

不安定因子の上限: : : : : : : : : :

336

8.4

Wilkinson さんの例題 : : : : : : : :

336

8.5

なぜ部分ピボット選択法が使われる

8.6

random matrix での実験

8.7

病理学的な例 : : : : : : : : : : : : : 340: : : : : : : : : : : : : : : :: : : : : :

338
339

8.7.1

Wright さんの例 : : : : : : :

340

8.7.2

Foster さんの例

341

: : : : : : :

317
317

5.10.2 MILUBCG 法 : : : : : : : :317: : : : : : : : : : : :

322

317
317

: : : : : : : : : : : : : : :: : : : : : : : : :

321

316

: : : : : : : : :5.12 BiCGStab 法

LU 分解との関係のか？

共役残差法 : : : : : : : : : : : : : : 316

5.11 QMR 法

消去過程 : : : : : : : : : : : : : : : 322

6.3

8 Gauss の消去法は安定か？

前処理付き共役勾配法 : : : : : : : : 315

GMRES 法6.2

: : : : : : : : : : :7 GECP のスピードアップ作戦

313

5.7.1

321

Gauss の消去法

313

共役勾配法の概要

319

6.1

6.8

307

309

反復法の概要: : : : : : :6 Gauss の消去法のアルゴリズム

306

分割法

: : : : : : :5.15 INSPEC を用いた実験

305

3.8

4.6

5

直交化法 : : : : : : : : : : : : : : : 305
3.5.1

illustration :

317
317

5.13 共役勾配法の復活: : : : : : : : : :318

5.14 反復解法の汎用性: : : : : : : : : :

319
292

Hidaki, Naoko

1 Introduction

0
BBB
B@

序章では、連立 1 次方程式が数値計算に登場する
背景について考えてみます。
なお、厳密に話を進めると定義がたくさん必要になり、読む人を退屈させますので、かなり荒っぽい説明になることをご勘弁願います。
また、以下特に断わらない限り、小文字は実数ま
たは実数のベクトル、大文字は実数の正方行列を表

x1
x2

.
.
.

n

x

1
CCC
CA

A
0!

0
BBB
B@

y1
y2

.
.
.

n

y

1
CCC
CA

線形作用素が逆作用素をもつとき、この作用素はします。

正則(non-singular) であるといいます。行列 A が正則であるとは、

1.1

01

線形作用素

A

連立 1 次方程式1 (linear equations) は、数値計算= I (単位行列)となる逆行列 (inverse matrix) A01 がただひとつ存を試みる人の前に頻繁に立ち現れる問題です。

X

01

A = AA

在することです。

方程式は連立 1 次方程式 (2) は、逆行列の値がきちんとわ

n

a

j =1

ij xj = bi ;i = 1; 1 1 1 ; n

かるならば、行列とベクトルの積

(1)

01

x = A

で表記されます。未知数は n 個の xj (j = 1; 1 1 1 ; n)b(3)で求めることができます。
ただし、実際の数値計算で逆行列 A01 を露骨にです。
(1) を行列形式で書けば計算にいくことはそんなに多くありません。 A01 が

Ax = b

きちんとわからなくても、 (2) をみたす x をもとめ

(2)

る効率のよい方法がたくさんあるからです。
それらの方法に比べ、逆行列を計算し、積 A01 bです。連立 1 次方程式の数値解法 (numerical soluequations)とは、行列 A とベクトの計算で x を求める方法は、多くの手間がかかり、ル b に対し、 (2) をみたす x を数値的に求める方法

精度もあまりよくない上、普通は記憶領域を大量にです。

消費します ([41] 参照)。

tion

of

linear

なお、この記事は全て実数の世界で考えますが、そのため、逆行列の値がすべて必要になるなどの基本的には複素数の世界にも拡張できます2 。

特別なことがない限り、 (3) の計算は行なわないの「1 次」とは「線形3 」のことです。ある作用が線が \常識" になっています。

形(linear) であるとは、それが比を保ち、さらに和連立 1 次方程式の数値解法は、次章以降で概観すに対する作用の結果がそれぞれに対する作用の和にることにします。なお \線形性" についてのイメー

一致することを意味します4 。

ジが欲しい方は [64] がお勧めです。

行列 A は、 n 次元実ベクトル空間 IRn から同じ

く n 次元実ベクトル空間 IRn への線形な作用素とし1.2て見ることができます。
1 「連立一次方程式」と書くこともあります。

近似モデル

微積分の概念は、自然現象や社会現象を数学の言

2 四則演算を実数の定義から複素数の定義に切り替えればほと葉でモデル化することを簡単にしました。

んどそのまま通用します。これは、複素数 z が平面上の点 (x; y )微分方程式や積分方程式のような無限次元5 の関にちょうど対応し、無理なく計算手順が拡張できるからです。
3 「線型」とも書きます。要は意味が通ればいいのです。

5 有限個の基底で表現できないという意味です。連続系ともい

4 記号で書くと T (v + v ) = T v + T v , T (av ) = aT (v )
1
2
1
2

います。

293

数方程式を数値的に解く場合は、数学モデルを有限次元の近似モデルに置き換える6 ことが必要です。

F (u) = 0

情報処理技術の急速な発達は、有限要素法・境界をみたす u を求めるという問題に書くことができま

要素法・差分法などの近似モデル化の研究を大いにす。 F は設定した空間上の非線形作用素です。

助けました。

設定した関数空間の要素 u(0) で F が微分可能な連立 1 次方程式は、この近似モデルを具体的な数とき、 F (u) を u(0) のまわりで Taylor 展開11 し、 3値にする段階でしばしば登場します。そして、他の項目以降を切り捨てて近似します12 。

計算7 に比べ、計算時間がたくさん必要なことと、(0)

F (u)  F (uちゃんと方程式が解けていないと後でたいへん困るという意味で、数値解析の中心部を占める大切な解F法のひとつです。

0) + F (u(0)(0))(u 0 u )(6)

0 (u(0) ) は F の u(0) での微分として決まる \線形"な作用素で、逆写像 F 0 (u(0) )01 を持つと仮定しま

もう少し詳しく見ていきましょう。

す13 。
ここで、 Taylor 展開で近似した (6) の右辺を 0 に1.3するような u(1) を考えます。つまり、線形方程式

非線形を線形にするからくりF (u立てた数学モデルが線形になることはあまりありません。たいていは非線形(nonlinear)な方程式にな

F

例えば、流体関係の学会などに顔を出すと、こん@t+ (v 1 r) v = 0grad p +div v = 0

u

1
Re

0 (0)
(1)
(0)
) + F (u )(u
0 u
)=0

(7)

0 (u(0) )01 を (7) の両辺に作用させ、変形するとな微分方程式をよく見かけます。
@v(0)を満足する u(1) を見つけることにします。

ります。

8>
><
>>:

(5)

(1)

(0)
0 (0) 01(0)=u
+ F (u ) F (u )で u(1) が求まります。この u(1) は、最初に適当に1v決めた u(0) よりも真の解 u に近いことが期待される

(4)

ものです。
そこで、今度は u(1) のまわりで同じように Tay-(v と p の境界条件)

lor 展開し、つぎの u(2) を(4) は、 Navier-Stokes 方程式8 とよばれる有名な方

程式です (詳しい説明はパスします。気になる方は

u[46] をお読み下さい)。

(2)

(1)
0 (1) 01(1)=u
+ F (u ) F (u )で求めます。

(4) などの非線形な方程式を直接解くことは、かなりの困難をともないます9 。そこで数値解析の常道

Newton-Raphson 法とは、このように適当な初期値 u(0) から出発し、 k = 1; 2; 1 1 1 について次々に近として、非線形方程式の代わりにある線形方程式の似解の列計算をくりかえすことで、近似解の収束列を作り出すということをよくやります。

u

その代表的な方法である Newton-Raphson 法

(k)

(k01)
0 (k01) 01(k01)=u+ F (u) F (u)(8)を作っていく方法です。

を例にとって説明します。

すると、いつかは u(k) と u(k01) が見わけがつか

Newton-Raphson 法

なくなるくらいに近くなるかも知れません。その時、 u(k) = u(k01) としてしまい、繰り返しの式 (8)うまい具合に関数空間10 と作用素 F を定義するこに放り込むととで、どんな非線形方程式も

(k )

u

6 この段階を「離散化」といいます。

(k )0 (k ) 01(k )=u+ F (u ) F (u )7 行列やベクトルを作ったり、収束性を調べたりメッシュを

11 イギリスの数学者 Brook Taylor さん (1685{1731) に由来

切ったりする作業。
8 Louis Marie Henri Navier さん (1785{1836) と George

します。
12 微分の定義や Taylor 展開の話、また、なぜ近似になるのか

Gabriel Stokes さん (1819{1903) の名前に由来します。

などの説明は省略します。要は大学の最初に習う 1 変数の実関数

9 解析的な解を求めるときも、数値的な近似解を求めるときもf (x) の延長だと思って下さい。

13 なぜ \線形" かというと、訳は簡単で、線形になるように決めです。
10 よくみかける空間を記号だけ紹介すると、 IRn (自分から自分ているからです。関数解析学では F 0 (u(0) ) を線形化作用素 (lin-への写像だと思って下さい), Sh , Lp , L1 , W m;p , H m , H0m ,H 01 , C m , C01 などなどです。

earizedoperator)

ともいいます。実際に定式化する場合、 F の微分可能性や関数空間の設定を注意深く説明しないと、関数空間論の人から怒られます。

294

つまり
F

0

(u

次元の関数空間の議論は、 [40], [44] などをお読み下(k) 01(k )) F (u ) = 0

さい)。

となります。この両辺に F 0 (u(k) ) を引っかけて
(k )

F (u

しかし、無限次元での定式化の場合は、数値計算のために、何らかの方法で無限次元を有限次元に\射影する" ことが不可避の作業ですので、計算機の)=0中では結局有限次元の問題、さらに最終的には IRnが得られます。つまり u(k) は F (u(k) ) = 0 の解に限での線形計算に落ち着きます17 。

nりなく近いに違いない、という作戦です。

(8) が有限次元空間、特に IR での反復だとすれNewton-Raphson 法は \F (u) はごく狭い範囲でば、線形作用素は行列の形で書くことができます。

は 1 次式で近似できる" という考え方をもとにした従って、方法で、出発点がもし真の解に十分近ければ、極めて急速に収束することが知られています。

u

この節の最初に登場した Navier-Stokes 方程式 (4)

(k )

(k01)
0 (k01) 01(k01)=u+ F (u) F (u)のアンダーラインの部分は、「ベクトル F (u(k01) )も、も ちろん (8) の 形に 記述す るこ とがで きま すに行列 F 0 (u(k01) ) の逆行列を左から引っかける」こ(具体的な方法は [16] を御覧下さい)。

とを意味しています。

Newton-Raphson 法は、簡単に Newton 法と呼ば

よって反復式 (8) は、連立 1 次方程式れることもあり、修正 Newton 法、準 Newton 法、
区間 Newton 法など、いろいろな変形版も加え、広F

く非線形方程式の数値解法に使われています14 。

0

(u

(k01)

(k 01)(k 01))^u= F (u)をu^(k01) について解き、この u^(k 01) を使って、次Newton-Raphson 法をはじめとする非線形方程式の u (k ) をの反復解法、近似解の精度的保証などの数学的な考
察は [47] をお読み下さい。

(k )

u

なお、 Newton-Raphson 法は、有名な Sir Isaac
Newton さん (1642{1727) と Joseph Raphson さん(k 01)
(k01)
=u
+u
^

で求めることと同じです。

(1648?{1715?) に由来します。

そのようなわけで、非線形方程式の数値解法では連立 1 次方程式を解く必要がしばしば生じます18 。

1.4 Newton-Raphson 法と連立
1 次方程式の関係

1.5

行列の次元

有限次元で例えば変数が 1 つのここで、本稿の主役の連立 1 次方程式が NewtonRaphson 法の何処に顔を出すのか探してみましょ3f (x) = xう。

先ほど考えていた非線形作用素 F は、特に IRn

2+x +x+1= 0

を Newton-Raphson 法で解く場合、反復列は上の作用素でなくとも構いません。空間が IRn の場

合、 F 0 (u(k) ) は F の Jacobi 行列となります。

(k )

x

Sobolev 空間15 などの無限次元の関数空間の場合(k 01)
0 (k01) 01(k 01)=x+ f (x) f (x)echet 微分16 とすることで、は、 F 0 (u(k) ) を F の Frです。このとき、変数は 1 個なので、 f 0 (xk01 ) はそのまま Newton-Raphson 法が適用できます (無限1 2 1 の行列です。よって、逆行列はいたって簡単で

14 昔の数値計算の教科書には「Newton-Raphson 法」と書か

0

f (x

れていましたが、現在は「Newton 法」とだけ書く本が増えました。しかし、 [50] によれば、 Newton さんは確かに方法は示し(k 01) 01
)
=

1

f

0 (x(k01) )で求まります。

ていたが、一般的な公式にまとめたのは Raphson さんだそうです。ここでは、 2 年程前に開催されたシンポジウムで、ある数値

17 つまり、すでに有限次元の非線形問題になっているものを計算の大御所が「Raphson 氏の偉大な業績を忘れて、皆 New-

Newton-Raphson 法で解くか、無限次元で定式化した Newton-

ton 法と呼ぶのはけしからん！ Newton-Raphson 法と呼べ！」

Raphson 法を離散化して解くかの違いです。手法は数値解に対しと吠えていたので、 Newton-Raphson と書いておきます。
15 Sergei L'vovich Sobolev さん (1908{?) は、 (当時) ソビエてどのレベルの誤差評価を期待するかで異なります。
18 この大きな流れを押さえていないばっかりに、修士論文の発トの数学者です。
16 Rene Maurice Fr
echet さん (1878{1973) は、フランスの問に絶句してしまい、著しく心証を悪くした人を知っています。

数学者です。

ご注意下さい。

表会で「結局それは何を計算するんですか？」という (温かい) 質295

演算を脳天気に信用するのはちょっと危険です21 。

一方、多次元の連立非線形方程式

少なくとも、学会や論文発表では「なんとなく実

f (x) = 0

験結果と似ているから大丈夫でしょう」という必殺Tn)T
f = (f1 f2 1 1 1 fn )x = (x1 x2 1 1 1 xの切札に加えて、「連立 1 次方程式を含む数値計算
が丸め誤差の混入によって致命的な精度損失を (多の場合19 、 f 0 (x(k01) ) は @ fi =@ xj を (i; j ) 成分とす分) 起こしていないであろう」ことを強弁するだけの根拠22 が必要です。なぜなら、る n 2 n の Jacobi 行列になります (詳しくは [60])。
もともとの問題が無限次元の場合、方程式の次元

もし、連立 1 次方程式を含む数値計算が丸め誤差によって致命的な精度損失を起こしていた場合、モデルを設定し、シミュレーションした意
味がなくなってしまうは離散化する近似空間の基底の数に一致します。差分法や有限要素法などを用いる場合、未知数は領域にばらまいた node の数に対応しますので、この総数が連立 1 次方程式の次元になります20 。

からです。

1.6

連立 1 次方程式の解の信用性については、最後の精度の問題

章で列変換をともなう Gauss の消去法の安定性の話
題を切り口にして、もっと掘り下げて考えたいと思

さて、いざ連立 1 次方程式を解き、数値解を求めいます。

る段階にくると、「誤差」の問題が浮上します。つまり、計算者はどれくらいの精度を数値結果に要求するかという問題です。
数学モデルとしての非線形方程式の解の存在や分岐現象などの挙動に興味のある人でしたら、収束性や誤差解析も含めて、数値解に高い精度を要求する
はずです。
一方、もともと自然現象や社会現象を数学モデル
化し、最終的な数値解をシミュレーションとして利用したいと思う人の立場から見れば、既に方程式を立てた段階である程度の誤差が入りこんでいるので、連立 1 次方程式をはじめとする数値計算を完璧な精度で求める必要はない
し、また意味もない。
と思うでしょう。
微分方程式などのモデルから出発したとしても、計算にかかるコストの問題から、どうせ離散化の段階で誤差が入るので、連立 1次方程式の解法に高い精度を追求して計算時間を無駄にするより、多少の誤差の混入には目をつむって、高速なプログラムを組もう。
と思う人もいるでしょう。
しかしながら、自分のプログラミング技術はとも線形作用素 A =かく、数値計算の普通の演算方法である浮動小数点
19 記号 \T " は、転置を意味します。

1=2

01=2

1=2

1=2

!
によって定義される 2 次元ベクトル場 (05  x1 ; x2  5)

20 テレビのヒーローもので、 4 次元や 5 次元の世界からやってくるのがありますが、まだまだ甘いといえるでしょう。スーパー

21 興味のある方は [62] を御覧下さい。

コンピュータでは数万次元の世界を日常的に扱っています (もっと22 例えば、問題自体の安定性、行列の性格の良さ、数値計算ラも、設定している空間の意味が全然違いますが…)。ただし、 2 次

イブラリの信用性、拡張精度との比較、誤差評価、優秀な学生が

元世界 (鏡の国だそうです) からやってきた正義の味方「ミラーマプログラムしたから大丈夫、この私が言うのだから間違いない！

ン」はすごいと思います。

etc.

296

2

連立 1 次方程式

0
BBB
BBB
BBB
@

前章では、数値計算の中で連立 1 次方程式が出現
してくる手順を簡単な例に沿って見てきました。
連立 1 次方程式の係数行列 A は、 n 次元ベクトル空間 IRn 上の線形作用素であり、 A の形は非線形問題の線形化の過程で決定されます。そして、モデルや離散化の方法によって行列 A の形は異なります。
連立 1 次方程式の数値解を求める段階になると、

疎行列

* *
*
*
*
*
*

*

*

*

*

*

*

*

2.1.3 帯行列

行列 A の形にあわせた様々な解法が登場します。これまでに提案され、現在も誰かが利用している手法

*

*

*
*

*

1
CCC
CCC
CCC
A

帯行列 (band matrix) は、疎行列であって、かつだけでも 100 種類以上あるのは確実でしょう。
この章では、線形作用素である行列 A の性格を分

0 でない成分が行列の対角線上と対角線の比較的近類した後、連立 1 次方程式の数値解法を概観したい

くに分布している行列です。これもよく見かけまと思います。

す。

0
BBB
BBB
BBB
@

個別の解法の説明は 3, 4, 5 章で行います。

2.1

行列の種類

連立 1 次方程式の解法は、行列 A の形に大きく依存します。そのための基礎知識として、行列用語を簡単に説明します。

2.1.1 密行列

帯行列

** **
* * *
0
** **
*** *
* *** **** *** *** *** *
*
**
**
*
*** *
** ** *** ***
* ***
0
* ***
*
**
*
*

1
CCC
CCC
CCC
A

2.1.4 三重対角行列

密行列 (dense matrix) は、とくに規則性や特徴のない一般の行列のことです。「密」とは、 0 以外の三重対角行列 (tridiagonal matrix) は、帯行列の要素が行列の中に充満していることを意味します。

特別な場合で、 ji 0 j j > 1 の成分がすべて 0 な行列

0
BBB
BBB
BBB
@

密行列

*********************************************************************************************************************************************************************************************************************************

1
CCC
CCC
CCC
A

です。つまり、 0 でない成分が行列の対角線上と、対角線の前後だけに分布している行列です。三項行列ともいいます。
A

値問題があっという間に解けます。

0
BBB
BBB
BBB
@

2.1.2 疎行列
疎行列 (sparse

matrix)が三重対角行列の場合、連立 1 次方程式や固有は、成分のほとんどが 0の行列をいいます。偏微分方程式を離散化して連立
1 次方程式に帰着させる場合、だいたい疎行列にな

ります。
297

三重対角行列

**
***
***
***
0
***
***
***
***
***
***
***
0
***
***
***
**

1
CCC
CCC
CCC
A

0
BBB
BBB
BBB
@

2.1.5 対称行列
対称行列 (symmetric matrix) は、 aij = aji すな

わち A = AT となる行列です。

行列が対称だと、解法の速度、メモリー、安定性などいろいろいいことがあります。

2.1.6 対角行列は、 \diagonal" の頭文字をとって \D" と書きます。
対角行列は大変素晴らしい行列です。なぜなら、

Hessenberg 行列23 は、 i > j + 1 ならば aij = 0となる行列です24 。

つまり Hessenberg 行列は、上三角行列に逆行列は各成分の逆数をとって並べ直したものです
し、固有値にいたっては最初から対角成分に並んで

0
BBB
BBB
BBB
@

*

対角行列

*
*

*

*
*

0
*
*

0

*
*

*

*
*

2.1.7 下三角行列

*
*

nn01

a21 ; a32 ; a43 ; 1 1 1 ; a

1
CCC
CCC
CCC
A

の部分を加えたものです。

0
BBB
BBB
BBB
@

Hessenberg 行列***********************************************************************************************************
0*************
*****
****
***
**

1
CCC
CCC
CCC
A

2.1.10 正定値行列

下三角行列 (lower triangular matrix) は、対角成正定値行列 (positive de nite matrix ) は、固有値分より上の部分

が全部プラスの値をとるような行列です。いい換え

ij j i < j gfa

ると、任意のベクトル x 6= 0 に対し、 2 次形式

がすべて 0 の行列です。

x

以下、下三角行列を表す場合は、 \lower" の頭文
字をとって \L" と書きます。

0
BBB
BBB
BBB
@

1
CCC
CCC
CCC
A

2.1.9 Hessenberg 行列

対角行列 (diagonal matrix) は、対角成分以外がすべて 0 の行列です。以下、対角行列を表す場合いるからです。

上三角行列***************************************************************************************************
0******
*****
****
***
**
*

下三角行列

*
**
***
****
*****
******
0***************************************************************************************************T Ax > 0

を満足する行列です。正値行列、正定符号行列とも

1
CCC
CCC
CCC
A

いいます。
A

が正定値で、しかも対称行列だったなら、素晴らしくいいことがおきます25 。

2.1.11 対角優位行列

X

対角優位行列 (diagonal dominant matrix) は、
ja

2.1.8 上三角行列ii j >

i6=j

ij j;

ja

i = 1; 1 1 1 ; n

をみたす行列です。つまり、対角成分が (絶対値の上三角行列 (upper triangular matrix) は、対角成意味で) 他の成分を圧倒して \でかい" 行列のことで分より下の部分す。
優対角行列ともいいます。

ij j i > j g

fa

23 K. Hessenberg さん (1874{1925) に由来します。

がすべて 0 の行列です。

24 上

Hessenberg

行列 という場合もあります。
25 ただし、後ほど出てくる Hilbert 行列のように、対称かつ正以下、上三角行列を表す場合は、 \upper" の頭文
字をとって \U " と書きます。

定値であっても、どうしようもなく性格が悪い行列もあります。

298

2.1.12 M 行列

0
BBB
BBB
BB@

次に、 P を A の右からかけると

M 行列 (M-matrix) は、正則かつすべての i 6= j
に対し aij  0 であり、かつ逆行列の成分すべてが0 以上の行列のことです。 M 行列は、物理の世界にAP

よく出てきます。対称な M 行列は正定値行列になります26 。

2.1.13 置換行列

=

a5

a2

b5

b2

c5

c2

d5

d2

e5

e2

f5

f2

a1

a6

a4

a3

b1

b6

b4

b3

c1

c6

c4

c3

d1

d6

d4

d3

e1

e6

e4

e3

f1

f6

f4

f3

1
CCC
CCC
CCA

となります。今度は、行はそのままで列が入れ換

置換行列 (permutation matrix) は、どの行にも、わっています。順番に、行列 A の《1 列と 5 列》《5またどの列にも、ちょうど 1 個ずつ 1 があり、他は列と 4 列》《4 列と 6 列》《6 列と 3 列》を入れ換

すべて 0 となる行列です。

えると AP になります28 。 2 列目はそのままです。

以下、置換行列を表す場合は、 \permutation" の置換行列は、行列の分解の際に大変重要な役割を

頭文字をとって \P " と書きます。
A果たしますが、 \単位行列もどき" の形をしている

置換行列 P を A の左からかけるということは、

ため、プログラムの中で配列として登場することはの行の順序を並べかえることを意味し、右からかなく、あくまでも裏方に徹します。

けるということは、 A の列の順序を並べかえることを意味します。

0
BBB
BBB
BB@

具体的な例をあげます。行列 A をA =置換行列 P を

P

=

0
BBB
BBB
BB@

a1

a2

a3

a4

a5

a6

b1

b2

b3

b4

b5

b6

c1

c2

c3

c4

c5

c6

d1

d2

d3

d4

d5

d6

e1

e2

e3

e4

e5

e6

f1

f2

f3

f4

f5

f6

0
BBB
BBB
BB@

0

0

1

0

0

0

0

1

0

0

0

0

0

0

0

0

0

1

0

0

0

0

1

0

1

0

0

0

0

0

0

0

0

1

0

0

1
CCC
CCC
CCA

直交行列 (orthogonal matrix) は、転置行列が逆行列となるような行列のことです。つまり
01

A

c1

c2

c3

c4

c5

c6

b1

b2

b3

b4

b5

b6

f1

f2

f3

f4

f5

f6

e1

e2

e3

e4

e5

e6

a1

a2

a3

a4

a5

a6

d1

d2

d3

d4

d5

d6

=A

T

となる行列です。
以下、直交行列を表す場合は \Q" と書きます。

(\直交" の意味は、 A の第 i 行と第 j 列の内積が

1
CCC
CCC
CCA

とします。すると P A は

PA =

2.1.14 直交行列



ij =

1

i = j

0

i 6= jとなることからきています。
直交行列を複素数に拡張した行列がユニタリ行列(unitary matrix) です。

直交行列は固有値問題の解法で主役をはります
が、連立 1 次方程式では単なる脇役の一人に過ぎま

1
CCC
CCC
CCA

せん。

2.2

解法の大まかな分類

連立 1 次方程式の数値解法は、次のように大別で
きます。
直接法
反復法

となります。行列 A の《1 行と 3 行》《3 行と 6 行》

連立 1 次方程式の数値解法《6 行と 4 行》《4 行と 5 行》を順番に入れ換える

共役勾配法と P A になります27 。 2 行目はそのままです。

各方法の組合せ

26 対称であることを M 行列の定義に加えている本もあります。
27 巡回置換 (cycle) の記号で書くと (13645) です。

以下、各解法の特徴をあげます。
28 こちらは (15463) です。

299

また、なかなか収束しないような \たち" の悪い

具体的な手法は、章を改めて概観します。

行列もあります。そのため、実際に計算をやってみ

2.3

ないと演算回数はわかりません。

直接法

2.5

直接法 (direct method) は、もしも丸めの誤差がなければ29 有限回の演算で厳密な解が得られる解法

共役勾配法

共役勾配法 (conjugate gradient method ) は、直のことです。

接法と反復法の長所をあわせ持った手法として、最直接法に必要な演算回数は、行や列の入れ換え操近とみに注目されています。

作が多少入るので不確定なところがありますが、計原理的には、丸めの誤差がなければ高々 n 回の反算する前からわかります。

復で厳密解が得られます。そのため、直接法の一種

そのため、行列の性質や精度に関係なく、ほぼ同といえます。

じ計算時間で答が出てきます。
A, b 0!

有限回の演算

0

0! x

x

1

! x

! 111 ! xn=x

行列が密行列であったり、都合のいい性質を持つ

が、実際は丸め誤差のためそうはいかず、欲しい精

かどうかわからない得体の知れないものであって度になるまで反復を続けます。従って、反復解法にも、計算が失敗しない限り、とりあえず解を出して分類されることもあります。
共役勾配法は、行列の形を整える \前処理"をほどくれます。そのため、直接法は標準的な解法とされこすことにより、格段の収束率の向上が得られ、特ています。
直接法の弱点としては、行列が疎行列であってもに正定値対称な疎行列の場合に威力を発揮します。

\0 が多い" という利点を生かしにくく、計算機の記また、非対称行列にも拡張が行われています。

憶容量が反復法に比べたくさん必要なことです。ただし、この困難を克服するためのいろいろな工夫も2.6

提案されています。

各方法の組合せ

また、演算回数は普通次数 n の 3 乗に比例しま

浮動小数点演算において発生する丸め誤差を考えす。そのため、次数が多くなればなるほど、計算機

れば、直接法を用いたとしても厳密解が得られるわの性能をうまく引き出したプログラムを書いておかけではありません。

ないと、 \有限回で終了する" といいいながらいつ出てくる数値結果からできるだけ誤差をなくすたまで経っても計算が終らない状況になります。

めに、直接法で得られた結果を反復法で改良する方法もあります。

2.4

また、反復法や共役勾配法の過程で、方程式の一反復法

部を直接法を用いて解く方法などもあります。

反復法 (iterative method) は、適当に選んだ初期

その他、どんな方法でもいいですが、得られた近値から出発して、真の解に収束していく近似解の列

似解をもとに、ある集合の中に厳密な解が存在するを逐次作成していく手法です。

ことを誤差限界とともに証明する精度保証付き数値計算法 (numerical methods with guaranteed accu-

反復の各段階では、それまでに得られた近似解のracy )

情報を駆使し、より良い近似解を生成する工夫が施

metic) と不動点定理をもとにした大変興味ある数値

されています。
x0もあります。これは区間解析 (interval arith-

! x

1

2

! x

計算法です。
! 111 ! x

最近、有理数演算を含むいくつかの精度保証付
きソフトウェアが利用できる環境になりましたの反復法の計算は、行列が 0 でない部分を格納する

で30 、機会があれば広報で紹介したいと思います。

記憶領域があれば可能です。反復法は、理屈の上でとりあえずの文献として [25], [65] をあげておき

は無限解の反復をしない限り厳密解を得ることはでます。

きません。そのため、適当なところで反復を打ち切る必要があります。
29 つまり、すべてを有理数演算で処理するか、無限桁の精度を

30 著作者の許可が得られるなら、センターでも公開したいと考持つ夢のような計算機で浮動小数点演算が実行できたならばえています。是非使ってみたいとお思いの方は、要望をrequest@cc.kyushu-u.ac.jp まで、どしどしお寄せ下さい。
300

3

直接法の概要3 章から 5 章までを使って、前章で大まかに分類

そのため、数値計算の手段としてはほとんど役にした連立 1 次方程式の数値解法をさらに細かく見て立ちません。しかし、線形代数のプログラミングのいきましょう。まずは直接法です。

題材にはいいかも知れません。

分類は Westlake さんのハンドブック [34] を参考

なお、 [48] に Cramer の公式の Fortran プログラに、最近の数値計算の本によく載っているものだけ

ムが載っています。

を紹介します31
これから紹介する方法は、プログラミングの演習

3.2

問題としては面白いかもしれませんが、実際の研究の手助けとするには心細い方法もあります。

Hamilton-Cayley の定理に基礎をおく方法

William Rowan Hamilton さん (1805{65) と Arthur

3.1 Cramer の公式

Cayley さん (1821{95) の名前がついたこの定理も、

これまた線形代数の世界ではよく知られています33 。

名前は大変よく知られていますが、数値計算で使われることは滅多にありません。 Gabliel Cramer

Hamilton-Cayley の定理

さん (1704{52) が 1750 年に発表した次の定理によっ正方行列 A の固有多項式 (characteristic poly-ています32 。

nomial) F () = jI 0 Aj をとれば, F (A) = 0

Cramer の公式
Aとなる.

が正則であれば, Ax = b は唯一の解を持ち,

解は
x

i

jA j

i=

jAj

この定理を用いて、 A01 を

i = 1; 1 1 1 ; n

01

A

ただし Ai は A の第 i 列を b でで与えられる.

X
n

=

i=1

i

p A

i01

の形で求めることができます。なお pi は Ai の対角置き換えたもの, j 1 j は行列式を表す.

要素の和を用いて順番に求めることができる実数です。

むしろ Cramer の公式といえば、上の公式より導!0

かれる 2 2 2 の逆行列の公式

a b
c d

1

= ad 01 bc 0dc 0ab

この方法も、 A のベキを記憶する必要があること!

と演算回数の問題などから実用的とは言えません。
が、線形代数プログラミングの題材にはいいかも知れません。詳しくは [3] を参照下さい。

の方が有名でしょう。

3.3

しかし、 Cramer の公式を数値計算で使おうとす

三角化に基礎をおく方法

ると、 n + 1 個の n 次行列式を計算せねばなりませ

直接法で現在もっともよく使われるのが、三角化ん。これは、他の解法と比べて反則のように手間がに基礎をおく方法です。プログラムが容易に書けるかかる上に、高精度の演算が要求されます。

上、一般の行列に広く適用でき、速度と精度もなか

31 連立 1 次方程式の分類および文献としては、 Forsythe さん

なかのものです。

が 1953 年に発表した有名な論文 [9], [10] があります。その中でこの手法のもとになる定理は、行列の分解に関す

紹介してある解法だけでも既に 400 を超えています。
32 Cramer さんはスイスの数学者です。 Cramer の公式は発表

る次の定理です。
33 Cayley-Hamilton の定理ともいいます。

の 20 年後、 Pierre Simon Laplace さん (1749{1827) の紹介で世に広まりました ([17])。

301

LDU 定理
A

を y について解きます。ここで L は下三角行列なので、前進代入 (forward substitution)35 と呼ばれるが正則ならば, 適当な置換行列 P を用いて簡単な計算方法ですぐに y が求まります。次にこの

P A = LDUyを右辺としてUx = yの形に分解できる. また, L, U の対角要素が 1のとき, 分解は一意に定まる.

を解きます。 U は上三角行列なので、今度は後退代
入 (backward substitution)と呼ばれる計算方法で x

この分解を LDU 分解 (LDU decomposition) とが求まります。

呼びます。 \三角化" とは、行列 A を L や U などの

本当に前進代入と後退代入でスラスラと解が求ま

三角行列の積に分解することをいいます。

るのか疑問の方は [45] を御覧下さい。これらは 6 章

3 2 3 行列で簡単な例を作ってみましょう。 P は!

!

単位行列となるので省略しています。
A

2 2 2
2 3 3
4 5 6

L

=

1 0 0
1 1 0
2 1 1

D

2 0 0
0 1 0
0 0 1

!

U

1 1 1
0 1 1
0 0 1

でもう一度登場します。

!

さて、この x は何ものかといいますと、先ほど解いた 2 つの方程式から y を消去することでLU x = P b

このように、任意の正則な行列は、適当な行の入の解であることがわかります。ところで、 A は P A =れ換え34 によって、下三角行列と対角行列と上三角

LUと分解してましたので、 x は行列の積の形に分解できます。
P Ax = P b

また、対角行列 D と上三角行列 U の積を改めてUだと考えれば、直ちに次のことがわかります (証の解です。 P は転置行列で正則ですので、結局 x は明は [49])。

Ax = b

LU 定理
A

基本的な LU (LDU ) 分解の方法は Gauss の消去が正則ならば, 適当な置換行列 P を用いて法ですが、行列の性格に適応させたり、スピードを追求するなどの目的に合わせてたくさんの変形が考P A = LU

案されています。

3.3.1 Gauss の消去法の形に分解できる. また, L の対角要素が 1 のとき, 分解は一意に定まる.

本来の意味での Gauss の消去法 (Gaussian elim-

これが有名な LU 分解(LU decomposition) で

ination)

す。

は、行列 A を成分の左端から結果が上三角行列になるように順番に消していく過程で、 b の処

LU

分解の形は、先ほどの LDU 分解で DU が一!

緒になったものです。
A

2 2 2
2 3 3
4 5 6
P

の解であることがわかります。

L

=

1 0 0
1 1 0
2 1 1

!

U

2 2 2
0 1 1
0 0 1

理もあわせて行ないます。つまり、 LU 分解の過程!

で同時に前進代入もやってしまいます。
この手順は、 \行列の LU 分解が完了したあとでbについて前進代入をする" 形に簡単に変更できます36 。
従って、最近 Gauss の消去法は「LU 分解法」とは行の順序を入れ換えるだけの置換行列ですの呼ばれることも多くなりました ([54])。

で、「適当に行を入れ換えた行列 A を考える」こと

Gauss の消去法は、 LU 分解の手法によって、まで、正則行列 A はた行列の種類により、さまざまな変種があります。
A = LU

内積形式に分解できると考えることもできます。

Gauss の消去法

内積形式 (inner-product form)では、分解処理を

行列を LU 分解して、それが何になるかと言え行列の左の列から右の列に、さらに列の上から下にば、この形からスラスラと解が求まるのです。

順次行ないます。

まず、連立 1 次方程式

35 前進消去(forward elimination) ともいいます。

Ly = P b

36 この方が、 A はそのまま b を変えて何度も計算する場合に便34 置換行列 P を左から引っかけることと同じです。

利です。

302

Crout 法との本質的な違いは、 Crout 法が U のデータがすべて主記憶に入りきらない大規模な
LU

分解で威力を発揮します ([49])。

対角成分を 1 とするのに対して、 Doolittle 法は L

Gauss の消去法の対角成分を 1 とするところだけです ([12], [34])。

外積形式

最近は Crout 法の方が有名になってしまい、影がうすくなっています。

外積形式 (outer-product form)は、 k = 1; 1 1 1 ; n

3.3.4 Cholesky 法に対し、 k 列の消去にひき続き、 k 行の分解を行ないます。数値計算の本に書いてあるアルゴリズムはCholesky 法38 は、次の分解定理をもとにしていだいたいこの方法です。

ます (証明は [33] 参照)。

この方法を応用すると、ベクトル計算機向きのプ

Cholesky 分解定理

ログラムが組めることが知られています。例えば、
行列をいくつかのブロックに分割した上で、ブロッA

クごとまとめて分解していく方法37 などです。

が対称かつ正定値ならば, 下三角行列 L が存在し

外積形式 Gauss 消去法のアルゴリズムとプログラ

A = LLムの注意点は、 6 章で詳しく考えます。

Tの形に分解できる. さらに, 分解は対角要素を正になるように定めることで一意に定まる.

縁どり法
Aが対称行列のときに使われることがあります。

この分解定理は、 [20] によれば F. R. Cholesky
さんが 1916 年に発表した方法です39 。

分解は外積形式と反対方向に行ないます。
行と列をどこかに追加した新しい行列に対して
LU

Cholesky 法は、この分解定理を利用し、 Gauss

分解する時に便利です ([59])。

対称の消去法を経由せずに直接 L を積和で計算する手法

Gauss 法です。分解ができてしまえば、あとは Gauss の消去法と同様に Ly = b と LT x = y を解くことで解が得られます。

対称 Gauss 法 (symmetric Gaussian method)は、Cholesky 法の計算は、分解の過程で対角成分の行列が対称な場合、その性質を利用し計算の手間

平方根を計算することから、平方根法(square rootと記憶容量を半分にするように工夫したものですmethod)([50])。

ともいいます。分解は L を求めるだけなので、演算回数が一般の Gauss の消去法の半分で済み

3.3.2 Crout 法ます。
Cholesky 分解は、 A が正定値でなくても (適当Crout 法は、 Prescot D. Crout さんが 1941 年

な転置行列を左からひっかけることで) 可能ですが、に発表した Gauss の消去法の変形版です。

L

手順は、 Gauss の消去法に出てくる中間結果をの対角要素が複素数になる場合があります。また、平方根の計算によって、精度が悪くなる場合が記録することを省略し、積計算を順に足し合わせてあります ([51])。

いくことで LU 分解を行います。順序は異なります

3.3.5 修正 Cholesky 法

が、計算は Gauss の消去法と同じになります。
Crout 法は、計算機の性能をうまく引き出すよう

修正 Cholesky 法 (modi ed Cholesky method)は、にプログラムすると Gauss の消去法以上のスピード

別名がいろいろあり、変形 Cholesky 法、改訂 Choleskyが出る場合があり、 LU 分解の有力な方法として幅38 本によっては \Choleski" と書いてあります。

広く使われています ([45], [49], [59])。

39 [72] には \Andre-Louis Cholesky" と書いてありました。

3.3.3 Doolittle 法

どちらが正しい のかはよくわかりませんでした。 Cholesky さんの連立 1 次方程式に関する寄与を紹介した Benoit さんの論文

Doolittle 法は、 M. H. Doolittle さんが 1878 年[1](1924 年) によれば、 \Commandant" Cholesky はフランスに発表した論文から名前がつきました。 Black 法と陸軍地誌部に属していて、第一次世界大戦 (大戦争と書いてありまもいいます。 Doolittle 法も Crout 法と同じようにした) 中の 1918 年に戦死されたそうです。また、 Mathematical
Tables and other Aids to Computation Vol.5, pp.112 (1951)

積和の形で LU 分解します。

を見ると、どうやら Rainsford[24] さんという人による Benoit

37 2 段同時消去法 ([45])、ブロッキング法などがあり、 SSL IIさんの論文 (仏語) の英語訳も存在するみたいです。

の VALU はブロッキング法を用いています ([70])。

303

3.3.6 Bunch の方法

法、対称 Cholesky 法、 square-root-free Cholesky法ともいいます。

修正 Cholesky 法のところで、行列 A が対称なら\みそ" は、 Cholesky 分解で生じる平方根の計算

ば原理的には LDLT 分解できるが、精度が悪くなを回避するところにあります。分解は次の定理が保る場合があることを述べました。

証します。

Bunch の方法は、 J. R. Bunch さんと L. Kauf-

修正 Cholesky 分解A

man さんが 1976 年に発表した手法 ([4], [5]) で、数が対称かつ正定値ならば,値的に安定した行列の分解を与える方法です。
分解は、次の定理に基づいています。

T

A = LDL

M DMの形に分解できる. さらに, 分解は対角行列 Dが対称ならば, 適当な置換行列 P を用いてAの要素を正にすることで一意に定まる.
P AP

修正 Cholesky 分解は、分解の形から LDLT 分解

T 分解

T = LDLTの形に分解できる. ただし, L は対角成分がすとも呼ばれます。修正 Cholesky 分解も、 Gauss のべて 1 の下三角行列, D は高々次数 2 の対称な消去法を経由しません。

ブロックから構成される対称ブロック対角行列.

なお、上の定理は、 A が対称であれば (適当な置
換行列を左から引っかけることで) 成立するのです

このような分解は数値的に安定となります。

が、対称正定値行列でないと分解の過程で非常に精Bunch さんの方法は、ブロック対角ピボッティン

度が悪くなることがあります。そこで、普通は A がグ手法 、あるいは \M DM T 分解"ともいいます。

対称かつ正定値の場合のみに使用されます。

分解の後は、順番に Ly1 =分解ができてしまえば、あとは

T
L y

3

P b,

Dy2

=

y1 ,

= y , x = P T y を計算すれば解が求まりま

2

3

す。

Ly = bSSL II, NUMPAC は Bunch さんの方法を用いたとサブルーチンをサポートしています ([69], [71])。

Tx = yDL

3.3.7 越境法を解くことで簡単に解が得られます。
修正 Cholesky 法も、 Gauss の消去法のように分越境法 (below-the-line method)は、方程式系 Ax =解の手順でバリエーションが幾つかあります ([59])。

外積形式修正b

Cholesky 法

を拡張して計算を行う方法です。
Aを拡張した行列 M をM =分解列より後の列を常に更新し、計算済みの列は影響を受けません。 A が帯行列、疎行列の場合によく使われます。

内積形式修正

A

K

H

0

!

で定義します。この後は 2 つの方法に分かれます。

Cholesky 法

Crout の越境法

外積形式と反対に、計算済みの行データを使ってM分解列を計算します。小規模行列に用いられます。

の n 行 n 列まで Crout 法を適用して A を LDU

分解する方法です。すると、 H と K の形をうまく決めることで、右辺ベクトル b を変化して作った行
列 B に対して A01 B が求まります。

縁どり法が、はっきり言って記憶領域の無駄使いです。

左の行から順番に分解行を求める手法です。 A がAitken の越境法

帯行列、疎行列の場合に使われます。

スカイライン(skyline) 法

こちらは拡張行列 M の A と H の行に Gauss の消去法を適用します。 H と K の形により、拡張行

envelope 法ともいいます。縁どり法の変形で、疎列 M の右下に結果として A01 や A01 B を作り出す

行列の場合に用います。非対称行列にも拡張がなさのが目的です。

れています。
304

名前は Alexander Craig Aitken さん (1895{1967)

まくいきません ([34]) し、連立 1 次方程式の解が欲に由来します。興味のある方は [34] を参照下さい。

しい場合ならば他の方法が有利です。

3.4

3.5

対角化に基礎をおく方法

直交化法\任意の行列は行と列に有限回の基本操作をほど直交化法 (orthoginal method) は LU 分解法などこすことによって、同値な対角行列に変形できる"に比べて演算回数がかなりかかるため、連立 1 次方という線形代数の基本的な定理にもとづいた方法で程式の数値解法で実際に用いられることはそんなにす。

ありません。

3.4.1 Gauss-Jordan の消去法

定性が信用できないときに、精度の比較で引っ張り

ただし、行列の性質が悪い場合や、他の解法の安出されることがあります。

Gauss の消去法に似た方法ですが、対角要素の上

0
BBB
B@

10
CCC BBB
CA B@の成分も消去することで、最終的に
^
a11

0
a
^22

..

.
^
ann

0

x1
x2

.
.
.

n

x

1 0 1
CCC BBB CCC
CA B@ CA=

3.5.1 直交ベクトル法

^
b1
^
b2
.
.
.
^
b

英語で書くと \method of orthogonal vectors" です。 A が対称なとき、 A に直交するベクトルを順次作ることでn

LALの形まで持っていきます。

の形に分解します。 L は下三角な直交行列です。こ

結果として、 Gauss の消去法のように後退代入をれよりやる必要がありません ([51])。

A

Gauss-Jordan40 の消去法は、掃き出し法 (sweepout

T =D

01

T D01 L=L

が簡単にわかり、従って x が計算できます。この方ともいわれ、昔は盛んに使われていま

method)法は一般の行列にも拡張できます ([34])。

したが、演算回数が Gauss の消去法の 1:5 倍もかか
ところが、行消去の演算が 1 行目から n 行目まで

3.5.2 直交行列法「せーの」で、しかも独立に行われるという性格か

直交行列法 (methodることから、次第に使われなくなりました41 。

ら、並列計算機向きの解法ではないかと指摘され、

tion)

最近にわかに注目されつつあります ([60])。

解から導かれます。

ofmatrixorthogonaliza-は、固有値の代表的な計算手順である QR 分

QR 分解

3.4.2 逆行列の積表現

Aが正則ならば, 直交行列 Q と上三角行列 R が存在して積行列法ともいわれます。 Gauss-Jordan の消去

A = QR法に関連した手法で、 \これとこれをかけあわせれば A01 になる" という乗数表を作り、精度と記憶場の形に分解できる. また, 分解は R の対角成分

所を節約しようという方法です ([34])。

をすべて正になるように定めれば一意である.

3.4.3 合同変換による方法

証明は [51] を御覧下さい。 QR 分解で方程式はAが対称行列ならば、ある直交行列 R によって
D = R

QRx = bT ARとなるので、まずと対角化できる、という性質を用いて A01 を生成す

Qy = bる方法です。しかし、一般に A が正定値でないとうを解きます。 Q は直交行列なので転置すれば逆行列

40 \Jordan" は、有名な数学者 Camille Jordan さん (1838{になり、1922) でなく、ドイツの測地学者 Wilhelm Jordan さん (1842{

y = Q1899) の名前からきています。従ってドイツ語読みをすれば「ヨルダン」が正しいそうです ([72])。
41 最近の数値計算の本で紹介されないことすらあります。このTb

で求まります。次に R は上三角行列なので、Rx = y世界では速いことが絶対の正義なのです。

305

3.6.4 階数零化法

を後退代入でスラスラ解けば終りです。
QR

分解を行なう主な算法は

階数零化法 (rank annihilation method) は、 H.

42 修正 Gram-Schmidt 直交化法
Householder 変換による方法43Givens 変換による方法44

S. Wilf さんが 1959 年に発表した方法で、逆行列が既にわかっている B とベクトル u, v を用いて A =B + uv

01

Aの 3 つです。これらは、どちらかといえば固有値計

01

0(B 01 u)(v T B 01 )1 + v T B 01 u

るかもしれませんが、連立 1 次方程式の解だけが欲れ以上の深入りは避けます。詳しくは [61] を御覧下しい場合には、計算時間の関係からあまりお勧めでさい。また、 Householder さんの簡潔にまとまったきません。

原論文 [19] もおすすめです。

3.7

分割法

三重対角行列専用の方法

偏微分方程式を差分法で解いたりすると、しばし

分割法 (partitioning method) はブロック分解法とば三重対角行列が現れます。 A が三重対角行列の場

もいい、行列が記憶領域に入りきれない場合などに合、その性質を存分に生かした専用の手法がありま使います。プログラム演習の題材にできそうな方法す。

ばかりです。詳細は [34] をお読み下さい45 。

3.7.1 Gauss の消去法, Cholesky 法

3.6.1 分割計算法

一般の Gauss の消去法, Cholesky 法を三重対角

正則な実行列の逆行列を、低い次数の逆行列を用

行列向きに書き直した方法です。

いて計算する方法です。計算時間の節約には全くな

通常の問題が O (n3 ) の計算時間を要するのに対りません。

し、三重対角の性質を使えば O(n) ですみます。

3.6.2 縁どり法

3.7.2 cyclic reduction 法

縁どり法 (bordering method)は分割計算法の変もとの方程式に適当な消去処理をすることによっ形で、 n 0 1 次の行列 An01 の逆行列がわかっていて、元数が半分の三重対角方程式に縮小 (reduction)

る場合、 An01 に 1 つの行と列をつけ加えた (縁どさせます。この方程式を更に半分に縮小することをりした) n 次の行列 An の逆行列を求める手法です。

繰り返し (cyclic)、最終的に次数を 1 次まで落して

3.6.3 エスカレータ法

しまって得られる解を逆にたどり、すべての解を求める方法です。
縮小処理はすべて並列に処理できるので、ベク

エスカレータ法 (escalator method) は、 J. Mor-

トル計算機、並列計算機にとても適した手法です

ris さんが 1946 年に発表した方法です。
k=B

で求める方法です。逆行列の値が必要な場合は使え

算や最小自乗法で大きい顔をする手法ですので、こ3.6

T と書けるとき、 A の逆行列を([70])。

を行列の次数とすると、低い次数から順番に

k k = bkA x

k = 1; 2; 1 1 1 ; n

3.8として方程式を解いていく方法です。

疎行列, 帯行列専用の方法

直接法の最大の弱点は、もとの行列 A を何らかの

42 J. P. Gram さん (1850{1916) と E. Schmidt さん (1876{

方法で変形するため、もともと 0 だったところがど

1959) の名前から来ています。

43 A. S. Householder さん (1904{?) が考案した方法です。

んどん汚されていく点です。

45 速さや精度はともかく、これらの連立 1 次方程式の数値解法

確保の問題が生じます。

そのため、大次元の行列の場合には、記憶領域の44 J. W. Givens Jr. さん (1910{?) が考案した方法です。

が提案された背景や数学的な証明を丁寧にフォローし、プログラ

LU 分解による行列の変化ムを組んで性能や問題点を調べるのは、結構面白いと思います。

次ページの図は、 2 次元正方領域で定常 NavierStokes 方程式の近似解を有限要素法で求めたときの306

行列の分布図です46 。□は成分が 0 でない場所を表します。■は成分が 0 のところです。有限要素法の分割を細かくすると、行列はどんどん \疎" の度合を増します。つまり■が多くなります。

LU

分解前の行列の形 (帯行列)

この帯行列に Gauss の消去法による LU 分解を
かけます。
LU

分解前の行列の形

この行列 A に Gauss の消去法による LU 分解をかけます。 L は A の下三角部分に、 U は上三角部
分にそれぞれ格納されます。

LU

分解後の行列の形 (帯行列)

帯行列の場合、 Gauss の消去法の過程で、対角線からある程度離れた成分 (帯の外側) は、ずっと 0 のまま変化しないことがわかります。
LU

LU

この性質を生かすと、計算時間と記憶場所を節約

分解後の行列の形

する帯行列専用のプログラムが書けます。アルゴリ
ズムは [60] に詳しくのっていますので参照下さい。

分解の結果、□の部分がほとんどを占めるよ

うになりました。このため、□の部分に対応した記Aが疎行列の場合も、 wavefront 法やスカイライ

憶場所を計算機内に確保しておく必要があります。

ン法、 DM 分解に基づく方法など、行列の特殊性

大規模な行列に何も考えずに直接法を適用するを生かした解法が開発されています。詳しくは [6],[51], [59], [60] などを参照下さい。

と、領域確保の問題で悩むことになります。
しかし、 A が帯行列ならば話が別です。次の図

3.9は、先ほどの LU 分解前の図から、対角付近の一部を抜き出した行列です47 。

現在人気の解法

46 同次境界条件で領域は四角形で等分割しています。近似は、LINPACK, LAPACK, IMSL などの著名な数値

流れ場が区分 2 次、圧力場が区分 1 次です。 Graphics は Math-

計算プログラム・パッケージの多くは、直接法のサ

ematicaの ListDensityPlot を使いました。

ブルーチンを主にサポートしています。

47 Laplace 作用素が 2 つ並んでいるみたいなものです。

307

SSL II

行列の形
一般の行列

NUMPAC

Gauss の消去法での LU 分解

Gauss の消去法での LU 分解

Crout 法での LU 分解

Doolittle 法での LU 分解

Householder 変換による QR 分解

Householder 変換による QR 分解

特異値分解法による最小二乗最小ノルム解
Bunch の方法での M DM T 分解

特異値分解法による最小二乗最小ノルム解

対称帯

Gauss の消去法での LU 分解

Gauss の消去法での LU 分解

三重対角

Gauss の消去法での LU 分解

Gauss の消去法での LU 分解

cyclic reduction 法

対称・帯

Bunch の方法での M DM T 分解

正定値・対称

修正 Cholesky 法での LDLT 分解

正定値・対称・帯

修正 Cholesky 法での LDLT 分解
修正 Cholesky 法での LDLT 分解

正定値・対称・三重対角

Bunch の方法での M DM T 分解

Bunch の方法での M DM T 分解
Cholesky 法での LLT 分解(修正)Cholesky 法での LLT (LDLT ) 分解
(修正)Cholesky 法での LLT (LDLT ) 分解
Cholesky 法での LLT 分解

SSL II, NUMPAC が採用している連立 1 次方程式の直接解法 (1995 年現在)

msec.
50000

その理由は、直接法が他の手法に比べて優れているからではなくて、ある最低の条件を満たせばサク
サクと答を出してくれる \汎用性" にあります。

DLSIX
DLAX
DVLAX

40000

反復法や共役勾配法は、解こうとする近似モデル

30000に依存して収束が速かったり遅かったりするので、個別撃破48 が要求される場合が多く、プログラム・

20000

パッケージとしてまとめるのは大変です49 。

10000

上の表は、九州大学大型計算機センターで現在公開している数値計算ライブラリの SSL II, NUMPAC00の中から、実係数の連立 1 次方程式に関する直接解法をまとめたものです。

1000

2000
3000
matrix dimension

4000

5000

なんと、記憶領域節約型の DLSIX が一番遅くなっ行列の形に対応していろいろな算法が提案されてています。次に対称正定値の場合を調べましょう。

いますが、やはり LU 分解などの \三角化に基礎を

DLSX, DVLSX は対称正定値専用のサブルーチンでおく方法" が圧倒的です。

す。 DVLSX はベクトル計算機向きにチューニング

この 中で どれ を 選ぶ かは (行 列の 形が 第一 で すされています。

が)、精度が極端に悪くなることを考えなければ、

msec.
80000

速い、または、記憶領域が節約できる

70000
DLAX
DVLAX
DVSLDL
DLSX

60000

ことが基準になります。

50000

例えば、対称行列専用のサブルーチンは、一般の40000

行列用のサブルーチンに比べ、半分の記憶領域 (行30000

列について) で実行できます。だからといって、実20000

行時間が半分になるとは限りません。

10000

ベクトル計算機 FUJITSU VP2600/10 で SSL II

0のサブルーチンを用いて実行速度を調べてみます。

0

1000

2000
3000
matrix dimension

4000

5000

DLSIX は対称行列専用、 DLAX, DVLAX は一般
DVLSX は期待通りの性能を見せましたが、 DLSX

行列用で、 DVLAX はベクトル計算機向きにチューはダントツの最下位におわりました。

ニングされています。

このことから、記憶領域の節約だけでなく、計算

48 要するに自分でプログラムを書け、ということです。

機の性能にあわせたチューニングも大切なことがわ

49 もっとも最近は、反復解法専用のプログラム・パッケージもかります。

出回るようになりました。

308

4

反復法の概要「反復法」と名がついているからには、何かの式

You will hardly ever again eliminate directly, atleast not when you have more than 2 unknowns.
The indirect procedure can be done while one ishalf asleep, or while thinking about other things.

を反復しながら収束する近似解を作るんだろうといGauss さんから Gerling さんへの手紙 (1823)

反復法は、直接法と比べむやみに種類が多く、とても分類する気になりません。ここでは、よく知られた反復法に限定して紹介します。

うことはすぐに察しがつきます。
Gauss さんの反復手法は、手紙をもらった Gerling

反復法でもっともよく使われるのは線形定常反復

さんからドイツの数学者のサークルに広まり、そ法 (linear stationary iterative method ) です。反復の成果に興味を持った Jacobi さんと、その弟子の列の作成は以下の通りです。

Seidel さんに引き継がれます51 。

線形定常反復法

なお、反復法についてのさらに詳しい歴史は [2]

適当な初期値 x0 から出発して,x

k = xk01 + R(b 0 Axk01 );を参照下さい。
k  1

(9)

4.2

を気がすむまで繰り返す. ただし R は A01 のJacobi 法は、ドイツの数学者 Carl Gustav Ja-近似であることが望ましい.

cob Jacobi さん52 (1804{1851) が 1845 年に発表したもので、 (9) で R = D 01 と選ぶ方法です。た(9) の \b0Axk 01 " は、 xk01 が解に近ければ小さくなることが期待される \残差" です。

だし、 D は A の対角成分を要素とする対角行列で実際、 R = A01 であれば、一発で反復は収束しす。
(9) を R = D 01 として変形すると、 Jacobi 法のます。もっとも、そんなことが最初からできるなら
苦労はしません。

反復列は線形定常反復法では、「どんな R を選べばいいの

k = (I 0 D01 A)xk01 + D 01 b

か」が最大のポイントになります。

4.1

Jacobi 法

xとなります。

反復法の創始者

Jacobi 法の収束については、以下のことが知られています ([60])。

反復法の歴史を遡ると、おなじみ Carl Friedrich

Jacobi 法の収束条件

Gauss さん (1777{1855) までたどりつきます。
Gauss さんは 1823 年に弟子の Christian Ludwig

Aが対角優位行列 ならば Jacobi 法は収束する.

Gerling さんにあてた手紙で、連立 1 次方程式を解
51 Gauss-Seidel 法の由来はここからきています。もっとも、くために反復法を使うことを推奨しています。

Forsythe さんによれば、現在 \Gauss-Seidel 法" とよばれる手

Gauss さんの考えた反復法は、次数が大きくなり

順 に つい て Gauss さ ん が 言及 し た 形 跡は ど こ に もな く ([8])、すぎて消去法での手計算が難しくなった場合、方程Sidel さんにいたってはこの方法についてふれたあと「この方法は式を 2 つ以上の系に分け、その系を順次解きながら近似解を生成し、望むまで繰り返すといった手法で使わない方がいい」と書いてあるそうです ([10])。
52 何といっても理科系の学生を悩ませる \Jacobian" で広く知した ([17])。 Gauss さんは手紙の中にこんなことをられています。強烈な性格だったらしく、人の反感を買うことが書いています50 :多かったそうです。晩年は財産を失い、健康を損ね、業績も振るわず、不遇のうちに亡くなりました ([73])。

50 もちろん、原文はドイツ語です。英語訳は [8] から抜粋しま

微積分学の勉強の際、 Jacobian でつまづく人が多いのはそのした。

ためかも知れません。

309

実際の計算において、 Jacobi 法は収束がかなり4.5

遅いため ([61])、収束速度を速める工夫が必要です。
なお、 Jacobi さんの連立 1 次方程式に対する取SOR 法

SOR 法 (successive over-relaxation method) は、り組みが [17] に詳細にまとめられていますので、興過剰緩和法、逐次過緩和法ともいいます。

味のある方はお読み下さい。

SOR 法は、 S. P. Frankel さんと David M. Young

さんがほぼ同時に発表した方法で、 Gauss-Seidel 法

4.3 Richardson 法の収束を加速させるため、加速係数 1 < ! < 2 を用いて R = (L + ! 01 D )01 とするものです。

Lewis Fry Richardson さん (1881{1953) は 1910

これを (9) に放り込むと、 n

年、 (9) の特別な場合として、 R を反復によって変化するスカラー k01 とおき、この手順を最適化す
る k01 の決定方法を提案しました ([32])。
また、適当な実数

でR =k = (D + !L)01 (!b +((10! )D 0!U )xk01 ) (10)x

T とおく手法もAとなります。 !

Richardson 法といいます。

4.4 Gauss-Seidel 法(10) を書き直すと、次の反復と同じです。
yvon Seidel さん (1821{96) によって発表され、 1874

SOR 法の収束については、 A. M. Ostrowski さんが 1954 年に示した有名な定理があります ([32])。

A = L +D +U

UOstrowski の定理
Aが正定値かつ対称で D + !L が正則ならば,0 < ! < 2 の範囲で SOR 法は収束する.

: A の上三角部分, ただし対角部分は 0

Gauss-Seidel 法は R = (L + D )01 と選びます。

最適な ! となる理論値はわかっていますが、実際これを (9) に放り込むと、にどのように選んだら収束が加速されるかは多分にケースバイケースで、これに関して多くの研究がな

k = D 01 (b 0 Lxk 0 U xk01 )x

されています。詳細は [32], [51], [61] を御覧下さい。

となります。 L が下三角行列なので、この式は、SOR 法の最大の利点は、プログラムが簡単に書k の i 番目を計算するとき、前回の xk01 を使う代
わりに、既に i 0 1 まで計算されている新しい xk のx

けることで、現在も流体解析の分野を中心に根強い人気を保っています ([55])。

要素を使用することを意味しています。

以下、 SOR 法の仲間を紹介します。

Gauss-Seidel 法の収束について、次のことが知ら

4.5.1 SLOR 法

れています ([54])。

Gauss-Seidel 法の収束条件
A

(11)

x

まず、行列 A を分離します。

の対角部分

k = D01 (b 0 Lxk 0 U xk01 )k = xk01 + ! (y k 0 xk01 )年にきちんとした形で著作となりました ([2])。

D : A

1 のとき Gauss-Seidel 法に比べ \過剰" に修正されるため、過剰緩和法の名前がつきました。

Gauss-Seidel 法は、 1862 年に Philipp Ludwigの下三角部分, ただし対角部分は 0

1 のときが Gauss-Seidel 法で

す。
! >

L : A

=

SOR の反復式 (11) の計算では、 D01 が必要でが対角優位行列、または対称かつ正定値行列す。しかし、 D は対角行列なので、逆行列は 1 をならば, Gauss-Seidel 法は収束する.

Dの各要素で割ったものとして簡単に求まります。

Gauss-Seidel 法は、 Jacobi 法に比べいくぶん収では、 SOR 法において、この D の部分が三重対束性がよくなっていますが、まだまだ改善の余地が角行列になるように A = L + D + U と分解したらあります ([32])53 。

どうなるでしょうか。ただし、分解された D は少な
くとも正則でないとだめです。
式は (11) と同じです。違うのは、三重対角行列に関する連立 1 次方程式

53 5.7 節に出てくる前処理技法は反復法にも適用でき、現在もD

盛んに研究されています。問題は並列化との相性です。

310

01

k 0 U x k 01 )(b 0 Lx

を直接解法で解く必要が生じるところです54 。
この 方法を SLOR 法(successive

line

もし、
A = A1 + A2 + Dover-55 ) または過剰線緩和法、逐次線
relaxation methodと分解できるとき、適当な行列 E1 , E2 をとれば、過緩和法といいます。

Ax = bSLOR 法は、直接法と反復法を組み合わせた方(A1 + D + E1 )x = b 0 (A2 0 E1 )x法で、収束性の向上、精度の安定などが報告され、(A2 + D + E2 )x = b 0 (A1 0 E2 )x

1950 年代終りから 1960 年代にかけてアメリカで大のどちらの式とも同値です。

いに使われました。

ADI 法は、 xk から一方の式を使い中間ベクトルkx^ を求め、次にもう一方の式を使い xk+1 を作る方

この成功に気を良くし、 SLOR 法を更に拡張した56 方法が S2SOR 法 (two-line successive overrelaxation method)

は

法です。

です。

ポイントは、うまい具合に E1 , E2 を選んで行列

なお、詳細な解説が [32] にありますので、興味のある方は御覧下さい。

A1 + D + E1 , A2 + D + E2

4.5.2 SSOR 法

い三重対角行列にするところです。
条件の良い三重対角行列を係数行列とする連立 1次方程式は、直接法でスラスラと解けます。このよ

SOR 法の反復行列は一般に非対称です。反復行

うに、 ADI 法58 は直接解法を一部採り入れた反復解列を対称に保つ方法が SSOR 法 (symmetric SOR
method)法です。

です。

なお、 [32] に 1 章まるまる割いた解説があります

後で出てくる Chebyshev 加速法と組み合わせてので、詳しくはこちらを参照下さい。

収束性を高めると、なかなかの性能を発揮します。
また、よく似た方法に USSOR 法 (un-symmetric
SOR method)

4.7があります ([51])。

4.5.3 SBOR 法
A

Aitken の  2 法

Alexander Craig Aitken さん (1895{1967) は、

数列 fxn g が極限に 1 次収束するとき、 fxn g よりを小行列のブロックに分割して、それぞれにSOR 法を適用する方法が SBOR 法 (successive block
over-relaxation method )もさらに速く収束する数列 fyn g を作り出す方法を提案しました。この変換方法を Aitken の  2 法と

いいます。

です。

n の作り方はいろいろありますが、例えば係数行列が対角成分に集まっている場合に威力をfy g

発揮します ([51])。
y

4.6 ADI 法

i = xi+2 0(xi+2 0 xi+1 )2

i

i

i

x +2 0 2x +1 + xなどで fyn g を求めます。これで反復法の収束を加速させることを狙います ([34])。

ADI 法 (Alternating Direction Implicit method) は、 D. W. Peaceman さんと H. H. Rachford Jr.

4.8

さんが 1955 年に発表した方法で、交互方向法ともいいます。

Chebyshev 加速法

Chebyshev 加 速法 も、反 復法の 収束 を速 める

ADI 法は、 2 次元の Laplace 方程式、 Poisson 方
程式57 、熱伝導方程式の陰解法などの差分法でよく有名な手法です。 Chebyshev 準反復法ともいいます。

Chebyshev 加速法の基本方針は、 \xk の計算に直前の xk01 しか使わないのはもったいないので、用いられます ([53])。
A1 , A2

を正定値で条件のいを対称かつ正定値かつ対角要素が正でそれそれより前の情報 xk02 ; 1 1 1 ; x0 も活用しよう" とい以外が非正の行列、 D を非負の対角行列とします。
54 厳密には直接法で解く必要はどこにもありません。が、直接うものです。

法がいちばん速いでしょう。
55 \線" の意味は、差分化した未知数の格子点の線を単位としてX

情報の使い方としては、反復列の一次結合

k
x
^ =

あつかうことに由来しています。
56 D のバンド幅を大きくした

57 Simeon Denis Poisson さん (1781{1840) に由来します。

k

j =0

jx

j

58 「交互方向法」という名称は、分割した格子を、まず水平方向に解き、次に垂直方向に解くところからきています。

311

に対し、 x^k の精度がもっともよくなるような j を決定することを考えます。

(

QR

この j は、 k 次の Chebyshev59 多項式T

k (t) =

分解などでも可能です。 A の条件がものすごく

悪くなければ、よりよい精度に改良されます ([60])。
また、 R の決め方は自由ですので、 A を LU 分
解に都合のいい行列と、それ以外 E に積極的に分け

cos(k arccos t);

jtj  1cosh(k arccos t);

jtj > 1て反復する方法もあります ([59])。
SSL II には、 LU 分解で得られた近似解の反復改良サブルーチン DLAXR がサポートされています。

の係数から導かれることが示せます ([61])。

Gauss の消去法での解法 DVLAX と組み合わせ、

この他、 j の決め方には CG 加速法などもあります。

精度がどれくらい改善されるかを調査しました。
計算機は VP2600/10、コンパイラは FORTRAN77
EX/VP V12L10 です。

4.9

8>
<
>:

モンテカルロ法

方程式はモンテカルロ法 (Monte Carlo method) は、逆行列の要素または連立 1 次方程式の解の成分を、期待値が望み通りの解になるような統計的変数のランダ
ム・サンプリングを行なうことにより、統計的に推X

ij = aji = n + 1 0 i
n0i+1
n0 k +1
bi =
n
k=1
a

ります。なお、誤差は相対誤差

反復法の出発値としてのごく粗い推定値を得たい場合などに用います ([34])。

max jnxi 0 1j
1inで評価しました。

反復改良法

relative error
10-7

さてここで、 (9) の行列 R は「A01 であることが望ましい」という、下手な仕様書みたいな条件を

DVLAX

10-8

思い出してみましょう。

10-9

浮動小数点演算の有界性より、直接法でもとめたAi  j)

で定義します。解は正確に x = (1=n; 1 1 1 ; 1=n) とな定する方法です。

4.10

(if

DVLAX + DLAXR

10-10の LU 分解は、ある \誤差" を含んでいます。こ10-11

れを行列の形で

10-12

A = LU + E

10-13
0

と書きます。
LU

分解によって求めた近似解 x^ は、反復法の絶

2000
3000
matrix dimension

4000

5000

反復改良により、約 1 桁精度が改良されていま

好の初 期値と いえま す。また、 LU の 逆行列 は、す。一方、実行時間は下の通りです。

「A01 であることが望ましい」という条件に適応しmsec.

た R の絶好の候補といえるでしょう。

50000

そこで、直接法と反復法を組み合わせた反復改良^、 R =
法 (iterative improvement) が、初期値を x

40000(LU )01 として次のように定義できます。
x

1000

DVLAX + DLAXR

30000

k = xk01 + (LU )01 (b 0 Axk01 )

20000
DVLAX

(LU )01 (b 0 Axk 01 ) の計算は、 c = (b 0 Axk01 )10000とすると、 LU x = c をみたす x を求めればいいの

00

で、前進代入・後退代入により求まります。

1000

2000
3000
matrix dimension

4000

5000

直接法によって求められた近似解の反復改良は
LU

分解に限らず、 LDU , LDLT , LLT , M DM T ,反復改良を追加すると、約 3 倍の実行時間がかかります。より精度のいい近似解が欲しいならば、反59 Pafnutii L'vovich Chebyshev さん (1821{1894) の名前か

復改良法を使うべきですが、その分コストがかかる

らきた有名な多項式です。

ことを覚悟すべきでしょう。
312

5共役勾配法の概要

共役勾配法の系統は、略字だらけで大いに混乱し脚光をあびるようになりました。
Hestenes さんと Stiefel さんも喜んでいることでますので、章の終りに一覧表をつけました。

しょう62 。

5.1

鮮烈なデビューとその後

共役勾配法60 は、 M. R. Hestenes さんと E. Stiefel
さんによって 1952 年に発表されました。

5.2

基本的な考え方共役勾配法の考え方は、次の性質が出発点です。
最小化問題

当時は電子計算機のための数値計算技術の確立が急務とされている時期で、共役勾配法は、反復法でAありながら有限回のステップで厳密解に到達するとが正定値かつ対称のとき, Ax = b を解くことは,いう性質と巧妙なアイディアが熱狂的に支持され、

f (x) =

「科学史上に残る大発明」として華やかなデビューを飾りました (当時の経緯は [52] に詳しく書いてあ12(x; Ax) 0 (x; b)を最小にする x を求めることと同値.

ります)。
( 1 ; 1 ) はベクトル内積

しかし、いざ共役勾配法を実際の問題に適用しようとすると、丸め誤差の影響を受けやすく、何回反(x; y ) = x

復してもいっこうに収束しない例も出てきました。

Ty =また、計算機の大容量化にあわせ、ライバルの直接法と SOR 法61 が理論と技法の両面から頑張ったとして定義します。

こともあり、華々しくデビューしたわりには、当時

0
BBB
B@

X
n

i=1

1
CCC
CA

i i

x y

同値性を確認しておきましょう。 f を微分すれば期待したほど普及しませんでした。
1960 年代の数値計算の本を見ると、連立 1 次方

程式の解法として Gauss の消去法、 Gauss-Jordan

rf (x) =法、 Jacobi 法、 Gauss-Seidel 法、 SOR 法は必ず紹介されていますが、共役勾配法をきちんと扱った本

@f
@ x1

.
.
.

@f
@x

はそれほど多くありません。しばらく不遇の時代がn= Ax 0 bとなります63 。

続きます。

f (x)

共役勾配法は、 1960 年代に非線形最適化問題のの最小点 x3 が見つかったとすると、 rf (x)解法として再登場し、 1970 年代には疎行列用の反での傾きはなくなるので、 x3 での f の微分は 0 と復解法としての研究が行われました。そして 1980

なり、つまり Ax3 0 b = 0 となります。
また、任意のベクトル h 6= 0 に対し、 A の対称性

年代に入ると、 \前処理"という新しい手法を付け加を使ってえることで、大規模行列に対する反復法として再び
60 「共役」は「きょうやく」と読みます。もともとは軛 (くびf (x + h) = f (x) + (h; Ax 0 b) +き) をともにして車を引く意味の「共軛」からきています。昔、数学の講義で、とある教授が「これを《きょうえき》と読むのは数

12(h; Ah)62 よく歌や映画にある筋書きで、売れない歌手の熱狂的なファ学が全くわかっとらん証拠じゃ〜」とおっしゃっていたのを記憶ンだった女の子が、華やかなスポットライトを浴びてどんどんメ

しています。数学者の前で発表するときは気をつけましょう。
61 一般論ですが、「理論」や「手法」などは、科学とか論理のジャーになっていく彼を悲しそうに見つめながら離れていくというのがあります。 1960 年代に共役勾配法を一所懸命に研究した

土俵をはみ出して、政治的学閥的立場、個人的好き嫌い、宗教観

人も、今はそんな気持ちなのでしょうか。
63 ピンとこない方は、 f (x) を全て成分の形に書き下して、一などで争われることがあります。しかし、動機はともあれ、結果
がみんなのためになればそれでいいのです。

つ一つ偏微分してみて下さい。

313

となることが簡単にわかります。
は正定値でしたので (h; Ah) > 0 となって64 、
3f (x f (x+ h) > f (x )3f (x) の値は f (x ) より大きくなり、このことから
3 が f の最小値を与えることがわかります。

!

具体的な例をあげましょう。
A =

2

01

01

2

b =

;

1

!

1

0+

p

0 ) が最小になるようなスカラーを計算し、それを 0 とする。

3となります。つまり、 x3 から少しでもはずれると、x

初期値 x0 と探索方向 p0 を決める。



ここで、 x3 が Ax3 = b をみたすと考えれば、 A



近似解を x1 = x0 + 0 p0 とする。



次の探索方向 p1 を決める。



以下、気が済むまで繰り返す。

最小値が欲しいは、 2 次式の最小化と同じになるので、簡単に求めることができます ([61])。

;

重要なのは、探索方向となるベクトル pk を如何に上手に決定するかです。この選び方によって、収とすると、 A は対称かつ正定値で、 Ax = b の解はT です。

束率が決定的に変わってきます。

x = (1; 1)

また、 f (x) は
2

5.4

2

f (x) = x1 + x2 0 x1 x2 0 x1 0 x2

最急降下法 (steepest descent method) は、探索方となります。 f をごちゃごちゃ変形すると
f (x) =

(x1 0 1)2 + (x2 0 1)2 + (x1 0 x2 )22

向 pk を xk における勾配方向kkk = 0rf (xk ) = ( @ f (x ) ; 1 1 1 @f (x ) )T
1
n
@x
@x

0 1

p

ですので、最小値は x1 = x2 = 1 のとき 01 となるにとる手法です。 SD 法ともいいます。

ことがすぐにわかります。
f (x1 ; x2 )

最急降下法の (1; 1)T のまわりの形を書くと以下の最急降下法は、 A の性質が悪い場合、収束がかな
り遅くなることから、あまり使われません。

ようになります。

5.5

共役方向法

共役方向法 (conjugate direction method) は、さらに上手な pk の設定方法です。 CD 法ともいいます。
探索方向の選び方は、これまでの探索方法と \共役" となる方向にとっていきます。 \共役" とはこのf (x)場合、 A に関して直交する方向を指します。ちゃんの形状と等高線図と書くと、この図で最小値となる x = (1; 1)T が Ax = b の共役の定義解です65 。

5.3

ベクトル x と y に対し(x; Ay ) = (y; Ax) = 0

計算の手順が成り立つとき, x と y は A に関して共役であ共役勾配法とは、適当な初期値を出発し、順番にる.

探索方向を定め、その方向に沿って最小値を探しに行く手法です。

となります。

具体的な手順は次の通りです。

つまり k 番目までの方向 pk が決まっているとき、

664 正定値の定義より h = 0 に対し hT Ah > 0 です。

k + 1

65 連立 1 次方程式を同値な最小問題に持ち込むこの方法は、共番目の方向 pk+1 は、その全てと共役になるよ

うに、すなわち役勾配法に代表されるような巧妙な理論を生みました。しかし、

k+1 ; Api ) = 0;

直接法や反復法に比べて最小化する近似解に「副産物」的な印象(pを受けるためか、敬遠する人もいます。

314

i = 0; 1 1 1 ; k

で決めます。

からです。

具体的な方法は [52], [61] を参照下さい。

従って、共役勾配法は高々 n 回の反復で厳密解に探索方向 pk を共役方向に定めるということは、収束します。

理論的に大変重要です。それは、次の定理が成り立といっても、これは理論上の話です。実際の計算では丸め誤差が混入するため、必ずしも n 回の反復

つからです。
共役方向法の収束定理
Aが対称かつ正定値ならば, 共役方向法は高々n回の反復で厳密解に到達する.

で解が求まるとは限りません。
共役勾配法は共役傾斜法、 CG 法ともいいます。

5.7

この、有限回の反復で厳密な解が得られるという前処理付き共役勾配法

前処理付き共役勾配法 (preconditioned conjugate

理論的保証が、共役方向法の最大の魅力です。

gradient method )5.6は、 PCG 法ともいいます。

手法は、適当な正則行列 C によって Ax = b を

共役勾配法

(C

真打ち共役勾配法 (conjugate gradient method)

01

AC

0T

)(C

T x) = C 01 b(12)は、共役方向法の一種で、次のように探索方向 pkと同値変形し67 、 (12) に対し共役勾配法を適用するを決定する手順を指します :ものです。

最急降下法で選んだ探索ベクトル pk は
0rf (x(12) の処理を前処理(preconditioning)、行列 C を
前処理行列(preconditioner)68 といいます。

k ) = b 0 Axk

でした。 b 0 Axk は、方程式の k 回目の残差と考え前処理をする理由られます。この残差ベクトルをr

共役勾配法の収束性は、解を構成する固有ベクトルに対応した固有値の分布によって大きく変わりまk = b 0 Axk

す。
係数行列の固有値が重複したり、狭い範囲に密集とおきます。

共役勾配法の探索方向 pk は、この勾配方向 r k か

していれば、収束が速くなることが知られていますら Apk01 の成分を除去し (pk ; Apk01 ) = 0 となるよ([56])。つまり、

うに決めます。つまり、pk = rk + k01 pk01

前処理は、行列の固有値分布を改善し、収束性

kk 01 )
k01 = 0 (r ; Apk01
k 01 )(p; Apを向上することを目的となるように決めます66 。

とした手法です。

このようにして選んだ pk には共役直交性

ij(p ; Ap ) = 0;i 6= jが、また rk には直交性

i j(r ; r ) = 0;

扱いにくそうな行列i 6= j+が帰納法を使って示せます ([54])。

これにより n 回反復すれば r n = 0 となる、すな

わち xn が Ax = b の解になることがわかります。
その理由は、n + 1 個の 0 でない n 次元ベクトルが互いに直交することはありえない。

前処理によって改善された行列
67 C 0T は (C T )01 = (C 01 )T の意味です。

66 本当に直交するかどうか計算してみればわかります。

68 お肌や髪の状態を整えるためにも使うそうな...

315

前処理付き共役勾配法は、 1980 年以降、偏微分

多くの問題で ICCG 法よりも収束が速くなること

方程式の離散近似から導かれる連立 1 次方程式の解が報告されています ([59])。

法に用いると効率的な近似解が得られることから、盛んに用いられるようになりました。

5.8

が、前処理技法の将来には並列計算機が立ち塞がります。どうやら、前処理技法と並列化技法はあま

共役残差法

共役勾配法は A が対称かつ正定値行列であること

り相性がよくないようです。

を前提とした方法でした。

5.7.1 ICCG 法が数多く提案され、各分野で大きな成果をあげてい

近年、共役勾配法を非対称行列へと拡張する手法
ます。以下、名前だけ簡単に紹介します。

前処理行列の選び方は、個々の問題によっていろ

共役残差法 (conjugate residual method) は、最小いろ考えることができます。

化する関数を残差ベクトルの平方和

その中でも、とくに一般的な効果が期待できる手法に不完全 Cholesky 分解 (incomplete Cholesky
decomposition)f (x) =があります。

これは、疎で正値対称な行列 A を、 \疎" の性質

1
2

として、今度はをあまり失わないように、適当な下三角行列 L と残差行列 E で

(Ax 0 b; Ax 0 b)

(Ap

T +E
A = LLk+1 ; Apk ) = 0となるように探索方向 pk を決定する方法です。

または共役残差法は CR 法ともいいます。詳しくは [59],
^D
^L
^
A = L

T + E;^
L = LD

1=2[61] を御覧下さい。

5.8.1 ILUCR 法と分解し、 L を前処理行列 C とする手法です。
まともな Cholesky 分解、修正 Cholesky 分解と違共役残差法にも当然前処理を施すことを考えまい、あらかじめ指定した成分以外は 0 とし、それ以す。前処理付き共役残差法は、頭に \preconditioned"

外は誤差行列 E に押しつけるという身勝手さのたがついて PCR 法ともいいます。

め69 、 \不完全 (incomplete)" と呼ばれます70 。

今度は A が対称でないので、 Cholesky さんの出

特に LLT , LDLT 分解した行列から (12) を構成

番はなく、 LU 分解の登場です。

する方法を ICCG 法 (incomplete Cholesky conjugate gradient method)

非対称行列 A を、 \ぎっちり"LU 分解しないで、、または不完全 Cholesky

都合のいい下三角行列 L と上三角行列 U を使って

分解付き共役勾配法 といいます。
ICCG 法は 1977 年に J. A. Meijerink さんと H.

A = LU + E

A. van der Vorst さんによって発表され ([22])、物と分解します。 E は誤差行列です。

理・工学などの分野で広く用いられるようになりま

この分解を不完全 LU 分解 (incomplete LU de-

した。

composition)

5.7.2 MICCG 法(LU )このあたりから記号が錯綜してきます。

MICCG 法 (modi ed ICCG method) は、行列
Aといいます。

ここで、
01

01

Ax = (LU )bに対し共役残差法を適用します71 。

の対角要素に重みづけをし、修正項を付加してこの手法を ILUCR 法 (incomplete LU

ICCG 法の加速を狙ったものです。

この行列の分解方法を修正不完全 Cholesky 分

gate residual method )conju-、または不完全 LU 分解付き共役残差法と呼びます。

解 (modi ed incomplete Cholesky decomposition)

71 A が対称の場合は対称性を保つために C 01 A(C 01 )T とすといいます。

る必要がありましたが、この場合 A はもともと非対称なので、左から (LU )01 をかけるだけです。

69 仕事でこんなことをやったら、間違いなく嫌われます。
70 Cholesky さんの名前は、数値解析の世界で知らない人はいないくらい有名になりましたが、「修正」「変形」「不完全」な
ど、ちょっと気の毒な言葉が頭にくっつきます。

316

5.8.2 MILUCR 法

5.10.1 ILUBCG 法

MILUCR 法 (modi ed ILUCR method) は、

ILUBCG 法 (incomplete

LU

BCG

method)共役勾配法のときの MICCG 法と同じく、行列 A

は、不完全 LU 分解で前処理を施した双対共役勾配の対角要素を重みづけることで修正項を付加し、法です。
同様に自乗共役勾配法に対する ILUCGS 法 (in-ILUCR 法の加速を狙ったものです。

この分解を修正不完全 LU 分解 (modi ed incomplete LU decomposition )complete LU CGS method)といいます。

5.10.2 MILUBCG 法

5.9 GMRES 法

修正不完全 LU 分解によって ILUBCG 法の加速を狙ったのが MILUBCG 法 (modi ed ILUBCG

行列 A が対称であっても必ずしも正定値でない場

method)合、 Lanczos 法と共役勾配法との関係をうまく利用した方法に MINRES 法 (minimal residual method)という方法があります72 。

この方法 を非対称行 列に拡張し た方法が GMRES 法 (generalized minimal residual method) でです。

ILCGS 法の改良版の MILUCGS 法 (modiILUCGS method)

5.11

す。具体的な手法は、最初に提案された論文 [26] をedもあります。

QMR 法

双対共役勾配法 (BCG 法, BiCG 法) は演算の途御覧下さい73 。

5.10もあります。

中でゼロ割りが発生して計算が破綻 (breakdown) したり、数値的に不安定となる可能性があります。
QMR 法 (quasi-minimal residual method)は、残差双対共役勾配法

ノルムの最小化により双対共役勾配法の不安定性を双対共役勾配法 (biconjugate gradient method )克服し、高い収束性を狙った手法です。詳細は [15]は、対称でない行列 A に対し、対称行列用の共役を御覧下さい。

勾配法を修正して適用するもので、 BCG 法またはBiCG 法ともいいます。

5.12

原理は、

i j

(Ap ; p^ )=0i i

(Ap ; p^ ) 6= 0(r ; r^ ) = 0;(r ; r^ ) 6= 0;

i

j

i

i

自乗共役勾配法 (CGS 法) の収束の不規則性を改i 6= j

善する目的で提案された手法に BiCGStab 法があ

i = j

ります。

となるベクトル列 fpk g, fr k g と fp^k g, fr^k g を反復\BiCGStab 法" は無理に展開すれば \biconjugate
gradient stabilized method" でしょうが、最初に提過程で作ることにあります。

案した H. A. van der Vorst さんによれば \CGS 法

なお、詳細は [59] を御覧下さい。

によく似ていて、安定性 (stability) が得られ、しか双対共役勾配法の変種に自乗共役勾配法 (conju-

も BiCG 法とも関係するので、 BiCGStab と名付

gate gradient squared method) があります。自乗

けた" そうです。

共役勾配法は CGS 法ともいいます。

BiCGStab 法にも前処理が適用できて、 van der

BCG 法と CGS 法は本質的には同じものですが、

Vorst さんは BiCGStab-P と命名しています。詳BCG 法よりも CGS 法の方が収束が速いことが報告

細は [31] をお読み下さい。

されています ([55])。

以上、共役勾配法に属する手法を用語の意味だけ

この二つに前処理を施すと、先頭に \precondi-

概観しました。ここで紹介したものは、雑誌や論文

tioned" がついて、 PBCG 法、 PCGS 法になりまでよく見かけるもののほんの一部です74 。

す。

特に非対称行列に対する手法は、現在盛んに研究

72 もう少し具体的に言うと、 A の固有値問題を解く方法でよく知られている Lanczos 法 ([72]) は、 v1
とき Krylov 空間 Kk=fが行なわれており、新しい手法が次々に提案されて

r0 =jjr0 jj とした
1 1 1 ; Ak01 v1 g 上の=span v1 ; Av1 ;BiCGStab 法います。
74 残念ながら、共役勾配法系統で現在 (1995 年 11 月) センター

Galerkin 法であり、この方法は Ax = b を解く共役勾配法に対応している ([26]) という関係です。
73 その他 FGMRES 法, GMRESR 法などがあるそうですが、で公開しているサブルーチンライブラリはありません。

ちゃんと調べてません。

317

\M" は \modi ed" の M， \IC" は \imcomplete Cholesky" の \IC"， \ILU" は \imcomplete LU" の \ILU"，\P" は \preconditioned" の \P" だと覚えておけば、共役勾配法関連のややこしい略語が整理できます。

略語

英語名

日本名

SD 法
CD 法
CG 法
PCG 法
ICCG 法MICCG 法steepest descent method
conjugate direction method
conjugate gradient method
preconditioned conjugate gradient method
imcomplete Cholesky conjugate gradient methodmodi ed imcomplete Cholesky conjugate gradient method

最急降下法

CR 法
PCR 法
ILUCR 法MILUCR 法

conjugate residual method
preconditioned conjugate residual method
incomplete LU conjugate residual methodmodi ed incomplete LU conjugate residual method

共役残差法

BCG 法
PBCG 法ILUBCG 法
MILUBCG 法

biconjugate gradient method
preconditioned biconjugate gradient method
incomplete LU biconjugate gradient methodmodi ed incomplete LU biconjugate gradient method

双対共役勾配法

CGS 法
PCGS 法ILUCGS 法
MILUCGS 法

conjugate gradient squared method
preconditioned conjugate gradient squared method
incomplete LU conjugate gradient squared methodmodi ed incomplete LU conjugate gradient squared method

自乗共役勾配法

QMR 法
MINRES 法
GMRES 法BiCGStab 法

quasi-minimal residual methodminimal residual method
generalized minimal residual methodbiconjugate gradient stabilized method3

準最小残差法 3

共役方向法
共役勾配法
前処理付き共役勾配法
不完全 Cholesky 分解付き共役勾配法
修正不完全 Cholesky 分解付き共役勾配法
前処理付き共役残差法

不完全 LU 分解付き共役残差法

修正不完全 LU 分解付き共役残差法
前処理付き双対共役勾配法

不完全 LU 分解付き双対共役勾配法

修正不完全 LU 分解付き双対共役勾配法
前処理付き自乗共役勾配法

不完全 LU 分解付き自乗共役勾配法

修正不完全 LU 分解付き自乗共役勾配法
最小残差法 3
一般化最小残差法 3
双対共役勾配安定化法 3【注意】 * は一般的な呼び方がわからなかったので適当に名付けたものです。

響を受けやすい性格が災いして、何回反復しても少上に、やたらと出てきた略語の一覧をのせます。

しも先に進まないという事態に陥ることでも有名で

5.13

した75 。

共役勾配法の復活

収束が遅い場合の対策として、反復が停滞してし1952 年に登場した共役勾配法は、当初の熱狂的な

まったら一旦反復を打ち切って再スタートする方法

ブームが去ると、一時は全然顧みられない時期があや、初期値をうまく改良する方法などが工夫されまりました。

したが、なかなか決め手を欠きました。原因は、丸め誤差に対する繊細さにありました。

たとえば、偏微分方程式についての古典的名著といわれる [11] では、次のように扱われています。

共役勾配法は丸め誤差に弱い。それは「n 回反復したのに残差が 0 にならない」という形で検出される。また、共役勾配法の計算を、単精度と倍精度で行なって比較してみると、単精度の場合には丸め誤差のために、収束 (事実上の) がかなり遅くなっていることがわかる。少ない桁数で計算すれば丸め誤差が入るのは当然であるが、共役勾配法は他の解法に比較して、丸め誤差の影響が大きいようである。

Hestenes と Stiefel (1952) の共役勾配法は有限回反復の族の 1 つで、構造は反復的でそれにもかかわらず有限回のステップで終了する。比較的大容量の記憶装置を必要とするのと、各反復段階が比較的複雑な構造のためにいままでのところ偏差分方程式には広く採用されなかった。
G. E. Forsythe and W. R. Wasow (1960)

戸川 隼人 (1975)『偏微分方程式の差分法による近似解法』『共役勾配法』

ただし、かなりの数の数値実験より、共役勾配法は\うまくいけば"SOR 法や Chebyshev 加速法と比較

75 倍精度演算が桁違いに遅かった時代は、やはり単精度演算がして格段に速く収束することが知られていました。

多く使われました。それが共役勾配法の評判を落した原因でもあります。理論的には高々 n 回の反復で収束する \はず" なのに全そのかわり、うまくいかない場合は、丸め誤差の影

然収束しなければ、当然イライラするでしょう。

318

うまく働くとされている反復法が成功する理由は、ほとんどの場合、それが解こうとしている実際の問題と密接に結び付いた方法になっているからである。したがって、反復解法のサブルーチンをライブラリの中で探しても無駄である。反復法自体は数学的に単純であるが、行列 A の構造は問題によって複雑になり、また、問題に特有の形をしていることが多い。

結局のところ、共役勾配法の持つ収束性のバラ
ツキを克服するためには、もとの方程式を適当に変数変換する、つまり何らかの行列を引っかけること
で、性根を叩き直すしかないという結論に達し、いかに叩き直すかの研究が行なわれました。
Meijerink さんと van der Vorst さんが 1977 年

G. E. Forsythe, M. A. Malcolmに発表した前処理の手法 ([22]) は、共役勾配法の収and C. B. Moler (1977)

束性を劇的に改善することが明らかとなり、その後『計算機のための数値計算法』の理論的な研究の発展と前処理技法との組合せによって、共役勾配法は大型疎行列を係数とする連立
1 次方程式の近似解法として、鮮やかな復活を遂げ

5.15

ます。
そして現在、前処理付き共役勾配法は次のような

INSPEC を用いた実験

最後に一つだけ、共役勾配法の復活を示すデータ

評価を得るまでになりました ([57])。

をあげます。
九州大学大型計算機センターでは、 1969 年以降

前処理の導入により、 CG 法はガウス消去法などの直接解法をメモリ使用量、計算速度の面で凌いだ。以来、 CG 法系の行列解法は各方面で研究され、正定値が保証されない、弱い非対称行列を解くことが可能な解法まで開発されている。現在ガウス消去法の優れている点は、よほどの悪条件行列でなければ解ける点である。一方、 CG 系の計算方法はパイプライン型スーパーコンピュータにも適しており、今後ますます使用されることであろう。

の物理、電気・電子工学、計算機科学、制御工学、情報工学関係の文献のデータベース INSPEC をサービスしています ([66])。
そのデータベースに、以下のキーワードで検索をかけました (1995 年 10 月現在)。

羽根 邦夫 (1987)『共役勾配法の計算方法を利用した\conjugate gradient" \linear equation"

\SOR" \linear equation"

\Gauss" \linear equation"

半導体デバイスシミュレータ』

実際には \equation" は \equations" を含みます。
また \Gauss" は \Gauss" と名前がつくすべてのも5.14のが対象です。従って、 \Gaussian" や \Gauss-Seidel"

反復解法の汎用性

や \Gauss-Jordan" も含んでいます。
このキーワードでヒットした件数を年代別にまただし、上の評価は、 \特定の偏微分方程式を特とめました。なお、年代別の総データ数が異なるた定の手法で離散化して得られた大規模疎行列に対しめ、 100 万件に対して何件の割合かという値に変換て" という付帯条件をつけるべきでしょう。現在のしています。

ところ、共役勾配法がすべての直接法と反復法に対しあらゆる面で優れているとはいえません。
結局のところ、設定した数学モデルと近似モデルから導かれる連立 1 次方程式の性格によって、適切``conjugate gradient''``SOR''``Gauss''400

な手法を選択するしかありません。

328

316

その場合、反復法と共役勾配法の難点76 は、特に

394

300

スーパーコンピュータの性能を極限まで引き出すには、専門家が作ったライブラリが充実している直

193

200

接法と違い、何らかの形でチューニングしなければならないという点です。理由は以下の通りです

106

100

([13])。

25.5

52.1

52.5

1980s

1990s

25.5

0

76 ということは論文のネタになりやすいところ

1970s

INSPEC のヒット件数の推移 (100 万件あたり)319

1970 年代は \SOR" と \conjugate

gradient" は同じ数でした。 1980 年代には共役勾配法に言及した論文の数が急速に増加しています。 1990 年代になると、 SOR 関係は 80 年代と大差ないですが、共役勾配法関係はさらに数が増えています。

1979 年当時のセンターオープン室

320

6 Gauss の消去法のアルゴリズム
ここまで、連立 1 次方程式の解法によく使われる

行列 A は一般の非対称な密行列で、もちろん正則

方法を直接法、反復法、共役勾配法の順にざっとみだとします。
60 と仮定します。第 1 式に
まず、 (13) で a11 =てきました。

(1)i = ai1 =a11 をかけ、 x1 の係数である a21 から
n1 までの成分を i  2, j  2 について

当然、これらの中でどれが最も優れた解法なの

m

か知りたくなります。しかし、いざ比較するときaは、行列 A の条件 (次数や形状など) や、計算機の

(1)

(1)
= aij 0 mi a1j ;

(1)

(1)
= bi 0 mi b1

a

性能、欲しい精度など、いろいろな条件を加味しなくてはならないので、誰もが納得するような比較を

ij

i

b

するのはそんなに簡単ではありません。

の過程で消去します。

逆に、これだけたくさんの手法が淘汰されずに

0
BBB
B@

残っている77 という事実は、連立 1 次方程式の決定

10
CCC BBB
CA B@

1 0
CCC BBBCA B@

すると方程式は、 (13) と同値な的な解法が未だ提案されていないことの証明でもあります。
今後は、超並列計算機・ベクトル計算機・ベクトル並列計算機などの、次世代スーパーコンピュータと親和性のある解法が生き残ることでしょう ([43],[45], [56])。

n

a11

a12

111

a1

0
.
.
.

a22

111

a2

0

(1)

.
.
.
(1)
an2

(1)

n

.
.
.
(1)

:::
111

nn

a

x1
x2

=

.
.
.

n

x

b1

(1)

b2

.
.
.
(1)

n

b

1
CCC
CA

の形になります。

さてこれから先は、最も有名、かつ、行列データ

一般の形で書くと、第 k 段目の消去は演算は、を全て主記憶に格納できた場合に最も汎用性があるaといわれている Gauss の消去法について考えたい

(k01)

kk

6= 0

と仮定するとき、 i  k + 1, j  k + 1

に対しと思います。

(k )(k 01) (k01)i = aik =akk(k )(k 01)(k) (k 01)
aij = aij
0 mi akj(k)(k01)(k) (k 01)
bi
= bi
0 mi
b
k

Crout 法も Cholesky 法も、結果的には LU 分解に帰着されるので、 Gauss の消去法78 の仲間としてm

考えることにします。
なお、この章は、以降で使用する Gauss の消去法のアルゴリズムと用語のための解説です。本題はで計算します。

最後の章なので、ご存知の方は読み飛ばして結構で(k01)

消去の段階における akk

す。

をピボット(pivot)または枢軸といいます。

6.1 Gauss の消去法

X

0
BBB
B@

n

a

ij xj = bi ;i = 1; 1 1 1 ; n

10
CCC BBB
CA B@

最終的に次の式になります。

連立 1 次方程式を

j =1

この消去を第 1 式から第 n 0 1 式まで行なうと、

(13)

a11

n

a12

111

a1

a22

111

a2

(1)

(1)
.
.
.

:::

0

n

(n01)
nn

a

x1
x2

.
.
.

x

n

1 0
CCC BBB
CA B@

b1

(1)

b2

=

.
.
.

b

(n01)
n

とします。
この方程式は (13) と同値です。

77 「速い」ことが正義だとすれば、既に淘汰された手法もかな(n01)

従って、 ann

6=0 ならば、ただちに n 番目のりあります。
78 Gauss-Jordan の掃き出し法は、 LU 分解法でないので、仲未知数 xn が求まり、これを直前の式に放り込むこ

間に入れてあげません。

とで xn01 がわかります。
321

1
CCC
CA

0
@

1
A

以下、順番に k = n 0 1; 1 1 1 ; 1 に対し

k=

x

b

(k01)

k

X
n

0

j =k+1a(k 01)kjjx

これが最後の消去過程です。
(k 01)kk=aで解が求まります。
以上が、基本的な Gauss の消去法の手順です。

6.2

消去過程

Gauss の消去法による行列 A の変化を図で見ま

しょう。

消去完了
(n01)

消去過程が終了すると右下に ann

6.3が出現します。

LU 分解との関係

先ほどの消去過程では、右辺 b の値も A と一緒に変形しました。しかし数値計算では、同じ A に対して何度も違う b を使って連立 1 次方程式を解くことがよくあります。その場合、先のアルゴリズムをそのままプログラムとして書くと、 b にあわせて何度
2 段目の消去過程

も同じ A を変形することになって、効率が悪くなります。

白い部分は既に消去が完了した部分です。残りのグそのため、まず行列 A を LU 分解し、三角行列

レーの部分は参照・代入があります。
Lと U を用いて b を変形するプログラムに書き直します。このようにすると、異なる b が与えられた
場合、以前の LU 分解で蓄積した情報を駆使してスムーズに解を計算することができます79 。
実は、 LU 分解は Gauss の消去法の計算で既に完了しています80 。

(k )

消去過程で出てきた miを、消去における乗数(multipliers) といいます。 LU 分解の下三角行列 L

0
BBB
BBB
B@

は、消去における乗数を順番に並べた行列
3 段目の消去過程

段が進むにつれて、参照される小行列が縮小してい
きます。

1
(1)
m2
(1)
m3
.
.
.
(1)
m

n

1
(2)

m3

.
.
.
(2)

m

n

0
1
..

111

.
(n01)

n

m

1

1
CCC
CCC
CA

となることが証明できます ([50], [60])。
一方、上三角行列 U は、消去法によって最終的に変形される行列です。従って、 Gauss の消去法によって L と U の値がすべてわかります。ということは、 LU 分解ができるわけです。
79 世の中に出回っている著名な Gauss の消去法関係のライブラ
リはほとんどこのようになっています。
80 そのため、 Gauss の消去法と LU 分解法は最近同じものとし
n0 1

段目の消去過程て扱われます。

322

0
BBB
B@

プログラムを組む場合は、記憶領域の節約のため、最初に A を記憶した場所に L と U を上書きします ([69])。
LU
z

分解が完了すれば、あとは作業用のベクトルを用意し、Lz = bを z について解いたあと、

a12

111

a1

0
.
.
.

a22

111

a2

0

(k01)
max jajki;j n ijの順に、次を計算する。

Xがピボットとなるように選ぶ方法です。

i01

この方法を完全ピボット選択法 (complete pivot-

ik zk

l

k=1

ing strategy)!

の順に、次を計算する。

X
n

i0

k=i+1

ik xk

u

、または、完全軸選択法 といいます。

0
BBB
B@

後退代入計算znna(前節のグレーの部分) すべての成分を調べて、その前進代入計算i=

111

中で絶対値最大の具体的な計算方法は、次の通りです。

x:::

2 つ目 の方 法は、 残って いる 方程式 の係 数行 列

後退代入 (backward substitution) と呼びます。

i = n; 1 1 1 ; 1

n

.
.
.
(1)

置換行列 P を引っかけるということです。

前者を前進代入 (forward substitution)、後者を

i = bi 0

.
.
.
(1)
an2

(1)

行を入れ換えるということは、 A に左から適当なを x について解けば、解が得られます。

z

(1)

1
CCC
CA

部分ピボット選択の対象となる成分 (2 段目)

Ux = zi = 1; 1 1 1 ; n

n

a11

ii

=u

n

a11

a12

111

a1

0
.
.
.

a22

111

a2

0

(1)

.
.
.
(1)
an2

(1)

n

:::
111

.
.
.
(1)

nn

a

1
CCC
CA

完全ピボット選択の対象となる成分 (2 段目)

6.4

完全ピボット選択法のためには、行の入れ換えにピボット選択プラスして列の入れ換えも必要です。これは、 A に(k01)

前節のアルゴリズムでは、ピボット akk

適当な置換行列 P1 を左から、また、置換行列 P2 をが0で右から引っかけることを意味します。

ないことを仮定していました。しかし、ピボットが

この二つのピボット選択法の比較は、最終章の中0 となる場合も十分ありえます。また、ピボットの心の話題です81 。

絶対値が小さい場合も、浮動小数点演算でオーバー

2 章で紹介した LU 定理は、 A が正則なとき、こフローが起き、次の段階に進めないことも考えられのような分解が可能なことを保証しています。

ます。
その場合、何らかの方法で 0 でない、ついでに、

6.5

数値的に最良のピボットを探す工夫が必要です。
数値解析の分野で多くの創造的な仕事をされた

ピボット選択の利点(k01)

James Hardy Wilkinson さん (1919{1986) は、 1961

ピボット選択を行う最大の理由は、ピボット akk

年に発表した論文 [35] で 2 つのピボット捜索法を提が 0 になり、 Gauss の消去法の処理が破綻するためのを防ぐためです82 。

案しています。

しかし、 Wilkinson さんの提案したピボット選択

1 つ目は、 k 段目の消去を行なう前に、 k 列の第
k

(k01)

法を採用した場合、たとえ akk

行から n 行までを調べ、その中で絶対値最大となるような(k 01)

max jakin ik

6= 0であっても、

行や列の入れ換えが生じる可能性があります。
その場合、この操作は無駄なのでしょうか？

j

81 ただし、 Wilkinson さんは、最大の要素をピボットに選ぶ方がピボットとなるように行を入れ換えてから、次の法が最良の方法であるとは言っていません。しかし、現在のとこ

段階の消去を行なう方法です。

ろ、これに代わりうる実用的なピボット選択法は提案されていまこの方法を部分ピボット選択法 (partial pivoting

せん ([34])。
82 実数を 0 でわる行為は、実数体に対する挑戦です。

strategy)、または、部分軸選択法 といいます。
323

また、対称かつ対角優位で対角要素が正の場合、

そうではありません。ピボット選択法は、次節で正定値であることが知られています83 。

述べるような、大変都合のいい行列であることがあ
らかじめわかっている場合を除けば、是非組み込むべき手法です。

6.7

その理由は次の通りです ([72])。



ピボット選択によって、 Gauss の消去法が破綻す(k )ることなく、しかも丸め誤差の増大が抑えられるこ

ピボットの絶対値が小さいと、乗数 miとがわかりました。

の計算で桁あふれを起こす危険がある。
 ja

(k01)

kk

j

しかし、これには「各行列成分の絶対値の大きさの値が小さいのは、前の消去の段がだいたい同じ程度であれば」という前提条件がつ階で桁落ちをした可能性が大きい。
 jm

(k)

i

j

(k 01)

を小さく、つまり jakk
(k01)

選択するほうが、 aij

スケーリングきます。

を大きく実際、各方程式に 0 でない定数をかけても、数学の情報を生かすこ

的な意味は不変です。しかし、ピボット選択の結果jとができる。

は大きくちがってきます。
一般に、要素の大きさが極端に違うことは好ましくありません。ピボット選択を用いても見当違いのGauss の消去法における理論的な丸め誤差解析に答がでる例も簡単に作れます ([61])。

よれば、ピボット選択によって丸め誤差の増大が著

そこで、消去法に入る前に各要素の絶対値の大き

しく抑えられることがわかっています ([12])。

6.6

さがだいたい同じ程度になるように調整し、ピボッ
ト選択に意味を持たせる必要があります。
この操作を行列のスケーリング (scaling) またはピボット選択が必要ない行列

正規化(normalization)、または均等化 (equilibrium)

一般の行列に Gauss の消去法を適用する場合は、といいます。

前節にあげた理由よりピボットの部分選択または完
全選択を行う必要があります。
ただし、行列 A が以下のどれかの場合は、ピボット選択をしなくも LU 分解可能で、丸め誤差の影響も少ないことがわかっています ([12], [74])。


正定値対称行列

対角優位行列

 M

5
2.5
0
-2.5
-5

行列

20
15

行列 A が正定値であるかどうかを判定するには、

105

定義通り固有値を計算して、符号を調べればいいの10ですが、次数が大きくなるとなかなかそうもいきま

5
15

せん。

20

その場合、行列を導いた \意味" を考えるとわか

スケーリング前の行列ることがあります。
例えば、有限要素法などではある基底で張られた

現在のところスケーリングの決定版は存在しませ関数 uh のノルム式がん84 が、よく使われるものは次の方法です ([52])。
2
Th L2 ( ) = z Az

6.7.1 部分ピボット選択の場合

jjgrad u jjなどと 2 次形式の形で書き表せる場合があります。

i = 1; 1 1 1 ; n

もしノルムと基底の性質から、 z 6= 0 の場合 2 次
形式が正ということがわかるなら、 A が正定値行列に対し、 i 行の成分の絶対値最大
jaij ji = 1maxj nであることが判定できます。

83 厳密には「既約」という条件も必要です。
84 [13] にスケーリングによって失敗する人工的な例があります。また、 [7] も参照下さい。

324

でもって、ab

3 =a =
ij
ij

3 =b =iii;

j = 1; 1 1 1 ; ni

ピボット選択を必ずする



行列のスケーリングもあわせて行うと、一応の精度が期待できます。
以下、 Gauss の消去法のプログラムの例をあげまと操作します。式の内容は変わりません。

す。部分ピボット選択のプログラムは数値計算の本

この変形の結果、すべての係数が絶対値 1 以下とによく載っていますので ([45], [49] の解説がお勧めなり、少なくとも 1 つは絶対値 1 の係数がどこかにです)、ここでは省略し、あまりお目にかかることがある行列になります。

ない完全ピボット選択法の Gauss の消去法プログラ

6.7.2 完全ピボット選択の場合

ムを紹介します。
言語は Fortran で書きましたが、 FORTRAN77と Fortran90 の中間みたいになってしまいました85 。

完全ピボット選択の場合は、式だけでなく未知数の、すなわち列のバランスもとる必要があります。

なお、簡単なチューニングを次章で解説します。

そのためには、各行をその行の絶対値最大の要素

サブルーチン名は \Gaussian Elimination with

で割ったのちに、各列をその列の絶対値最大の要素

Complete Pivoting" の頭文字をとって、さらに、\Ver.0" の意味をつけて \GECP0" としておきます。

3jaij jj = 1maxin! =============================================!
完全ピボット選択 Gauss の消去法プログラム!
<<注意>> チューニング前なので遅いです！
! =============================================SUBROUTINE GECP0(A,NX,N,B,EPS,IW,P,Q,VP,&
VS,VW,ICON)DOUBLE PRECISION A(NX,1),B(1),VW(1),VS(1),&VP(1),D,EPSINTEGER N,NX,I,J,K,L,C,P(1),Q(1),IW,ICONICON=0IF(N.LT.1.OR.N.GT.NX.OR.EPS.LE.0.0D0)&
GO TO 8000
IF(IW.EQ.2) GO TO 6000
! ----------------- 各式の正規化 ---------------DO 10 I=1,N
D=0.0D0
DO 20 J=1,N
IF(ABS(A(I,J)).GT.D) D=ABS(A(I,J))
20
CONTINUE
IF(D.EQ.0.0D0) GO TO 7000
DO 30 J=1,N
A(I,J)=A(I,J)/D
30
CONTINUE
VP(I)=D
10
CONTINUE! ----------------- 未知数の正規化 -------------DO 40 J=1,N
D=0.0D0
DO 50 I=1,N
IF(ABS(A(I,J)).GT.D) D=ABS(A(I,J))
50
CONTINUE
IF(D.EQ.0.0D0) GO TO 7000
DO 60 I=1,N
A(I,J)=A(I,J)/D
60
CONTINUE
VS(J)=D40
CONTINUE
DO 110 K=1,N
P(K)=K
Q(K)=K110
CONTINUEでもってa

33 = a3 =

ij

ij

3 =

j xj

j

x

j;

i = 1; 1 1 1 ; nと操作します。 i の情報はあとで使うので記録しておく必要があります。

5
2.5
0
-2.5
-5

20
15
10

5
10

5
15
20

スケーリング後の行列
このようなスケーリングによって、新しい行列のどの成分も絶対値 1 以下となり、しかもどの行にも、また、どの列にも少なくとも 1 つの絶対値 1 の要素が含まれる結果となります。

6.8

プログラム例

以上のことから、行列が対称なのか正定値なのか

85 コメントの \!" は、最近の Fortran コンパイラならほとんど

優対角なのか性格が判然としない連立 1 次方程式に解釈してくれると思いますが、 \C" と同じです。

対応した Gauss の消去法のプログラムを書く場合、325

GECP0 では、ピボット選択の過程で生じる方程! ------------------ Pivot の捜索 --------------DO 120 K=1,N
D=0.0D0
DO 130 I=K,N
DO 130 J=K,N
IF(ABS(A(P(I),Q(J))).GT.D) THEN
D=ABS(A(P(I),Q(J)))
L=I
C=J
END IF
130
CONTINUE
IF(D.EQ.0.0D0) GO TO 7000
! ----------- index のつけ変え ------------------IF(K.NE.L) THEN
J=P(K)
P(K)=P(L)
P(L)=J
END IF
IF(K.NE.C) THEN
J=Q(K)
Q(K)=Q(C)
Q(C)=J
END IF! ----------- LU Decomposition ----------------IF(ABS(A(P(K),Q(K))).LT.EPS) GO TO 7000
A(P(K),Q(K))=1.0D0/A(P(K),Q(K))
DO 140 I=K+1,N
A(P(I),Q(K))=A(P(I),Q(K))*A(P(K),Q(K))
DO 150 J=K+1,N
A(P(I),Q(J))=A(P(I),Q(J))&-A(P(I),Q(K))*A(P(K),Q(J))
150
CONTINUE
140
CONTINUE
120
CONTINUE! ----------- Forward Substitution ------------6000 DO 160 I=1,N
VW(I)=B(P(I))/VP(P(I))
DO 170 J=1,I-1
VW(I)=VW(I)-A(P(I),Q(J))*VW(J)170
CONTINUE
160
CONTINUE! ---------- Backward Substitution ------------DO 180 I=N,1,-1
B(I)=VW(I)DO 190 J=I+1,N
B(I)=B(I)-A(P(I),Q(J))*B(J)190
CONTINUE
B(I)=B(I)*A(P(I),Q(I))
180
CONTINUE! ---------- 解の index つけ換え ----------------DO 200 I=1,N
VW(I)=B(I)
200
CONTINUE
DO 210 I=1,N
B(Q(I))=VW(I)/VS(Q(I))
210
CONTINUE
IC=0
RETURN! --------- ERROR が出た場合の処理 -------------7000 ICON=2000
RETURN
8000 ICON=3000
RETURN
END

式と変数の変換情報を記録する整数のインデックス
ベクトル
p = (p(1); p(2); 1 1 1 ; p(n))
q = (q (1); q (2); 1 1 1 ; q (n))

を用意し、この p(k) と q (k ) の値を交換することで正しい結果を得るようにプログラムを組んでいます。
分解のプログラムでは、 i の情報をそのま
ま部分ピボット選択に利用する方法を採用する場合
LUが多いのですが (具体的なプログラムは [49] 参照)、
LU

分解のコストに比べスケーリングのコストは大したことないことと、プログラムを見やすくするた
め、陽に行列をスケーリングしています。
次の章では、サブルーチン GECP0 をできるだけ楽をしてスピードアップする方法を考えます。

なお、アルゴリズムは [60] を参考にしました。

ちょっと休憩

326

7 GECP のスピードアップ作戦
この章では、完全ピボット選択付き Gauss の消去

2. ピボット選択と行・列の交換法プログラム \GECP0" を、スカラー計算機、およ

3. LU 分解びベクトル計算機でとりあえずのスピードが出るよ4. 前進代入／後退代入、および解の並べ換え

うにチューニングします。
計算機の特性にあわせたチューニングのためにすぐれたプログラミング言語には、実行時の振舞

は、計算機内部の作りに立ち入った知識が多少必要

いを解析するソフトウェアが付属しています87 。

になります。しかし、最近のプログラム言語の最適

しかし、ここではもっと具体的に、プログラム内化技術は優秀なので、ある程度の筋道を示すと、あに CPU 時間を計測する関数を埋め込んで、 4 つのとは勝手に仕事をしてくれます86

処理の経過時間を測定します。

ここでは、富士通株式会社のコンパイラ部隊が測定は 1995 年 10 月 6 日におこないました。計算総力を結集して作り上げた Fortran コンパイラの機は FUJITSU VP2600/10、コンパイラは FOR-

\FORTRAN77 EX/VP" を用いて、 GECP0 のス

TRAN77 EX/VP です。最適化オプションはOPT(E)ピードアップ作戦を展開します。

7.1

として、スカラーモードとベクトルモードで方程式の次数を変えてデータを集めました88 。
以下は、 GECP0 における各処理過程の全体の実実行時間の調査

行時間に占める割合です。まずはスカラーユニット

まず [45] の指摘を紹介します。

のみを使った実行結果です。

経験的事実として、 CPU 時間が全体の中で支配的なプログラムでは、 CPU 時間の大部分 (90% 以上) は
ソースプログラムの数 % 以下の部分で集中して消費される場合が多い。したがって、プログラムの改良を
効率よく行なうためには、計算時間の長いプログラムの部分を検出することが重要である。

各処理に要するコスト (スカラーモード)
処理過程
n = 100
n = 1000
行列の正規化

2.64%

0.23%

ピボット選択処理

49.53%

38.25%

分解

46.39%

61.40%

前進消去／後退代入

1.44%

0.12%

LU

島崎 眞昭『スーパーコンピュータとプログラミング』

データからすぐわかるように、ピボット選択とプログラム全体で考えると、連立 1 次方程式や固LU分解にかかるコストが他を圧倒しています。

有値計算などが「CPU 時間の大部分」を消費する
各処理に要するコスト (ベクトルモード)
処理過程
n = 100
n = 1000

ことはよく知られています。
つまり、 GECP0 の計算が遅いと、全体の計算時間が遅くなってしまいます。その意味で GECP0 のチューニングは非常に大事になります。
第一にやるべきことは、 [45] の指摘通り、

行列の正規化

0.66%

0.00%

ピボット選択処理

20.57%

0.80%

分解78.03%

98.20%

前進消去／後退代入

0.73%

0.00%

LU

各処理過程で最もコストがかかる部分を特定する

対してベクトルモードでは、 LU 分解に大半のコスト (n = 1000 で 98% 以上) が集中しています。

ことです。

これは極めて重要なことです。なぜこのような違GECP0 は 4 つの部分に分けることができます。

いが生じるのでしょうか？
1. 正規化 (スケーリング) 処理

87 FORTRAN77 EX/VP には ANALYZER というツールが86 こういった優秀なのが研究室にたくさんいればいいのですあります。
88 このあたりの用語がわからない方は、 [63] をごらん下さい。

が、現実はそうもいきません。

327

ベクトル化できない原因
n = 1000

P(I) と Q(I) は 1 から n までのどれかの値をとり、同じ値となることは絶対にないでの実行時間を見ると、

ピボット選択処理
LU

分解

スカラー

ベクトル

ことを知っています。よって、I を K+1 から N まで

149 秒

1.9 秒

「せーの」でベクトル処理して欲しいと期待するの

239 秒

237 秒

は当然です。
ところが、コンパイラはそんなことはつゆ知りまです。ベクトルモードでは、ピボット選択処理の実せんので行時間が 149 秒から 1.9 秒と大幅に短縮できています。対して LU 分解の方は全く改善されていませ「もしかしたら P(I) は 1, 2, 2, 1,... な
どのように、重なりがある回帰データかも知れないな。もしそうなら、ベクトル処理によってA が変な値に書きかわってしまうかもしれないぞ。これはまずい。スカラーユニットで順番に処理しようっと。」

ん。これは、 LU 分解の部分が全くベクトル化できてないことを意味します。ベクトル計算機にとって、この事態は極めて遺憾なことです。
従って、ベクトルモードの処理において LU 分解に大量のコストがかかる理由がわかります。つまり、他の処理がそこそこベクトル化できているのにと考え、この部分をスカラーモードで実行します。

対して、 LU 分解の部分だけ全くベクトル化できてこれは、コンパイラにとって当然の選択です。

いないため、相対的なコストが増大したわけです。

その結果、本来はベクトルユニットで実行されるでは、この部分の改良を考えましょう。

7.2

べき DO ループが非常に遅くなってしまいます91
では、どうすべきか？ 対策は次のどちらかです。

行列の値の交換《対策 1》 ベクトル化を推進する最適化指示行を

7.2.1 親切なコンパイラ

ソースプログラムに埋め込み、配列に回帰参照
がないことをコンパイラに指示する。

GECP0 のプログラムでピボット選択によって生じる行と列の交換は、行の情報を記憶するインデッ《対策 2》 おしゃれなインデックスベクトルに頼る

クスベクトルの値を交換することですませていまのをやめ、ピボット選択時に本当に行と列の値す。

を交換することでコンパイラに自分の意図を思い知らせる。

このようにすれば、実際に行列の値を変換をする

7.2.2 最適化指示行の挿入手間が省け、計算時間の節約になる、と思うのが普通です。インデックスベクトルを添字に用いる方法
は George E. Forsythe さんと Cleve B. Moler さん《対策 1》は、 DO ループの前に「この配列にはの書いた有名な数値計算の教科書 [12] などでも推奨

回帰参照はないよ」という意味の行を 1 つ追加しま

されています89 。が、このインデックスベクトルをす。 FORTRAN77 EX/VP の場合、次のようにな使ったプログラムは、ベクトル計算機に適さないこります。

とが知られています。
*VOCL LOOP,NOVREC(A) ! この行を追加
DO 100 I=K+1,N
A(P(I),Q(K))=A(P(I),Q(K))*A(P(K),Q(K))
100
CONTINUE

原因は、プログラムを組む人のベクトルコンパイ
ラに対する説明不足にあります90 。例えば、

140

DO 140 I=K+1,N
A(P(I),Q(K))=A(P(I),Q(K))*A(P(K),Q(K))
CONTINUENOVREC は \no vector recurrence" の略です。詳しくは [68] を御覧下さい。

の部分をコンパイラが見た場合、添字の P(I) と

ただしこの方法は、ベクトル計算機とコンパイラ

Q(I) の値は実行するまでわからないと判断します。
実際、ピボット選択が済むまで P(I) とQ(I) の値の種類によって異なりますので、あくまでも個別撃破の作戦となります92 。

はわかりません。しかし、プログラムを書く人は、

91 [63] の実験でわかるように、ベクトル計算機からベクトルユ

89 もっとも、これは部分ピボット選択のプログラムに関するこニットをとれば、汎用機どころかパソコンにも負ける可能性があとで、この場合 A(P(I),J) に対しての話です。完全ピボット選

ります。
92 また、この例でのベクトル化の効果は《対策 2》と比べて半択の場合 A(P(I),Q(J)) なので状況は異なります。
90 あるいはコンパイラの認識不足ともいえるのですが、それは分以下です。

酷というものです。

328

7.2.3 本当に行と列の値を交換する

コードをはきだしたため、スカラーモードに比べて何と 498 倍のスピードがでました。こうでなくちゃ

汎用性を追求するならば、対策 2 です。さっそくいけません。

プログラムを修正します。

また、行と列を交換したプログラムの方が、オリDO 120 K=1,N!-------------- Pivot の捜索 ------------------ !
D=0.0D0
DO 130 I=K,N
DO 130 J=K,N
IF(ABS(A(I,J)).GT.D) THEN
! 変更 1
D=ABS(A(I,J))! 変更 2L=I
C=J
END IF
130
CONTINUE!------------- index のつけ変え ---------------- !
IF(K.NE.L) THEN
J=P(K)
P(K)=P(L)
P(L)=JDO 131 J=1,N! 追加 1
D=A(K,J)! 追加 2A(K,J)=A(L,J)! 追加 3A(L,J)=D! 追加 4131
CONTINUE! 追加 5END IF
IF(K.NE.C) THEN
J=Q(K)
Q(K)=Q(C)
Q(C)=JDO 132 I=1,N! 追加 6
D=A(I,K)! 追加 7A(I,K)=A(I,C)! 追加 8A(I,C)=D! 追加 9132
CONTINUE! 追加 10END IF

ジナルのインデックスベクトルを使ったプログラム
よりスカラーモードでも速くなっています。行と列の交換部分を新規に追加するため、その分計算時間が余計にかかるのでは、と心配されたピボット選択
処理部分も、 37 秒速くなりました。その理由は、


ピボット選択の実行時間の大部分は、最大ピボット探索の DO ループが消費する



その DO ループ内で A(P(I),Q(J)) が A(I,J)とスッキリしたため、実行効率があがった

ためです。

7.3

添字の動き

続いては、 Fortran 言語の特徴を利用します。次の文をお読み下さい ([60])。
一般に、 FORTRAN では、 2 次元以上の配列の配列
要素は、実際のメモリ上では、第 1 の添字の順に連続して並び、次に第 2 の添字が一つ増えて再び第 1 の添
字の順に連続して並ぶようにコンパイルされることになっている。（中略）したがって、たとえば 2 次元配列の添字のそれぞれを二重の DO ループで動かすような場合には内側の DO ループで第 1 の添字が動き、外側の DO ループで第 2 の添字が動くようにプログラムを書くべきである。

以上はピボット選択の部分です。
\追加" の部分で行と列の入れ換えを行なうように

森 正武『数値解析法』

しました。なお、入れ換えの情報 P, Q は、前進／後退代入で使うのでそのまま残します。

これに従ってプログラムを修正します。まず、ピ! -------- LU Decomposition ---------------- !
A(K,K)=1.0D0/A(K,K)
DO 140 I=K+1,N
A(I,K)=A(I,K)*A(K,K)
DO 150 J=K+1,N
A(I,J)=A(I,J)-A(I,K)*A(K,J)
150
CONTINUE
140
CONTINUE

ボット捜索の部分の以上は LU 分解の部分です。インデックスベクト

130

ル P と Q が追放され、すっきりしました。

の I と J を入れ換え

さっそく実行時間を計測します。さきほどと同じ
n = 1000でのデータです。

ピボット選択処理
LU

分解

スカラー

ベクトル

112 秒

0.88 秒

224 秒

0.45 秒

DO 130 I=K,N
DO 130 J=K,N
IF(ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
L=I
C=J
END IFCONTINUE

うまくベクトル化できました。コンパイラは今

130

度は LU 分解をすべてベクトルユニットで処理する
329

DO 130 J=K,N!
<--- 修正DO 130 I=K,N!
<--- 修正
IF(ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
L=I
C=J
END IFCONTINUEとします。次に、 LU 分解の部分

150
140

! =========================================== !
! ------------------------------------------- !
!
SUBROUTINE GECP(A,NX,N,B,EPS,ISW,P,Q,VP,&
VS,VW,ICON)! ------------------------------------------- !
! <<概要>>!
実行列の連立 1 次方程式 Ax=b を complete pivoting!
を用いた Gauss の消去法で解く倍精度用サブルーチン! <<引数>>!
A : 係数行列．演算後, 内容は保存されない．
!
NX : 配列 A の整合寸法!
N : 係数行列 A の次数 (N>1)!
B : 定数ベクトル B．解ベクトル x として出力!
EPS : 係数行列の特異性の判定定数!
ISW : ISW=1 のとき -- 最初から処理を行う!
ISW=2 のとき -- LU 分解のデータはそのまま!
B の値だけを換えて実行!
P : pivoting の列の入れ換え情報ベクトル!
Q : pivoting の行の入れ換え情報ベクトル!
VP : 式のスケーリング情報ベクトル!
VS : 方程式のスケーリング情報ベクトル!
VW : 作業用実ベクトル! ICON : condition code!
0 -- 正常終了!
2000 -- 係数行列のある行か列がゼロ!
または pivot 判定が EPS 以下!
3000 -- 引数が間違っている! ------------------------------------------- !
IMPLICIT NONEDOUBLE PRECISION A(NX,1),B(1),VW(1),VS(1),&VP(1),D,EPSINTEGER N,NX,I,J,K,L,C,P(1),Q(1),ISW!
! ------------------------------------------- !
引数のチェック!
! ------------------------------------------- !
ICON=0IF(N.LT.1.OR.N.GT.NX.OR.EPS.LE.0.0D0)&
GO TO 8000
IF(ISW.EQ.2) GO TO 6000!
! =========================================== !
!
スケーリング (正規化処理)! ------------------------------------------- !
! * 処理内容 *!
Pivoting の前段階として，各係数の大きさがほぼ!
等しくなるように，個々の式をその式の絶対値最大の!
係数で割り調整する．
!
さらに各未知数間のバランスをとるため，係数行列の!
各要素をその列の絶対値最大の係数で割る.
! * 注意 *式のスケーリング情報を VP に保存．また,!
!
未知数の情報を VS に保存する! =========================================== !
!
! ------------- 各式の正規化 ------------------ !
DO 10 I=1,N
D=0.0D0
DO 20 J=1,N
IF(ABS(A(I,J)).GT.D) D=ABS(A(I,J))
20
CONTINUE
IF(D.EQ.0.0D0) GO TO 7000 ! 行列が非正則
DO 30 J=1,N
A(I,J)=A(I,J)/D
30
CONTINUE
VP(I)=D
10
CONTINUE

DO 140 I=K+1,N
A(I,K)=A(I,K)*A(K,K)
DO 150 J=K+1,N
A(I,J)=A(I,J)-A(I,K)*A(K,J)
CONTINUE
CONTINUE

で二重ループの内側で第 1 添字の I が動くように

140

150
145

DO 140 I=K+1,N
A(I,K)=A(I,K)*A(K,K)! <-- ループ分割
CONTINUE! <-- 修正DO 145 J=K+1,N! <-- 修正DO 150 I=K+1,N! <-- 修正A(I,J)=A(I,J)-A(I,K)*A(K,J)
CONTINUE
CONTINUE! <-- 修正と修正します。その他、前進／後退代入計算も修
正することができますが、全体に占めるコストが微々たるものなので、何もしないことにします。
実行結果を見ましょう。

ピボット選択処理
LU

分解

スカラー

ベクトル

38 秒

0.87 秒

29 秒

0.45 秒

ベクトルモードで実行時間が変わらないのは、既に自動ベクトルコンパイラによってチューニングが完了しているためです ([45], [49])。
一方、スカラー処理では、たったこれだけの修正なのに、ピボット選択で約 3 倍、 LU 分解で 7.7 倍もの性能がでました。以上 2 つの改良方法は、ベク
トル計算機だけでなく、汎用機やワークステーション、パーソナルコンピュータなどにも共通する技法です。

7.4 GECP プログラム
プログラムのチューニングの余地はまだまだあります。例えば、前進／後退代入とスケーリングの部
分はまったく手つかずです。また \ループアンロー
リング" という有名な手法をとり入れるのも効果的でしょう。さらに、記憶領域の節約の面では、引数で使うスケーリング情報を格納するベクトルは、宣言を不要にすることができます。
しかし、最低限のチューニングをするだけでも、(もとのプログラムに比べて) 格段の性能アップができましたので、一応の完成にします。
出来上がりの完全ピボット選択サブルーチン名を\GECP" とします。

以下はプログラムリストです。

330

! =========================================== !
!
LU Decomposition! ------------------------------------------- !
IF(ABS(A(K,K)).LT.EPS) GO TO 7000
A(K,K)=1.0D0/A(K,K)
DO 140 I=K+1,N
A(I,K)=A(I,K)*A(K,K)
140
CONTINUE
DO 145 J=K+1,N
DO 150 I=K+1,N
A(I,J)=A(I,J)-A(I,K)*A(K,J)
150
CONTINUE
145
CONTINUE
120 CONTINUE!
! =========================================== !
!
Forward Substitution! ------------------------------------------- !
6000 DO 160 I=1,N! ISW=2 の場合, ここから
VW(I)=B(P(I))/VP(P(I))
DO 170 J=1,I-1
VW(I)=VW(I)-A(I,J)*VW(J)170
CONTINUE
160 CONTINUE!
! =========================================== !
!
Backward Substitution! ------------------------------------------- !
DO 180 I=N,1,-1
B(I)=VW(I)DO 190 J=I+1,N
B(I)=B(I)-A(I,J)*B(J)190
CONTINUE
B(I)=B(I)*A(I,I)
180 CONTINUE!
! =========================================== !
解の index つけ換え!
! ------------------------------------------- !
DO 200 I=1,N
VW(I)=B(I)200 CONTINUE!
DO 210 I=1,N
B(Q(I))=VW(I)/VS(Q(I))
210 CONTINUEICON=0! ここまで来れば正常終了
RETURN!
! =========================================== !
!
ERROR が出た場合の処理! ------------------------------------------- !
7000 ICON=2000
RETURN
8000 ICON=3000
RETURN
END! ------------- 未知数の正規化 ---------------- !
DO 40 J=1,N
D=0.0D0
DO 50 I=1,N
IF(ABS(A(I,J)).GT.D) D=ABS(A(I,J))
50
CONTINUE
IF(D.EQ.0.0D0) GO TO 7000 ! 行列が非正則
DO 60 I=1,N
A(I,J)=A(I,J)/D
60
CONTINUE
VS(J)=D40
CONTINUE!
! =========================================== !
!
Complete Pivoting! ------------------------------------------- !
! * 処理内容 *!
LU 分解の過程の小行列で絶対値最大の要素を特定!
対応する行と列の番号をベクトル P, Q に保存! =========================================== !
DO 110 K=1,N
P(K)=K! 初期値の設定
Q(K)=K
110 CONTINUE! ------------- pivot の捜索 ----------------- !
DO 120 K=1,N
D=0.0D0
DO 130 J=K,N
DO 130 I=K,N
IF(ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
L=I! pivot となる行番号
C=J! pivot となる列番号
END IF
130
CONTINUE
IF(D.EQ.0.0D0) GO TO 7000 ! 行列が非正則! ----------- index のつけ変え ---------------- !
IF(K.NE.L) THEN
DO 131 J=1,N
D=A(K,J)
A(K,J)=A(L,J)! 行の入れ換えA(L,J)=D131
CONTINUE
J=P(K)
P(K)=P(L)! 行情報の入れ換え
P(L)=J
END IF!
IF(K.NE.C) THEN
DO 132 I=1,N
D=A(I,K)
A(I,K)=A(I,C)! 列の入れ換えA(I,C)=D
132
CONTINUE
J=Q(K)
Q(K)=Q(C)! 列情報の入れ換え
Q(C)=J
END IF

331

では、 GECP サブルーチンの処理速度を他の連

7.5

立 1 次方程式のサブルーチンと比較します。計算は実行時間の比較

ベクトル計算機 VP2600/10 で行ないました。連立完成したプログラム \GECP" で各処理がどれく

1 次方程式のサブルーチンは次の 4 つです。

らいのコストを占めるかをベクトル・スカラーモー
ド別に表にします。
サブルーチン名
各処理に要するコスト (スカラーモード)
処理過程
n = 100
n = 1000

Gauss の消去法 (部分ピボット選択)

LEQLUW

Crout 法 (部分ピボット選択)
Householder 変換による QR 分解
Gauss の消去法 (完全ピボット選択)

行列の正規化

4.30%

1.33%

DLAXL

ピボット選択処理

56.74%

55.88%

GECP

分解

37.14%

42.29%

前進消去／後退代入

1.82%

0.52%

LU

解法

DVLAX

DVLAX と DLAXL は SSL II のサブルーチンです。
DLAXL は、もともと最小自乗法用なのですが、行と列の数を揃えると連立 1 次方程式の解法としても相変わらずピボット処理と LU 分解の部分にコストが集中しています。うち、ピボット操作の実行時使えます。 LEQLUW は NUMPAC のサブルーチン

間が半分以上を占めています。

です。
この中で VP2600/10 の性能を最高に引き出すサ

各処理に要するコスト (ベクトルモード)
処理過程
n = 100
n = 1000

ブルーチンは DVLAX です。 [63] での調査によれ
ば、 VP2600/10 のピーク性能に迫る数値を出して行列の正規化

2.61%

0.63%

ピボット選択処理

64.50%

65.40%

よって、 DVLAX の 処理時 間を 1 とし て、各 サ

分解30.94%

33.59%

ブルーチンがその何倍の処理時間になるか、次数を前進消去／後退代入

1.94%

0.38%

変えながら測定しました。

LUいるので、他はとてもたちうち出来ません。

ベクトルモードでは、ピボット処理にかかる割合
がさらに増加して 60% をこえています。さらに詳し
ratio to DVLAX

く調べると、ピボット選択処理の 95% 以上が130

LEQLUW
DLAXL
GECP

100

DO 130 J=K,N
DO 130 I=K,N
IF(ABS(A(I,J)).GT.D) THEN
D=ABS(A(I,J))
L=I
C=J
END IFCONTINUE

10

1
0

1000

2000

3000
4000
matrix dimension

5000

6000

7000

サブルーチンの処理性能の部分に集中していることがわかります。コンパ
イラの情報を見ると、この部分は既にベクトル化さ

QR

分解の DLAXL はもともと演算数が LU 分解

れています。従って、数倍の order はあったとしてに比べて多いことと、ベクトル計算機にうまく対応

も、数十倍、数百倍の order での改善は難しいと思していないこともあって、 N = 4000 で DVLAX

われます93 。

の 300 倍以上の実行時間がかかってしまいました。

この部分のコストは、部分ピボット選択ではほと

LEQLUW は N

んど無視できる値です。つまり、どんなに頑張って= 2000 までほぼ 2 倍の処理速度でしたが、 N が増えるとだんだん性能が低下して、

GECP をチューニングしても、絶対に部分ピボット選択法より速くなることはありません94 。演算回数N
N=

7000 では 2.8 倍になりました。 GECP は= 500 のときに最大 7.6 倍となりましたが、 Nを比較しても、計算時間が 2 倍程度なら素晴らしいが大きくなるとじわじわと差を詰め、 N = 7000 で速さだといえます。

は 6 倍の実行時間となりました。

93 はっきり言って、この部分のチューニングをどうやっていい

LEQLUW と GECP の比較では、目標の \2 倍"のかよくわかりません。ループアンローリングくらいでしょう

から見るとまあまあの結果です。ただし DVLAX のか？
94 部分ピボット選択のプログラムも同様にチューニングした場

レベルから見れば、まだまだと言えるでしょう。

合の比較です。

332

8 Gauss の消去法は安定か？
そのため、紙と鉛筆で連立 1 次方程式を解く場合

最後の章は、電子計算機が世の中に登場して以の有力な手法として、 Gauss の消去法は揺るぎない

来、ずっと議論の種になっている古典的命題地位を獲得しました。
部分ピボット選択の Gauss の消去法は、はたし1940 年代になると、電子計算機が本格的な実用段て安定な数値解法か？

階に入り、それにともなって、丸め誤差の研究もさかんになります。当然、連立 1 次方程式の簡単な解法である Gauss の消去法も格好の研究対象となりまをとりあげたいと思います。

した。

この問題は、理論的には既に解決済みです。ただ最初になされた研究では、丸め誤差が Gauss の消し、経験的・実用的には \もやもや" としたものが去法にあたえる影響はきわめて悲観的でした。 [30]

残ります。

8.1

から引用します。
電子計算機時代の幕開けの頃、 Gauss の消去法は連立 1 次方程式の解法として役に立たないのでは、と懸念された。 Hotelling は 1943 年の論文で、 AT Ax =b の形をした n n 方程式の解において誤差は 4n01程度になるであろうと予言した。彼はこう指摘する。
「78 行の行列では、小数点第 1 位の近似精度を保証するのですら、 46 桁が必要となるだろう。」

消去法の歴史Gauss さんは、 1809 年に出版された『天体運動2論』という本の中で消去法を紹介し ([50])、この方法は広く世に知られるところとなりました。
しかし、 Gauss さんが消去法のアルゴリズムの発明者というわけではないようです。

L. N. Trefethen and R. S. Schreiber

18 世紀後半にフランスで活躍した数学者 Joseph
Louis Lagrange さん (1736{1813)95 は、 Gauss さ\Average-Case Stability of Gaussian Elimination"

んが『天体運動論』を出す数十年前に、斉次な系を最後の \桁数云々" の部分は、行列の次数を n =解くための消去法のアイデアを既にほのめかしてい78 として 10 進数の計算機を用意したとき、誤差が4n01 ですから、ました ([17])。
また、次の事実もあるそうです。

x = 47710

遠く中国では『九章算術』の第 8 章「方程」に、ガウスの消去法とまったく同じアルゴリズムの記述がある。『九章算術』は紀元前 2 世紀の頃に編纂されたというから、すでに 2000 年も前から知られていたわけである。

の解として、桁数が 46 以上絶対に必要なことを意味します。
当初の研究での結論は、 100 個程度の未知数を持
つ方程式系を解く場合でも、不可能なほどの多くの

戸川 隼人『数値計算』

有効桁を保って演算しない限り、必要な精度はとても得られまい、というものでした。

みなさんも、高校や大学で問題を解いた96 ことがそのため、一時期 Gauss の消去法は計算機を用いあると思いますが、消去法や掃き出し法は、理論がた数値計算の舞台から姿を消し、しばらく反復法の明解な上、四則演算さえ知っていれば、あとは機械天下になります97 。

的な行列計算で答が出てきます。

97 「半分眠りながらでも解けるよ」と弟子への手紙で反復法を

95 変分法、解析力学の創始者です。また、代数の世界でも偉大勧めていた Gauss さんも、生きていれば喜ばれたでしょう。消去

な足跡を残しました。メートル法の制定における委員長でもあっ法に \Gauss の" という言葉がくっついていることを当人が知った人です。
96 正確には解こうと努力したていたかどうかよく解りません。

333

8.2.2 条件数

浮動小数点演算と固定小数点演算98 における誤差解析を理論的に発展させ、 Gauss の消去法を数値計

誤差を測るためには、何らかの形の「距離」の概

算の表舞台に復活させた人たちは、 Alan Mathison

念が必要です。その尺度をノルム(norm) といいま

Turing さん (1912{1954)、 John von Neumann さす。

ん (1903{1957) と Herman Heine Goldstine さん

ベクトルのノルムを(1903{?)、そしてとりわけ James Hardy Wilkinson

さん (1919{1986) です ([36])。

8.2

ijjxjj1 = max jx ji

また、行列のノルムを後退誤差解析

Xn

jjAjj1 = max

Wilkinson さんの誤差評価手法は、後退誤差解析

i(backward error analysis)といわれる手法です。読

ja

j =1

ij j

で定義します101 。

んで字のごとく、解をたどって定式化を行ない、誤差を見積もる方法です。

ノルムの性質をより、直ちに次のことが言えま

後退誤差解析の方法と成果は [35], [36], [37] に詳す。

細にまとめられています99 。これに対して、前進誤

条件数

差解析 (forward error analysis)もあります。こちらA

は [27], [28] を御覧下さい。

が 1A だけ変化したとき x が 1x だけ変わったとすると後退誤差解析より得られた知見の中から、 Gaussの消去法に関する誤差評価を紹介します。まずは、

jj1xjj1

用語の説明からはじめます。

01

jjx + 1xjj1

8.2.1 マシンイプシロン

 jjAjj1 jjA

jj1

jj1Ajj1
jjAjj1をみたす. また, b が 1b だけ変化したときは
jj1xjj1

マシンイプシロン (machine epsilon)は、浮動小jjxjj1数点数の相対誤差の上限をあらわす値で、計算機と浮動小数点体系によって決定されます ([60], [54])100 。

01

 jjAjj1 jjA

jj1

jj1bjj1
jjbjj1となる.

マシンイプシロンを "M とすれば、
"

です。

M =c

この二つの不等式は、 x の変化が A や b の変化の

10t

(14)

何倍に拡大されるかをあらわしています。
この拡大を表す値を

は基数、 t は仮数部の桁数、 c は丸めの方法によって決まる定数で、四捨五入の場合 1=2、切1 (A) = jjAjj1 jjA

01

(15)

jj1

捨ての場合 1 となります。
で書き、 1 (A) を行列 A の条件数 (condition num-

現在よく使われている数値表現での "M の値は以ber)

下の通りです ([55], [67])。

といいます。

一般に、条件数が大きくなると、係数のわずかなt

IBM 形式 単精度

16

6

倍精度

16

14

4 倍精度

16

28

IEEE 形式 単精度

2

24

"

M

変化で解の変化に大きな影響がでます ([58])。

9:54 2 1007
2:22 2 10016

条件数がものすごく大きくなるような行列は、相対誤差が拡大する可能性が高く、不安定性を持ちま3:10 2 10033
5:96 2 1008

す。このような行列を悪条件 (ill-conditioned) な行
列といいます。

倍精度
253
1:11 2 10016《注意》 "M は 10 進数に換算した近似値

実際の計算で jjA01 jj1 の値を少ない計算量で正確に見積もるのはなかなか難しい問題で、感度解析とよばれる分野の中心問題です。詳しくは [55] をM は、反復法の収束判定などによく利用される"

御覧下さい。

基本的な数値です。

101 有限次元空間の場合、 \ノルムの定義" をみたすノルムはすべ

98 現在はあまり使われなくなりました。詳しくは [36] を御覧下て同じと考えることができます。ノルムには他にもさい。
99 今はもう古典の部類に属してるといえるかも知れません。

などがあります (詳しくは [32], [60] 参照)。

100 「計算機イプシロン」ともいいます。

334

jj1jj1 や jj1jj2

8.2.3 誤差評価

8.2.5 Hilbert 行列
悪条件となる行列の有名な例として、 Hilbert

以上の数値を用いて誤差評価を行ないます。
x^を浮動小数点演算による Gauss の消去法を用い行列があります。 Hilbert 行列は、 David Hilbert
さん103 (1862{1943) が 1894 年に発表した論文よりて求めた近似解、また  を =(k)maxi;j;k jaij j
max i; jjaij j

名付けられたもので、(16)で定義します。このとき Wilkinson さんはいくつ

有名な行列のわりには、どうやって導かれたのか([37])。

数値解析の本にはあまり書いてないので、経緯をWilkinson の定理[12] から簡単に説明します。
f (x)き、その相対誤差は次で評価される。

jjxjj1

 4n 1 (A) "(17)

E (c) = jjf

n

i

c x

i=1

i01

n (x) 0 f (x)jj 2
L( )を最小にするように係数 ci を決定するものです。なおいては不可避の値です。しかし、倍精度演算をすお、 L2 ノルムとはれば、マシンイプシロンの頑張りでなんとか小さくなりそうです102

jjg (x)jj

結局、 Wilkinson さんの与えた定理より次のことがわかります。

2

L2 ( ) =Z

10g (x )

2dx

です。 E (c) の微分は 0 となるので、 ci でひとつず

X Z

 Z

つ微分して変形すれば、Gauss の消去法は、条件数と  がめちゃくちゃn

大きくならない限り安定である。

0

j =1

tor)

n (x) =の形で近似したときの L2 ノルムの誤差

また、 " も計算機と精度により決まります。
4n2 は次数 n によって決まり、浮動小数点演算にX

を n 0 1 次までの多項式でf

条件数 1 (A) は行列 A が決まれば固定されます。

= [0; 1] で定義される連続関数

問題は、区間連立 1 次方程式を Gauss の消去法で求めたと

2

(18)

i+ j 01となる行列のことです。

かの仮定のもと、次の有名な定理を証明しました

jj x^ 0 xjj1

1

ij =

a

1

i+j 02 dx
x

c

j=

1

0

i01 dxf (x)xを得ます。

は、三角行列化における不安定因子(growth fac-

これは、 n 個の未知数に対する連立 1 次方程式に、または増加因子と呼ばれます ([38])。

Z

なっています。行列 A の成分は8.2.4 悪条件の場合の対策a

条件数 1 (A) がどれくらいの大きさになれば \悪

ij =

1

0

x

i+j 02 dx =1
i +j 0 1となって、 (18) が出てきます。

条件" と言うか、基準はいろいろです。少なくともnn に関して指数関数的な勢いで、たとえば 2 とかne とかで増加してしまうと、もうお手上げです。

Hilbert 行列は、対称かつ正定値行列なのですが、その条件数は指数関数的に増大します (詳しい値が知りたい方は [12], [18] を参照下さい)。

[55] から引用します。

従って、ほんのわずかな誤差の混入が著しく解の良条件の問題は、誤差が増大しない安定な方法を用いれば精度よく解くことができる。悪条件の問題を精度よく解くためには、いくら解法を工夫してもだめである。唯一の対策は計算桁数を多くすることである。

精度を損ないます。対策としては、計算桁数を多くとるしかありません。
Hilbert 行列は、そのあまりの条件数の悪さ、即ち解きにくさから、 1940 〜 50 年代にかけて盛んに

名取 亮『数値解析とその応用』研究されました104
103 20 世紀前半の最も偉大な数学者の一人です。 1900 年パリの102 逆にこの定理は、単精度演算で n = 10000 や n = 100000

国際数学者会議における『23 の数学の問題』は特に有名です。
104 n = 10 までの逆行列を与えるだけの論文というのもありまなどの計算をすることの無謀さを語っています。ご注意下さい。

す。

335

ただし、 Hilbert 行列の構造を研究するのも大事

これによって、 A が悪条件でない場合、対称正定なのですが、 George E. Forsythe さんの次の指摘も値行列106 と優対角行列はピボット選択をしなくても大事です ([12])。

安定、また Hessenberg 行列でも部分ピボット選択
によってほぼ安定といえます。

Hilbert 行列の条件の悪さは、この行列の導入に用いた近似の問題にさかのぼることができる。区間 0 x  1 の上での関数 xi (i = 0; 1 1 1 ; n 0 1) はほとんど
1 次従属である。このことは Hilbert 行列の行ベクト
ルがほとんど 1 次従属であること、すなわちこの行列

また、 Wilkinson さんによれば、完全ピボット選
択の上限 (20) は極めて「ゆるい」限界であって、証明した本人も「ひどい過大評価だ」と認めています。実数行列において今まで知られている最大の 
は n です。任意の正則行列に対して   n が成り立がほとんど特異であることを意味する。（中略）要す
るに、この近似の問題は、 Hilbert 行列のような行列を導く形をしている以上、問題の \たて方が悪い"。

つのではないか、という Wilkinson さんの有名な予想があり、特定の n では証明ができていますが、一G. E. Forsythe and C. B. Moler

般の証明はまだのようです ([23], [30])。

『計算機のための線形計算の基礎』

ここでは、その予想を信用することにします。
完全ピボット選択の Gauss の消去法は、条件数8.3が指数関数的に増大しない限り安定不安定因子の上限

Gauss の消去法に おける Wilkinson さんの 定理

問題は、部分ピボット選択の上限 2n01 です。 (19)(17) のもう一つの主役は、不安定因子  です。

と (20) を比べて見ましょう。

条件数が純粋に行列の性格、つまり問題の設定そのものに起因するのに比べ、不安定因子は行列の性格にプラスして、算法の方法によって大きく変わります。
Wilkinson さんは、誤差定理を証明するだけでなく、さらに [35] において、  の理論的な上限値を与partial pivoting

512524; 28875
5:6 2 1014

8.4
(19)

Wilkinson さんの例題

部分ピボット選択を用いた Gauss の消去法における不安定因子の上限 2n01 を証明した Wilkinson さ

完全ピボット選択を用いたとき



50

15

になりそうです。

の上限 (Wilkinson)

部分ピボット選択を用いたとき

 

20

い勢いで増加します。これは、ちょっとまずいことると仮定します。

n01
  2

10
6:1

部分ピボット選択の不安定因子の上限は、ものすご

えました105 。以下、方程式はスケーリングされていncomplete pivoting

んは、同じ論文で正確な上限を与える行列を披露し1=2 (21 31=2 41=3 1 1 1 n1=(n01) )1=2
nています ([35])。

1:8n1=4 log n

n

(20)

次行列を An とします。行列の作り方は以下の通り、いたって簡単です。

対称正定値行列のとき (ピボット選択なし)

k = 1; 2; 1 1 1 ; n 0 1

  1

対角優位行列のとき (ピボット選択なし)
  2

に対し、 An の第 k 列を

上から k 0 1 番目までが 0

そのあとに 1; 1; 01; 1; 01; 1; 1 1 1 が続く



最後の列だけは 01; 1; 01; 1 1 1として作成します。 A の条件数はそれほど大きくあ

Hessenberg 行列に部分ピボット選択を用いた

りません ([35])。従って、解の安定性は不安定因子

とき



  n

の大きさにかかってきます。

106 非対称正定値でも  1 が示せますが、その場合、 L の成分が大きくなる可能性があります。

105 1961 年のこの論文の中で初めて \complete pivoting" と\partial pivoting" の用語が用いられました。

336

n = 5

0
BBB
BB@

のときは
1

0

0

0

01

1

1

0

0

1

01

1

1

0

01

1

01

1

1

1

01

1

01

1

01

1
CCC
CCA

で定義します。
n

測定しました。計算は倍精度でおこなっています。
計算機は FUJITSU M-1800/20、 FORTRAN77 EX
コンパイラです。スケーリングつき部分ピボット選
択の Gauss の消去法サブルーチンとして、 SSL II
拡張機能の DVLAX を使いました ([70])。

です。これを Gauss の消去法、または Crout 法でLU

を 10 から 180 まで 10 きざみで動かし、誤差をrelative error

分解します。部分ピボット選択での行の変換は起きません107 。下三角行列 L は、もとの A5 の下三1036

角成分と完全に一致します。一方 U は、

1031

0
BBB
BB@

1

0

0

0

01

1

0

0

2

1

0

04

1

8

0

016

1
CCC
CCA

1026
1021
1016
1011
106

となり、右下の成分の絶対値が 24 = 2n01 となって

10
20

います。

60

100
matrix dimension

140

180

不安定因子が  = 2n01 の場合の解の振舞い行列の形を A50 で描いてみます。白い部分が 0、黒が 1、グレーが 01 をとります。

n = 60になる前から真の解から外れ、見当違いな方向に飛んでいってしまいました。
この原因を調べてみます。解が飛んで行ってしまう 55 から 60 までの ann の値を書き出します。
n

以下同様に、一般の An を部分ピボット選択で LU

36028797018963968

57

72057594037927936

58

144115188075855870

59

288230376151711740

60

576460752303423490

は、情報をため込む領域ははじめから決まっている
ので、計算機は仕方なく下線部の情報をカットしま

このことより、 n が大きくなると、部分ピボット

す。

選択の Gauss 消去法は非常に不安定になり、 U のその結果、 ann 0 1=2 などの計算は、意味をなくわずかな違いが見当違いの LU の値を生む危険性をし、その演算を繰り返すことによって、必要な情報

はらんでいます。

が次々に落ちてしまい、結果的に答がどんどん正解

なぜ危険かの理論的な説明は [33], [35], [38] に任から離れていきます。

せることにして、さっそく数値実験をします。

この反例より、次の悲観的な結論が導かれます。

Wilkinson さんの定義に従い、行列 An を作ります。 b は、解が x = (1; 1; 1 1 1 ; 1)T となるように

部分ピボット選択の Gauss の消去法は、条件数

n

j =1

56

度以上の情報が必要になります。浮動小数点演算でとなることが確認できます。

i=

18014398509481984

17 桁です。 n を 58 まで大きくすると、仮数部の精n n01(01) 2bnn

55

IBM の浮動小数点の精度は、 10 進数に直すと約分解した場合、上三角行列 U の右下成分が

X

a

a

がよい場合でも不安定となる場合がある。

ij

107 ここがミソです。完全ピボット選択の場合、ピボットが変化
して、不安定因子は常に 2 以下に抑えられます。

337

おそらく、電卓や数式処理ソフトの普及が、紙と8.5鉛筆を使って自力で計算する能力を減退させている

なぜ部分ピボット選択法が使のは間違いないでしょう。

われるのか？
Wilkinson さんの後退誤差解析の結論を見れば、ピボットの部分選択法に関して、 1943 年の Hotellingさんの予想の \4n01 " が \2n01 " になっただけのよ
うに見えます。

プログラム 5.2 は Gauss の消去法を使って連立 1 次方程式を解く。このプログラムの算法が単純である
ことは教育の目的には理想的である。一般に優れた
Gauss-Jordan 法の算法はプログラム 9.7 で示す。

では、理論的にすぐれた手法108 である完全ピボッJ. M. McCormick and M. G. Salvadori (1964)『FORTRAN による数値計算プログラミング』

ト選択法や、さらにすぐれた手法である QR 分解
法109 が、なぜ部分ピボット選択に代わって台頭しな

この本 ([21]) では、 Gauss-Jordan 法に部分ピボッ

かったのでしょうか。

ト選択を施すプログラムをこの後で紹介してます。

それには、 \実際には" という、数値計算におけGauss-Jordan 法が、 Gauss の消去法に比べて「一る膨大な経験則が関わってきます。

般に優れている」理由は、どうやら逆行列の計算がいくつかの本と論文から、 Gauss の消去法に関す簡単にできるためらしいです。

るコメントを年代順に拾っていきます。
疑いもなく密な行列を持つ連立 1 次方程式を解くための最もよく知られた方法は Gauss の消去法およびその種々の変型である。 N 個の未知数に関する N 個の密な方程式は N 3 =3 程度の回数の乗算とそれに通常あ
らい時間評価をする際には無視される他の算術演算で解くことができる。
von Neumann と Goldstein による完全な誤差解析か
らディジタル型計算機で連立 1 次方程式の固定小数点方式の解における積み重ねられた丸めの誤差に対する
くわしい誤差限界が与えられる。

A が対角優位でなく、正定値対称でも Hessenberg でもなく、 U の増大について何も情報をもっていない
場合には、完全ピボット選択を行なう必要がある。
B. Wendro(1966)『理論数値解析』

理論的な立場からは、この結論になるでしょう。
肝心の Wilkinson さんの意見はどうでしょうか。

G. E. Forsythe and W. R. Wasow (1960)『偏微分方程式の差分法による近似解法』

von( = 2n01 となるような) どちらかというと特殊な行列も存在する。しかし、実際には、そのように大きくなるのは非常にまれである。任意の要素が 8 の大きさになることも、ごく珍しいことであり、行列がまっ(k)
たく悪条件であるときは、次々と aij は減少する場
合のほうがずっとありそうなことである。さらに、その成長がもっとゆるやかな因子で制限されていることが \事前的に" わかっているようないくつかの重要な類がある。

Neumann さん と Goldstein さ んは、 1947

年に発表した論文の中で、正定値対称行列の場合に、丸め誤差の影響は、それまで考えられていたほど大きくないことを示しました ([38])。

（完全ピボット選択法を説明した後）消去の進行を最後まで円滑に進行させるには、このようにするのが最も望ましいことである。人間の眼は状況を 2 次元で認識するので、上のような絶対値最大なものの捜索は人手による計算ではむしろ簡単である。機械計算では 1次元の過程でこれを走査していかなれけばならない面倒があるが、しかし機械の高速性を考慮すれば大したことではないので、ぜひこの手つづきを経由するのが望ましい。

J. H. Wilkinson (1963)『基本的演算における丸め誤差解析』\いくつかの重要な類" とは、先ほどの定理であげたように、 A が対称正定値や対角優位や Hessenberg 行列な場合のことです。
つまり、  = 2n01 となるようなものは \例外" で宇野 利雄 (1963)

あり、人工的に意地悪して作らないかぎり、不安定『計算機のための数値計算』

因子が増大することはない、と Wilkinson さんは主張します。

108 誤差が少なくて安定しているという意味。
109 A を Householder 法で QR 分解した場合、不安定因子は n1=2 になります ([38])。

さらに、完全ピボット選択法の問題点も指摘して



338

います。

計算に誤差がないならば、有限回の 4 則演算で解が得られるのが直接法である。この中で、ピボットの選択を伴う Gauss 消去法が計算の手間と精度の上で最も優れた方法とされている。一時期、その精度について
悲観的見方があったが、誤差の後退解析によって、この方法の信頼性が保証された。（中略）完全選択法を(k )用いた場合の方が aij の増大率について、より小さい上界が得られているが、経験上、部分選択法で十分である。

こうして理論的に優れていることが示されても、完全ピボット選択法の利用が有効であるかどうかは疑わしい。とくに、行列の各成分の数値の記憶に磁気テープを使うときには、実際の管理上の問題が完全ピボット選択法においては、重要なことになってくるのである。完全ピボット選択法の r 段階において、最大成分を定めるためには、 (n r + 1)2 回の減算をしなければならないが、浮動小数点方式の演算では、減算の所
要時間は、乗算に比較して無視できるというわけには
いかない。

0

『岩波数学辞典』 (1985)

J. H. Wilkinson (1967)『性質の悪い一次方程式の解法』

そうはいっても、根本的な疑問が残ります。 Trefethen さんの言葉を借りましょう ([30])。

ピボット選択のために必要な比較の回数を調べると
部分選択

n(n + 1)=2

回

完全選択

n(n + 1)(2n + 1)=6

回

Wilkinson による Gauss の消去法についての標準的誤差解析は、ただちに安定性の議論を消去過程で成分のサイズがどのくらい増大するのかという問題に帰着させる。つまり、 P A = LU の作成段階で、 Uの成分が A のそれと比べて極端に大きく成長することがありうるか、ということである。原理的には答えは yes である。簡単な例により 2n01 倍もの大きさにすることができる。しかし現実には決して起きない。
この経験的な結論はあまりによく確立してしまった
ため、次の注目すべき状況を忘れがちである。すなわ
ち、なぜ (ピボットの部分選択付き)Gauss の消去法が安定であるのか、誰も知らないのである！

です。 LU 分解にかかる演算回数 (2=3)n3 回 ([41])と比較して、部分選択の場合は無視できる回数ですが、完全選択の場合、ほぼ同じ回数がピボットを選ぶ作業で必要になります。また、実際にこの点が完全選択法の最大のネックになることが、前章の数値
例を見ることで納得できます。
つまり、ピボットの完全選択法の方が、数少ない(と Wilkinson さんは言っています) 特別な場合によりよい精度の結果を得ることができることは理論

L. N. Trefethen (1985)

的にわかっているのですが、 \経験的に" 部分選択

\Three Mysteries of Gaussian Elimination"法でも安定な結果が得られているので、計算時間の節約のため、部分選択で十分だと Wilkinson さんはおっしゃっているわけです。

8.6

この、 \経験的に安定" と言い切れるだけの数値
実験の積み重ねが、ピボットの部分選択付き Gauss

random matrix での実験の消去法を「最も有名かつ信頼できる (密行列に関L. N. Trefethen さんは、この問題を 5 年間あたする) 連立 1 次方程式の数値解法」の地位に押し上ためた結果、自分の問いかけに対して自分で答を出げることになります。

そうとしました。 R. S. Schreiber さんとの共著で1990 年に発表した論文 [30]110 において、 Trefethen

1970 年代に入ると、数値計算の専門家は次のようさんは、次数が 1024 以下の行列に対し、適当な乱な結論を下します ([13])。

数を発生させて作った random matrixで Gauss の消去法を片っぱしから解き、不安定因子の平均値を行列の計算にとり組んできた人々が過去 15 ないし 20
年間に知った唯一のしかもきわめて重要なことは、次の事実であろう。ピボットの部分的選択を行う Gauss
消去法では残差が小さいことが保証される。

算定しました。
それによれば、平均値はG. E. Forsythe, M. A. Malcolm and C. B. Moler(1977) 『計算機のための数値計算法』

2=3

部分ピボット選択のとき

  n

完全ピボット選択のとき

  n

1=2

となりました111 。
これをグラフにします。

1980 年代に入っても、ピボットの部分選択付き
Gauss の消去法に与えられた地位はますます堅固な110 副題は \Dedicated to the memory of Jim Wilkinson."

ものになります。 [73] を引用します。

です。
111 random matrix がどれくらい有効性を持つのかは、意見が分かれます。 random matrix を用いたテストは大変重要だが、十分ではない、という意見もあります ([14])。

339

最近、 Wilkinson さんが作った人工的な行列では

105

なく、 \実際のモデル" から導かれたそれほど悪条
bound for GECP

bound for GEPP

件でない行列に対して、不安定因子が指数関数的にgrowth factor

104

増加する例がいくつか報告されてましたので、簡単な紹介と行列の作り方をお教えします。

1000

100

average for GEPP

10average for GECP

1
10

100

8.7.1 Wright さんの例
1 つめは、 S. J. Wright さんが 1993 年に [39] で発表した例です。 IRn の 2 点境界値問題

n
1000

y

0

= M (t)y + q (t);y (a) = y (b) =を考えます。求めたいのは y (t) 2 IRn です。 M (t)不安定因子の増大の様子

は t についての n 2 n の関数行列、 q (t) はベクトル上の二つの曲線は、 Wilkinson さんが証明した不です。この問題を標準的な shooting method112 を適安定因子  の上限 (19), (20) です。下の二つは Tre-

用し離散化します。 Wright さんは特別な場合としfethen さんが random matrix を用いて計算した不てn=2の安定因子の平均値です。
L. N. Trefethen さんと R. S. Schreiber の実験がM (t) =

実際に数値計算で使用される一般の行列をうまく反映していると信じるならば、このグラフによって不

0
BBB
BBB
B@

01=6

1

1

01=6

!

a = 0;

;

1
CCC
CCC
CA

を考えました。係数行列 A は次のようになります。

安定因子の増加は最悪の理論値よりも相当低いことになります。
不安定因子が指数関数的に増加する行列は、平均値から遠く離れた極めて特異なものであり、もし行A =列の条件数がめちゃくちゃに大きくなければ、部分
ピボット選択付き Gauss の消去法は安定であると結

I

Mh

0e

I

0
.
.
.

I

Mh
0e

I

..

..

.

I

この頼もしい結果も取り込んで、現在数値計算をする人の共通の認識は次のようになっています。

0I

は 2 次の単位行列、 h は区間 [a; b] の分割幅、Mh = I +M h+O(h2 ) 0e

部分的ピボッティングのコストは、計算全体のコストに比して無視できるくらい軽微なものであるが、完全ピボッティングのコストは大きく、計算全体そのものと同程度である。一方、実用上は部分的ピボッティン
グで十分であることがわかっているので、通常は部分的ピボッティングを採用する。部分的ピボッティングが有効に働かないで誤差の増大を起こすような例題も作れるが、それは病理学的な例で、実際上、めったに出会わないものである。

.

Mh
0e

論づけられます。

h=6 0 10h0h

h=6 0 1!

です。 A は h が小さいなら、それほど条件数は悪くなりません。行列の作り方は次のように極めて簡単です。

津田 孝夫 (1988)

10『数値処理プログラミング』

30

8.7

b = 60

病理学的な例「めったに出会わない」ということは、たまには40

出会うということです。その意味で Trefethen さん
が「現実には決して起きない」と言うのは行き過ぎ

read(5,*) n! 分割数
h = 60.0D0/n! 分割幅
m=2*n+2! 行列の次数do 10 i=1,m
do 10 j=1,mA(i,j) = 0.0D0 ! 初期化
continue
do 30 i=1,mA(i,i) = 1.0D0 ! 対角成分
continue
do 40 i=1,nA(2*i+1,2*i-1) = h/6-1.0D0
A(2*i+2,2*i-1) = -hA(2*i+1,2*i)= -hA(2*i+2,2*i)= h/6-1.0D0
continueA(1,m-1) = 1.0D0 ! 残りA(2,m)= 1.0D0です。
112 手法は [39] の参考文献を御覧下さい。

340

ところが、この行列 A を係数行列にする連立 1

する 2 つ目の例として、 L. V. Foster さんが 1994
年に [14] に発表した論文115 を紹介します。

次方程式を部分ピボット選択付き Gauss の消去法

問題は 2 階の線形 Volterra116 型積分方程式から

で LU 分解すると、不安定因子が指数関数的に増大してしまいます113 。 Wilkinson さんの行列と同じ導かれる

く、解が x = (1; 1; 1 1 1 ; 1)T となるように b を決め

x(s) 0て数値実験します。計算機は FUJITSU M-1800/20
です。比較のサブルーチンは

GECP(完全ピボット選択)

GEPP(部分ピボット選択)



QR 分解

Z

0

s
k (s; t)x(t) dt +(s)x(L) = G(s)を一般的な求積法で近似するものです。求めたいのは区間 [0; L] 上の x(s)、その他の関数は既知です。
Foster さんは、既知の関数を連続問題の解が具体的にわかるように選び、その離散近似を行いました。
積分の近似は Newton-Cotes117 公式と Simpson118 ので行いました114 。

公式を組み合わせて用います。
10

0

GEPP
GECP
QR_decomposition

relative error

10-2
10-4

この近似から導かれる行列 A は悪条件ではあり

partial pivotingません。しかし、部分ピボット選択による LU 分解をすると、ある n 以上ではピボット選択をしなくな

10-6

り、そのため不安定因子が急激に増大します。

10-8

Foster さんは素晴らしいことに、

10-10
10-12

ワーク上で公開していました119 。そこで、プログラ

QR decomposition

10-14

ムをとり寄せて Fortran にプログラムを組み直し、

complete pivoting

10-16
0

QR

MATLAB という行列演算ソフトを使ったプログラムをネット

50

100
matrix dimension

150

200

数値実験しました。

分解と完全ピボット選択は、期待通りの精度を
DLAX
DGEFA/SL
GECP
GENP

104

示しています。一方、部分ピボット選択では、行列の次数が 120 くらいで見当違いな解に発散します。

102

原因は、 Wilkinson さんの例題と同じで、不安定因
relative error

子の指数関数的な増加のためです。
参考のために、行列 A の形を graphics にします。
□が 0、■が 1、残りはグレーで書きます。

100
10-2
10-4
10-6
0

50

100
matrix dimension

150

200

これは [14] の最初の例題です。誤差は、連続問題の真の値と近似解の相対誤差です。 DLAX は部分
ピボット選択付きの Crout 法の LU 分解、 DGEFA/SL
は LINPACK のサブルーチンです。また、 GENP
はスケーリングもピボット選択もしない Gauss の消去法の解法、 GECP はおなじみスケーリングと完

Wright さんの作った行列の形

全ピボット選択を行なったサブルーチンです。

8.7.2 Foster さんの例

115 \Gaussian Elimination with Partial Pivoting Can Failin Practice" という、そのものズバリの表題です。

116 Vito Volterra(1860{1940) さんはイタリアの数学者です。

実際の数学モデルから導かれる連立 1 次方程式に117 Roger Cotes さん (1682{1716) の名前 が半分入っていま

おいて部分ピボット選択付き Gauss の消去法が失敗

2す。
118 Thomas Simpson さん (1710{1761) から来ています。

1021 くらいに増大します。
114 GEPP は SSL II の DVLAX、 QR 分解は同じく SSL II

113 n = 200 程度で 2:59

119 アクセス先は論文に書いてあります。

の Householder 変換の DLAXL を用いました。

341

のはかなり難しいでしょう121 。

GECP が順調に収束していくのに対し、他のサブ

逆に、行列の形から部分ピボット選択の Gauss のルーチンは、ある n から急に力つきて飛んで行ってしまう様子がよくわかります。

消去法が失敗することを事前に判断する基準ができるならば、これは素晴らしいことです122 。

1017

difference of elements

1015
DLAX
DGEFA/SL
GECP
GENP

1013
10

11

109
107
105
103
101
10-1
0

50

100
matrix dimension

150

200

Foster さんの作った行列の形

こちらは、不安定因子の値をプロットしたものでFoster さんは、論文の最後の章の「なぜ部分ピボッ

す。完全ピボット選択以外のプログラムが真の解か

ト選択が失敗する実際の例がこれまで報告されなら突然離れていく原因が、不安定因子の急激な増大

かったのか？」という項で、次のように述べていまによる丸め誤差の累積によることがわかります120 。

す ([14])。

[14] にはもう一つ例題があり、上と同様な挙動を我々は (失敗する例がこれまで報告されなかった) 本
質的な理由は、このような例が極めて \まれ" であるからだと期待する。問題の設定と解法へのアプローチ
は慎重に選択されたものである。しかしながら、数値線形計算のソフトウェアパッケージのほとんどが不安定因子について何の情報も提供しておらず、そして不安定因子の増大が、時々ではあるが、実際に発生している可能性があることを指摘すべきであろう。

示します。こっちは A がすぐに作れます。 Fortran
で書けば、

20

30

40

read(5,*) n! 行列の次数k = 1.0D0! 定数 1C = 6.0D0! 定数 2
L =40.0D0! 区間幅
h = L/dble(n-1)! 分割幅
do 20 i=1,ndo 20 j=1,nA(i,j) = 0.0D0 ! 初期化
continue
do 30 i=3,ndo 30 j=2,i-1
A(i,j) = -k*h
continue
do 40 i=2,nA(i,1) = -k*h/2.0D0A(i,i) = -k*h/2.0D0 + 1.0D0
A(i-1,n) = -1.0D0/C
continueA(1,1) = 1.0D0
A(n,n) = 1.0D0 - 1.0D0/C - k*h/2

L. V. Foster (1994)\Gaussian Elimination with Partial
Pivoting Can Fail in Practice"\速さ" の名のもとに疎外されている完全ピボット選

択付き Gauss の消去法が \精度" の名のもとに頼りにされるのは、このような場面です。
おそらく、ほとんどの専門家が指摘する通り、部
分ピボット選択が実際の数値計算において致命的な不安定因子の増大を起こすことは \まれ" なのでしょう123 。しかし、プログラムのバグ取りも兼ねて
数値計算の信頼性をあげるために複数の解法を試しで部分ピボット選択付きの Gauss の消去法が倍精てみることは、大変効果があります124 。そんなと度演算で失敗する例ができます。

き、完全ピボット選択付き Gauss の消去法は (2 倍

この行列の形状を次のページに描きます。□が以上の実行時間がかかりますが) うってつけの解法

0、グレーが 0kh=2 の部分です。

だといえます。そして、もしかしたら皆さんを \致Wilkinson さん、 Wright さん、 Foster さんの行

命的な失敗" から救ってくれるかも知れません。

列の形を見比べると、なんとなく似ていることに気121 条件を限定した random matrix を作れば、違う結果が出るづきます。このような形の行列を乱数で発生させるかもしれません。
122 ただし、 \この形以外は安定である" というのを証明するのは120 これだけ不安定因子が増大すると、反復改良をかけても、行列が \特異" だと判断されてしまいます。

大変だと思います。
123 なんせ見つけただけで論文になるのですから。

124 ついでに精度や定数もいろいろ変えてみましょう。

342

参考文献と索引

あとがき

参考文献

この原稿は、もともと連立 1 次方程式の数値解法[1] Benoit, Commandant : Note sur une methode deresolution des equations normales provenant del'application de la methode des moindres carres a
un systeme d'equations lineaires en nombre inferieura celui des inconnues |Application de la methodea la resolution d'un systeme de ni d'equations
Lineaires (Precede du Commandant Cholesky),
Bulletin Geodesique, No.2, pp.67{77 (1924).

の一覧表として (自分のために) 書き始めした。現
在の研究のいきがかりから、連立 1 次方程式の解の厳密な評価 (浮動小数点演算を使う以上、ある幅を持った区間の中に解を包み込むことになります) が必要になったからです。ついでに最大固有値の厳密な評価もいるのですが、こっちは手つかずです。
解の評価のためには精度良い近似値が必要です。

| \Commandant"Cholesky の業績を紹介した長近似解を得る手段は \どんな方法を用いても" 構いいタイトルの論文です。おそらくここから \Choleskyません。そこで、世間にはどんな解法があるのだろ分解" の名前が普及したものと思われます。東京大学

うと思い、数値解析の教科書などを見ながらメモしに所蔵されています。

[2] Bodewig, E. : Matrix Calculus, North-Holland
Publishing Company-Amsterdam (1956).

ていった結果がこの記事となりました。
もともとが \精度がいちばん、速度は二の次" という妙な視点で見ていったため、 6 章以下、特に 8| 反復法の歴史が詳しく書いてあります。

[3] Booth, A. D. : Numerical Methods, Butterworths,
London (1957).

章は部分ピボット選択付き Gauss の消去法が失敗する可能性ばかりを探ってしまいました。

| Hamilton-Cayley の定理に基礎をおく A01 の

しかし、 Gauss の消去法の歴史や丸め誤差についての議論など、ご存知なかった方にはなかなか面白求め方が載っています。

[4] Bunch, James R., Kaufman, Linda and Parlett,Beresford N. : Decomposition of a Symmetric
Matrix : Numerische Mathematik, Vol.27, pp.95{
109 (1976).

かったのではないでしょうか。また、 Wilkinson さんや Forsythe さんなどの名前が出てきて懐かしく思われた方もいらっしゃると思います。
さて、肝心の厳密な解の評価ですが、こちらは[5] Bunch, James R. and Kaufman, Linda : Some
Stable Methods for Calculation Inertia and Solving Symmetric Linear Systems, Mathematics ofComputation, Vol.31, pp.163{179 (1977).

\区間解析" と \不動点定理" という数学的な側面と、アルゴリズムを実現する精度保証付き数値演算ソフ
トの使い方という \利用の手引" の側面があって、今回は解説する時間 (および能力) がありませんでし| [4], [5] では、対称行列に対する安定した数値解た。近いうちに、この原稿の続編として紹介できれ

法 (M DM T 分解) を提案しています。具体的なアルば思います。

ゴリズムおよび ALGOL プログラムがついています。

[6] Du , I.S., Erisman, A. M. and Reid, J. K. : Direct
Methods for Sparse Matrices, Clarendon Press,
Oxford (1986).

参考文献について以下の参考文献の中には、『古典的名著』の誉れ| 疎行列の形状の graphics が付録に載っています。

を受けながら絶版になっている本も少なくありませ

行列の形状の図はこれを真似ました。

[7] Dyksen, Wayne. R. and Rice, John R. : The Importance of Scaling for the Hermite Bicubic Collocation Equations, SIAM Journal on Scienti c
and Statistical Computing, Vol.7, No.5, pp.707{
719 (1986).

ん。一応、孫引きは避け、原典を極力見るようにしましたので、参考文献のほとんどは九州大学中央図書館または大型計算機センター図書室にあるか、または書店で注文すれば手に入ります (1995 年 11 月
現在)。文献に簡単な紹介をつけました。何かの参考| 数値解析におけるスケーリングの必要性を豊富なになれば幸いです。

データを駆使して訴えています。

343

[8] Forsythe, George E. : Gauss to Gerling on Relaxation, Mathematical Tables and other Aids toComputation (現在 Mathematics of Computation),
Vol.5, pp.255{258 (1951).

[15] Freund, Roland W., Nachtigal, Noel M. : An Implementation of the QMR Method Based on Coupled Two-Term Recurrences, SIAM Journal on
Scienti c Computing, Vol.15, No.2, pp.313{337(1994).

| Gauss さんが弟子の Gerling さんにあてて書いた手紙の英訳です。

| QMR 法を提案した本人が解説しています。

[9] Forsythe, George E. : Tentative Classi cation of
Methods and Bibliography on Solving Systems ofLinear Equations, in Simultaneous Linear Equations and the Determination of Eigenvalues Editedby L. J. Paige and Oliga Taussky, National Bureau of Standards Applied Mathematics Series 29(1953).

[16] Girault, Vivette and Raviart, Pierre-Arnaud : Fi-

nite Element Methods for Navier-Stokes equations,
Theory and Algorithms, Springer Series in Com-

putational Mathematics 5, Springer-Verlag (1986).
| Navier-Stokes 方程式に対する有限要素法の解法を解説した、その方面では有名な本です。

| 連立 1 次方程式の解法についての総合分類および当時の最新の参考文献の一覧です。あきれるくらい完璧に分類されています。

[10] Forsythe, George E. : Solving Linear Algebraic
Equations Can be Interesting, Bulletin of the American Mathematical Society, Vol.59 (1953).
| こちらは [9] から興味ある解法をピックアップして、適切な解説を加えた論文です。

[11] Forsythe, George E. and Wasow, Wolfgang R., :Finite-Di erence Methods for Partial Di erential
Equations, John Wiley & Sons, New York (1960).

[17] Goldstine, Herman H. : A History of Numerical
Analysis from the 16th through the 19th Century,
Studies in the History of Mathematics and Physical Sciences 2, Springer (1977).
| 16 世紀から 19 世紀にかけての数値解析の歴史が書いてあります。

[18] Gregory, R. T. and Karney, D. L. : A collectionof matrices for testing computational algorithms,
John Wiley & Sons, New York (1969).
| 連立 1 次方程式、逆行列、固有値問題などのテス

ト行列を集めた本です。

[19] Householder, Alston S. : Unitary Triangularization of a nonsymmetric Matrix, Journal of the
Association for Computing Machinery, Vol.5, pp.339{
342 (1958).

［邦訳］ 藤野 精一 訳 : 偏微分方程式の差分法による近似解法 (上),(下), 吉岡書店 (1968).
| 偏微分方程式の数値解法についての古典的名著として知られています。訳者の藤野先生 (我々は \大御大" とお呼びしています) には、修士の頃に御指導い| 後に Householder 変換と呼ばれる手法が最初にただきました。

載った論文です。

[20] Macon, Nathaniel : Numerical Analysis, John
Wiley & Sons, New York (1963).
| LLT 分解が \geodetic surveyor" の Cholesky[12] Forsythe, George E., Moler, Cleve B., : Computer
Solution of Linear Algebraic Systems, PrenticeHall, (1967).
［邦訳］ 渋谷 政昭, 田辺 国士 : 計算機のための線形計算の基礎 | 連立 1 次方程式のプログラミング |,培風館 (1969).

さんによって 1916 年に示された、と書いてある本です。

[21] McCormick, John M. and Salvadori, Mario G. :Numerical Methods in FORTRAN, Prentice-Hall(1964).
［邦訳］清水 留三郎 訳 : FORTRAN による数値計
算プログラミング, サイエンスライブラリ 情報電算機
1, サイエンス社 (1970).

| 連立 1 次方程式の数値解法に関する古典的名著といわれています。序文で「この本を J. H. Wilkinson に捧げたい」と書いてあります。

[13] Forsythe, George E., Malcolm, Michael A. andMoler Cleve B. : Computer Methods for Mathematical Computations, Prentice-Hall (1977).
［邦訳］ 森 正武 訳 : 計算機のための数値計算法, 日| Fortran90 の廃止予定事項に入ってしまった部分

がかなりあります。

本コンピュータ協会編 コンピュータ・サイエンス研
究所シリーズ, 科学技術出版社 (1978).

[22] Meijerink, J. A. and van der Vorst, H. A. : An
Iterative Solution Method for Linear Systems ofwhich the Coecient Matrix is a Symmetric MMatrix, Mathematics of Computation, Vol.31, pp.148{
162 (1977).

| 完成したプログラムに重点をおいた解説がされています。 Forsythe さんの遺作です。

[14] Foster, Leslie V. : Gaussian Elimination with
Partial Pivoting Can Fail in Practice, SIAM Journal on Matrix Analysis and Applications, Vol.15,
No.4, pp.1354{1362 (1994).

| 不完全 LU 分解、不完全 Cholesky と共役勾配法の相性のよさを指摘して共役勾配法を復活させた画期的な論文です。ちなみに Bunch さんの論文がこの| 部分ピボット選択付き Gauss の消去法が不安定

すぐあとに続きます。

となる例を実際の問題から導いて紹介しています。

344

[32] Varga, Richard S. : Matrix Iterative Analysis,
Prentice-Hall, (1962).
［邦訳］ 渋谷 政昭, 棚町 芳弘, 金子 正久, 野田 隆
訳 : 計算機による大型行列の反復解法, サイエンスラ
イブラリ 情報電算機 10, サイエンス社 (1972).

[23] Puschmann, Heinrich and Cortes, Joaquin : The
Coordinex Problem and Its Relation to the Conjecture of Wilkinson, Numerische Mathematik,
Vol.42, pp.291{297 (1983).
| 完全ピボット選択法についての Wilkinson さん| 連立 1 次方程式の反復解法のアルゴリズムと数学の推測を扱った論文です。

[24] Rainsford, Hume F. : Survey Adjustments andLeast Squares, Constable & Company (1957).

的背景を詳しく論じた古典的名著です。

[33] Wendro , B., : Theoretical Numerical Analysis,
Academic Press, New York (1966).
［邦訳］ 戸川 隼人 訳 : 理論数値解析 | 計算機による数値計算の基礎 |, サイエンスライブラリ 情報電算機 19, サイエンス社 (1973).

| pp.74 に、「\Commandant Cholesky" の提案した方法が最初に公なった 1924 年の論文は私が英語に訳した」と書いてある本です。そこ以外は読んでないので、内容は知りません。

Veri cation methods for
dense and sparse system of equations, Berichte| 数値解析における理論的な側面を平易に解説した[25] Rump, Siegfried M. :本です。

[34] Westlake, Joan R. : A Handbook of Numericaldes Forschungsschwerpunktes Informations- und
Kommunikationstechnik, TUHH (1993).

Matrix Inversion and Solution of Linear Equations, John Wiley & Sons, New York (1968).

| 精度保証付き数値計算を用いた連立方程式の解法［邦訳］ 戸川 隼人 訳 : コンピュータのための線形計算ハンドブック, 培風館 (1972).

を扱っています。

[26] Saad, Youcef, Schultz Martin H. : GMRES: A
Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems, SIAM Journal on Scienti c and Statistical Computing, Vol.7,
No.3, pp.856{869 (1986).

| 連立 1 次方程式の数値解法をきめ細かく分類した

ハンドブックです。

[35] Wilkinson, J. H., : Error Analysis of Direct Methods of Matrix Inversion, Journal of the Association for Computing Machinery, Vol.8, pp.281{
330 (1961).

| GMRES 法の原論文です。

[27] Stummel, Friedrich : Forward Error Analysis ofGaussian Elimination Part I : Error and Residual Estimates, Numerische Mathematik, Vol.46,
pp.365{395 (1985).

| Gauss の消去法によって発生する浮動小数点演算の誤差を、後退誤差解析によって評価した古典です。

[36] Wilkinson, J. H., : Rounding Errors in Algebraic
Processes, Prentice Hall (1963).
［邦訳］一松 信, 四条 忠雄 訳 : 基本的演算における丸め誤差解析, 培風館 (1974).

[28] Stummel, Friedrich : Forward Error Analysis ofGaussian Elimination Part II : Stability Theorems, Numerische Mathematik, Vol.46, pp.397{
415 (1985).

| 誤差解析の基礎理論を解説した本です。

| [27], [28] は前進誤差解析を扱っています。 Part
I, Part II に分かれています。

[29] Trefethen, Lloyd N. : Three Mysteries of Gaussian Elimination, ACM SIGNUM Newsletter, Vol.20,
No.4, pp.2{5 (1985).

[37] Wilkinson, J. H., : The Algebraic Eigenvalue Problem, Oxford University Press, London (1965).
| それまでの Wilkinson さんの成果を結集して書かれた 600 ページを越える大作です。一時期、数値解析の論文では必ずこの本が参考文献に載せられて| 計算量、直交変換法との関係、安定性について論いました。

じています。

[30] Trefethen, Lloyd N., Schreiber, Robert S. : AverageCase Stability of Gaussian Elimination, SIAM
Journal on Matrix Analysis and Applications, Vol.11,
No.3, pp.335{360 (1990).
| random matrix を用いて不安定因子を計算することで Gauss の消去法の平均的な安定性を主張する[38] Wilkinson, J. H., : 性質の悪い一次方程式の解法, in
Mathematical Methods for Digital Computers II,
Edited by Anthony Ralson, Herbert S. Wilf, John
Wiley & Sons, New York (1967).
［邦訳］片岡 信二 監訳 : 電子計算機のための数学的方法 II, 鹿島出版会 (1975).
| Wilkinson さんの方法のダイジェスト版です。

論文です。

[31] van der Vorst, H. A. : Bi-CGStab: A Fast andSmoothly Converging Variant of Bi-CG for the
Solution of Nonsymmetric Linear Systems, SIAM
Journal on Scienti c and Statistical Computing,
Vol.13, No.2, pp.631{644 (1992).

[39] Wright, Stephen J. : A Collection of Problems
for which Gaussian Elimination with Partial Pivoting is Unstable, SIAM Journal on Scienti c
Computing, Vol.14, No.1, pp.231{238 (1993).
| 部分ピボット選択付き Gauss の消去法が失敗する例をあげた論文です。

| BiCGStab 法の原論文です。

345

[40] Zeidler, E. : Nonlinear Functional Analysis and
Its Applications, Springer-Verlag, Berlin (1986).

[62] 渡部 善隆 : 数値計算と誤差の話, 九州大学大型計算機
センター広報, Vol.27, No.5, 556{580 (1994).
[63] 渡部 善隆 : 数値計算と速度の話, 九州大学大型計算機
センター広報, Vol.28, No.3, 207{231 (1995).

| 関数解析の世界では有名な本です。たまに証明が間違っていたりします。

::::::::::::::::::::::::::::::::::::::::::::::[64] 線型代数ア・ラ・カルト, 数学セミナー・リーディン
グス, 数学セミナー増刊, 日本評論社 (1990).

以下は国内の数値計算に関する文献です。

[41] 伊理 正夫, 藤野 和建 : 数値計算の常識, 共立出版
(1985).

[65] 精度保証付き数値計算とその応用, 情報処理 Vol.31,
No.9 (1990).

[42] 宇野 利雄 : 計算機のための数値計算, 応用数学力学講座 14, 朝倉書店 (1963).

[66] 情報検索システム AIR |AIR を用いた文献情報の検索 |, 九州大学大型計算機センター (1995).

[43] 唐木 幸比古 : スーパーコンピュータと行列計算, 情報
処理, Vol.28, No.11, pp.1441{1451 (1987).

[67] OS IV/MSP FORTRAN77 EX 使用手引書 V12
用, 79SP-5031，富士通株式会社 (1991).

[44] 柴垣 和三雄, 関数解析と数値計算の基礎, 森北出版
(1985).

[68] OS IV/MSP FORTRAN77 EX/VP 使用手引書 V12
用, 79SP-5041，富士通株式会社 (1991).

[45] 島崎 眞昭 : スーパーコンピュータとプログラミング,
計算機科学／ソフトウェア技術講座 9, 共立出版 (1989).

[69] SSL II 使用手引書 (科学用サブルーチンライブラリ),
99SP{4020, 富士通株式会社 (1987).

[46] 数値流体力学編集委員会 編 : 非圧縮性流体解析, 数値流体力学シリーズ 1, 東京大学出版会 (1995).

[70] SSL II 拡張機能使用手引書 (科学用サブルーチンライブラリ), 99SP{4070, 富士通株式会社 (1991)．

[47] 杉原 正顕, 室田 一雄 : 数値計算法の数理, 岩波書店
(1994).

[71] NUMPAC 利用手引書, 富士通株式会社 (1994)．

[48] 谷口 博 : 機械工学における電子計算機の応用, サイエ
ンスライブラリ情報電算機 34, サイエンス社 (1976).

[73] 岩波数学辞典 (第 3 版), 岩波書店 (1985).

[72] 現代数理学事典, 大阪書籍 (1991).
[74] アルゴリズム辞典, 共立出版 (1994).

| Cramer の公式の Fortran プログラムが載っています。

[49] 津田 孝夫 : 数値処理プログラミング, 岩波講座ソフト
ウェア科学 9, 岩波書店 (1988).
[50] 戸川 隼人 : 数値計算, 情報処理入門コース 7, 岩波書
店 (1991).
[51] 戸川 隼人 : マトリクスの数値計算, オーム社 (1971).
[52] 戸川 隼人 : 共役勾配法, シリーズ新しい応用の数学
17, 新曜社 (1977).
[53] 戸川 隼人 : 微分方程式の数値解法 | 有限要素法と差分法 |, オーム社 (1973).
[54] 名取 亮 : 線形計算, すうがくぶっくす 12, 朝倉書店
(1993).
[55] 名取 亮 : 数値解析とその応用, コンピュータ数学シ
リーズ 15, コロナ社 (1990).
[56] 名取 亮, 野寺 隆 : 大規模行列計算における反復解法,
情報処理, Vol.28, No.11, pp.1452{1459 (1987).
[57] 名取 亮, 野寺 隆 編: スーパーコンピュータと大型数
値計算, bit 臨時増刊, 共立出版 (1987).
[58] 一松 信, 戸川 隼人 編 : 数値計算における誤差, bit 臨
時増刊, 共立出版 (1975).
[59] 村田 健郎, 小国 力, 唐木 幸比古 : スーパーコンピュータ | 科学技術計算への適用 |, 丸善 (1985).
[60] 森 正武 : 数値解析法, 朝倉現代物理学講座 7, 朝倉書店 (1984).
[61] 森 正武, 杉原 正顕, 室田 一雄 : 線形計算, 岩波講座
応用数学 8, 岩波書店 (1994).

Cholesky 法を紹介した 1924 年の論文 ([1])

346