概要
数学（＝理論）と計算（＝実践）は全く性質の異なるもので，計算の効率性や正確性を追及すると数学とは別の方法論が必要になります。この講義では２次方程式や平均の計算を例に，「高性能な計算」をコンピューター上で行う手法を考えていきます。
１．数学＝理論，計算＝実践
私の専門は「高性能計算(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)」と呼びます。
私が取り組んでいる「高性能計算」は，式を計算するためのアルゴリズムをどのように決定すれば効率が良くなるか，ということを考える学問分野なのです。
２．コンピューターの構造
コンピュータは，全ての情報を電気信号として保持し，その信号を伝える回路（バス）と信号のやり取りを通じて何らかの処理を行う機器（デバイス）によって構築された金物（ハードウェア）と，ハードウェア内を伝わる電気信号の膨大な積み重ねによって構築された処理内容の蓄積，即ち軟体（ソフトウェア）によって構成されています。ハードウェアは堅い金属やプラスティックのカタマリですから，目で見ることができたり（可視物），触ってその存在を確認することができますが，ソフトウェアは単なる信号ですから，我々人間が感知するためにはデバイスを介して可視化する必要があります。処理を支持する信号を人間や外部の危機とやり取りするためのデバイスを入出力装置（I/O = Input & Output）と呼びます。こうした入出力装置を介しつつ，コンピューターに指令を与えるソフトウェアを作って初めてコンピューターを自在に操ることができるようになります。
この原則は，先にあげた２次方程式の解の計算でも，平均値の計算でも変わりません。例えばC++(シープラスプラス)というコンピューター言語でこの計算を書くと次のようになります。これは計算部分のみの抜粋です（プログラム全体→quadratic_eq.cc)。
// print the given equation
cout << a << " * x^2 + (" << b << ") * x + (" << c << ") = 0 " << endl;
// Calculate solutions with analytic formulas
d = b * b - 4.0 * a * c;
sqrt_d = sqrt((complex<double>)d);
a2 = 2.0 * a;
x[0] = (- b + sqrt_d) / a2;
x[1] = (- b - sqrt_d) / a2;
// Print solutions
cout << "x[0] = " << setprecision(15) << x[0] << endl;
cout << "x[1] = " << setprecision(15) << x[1] << endl;
解は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“は次のようになります。
// print the given datum for(i = 0; i < num_c; i++) cout << c[i] << " "; cout << endl; // Calculate average average_mu = 0.0; for(i = 0; i < num_c; i++) average_mu += c[i]; average_mu /= num_c; // Print solutions cout << "average = " << setprecision(15) << average_mu << endl;
このプログラムをコンパイルして，\((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 ←平均値
４．浮動小数点数と桁落ち
さて，例として挙げた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次方程式のほどの抜本的効果はありませんので，小数の桁数を増やすことが一番の解決策となります。
５．マルチコアと並列化
コンピュータに求められる能力は増えることはあれ，減ることはありません。人間の利便性を高めるために高性能化は常に求められています。
コンピュータの性能を高める方法は，CPUの処理能力を上げることです。その一つの方法として，動作周波数(clock)を高めるという方法があります。実際に，20世紀最後の10年ほどで飛躍的に動作周波数は増加しました。
しかし，21世紀に入ると打ち止め状態になっています。動作周波数を高めると電力を食うようになります。一般家庭でも使用できるようなパソコンが何千ワットも食うようでは使い物になりません。また，電力消費に伴う発熱も深刻で，安定的に動作させるための空冷ファンは大きくなるばかり。余計に電力消費が増えていきます。
性能は上げなくてはならない，しかし，電力消費量の増加は極力抑えたい・・・という二律背反な要求に対し，CPUメーカーは電子回路の改良と共に，「マルチコア」という方向にかじを切りました。簡単に言うと，CPUの処理の中核を担う回路（コア）を複数（マルチ）持たせて，あたかも複数のパソコンを同時に使えるようにする技術です。もし6つのコアを持つCPUがあれば，6つの処理を同時にこなすことができ，短時間で終了することで多少電力消費は増えてもトータルでは少なくすることができるはずです。
このように，複数の処理を同時並行にこなすことを「並列処理」と呼びます。対して，一つの処理だけを順次行う処理を「逐次処理」と呼びます。現在のパソコンやスマートフォン（スマホ）は2～4コアCPUが普通ですし，GPUのようにもっと大量のコアを備えてグラフィックスや計算を高速にこなす「メニー(many, 多数)コア」ハードウェアも登場しています。
では，どうやったら複数のコアを当時に使うことができるのでしょうか？ 並列処理を行うプログラム技法はいくつかありますが，一つのプログラムの中の処理を「スレッド」という単位に分割して並列に実行可能にさせる，というマルチスレッド化が良く用いられます。例えば，OpenMPという技法を使って2次方程式の解を同時に計算させるプログラムは次のようになります。
#pragma omp parallel sections
{
// thread 0
#pragma omp section
{
x[0] = (- b + sqrt_d) / a2;
}
// thread 1
#pragma omp section
{
x[1] = (- b - sqrt_d) / a2;
}
}
\(x_1, x_2\)の計算は相互に依存していませんので，並行して計算させることができるわけです。逆に，桁落ちが起きないように改良してしまうと，片方の解が出るまでは次の解の計算に進めず，このような並列処理はできなくなります。
しかし，2次方程式の計算は大変軽いものなので，並列化しても大したメリットはありません。桁落ちを防ぐ処理を行う方が実用的には有用と言えるでしょう。
６．並列処理の効果のほどは？
では平均の計算はどうでしょう？ 大量のデータ（数万，数億，・・・）の平均値の計算は普通に行われるものなので，並列処理するメリットはかなりありそうです。ではどのように行えばいいのでしょうか？ 前述した
\[ \frac{1 + 2 + \cdots + 9 + 10}{10} = 5.5 \]
の計算を例に考えてみましょう。
2コア（デュアルコア）のCPUでは，同時に２スレッドの処理が可能です。その際には，半分ずつ和の計算を行い，次にその部分和を足し合わせて1スレッドで合計値を求め，最後に10で割ります。
4コア（クアッドコア）が使える場合はどうでしょう？ 最初の部分和は4スレッドで計算し，次の部分和の計算は2スレッドで分担します。何故なら，最初の部分和の計算よりも計算量が増えると，ここで時間が余計にかかってしまうからです。
最後は2個の和の計算を1スレッドで行って10で割ります。
この場合，和の計算は最大10回で済みますから，5スレッド以上使える環境があっても意味がありません。
こうして2スレッド，4スレッドそれぞれ並列実行できたとすると，計算時間は4スレッドの方が短くなることが期待されます。
では更に大量のデータの平均の計算を実行しようとするとどうでしょうか？ データ数\(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)
// Calculate average num_threads = 4; if(num_threads > num_c) num_threads = num_c / 2; omp_set_num_threads(num_threads); cout << "#threads: " << omp_get_max_threads() << endl; cout << "num_c : " << num_c << endl; #pragma omp parallel for reduction(+:average_mu) for(i = 0; i < num_c; i++) average_mu += c[i]; average_mu /= num_c;
使用するスレッド数は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桁程度で済んでます。
７．まとめ
以上，駆け足でしたが，数学と計算の違い，そして計算そのものもきちんと考えることは意外と面倒なことが理解できたことでしょう。今まで述べてきたことのポイントを下記にまとめておきます。
- 理論式をそのまま逐次処理的に計算するにしても，桁落ちが起きないような工夫が可能であれば行う必要があります。
- マルチコアCPUでは並列処理が可能なので，大量のデータの処理を行うときには威力を発揮できます。しかしそれは並列処理可能な部分に限られます。
- 並列処理を行って性能向上を図るためには，処理の流れ，即ちアルゴリズムを工夫する必要があります。
現在の科学技術計算は，シミュレーションの理論化，計算処理，それぞれ分担して専門家が担うようになっています。ここで述べたのはごく簡単な例に基づいた計算処理の部分だけですが，数式をそのまま計算するというだけでは済まない面倒さがあることを覚えて頂ければ幸いです。