プログラマの為の
数学勉強会
第8回

謝辞

この会の企画・会場設備の提供をして頂きました
㈱ ワークスアプリケーションズ様
にこの場をお借りして御礼申し上げます。

この資料について

• http://nineties.github.com/math-seminar に置いてあります。
• SVGに対応したブラウザで見て下さい。主要なブラウザで古いバージョンでなければ大丈夫だと思います。
• 内容の誤り、プログラムのバグは@9_tiesかkoichi.nakamur AT gmail.comまでご連絡下さい。
• サンプルプログラムはPythonで記述しています。

行列式・逆行列

　前回,3つの行基本変形

• \(P_{i,j}\): 第\(i\)行と第\(j\)行を入れ替える
• \(Q_{i,c}\): 第\(i\)行を\(c \neq 0\)倍する
• \(R_{i,j,c}\): 第\(i\)行を\(c\)倍したものを第\(j\)行に加える

を組み合わせて連立一次方程式を解くことが出来るという話をしました。

行列式とは,このような行列の変形に関して良い性質を持った量です。 連立一次方程式の可解性や,(来週やる)線型写像の性質など様々な分析に利用出来る便利な道具であると言えます。

交代性より,同じ列があったら行列式の値は\(0\)となります。 \[|\cdots, \mathbf{a}, \cdots, \mathbf{a}, \cdots| = -|\cdots, \mathbf{a}, \cdots, \mathbf{a}, \cdots|\ \Leftrightarrow\ \color{yellow}{|\cdots, \mathbf{a}, \cdots, \mathbf{a}, \cdots| = 0} \]

この行列式の3つの性質から,行列式の値を一意に計算する事が出来ます。

\[ \begin{vmatrix} 1 & 2 \\ 3 & 4 \end{vmatrix} \] を例に考えて見ましょう。

まず第1列を \[ \begin{pmatrix} 1 \\ 3 \end{pmatrix} = \color{yellow}{\begin{pmatrix} 1 \\ 0 \end{pmatrix} + 3\begin{pmatrix} 0 \\ 1 \end{pmatrix}} \] と分解して線型性を用いれば \[ \begin{vmatrix} 1 & 2 \\ 3 & 4 \end{vmatrix} = \begin{vmatrix} \color{yellow}{1} & 2 \\ \color{yellow}{0} & 4 \end{vmatrix} + \color{yellow}{3}\begin{vmatrix} \color{yellow}{0} & 2 \\ \color{yellow}{1} & 4 \end{vmatrix} \] となります。(これは第一列に関する余因子展開と呼ばれます。)

続いて第2列を \[ \begin{pmatrix} 2 \\ 4 \end{pmatrix} = \color{yellow}{2\begin{pmatrix} 1 \\ 0 \end{pmatrix} + 4\begin{pmatrix} 0 \\ 1 \end{pmatrix}} \] と分解して線型性を用いれば \[\begin{aligned} \begin{vmatrix} 1 & 2 \\ 3 & 4 \end{vmatrix} &= \begin{vmatrix} 1 & 2 \\ 0 & 4 \end{vmatrix} + 3\begin{vmatrix} 0 & 2 \\ 1 & 4 \end{vmatrix} \\ &= \left\{\color{yellow}{2}\begin{vmatrix} 1 & \color{yellow}{1} \\ 0 & \color{yellow}{0} \end{vmatrix}+ \color{yellow}{4}\begin{vmatrix} 1 & \color{yellow}{0} \\ 0 & \color{yellow}{1} \end{vmatrix}\right\} + 3\left\{\color{yellow}{2}\begin{vmatrix} 0 & \color{yellow}{1} \\ 1 & \color{yellow}{0} \end{vmatrix}+ \color{yellow}{4}\begin{vmatrix} 0 & \color{yellow}{0} \\ 1 & \color{yellow}{1} \end{vmatrix}\right\} \end{aligned} \] となりますが,交代性より列が重複すると0なので \[ \begin{vmatrix} 1 & 2 \\ 3 & 4 \end{vmatrix} = 4\begin{vmatrix} 1 & 0 \\ 0 & 1 \end{vmatrix} + 6\begin{vmatrix} 0 & 1 \\ 1 & 0 \end{vmatrix} \] となります。

次に列の並べ替えをして\(E\)に形に揃えます。交代性より,列を交換する度に符号が変化します。 \[ \begin{vmatrix} 1 & 2 \\ 3 & 4 \end{vmatrix} = 4\begin{vmatrix} 1 & 0 \\ 0 & 1 \end{vmatrix} - 6\begin{vmatrix} 1 & 0 \\ 0 & 1 \end{vmatrix} \] 最後に\( \det E = 1\)を利用すれば \[ \begin{vmatrix} 1 & 2 \\ 3 & 4 \end{vmatrix} = 4 - 6 = \color{yellow}{-2} \] が求める結果となります。

このように,行列式を計算する為には行列の列を \[ \begin{aligned} \begin{pmatrix} a_1\\a_2\\ \vdots \\a_n \end{pmatrix} &= a_1\begin{pmatrix} 1 \\ 0 \\ \vdots \\ 0 \end{pmatrix} + a_2\begin{pmatrix} 0 \\ 1 \\ \vdots \\ 0 \end{pmatrix} + \cdots + a_n\begin{pmatrix} 0 \\ 0 \\ \vdots \\ 1 \end{pmatrix} \\ &= \color{yellow}{a_1\mathbf{e_1}+ a_2\mathbf{e_2}+\cdots+ a_n\mathbf{e_n}} \end{aligned} \] と分解し線型性を用いて展開するという事を繰り返していきます。但し\(\mathbf{e}_i\)は第\(i\)成分が\(1\),それ以外が\(0\)のベクトルです。

展開結果には\(\mathbf{e}_1,\cdots,\mathbf{e}_n\)を１つずつ列に持つ項だけが残り,\(E\)の形にする為に必要な列の交換の回数だけ\(-1\)倍されるという事になります。

【答え】
\[ \begin{aligned} \begin{vmatrix} a & b \\ c & d \end{vmatrix} &= | a\mathbf{e}_1 + c\mathbf{e}_2, b\mathbf{e}_1 + d\mathbf{e}_2 |\\ &= ad|\mathbf{e}_1,\mathbf{e}_2|+bc|\mathbf{e}_2,\mathbf{e}_1| \\ &= ad|\mathbf{e}_1,\mathbf{e}_2|-bc|\mathbf{e}_1,\mathbf{e}_2| \\ &= \color{yellow}{ad-bc} \end{aligned} \]

【答え】
\[ \begin{aligned} &\begin{vmatrix} a & b & c \\ d & e & f \\ g & h & i\end{vmatrix} = |a\mathbf{e}_1+d\mathbf{e}_2+g\mathbf{e}_3,b\mathbf{e}_1+e\mathbf{e}_2+h\mathbf{e}_3,c\mathbf{e}_1+f\mathbf{e}_2+i\mathbf{e}_3| \\ =& aei|\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3| + ahf|\mathbf{e}_1,\mathbf{e}_3,\mathbf{e}_2| + dbi|\mathbf{e}_2,\mathbf{e}_1,\mathbf{e}_3| +\\ &dhc|\mathbf{e}_2,\mathbf{e}_3,\mathbf{e}_1| + gbf|\mathbf{e}_3,\mathbf{e}_1,\mathbf{e}_2| + gec|\mathbf{e}_3,\mathbf{e}_2,\mathbf{e}_1| \\ =& aei|\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3| -ahf|\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3| -dbi|\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3| +\\ &dhc|\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3| + gbf|\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3| -gec|\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3| \\ =&aei-afh+bfg-bdi+cdh-ceg \end{aligned} \]

