模擬講義「コンピュータはどうやって計算しているのか?」(改題「数学と計算 ~理論と実践~」)


概要

 数学(=理論)と計算(=実践)は全く性質の異なるもので,計算の効率性や正確性を追及すると数学とは別の方法論が必要になります。この講義では2次方程式や平均の計算を例に,「高性能な計算」をコンピューター上で行う手法を考えていきます。

1.数学=理論,計算=実践

 私の専門は「高性能計算(High Performance Computing)」です。読んで字のごとく,「高性能」な「計算法」を追求する学問の一分野ということになります。まだコンピューターが開発される前,20世紀までは人間があらゆる計算を行っていましたが,21世紀の現在,計算のうち単純なものはコンピューターが担うようになっています。従って,「高性能計算」も「コンピューターにとって」好都合なように計算方法を工夫するものになっています。この模擬講義では現在のコンピューターがどのような工夫のもとに「高性能計算」を実現しているかを具体的な例を挙げて解説します。

 さて,皆さんが思い浮かべる「計算」とはどのようなものでしょう? 「計算高い奴」と言うと,ずるがしこく立ち振る舞うという意味になりますが,ここで取り扱う計算とは「ごく単純な数の操作」とします。一番簡単な「数の操作」は,小学校で習う四則演算(加減乗除)ですね。ここで\(a\), \(b\)は実数とします。

  • たし算(加算)・・・\(a + b\)
  • ひき算(減算)・・・\(a – b\)
  • かけ算(乗算)・・・\(a \times b = ab\)
  • わり算(除算)・・・\(a / b\)

 四則演算だけで実行できる計算といえば,平均値の計算があります。例えば\(n\)個の実数\(c_1, c_2, …, c_n\)の平均値\(\mu\)は,\(n-1\)の加算と1回の除算を組み合わせて
\[ \mu = \frac{1}{n}\sum^n_{i=1} c_i = \frac{c_1 + c_2 + \cdots + c_n}{n} \]
と計算することができます。

 もう少し複雑な計算になると,平方根(\(\sqrt{a}\)),三角関数(\(\sin a, \cos a, \tan a \)),指数関数(\(\exp(x) = e^x\)),対数関数(\(\log a\))なども使う必要が出てきます。平方根を使う計算では,例えば2次方程式
\[ ax^2 + bx + c = 0 \]
の二つの解\(x_1\), \(x_2\)は
\[ x_1 = \frac{-b + \sqrt{b^2 – 4ac}}{2a}, x_2 = \frac{-b – \sqrt{b^2 – 4ac}}{2a} \]
として計算することができます。
 こういう計算は今では普通に電卓やパソコンソフト(電卓ソフトや表計算)で簡単に実行することができます。

 ところで,こういう計算はどうしてこういう式として表わすことができるのでしょうか? 「平均値」の計算はどうして全ての数を足してその個数で割るのでしょうか? 2次方程式の解はどうしてこのような形になるのでしょうか?

 そのような疑問に答えるための学問体系が「数学(Mathematics)」です。

 例えば平均値。解釈の仕方はいろいろ考えられますが,例えばこう考えてはどうでしょうか? いま,背丈の違う二人の人がいるとします。そのうち背の高い方から低い方へ背丈を分けて同じ身長にならしたとすれば,それがこの二人の「身長の平均値」を意味することになります。これが3人でも4人でも5人でも・・・\(n\)人でも,大きい人から小さい人へ背丈と分けて同じ身長に均した時の背丈が「\(n\)人の身長の平均値」であると言うことができます。

 2次方程式の解はどうでしょう? 教科書にはどうしてこういう解の公式が出てくるのか,ちゃんと解説が書いてあったはずですが皆さん覚えていますか?

 2次方程式の両辺を整理して,左側に\(x\)を含む1次式の2乗が来るようにすると
\[ \left(x + \frac{b}{2a}\right)^2 = \frac{b^2 – 4ac}{4a^2} \]
となりますので,
\[ x + \frac{b}{2a} = \pm \frac{\sqrt{b^2 – 4ac}}{2a} \]
から
\[ x = \frac{-b \pm \sqrt{b^2 – 4ac}}{2a} \]
となる訳です。

 こういう,「平均値の解釈」や「2次方程式の解の導出法」も「数学」の一部です。もちろん「平均値の計算」や「2次方程式の解の計算」も「数学」の一部ですが,これらの計算は「解釈」は「導出」の結果得られる「手続き」に過ぎません。計算の仕方が分かったとしても,どういう理由でそのような計算式になったのか,どういう理屈で計算式が導かれたのかは別に考える必要がある訳です。

 現在,コンピューターが担っている「計算」は数学的に導かれた解釈や導出法から得られる単純処理になっています。しかしこの先にもう一つ難しい問題が残っているのです。

 コンピューターはあらかじめ決められた手続きを教えてもらわなければ何も実行しません。手続きを記述した文書はプログラム(ソフトウェア)として人間が書いてやらなければなりません。平均値の計算や2次方程式の解の計算も,式に基づいてどういう手順で計算するのかを,四則演算や関数の組み合わせ方法まで指示しないと,効率的な計算ができないことがあります。このような手続きを「アルゴリズム(algorithm)」と呼びます。

 私が取り組んでいる「高性能計算」は,式を計算するためのアルゴリズムをどのように決定すれば効率が良くなるか,ということを考える学問分野なのです。

