2. 連立一次方程式

2.1 概要

例えば, 未知数が 2 つの連立一次方程式 (2 元連立一次方程式) は次のようなものでした.
\[
a_{11}x_1 + a_{12}x_2 = b_1 \\
a_{21}x_1 + a_{22}x_2 = b_2 \\
\] \(x_1\) と \(x_2\) は未知数です.

これを行列を使って表すと次のようになります.
\[
\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}
\] ただし,
\[
\boldsymbol{A} =
\begin{pmatrix}
a_{11} & a_{12} \\
a_{21} & a_{22} \\
\end{pmatrix},
\boldsymbol{b} =
\begin{pmatrix}
b_1 \\
b_2 \\
\end{pmatrix},
\boldsymbol{x} =
\begin{pmatrix}
x_1 \\
x_2 \\
\end{pmatrix}
\] \(\boldsymbol{A}\) を係数行列, \(\boldsymbol{b}\) を右辺ベクトル, \(\boldsymbol{x}\) を解ベクトルといいます. 未知数が n 個の連立一次方程式 (n 元連立一次方程式) は行列とベクトルを使った表記は同じになります.

逆行列 \(\boldsymbol{A^{-1}}\) がわかれば, 左から方程式の両辺に掛けてやると,
\[
\boldsymbol{A^{-1}Ax = A^{-1}b}
\] となり, \(\boldsymbol{A^{-1}A}\) が消えて次のように \(\boldsymbol{x}\) を求めることができます.
\[
\boldsymbol{x} = \boldsymbol{A^{-1}}\boldsymbol{b}
\] Excel では MINVERSE という逆行列を求めるワークシート関数があるので, これを使って逆行列を求め, 右辺ベクトル \(\boldsymbol{b}\) に左から掛けてやることにより解ベクトル \(\boldsymbol{x}\) を求めることができます. しかし, この方法は計算量が多くなることと精度が悪くなりやすいことから, 逆行列を求めないで直接計算する方がよいとされています.

2.1 連立一次方程式の解法

コンピュータで連立一次方程式を解く一般的な方法は本質的には手計算を行うときと同じです. すなわち, 2 本の式の片方に定数を掛けてもう一方から引くことにより変数を 1 つずつ消去していくことを繰り返します. これはガウスの消去法と呼ばれるものです.

これを行列の式で表すと係数行列を三角行列の積に分解していく過程と同じになります. 三角行列を係数とする連立一次方程式は代入操作だけで簡単に解くことができます.

2.1.1 一般行列の場合

まず, 係数行列 \(\boldsymbol{A}\) を次のように分解します. これを LU 分解と呼びます.
\[
\boldsymbol{A} = \boldsymbol{LU} \space (ただし, \boldsymbol{L} は下三角行列, \boldsymbol{U} は上三角行列)
\] LU 分解の計算手順はいくつかのバリエーションがありますが基本的にはガウスの消去法に基づいたものです. なお, この過程を前進消去といいます.

LU 分解の計算中に対角要素が 0 になると計算を進めることができなくなります. このような場合には, 対角要素が 0 でない残りの行のどれかと入れ替えを行えばよいことになります. なお, 誤差をできるだけ少なくするという意味では対角要素の値が最大の行と入れ替えるのがよさそうです. このあたりも種々のやり方が考えられます. このような入れ替えを行うことをピボットの選択といいます.

次に, 以下のように \(\boldsymbol{y} = \boldsymbol{Ux}\) とおいて, \(\boldsymbol{Ly} = \boldsymbol{b}\) から \(\boldsymbol{y}\) を求め, 最後に \(\boldsymbol{Ux} = \boldsymbol{y}\) から \(\boldsymbol{x}\) を求めます. \(\boldsymbol{L}\) と \(\boldsymbol{U}\) は三角行列なので代入操作だけで容易に解を計算することができます. これを後退代入といいます.
\[
\begin{align}
&\boldsymbol{LUx} = \boldsymbol{b} \\
&\boldsymbol{Ly} = \boldsymbol{b} \\
&\boldsymbol{Ux} = \boldsymbol{y} \\
\end{align}
\] なお, 分解(前進消去)と解の計算(後退代入)を分けておくことにより, 同じ係数行列で右辺だけ異なる連立一次方程式を解くときに後退代入の計算だけ行えばよいという利点があります.

2.1.2 対称行列の場合

まず, 係数行列 \(\boldsymbol{A}\) を次のように分解します. これをコレスキー分解と呼びます.
\[
\boldsymbol{A} = \boldsymbol{LL^T} または \boldsymbol{A} = \boldsymbol{U^TU} \space (ただし \boldsymbol{L} は下三角行列, \boldsymbol{U} は上三角行列)
\] コレスキー分解の計算手順は LU 分解と同様ですが, 対角要素が 0 になる対称行列が現れることは応用上ないので, 通常はピボットの選択は行いません.

分解が完了したら一般行列の場合と同様に三角行列の代入操作により解を求めます.

2.1.3 条件数について

連立一次方程式を解く際には係数行列の条件数を求めておくと有益です.

条件数は次のように定義されます.
\[
cond(\boldsymbol{A}) = ||\boldsymbol{A}||\cdot||\boldsymbol{A^{-1}}||
\] 条件数の値は, 最も条件が良い単位行列 \(\boldsymbol{I}\) の場合 1, 最も条件が悪い特異行列の場合 \(\infty\) で, 普通はその間の値となります. 条件数は定義どおりに計算すると逆行列を求めなければならず計算量が多くなるため, 通常は簡便に計算できる推定値が使わます.