以上で見たことから以下の公式が成立します。

ここで,\(\mathrm{sign}(i_1, \cdots, i_n)\)は,\((i_1,\cdots,i_n)\)を\((1,\cdots,n)\)に並べ替える為に必要な二数の入れ替え(互換)の回数が偶数なら1,奇数なら-1です。

例えば\( \mathrm{sign}(2,3,1) = 1,\ \mathrm{sign}(3,2,1) = -1 \)といった具合です。 \[ \begin{aligned} & (2, 3, 1) \ \xrightarrow{\text{1,2入れ替え}}\ (1, 3, 2)\ \xrightarrow{\text{2,3入れ替え}}\ (1, 2, 3) \qquad \text{(偶数回)}\\ & (3, 2, 1) \ \xrightarrow{\text{1,3入れ替え}}\ (1,2,3) \qquad\text{(奇数回)} \end{aligned} \]

\(\{1,2\}\)の置換は\((1,2)\)と\((2,1)\)ですから \[ \begin{aligned} \det A &= \mathrm{sign}(1,2)a_{11}a_{22} + \mathrm{sign}(2,1)a_{21}a_{12} \\ &= a_{11}a_{22}-a_{12}a_{21} \\ &= ad-bc \end{aligned} \] となります。

【答え】
\[ \begin{aligned} & \det A\\ = & \mathrm{sign}(1,2,3)a_{11}a_{22}a_{33} + \mathrm{sign}(1,3,2)a_{11}a_{32}a_{23} + \mathrm{sign}(2,1,3)a_{21}a_{12}a_{33} + \\ & \qquad \mathrm{sign}(2,3,1)a_{21}a_{32}a_{13} + \mathrm{sign}(3,1,2)a_{31}a_{12}a_{23} + \mathrm{sign}(3,2,1)a_{31}a_{22}a_{13} \\ = & a_{11}a_{22}a_{33} - a_{11}a_{32}a_{23} - a_{21}a_{12}a_{33} + a_{21}a_{32}a_{13} + a_{31}a_{12}a_{23} - a_{31}a_{22}a_{13} \\ \end{aligned} \]

公式 \[ \det A = \sum \mathrm{sign}(i_1,\cdots,i_n) a_{i_1,1}a_{i_2,2}\cdots a_{i_n,n} \] より行と列を入れ替えた \[ \det A = \sum \mathrm{sign}(i_1,\cdots,i_n) a_{1,i_1}a_{2,i_2}\cdots a_{n,i_n} \] が成立する事も容易に示せます。\((i_1,\cdots,i_n)\)を並べ替えて\((1,\cdots,n)\)にする互換の回数と, \((1,\cdots,n)\)を並べ替えて\((i_1,\cdots,i_n)\)にする互換の回数は当然等しいからです。

従って

が成り立ちます。また,今まで列に対して行っていた全ての変形は行についても同様に行えるという事になります。

行基本変形と行列式

ということで,3つの行基本変形と行列式の関係を述べる事が出来ます。

• \(P_{i,j}\)を行うと\(-1\)倍 \[|\cdots, \mathbf{a}_i, \cdots, \mathbf{a}_j, \cdots| = - |\cdots, \mathbf{a}_j, \cdots, \mathbf{a}_i, \cdots|\]
• \(Q_{i,c}\)を行うと\(c\)倍 \[|\cdots, c\mathbf{a}_i, \cdots| = c |\cdots, \mathbf{a}_i, \cdots| \]
• \(R_{i,j,c}\)に関して不変 \[ \begin{aligned} &|\cdots, \mathbf{a}_i,\cdots,c\mathbf{a}_i+\mathbf{a}_j,\cdots| \\ = &c|\cdots, \mathbf{a}_i,\cdots,\mathbf{a}_i,\cdots| + |\cdots, \mathbf{a}_i,\cdots,\mathbf{a}_j,\cdots| \\ = &|\cdots, \mathbf{a}_i,\cdots,\mathbf{a}_j,\cdots| \end{aligned} \]

例えば,

\[ \begin{aligned} \left[\begin{array}{cc} 3& 1 \\ 2& 4 \\ \end{array}\right] &\xrightarrow{P_{1,2}} \left[\begin{array}{cc} 2& 4 \\ 3& 1 \\ \end{array}\right] \xrightarrow{Q_{1,\frac{1}{2}}} \left[\begin{array}{cc} 1& 2 \\ 3& 1 \\ \end{array}\right] \xrightarrow{R_{1,2,-3}} \left[\begin{array}{cc} 1& 2 \\ 0&-5 \\ \end{array}\right] \\ &\xrightarrow{Q_{2,-\frac{1}{5}}} \left[\begin{array}{cc} 1& 2 \\ 0& 1 \\ \end{array}\right] \xrightarrow{R_{2,1,-2}} \left[\begin{array}{cc} 1& 0 \\ 0& 1 \\ \end{array}\right] \end{aligned} \] という行基本変形の列に関して \[ \begin{vmatrix}3 & 1 \\ 2 & 4\end{vmatrix} \times (-1) \times \frac12\times \left(-\frac15\right) = \begin{vmatrix} 1 & 0 \\ 0 & 1\end{vmatrix}\] つまり \[ \begin{vmatrix}3 & 1 \\ 2 & 4\end{vmatrix} = 10 \] が成り立っています。

続いて,正方行列\(A,B\)の積の行列式 \[ \det (AB) \] がどうなるか説明します。

二次の場合を考えてみましょう。 \[ \det \left(\begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} e & f \\ g & h \end{pmatrix}\right) \] ここで, \[ \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} e & f \\ g & h \end{pmatrix} = \left( e\begin{pmatrix} a \\ c\end{pmatrix} + g\begin{pmatrix} b \\ d\end{pmatrix}, f\begin{pmatrix} a \\ c\end{pmatrix} + h\begin{pmatrix} b \\ d\end{pmatrix} \right) \] という分解が出来るので,

線型性・交代性を使えば \[ \begin{aligned} \left|\begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} e & f \\ g & h \end{pmatrix}\right| &= \left| e\begin{pmatrix} a \\ c\end{pmatrix} + g\begin{pmatrix} b \\ d\end{pmatrix}, f\begin{pmatrix} a \\ c\end{pmatrix} + h\begin{pmatrix} b \\ d\end{pmatrix} \right| \\ &= eh\begin{vmatrix}a&b\\c&d\end{vmatrix} +fg\begin{vmatrix}b&a\\d&c\end{vmatrix} \\ &= eh\begin{vmatrix}a&b\\c&d\end{vmatrix} -fg\begin{vmatrix}a&b\\c&d\end{vmatrix} \\ &= \begin{vmatrix}a&b\\c&d\end{vmatrix}\times(eh-fg) \\ &= (ad-bc)(eh-fg) \\ \end{aligned} \] となり, \[ \det(AB) = (\det A)(\det B) \] が成立しています。

