47

第 5 章 計算量について
それではあんまりだということで，二十年a ほど前から“計算複雑度 (computational complexity)”あるいは “計算量”の理論というのが登場した。(略)「この問題をこのプログラムで得にはこのくらいの手間がかかる」というような形の議論をきちんとやろうというものである。
この計算複雑度の理論の中では数値計算の技術に直接関係するような話題が多数取り上げられ，それらに対して，少なくとも理論的には，大変興味ある成果が多く得られている。そうなれば，数値計算に関心のある者としては，そのような成果のうち重要なものは常識として知っていなければならないということになる。
伊里正夫・藤野和建「数値計算の常識」(共立出版)
a この本の初版は 1985 年であるから，それよりも 20 年前という意味であろう。

現在のコンピュータは bit 単位の論理演算を多数組み合わせることで，複雑な浮動小数点演算を実現している。そして，ユーザが直接扱うことの出来る機械語命令は• byte(or word) 単位の論理演算及び bit 列操作• 符号付き (signed) もしくは 符号なし (unsigned) 整数演算
• 浮動小数点数演算（四則演算，初等関数）に大別される [45]。数値計算のアルゴリズムの優劣を論じる時，一つの指標として演算回数の大小がよく使用される。特に 1CPU の PC や WS で，アルゴリズムを逐次的に実行する場合は計算時間の増大を演算回数の order で予測することがある程度は可能である。本章では，数値計算で重要な役割を果たす浮動小数点演算に焦点を絞って，以下の章で紹介する各種アルゴリズムの計算量を考察する基礎を提供する。

5.1

浮動小数点数の四則演算と Landau の O 記号

小学校で習う小数の加減乗除は整数のそれと本質的には同じものである。人間の感覚で言う「面倒くさい」計算は，そのままコンピュータについても当てはまる。面倒な計算は時間を要する。従って四則演算は前章のベンチマークテストからも分かるように，T (加算 (FADD)) = T (減算 (FSUB)) ≤ T (乗算 (FMUL)) < T (除算 (FDIV))(5.1)

第5章48

計算量についてという順に計算時間を要すると考えて良い。後で述べる初等関数はこれらの演算を組み合わせて実行されるため，更に時間を要するのが普通である1 。
従って，数値計算のアルゴリズムの計算時間は加減算の実行回数以上に，乗除算や初等関数の実行回数に左右される。「アルゴリズムの演算回数」という言葉がしばしば後者の意味で使われるのはそのためである。
計算回数に限らず，様々な場面で使用される言葉としてオーダ (order) がある。これは以下に示す Landau の O(ラージオー) と同義である。
定義 5.1.1 (Landau の O 記号)
ある一変数実関数 f (x), g(x) ∈ R に対して， f (x) = O(g(x)) とはlim

f (x)x→α g(x)= 定数 (, 0)となることを意味し，このような f (x) は g(x) のオーダ (order) であると呼ぶ。この O(g(x)) を Landauの O(ラージオー) 記号という。α としては 0 もしくは ±∞ がよく使用される。
またlim

f (x)

x→α g(x)=0であるときは特に f (x) = o(g(x)) と書き，これを o(スモールオー) 記号と呼ぶ。
以降，オーダという言葉は O(ラージオー) の意味で使用する。g(x) としてよく使用されるのは xの多項式であり，特に x2 , x3 , ...xn である。O(xn ) はほぼ xn に比例していることを表しており，直感的に理解しやすいため，様々な場面で使用される。

5.2

複素数の四則演算

本書は実数の演算が主体であるが，複素数の演算が必要となる場面に遭遇することもある。ここ
で復習も兼ねて，複素数の演算とその演算量について若干の考察を行う。
√
任意の複素数 c = Re(c) + Im(c) −1 ∈ C は 2 つの実数の組 (Re(c), Im(c)) として表現できる。従って，複素数の四則演算は全て実数のそれを組み合わせることによって実現できる。

|a|=a±b =ab=√(Re(a))2 + (Im(a))2√(Re(a) ± Re(b)) + −1(Im(a) ± Im(b))√(Re(a)Re(b) − Im(a)Im(b)) + −1(Im(a)Re(b) + Re(a)Im(b))

ab|b|2√
ここで b = Re(b) − Im(b) −1 である。
a/b =