これを使って, 求められた解 \(\boldsymbol{x}\) の相対誤差を次のようにして推定することができます.
\[
||\boldsymbol{x^*} – \boldsymbol{x}||/||\boldsymbol{x}|| \leq cond(\boldsymbol{A})\cdot \epsilon
\] ここで, \(\boldsymbol{x^*}\) は真の解, \(\epsilon\) は計算機イプシロン(*)です. 条件数が非常に大きい場合を悪条件であるといい, 例えば条件数が \(10^{15}\) 程度であればほとんど有効桁がない可能性があることがわかります.

(*) 計算機イプシロンは実数を計算機の浮動小数で近似したための誤差で, Excel (IEEE形式の倍精度) ではおおよそ \(1.11 \times 10^{-16}\) となります

2.2 XLPack を使った連立一次方程式の解き方

LU 分解およびコレスキー分解の XLPack の VBA サブルーチンはそれぞれ Dgesv および Dposv です. ワークシート関数はそれぞれ, WDgesv, WDposvです.

VBA プログラムで条件数の推定値を求めるためには Dlange と Dgecon および Dlansy と Dpocon が使用できます. ワークシート関数の場合には, WDgesv および WDposv の計算結果の一部として条件数の推定値も返されます.

これらは, 線形計算ライブラリ LAPACK のサブルーチン DGESV, DLANGE, DGECON (LU 分解) および DPOSV, DLANSY, DPOCON (コレスキー分解) を使用しています.

例題

次の連立一次方程式を解く.
\[
\begin{align}
& 10x_1 – 7x_2 = 7 \\
& -3x_1 + 2x_2 + 6x_3 = 4 \\
& 5x_1 – x_2 + 5x_3 = 6 \\
\end{align}
\] これは, 行列で表すと次のようになります.
\[
\boldsymbol{A}\boldsymbol{x} = \boldsymbol{b}
\] ただし,
\[
\boldsymbol{A} =
\begin{pmatrix}
10 & -7 & 0 \\
-3 & 2 & 6 \\
5 & -1 & 5 \\
\end{pmatrix},
\boldsymbol{b} =
\begin{pmatrix}
7 \\
4 \\
6 \\
\end{pmatrix},
\boldsymbol{x} =
\begin{pmatrix}
x_1 \\
x_2 \\
x_3 \\
\end{pmatrix}
\] ここで, \(\boldsymbol{A}\) は係数行列, \(\boldsymbol{b}\) は右辺ベクトル, \(\boldsymbol{x}\) は解ベクトルです.

2.2.1 ワークシート関数を使用した解き方

まず, ワークシートの適当な場所に係数行列 \(\boldsymbol{A}\) と右辺ベクトル \(\boldsymbol{b}\) のデータを入力します (オレンジ色のセル). そして, 解ベクトルを入れる場所として(解の要素数 + 2)個のセルを選択してワークシート関数 WDgesv を入力します (緑色のセル).

WDgesv の必要なパラメータは N, A, B, Nrhs です.
N: 要素数 (この例では 3).
A: 係数行列 \(\boldsymbol{A}\) のセル範囲.
B: 右辺ベクトル \(\boldsymbol{b}\) のセル範囲.
Nrhs: 右辺ベクトルの本数 (\(\boldsymbol{b}\) だけ変えて複数回計算する場合に入力. 省略時は 1 とみなす).

入力が終了したら Ctrl + Shift + Enter を押します.

これで, 解 \(\boldsymbol{x} = \begin{pmatrix} 0 & -1 & 1\end{pmatrix}^T\) が求められました. なお, 表示されている 0.078283 は条件数の推定値の逆数, 一番下の 0 はリターンコード (正常終了を示す) です.

2.2.2 VBA プログラムを使用した解き方

上と同じ例題を VBA プログラムにより解きます. Dgesv を使ったプログラム例を示します.

Sub Start()
    Const NMax = 10
    Dim N As Long, A(NMax, NMax) As Double, B(NMax) As Double
    Dim Anorm As Double, RCond As Double
    Dim Ipiv(NMax) As Long, Info As Long, I As Long, J As Long
    '--- Input data
    N = 3
    For I = 0 To N - 1
        For J = 0 To N - 1
            A(I, J) = Cells(5 + I, 1 + J)
        Next
        B(I) = Cells(5 + I, 4)
    Next
    '--- Solve equation
    Anorm = Dlange("1", N, N, A())
    Call Dgesv(N, A(), Ipiv(), B(), Info)
    If Info = 0 Then
        '---  Compute condition number
        Call Dgecon("1", N, A(), Anorm, RCond, Info)
    End If
    '--- Output solution
    For I = 0 To N - 1
        Cells(5 + I, 5) = B(I)
    Next
    Cells(8, 5) = RCond
    Cells(9, 5) = Info
End Sub

プログラム中, Dlange と Dgecon は条件数を求めるためのもので方程式を解くだけならば不要です.

所定の位置(下のオレンジ色のセル)に係数行列 \(\boldsymbol{A}\) と右辺ベクトル \(\boldsymbol{b}\) の値を入力し, マクロ Start としてプログラムを実行すると次の結果が得られます. ワークシート関数を使用したときと異なり, 係数行列 \(\boldsymbol{A}\) と右辺ベクトル \(\boldsymbol{b}\) の値を入力しただけでは結果が得られず, VBA プログラムを実行してやる必要があります.

2.2.3 係数行列が対称な場合

係数行列 \(\boldsymbol{A}\) が正定値対称行列の場合には WDposv と Dposv を使った方が効率よく計算を行うことができます. 手順は上と同様ですが, A のデータは上または下三角部分(対角要素を含む)だけ設定すればよくなります. 残りの部分は参照されません.

ここには掲載しませんが, ダウンロードできる例題ワークシートには WDposv と Dposv を使った例が入っています.