一般の場合も同様で,以下の性質が成立します。

余因子展開

行列式の計算に利用した \[ \begin{vmatrix} 1 & \color{yellow}{2} & 3 \\ 4 & \color{yellow}{5} & 6 \\ 7 & \color{yellow}{8} & 9\end{vmatrix} = \color{yellow}{2}\begin{vmatrix} 1 & \color{yellow}{1} & 3 \\ 4 & \color{yellow}{0} & 6 \\ 7 & \color{yellow}{0} & 9\end{vmatrix} + \color{yellow}{5}\begin{vmatrix} 1 & \color{yellow}{0} & 3 \\ 4 & \color{yellow}{1} & 6 \\ 7 & \color{yellow}{0} & 9\end{vmatrix} + \color{yellow}{8}\begin{vmatrix} 1 & \color{yellow}{0} & 3 \\ 4 & \color{yellow}{0} & 6 \\ 7 & \color{yellow}{1} & 9\end{vmatrix} \] といった展開を余因子展開と言います。上は第2列に関する余因子展開です。

ここで\(\mathbf{e}_i\)で行列\(A\)の第\(j\)列を置き換えた物の行列式を第\((i,j)\)余因子と言い,\(A_{ij}\)と書くことにします。すると,上の式は \[ |A| = a_{12}A_{12} + a_{22}A_{22} + a_{32}A_{32} \] と書くことが出来ます。

ところで,余因子に関して \[ \begin{aligned} A_{23} &= \begin{vmatrix} 1 & 2 & \color{yellow}{0} & 4\\ \color{yellow}{5} & \color{yellow}{6} & \color{yellow}{1} & \color{yellow}{8} \\ 9 & 10 & \color{yellow}{0} & 12\\ 13 & 14 & \color{yellow}{0} & 16 \end{vmatrix} = (-1)^{2+3}\begin{vmatrix} 1 & 2 & 4 &\color{yellow}{0}\\ 9 & 10 & 12 &\color{yellow}{0}\\ 13 & 14 & 16 &\color{yellow}{0}\\ \color{yellow}{5} & \color{yellow}{6} & \color{yellow}{8} &\color{yellow}{1} \end{vmatrix} \\ &= (-1)^{2+3}\begin{vmatrix} 1 & 2 & 4 \\ 9 & 10 & 12 \\ 13 & 14 & 16 \\ \end{vmatrix} \end{aligned} \] といった行列式のサイズの縮小が可能です。

\((-1)^{2+3}\)は第2行,第3列を第4行,第4列に移す際の符号の変化です。移した後は,第\((4,4)\)成分の1を含む項しか残らず,さらに\(\mathrm{sign}(i,j,k,4)=\mathrm{sign}(i,j,k)\)である事から最後の等式が成立します。

転置行列を考えれば,第\(i\)行に関する余因子展開\( \det A = \sum_i a_{ij} A_{ij} \)も成立します。

第3列に関して展開すれば \[ \begin{aligned} \det A &= a_{13}A_{13} + a_{23}A_{23} + a_{33}A_{33} \\ &= A_{13} + 4A_{33} \quad(\because a_{23} = 0) \\ &= (-1)^{1+3}\begin{vmatrix} 3 & 2 \\ 3 & 1 \end{vmatrix} + (-1)^{3+3}\times 4 \begin{vmatrix} 2 & 4 \\ 3 & 2 \end{vmatrix} \\ &= -3 + 4\times (-8) \\ &= -35 \end{aligned} \] となります。別の列・行での展開もやってみて下さい。

【答え】
\[ \begin{aligned} & \begin{vmatrix} 2 & 0 & 1 & -1 \\ 0 & 2 & 3 & 0 \\ 2 & 0 & 1 & 2 \\ -2 & 1 & 0 & 3 \end{vmatrix} = 2\begin{vmatrix} 2 & 1 & -1 \\ 2 & 1 & 2 \\ -2 & 0 & 3 \end{vmatrix} + \begin{vmatrix} 2 & 1 & -1 \\ 0 & 3 & 0 \\ 2 & 1 & 2 \end{vmatrix} \\ = & 2\left\{ -\begin{vmatrix}2 & 2 \\ -2 & 3 \end{vmatrix} + \begin{vmatrix}2 & -1 \\ -2 & 3 \end{vmatrix} \right\} + 3\begin{vmatrix}2 & -1 \\ 2 & 2 \end{vmatrix} \\ = & 2\times(-10+4)+3\times 6\\ = &6 \end{aligned} \]

逆行列

【一意性の証明】
\( AX=XA=E \)と\( AY=YA=E \)が成立するならば \[ AX=E\Rightarrow YAX=Y\Rightarrow EX=Y\Rightarrow X=Y\] □

行列式の理論を用いて逆行列の公式を導きましょう。

例えば, \[ A = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix} \] を第1列について余因子展開すると \[ \begin{vmatrix} \color{yellow}{1} & 2 & 3 \\ \color{yellow}{4} & 5 & 6 \\ \color{yellow}{7} & 8 & 9\end{vmatrix} = \color{yellow}{1}\begin{vmatrix} \color{yellow}{1} & 2 & 3 \\ \color{yellow}{0} & 5 & 6 \\ \color{yellow}{0} & 8 & 9\end{vmatrix} + \color{yellow}{4}\begin{vmatrix} \color{yellow}{0} & 2 & 3 \\ \color{yellow}{1} & 5 & 6 \\ \color{yellow}{0} & 8 & 9\end{vmatrix} + \color{yellow}{7}\begin{vmatrix} \color{yellow}{0} & 2 & 3 \\ \color{yellow}{0} & 5 & 6 \\ \color{yellow}{1} & 8 & 9\end{vmatrix} \] つまり,\( a_{11}A_{11} + a_{21}A_{21} + a_{31}A_{31} = |A|\)が成り立ちます。第2列,第3列に関して展開しても同様なので \[ \begin{aligned} &a_{11}A_{11} + a_{21}A_{21} + a_{31}A_{31} = |A|\\ &a_{12}A_{12} + a_{22}A_{22} + a_{32}A_{32} = |A|\\ &a_{13}A_{13} + a_{23}A_{23} + a_{33}A_{33} = |A|\\ \end{aligned} \] が成立します。

では,列の番号が揃っていない \[ a_{1\color{yellow}{2}}A_{1\color{red}{1}} + a_{2\color{yellow}{2}}A_{2\color{red}{1}} + a_{3\color{yellow}{2}}A_{3\color{red}{1}} \] などはどうなるでしょうか？