(5.2)
(5.3)
(5.4)
(5.5)

1 現在の CPU は内部に積み込んだ高速転送可能なキャッシュ(cache) メモリを持っており，一度メインメモリから読み込んだ値をそこに記憶しておき，2 度目以降のアクセスはそれを取り出すだけで済む。従って，このキャッシュメモリをうまく利用できるようにした線型計算プログラムは，素朴に組んだものより高速になる。よって単純に計算量だけでは計算時間を推定できないこともある。

5.3. 基本線型計算49

表 5.1: 複素数演算の計算回数
加減算 乗算 除算 平方根|a| ((5.2) 式)

120

1|a| ((5.6) 式)

1211a±b2000

ab2400a/b ((5.5) 式)

3620a/b ((5.7) 式)

3330

但し，オーバーフローを防止するため，|a| と a/b は次のように計算するのが良いとされている

[13]。




|a| = 

















a/b = 



Re(a)Im(a)√( Im(a) )2|Re(a)| 1 + Re(a)√( Re(a) )2|Im(a)| 1 + Im(a)(if Im(a) = 0)
(if Re(a) = 0)
(if |Re(a)| ≥ |Im(a)| > 0)
(if |Im(a)| > |Re(a)| > 0)

s

計算不能
()
−Re(a) · Im(b)Re(b) + Im(a) √+−1s( Re(b) )Re(a) · Im(b)+ Im(a))(−Re(a) + Im(a) · Re(b)Im(b) √( Im(b) )Re(a) + Im(a) · Re(b)s+s(5.6)−1(if b = 0 (即ち Re(b) = Im(b) = 0))(if |Re(b)| ≥ |Im(b)| ≥ 0)
ここで s = Re(b) + Im(b) ·( Im(b) )
Re(b)

(5.7)

(if |Im(b)| ≥ |Re(b)| ≥ 0)()
ここで s = Re(b) · Re(b)Im(b) + Im(b)以上の複素数演算の計算回数を表 5.1 にまとめておく。対応する実数の演算と比べて，2〜3 倍の演算量を必要とすることが分かる。従って，複素数の演算は実数のそれに比べてかなり「高くつく」ことを認識しておく必要がある。当然のことながら，必要となるメモリ量も実数の 2 倍になる。
問題 5.2.1

1. 式 (5.6)，(5.7) がそれぞれ |a| と a/b を計算していることを確認せよ。
2. 前章の多倍長浮動小数点数の四則演算のベンチマークテストの結果を用いて，複素数演算の性能評価を行え。また実際にベンチマークテストを行った結果と比較せよ。

5.3

基本線型計算

現在の数値計算は大規模化が進んでおり，そこではベクトル及び行列の基本線型計算 (linearcompuation) が多用される。基本線型計算は次元数が上がるにつれて莫大な計算量を必要とすることを認識しておく必要がある。ここではその一端に触れることにする。

第5章

50

計算量について実ベクトル a = [a1 a2 · · · an ]T ∈ Rn の基本線型計算はそれぞれ次のように行う。

 αa1 
 αa2 

スカラー積 αa =  .  (ここでα ∈ R)
 .  . 


αan

 a1 ± b1 
 a2 ± b2 

加減算 a ± b = ..


.


an ± bnn∑
内積 (a, b) =ai bi

(5.8)

(5.9)

(5.10)

i=1

同様に，実正方行列 A = [ai j ] ∈ Mn (R) (i, j = 1, 2, ..., n) の基本線型計算は次のようにして行われる。

 ∑n
 j=1 a1 j b j  ∑n j=1 a2 j b j 

ベクトルとの積 Ab = ..
.
 ∑nabnjjj=1
 αa11 αa12 · · · αa1n  αa21 αa22 · · · αa2n 
スカラー積 αA =  .
..
..  (ここでα ∈ R)
 .
.
.
.


αan1 αan2 · · · αann a11 ± b11 a12 ± b12 · · · a1n ± b1n 
 a21 ± b21 a22 ± b22 · · · a2n ± b2n 
加減算 A ± B = ..
..
..


.
.
.


an1 ± bn1 an2 ± bn2 · · · ann ± bnn
 ∑n
∑n
∑n
 j=1 a1 j b j1···
j=1 a1 j b j2
j=1 a1 j b jn  ∑n∑n∑n j=1 a2 j b j1···
j=1 a2 j b j2
j=1 a2 j b jn 
乗算 AB = ..
..
..


.
.
.

∑n
∑n
 ∑n···
j=1 an j b j1
j=1 an j b j2
j=1 an j b jn
ベクトル演算，行列演算の計算回数を表 5.2，5.3 にまとめておく。
問題 5.3.1

1. α, β ∈ R，a, b ∈ Rn であるとき，αa ± βb の計算量を求めよ。
2. α, β ∈ R, A, B, C ∈ Mn (R) である時，(αA ± βB)C の計算量を求めよ。

(5.11)(5.12)(5.13)(5.14)

5.4. 収束判定と停止則について51

表 5.2: ベクトル演算の計算回数

表 5.3: 行列演算の計算回数

加減算

乗算

加減算

乗算

αa0nAbn(n − 1)

n2

a±bn0

αA0n2

(a, b)n−1

ABn (n − 1)n3n3. a, b ∈ Cn である時，内積 (a, b) は(a, b) =n∑2ai bii=1と計算する。この計算量を求めよ。

5.4

収束判定と停止則について

多くの数値計算アルゴリズムは，最終的には収束する無限数列もしくはベクトル列を漸化式に基づいて反復的に生成するという形で表現される。本書で取り上げるアルゴリズムのうち，有限回の演算で終わるものは連立一次方程式の直接法 (第 7 章) だけであり，実用的な観点から，これをわざわざ無限ベクトル列を生成する形に直したアルゴリズム (第 9 章，第 10 章2 ) も使用される。
その際，どこで計算を打ち切るべきか，その機構をプログラム中に予め埋め込んでおく必要がある。この収束を判定する機構を停止則 (stopping rule) と呼ぶ。実際には収束しない計算を行うことも考えられるため，停止則に加えて最大反復回数を指定し，2 重に保険をかけておくのが普通である。
初等関数の近似計算のように，事前に反復回数が分かっているケースもあるが，多くは数列 (ベクトル列) の値そのものを見て判別する必要がある。どのように判別すべきかは，問題やアルゴリズムの数学的，数値的性質や，実用上の必要に応じてさまざまな方法が考えられる。

5.4.1

数列の場合

まず， sc (, 0) へ単調に収束する無限数列 s0 , s1 , ..., sn , ... を考える。このとき一般項 sn を

sn = sc + S (n) (ここで lim S (n) = 0)n→∞(5.15)と書く。この時 S (n) は単調に 0 に収束することになる。この S (n) を sn の理論誤差 (打ち切り誤差)と呼ぶ。有限桁の浮動小数点数演算においては，打ち切り誤差が ε M 以下になれば，誤差全体としてはそれ以下にすることは不可能と考えてよい。よって，S (n) ≈ ε M になった時には数列の生成を停止する必要がある。この時，|S (n)| > ε M となる n の範囲を打ち切り誤差領域，|S (n)| ≤ ε M となるn の範囲を丸め誤差領域と呼ぶ。ちょうどこの二つの領域の境目となる n において，最小の相対誤差が得られることになる。この n を最適点と呼ぶ。
2 CG 法は理論的には有限のベクトル列で解に収束するが，丸め誤差の混入する近似計算では無限ベクトル列と同じ扱いが要求される。

第5章

52

計算量について一般には最適点における相対誤差になるまで収束させる必要がないことも多く，その場合には要求される相対誤差 (相対許容度 (relative tolerance))εR > 0 以下になった時に停止すればよい。さすれば，十分大きな n に対しては|sn+1 − sn | = |S (n + 1) − S (n)| ≤ 2 max(|S (n + 1)|, |S (n)|) ≈ εR |sc |が成立した時ということになるので，|sn | ≈ |sc | であるから，実用的には|sn+1 − sn | ≤ εR |sn |(5.16)となった時に停止すれば，rE(sn+1 ) ≈ εR と考えられる。
例題 5.4.1

IEEE754 倍精度計算でsn =)in (∑1−2i=0(5.17)を計算する例を考える。この時 sc = 2/3 = 0.666 · · · である。横軸に n を取り，|sn+1 − sn |, |sn |, rE(sn )をプロットしたのが図 5.1 である。ちょうど |sn+1 − sn | が rE(sn ) と連動して小さくなっていること
Example 1
1.E+01
1.E-01
1.E-03
1.E-05
1.E-07
1.E-09
1.E-11
1.E-13
1.E-15
1.E-17打ち切り
誤差領域

丸め誤差領域|s_{n+1} - s_n||s_n|Rel. Error

0102030n405060

図 5.1: 数列 (5.17) の収束状況が分かる。また，十分大きな n に対しては |sn | ≈ |sc | であるから，|sn+1 − sn | ≤ εR |sn | という停止則は，rE(sn ) ≈ εR である所でうまく停止させることができることが分かる。また，εR < ε M となる丸め誤差領域では，ε M 以下の相対誤差となる収束値を得ることが不可能なことも理解できよう。
次に， sc = 0 の場合を考える。この場合は |S (n)| ≈ ε M となる所が最適点となるので，絶対許容度 εA > ε M を使用して|sn+1 − sn | ≤ εA(5.18)となるところで停止する必要がある。これを例で示す。
例題 5.4.20 へ単調に収束する数列)in (∑21−−sn =23i=0(5.19)

5.4. 収束判定と停止則について53を考える。前述の例題と同様に，横軸に n を取り，|sn+1 − sn |, |sn |, rE(sn ) をプロットしたグラフを
図 5.2 に示す。
Example 2
1.E+00
1.E-02打ち切り
誤差領域

1.E-04
1.E-06

丸め誤差領域

1.E-08
1.E-10|s_{n+1} - s_n||s_n|Rel. Error

1.E-12
1.E-14
1.E-160102030n405060

図 5.2: 数列 (5.19) の収束状況

0 に収束する数列であるから，当然十分大きな n に対しては |sn | ≈ 0 となり，他の全ての値も 0に近づいていく。よって絶対許容度 εA が相対誤差の指標となる。これは真値が 0 の時の相対誤差が絶対誤差と同じであることと同義である。
以上の例により，0 を含む収束値を持つ単調数列を停止させるには，(5.16), (5.18) の両方の条件式の機能を持たせる必要がある。よって，ユーザには相対許容度 0 < εR < ε M と絶対許容度 0 < εAの両方を指定させ，|sn+1 − sn | ≤ εR |sn | + εA

(5.20)

が満足した時に停止させるようにすることが望ましい。εA は一般的には前述の通り εA > ε M とすべきだが，0 , sc < ε M であるケースも考えられるので，問題やアルゴリズムの持つ性質を十分見極めた上で，εA < ε M とする場合も考えられる。
最後に，単調でない収束列の場合を述べておく。一般には収束値も，収束の状況も不明であるとあり得る訳で，その際には (5.20) 式の停止則では真に相対誤差に近い値が望めないこともある。
例えば，方程式 f (x) = 0 を満足する解に収束する数列 x0 , x1 , ..., xn , ... が得られた時は，数列の値そのものを見るよりは，いかにその値が方程式を満足しているか，すなわち残差 (residual) f (xn ) がf (xn ) ≈ 0 であるかを調べる方が合理的である。この場合は| f (xn )| < εA(5.21)

を満足していればそこで停止する，という条件を使うべきであろう。あるいは，収束値にある程度近い出発値 s0 が取れるとすれば，出発値よりは εR 程度の残差の改善が見込めるときに停止させる| f (xn )| < εR | f (x0 )|(5.22)という条件も考えられる。もちろんこの場合も，(5.20) と同様にして| f (xn )| < εR | f (x0 )| + εAという形にしておいた方が汎用性が高まることになる。

(5.23)

第5章

54

計算量について5.4.2 ベクトル列，行列列の場合収束する m 次元ベクトル列 s0 , s1 , ..., sn , ... あるいは行列の列についても前述の停止則の議論がそのまま適用できる。行列についてはベクトルの場合と同様に議論できるので，以下ではベクトル列のみ扱う。
べクトルを一つの塊と見てよい場合は，ベクトルノルム (第 8 章参照) を用いて，単調収束の場
合は (5.20) 式より∥sn+1 − sn ∥ ≤ εR ∥sn ∥ + εA

(5.24)となるし，方程式 f(x) = 0 の解に収束する場合に残差を用いるのであれば，(5.23) 式より∥f(xn )∥ ≤ εR ∥f(x0 )∥ + εA

(5.25)としてよい。
(n)(n) T
しかし実用上，sn = [s(n)1 s2 ...sm ] の特定の要素，あるいは全要素が収束条件 (5.20)，あるいは(5.23) を満足していなければならないこともあり得る。ことに要素の絶対値の大きさが格段に異なる場合，一番絶対値の大きな要素の収束だけをチェックしただけでは不十分なこともある。
繰り返しになるが，停止則の選択は，ことにベクトル列の収束判定においては，問題の数学的性質，アルゴリズムの数値的性質，実用上の要求 (精度) の 3 要素を全て吟味して，それに応じたものを選択する必要がある。
問題 5.4.1
行列 A が 1/2 0  ∈ M2 (R)A = 
1 1/3である時，行列の列 S n , Pn , Qn を次のように定義する。

1. S n = I +2. Pn = I +3. Qn =

∑ni=1 A∑nii=1 (1/i!)Ai∑niii=1 ((−1) /i)A

この時，S = limn→∞ S n , P = limn→∞ Pn ，Q = limn→∞ Qn の近似値を，次の条件を満足するように求めよ。また，その停止則に必要な計算量も算出せよ。

1. 全成分の相対誤差が 10−5 以下
2. Frobenius ノルム ((8.10) 式参照)∥ · ∥F に基づく相対誤差が 10−5 以下
3. 1 ノルム ∥ · ∥1 に基づく相対誤差が 10−5 以下5.4. 収束判定と停止則について

55

演習問題
√√
√
1. 複素数 a = π + e −1, b = 2 + −3 について次の問いに答えよ。
(a) 10 進 5 桁の浮動小数点数を用いて，a, b を表わせ。その際発生する実数部，虚数部の相対丸め誤差をそれぞれ求めよ。

(b) 上で求めた近似値をそれぞれ ea, eb とする時，ea +eb, eaeb, ea/eb をなるべく最小の計算量で済むようなアルゴリズムで計算し，10 進 5 桁で表示せよ。その際，計算において発生する相対丸め誤差もそれぞれ求めよ。

(c) 上の計算で必要な計算量を求めよ。
(d) (b) の計算をソフトウェアを用いて実行し，その手順と結果を説明せよ。
√ √2. ベクトル a = [π e 2 3]T ，b = [1/3 1/7 1/9 1/11]T ∈ R4 とする時，次の問いに答えよ。
(a) 10 進 5 桁の浮動小数点数を用いて，a, b を表わせ。その際発生する各成分ごとの相対丸め誤差を求めよ。

(b) 上で求めた近似値をそれぞれ ea, eb とする時，ea +eb, (ea, eb) をなるべく最小の計算量で済むようなアルゴリズムで計算し，10 進 5 桁で表示せよ。その際，計算において発生する相対丸め誤差も各成分ごとに求めよ。

(c) 上の計算に必要な計算量を答えよ。
(d) (b) の計算をソフトウェアを用いて実行し，その手順と結果を説明せよ。
3. 行列 A， B ∈ M2 (R) が次のような成分を持つとき，次の問いに答えよ。
 √
 − 2 π  e π 1√  , B =A = √ 
π2π−eπ2 − e 2(a) 10 進 5 桁の浮動小数点数を用いて， A, B を表わせ。その際発生する各成分ごとの相対丸め誤差を求めよ。
eBe とする時，Ae + B,e A
eB
e をなるべく最小の計算量で済む(b) 上で求めた近似値をそれぞれ A,ようなアルゴリズムで計算し，10 進 5 桁で表示せよ。その際，計算において発生する相対丸め誤差も各成分ごとに求めよ。

(c) 上の計算に必要な計算量を答えよ。
(d) (b) の計算をソフトウェアを用いて実行し，その手順と結果を説明せよ。

参考図書
数値計算にはある程度，アルゴリズムの計算量という観点が必要である。このことは伊理 [13]も自著の中で触れている。大学初年度向けのテキストは多数出版されているが，例えば

第5章

56
アルゴリズムの設計と解析 I, II

A.V.Aho et al./野崎昭弘・野下浩平 訳
サイエンス社
1977 年などを，結果だけでも良いから通読しておくと参考になる。
数値計算も含めた各種のアルゴリズムをスマートに実装したプログラム満載のJava によるアルゴリズム事典
奥村晴彦ほか技術評論社

2003 年は，座右の書として便利である。

計算量について