これは模擬講義「数学と計算 ～理論と実践～」の続きになります。まずそちらを読んで頂いてから，こちらの文章に取り組んで下さい。
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つの形の連立一次方程式を導いて，次のようにして解きます。
- [LU分解] 係数行列\(A\)から下三角行列\(L\)と上三角行列\(U\)を求める(実は\(A = LU\)という関係になります)
- [前進代入] \(L\mathbf{y} = \mathbf{b}\)を\(\mathbf{y}\)について解く
- [後退代入] \(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]\]
となります。この過程を図示すると下記のようになります。
黒く塗りつぶした四角部分は
\[ 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 \]個となり，それぞれ１回ずつの割算と掛算が入っていますので，この割算と掛算の回数だけでも
\[ \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分解は下記のように並列化します。色分けした部分を各スレッドに分配して並列に計算しています。
前進代入，後退代入も下記のように並列化します。
LU分解も前進代入，後退代入も，計算が進むにつれて計算量が減り，結果的にだんだん並列化効率が落ちていきますが，これ以上の良い手段はなかなか見つからないというのが現状です。しかしこのように並列化することで，マルチコアCPU環境・GPU環境では高速に計算することができるようになります。実際に12コアCPUとTesla K20というGPUを積んだワークステーション環境で解いた結果が下記になります。
最大\(n = 8192\)までの規模の連立一次方程式を説いた時の計算時間と計算に伴う丸め誤差をプロットしています。GPUの方が圧倒的に高速に計算できていることが分かります。
5. まとめ
駆け足で解説してきましたが，ここで説明した大規模な連立一次方程式をLU分解・前進代入・後退代入を使って解く方法は直接法(direct method)と呼ばれ，Linpack（リンパック）ベンチマークテストと呼ばれる，スーパーコンピュータの速度比較プログラムに使われています。Top500というサイトには，毎年半期ごとに世界中のスーパーコンピュータのLinpackテストの結果に基づいたランキングが掲載されており，2015年9月現在のナンバーワンは中国にあります。日本では神戸にある「京（ケイ）」というスーパーコンピュータが最高速で，現在は世界第4位に後退しました。
コンピューターの性能評価のためのLinpackテストは，皆さんが中学校で習う連立一次方程式の拡大版で行われており，それらは\(O(n^3)\)の処理時間を要するものである，ということを理解してもらえば幸いです。