連立一次方程式の数値計算
4．連立方程式の計算方法
4．1．連立方程式の解法
連立一次方程式の解を行列演算によって求める方法はいろいろあります。 以下に主な方法を紹介します。（1）逆行列を使って求める方法
行列\( A \)の逆行列\( A^{-1} \)を求めることで、変数ベクトル\( \bf{x} \)を求めます。高校で習うのは主にこの方法です。
\[
A \bf{x}
=
\bf{b}
\rightarrow
\bf{x}
=
A^{-1} \bf{b}
\]
しかしながら数値計算で逆行列を求めることは、計算精度、安定性やメモリ占有の観点からほとんどありません。
また、係数行列が正方行列の場合にしか対応できず、最適化計算等への対応もできません。
従って連立方程式の解法としては次以降に示す行列分解による方法が用いられます。
（2）LU分解を使って求める方法
行列\( A \)を、基本変形変換を用いて下三角行列\( L \)と上三角行列\( U \)に分解することで、変数ベクトル\( \bf{x} \)を求めます。\( A = LU \)とできれば、連立方程式は次のように書き換えられます。
\[
LU \bf{x} = \bf{b}
\ \rightarrow \
U \bf{x} = L^{-1} \bf{b}
\]
これを成分表示すると
\[
\begin{pmatrix}
u_{11} & \cdots & u_{1n} &
\\
& \ddots & \vdots &
\\
& & u_{nn} &
\end{pmatrix}
\left(
\begin{array}{c}
x_1
\\
\vdots
\\
x_n
\end{array}
\right)
=
\left(
\begin{array}{c}
b'_1
\\
\vdots
\\
b'_n
\end{array}
\right)
\]
となるので、n行目から順にn-1行目、・・・、1行目と計算していけば変数ベクトル\( \bf{x} \)は求まります
（この方法を後退代入法と呼びます）。
また、係数行列\( A \)は正方行列でなくても計算できます。
LU分解の詳細はLU分解ページで説明します。
（3）QR分解を使って求める方法
行列\( A \)をユニタリ行列\( Q \)と上三角行列\( R \)に分解することで、変数ベクトル\( \bf{x} \)を求めます。\( A = QR \)とできれば、連立方程式は次のように書き換えられます。
\[
QR \bf{x} = \bf{b}
\ \rightarrow \
R \bf{x} = Q^{-1} \bf{b}
\]
これを成分表示すると
\[
\begin{pmatrix}
r_{11} & \cdots & r_{1n} &
\\
& \ddots & \vdots &
\\
& & r_{nn} &
\end{pmatrix}
\left(
\begin{array}{c}
x_1
\\
\vdots
\\
x_n
\end{array}
\right)
=
\left(
\begin{array}{c}
b'_1
\\
\vdots
\\
b'_n
\end{array}
\right)
\]
となるので、LU分解同様、n行目から順に計算していけば変数ベクトル\( \bf{x} \)は求まります。
また、係数行列\( A \)は正方行列でなくても計算できます。
QR分解は計算安定性が高いので、最小二乗法の計算等に用いられます。
QR分解の詳細はQR分解ページで説明します。
（4）反復法を使って求める方法
大規模な連立方程式を解く必要がある場合、上記（2）、（3）の方法では計算負荷が非常に高くなり、実用性がなくなります。 このような場合、適当な初期解を与えて繰り返し計算を行うことで、真の解に収束させていく方法をとります。 この計算方法が反復法と呼ばれるゆえんです。 それに対し、（2）、（3）は反復法に対して直接法と呼びます。反復法は、必ずしも解が得られる（収束する）わけではありません。 そのため、解が収束するための条件を満たす必要があります。 それは、係数行列の対角成分の絶対値が、同じ行の非対角成分の絶対値の和より大きいこと、つまり対角優位であることが条件です。
4．2．後退代入
連立方程式を解く際に係数行列を上三角化した場合（LU分解、QR分解など）、実際に計算する部分を特定すると、次のようになります。 なお、係数行列\( A \in R^{m \times n} \ \)、定数項ベクトル\( \bf{b} \in R^m \)、変数ベクトル\( \bf{x} \in R^n \)、行列\( A \)の階数（ランク）を\( r \)とします。 また、上三角化した行列を\(U\)で表します。-
\( r = n \leq m \land \bf{b_{m-n}} = \bf{o} \)の場合
\[ \begin{pmatrix} U_{nn} & \\ O_{m-n,n} & \end{pmatrix} \bf{x_n} = \begin{pmatrix} \bf{b_n} \\ \bf{b_{m-n}} \end{pmatrix} \ \rightarrow \ U_{nn} \bf{x_n} = \bf{b_n} \]計算する部分だけを成分表示すると、次のようになります。
\[ \begin{pmatrix} u_{11} & \cdots & u_{1n} \\ & \ddots & \vdots \\ & & u_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ \vdots \\ b_n \end{pmatrix} \tag{4.2-1} \] -
\( r \leq m \lt n \)または\( r \lt n \lt m \land \bf{b_{m-r}} = \bf{o} \)の場合
\[ \begin{pmatrix} U_{rr} & O_{r,n-r} & \\ O_{m-r,r} & O_{n-r,n-r} & \end{pmatrix} \begin{pmatrix} \bf{x_r} \\ \bf{x_{m-r}} \end{pmatrix} = \begin{pmatrix} \bf{b_r} \\ \bf{b_{m-r}} \end{pmatrix} \ \rightarrow \ U_{rr} \bf{x_r} = \bf{b_r} \]このケースでは、3.2節でみたように定数項ベクトルから一般解の任意係数倍を加える必要があります。 このとき特殊解は、この任意係数をすべて0と置くことで求めることができます。
結局、計算する部分だけを成分表示すると、次のようになります。
\[ \begin{pmatrix} u_{11} & \cdots & u_{1r} \\ & \ddots & \vdots \\ & & u_{rr} \end{pmatrix} \begin{pmatrix} x_1 \\ \vdots \\ x_r \end{pmatrix} = \begin{pmatrix} b_1 \\ \vdots \\ b_r \end{pmatrix} \tag{4.2-2} \]
この場合、添字の大きい方から順番に計算することで\( \bf{x} \)を特定できます。
この計算方法を後退代入と呼びます。
ピボッティングとスケーリングは基本変形行列によって実行可能です。
例えば、1列目は\( |a_{31}| \)が最大をとる場合、3行目と1行目を入れ替えます。
例えば、1行目は\( |a_{13}| \)が最大をとる場合、3列目と1列目を入れ替えます。
例えば、2行3列成分を“1”とする2行目へのスケーリングは、次のようになります。
スケーリングやピボッティングを実施する必要のない場合は単位行列\( E \)を適用します。
4．3．ピボッティングとスケーリング
前節の後退代入では、上三角行列の対角成分は0でないことを前提に話を進めました。 しかしこれは保証されず、対角成分=0のとき計算不能となります。 また、対角成分の絶対値が非常に小さい場合、割り算による計算精度上の問題が生じます。 これらの問題を回避するために、“対角成分≠0”となるよう、行または列の入れ替え操作を行います。 この入れ替え操作をピボッティングと呼び、以下に示す方法で行います。- 行ピボッティング：\( k \)列目\( k+1～n \)行目の中から絶対値最大となる成分を含む\( i \)行と\( k \)行を交換
- 列ピボッティング：\( k \)行目\( k+1～n \)列目の中から絶対値最大となる成分を含む\( j \)列と\( k \)列を交換
- 完全ピボッティング：行ピボッティングと列ピボッティングを両方行う
ピボッティングとスケーリングは基本変形行列によって実行可能です。
（1）行ピボッティングの場合
\( i \)行と\( j \)行を入れ替えるために、 交換行列\( P(i,j) \)を左から掛けます。例えば、1列目は\( |a_{31}| \)が最大をとる場合、3行目と1行目を入れ替えます。
\[
\begin{eqnarray}
& P(i,j) A \quad
& =
\begin{pmatrix}
0 & 0 & \txc{red}{1} & 0
\\
0 & 1 & 0 & 0
\\
\txc{red}{1} & 0 & 1 & 0
\\
0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
a_{11} & a_{12} & a_{13} & a_{14}
\\
a_{21} & a_{22} & a_{23} & a_{24}
\\
a_{31} & a_{32} & a_{33} & a_{34}
\\
a_{41} & a_{42} & a_{43} & a_{44}
\end{pmatrix}
\\
\\
& & =
\begin{pmatrix}
\txc{red}{a_{31}} & \txc{red}{a_{32}} & \txc{red}{a_{33}} & \txc{red}{a_{34}}
\\
a_{12} & a_{22} & a_{23} & a_{24}
\\
\txc{red}{a_{11}} & \txc{red}{a_{12}} & \txc{red}{a_{13}} & \txc{red}{a_{14}}
\\
a_{41} & a_{42} & a_{43} & a_{44}
\end{pmatrix}
\end{eqnarray}
\]
（2）列ピボッティングの場合
\( i \)列と\( j \)列を入れ替えるために、 交換行列\( P(i,j) \)を右から掛けます。例えば、1行目は\( |a_{13}| \)が最大をとる場合、3列目と1列目を入れ替えます。
\[
\begin{eqnarray}
& P(i,j) A \quad
& =
\begin{pmatrix}
a_{11} & a_{12} & a_{13} & a_{14}
\\
a_{21} & a_{22} & a_{23} & a_{24}
\\
a_{31} & a_{32} & a_{33} & a_{34}
\\
a_{41} & a_{42} & a_{43} & a_{44}
\end{pmatrix}
\begin{pmatrix}
0 & 0 & \txc{red}{1} & 0
\\
0 & 1 & 0 & 0
\\
\txc{red}{1} & 0 & 1 & 0
\\
0 & 0 & 0 & 1
\end{pmatrix}
\\
\\
& & =
\begin{pmatrix}
\txc{red}{a_{13}} & a_{12} & \txc{red}{a_{11}} & a_{14}
\\
\txc{red}{a_{23}} & a_{22} & \txc{red}{a_{21}} & a_{24}
\\
\txc{red}{a_{33}} & a_{32} & \txc{red}{a_{31}} & a_{34}
\\
\txc{red}{a_{43}} & a_{42} & \txc{red}{a_{41}} & a_{44}
\end{pmatrix}
\end{eqnarray}
\]
（3）スケーリングの場合
\( i \)行目を\( 1/a_{ij} \)倍するために、 スカラー倍行列\( M(i:1/a_{ij}) \)を左から掛けます。例えば、2行3列成分を“1”とする2行目へのスケーリングは、次のようになります。
\[
\begin{eqnarray}
M(i:1/a_{ij}) A
& =
\begin{pmatrix}
1 & 0 & 0 & 0
\\
0 & \txc{red}{1/a_{23}} & 0 & 0
\\
0 & 0 & 1 & 0
\\
0 & 0 & 0 & 1
\end{pmatrix}
\begin{pmatrix}
a_{11} & a_{12} & a_{13} & a_{14}
\\
a_{21} & a_{22} & a_{23} & a_{24}
\\
a_{31} & a_{32} & a_{33} & a_{34}
\\
a_{41} & a_{42} & a_{43} & a_{44}
\end{pmatrix}
\\
\\
& =
\begin{pmatrix}
a_{11} & a_{12} & a_{13} & a_{14}
\\
\txc{red}{a_{12}/a_{23}} & \txc{red}{a_{22}/a_{23}} & \txc{red}{a_{23}/a_{23}} & \txc{red}{a_{24}/a_{23}}
\\
a_{31} & a_{32} & a_{33} & a_{34}
\\
a_{41} & a_{42} & a_{43} & a_{44}
\end{pmatrix}
\end{eqnarray}
\]