2.コンピューターの構造

 コンピュータは,全ての情報を電気信号として保持し,その信号を伝える回路(バス)と信号のやり取りを通じて何らかの処理を行う機器(デバイス)によって構築された金物(ハードウェア)と,ハードウェア内を伝わる電気信号の膨大な積み重ねによって構築された処理内容の蓄積,即ち軟体(ソフトウェア)によって構成されています。ハードウェアは堅い金属やプラスティックのカタマリですから,目で見ることができたり(可視物),触ってその存在を確認することができますが,ソフトウェアは単なる信号ですから,我々人間が感知するためにはデバイスを介して可視化する必要があります。処理を支持する信号を人間や外部の危機とやり取りするためのデバイスを入出力装置(I/O = Input & Output)と呼びます。こうした入出力装置を介しつつ,コンピューターに指令を与えるソフトウェアを作って初めてコンピューターを自在に操ることができるようになります。

basis_computing1

basis_computing2

 この原則は,先にあげた2次方程式の解の計算でも,平均値の計算でも変わりません。例えばC++(シープラスプラス)というコンピューター言語でこの計算を書くと次のようになります。これは計算部分のみの抜粋です(プログラム全体→quadratic_eq.cc)。

 解はx[0], x[1]という名前のメモリ上の領域に計算され格納されます。このプログラムはコンピュータが直接解釈できる形式に変換され(コンパイル),実行時に2次方程式の係数\(a\), \(b\), \(c\)を入力することで解\(x_1\), \(x_2\)のをディスプレイに表示(出力)します。例えば
\[ x^2 + 3x + 2 = 0 \]
の解\(-1, -2\)を計算する場合は次のように実行します。

$ ./quadratic_eq 1 3 2 ← プログラム名,係数の入力
1 * x^2 + (3) * x + (2) = 0 ← 方程式の確認
x[0] = (-1,0) ←解1
x[1] = (-2,0) ←解2

 また,平均を計算するC++プログラム”average.cc“は次のようになります。

 このプログラムをコンパイルして,\((1 + 2 + \cdots + 10) / 10 = 5.5 \)を計算させると,次のような結果が得られます。

$ ./average 10 1 2 3 4 5 6 7 8 9 10 ←プログラム名 データ数 データ1 ・・・ データ10
12
2
3
4
5
6
7
8
9
10
11
1 2 3 4 5 6 7 8 9 10 ←データの確認
average = 5.5 ←平均値

4.浮動小数点数と桁落ち
 さて,例として挙げた2次方程式の解計算と平均値の計算,実はどちらも正確なものではありません。コンピュータ内部での実数は,有限桁の小数でしか表現できないため,どうしても真の値とのズレ(誤差)が発生してしまいます。

 例えば,円周率\(\pi = 3.1415926535…\)は,同じ数のパターンの繰り返しが現れない無理数として良く知られたものですが,正確に表現するためには無限の桁数が必要になります。しかしコンピュータのメモリは有限なので,どこかで打ち切って丸めなければ格納できません。例えば15桁で表現するためには
\[ \pi \approx 3.14159265358979 \]
としなければならず,
\[\pi – 3.14159265358979 = 0.0000000000000032384626433832795… \]
という誤差が生じます。これを丸め誤差と呼びます。例として挙げた二つのプログラムではどちらも約15桁程度の正確さを持っておりますので,16桁目以降は不正確ということになります。

 実用的には15桁も正確であればさほど問題発生しませんが,稀に「桁落ち」という現象が起き,結果が著しく不正確になることがあります。例えば2次方程式の例では
\[ x^2 – 12345x + 0.00001 = 0\]
を計算してみると

$ ./quadratic_eq 1 12345 0.00001
1 * x^2 + (12345) * x + (1e-05) = 0
x[0] = (-8.10359779279679e-10,0)
x[1] = (-12344.9999999992,0)

