この記事は機械学習に必要な高校数学やり直しアドベントカレンダー Advent Calendar 2016の17日目の記事です。
技術系の開発をしていると、携わる分野によって使う数学はだいたい決まってしまいます。仕事で使わない数学は時間とともにだんだん忘れ去っていきます。そんなシニアなエンジニアが、高校で習った数学をなんとか思い出しつつ、線形代数の初歩を書いてみます。
さて、高校数学で習う線形代数って何かと考えた時、連立方程式やベクトル、行列ではないかと思います。この中で、連立一次方程式の解法について書いてみます。なお、ここでは解があり、簡単に解ける問題を例にします。
消去法
まず、最初に習うのが消去法で解く方法だと思います（もしかしたら中学かも）。消去法は連立方程式に対して以下の操作を繰り返し行って、解を計算します。
- 1つの式に0でない定数を掛ける
- 1つの式に式を加える
- 2つの式を入れ替える
この3つの操作を行基本変形と呼びます。では、例題を使って実際に計算してみます。
\left\{%
\begin{align}
\; x + y + 2z&=9&\cdots (1) \\
\; x + 2y + z&=11&\cdots (2) \\
\; 2x + y + z&=8&\cdots (3)
\end{align}
\right.
を使って計算してみます。まず、(2)-(1)を計算します。
\begin{equation}
y-z=2\;\;\;\;\cdots (4)
\end{equation}
次に(3)-2×(1)から
\begin{equation}
-y-3z=-10\;\;\;\;\cdots (5)
\end{equation}
そして、(5)＋(4)で
\begin{equation}
-4z=-8\;\;\;\;\cdots (6)
\end{equation}
x, y, zの係数が1になるように式(5), (6)に定数をかけると、元の連立方程式は
\left\{%
\begin{align}
\; x + y + 2z&=9&\cdots&(1)\\
\; y + 3z&=10&\cdots&(5)'\\
\; z&=2&\cdots&(6)'
\end{align}
\right.
となります。ここでzの値は求まってますので、(6)'と(5)'でyが求まり、さらに、y、zを(1)に代入してxが求まります。答えは以下になります。
\left\{%
\begin{eqnarray}
\; x=1\\
\; y=4\\
\; z=2
\end{eqnarray}
\right.
消去法は式(1), (5)', (6)'を見ていただければ分かるように、各式から変数を一つずつ消していき、次に、下から上に向かって変数に求まった値を代入していき、解を求めるものです。
ずっと昔、面倒だなぁと思いながら問題集を解いていたのを思い出しました。
行列
行列はもともと連立方程式を扱うために考え出されたものなので、この連立方程式も行列で表してみます。連立方程式を行列の乗算を使って表すと次のようになります。
\begin{equation}
\left(
\begin{matrix}
1 & 1 & 2 \\
1 & 2 & 1 \\
2 & 1 & 1
\end{matrix}
\right)\left(
\begin{matrix}
x \\
y \\
z
\end{matrix}
\right)
=
\left(
\begin{matrix}
9 \\
11 \\
8
\end{matrix}
\right)
\end{equation}
ここで、
\begin{equation}
A = \left(
\begin{matrix}
1 & 1 & 2 \\
1 & 2 & 1 \\
2 & 1 & 1
\end{matrix}
\right)
\;\;x = \left(
\begin{matrix}
x \\
y \\
z
\end{matrix}
\right)
\;\;b =\left(
\begin{matrix}
9 \\
11 \\
8
\end{matrix}
\right)
\end{equation}
とすると、元の連立方程式は
\begin{equation}
Ax = b
\end{equation}
と表されます。なお、Aは係数行列、xは未知数ベクトル、bは定数ベクトルと呼ばれます。
Gauss-Jordanの消去法（掃き出し法）
連立方程式を行列で表現したので、行列のままで解を計算してみます。まずは掃き出し法を使って計算してみます。方法は上で説明した消去法と同じですが、行列の時は掃き出し法と呼ばれることが多いです。
計算を行うために、まず、係数行列Aと右辺の定数ベクトルを一つに合わせた係数拡大行列
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 1 & 2 & 9\\
1 & 2 & 1 & 11\\
2 & 1 & 1 & 8
\end{array}
\right)
\end{equation}
を作ります。この係数拡大行列の縦線より左側が単位行列になれば、定数ベクトル部分が解になります。つまり、
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 0 & 0 & ?\\
0 & 1 & 0 & ?\\
0 & 0 & 1 & ?
\end{array}
\right)
\end{equation}
と変形できれば、?の部分が解になります。それでは計算をしていきます。
まず、1行目を使って2行目、3行目の1列目が0になるようにします。つまり、
\begin{align}
2行目&=2行目 - 1行目\\
3行目&=3行目 - 1行目\times2
\end{align}
の計算を行います。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 1 & 2 & 9\\
0 & 1 & -1 & 2\\
0 & -1 & -3 & -10
\end{array}
\right)
\end{equation}
同様に2行目を使って1行目と3行目の2列目を0にします。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 0 & 3 & 7\\
0 & 1 & -1 & 2\\
0 & 0 & -4 & -8
\end{array}
\right)
\end{equation}
3行目の-4が1になるように3行目を-4で割ります。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 0 & 3 & 7\\
0 & 1 & -1 & 2\\
0 & 0 & 1 & 2
\end{array}
\right)
\end{equation}
最後の3行目を使って、1行目と2行目の3列目を0にします。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 0 & 0 & 1\\
0 & 1 & 0 & 4\\
0 & 0 & 1 & 2
\end{array}
\right)
\end{equation}
これで縦棒より右側に解が求まりました。元の式の形に戻すと、
\begin{equation}
\left(
\begin{matrix}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1
\end{matrix}
\right)\left(
\begin{matrix}
x \\
y \\
z
\end{matrix}
\right)
=
\left(
\begin{matrix}
1 \\
4 \\
2
\end{matrix}
\right)
\end{equation}
となります。これが高校数学での標準的な解き方ということです。
Gaussの消去法
高校数学ではまた別の解き方も習うはずです。それがGaussの消去法です。先ほどのGauss Jordanとよく似ていますが、先ほどの問題をGaussの消去法でも解いてみます。
Gauss Jordanと同様に係数拡大行列を作ります。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 1 & 2 & 9\\
1 & 2 & 1 & 11\\
2 & 1 & 1 & 8
\end{array}
\right)
\end{equation}
Gauss Jordanと違うのは、対角要素より下の要素のみを0にしていくところです。実際にやってみます。
まず、1行目を使って、2行目と3行目の1列目が0になるようにします。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 1 & 2 & 9\\
0 & 1 & -1 & 2\\
0 & -1 & -3 & -10
\end{array}
\right)
\end{equation}
次に2行目を使って3行目の2列目を0にします。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 1 & 2 & 9\\
0 & 1 & -1 & 2\\
0 & 0 & -4 & -8
\end{array}
\right)
\end{equation}
最後に3行目の3列目が1になるよう、3行目を-4で割ります。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 1 & 2 & 9\\
0 & 1 & -1 & 2\\
0 & 0 & 1 & 2
\end{array}
\right)
\end{equation}
これでまずzの値が戻りました。ここで縦線より左側の係数行列を見ると以下のように対角線より下は全部0の行列なっています。
\begin{equation}
\left(
\begin{matrix}
1 & 1 & 2 \\
0 & 1 & -1\\
0 & 0 & 1
\end{matrix}
\right)
\end{equation}
こののような行列を上三角行列と呼び、ここまでの操作を前進消去と呼びます。さて、計算に戻りましょう。今度は逆に下から上に向かって、対角線上の要素だけが1となり、それ以外は0となるように操作していきます。まず、2行目に3行目を足します。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 1 & 2 & 9\\
0 & 1 & 0 & 4\\
0 & 0 & 1 & 2
\end{array}
\right)
\end{equation}
最後に1行目から2行目と3行目を2倍したものを引きます。
\begin{equation}
\left(
\begin{array}{ccc|c}
1 & 0 & 0 & 1\\
0 & 1 & 0 & 4\\
0 & 0 & 1 & 2
\end{array}
\right)
\end{equation}
これでGauss Jordanと同じ解が求まりました。
数十年前にせっせとこういう計算をしていたのを思い出してきました。すっかり忘れたと思ってましたが、意外と思い出すものです。
プログラム
ここからは高校数学の範囲を（たぶん）超えますが、アルゴリズムが分かるとプログラムしてみたくなるので、JavaScriptでコードを書いてみました。
Gauss Jordanの消去法
コードは以下のようになります。繰り返しますが、単純に解ける場合しか考慮していません。実用にするにはpivotが0の場合の処理とか、解が求まらない場合の対処、誤差を減らす工夫などが必要になります。
function gaussJordan(a, c) {
var i, j, k, q;
var n = a.length - 1 ;
for (k = 0; k <= n; k++) {
for (i = 0; i <= n; i++) {
if (i === k) continue;
q = a[i][k] / a[k][k];
for (j = k + 1; j <= n; j++) {
a[i][j] -= q * a[k][j];
}
c[i] -= q * c[k];
}
}
for (k = 0; k <= n; k++) {
c[k] /= a[k][k];
}
}
手計算のやり方をそのままコードにしているので、手計算と比べていただければわかると思います。
Gaussの消去法
次はGaussの消去法のプログラムです。
function gauss(a, c) {
var i, j, k, q, s;
var n = a.length - 1 ;
// Forward elimination
for (k = 0; k < n; k++) {
for (i = k + 1; i <= n; i++) {
q = a[i][k] / a[k][k];
for (j = k + 1; j <= n; j++) {
a[i][j] -= q * a[k][j];
}
c[i] -= q * c[k];
}
}
// Back subsitution
c[n] /= a[n][n];
for (k = n - 1; k >= 0; k--) {
s = c[k];
for (j = k + 1; j <= n; j++) {
s -= a[k][j] * c[j];
}
c[k] = s / a[k][k];
}
}
前進消去と後退代入の2ステップあるので、コードはGauss Jordanより長くなっています。
計算量
せっかくプログラムまで作ったので、それぞれの計算量を見積もって比較してみます。それによってプログラムに使うにはどちらの解法が良いかを知る一つの基準になります。
Gauss Jordanの計算量
まず一番内側のループは、
for (j = k + 1; j <= n; j++) {
a[i][j] -= q * a[k][j];
}
です。このループの中で減算が1回、乗算が1回行われます。このループでの演算回数は、
\begin{equation}
乗除算:
\sum_{j=k+1}^{n} 1 = (n - k)\\
加減算:
\sum_{j=k+1}^{n} 1 = (n - k)
\end{equation}
になります。次に、その一つ外側のループをみます。
for (i = 0; i <= n; i++) {
if (i === k) continue;
q = a[i][k] / a[k][k];
for (j = k + 1; j <= n; j++) {
a[i][j] -= q * a[k][j];
}
c[i] -= q * c[k];
}
このループでは、先ほどのループに加え乗算を1回、除算を1回、減算を1回行います。したがって、
計算量は、
\begin{equation}
乗除算:
\sum_{i=1}^{n-1} (n - k + 2) = (n - 1)(n - k + 2)\\
加減算:
\sum_{i=1}^{n-1} (n - k + 1) = (n - 1)(n - k + 1)
\end{equation}
一番外側のループは、内側のループが入っているだけなので、計算量は、
\begin{equation}
乗除算:
\sum_{k=1}^{n} (n - 1)(n - k + 2) \\
加減算:
\sum_{k=1}^{n} (n - 1)(n - k + 1)
\end{equation}
となります。さらに、もう一つ
for (k = 0; k <= n; k++) {
c[k] /= a[k][k];
}
のループがあります。こののループでは除算が1回行われますので、このループの計算量は、
\begin{equation}
乗除算:
\sum_{k=1}^{n} 1 = n
\end{equation}
となります。まず、乗除算分を計算します。l=n-kとすると
\begin{align}
\sum_{k=1}^{n} (n - 1)(n - k + 2) + n
&= (n - 1)\sum_{l=0}^{n-1} (l+2) + n\\
&= (n -1)(\frac{n(n-1)}{2} +2n) + n \\
&= \frac{n(n-1)^2}{2} + 2n(n-1)+n \\
&= \frac{n(n-1)^2}{2} + n(2n-1)\\
&= n\frac{(n-1)^2 + 4n - 1}{2} \\
&= \frac{n^2(n+2)}{2}
\end{align}
また、加減算の計算量は
\begin{align}
\sum_{k=1}^{n} (n - 1)(n - k + 1)
&= (n-1)\sum_{k=1}^{n}(n - k + 1) \\
&= (n -1)\sum_{l=0}^{n-1} (l + 1) \\
&= (n - 1)\biggl(\sum_{l=0}^{n-1} l + \sum_{l=0}^{n-1} 1\biggr) \\
&= (n - 1)\biggl(\frac{n(n-1)}{2} + n\biggr)\\
&= (n - 1)\frac{n(n+1)}{2} \\
&= \frac{n(n^2-1)}{2}
\end{align}
したがって、全体の計算量は
\begin{align}
\frac{n^2(n+2)}{2} + \frac{n(n^2-1)}{2}
&= \frac{n(n^2+2n+n^2-1)}{2}\\
&= \frac{n(2n^2+2n-1)}{2}\\
&\sim \frac{2n^3}{2} \\
&= n^3
\end{align}
となります。
Gaussの消去法の計算量
前進消去の計算量
まず、前進消去(forward elimination)の部分を計算します。一番内側のループは、
for (j = k + 1; j <= n; j++) {
a[i][j] -= q * a[k][j];
}
です。このループの中では減算が1回、乗算が1回行われてますので、このループ全体の合計は、
\begin{equation}
乗除算:
\sum_{j=k+1}^{n} 1 = (n - k)\\
加減算:
\sum_{j=k+1}^{n} 1 = (n - k)
\end{equation}
となります。その外側のループを見ると
for (i = k + 1; i <= n; i++) {
q = a[i][k] / a[k][k];
for (j = k + 1; j <= n; j++) {
a[i][j] -= q * a[k][j];
}
c[i] -= q * c[k];
}
一番内側のループの他に乗除算が2回、加減算が1回行われています。したがって、
\begin{equation}
乗除算:
\sum_{i=k+1}^{n} (n - k + 2) = (n - k)(n - k + 2)\\
加減算:
\sum_{i=k+1}^{n} (n - k + 1) = (n - k)(n - k + 1)
\end{equation}
となります。最後に一番外側のループを見ると2番目のループが入っているだけなので、
\begin{equation}
乗除算:
\sum_{k=1}^{n-1} (n - k)(n - k + 2) \\
加減算:
\sum_{k=1}^{n-1} (n - k)(n - k + 1)
\end{equation}
となります。この式はΣの内側にkがあるので先ほどのように簡単には計算できまん。まず、乗除算分を計算してみます。l = (n - k)とおくと
\begin{align}
\sum_{k=1}^{n-1} (n - k)(n - k + 2)
&= \sum_{l=1}^{n-1} l(l+2)\\
&= \sum_{l=1}^{n-1} (l^2+2l)\\
&= \sum_{l=1}^{n-1} l^2 + 2\sum_{l=1}^{n-1} l \\
&= \frac{n(n-1)(2n-1)}{6} + n(n-1) \\
&= \frac{n(n-1)(2n+5)}{6}
\end{align}
加減算分は同様に計算すると、
\begin{align}
\sum_{k=1}^{n-1} (n - k)(n - k + 1)
&= \sum_{l=1}^{n-1} l(l+1)\\
&= \sum_{l=1}^{n-1} (l^2+l)\\
&= \sum_{l=1}^{n-1} l^2 + \sum_{l=1}^{n-1} l \\
&= \frac{n(n-1)(2n-1)}{6} + \frac{n(n-1)}{2} \\
&= \frac{n(n-1)(2n+2)}{6}\\
&= \frac{n(n-1)(n+1)}{3}
\end{align}
となります。
後退代入の計算量
次に後退代入の計算量を見積もります。後退代入のコードは、
c[n] /= a[n][n];
for (k = n - 1; k >= 0; k--) {
s = c[k];
for (j = k + 1; j <= n; j++) {
s -= a[k][j] * c[j];
}
c[k] = s / a[k][k];
}
です。内側のループの演算回数は、
\begin{equation}
乗除算:
\sum_{j=k+1}^{n} 1 = (n - k)\\
加減算:
\sum_{j=k+1}^{n} 1 = (n - k)
\end{equation}
となります。プログラム全体で、乗除算の回数は
\begin{align}
\biggl(\sum_{k=1}^{n-1} ((n - k) + 1)\biggr)+1
&= \sum_{k=1}^{n-1}n - \sum_{k=1}^{n-1} k + \sum_{k=1}^{n-1} 1 + 1 \\
&= n(n-1) - \frac{n(n-1)}{2} + (n-1) + 1\\
&= \frac{n(n-1)}{2}+n\\
&= \frac{n(n+1)}{2}
\end{align}
となります。加減算の回数は
\begin{align}
\sum_{k=1}^{n-1} (n - k)
&= n(n-1) - \frac{n(n-1)}{2}\\
&= \frac{n(n-1)}{2}
\end{align}
です。
全体の計算量
乗除算の全体が
\begin{align}
\frac{n(n-1)(2n+5)}{6} + \frac{n(n+1)}{2}
&= \frac{n((n-1)(2n+5)+3(n+1))}{6} \\
&= \frac{n(2n^2 + 3n - 5 + 3n + 3)}{6} \\
&= \frac{n(2n^2 + 6n - 2)}{6} \\
&= \frac{n(n^2 + 3n - 1)}{3}
\end{align}
となり、加減算の全体が
\begin{align}
\frac{n(n-1)(n+1)}{3} + \frac{n(n-1)}{2}
&= \frac{n(2(n-1)(n+1)+3(n-1))}{6} \\
&= \frac{n(2n^2 - 2 + 3n - 3)}{6} \\
&= \frac{n(2n^2 + 3n - 5)}{6} \\
&\end{align}
となります。したがって、全体の計算量は
\begin{align}
\frac{n(n^2 + 3n - 1)}{3} + \frac{n(2n^2 + 3n - 5)}{6}
&= \frac{n(2n^2+6n-2 + 2n^2+3n-5)}{6} \\
&= \frac{n(4n^2+9n-5)}{6} \\
&\approx \frac{4}{6}n^3 \\
&= \frac{2}{3}n^3
\end{align}
Gauss Jordanの計算量が
n^3
でしたので、Guassの消去法の方が計算量が少ないことになります。
まとめ
Gaussの消去法とGauss Jordanでは操作がよく似ていましたが、計算量には明らかに差がありました。Gaussの消去法の計算量は、計算が簡単な割にはなかなか優秀です。今回は解が簡単に求まる例を使いましたので、サンプルのプログラムも簡単なものになっています。連立1次方程式は必ず解が求まるわけではなく、解が不定な場合や解がない場合があります。また、操作の途中で対角要素が0になると行の交換を行わなければ操作を続けられません。実用的なブログラムにするためには、それらを考慮する必要があります。また、除算を多用するため誤差の考慮も必要です、
また、高校数学の範囲でなければ、LU分解など他の解法も幾つかあります。連立一次方程式は中学ぐらいから習いますが、結構奥が深いものです。