これは \[ \color{yellow}{2}\begin{vmatrix} \color{yellow}{1} & 2 & 3 \\ \color{yellow}{0} & 5 & 6 \\ \color{yellow}{0} & 8 & 9\end{vmatrix} + \color{yellow}{5}\begin{vmatrix} \color{yellow}{0} & 2 & 3 \\ \color{yellow}{1} & 5 & 6 \\ \color{yellow}{0} & 8 & 9\end{vmatrix} + \color{yellow}{8}\begin{vmatrix} \color{yellow}{0} & 2 & 3 \\ \color{yellow}{0} & 5 & 6 \\ \color{yellow}{1} & 8 & 9\end{vmatrix} = \begin{vmatrix} \color{yellow}{2} & 2 & 3 \\ \color{yellow}{5} & 5 & 6 \\ \color{yellow}{8} & 8 & 9\end{vmatrix} = 0 \] となり,第2列が重複してしまうので\(0\)となります。

この様に \[ a_{1i}A_{1j} + a_{2i}A_{2j} + a_{3i}A_{3j} \] は\(i = j\)(列番号が同じ)なら\(|A|\),\(i \neq j\)(列番号が違う)なら\(0\)となります。

すると \[ \begin{pmatrix} A_{11} & A_{21} & A_{31} \\ A_{12} & A_{22} & A_{32} \\ A_{13} & A_{23} & A_{33} \\ \end{pmatrix} \begin{pmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \\ \end{pmatrix} = \begin{pmatrix} |A| & 0 & 0 \\ 0 & |A| & 0 \\ 0 & 0 & |A| \\ \end{pmatrix} \] となります。つまり,\(|A| \neq 0\)であるならば \[ \frac{1}{|A|} \begin{pmatrix} A_{11} & A_{21} & A_{31} \\ A_{12} & A_{22} & A_{32} \\ A_{13} & A_{23} & A_{33} \\ \end{pmatrix} A = E \] となっています。行に関する余因子展開を考えれば,全く同様にして \[ A \cdot \frac{1}{|A|} \begin{pmatrix} A_{11} & A_{21} & A_{31} \\ A_{12} & A_{22} & A_{32} \\ A_{13} & A_{23} & A_{33} \\ \end{pmatrix} = E \] も示されるので,逆行列が得られた事になります。

十分性は前ページの通りです。また\(A\)が正則ならば \[ AA^{-1} = E\ \Rightarrow\ (\det A)(\det A^{-1}) = \det E \neq 0 \] なので\(\det A \neq 0\)である事が必要です。

まず,\(\det A = ad-bc\)です。そして \[ \begin{aligned} A_{11} &= (-1)^{1+1}|d| = d\\ A_{12} &= (-1)^{1+2}|c| = -c\\ A_{21} &= (-1)^{2+1}|b| = -b\\ A_{22} &= (-1)^{2+2}|a| = a \end{aligned} \] ですから\(ad-bc\neq 0\)の時 \[ A^{-1} = \frac{1}{ad-bc}\begin{pmatrix} A_{11} & A_{21} \\ A_{12} & A_{22} \end{pmatrix} = \color{yellow}{\frac{1}{ad-bc}\begin{pmatrix} d & -b \\ -c & a \end{pmatrix}} \]

まず\(\det A = -3\)で, \[ \begin{aligned} A_{11} &= 1\cdot1 - 3\cdot 0 = 1,\quad A_{12} = -(2\cdot 1-3\cdot 0) = -2,\quad A_{13} = 2\cdot 0 -1\cdot 0 = 0 \\ A_{21} &= -(2\cdot 1-0\cdot 0) = -2,\quad A_{22}=1\cdot 1-0\cdot 0=1,\quad A_{23} = 1\cdot 0-2\cdot 0 = 0 \\ A_{31} &= 2\cdot 3-0\cdot 1=6,\quad A_{32} = -(1\cdot 3-0\cdot 2)=-3,\quad A_{33} = 1\cdot 1-2\cdot 2 = -3 \end{aligned} \] となりますから \[ A^{-1} = -\frac{1}{3}\begin{pmatrix} 1 & -2 & 6 \\-2 & 1 & -3 \\0 & 0 & -3 \end{pmatrix} =\color{yellow}{\begin{pmatrix} -\frac13 & \frac23 & -2 \\ \frac23 & -\frac13 & 1 \\ 0 & 0 & 1 \end{pmatrix}} \]

【証明】
\(AX=E\)ならば\((\det A)(\det X)=1\neq 0\)より\(\det A\neq 0\)。よって\(A^{-1}\)が存在するので \[ AX=E\Rightarrow A^{-1}AXA=A^{-1}A \Rightarrow XA=E \] 逆も同様。
2つ目の命題は自明。3つ目は \[ ABB^{-1}A^{-1} = AA^{-1} = E \] であることより成立する。

連立一次方程式の可解性

連立一次方程式はいつでも唯一の解を持つわけではありません。例えば \[ \left\{\begin{array}{c} x + y = 1 \\ 2x + 2y = 2 \end{array}\right. \ \Leftrightarrow\ x + y = 1 \] は無数の解を持ちますし, \[ \left\{\begin{array}{c} x + y = 1 \\ 2x + 2y = 3 \end{array}\right. \ \Leftrightarrow\ x + y = 1 = \frac{3}{2} \] は解を持ちません。

この様に,連立一次方程式は

• 唯一の解を持つ
• 無数の解を持つ
• 解を持たない

の3つの場合があり,特に唯一解を持つ条件の判定は重要です。

まず,連立一次方程式 \[ A\mathbf{x} = \mathbf{b} \] について\(A\)が正則ならば,両辺に\(A^{-1}\)を左から掛ければ \[ \color{yellow}{\mathbf{x} = A^{-1}\mathbf{b}} \] となります。つまり, \[ \color{yellow}{\text{$A$が正則}\ \Rightarrow\ \text{$A\mathbf{x}=\mathbf{b}$が唯一解を持つ}} \] が言えます。

次に\(A\mathbf{y}=0\)を満たす\(\mathbf{y}\neq 0\)が存在するとしましょう。すると \[ A\mathbf{x} = \mathbf{b},\ A\mathbf{y} = 0\ \Rightarrow A(\mathbf{x}+\mathbf{y}) = \mathbf{b}\] であり,\(\mathbf{x} \neq \mathbf{x}+\mathbf{y}\)ですから \[ A\mathbf{y} = 0\text{が$0$でない解を持つ}\ \Rightarrow\ \text{$A\mathbf{x}=\mathbf{b}$が唯一解を持たない} \] が成立します。さらに,次ページで示すように \[ \text{$A$が正則でない}\ \Rightarrow\ A\mathbf{y}=0\text{が$0$でない解を持つ}\qquad\cdots(1) \] なので,対偶を考えて \[ \color{yellow}{\text{$A\mathbf{x}=\mathbf{b}$が唯一解を持つ}\ \Rightarrow\ \text{$A$が正則}} \] も成立します。

数学的帰納法で証明します。

【証明】
\((1,1)\)型行列については明らかに成立する。
そこで\(n-1\)元連立一次方程式について(1)が成り立つと仮定し,\(n\)元連立一次方程式 \[ \left[\begin{array}{cccccc} a_{11}&a_{12}&\cdots&a_{1n}&:& b_1 \\ a_{21}&a_{22}&\cdots&a_{2n}&:& b_2 \\ \vdots&\vdots&\ddots&\vdots&:&\vdots\\ a_{n1}&a_{n2}&\cdots&a_{nn}&:& b_n \\ \end{array}\right] \] について考える。ここで第1列が全て\(0\)なら\(x_1\)を自由に取れるので(1)が成立し,そうでないならばピボット行の交換によって\(a_{11}\neq 0\)と出来,前進消去によって \[ \left[\begin{array}{cccccc} a_{11}&a_{12}&\cdots&a_{1n}&:& b_1 \\ 0 &a'_{22}&\cdots&a'_{2n}&:& b'_2 \\ \vdots&\vdots&\ddots&\vdots&:&\vdots\\ 0 &a'_{n2}&\cdots&a'_{nn}&:& b'_n \\ \end{array}\right] \] とできる。この際,係数行列の行列式は符号を除いて変化しない。従って\(\det A = 0\)ならば \[ \det\begin{pmatrix} a'_{22}&\cdots&a'_{2n}\\ \vdots&\ddots&\vdots \\ a'_{n2}&\cdots&a'_{nn}\\ \end{pmatrix} = 0 \] となるが,この時帰納法の仮定より\((x_2,\cdots,x_n)^T\neq 0\)なる解が存在する。 □

ここで\(A^{-1}\)の公式を代入すればクラメルの公式という公式が得られますが,あまり使わないので省略します。

行列式の計算法

計算機で行列式を計算する為には,前進消去を用いる事が出来ます。

前進消去で得られた上三角行列では対角成分の積が行列式の値となります。(第1列からは\(a_{11}\)を取る他ない,すると第2列からは\(a_{22}\)を取る他ない\(\cdots\)) \[ \det\begin{pmatrix} a_{11}&a_{12}&\cdots&a_{1,n-1}&a_{1n} \\ 0 &a_{22}&\cdots&a_{2,n-1}&a_{2n} \\ 0 &0 &\ddots&\ddots&a_{3n} \\ \vdots&\vdots&\ddots&\ddots&\vdots \\ 0 &0 &\cdots&0&a_{nn} \\ \end{pmatrix} = a_{11}a_{22}\cdots a_{nn} \]

前進消去では行基本変形\(P_{i,j},R_{i,j,c}\)を使いますので,\(P_{i,j}\)による符号の変化だけ気をつければ良いです。また,行列式が\(0\)となる条件はピボットのどれかが\(0\)となる事です。

細かい話になりますが,\(n\)が大きいと\(a_{11}a_{22}\cdots a_{nn}\)がオーバーフロー・アンダーフローを起こします。そういう時は 行列式の絶対値の対数\( \ln|\det A|=\sum_i\ln|a_{ii}|\)と符号に分けて計算するルーチンを使います。

コーディング例です。前回やったガウス消去法のコードをちょっといじるだけです。 また,計算量は前進消去の計算量で決まるので\(\mathcal{O}(n^3)\)となります。

# -*- coding: utf-8 -*-
import numpy as np

def determinant(A):
    # 前進消去
    A = np.copy(A)
    n = A.shape[0]
    p = np.arange(n)    # [0,1,2,...,n-1] 
    det = 1.0
    for k in xrange(n-1):

        # ピボット選択
        pivot_idx = p[k]
        pivot_max = 0
        for i in xrange(k, n):
            v = abs(A[p[i], k])
            if v > pivot_max:
                pivot_max = v
                pivot_idx = i

        # 実際には誤差があるので pivot_max == 0.0よりは
        # pivot_max < (小さい値) とした方が良いかも。
        if pivot_max == 0.0:
            return 0.0

        # ピボット行の交換
        if p[k] != pivot_idx:
            p[k], p[pivot_idx] = p[pivot_idx], p[k]
            det *= -1 # ピボット交換では符号を変える

        pivot = A[p[k],k]
        det *= pivot # 対角成分を掛け合わせる
        for i in xrange(k+1, n):
            l = A[p[i], k]/pivot

            for j in xrange(k+1, n):
                A[p[i], j] -= l * A[p[k], j]
    det *= A[p[n-1], n-1]   # これを忘れずに
    return det

A = np.array([[3,2,5,1],[-2,3,4,1],[1,2,5,6],[0,7,1,3]],dtype=float)
print determinant(A)

逆行列の計算法

逆行列だけを単独で求めたい場面は少ないと思います。線型代数の殆どの問題は連立一次方程式 \[ A\mathbf{x} = \mathbf{b} \] を解くことであり,これは直接解くアルゴリズム(ガウス消去他)があります。わざわざ\(A^{-1}\)を求めてから\(\mathbf{x}=A^{-1}\mathbf{b}\)を計算する理由はありません。

しかし,もし逆行列が必要であれば \[ AX = E \] という連立一次方程式の解が\(X = A^{-1}\)となります。従って新しいアルゴリズムは不要です。

>>> import numpy as np
>>> from numpy import linalg as LA
>>> A = np.array([[1,5,2],[3,1,2],[4,3,2]], dtype=float)
>>> E = np.identity(3)
>>> A
array([[ 1.,  5.,  2.],
       [ 3.,  1.,  2.],
       [ 4.,  3.,  2.]])
>>> E
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
>>> X = LA.solve(A, E)    # AX=Eを解けば逆行列が求まる
>>> X
array([[-0.25  , -0.25  ,  0.5   ],
       [ 0.125 , -0.375 ,  0.25  ],
       [ 0.3125,  1.0625, -0.875 ]])
>>> np.dot(A, X)          # AX=Eの確認
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
>>> np.dot(X, A)          # XA=Eの確認
array([[ 1.,  0.,  0.],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
>>> LA.inv(A)             # もちろん,専用の関数も用意されている
array([[-0.25  , -0.25  ,  0.5   ],
       [ 0.125 , -0.375 ,  0.25  ],
       [ 0.3125,  1.0625, -0.875 ]])

一歩進んだ話

以上が行列式・逆行列・連立一次方程式の基本です。しかし,実戦的には係数行列の特徴を利用して効率の良い計算をする事が大切になります。

例えば,次のような左下の部分が全て\(0\)になっている係数行列を考えましょう。 ただし\(A,D\)は正方行列であるとします。 \[ \begin{pmatrix} A & B \\ \color{yellow}{O} & D \end{pmatrix} \mathbf{x} = \mathbf{b} \] この特徴を利用して行列式・逆行列の計算,連立一次方程式の求解を行うにはどうすれば良いでしょうか？

行列式を計算するには,前進消去法が使えるのでした。そこで \[ \begin{pmatrix} A & B \\ O & D \end{pmatrix} \rightarrow \begin{pmatrix} U_1 & B' \\ O & U_2 \end{pmatrix} \] と変形しましょう。\(U_1,U_2\)は上三角行列です。

この行列式は\(U_1,U_2\)の対角成分を掛けて,ピボット行の交換回数だけ符号を反転すれば良いのですが,行の交換回数は \[ A \rightarrow U_1 \] の変形に必要な交換回数と \[ D \rightarrow U_2 \] に必要な交換回数の和です。この事から \[ \color{yellow}{\begin{vmatrix} A & B \\ O & D \end{vmatrix} = |A||D|} \] が成立している事が判ります。

\(A^{-1}\)を求めるには\(AX=E\)を解けば良いのですから,先ほどの行列に対しては \[ \begin{aligned} &\begin{pmatrix} A & B \\ O & D \end{pmatrix} \begin{pmatrix} X & Y \\ Z & W \end{pmatrix} = \begin{pmatrix} E & O \\ O & E \end{pmatrix} \\ \Leftrightarrow & AX+BZ=E,\ AY+BW=O,\ DZ=O,\ DW=E \end{aligned} \] を解けば良いです。

求める逆行列が存在するならば前頁の話より\(|A||D|\neq 0\ \Leftrightarrow\ |A|\neq 0, |D|\neq 0\)ですから,\(A^{-1},D^{-1}\)が存在します。従って第3,4式から \( Z = O, W = D^{-1} \)となり,さらにこれを第1,2式に代入すれば \( AX=E,\ AY=-BD^{-1} \)となる事から\( X = A^{-1},\ Y = -A^{-1}BD^{-1} \)となります。

すなわち \[ \color{yellow}{\begin{pmatrix} A & B \\ O & D \end{pmatrix}^{-1} = \begin{pmatrix} A^{-1} & -A^{-1}BD^{-1} \\ O & D^{-1} \end{pmatrix} } \] が成立します。

従って連立方程式 \[ \begin{pmatrix} A & B \\ O & D \end{pmatrix} \begin{pmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \end{pmatrix} = \begin{pmatrix} \mathbf{b}_1 \\ \mathbf{b}_2 \end{pmatrix} \] の解を明示的に書けば \[ \begin{pmatrix} \mathbf{x}_1 \\ \mathbf{x}_2 \end{pmatrix} = \begin{pmatrix} A^{-1} & -A^{-1}BD^{-1} \\ O & D^{-1} \end{pmatrix} \begin{pmatrix} \mathbf{b}_1 \\ \mathbf{b}_2 \end{pmatrix} = \begin{pmatrix} A^{-1}(\mathbf{b}_1 - BD^{-1}\mathbf{b}_2) \\ D^{-1}\mathbf{b}_2 \end{pmatrix} \] です。

また,計算機で解く場合には逆行列の計算は不要で

1. \( D\mathbf{x}_2 = \mathbf{b}_2 \)を解く
2. \( A\mathbf{x}_1 = \mathbf{b}_1-B\mathbf{x}_2 \)を解く

とすればよいです。

こういった話題は,線型代数を利用したソフトウェア開発を始めてから徐々に覚えていけば良いと思います。 本日の内容をよく理解してもらえれば,後に専門的な文献を読む事も出来ると思います。

連立一次方程式の直接解法

前回,最も重要な直接解法であるガウス消去法を紹介しました。今回はより発展的な手法を紹介していきます。

計算量について復習すると\(A\)が\(n\)次正方行列のとき \[ A\mathbf{x} = \mathbf{b} \] をガウス消去法で解くためには,前進消去法で\(\mathcal{O}(n^3)\),後退代入で\(\mathcal{O}(n^2)\)の計算量が必要なのでした。

LU分解

行列\(A\)を下三角行列\(L\)と上三角行列\(U\)の積 \[ A = LU \] と分解する事をLU分解と言います。

すると \[ A\mathbf{x} = \mathbf{b}\ \Leftrightarrow\ LU\mathbf{x}=\mathbf{b}\ \Leftrightarrow\ \left\{\begin{array}{c} L\mathbf{y} = \mathbf{b} \\ U\mathbf{x} = \mathbf{y} \end{array}\right. \] となりますが,前回やった後退代入法と同様に三角行列に対して \[ L\mathbf{y} = \mathbf{b},\ U\mathbf{x} = \mathbf{y}\] はどちらも\(\mathcal{O}(n^2)\)で解くことが可能です。

例えば,連立方程式 \[ \left\{\begin{array}{l} 4x + 5y = 9\\ 2x + y = 3 \end{array}\right. \] を考えましょう。ここで係数行列は \[ \begin{pmatrix} 4 & 5 \\ 2 & 1 \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ \frac12 & 1\end{pmatrix} \begin{pmatrix} 4 & 5 \\ 0 & -\frac{3}{2}\end{pmatrix} \] と分解出来るので,まず下三角行列について \[ \begin{bmatrix} 1 & 0 & : & 9 \\ 1/2 & 1 & : & 3 \end{bmatrix} \rightarrow \begin{bmatrix} 1 & 0 & : & 9 \\ 0 & 1 & : & -3/2 \end{bmatrix} \] と解き,これを使って上三角行列について \[ \begin{bmatrix} 4 & 5 & : & 9 \\ 0 & -3/2 & : & -3/2 \end{bmatrix} \rightarrow \begin{bmatrix} 1 & 0 & : & 1 \\ 0 & 1 & : & 1 \end{bmatrix} \] と解きます。つまり\(x = y = 1\)が求める解となります。

次に,LU分解のアルゴリズムを紹介します。LU分解にもいろいろな種類があるのですが,ドゥーリトル法という\(L\)の対角成分を\(1\)にした \[ \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & & 0 \\ l_{21} & 1 & & 0 \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & 1 \\ \end{pmatrix} \begin{pmatrix} u_{11} & u_{12} & \cdots & u_{1n} \\ 0 & u_{22} & \cdots & u_{2n} \\ & & \ddots & \vdots \\ 0 & 0 & & u_{nn} \\ \end{pmatrix} \] 分解を求める方法を紹介します。(\(U\)の対角成分を\(1\)にする方法はクラウト法と言います。)

アルゴリズムは

\[ \left(\begin{array}{c|ccc} a_{11} & a_{12} & \cdots & a_{1n} \\ \hline a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \\ \end{array}\right) = \left(\begin{array}{c|ccc} 1 & 0 & & 0 \\ \hline l_{21} & 1 & & 0 \\ \vdots & \vdots & \ddots & \\ l_{n1} & l_{n2} & \cdots & 1 \\ \end{array}\right) \left(\begin{array}{c|ccc} u_{11} & u_{12} & \cdots & u_{1n} \\ \hline 0 & u_{22} & \cdots & u_{2n} \\ & & \ddots & \vdots \\ 0 & 0 & & u_{nn} \\ \end{array}\right) \]

と小行列に分解して考えれば,簡単に導けます。つまり \[ \begin{pmatrix} a_{11} & \mathbf{a}_1^T \\ \mathbf{a}_2 & A' \end{pmatrix} = \begin{pmatrix} 1 & \mathbf{0} \\ \mathbf{l} & L' \end{pmatrix} \begin{pmatrix} u_{11} & \mathbf{u}^T \\ \mathbf{0} & U' \end{pmatrix} = \begin{pmatrix} u_{11} & \mathbf{u}^T \\ u_{11}\mathbf{l} & \mathbf{l}\mathbf{u}^T+L'U' \end{pmatrix} \] と書けるので, \[ u_{11} = a_{11},\ \mathbf{u} = \mathbf{a}_1,\ \mathbf{l} = \frac{1}{u_{11}}\mathbf{a}_2 \] を計算した後, \[ A'-\mathbf{l}\mathbf{u}^T = L'U' \] という１つサイズの小さいLU分解を行えばい良いです。

以下がドゥーリトル法のプログラムです。三重ループがあるので計算量は\(\mathcal{O}(n^3)\)となります。

# -*- coding: utf-8 -*-
import numpy as np

def lu(A):
    A = np.copy(A)  # 作業用

    n = A.shape[0]
    L = np.identity(n, dtype=float)      # 単位行列
    U = np.zeros((n,n), dtype=float)     # 零行列

    for k in xrange(n):
        U[k,k] = A[k,k]   # u11 = a11
        for i in xrange(k+1,n):
            U[k,i] = A[k,i]         # u = a1
            L[i,k] = A[i,k]/U[k,k]  # l = a2/u11
        for i in xrange(k+1,n):
            for j in xrange(k+1,n):
                A[i,j] -= L[i,k] * U[k,j]   # A - lu^T
    return (L, U)

A = np.array([[4, 5, 1, 9], [2, 1, 3, 1], [3, 1, 5, 3], [-2, 1, 4, 3]], dtype=float)
L,U = lu(A)
print "A=\n",A
print "L=\n",L
print "U=\n",U
print "LU=\n", np.dot(L, U)

実装例は省略しますが,実際にはガウス消去法と同様,部分ピボット選択を行う必要があります。既存のライブラリを使う場合には,ピボット行交換の情報が一緒に返されるので注意して下さい。

例えば,SciPyの実装で\(\begin{pmatrix} 2 & 1 \\ 4 & 5\end{pmatrix}\)をLU分解すると \[ \begin{pmatrix} 2 & 1 \\ 4 & 5\end{pmatrix} = \color{yellow}{\begin{pmatrix} 0 & 1 \\ 1 & 0\end{pmatrix}} \begin{pmatrix} 1 & 0 \\ 1/2 & 1\end{pmatrix} \begin{pmatrix} 4 & 5 \\ 0 & -3/2\end{pmatrix} \] という物が帰ってきます。これは1行目と2行目を入れ替えたという事を表しています。

>>> A = np.array([[2,1],[4,5]])
>>> P,L,U = LA.lu(A)
>>> P
array([[ 0.,  1.],
       [ 1.,  0.]])
>>> L
array([[ 1. ,  0. ],
       [ 0.5,  1. ]])
>>> U
array([[ 4. ,  5. ],
       [ 0. , -1.5]])

結局LU分解に\(\mathcal{O}(n^3)\)の計算量が必要なのだから,速くなっていないのでは？と思った人もいるかもしれません。

しかし,係数行列が同一の連立方程式を大量に解く場合,一回だけ\(\mathcal{O}(n^3)\)で分解しておけば後は全て\(\mathcal{O}(n^2)\)で計算出来るということになります。

# -*- coding: utf-8 -*-
import numpy as np
import scipy
from scipy import linalg as LA
from time import time

N = 1000
A = np.random.randn(N, N)   # 係数行列
b = np.random.randn(N)      # 定数項

REPEAT = 100    # とりあえず定数項は固定して REPEAT 回 Ax=bを解いてみる

print "== LU分解不使用 =="
start = time()
for i in range(REPEAT):
    x = LA.solve(A, b)
print time()-start, "sec"

print "== LU分解使用 =="
start = time()
lu_A = LA.lu_factor(A)    # LU分解
for i in range(REPEAT):
    LA.lu_solve(lu_A, b)
print time()-start, "sec"

LU分解の発展

LU分解の特別なものにコレスキー分解という物があります。

これは固有値が全て正の対称行列に対して行う事が出来るもので,下三角行列\(L\)によって \[ \color{yellow}{ A = LL^T } \] と分解します。固有値が全て正の行列は正定値行列と呼ばれます。

\(L^T\)は上三角行列なのでこれはLU分解に他なりませんが,\(L^T\)の成分を計算する必要がないためメモリ使用量・計算量を減らす事が出来ます。

また,\(A\)が正定値でない場合には,下三角行列\(L\)と対角行列\(D\)で \[ \color{yellow}{ A = LDL^T } \] とする修正コレスキー分解という手法もあります。

更に,特殊な形の行列(帯行列と呼ばれる物など)にLU分解やコレスキー分解を行う事によって得られる公式も様々存在します。

これらの詳細は本日は省略しますが,是非勉強してみてください。

係数行列が正方でない場合

続いて,\(A\)が一般の\((m,n)\)行列の時の連立一次方程式 \[ A\mathbf{x} = \mathbf{b} \] について考えます。

\(m > n\)の時は方程式の数が変数の数より多い場合, \(m < n\)の時は方程式の数が変数の数より少ない場合となります。

求めたいパラメータの数に対して得られる情報が過多である,情報が足りないという場面は現実の問題では良くあります。

準備:ベクトルでの微分法

準備として微分学の知識が必要です。

多変数のスカラー値関数\(f(\mathbf{x}) = f(x_1,\cdots,x_n)\)の勾配を \[ \color{yellow}{\frac{\mathrm{d} f}{\mathrm{d}\mathbf{x}}} = \mathrm{grad} f = \left(\frac{\partial f}{\partial x_1}, \cdots, \frac{\partial f}{\partial x_n}\right)^T\] と書き表す事にします。

まず定数\(f(\mathbf{x}) = c\)に対しては明らかに \[ \color{yellow}{\frac{\mathrm{d} c}{\mathrm{d}\mathbf{x}} = \mathbf{0}} \] です。

定ベクトルと変数ベクトルの内積 \[ (\mathbf{a},\mathbf{x}) = \mathbf{a}^T\mathbf{x} = \mathbf{x}^T\mathbf{a} = a_1x_1 + a_2x_2 + \cdots + a_nx_n \] については, \[ \frac{\partial}{\partial x_i}(a_1x_1 + a_2x_2 + \cdots + a_nx_n) = a_i \] なので \[ \color{yellow}{\frac{\mathrm{d}}{\mathrm{d} \mathbf{x}}\mathbf{a}^T\mathbf{x} = \frac{\mathrm{d}}{\mathrm{d} \mathbf{x}}\mathbf{x}^T\mathbf{a} = \mathbf{a}} \] となります。

変数ベクトルのノルムの2乗 \[ ||\mathbf{x}||^2 = \mathbf{x}^T\mathbf{x} = x_1^2 + x_2^2 + \cdots + x_n^2 \] については \[ \frac{\partial}{\partial x_i}(x_1^2 + x_2^2 + \cdots + x_n^2) = 2x_i \] なので \[ \color{yellow}{\frac{\mathrm{d}}{\mathrm{d} \mathbf{x}}||\mathbf{x}||^2 = 2\mathbf{x}} \] となります。

また\(\mathbf{x}^TA\mathbf{x}\)という関数 \[ \mathbf{x}^TA\mathbf{x}=\mathbf{x}^T(\sum_k a_{jk}x_k)_{n,1} = \sum_j \sum_k a_{jk} x_j x_k \] については\(\partial x_j/\partial x_i = \delta_{ij}\)である事に注意すれば \[ \begin{aligned} \frac{\partial}{\partial x_i}\mathbf{x}^TA\mathbf{x} &= \sum_j \sum_k a_{jk} (\delta_{ij}x_k + \delta_{ik}x_j) \\ &= \sum_j \sum_k a_{jk} \delta_{ij}x_k + \sum_j\sum_k a_{jk} \delta_{ik}x_j \\ &= \sum_k a_{ik} x_k + \sum_j a_{ji} x_i \\ \end{aligned} \] であるので \[ \color{yellow}{\frac{\mathrm{d}}{\mathrm{d}\mathbf{x}}\mathbf{x}^TA\mathbf{x} = (A+A^T)\mathbf{x}} \] となります。

以上に明らかに成立する微分操作の線型性を加えてまとめると以下のようになります。

まず \[ \begin{aligned} f(\mathbf{x}) &= (A\mathbf{x}-\mathbf{b})^T(A\mathbf{x}-\mathbf{b}) = (\mathbf{x}^TA^T-\mathbf{b}^T)(A\mathbf{x}-\mathbf{b}) \\ &= \mathbf{x}^TA^TA\mathbf{x}-\mathbf{b}^TA\mathbf{x}-\mathbf{x}^TA^T\mathbf{b}+\mathbf{b}^T\mathbf{b} \end{aligned} \] なので \[ \begin{aligned} \frac{\mathrm{d}f}{\mathrm{d}\mathbf{x}} &= (A^TA + (A^TA)^T)\mathbf{x} - (\mathbf{b}^TA)^T - A^T\mathbf{b} \\ &= 2A^TA\mathbf{x} -2 A^T\mathbf{b} \\ &= \color{yellow}{2A^T(A\mathbf{x}-\mathbf{b})} \end{aligned} \] となります。\(((ax-b)^2)'=2a(ax-b)\)との類似性に注目しましょう。

方程式の数\( > \)変数の数の場合

話を元に戻して,方程式の数が変数の数より多い(\(A\)が縦長)場合の \[ A\mathbf{x}=\mathbf{b} \qquad\cdots(1)\] について考えます。

もちろん,一般に\((1)\)を満たす\(\mathbf{x}\)は存在しないので,誤差の平方 \[ E(\mathbf{x}) = ||A\mathbf{x}-\mathbf{b}||^2 \] が最小となる\(\mathbf{x}\)を考えます。つまり最小二乗法です。

\(E\)が極小値を取れば良いので,前頁の結果より \[ \mathrm{grad}E = \mathbf{0}\ \Leftrightarrow\ A^T(A\mathbf{x}-\mathbf{b}) = \mathbf{0}\ \Leftrightarrow\ A^TA\mathbf{x}=A^T\mathbf{b}\] を解けば良いという事になります。つまり\(A^TA\)が正則ならば \[ \color{yellow}{\mathbf{x} = (A^TA)^{-1}A^T\mathbf{b}} \] です。

>>> import numpy as np
>>> from numpy import linalg as LA
>>> A = np.array([[2,3,1],[1,2,2],[1,-1,3],[1,4,-1]],dtype=float)
>>> b = np.array([1,3,1,1],dtype=float)
>>> x = LA.solve(np.dot(A.T, A), np.dot(A.T, b))    # A^TAx = A^Tbを解く
>>> x
array([-1.76470588,  1.04705882,  1.30588235])
>>> np.dot(A, x)
array([ 0.91764706,  2.94117647,  1.10588235,  1.11764706])
>>> LA.norm(np.dot(A, x) - b)     # ||Ax-b||
0.18786728732554484

方程式の数\( < \)変数の数の場合

方程式の数が変数の数より少ない(\(A\)が横長)場合の \[ A\mathbf{x}=\mathbf{b} \qquad\cdots(1)\] について考えます。

この場合は一般に解が無数に存在し,一意に定まりません。そこで適当なスカラー値関数\(f(\mathbf{x})\)に対して \[ \color{yellow}{A\mathbf{x}=\mathbf{b}\text{を満たす$\mathbf{x}$のうち}f(\mathbf{x})\text{が最小であるもの}} \] を求めるという問題を考えます。

準備: ラグランジュの未定乗数法

このように \[ \text{条件}g(\mathbf{x}) = 0\text{の時の}f(\mathbf{x})\text{の極値} \] を考える場合にはラグランジュの未定乗数法という手法を使います。

元の問題に戻りましょう。例題として \[ A\mathbf{x}=\mathbf{b}\text{で}||\mathbf{x}||^2\text{が最小のもの} \] を求めてみます。

ラグランジュの未定乗数法より \[ F(\mathbf{x},\mathbf{\lambda}) = \frac{1}{2}||\mathbf{x}||^2-\mathbf{\lambda}^T(A\mathbf{x}-\mathbf{b}) \] の勾配が\(\mathbf{0}\)であれば良いので \[ \begin{aligned} &\frac{\mathrm{d} F}{\mathrm{d}\mathbf{x}} = \mathbf{0},\ \frac{\mathrm{d} F}{\mathrm{d}\mathbf{\lambda}} = \mathbf{0} \Leftrightarrow \mathbf{x}-A^T\mathbf{\lambda}=\mathbf{0},\ A\mathbf{x}-\mathbf{b}=\mathbf{0}\\ \Leftrightarrow& \mathbf{x}=A^T\mathbf{\lambda},\ A\mathbf{x}=\mathbf{b} \Leftrightarrow \mathbf{x}=A^T\mathbf{\lambda},\ AA^T\mathbf{\lambda}=\mathbf{b}\\ \end{aligned} \] となります。従って\(AA^T\)が正則ならば第２式より \[ \mathbf{\lambda} = (AA^T)^{-1}\mathbf{b} \] となるので,第1式より \[ \color{yellow}{\mathbf{x} = A^T(AA^T)^{-1}\mathbf{b}} \] が求める解となります。

>>> import numpy as np
>>> from numpy import linalg as LA
>>> A = np.array([[2,1,1],[1,3,2]],dtype=float)
>>> b = np.array([1,-1])
>>> lam = LA.solve(np.dot(A, A.T), b)   # lam = (AA^T)^{-1}b
>>> x = np.dot(A.T, lam)                # x = A^T lam
>>> x
array([ 0.82857143, -0.51428571, -0.14285714])
>>> np.dot(A, x)
array([ 1., -1.])

まとめましょう

問題によって「\( ||A\mathbf{x}-\mathbf{b}|| \)が最小」とか「\(||\mathbf{x}||\)が最小」という条件は変化します。よって,上の公式をただ覚えるのではなくて以上の議論をしっかり理解する事が大切です。

また,\(A^TA\)とか\(AA^T\)が正則じゃない時はどうするのか？という疑問があると思います。興味のある方はペンローズの擬似逆行列について調べてみましょう。

今回はここで終わります。

まだ,反復解法と解法の安定性の話をしていませんが,連立一次方程式の基本的な話題は概ね終了です。

次回は線形空間・線型写像・固有値・固有ベクトルの理論を中心に進めて行きます。