という結果が得られますが,この場合\(x_1 = -8.100445524504379… \times 10^{-10}\)なので(正確な計算はこちらで実行できます),最初の3桁しか合っていません。誤差を含む,ほとんど同じ大きさの有限小数の加減算を行うとこのような「桁落ち」が起きます。これを防ぐには小数の桁数を増やすか,計算方法を変えるしかありません。

 ちなみに,桁数を約32桁に増やして計算してみると(プログラムはこちら


$ ./quadratic_eq_dd 1 12345 0.00001
1.000000000000000000000000000000e+00 * x^2 + (1.234500000000000000000000000000e+04) * x + (1.000000000000000000000000000000e-05) = 0
x[0] = (-8.100445524504379240908904588516e-10,0.000000000000000000000000000000e+00)
x[1] = (-1.234499999999918995544754956208e+04,0.000000000000000000000000000000e+00)

となり,桁落ち分をある程度フォローできていることが分かります。とはいえ,桁数を増やすより,もう少し賢いやり方はないのでしょうか?

 2次方程式の場合,次のように桁落ちしない正確な解\(x_1’\)を使って次の解\(x_2’\)を求めるという手順を取ることで,桁落ちを回避できることが知られています。
\begin{eqnarray*}
x_1′ &=& \frac{-b – \mbox{sign}(b)\sqrt{b^2 – 4ac}}{2a} \\
x_2′ &=& \frac{c}{2x_1′}
\end{eqnarray*}

 このような改良を加えたプログラムを”quadratic_eq2.cc“として実行してみます。

$ ./quadratic_eq2 1 12345 0.00001
1 * x^2 + (12345) * x + (1e-05) = 0
x[0] = (-12344.9999999992,0)
x[1] = (-8.10044552450438e-10,-0)

 このように,正確な解を得ていることが分かります。

 平均値の計算でも桁落ちが発生します。絶対値の小さい数から足しこむことである程度,桁落ちを防ぐことはできますが,2次方程式のほどの抜本的効果はありませんので,小数の桁数を増やすことが一番の解決策となります。

5.マルチコアと並列化

 コンピュータに求められる能力は増えることはあれ,減ることはありません。人間の利便性を高めるために高性能化は常に求められています。
 コンピュータの性能を高める方法は,CPUの処理能力を上げることです。その一つの方法として,動作周波数(clock)を高めるという方法があります。実際に,20世紀最後の10年ほどで飛躍的に動作周波数は増加しました。

Intel CPUの動作周波数

 しかし,21世紀に入ると打ち止め状態になっています。動作周波数を高めると電力を食うようになります。一般家庭でも使用できるようなパソコンが何千ワットも食うようでは使い物になりません。また,電力消費に伴う発熱も深刻で,安定的に動作させるための空冷ファンは大きくなるばかり。余計に電力消費が増えていきます。
 性能は上げなくてはならない,しかし,電力消費量の増加は極力抑えたい・・・という二律背反な要求に対し,CPUメーカーは電子回路の改良と共に,「マルチコア」という方向にかじを切りました。簡単に言うと,CPUの処理の中核を担う回路(コア)を複数(マルチ)持たせて,あたかも複数のパソコンを同時に使えるようにする技術です。もし6つのコアを持つCPUがあれば,6つの処理を同時にこなすことができ,短時間で終了することで多少電力消費は増えてもトータルでは少なくすることができるはずです。

muticore

 このように,複数の処理を同時並行にこなすことを「並列処理」と呼びます。対して,一つの処理だけを順次行う処理を「逐次処理」と呼びます。現在のパソコンやスマートフォン(スマホ)は2~4コアCPUが普通ですし,GPUのようにもっと大量のコアを備えてグラフィックスや計算を高速にこなす「メニー(many, 多数)コア」ハードウェアも登場しています。

gpu

 では,どうやったら複数のコアを当時に使うことができるのでしょうか? 並列処理を行うプログラム技法はいくつかありますが,一つのプログラムの中の処理を「スレッド」という単位に分割して並列に実行可能にさせる,というマルチスレッド化が良く用いられます。例えば,OpenMPという技法を使って2次方程式の解を同時に計算させるプログラムは次のようになります。

\(x_1, x_2\)の計算は相互に依存していませんので,並行して計算させることができるわけです。逆に,桁落ちが起きないように改良してしまうと,片方の解が出るまでは次の解の計算に進めず,このような並列処理はできなくなります。

 しかし,2次方程式の計算は大変軽いものなので,並列化しても大したメリットはありません。桁落ちを防ぐ処理を行う方が実用的には有用と言えるでしょう。

6.並列処理の効果のほどは?

 では平均の計算はどうでしょう? 大量のデータ(数万,数億,・・・)の平均値の計算は普通に行われるものなので,並列処理するメリットはかなりありそうです。ではどのように行えばいいのでしょうか? 前述した

\[ \frac{1 + 2 + \cdots + 9 + 10}{10} = 5.5 \]

の計算を例に考えてみましょう。

 2コア(デュアルコア)のCPUでは,同時に2スレッドの処理が可能です。その際には,半分ずつ和の計算を行い,次にその部分和を足し合わせて1スレッドで合計値を求め,最後に10で割ります。

parallel2

 4コア(クアッドコア)が使える場合はどうでしょう? 最初の部分和は4スレッドで計算し,次の部分和の計算は2スレッドで分担します。何故なら,最初の部分和の計算よりも計算量が増えると,ここで時間が余計にかかってしまうからです。

 最後は2個の和の計算を1スレッドで行って10で割ります。

parallel4

 この場合,和の計算は最大10回で済みますから,5スレッド以上使える環境があっても意味がありません。

 こうして2スレッド,4スレッドそれぞれ並列実行できたとすると,計算時間は4スレッドの方が短くなることが期待されます。

parallel_time

 では更に大量のデータの平均の計算を実行しようとするとどうでしょうか? データ数\(n\)に対して,同時に実行可能なスレッド数を\(\alpha\) (\(n > \alpha\))とすると最初の部分和の計算は\(\lceil n / \alpha\rceil \)回の和で済みます。ここで\(\lceil x \rceil\)は \(x\)を超えない最大の整数を意味します。

 次の段階で足すべきデータは\( n/2 \)個になっています。これを繰り返すことで,最後は2個の和で済ますことができます。最終的に2個になったデータの和と\(n\)の除算は1スレッドで終了しますので,トータルでは和の計算の部分は最大
\[ \lceil \log_2 n \rceil \]
回の繰り返しで平均の計算が完了することになります。先にあげた10個のデータの場合は\(\lceil\log_2 10\rceil\) \(=\lceil 3.321928…\rceil = 4\)となり,和の計算は4段階必要なことになります。

 このような集約を行う計算をリダクションと呼びます。これをOpenMPで記述すると下記のようになります。(プログラム全体: average_omp.cc)

 使用するスレッド数はomp_set_num_threads関数で設定しています。データの数が少ない時はスレッド数を制限するようにしています。

 では実際に並列処理によってどの程度計算時間が軽減されるかを見ていきましょう。パソコンの環境は次の通りです。

  • CPU: Intel Core i7-3820 3.6GHz 4コア
  • メモリ: 32GB
  • OS: CentOS 6.5 x86_64
  • コンパイラ: Intel C++ 14.0.2
  • データ量: 16GB (\(= 8\mbox{Bytes} \times 2 \times 1024^3\))

4コアありますので,最大4スレッドを並列に実行できることになります。そこで,スレッド数を1, 2, 4として実行してみました。”Time(s)”の所が処理時間(秒)になります。


$ export OMP_NUM_THREADS=4;./average_omp 1
#threads: 1
num_c : 2147483648
average = 1073738907.87731
Time(s) = 0.998027801513672
$ export OMP_NUM_THREADS=4;./average_omp 2
#threads: 2
num_c : 2147483648
average = 1073738907.87725
Time(s) = 0.568727016448975
$ export OMP_NUM_THREADS=4;./average_omp 4
#threads: 4
num_c : 2147483648
average = 1073738907.87729
Time(s) = 0.453949928283691

1スレッドでは約1秒でしたが,2スレッドで約0.57秒,4スレッドで0.45秒で終了しています。2スレッドで大体半減しますが,4スレッドにすると更に半分に,とはなりません。スレッド数が多いとその準備の手間で,理論的な見積りよりも処理時間が余計にかかります。

 ちなみに,太線部分のところがスレッド数によって異なっていることにお気づきでしょうか? 並列化してしまうと,特段指定しない限り,計算の順番は全く保障されないので,丸め誤差の差異が発生しているのです。これは桁落ちしていないのでせいぜい末尾の2桁程度で済んでます。

7.まとめ

 以上,駆け足でしたが,数学と計算の違い,そして計算そのものもきちんと考えることは意外と面倒なことが理解できたことでしょう。今まで述べてきたことのポイントを下記にまとめておきます。

     

  • 理論式をそのまま逐次処理的に計算するにしても,桁落ちが起きないような工夫が可能であれば行う必要があります。
  • マルチコアCPUでは並列処理が可能なので,大量のデータの処理を行うときには威力を発揮できます。しかしそれは並列処理可能な部分に限られます。
  • 並列処理を行って性能向上を図るためには,処理の流れ,即ちアルゴリズムを工夫する必要があります。

 現在の科学技術計算は,シミュレーションの理論化,計算処理,それぞれ分担して専門家が担うようになっています。ここで述べたのはごく簡単な例に基づいた計算処理の部分だけですが,数式をそのまま計算するというだけでは済まない面倒さがあることを覚えて頂ければ幸いです。

京都大学のOffice365サービス移行記録

 本学は本年度(2016年度)よりMicrosoftのOffice365サービスの利用を開始しました。メールが世界どこにいても読み書き返信できるようになり,Office 2016の利用も手軽になるなど良いことづくめのようですが,反面,すべての情報を先方に預けるということでもあり,良いことづくめというわけではありません。

 京都大学は本学より早くOffice365とその前身となるサービスへの移行を開始しました。その貴重な記録が公開されています。ネットワーク管理やOffice365サービスに興味のある方は読む価値があります。

 ネットワーク管理者としては,サーバのメンテナンス業務が楽になる反面,こちらではどうしようもないトラブルへの対処も先方任せになるというジレンマに見舞われることになる訳で,日々胃が痛いことには変わりないなぁという感想を持ちました。

「効率的な手順とは?」 : 2/1 裾野高校模擬講義

1.初めに
 世の中にはいろいろな「作業」があります。楽しい作業,苦しい作業,大変な作業,楽な作業・・・傍から見ているのと実際にその作業を行ってみるのでは全く異なる感想を持つ,なんてことはよくあることです。「計算」なんてその最たるものでしょう。数学と聞くともう数式を見るのも嫌い!,という人は多いですが,小学校の時に習った四則演算(足し算,引き算,掛算,割算)であれば,楽しかったかどうかはともかく,少ない桁数の自然数どうしでしたらそれほど時間をかけることなく実行できるはずです。例えば\[2\times 2\times 2\times 2\times 2 = 2^5 = 32\]なんて計算ができないという人はまずいないでしょう。

 さて「作業」にもイロイロありますが,ここではその開始から終了までに要する「時間」を考えてみたいと思います。作業に要する時間は,どのぐらいの手間がかかるのか,つまり,その作業を行うために必要となる最も簡単かつ短時間で終わる小作業の回数をカウントすることで大体の見積もりを取ることができます。先に挙げた\(2^5\)を計算する作業は,2を掛けるという小作業を4回繰り返すことで32という結果を得て完了することができました。従って,これが\(2^6\), \(2^7\), …, \(2^n\)であっても,掛算をそれぞれ5回,6回,…, \(n\)回実行すれば,計算結果を得ることができ,作業を完了することができるわけです。

 しかしもっとこの作業を効率的にできないものでしょうか? 同じ作業を実行するにしても,初めて経験する人よりは,普通は,その作業を経験したことのある人の方が「効率的」=「短い時間で済む」=「少ない小作業の回数で済ませる」ことができますよね? それは同じ作業を繰り返すことで,「この小作業は不要,ここでこの小作業をすることで効率的になる」ということに気が付くからです。現在はどんな工場でも,どんな事務所でも,どんな会社でも,短い時間で効率的に作業を行うことが求められます。これができないとどうなるか? そう,同じ仕事量でも完了までにいたずらに時間がかかり,いつまでたっても帰れない「ブラック企業」と言われることになってしまうわけです。サービス残業させられて残業代も出ないようなら,短い時間で済ませてさっさと定時退社した方がいいですよね?

 さて,先ほど実行した\(2^n\)の計算ですが,もっと効率的に実行する方法がありますよね? そう,掛算の回数を減らすことはそんなに難しくないですよね? まぁ電卓が手元にあれば,という前提で考えてもらうと,例えば\[2^5 = (2\times 2)^2 \times 2 = 32\]という手順であれば,3回の掛算で済むことになります。\(2^6\)ですと,\[ (2\times 2)^3 = (4\times 4) \times 4 = 64\] となってこれも3回で済みます。・・・と考えていくと,一般に\(2^n\)の場合,\(n\)の2進表記が\(n = (p_k p_{k-1} \cdots p_1 p_0)_2\)であれば,\[2^n = 2^{2^k\cdot p_k} \times 2^{2^{k-1}p_{k-1}} \times \cdots \times 2^{2{p_1}}\times 2^{p_0} \]という手順で済むことが分かります。つまり\(2^{100}\)であれば,99回の掛算が必要になるところを,\(100 = (110010)_2\)より,\[ 2^{100} = 2^{2^6+2^5+2^2} = 2^{64 + 32 + 4} = 2^{64}\times 2^{32} \times 2^4 \]となり,\[ 2^{64} = (((((2^2)^2)^2)^2)^2)^2 \] ですから6回の掛算で済み,この途中結果を覚えておくとすると自動的に\(2^{32}\)と\(2^4\)が得られますので,合計8回の掛算で済むわけです。ちなみに答えは\[\begin{split} 2^{100} &= 18446744073709551616\times 4294967296\times 16\\ &= 1267650600228229401496703205376\end{split} \] となります。

 単純な掛算でも計算回数を減らすことができるわけです。これが膨大な数の掛算を行うとなると,効率的な「手順」=「アルゴリズム」を使わない手はありません。

2.ハノイの塔問題
 ではもう少しゲーム的な「作業」を行ってみましょう。アルゴリズムを考えるときには必ず題材になる有名なパズル「ハノイの塔」というものです。
 次のような積木があるとします。

hanoi_3units

 この積木を一番端の棒に移動させて下さい。但し次の条件を必ず守って下さい。

  • 円盤は1枚ずつしか移動できません。移動先は必ず他の棒になります。
  • 小さい円盤の上に大きな円盤を載せてはいけません。
  • 作業は全ての円盤が一番端の棒に,大きさ順に刺さったところで完了します。

 さて,ではこの3枚のハノイの塔問題を,なるべく「効率的」に解いて下さい。

 解答は下記のYouTube動画にあります。

3.トランプの並べ替え:クイックソート
 日常的に良くある作業として「並べ替え」があります。これを効率的にできないものでしょうか?

 例えばジョーカーの入っていない52枚のトランプをスペード,クローバー,ハート,ダイヤモンドの組に分け,それぞれ番号の小さい順に並べ替える問題を考えてみます。どうすればなるべく早くに並べ替えが完了するでしょうか?

 例えば,下記のYouTube動画に示す「クイックソート」という手順は如何でしょうか?

4.おわりに
 コンピューターは人間と違って莫大なデータ量(ビッグデータ)を扱うために使用されます。従ってちょっとでも非効率的なアルゴリズムを使うと途端に膨大な時間を無駄にすることになりかねません。効率的な手順を考えることは,コンピューターを効率的に使うことにもつながります。作業が早く完了すればそれだけ電気代が節約できますからね。
 我々コンピューターの専門家はこのような「効率的」なアルゴリズムを考えることが重要な仕事の一つなのです。

[模擬講義2] 大規模計算とCPU&GPU


 これは模擬講義「数学と計算 ~理論と実践~」の続きになります。まずそちらを読んで頂いてから,こちらの文章に取り組んで下さい。

1. 大規模計算とは?

 大規模計算とは,現在のコンピュータが実行するに際し,大量のデータをメモリに格納し,計算時間にも数秒,数分,数時間・・・数日・・・費やす必要のある一連の計算処理のことです。
 既にお話しした通り,コンピュータは,処理する対象となるデータを一度メモリに格納し,CPUでそれらを読み込んで処理し,処理結果をメモリに書き戻すという作業を繰り返します。この作業1回分の処理時間はごく短く,マイクロ秒(μs)\(=10^{-6}\)秒単位の世界です。しかし例えば5μsの処理を10回行えば50μs,100回行えば500μs,…と,繰り返しの数\(n\)が大きくなればなるほど,合計処理時間\(5n\) μsは大きくなります。このように\(n\)に比例して時間が増える処理のことを「\(O(n)\)のアルゴリズム」と称します。先に挙げた平均の計算
\[ \mu = \frac{1}{n} \sum^n_{i=1} c_i = \frac{c_1 + c_2 + \cdots + c_n}{n} \]
は典型的な\(O(n)\)のアルゴリズムです。

 しかし,実用的な処理の場合,規模,即ち\(n\)が大きくなるにつれてもっと大きな処理時間を要するものが存在します。少し話が難しくなりますが,以下では連立一次方程式(連立方程式の一種)を見ていくことにしましょう。

2. 連立一次方程式と大規模化

 実数の定数\(a\), \(b\), \(c\), \(d\)に対して,次の二つの方程式を同時に満足する解\(x\), \(y\)を求める,というのが,中学校で初めて習う連立方程式,正式には連立「一次」方程式です。

\[\left\{\begin{array}{ccc} ax + by &=& \alpha \\
cx + dy &=& \beta \end{array}\right. \]

 具体的には例えば

\[\left\{\begin{array}{ccc} 3x + 2y &=& 1 \\
-4x – 2y &=& -2 \end{array}\right. \]

のようになります。この場合,\(x = 1\), \(y = -1\)が解となります。

 連立一次方程式は,様々な場面で登場するとても重要な方程式です。そしてそれだけに,とても大規模なものを扱うことが良く行われています。そのため,大規模化しやすいよう,定数や変数は\(a\)と\(x\)という二つの文字に制限し,大規模化した時には定数や変数に添え字(下付き数字)をつけて区別するように変更します。例えば上の2元2次連立一次方程式は下記のように書き直します。

\[\left\{\begin{array}{ccc} a_{11}x_1 + a_{12}x_2 &=& b_1 \\
a_{21}x_1 + a_{22}x_2 &=& b_2 \end{array}\right. \]

 添え字が異なる文字は違うものであることを意味します。ちょっと目がチカチカしますが,我慢して下さい。

 さて,連立一次方程式の重要なところは,変数の個数と方程式の個数が同じであるところです。2元2次の場合は,それぞれの方程式が直線になっており,2つの直線の交わる交点が解となっている,ということを考えると,解が一つに定まるためには2つの方程式=2本の直線が必要であることは感覚的に理解できるでしょう。

 では変数が3つ,即ち\(x_1\), \(x_2\), \(x_3\)になった時はどうでしょう? この場合,方程式は
\[ ax_1 + bx_2 + cx_3 = d \]
という形になりますので,図形としては平面を意味することになります。つまり「板(いた)」ですね。並行でない「板1」と「板2」が二つ重なるところは直線になりますので,点としては一つに絞れません。更にもう一つ「板3」を持ってくると,既に存在している「板1」と「板2」が交わる「直線12」がこの「板3」とぶつかる交点ができますので,これが解として一つに定まることになります。つまり,一つの解=交点が決まるには,板=方程式は3つ必要になります。従って3つの変数に対しては下記のような3元3次連立一次方程式である必要がある訳です。

\[\left\{\begin{array}{ccc} a_{11}x_1 + a_{12}x_2 + a_{13}x_3 &=& b_1 \\
a_{21}x_1 + a_{22}x_2 + a_{23}x_3 &=& b_2 \\
a_{31}x_1 + a_{32}x_2 + a_{33}x_3 &=& b_3 \\
\end{array}\right. \]

 以下は同様です。変数が4つになれば方程式は4本必要になります。変数が5つになれば方程式は5本,・・・,つまり,変数が\(n\)個あれば方程式は\(n\)本必要になる訳です。それを下記のように表現します。

\[\left\{\begin{array}{cccc} a_{11}x_1 + a_{12}x_2 + \cdots + a_{1n}x_n &=& b_1 \\
a_{21}x_1 + a_{22}x_2 + \cdots + a_{2n}x_n &=& b_2 \\
\cdots & & \cdots \\
a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n &=& b_n \\
\end{array}\right. \]

 さらに簡略化して書きたいときには,係数\(a_{ij}\),変数(解)\(x_i\),右辺の定数項\(b_i\)をバラシて
それぞれのカタマリを,\(A\), \(\mathbf{x}\),\(\mathbf{b}\)と表記し

\[ A = \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{array}\right],\ \mathbf{x} = \left[\begin{array}{c}
x_1 \\
x_2 \\
\vdots \\
x_n
\end{array}\right],\ \mathbf{b} = \left[\begin{array}{c}
b_1 \\
b_2 \\
\vdots \\
b_n
\end{array}\right] \]

以下のように一行の数式で表現します。

\[ A\mathbf{x} = \mathbf{b} \]

ここで\(A\)を係数からなる「行列(matrix)」,\(\mathbf{x}\)を「解ベクトル(vector)」,\(\mathbf{b}\)を「定数ベクトル」と呼びます。最初に例示した2元2次方程式

\[\left\{\begin{array}{ccc} 3x_1 + 2x_2 &=& 1 \\
-4x_1 – 2x_2 &=& -2 \end{array}\right. \]

に対しては
\[ A = \left[\begin{array}{cc} 3 & 2 \\
-4 & -2 \end{array}\right],\ \mathbf{x} = \left[\begin{array}{c}
x_1 \\
x_2
\end{array}\right],\ \mathbf{b} = \left[\begin{array}{c}
1 \\
-2
\end{array}\right] \]
となります。

3.連立一次方程式の解き方

 さて,では大規模な連立一次方程式をどう解いていくか? その方法を説明していきます。詳細については大学学部の「線形代数」という科目で習う内容になるので,そのエッセンスだけを述べることにします。

 まず,簡単な2元2次の例で確認してみましょう。肝心なのは,「自動的に解ける形に方程式を変形する」→「自動的に解いてしまう」という手順を踏むということです。例えば\(u_{21} = 0\)であれば

 \[\left\{\begin{array}{rcl} u_{11}x_1 + u_{12} x_2 &=& y_1 \\
u_{22} x_2 &=& y_2 \end{array}\right. \]

と表現することができますので,下の式から順に\(x_2 = y_2 / u_{22} \)として\(x_2\)を求め,次にそれを上の式に代入して\(x_1 = u_{12} x_2 / u_{11}\)を求めればよいことになります。このように,一番下の式が変数が一つだけの方程式になった格好の連立一次方程式は,係数行列\(A\)が上三角行列\(U\),即ち
\[ U = \left[\begin{array}{cc}
u_{11} & u_{12} \\
0 & u_{22}
\end{array}\right] \]
を用いて
\begin{equation}
U\mathbf{x} = \mathbf{y} \label{eqn:upper_triangle}
\end{equation}
と表現することができます。このように上三角行列を係数として持つ連立一次方程式を下の方から解いていく方法を「後退代入」と呼びます。

 同様に\(l_{11} = l_{22} = 1\)かつ\(l_{12} = 0\)であれば

\[\left\{\begin{array}{rcl} y_1 \phantom{+ u_{12}y_2} &=& b_1 \\
l_{21}y_1 + y_2 &=& b_2 \end{array}\right. \]

と表現できますので,これも上から順に\(y_1 = b_1\),次に\(y_2 = b_2 – l_{21}y_1\)として解を求めることができます。このとき係数行列は下三角行列\(L\),即ち
\[ L = \left[\begin{array}{cc}
1 & 0 \\
l_{21} & 1
\end{array}\right] \]
を用いて
\begin{equation}
L \mathbf{y} = \mathbf{b}
\end{equation}
と表現することができます。このような下三角行列を係数として持つ連立一次方程式を上の方から解いていく方法を「前進代入」と呼びます。

 結論を先に述べると,上記2つの形の連立一次方程式を導いて,次のようにして解きます。

  1. [LU分解]  係数行列\(A\)から下三角行列\(L\)と上三角行列\(U\)を求める(実は\(A = LU\)という関係になります)
  2. [前進代入] \(L\mathbf{y} = \mathbf{b}\)を\(\mathbf{y}\)について解く
  3. [後退代入] \(U\mathbf{x} = \mathbf{y}\)を\(\mathbf{x}\)について解く

 実際に
\[\left\{\begin{array}{ccc} 3x_1 + 2x_2 &=& 1 \\
-4x_1 – 2x_2 &=& -2 \end{array}\right. \]
でLU分解→前進代入→後退代入を行っていきます。

 まずLU分解です。
\[ A = \left[\begin{array}{cc}
3 & 2 \\
-4 & -2
\end{array}\right] \]
を,行(横一列)をひとまとめに足したり引いたり定数倍したりして,上三角行列\(U\)にします。
\[ \left[\begin{array}{cc}
3 & 2 \\
-4 & -2
\end{array}\right] \Rightarrow \left[\begin{array}{cc}
3 & 2 \\
-4 – (-4)/3\times 3 & -2 – (-4)/3\times 2
\end{array}\right] = \left[\begin{array}{cc}
3 & 2\\
0 & 2/3 \\
\end{array}\right] = U \]
第1行目に\((-4)/3\)を掛けて,第2行目から引いているわけです。こうして\(a_{21}\Rightarrow 0\)にすることができました。実はこの時,下三角行列\(L\)も自動的に求められており,
\[ L = \left[\begin{array}{cc}
1 & 0 \\
-4/3 & 1
\end{array}\right] \]
になります。実は\(a_{21} / a_{11} = -4 / 3\)を計算しただけになります。

 こうして出てきた\(L, U\)を用いて
\[ L\mathbf{y} = \mathbf{b} \Longleftrightarrow \left[\begin{array}{cc}
1 & 0 \\
-4/3 & 1
\end{array}\right] \left[\begin{array}{c}
y_1 \\
y_2
\end{array}\right] = \left[\begin{array}{c}
1 \\
-2
\end{array}\right] \]
を解いて\[ \mathbf{y} = \left[\begin{array}{c}
1 \\
-2/3
\end{array}\right] \] を求め,次に
\[ U\mathbf{x} = \mathbf{y} \Longleftrightarrow \left[\begin{array}{cc}
3 & 2 \\
0 & 2/3
\end{array}\right] \left[\begin{array}{c}
x_1 \\
x_2
\end{array}\right] = \left[\begin{array}{c}
1 \\
-2/3
\end{array}\right] \]
を解いて\[ \mathbf{x} = \left[\begin{array}{c}
1 \\
-1
\end{array}\right] \] を求めます。

 この方法は\(n\)元\(n\)次連立一次方程式になっても変わりません。しかし,計算の手間が膨大に増えていきます。実際,LU分解は
\[ A = \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21} & a_{22} & \cdots & a_{2n} \\
\vdots & \vdots & & \vdots \\
a_{n1} & a_{n2} & \cdots & a_{nn}
\end{array}\right]\]
\[\Rightarrow \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a_{21}/a_{11} & a_{22} – a_{21}/a_{11}\times a_{12} & \cdots & a_{2n} – a_{21}/a_{11}\times a_{1n} \\
\vdots & \vdots & & \vdots \\
a_{n1}/a_{11} & a_{n2} – a_{n1}/a_{11}\times a_{12} & \cdots & a_{nn} – a_{n1}/a_{11}\times a_{1n}
\end{array}\right]\]
\[\Rightarrow \cdots \Rightarrow \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
a^{(1)}_{21} & a^{(2)}_{22} & \cdots & a^{(2)}_{2n} \\
\vdots & \vdots & & \vdots \\
a^{(1)}_{n1} & a^{(2)}_{n2} & \cdots & a^{(n-1)}_{nn}
\end{array}\right]\]
となります。この過程を図示すると下記のようになります。

lu

黒く塗りつぶした四角部分は
\[ a^{(i-1)}_{ji} / a^{(i-1)}_{ii}\ (j = i + 1, …, n) \]
の計算を行い,灰色の部分は
\[ a^{(i-1)}_{jk} – a^{(i-1)}_{ji} / a^{(i-1)}_{ii} \times a^{(i-1)}_{ik}\ (j, k = i+1, …, n) \]
の計算を行っています。

黒い四角部分は\[ (n-1) + (n-2) + \cdots + 2 + 1 \]個となり,灰色の部分は\[ (n-1)^2 + (n-2)^2 + \cdots + 2 + 1 \]個となり,それぞれ1回ずつの割算と掛算が入っていますので,この割算と掛算の回数だけでも
\[ \sum^{n-1}_{i=1} i(i + 1) = \frac{n(n^2 – 1)}{3} \]
となりますから,最終的には\(O(n^3)\)の手間がかかるようになります。\(n\)が10倍ずつ増えていくとき,\(O(n)\)の計算が10秒,100秒,1000秒・・・と増えていくとすると,\(O(n^3)\)は1000秒,1000000秒, 1000000000秒・・・のように増えていきます。計算時間の増え方が1000倍になっていくわけです。

 こうして\(L, U\)は
\[ L = \left[\begin{array}{cccc}
1 & 0 & \cdots & 0 \\
a^{(1)}_{21} & 1 & \cdots & \vdots \\
\vdots & \ddots & \ddots & 0 \\
a^{(1)}_{n1} & \cdots & a^{(n-1)}_{n,n-1} & 1
\end{array}\right],\ U = \left[\begin{array}{cccc}
a_{11} & a_{12} & \cdots & a_{1n} \\
0 & a^{(2)}_{22} & \cdots & a^{(2)}_{2n} \\
\vdots & \ddots & & \vdots \\
0 & \cdots & 0 & a^{(n-1)}_{nn}
\end{array}\right] \]
として得られますので,前進代入は\(L\mathbf{y} = \mathbf{b}\)を解いて
\begin{align*}
y_1 &:= b_1 \\
y_2 &:= b_2 – a^{(1)}_{21} y_1 \\
& \vdots \\
y_{n-1} &:= b_{n-1} – \sum^{n-2}_{j=1} a^{(j)}_{n-1,j}y_j \\
y_n &:= b_n – \sum^{n-1}_{j=1} a^{(j)}_{nj}y_j
\end{align*}
となり,後退代入は\(U \mathbf{x} = \mathbf{y}\)を解いて
\begin{align*}
x_n &:= \frac{y_n}{a^{(n-1)}_{nn}} \\
x_{n-1} &:= \frac{y_{n-1} – a^{(n-2)}_{n-1,n} x_n}{a^{(n-2)}_{n-1, n-1}} \\
& \vdots \\
x_2 &:= \frac{y_2 – \sum^{n}_{j=3} a^{(1)}_{2j} x_j}{a^{(1)}_{22}} \\
x_1 &:= \frac{y_1 – \sum^{n}_{j=2} a_{1j} x_j}{a_{11}}
\end{align*}
として計算することができます。

 従って,前進代入,後退代入の手間は\(O(n^2)\)になります。

4. 並列で計算するには?

 さて,このように大規模化した連立一次方程式を現在のコンピュータ上で高速に解くためには,LU分解,前進代入,後退代入の計算を並列化する必要があります。

 LU分解は下記のように並列化します。色分けした部分を各スレッドに分配して並列に計算しています。

parallel_lu

 前進代入,後退代入も下記のように並列化します。

parallel_subst

 LU分解も前進代入,後退代入も,計算が進むにつれて計算量が減り,結果的にだんだん並列化効率が落ちていきますが,これ以上の良い手段はなかなか見つからないというのが現状です。しかしこのように並列化することで,マルチコアCPU環境・GPU環境では高速に計算することができるようになります。実際に12コアCPUとTesla K20というGPUを積んだワークステーション環境で解いた結果が下記になります。

dsgesv_cpu_gpu

 最大\(n = 8192\)までの規模の連立一次方程式を説いた時の計算時間と計算に伴う丸め誤差をプロットしています。GPUの方が圧倒的に高速に計算できていることが分かります。

5. まとめ

 駆け足で解説してきましたが,ここで説明した大規模な連立一次方程式をLU分解・前進代入・後退代入を使って解く方法は直接法(direct method)と呼ばれ,Linpack(リンパック)ベンチマークテストと呼ばれる,スーパーコンピュータの速度比較プログラムに使われています。Top500というサイトには,毎年半期ごとに世界中のスーパーコンピュータのLinpackテストの結果に基づいたランキングが掲載されており,2015年9月現在のナンバーワンは中国にあります。日本では神戸にある「京(ケイ)」というスーパーコンピュータが最高速で,現在は世界第4位に後退しました。

 コンピューターの性能評価のためのLinpackテストは,皆さんが中学校で習う連立一次方程式の拡大版で行われており,それらは\(O(n^3)\)の処理時間を要するものである,ということを理解してもらえば幸いです。