関数LU_decomposition マニュアル

(The documentation of function LU_decomposition)
Last Update: 2021/12/6


◆機能・用途(Purpose)

正方行列をLU分解する。
Calculate the LU decomposition of a square matrix.

◆形式(Format)

#include <matrix/inverse.h>
inline struct LU LU_decomposition(const struct matrix A)


◆引数(Arguments)

A LU分解をしたい行列\(\myvector{A}\)。正方行列のみ可。
The matrix \(\myvector{A}\) for which the LU decomposition is to be calculated. This must be a square matrix.


◆戻り値(Return value)

行列\(\myvector{A}\)のLU分解。 厳密には\(\myvector{A}\)自体のLU分解ではなく、 \(\myvector{A}\)の行の交換(ピボット選択)を\(N-1\)回行って得られる行列 \(\myvector{A^{(N-1)}}\)(下記「計算式」参照)のLU分解を求める。 但し\(N\)は\(\myvector{A}\)のサイズを表す。
The LU decomposition of matrix \(\myvector{A}\). Strictly, not the LU decomposition of \(\myvector{A}\) itself but the LU decomposition of \(\myvector{A^{(N-1)}}\) obtained by swapping the rows of \(\myvector{A}\) (pivot selection) for \(N-1\) times; see “Formula” below for detail. Here \(N\) represents the size of \(\myvector{A}\).
戻り値のメンバの値は以下のようになる。
The values of members of the return value are as follows.
戻り値のメンバ
Member of the return value 		値
Value
L LU分解で得られた下三角行列\(\myvector{L}\)。
The lower triangular matrix \(\myvector{L}\) obtained by the LU decomposition.
U LU分解で得られた上三角行列\(\myvector{U}\)。
The lower triangular matrix \(\myvector{U}\) obtained by the LU decomposition.
各\(i\)に対するpivot[i]
pivot[i] for each \(i\) 		行列\(\myvector{A^{(N-1)}}\)の第\(i\)行成分が もともと行列\(\myvector{A}\)のときに代入されていた行番号。 行列\(\myvector{A}\)の\((i,j)\)成分を\(a_{ij}\)、 行列\(\myvector{A^{(N-1)}}\)の\((i,j)\)成分を\(a_{ij}^{(N-1)}\) としたとき、例えば\(a_{ij}^{(N-1)}=a_{Ij}\)であれば pivot[i]\(=I\)である。
The row index of \(\myvector{A}\) corresponding to the \(i\)th row of \(\myvector{A^{(N-1)}}\). For example, pivot[i]\(=I\) if \(a_{ij}^{(N-1)}=a_{Ij}\), where \(a_{ij}\) and \(a_{ij}^{(N-1)}\) represent \((i,j)\)th component of matrices \(\myvector{A}\) and \(\myvector{A^{(N-1)}}\), respectively.


◆使用例(Example)

struct matrix A;
#define FORCE_LU_DECOMPOSITION TRUE
struct LU A_LU=LU_decomposition(A);


◆補足(Additional notes)

LU分解が定義できるのは\(\myvector{A}\)が正則な場合のみである。 非正則行列の場合にはプログラムをエラーとして終了する。 しかしマクロFORCE_LU_DECOMPOSITIONをTRUEにした場合には 非正則行列であっても計算を続行する。 非正則行列の判定条件や計算を続行する場合に用いる値など、 詳しくは下記「4. 非正則行列とゼロ割りの問題の扱い」を参照のこと。
The LU decomposition can be defined only when \(\myvector{A}\) is a regular matrix. In case of an irregular matrix, the program finishes as an error. However, users can choose to continue the calculation by setting a macro FORCE_LU_DECOMPOSITION to TRUE. For more detail (e.g., judgement of whether the matrix is regular or irregular, and the values used when the calculation is forced to continue), see “4. Treatments of irregular matrices and zero-division problems” below.

◆計算式(Formula)

以下ではコンピュータ・プログラムと添字を揃えるために 行列の行番号・列番号は0から始まるものとして定式化する。
To be consistent with the computer program, the formula below will use the row and column indices of the matrices starting from zero.

1. LU分解の方程式 (Equation for the LU decomposition)

正方行列\(\myvector{A}\)を下三角行列\(\myvector{L}\)と 上三角行列\(\myvector{U}\)の積 \[\begin{equation} \myvector{A}=\myvector{L}\myvector{U} \label{eq.decomposition} \end{equation}\] に分解することを考える。 \(\myvector{A}\)の\((i,j)\)成分を\(a_{ij}\)、 \(\myvector{L}\)の\((i,j)\)成分を\(l_{ij}\)、 \(\myvector{U}\)の\((i,j)\)成分を\(u_{ij}\)として、 (\ref{eq.decomposition})式を成分で書き下すと \[\begin{equation} a_{ij}=\sum_{k=0}^{N-1} l_{ik}u_{kj} \label{eq.component.all} \end{equation}\] と書ける。\(\myvector{L}\)は下三角行列、\(\myvector{U}\)は上三角行列であるので \[\begin{equation} l_{ik}=0 \hspace{1em} (i<k) \label{eq.zero_l} \end{equation}\] \[\begin{equation} u_{kj}=0 \hspace{1em} (k>j) \label{eq.zero_u} \end{equation}\] である。したがって(\ref{eq.component.all})式の \(k\)についての和を取る範囲を狭めることができて \[\begin{equation} a_{ij}=\sum_{k=0}^{min\{i,j\}} l_{ik}u_{kj} \label{eq.component} \end{equation}\] と書ける。 基本的には(\ref{eq.component})を\(l_{ik}\), \(u_{kj}\)について解けば LU分解が求まることになる。
Let us consider to decompose a square matrix \(\myvector{A}\) into a lower triangular matrix \(\myvector{L}\) and an upper triangular matrix \(\myvector{U}\) as in the form of eq. (\ref{eq.decomposition}). This equation can be rewritten by components as eq. (\ref{eq.component.all}), where \(a_{ij}\), \(l_{ij}\), and \(u_{ij}\) represent the \((i,j)\)th components of the matrices \(\myvector{A}\), \(\myvector{L}\), and \(\myvector{U}\), respectively. Since \(\myvector{L}\) and \(\myvector{U}\) are lower and upper triangular matrices, respectively, eqs. (\ref{eq.zero_l}) and (\ref{eq.zero_u}) hold. Thus the range of \(k\) for the summation in eq. (\ref{eq.component.all}) can be reduced, giving eq. (\ref{eq.component}). Basically, the LU decomposition is obtained by solving eq. (\ref{eq.component}) with respect to \(l_{ik}\) and \(u_{kj}\).

2. 制約条件 (Constraints)

(\ref{eq.component})式を\(l_{ik}\),\(u_{kj}\)について解くことを考えると 式の数に比べ未知数の方が多くなっている。 式の数は(\ref{eq.component})で\(i,j=0,\cdots,N-1\)とした\(N^2\)個である。 未知数の数は 行列\(\myvector{L}\)の下三角成分\(l_{ik}\) (\(i\geq k\))が\(N(N+1)/2\)個、 行列\(\myvector{U}\)の上三角成分\(u_{kj}\) (\(k\leq j\))が\(N(N+1)/2\)個、 合わせて\(N(N+1)=N^2+N\)個である。 つまり未知数の方が\(N\)個多い。 そこで\(\myvector{L}\)の対角成分を1とおく。 すなわち \[\begin{equation} l_{ii}=1 \label{eq.constraint} \end{equation}\] とする。こうすれば未知数は\(N^2\)個となって式の数と一致する。
Eq. (\ref{eq.component}) has more unknowns than the number of equations. Indeed, the number of equations is \(N^2\), corresponding to \(i,j=0,\cdots,N-1\) in eq. (\ref{eq.component}). The number of unknowns is \(N(N+1)=N^2+N\), composed of \(N(N+1)/2\) unknowns for \(l_{ik}\) (\(i\geq k\)) corresponding to the lower triangular components of \(\myvector{L}\) and additional \(N(N+1)/2\) unknowns for \(u_{kj}\) (\(k\leq j\)) corresponding to the upper triangular components of \(\myvector{U}\). Thus the number of unknowns is by \(N\) larger than the number of equations. To solve this situation, the diagonal components of \(\myvector{L}\) are set to 1 (eq. \ref{eq.constraint}). Then there are \(N^2\) remaining unknowns, whose number is equal to the number of equations.

3. 方程式の解 (Solution of the equation)

制約条件(\ref{eq.constraint})のもとで 方程式(\ref{eq.component})の解を具体的に書き下すことを考える。 \(i \leq j\)のとき、(\ref{eq.constraint})を用いると(\ref{eq.component})式は \[\begin{eqnarray} a_{ij} &=& \sum_{k=0}^i l_{ik}u_{kj} \nonumber \\ &=& \sum_{k=0}^{i-1} l_{ik}u_{kj} + l_{ii}u_{ij} \nonumber \\ &=& \sum_{k=0}^{i-1} l_{ik}u_{kj} + u_{ij} \label{eq.component.right} \end{eqnarray}\] と書ける。これを\(u_{ij}\)について解けば \[\begin{equation} u_{ij}=a_{ij}-\sum_{k=0}^{i-1} l_{ik}u_{kj} \label{eq.u} \end{equation}\] が得られる。 次に、\(i>j\)のとき、(\ref{eq.component})式は \[\begin{eqnarray} a_{ij} &=& \sum_{k=0}^j l_{ik}u_{kj} \nonumber \\ &=& \sum_{k=0}^{j-1} l_{ik}u_{kj} + l_{ij}u_{jj} \label{eq.component.left} \end{eqnarray}\] と書ける。これを\(l_{ij}\)について解けば \[\begin{equation} l_{ij}=\frac{1}{u_{jj}}\left(a_{ij}-\sum_{k=0}^{j-1} l_{ik}u_{kj}\right) \label{eq.l} \end{equation}\] が得られる。(\ref{eq.u})(\ref{eq.l})式を見ると、 右辺に登場する\(l_{ik}\)や\(u_{kj}\)は 左辺の要素よりも左上の要素のみであることが分かる。 したがって、(\ref{eq.u})(\ref{eq.l})式を繰り返し用いれば \(l_{ik}\)や\(u_{kj}\)の全要素を左上から順に求めていくことができる。
We now consider to write down the solution of eq. (\ref{eq.component}) under the constraint (\ref{eq.component}). When \(i \leq j\), eq. (\ref{eq.component}) can be rewritten as (\ref{eq.component.right}), where the constraint (\ref{eq.constraint}) was used. Eq. (\ref{eq.component.right}) can be solved for \(u_{ij}\) as (\ref{eq.u}). When \(i>j\), eq. (\ref{eq.component}) can be rewritten as (\ref{eq.component.left}), the solution of which with respect to \(l_{ij}\) is (\ref{eq.l}). Note that all the \(l_{ik}\) and \(u_{kj}\) values that appear in the right hand sides of eqs. (\ref{eq.u}) and (\ref{eq.l}) are located at upper and left sides than those that appear in the left sides. This means that all the values of \(l_{ik}\) and \(u_{kj}\) can be calculated from left top to right bottom by repeatedly using eqs. (\ref{eq.u}) and (\ref{eq.l}).

4. 計算を安定させるための行の交換(ピボット選択) (Pivot selection for stabilizing the calculation)

(\ref{eq.l})式では分母に\(u_{jj}\)が登場するので これが小さな値になると計算誤差が増大する。 これを避けるためにLU分解では 行列\(\myvector{A}\)の行の交換(部分ピボット選択)を行う。 交換する行の具体的な選択方法については次節で述べるが、 基本的な考え方としては\(|u_{jj}|\)を最大化するように選ぶ。 以下、行列\(\myvector{A}\)から行の交換を \(m\)回行って得られる行列を\(\myvector{A^{(m)}}\)、 その\((i,j)\)成分を\(a_{ij}^{(m)}\)と書くことにする。 この関数では\(\myvector{A}\)自体のLU分解ではなく \(\myvector{A^{(N-1)}}\)のLU分解を求める。 すなわち、厳密には(\ref{eq.decomposition})ではなく \[\begin{equation} \myvector{A^{(N-1)}}=\myvector{L}\myvector{U} \label{eq.decomposition.pivot} \end{equation}\] を解く。したがって解も(\ref{eq.u})(\ref{eq.l})ではなく \[\begin{equation} u_{ij}=a_{ij}^{(N-1)}-\sum_{k=0}^{i-1} l_{ik}u_{kj} \label{eq.u.pivot} \end{equation}\] \[\begin{equation} l_{ij}=\frac{1}{u_{jj}} \left(a_{ij}^{(N-1)}-\sum_{k=0}^{j-1} l_{ik}u_{kj}\right) \label{eq.l.pivot} \end{equation}\] となる。
Since \(u_{jj}\) appears in the denominator of eq. (\ref{eq.l}), a small value of \(u_{jj}\) results in a large computational error. To suppress this problem, swappings of the rows of \(\myvector{A}\) (partial pivot selection) are applied in the LU decomposition. Basically, the row to be swapped is determined to maximize \(|u_{jj}|\); more detail on the selection of the row is described in the next section. Below, let \(\myvector{A^{(m)}}\) be the matrix obtained by \(m\)-times swappings of the rows from the original matrix \(\myvector{A}\), and \(a_{ij}^{(m)}\) be its \((i,j)\)th component. This function does not calculate the LU decomposition of \(\myvector{A}\) but calculates the LU decomposition of \(\myvector{A^{(N-1)}}\). Therefore strictly the equation to be solved is not (\ref{eq.decomposition}) but (\ref{eq.decomposition.pivot}). The solution is thus not given by eqs. (\ref{eq.u}) and (\ref{eq.l}) but given by (\ref{eq.u.pivot}) and (\ref{eq.l.pivot}).

5. ピボット選択において交換する行の選択 (Choice of the row to swap in the pivot selection)

行列\(\myvector{A^{(m+1)}}\)は \(\myvector{A^{(m)}}\)の第\(m\)行とそれよりも後の行との 交換によって生成する (但し\(\myvector{A^{(0)}}=\myvector{A}\)とする)。 交換する行は以下のように選ぶ。まず各\(i(\geq m)\)について \[\begin{equation} c_i^{(m)} \equiv a_{im}^{(m)}-\sum_{k=0}^{m-1} l_{ik}u_{km} \label{eq.c} \end{equation}\] を計算する。\(|c_i^{(m)}|\)を\(\myvector{A^{(m)}}\)の第\(i\)行の絶対値最大値で 規格化した値\(|\hat{c}_i^{(m)}|\)を求め、その最大値を与える行\(I\)を \(\myvector{A^{(m)}}\)の第\(m\)行と交換する。
The matrix \(\myvector{A^{(m+1)}}\) is generated by swapping the \(m\)th row and a later row of \(\myvector{A^{(m)}}\); here \(\myvector{A^{(0)}}=\myvector{A}\). The row to be swapped is selected as follows. First, \(c_i^{(m)}\) defined by eq. (\ref{eq.c}) is calculated for each \(i (\geq m)\). Then the values \(|\hat{c}_i^{(m)}|\), defined as \(|c_i^{(m)}|\) normalized by the absolute maximum of the \(i\)th row of \(\myvector{A^{(m)}}\), are calculated. The row \(I\) that gives the maximum of \(|c_i^{(m)}|\) is swapped with the \(m\)th row of \(\myvector{A^{(m)}}\).
この選び方では基本的に\(|u_{mm}|\)が最大になるように 交換する行を選んでいることになる。 このことは次のようにして理解できる。 \(c_i^{(m)}\)は「もしも\(\myvector{A^{(m)}}\)の第\(i\)行を 第\(m\)行と交換したとしたら\(u_{mm}\)になるはずの値」である。 というのも仮に\(\myvector{A^{(m)}}\)の第\(i\)行を第\(m\)行と交換したとしたら \[\begin{equation} a_{im}^{(m)}=a_{mm}^{(m+1)}=a_{mm}^{(N-1)} \label{eq.a.swap_im} \end{equation}\] が成り立つ (1つ目の等号は\(\myvector{A^{(m)}}\)の第\(i\)行と第\(m\)行を交換したものが \(\myvector{A^{(m+1)}}\)であるため、 2つ目の等号はそれ以降に第\(m\)行と他の行の交換を行わないため)。 (\ref{eq.a.swap_im})を(\ref{eq.c})の右辺に代入すると (\ref{eq.u.pivot})で\(i=j=m\)とした式の右辺と一致する。 以上により\(c_i^{(m)}\)は「もしも\(\myvector{A^{(m)}}\)の第\(i\)行を 第\(m\)行と交換したとしたら\(u_{mm}\)になるはずの値」であることが示され、 このことから\(|c_i^{(m)}|\)の最大値を与える\(i\)を選べば \(|u_{mm}|\)が最大化されるというわけである。
By this method, basically the row that maximizes \(|u_{mm}|\) is chosen, as is explained below. The quantity \(c_i^{(m)}\) represents “the value of \(u_{mm}\) which would be obtained if the \(i\)th and \(m\)th rows of \(\myvector{A^{(m)}}\) were swapped”. This is because if the \(i\)th and \(m\)th rows of \(\myvector{A^{(m)}}\) were swapped, eq. (\ref{eq.a.swap_im}) would hold; the first equation holds because \(\myvector{A^{(m+1)}}\) is obtained by swapping the \(i\)th and \(m\)th rows of \(\myvector{A^{(m)}}\), and the second does because the \(m\)th row will no more be swapped with other rows. Inserting eq. (\ref{eq.a.swap_im}) into the right hand side of (\ref{eq.c}) results in the right hand side of (\ref{eq.u.pivot}) with \(i=j=m\). Therefore it is verified that the quantity \(c_i^{(m)}\) represents “the value of \(u_{mm}\) which would be obtained if the \(i\)th and \(m\)th rows of \(\myvector{A^{(m)}}\) were swapped”; thus choosing \(i\) that maximizes \(|c_i^{(m)}|\) results in maximization of \(|u_{mm}|\).
\(|c_i^{(m)}|\)ではなく\(|\hat{c}_i^{(m)}|\)の最大値を与える\(i\)を選んでいるのは LU分解が連立方程式の解の計算にしばしば用いられることを踏まえたものである。 連立方程式においてはある行を丸ごと定数倍しても解は変わらない。 そのような定数倍の因子によってピボットが変わってしまうのは 好ましくないと考えられるため、 各行の絶対値最大値で規格化した上で大小を評価するようにしている。
The \(i\) that does not maximizes \(|c_i^{(m)}|\) but does maximizes \(|\hat{c}_i^{(m)}|\) is chosen because the LU decomposition is often used for solving simultaneous equations. The solution of a simultaneous equation does not change by multiplying a constant factor with a line. Therefore the pivod selection should not depend on the constant factor, so that each line is evaluated after normalization by the absolute maximum of the line.

◆実際の計算手順(Calculation steps)

1. 計算の仕組み (Outline of the calculation procedure)

もしも\(\myvector{A}\)の行の交換だけ先に済ませてしまって \(\myvector{A^{(N-1)}}\)を得てから (\ref{eq.u.pivot})(\ref{eq.l.pivot})式を用いて \(\myvector{L}\),\(\myvector{U}\)を求めることができるなら アルゴリズムは単純で分かり易いものになる。 しかし実際は\(\myvector{A}\)の行の交換と\(l_{ij}\),\(u_{ij}\)の導出を 同時進行で行わなければならず、 行を交換する途中段階の(まだ確定していない) \(\myvector{A}\)を用いて\(l_{ij}\),\(u_{ij}\)を一度求めたあとで \(\myvector{A}\)の行を交換するというややこしい手順になる。
If the rows of \(\myvector{A}\) can be swapped first to obtain \(\myvector{A^{(N-1)}}\), before applying eqs. (\ref{eq.u.pivot}) and (\ref{eq.l.pivot}) to obtain \(\myvector{L}\) and \(\myvector{U}\), then the algorithm would be simple and comprehensive. Indeed, the swappings of the rows of \(\myvector{A}\) and the calculations of \(l_{ij}\) and \(u_{ij}\) must be done simultaneously. This means that \(l_{ij}\) and \(u_{ij}\) must be calculated from \(\myvector{A}\) at an intermediate state (i.e., the row swappings did not complete), after then rows of \(\myvector{A}\) are swapped. This is the difficult aspect of the LU decomposition algorithm.
計算手順を理解する上でポイントとなるのは \(\myvector{A^{(m)}}\)の第\(m\)行を それよりも後の行と交換することによって \(\myvector{A^{(m+1)}}\)を求めるという点である。 これを\(m\)を増加させながら行うのであるから、 \(\myvector{A^{(m)}}\)が求まった時点で第\(m-1\)行までは確定している。 つまり\(\myvector{A^{(m)}}\)と\(\myvector{A^{(N-1)}}\)は 第\(m-1\)行までは同じである。 この性質を利用することで\(\myvector{L}\)と\(\myvector{U}\)の成分が 順々に求まる仕組みになっている。
The important aspect for understanding the algorithm is that \(\myvector{A^{(m+1)}}\) is obtained by swapping the \(m\)th and later lines of \(\myvector{A^{(m)}}\). This operation is repeated for increasing \(m\), meaning that the (\(m-1\))th and previous rows are fixed when \(\myvector{A^{(m)}}\) is obtained. That is, \(\myvector{A^{(m)}}\) and \(\myvector{A^{(N-1)}}\) have the same values in (\(m-1\))th and previous rows. The algorithm makes use of this nature to sequentially calculate the components of \(\myvector{L}\) and \(\myvector{U}\).
\(\myvector{U}\)については上記の性質を用いることで 各成分の値を1発で求めることができる。 一方、\(\myvector{L}\)については値の設定後に行の交換が必要になる。 \(\myvector{A^{(m)}}\)を用いてひとまず仮の値を与えておき、 \(\myvector{A^{(m+1)}}\)の計算に合わせて\(\myvector{L}\)の行の交換も行う。 そこで、\(\myvector{A^{(m)}}\)から計算した 仮の\(\myvector{L}\)を\(\myvector{L^{(m)}}\)、 その\((i,j)\)成分を\(l_{ij}^{(m)}\)と表す。 \(\myvector{L^{(m)}}\)の全成分の値が求まるわけではないが、便宜上このように表記する。 \(\myvector{L^{(N-1)}}\)が最終的な\(\myvector{L}\)である。 後述の通り、\(\myvector{L^{(m)}}\)の第\(m\)行とそれ以降の行の 第\(m-1\)列までの成分を入れ替えたものが\(\myvector{L^{(m+1)}}\)となる。 したがって\(\myvector{L^{(m)}}\)についても\(\myvector{A^{(m)}}\)と同様に 第\(m-1\)行までは確定している。 すなわち\(\myvector{L^{(m)}}\)と\(\myvector{L^{(N-1)}}\)は 第\(m-1\)行までは同じである。
This nature enables a direct calculation of \(\myvector{U}\), whereas it is needed to swap the rows of \(\myvector{L}\) after setting the values. Temporal values of the components of \(\myvector{L}\) are set using \(\myvector{A^{(m)}}\), and then the rows of \(\myvector{L}\) are swapped corresponding to the calculation of \(\myvector{A^{(m+1)}}\). Hereafter, the temporal \(\myvector{L}\) determined by \(\myvector{A^{(m)}}\) is denoted by \(\myvector{L^{(m)}}\), and its \((i,j)\) component is denoted by \(l_{ij}^{(m)}\), although not all the components of \(\myvector{L^{(m)}}\) are determined. The matrix \(\myvector{L^{(N-1)}}\) is the final \(\myvector{L}\). As is explained later, \(\myvector{L^{(m+1)}}\) is obtained by swapping \(m\)th and later rows of \(\myvector{L^{(m)}}\) within (\(m-1\))th and previous columns. Therefore, same as \(\myvector{A^{(m)}}\), the (\(m-1\))th and previous rows are fixed when \(\myvector{L^{(m)}}\) is obtained. That is, \(\myvector{L^{(m)}}\) and \(\myvector{L^{(N-1)}}\) have the same values in (\(m-1\))th and previous rows.

2. 計算の手順 (The calculation procedure)

\(m=0,\cdots,N-1\)について以下の処理を行う。
The following procedures are repeated for \(m=0,\cdots,N-1\).
  1. \(\myvector{U}\)の第\(m\)列の第\(m-1\)行までの成分の計算。 \(i=0,\cdots,m-1\)について、 \[\begin{equation} u_{im}=a_{im}^{(m)}-\sum_{k=0}^{i-1} l_{ik}^{(m)}u_{km} \label{eq.uim} \end{equation}\] を計算する。
    Calculate the components of \(\myvector{U}\) for (\(m-1\))th and previous rows of the \(m\)th column, by applying eq. (\ref{eq.uim}) for \(i=0,\cdots,m-1\).
  2. ピボット選択。 \(\myvector{A^{(m+1)}}\)を作成するために \(\myvector{A^{(m)}}\)の第\(m\)行と交換する行\(I\)を以下の手順で決定する。 まず\(i=m,\cdots,N-1\)について \[\begin{equation} c_i^{(m)} \equiv a_{im}^{(m)}-\sum_{k=0}^{m-1}l_{ik}^{(m)}u_{km} \label{eq.c.use} \end{equation}\] を計算する。 これを\(\myvector{A^{(m)}}\)の第\(i\)行の絶対値最大値で割ることによって \(\hat{c}_i^{(m)}\)を求める。 \(|\hat{c}_i^{(m)}|\)の最大値を与える\(i\)が求める\(I\)である。
    The pivot selection. The row \(I\) to be replaced with the \(m\)th row of \(\myvector{A^{(m)}}\) to create \(\myvector{A^{(m+1)}}\) is determined by the following steps. First, the value \(c_i^{(m)}\) is calculated by eq. (\ref{eq.c.use}) for each \(i=m,\cdots,N-1\). Then this value for each \(i\) is divided by the absolute maximum of the \(i\)th row of \(\myvector{A^{(m)}}\) to obtain \(\hat{c}_i^{(m)}\). The optimal \(i\) that gives the maximum \(|\hat{c}_i^{(m)}|\) is the row \(I\) needed here.
  3. \(\myvector{A^{(m)}}\), \(\myvector{L^{(m)}}\)の 第\(m\)行と第\(I\)行を交換することによって \(\myvector{A^{(m+1)}}\), \(\myvector{L^{(m+1)}}\)を作成する。 但し\(\myvector{L^{(m+1)}}\)については第\(m-1\)列までの範囲を対象とする (第\(m\)列以降はまだ求まっていない)。 式で書くと \[\begin{equation} a_{ij}^{(m+1)}= \begin{cases} a_{Ij}^{(m)} & (i=m) \\ a_{mj}^{(m)} & (i=I) \\ a_{ij}^{(m)} & (i\neq m,i\neq I) \end{cases} \label{eq.a.update} \end{equation}\] \[\begin{equation} l_{ij}^{(m+1)}= \begin{cases} l_{Ij}^{(m)} & (i=m,j<m) \\ l_{mj}^{(m)} & (i=I,j<m) \\ l_{ij}^{(m)} & (i\neq m,i\neq I,j<m) \end{cases} \label{eq.l.update} \end{equation}\] である。
    Matrices \(\myvector{A^{(m+1)}}\) and \(\myvector{L^{(m+1)}}\) are created by swapping the \(m\)th and \(I\)th rows of \(\myvector{A^{(m)}}\) and \(\myvector{L^{(m)}}\), respectively. Here, only the (\(m-1\))th and previous columns are used for \(\myvector{L^{(m+1)}}\) because \(m\)th and later columns have not been obtained. This operation is written by eqs. (\ref{eq.a.update}) and (\ref{eq.l.update}).
  4. \(\myvector{U}\)の\((m,m)\)成分の計算。 \[\begin{equation} u_{mm}=a_{mm}^{(m+1)}-\sum_{k=0}^{m-1} l_{mk}^{(m+1)}u_{km} \label{eq.umm} \end{equation}\] を計算する。
    Calculate the \((m,m)\) component of \(\myvector{U}\) by applying eq. (\ref{eq.umm}).
  5. \(\myvector{L^{(m+1)}}\)の第\(m\)列の第\(m+1\)行以降の成分の計算。 \(i=m+1,\cdots,N-1\)について、 \[\begin{equation} l_{im}^{(m+1)} \equiv \frac{1}{u_{mm}} \left(a_{im}^{(m+1)}-\sum_{k=0}^{m-1} l_{ik}^{(m+1)}u_{km}\right) \label{eq.lim} \end{equation}\] を計算する。
    Calculate the components of \(\myvector{L^{(m+1)}}\) for (\(m+1\))th and later rows of the \(m\)th column, by applying eq. (\ref{eq.lim}).


3. 計算手順の理解 (Understanding the calculation procedure)

上の手順でLU分解を実行できることは次のようにして理解できる。
Below, we explain why the LU decomposition can be conducted by the procedure above.

3-1. \(m=0\)での動作 (Procedures for \(m=0\))
\(m=0\)のときはステップ1は(\(i=0,\cdots,m-1\)についての処理なので)スキップされる。 またステップ2以降でも 計算式の右辺の\(k\)に関する和の項が(\(k=m-1\)までの和なので)登場しない。 したがって\(\myvector{L}\)や\(\myvector{U}\)の成分が 全く求まっていない状態から始めて 上記の計算手順を実行できる。
When \(m=1\), the step 1 is skipped (because this step is conducted for \(i=0,\cdots,m-1\)). After the step 2, the summation terms with respect to \(k\) in the right hand sides of the formulas do not appear (because the summations are taken until \(k=m-1\)). Therefore the procedures can be conducted with no value given to the components of \(\myvector{L}\) and \(\myvector{U}\).

3-2. 計算式の右辺には計算済みの量のみが登場することの確認 (Velification that the right hand sides of the formulas do not have unknowns)
ステップ1-5において求める行列成分と計算に必要な行列成分を整理すると 表1のようになる。
Table 1 shows the components of matrices to be calculated and those needed in individual steps of 1-5.
表1. 計算の各ステップで求める行列成分と必要な行列成分
Table 1. Components of matrices computed and needed in individual steps
ステップ
Step 		求める成分(\(\myvector{L}\))
Components computed (\(\myvector{L}\)) 		求める成分(\(\myvector{U}\))
Components computed (\(\myvector{U}\)) 		必要な成分(\(\myvector{L}\))
Components needed (\(\myvector{L}\)) 		必要な成分(\(\myvector{U}\))
Components needed (\(\myvector{U}\))
1   行(Row): \(<m\)
列(Column): \(=m\)		 行(Row): \(<m\)
列(Column): \(<m\) 1		 行(Row): \(<m\) 2
列(Column): \(=m\)
2     行(Row): \(\geq m\)
列(Column): \(<m\)		 行(Row): \(<m\)
列(Column): \(=m\)
4   行(Row): \(=m\)
列(Column): \(=m\)		 行(Row): \(=m\)
列(Column): \(<m\)		 行(Row): \(<m\)
列(Column): \(=m\)
5 行(Row): \(>m\)
列(Column): \(=m\)		   行(Row): \(>m\)
列(Column): \(<m\)		 行(Row): \(\leq m\)
列(Column): \(=m\)
  1. 但し行番号\(>\)列番号の成分のみ
    Only the components with the row index greater than the column index
  2. 但し求める成分よりも前の行のみ
    Only the rows before that for computation

まず\(\myvector{U}\)の成分がどのように計算式に登場するかを見てみる。 表1より、いずれのステップにおいても \(\myvector{U}\)の第\(m\)列のみが計算式の右辺に登場することが分かる。 そして行番号範囲としては、 ステップ1では求める成分よりも前の行のみ、 ステップ2と4ではステップ1で計算済みの第\(m-1\)行目までの成分のみ、 ステップ5では(\ref{eq.lim})ステップ1,4で計算済みの 第\(m\)行目までの成分のみが計算式の右辺に登場する。 したがってどのステップでもその時点で既に計算済みの\(\myvector{U}\)の成分のみが 用いられていることが分かる。
First, let us focus on the components of \(\myvector{U}\) that appear in the formulas. According to Table 1, only the \(m\)th column of \(\myvector{U}\) appears in the right hand sides of the formulas. The index ranges of the rows that appear in the right hand sides of the formulas are: those prior to the component to be calculated in the step 1, (\(m-1\))th and previous rows (which have already been calculated in the step 1) in the steps 2 and 4, and \(m\)th and previous rows (which have already been calculated in the steps 1 and 4) in the step 5. Therefore only the components of \(\myvector{U}\) that have already been calculated appear in each step of the calculation.
次に\(\myvector{L}\)について見てみる。 表1より、いずれのステップにおいても 行番号\(>\)列番号かつ列番号\(<m\)の成分のみが登場することが分かる。 ステップ5において第\(m\)列の第\(m+1\)行以降の全成分を計算するのであり、 これを\(m=0\)から始めて\(m\)の値を増やしながら繰り返すのであるから、 ある\(m\)においてステップ1が始まる以前の段階で 第\(m-1\)列目までの行番号\(>\)列番号の成分は全て計算済みである。 したがってステップ1-5においては 既に計算済みの\(\myvector{L}\)の成分のみが用いられていることが分かる。
Next, let us focus on \(\myvector{L}\). According to Table 1, only the column indices less than the row indices and less than \(m\) appear in each step. All the (\(m+1\))th and later row components in the \(m\)th column are computed in the step 5, which is repeated with increading \(m\) starting from \(m=0\), so that all the components in the (\(m-1\))th and previous columns having the row indices greater than the column indices have been calculated prior to the step 1 of each \(m\). Therefore only the components of \(\myvector{L}\) that have already been calculated appear in each step of the calculation.

3-3. 行の交換が計算結果に及ぼす影響範囲 (Extent of the effects of row swappings on the computation results)
ステップ3において\(\myvector{A}\)と\(\myvector{L}\)の第\(m\)行を それ以降の行と交換している。 行の交換が行われるのはこのステップのみである。 したがって各計算式の右辺に登場する\(\myvector{U}\)の成分は その後の行の入れ替え操作によって影響を受けない最終的な値である。 またステップ1,2における\(\myvector{A}\)と\(\myvector{L}\)の第\(m-1\)行までの成分、 ステップ4,5における\(\myvector{A}\)と\(\myvector{L}\)の第\(m\)行までの成分も その後の行の入れ替え操作によって影響を受けない最終的な値である。 以降の解説ではこのことを前提とする。
In the step 3, the \(m\)th and later rows are swapped for \(\myvector{A}\) and \(\myvector{L}\). The row swappings appear only in this step. Therefore the components of \(\myvector{U}\) in the right hand side of each equation are the final values which will not change by later row swapping operations. Also, the components of \(\myvector{A}\) and \(\myvector{L}\) in the (\(m-1\))th and previous rows at steps 1 and 2, as well as those in the \(m\)th and previous rows at steps 4 and 5, are the final values which will not change by later row swapping operations. This nature is fully used in the description below.

3-4. ステップ1における\(\myvector{U}\)の計算式(\ref{eq.uim})が 本来の計算式(\ref{eq.u.pivot})と等価であることの確認 (Confirmation of the equivalence between eq. (\ref{eq.uim}) in the step 1 and the original formula (\ref{eq.u.pivot}) for the computation of \(\myvector{U}\))
(\ref{eq.uim})式の右辺に登場するのは \(\myvector{A}\)と\(\myvector{L}\)の第\(i\)行成分のみである。 \(i=0,\cdots,m-1\)について計算するのであるから、 右辺に登場する量はすべて その後の行の入れ替え操作によって影響を受けない最終的な値であることが分かる。 すなわち(\ref{eq.uim})式において \(a_{im}^{(m)}=a_{im}^{(N-1)}\), \(l_{ik}^{(m)}=l_{ik}^{(N-1)}=l_{ik}\)であり、 ステップ1で用いられる(\ref{eq.uim})式は 本来の\(u_{ij}\)の計算式(\ref{eq.u.pivot})で\(j=m\)とした式になっている。
Only the \(i\)th row of \(\myvector{A}\) and \(\myvector{L}\) appears in the right hand side of eq. (\ref{eq.uim}). Since this calculation is conducted for \(i=0,\cdots,m-1\), all the quantities in the right hand side of this equation are the final values which will not change by later row swapping operations. Namely, relations \(a_{im}^{(m)}=a_{im}^{(N-1)}\) and \(l_{ik}^{(m)}=l_{ik}^{(N-1)}=l_{ik}\) hold in eq. (\ref{eq.uim}), so that this equation is equivalent to the original formula (\ref{eq.u.pivot}) for \(u_{ij}\) with \(j=m\).

3-5. ステップ2における行の評価式(\ref{eq.c.use})と 元々の式(\ref{eq.c})の関係 (Relation between the row evaluation formula of eq. (\ref{eq.c.use}) in the step 2 and the original formula of eq. (\ref{eq.c}))
ステップ2では\(i=m,\cdots,N-1\)について (\ref{eq.c.use})式を計算するのであるから、 右辺の\(a_{im}^{(m)}\)や\(l_{ik}^{(m)}\)は後から変化しうる値である。 しかし、仮に直後のステップ3で第\(i\)行を第\(m\)行と交換したとしたら それ以上は変化しない。 すなわち\(a_{im}^{(m)}\), \(l_{ik}^{(m)}\)は 「仮に第\(m\)行と第\(i\)行を交換したとしたら 最終的に\(a_{mm}^{(N-1)}\), \(l_{mk}\)となる値」である。 (\ref{eq.c})式では行の交換に伴い \(\myvector{L}\)が変化することを考慮していなかった。 これを考慮に入れた上で行の交換後の\(u_{mm}\)を評価する式が(\ref{eq.c.use})であり、 この式に基づいて計算した\(|\hat{c}_i^{(m)}|\)が最大になるように\(I\)を選べば \(u_{mm}\)が最大化される。
Since eq. (\ref{eq.c.use}) is applied to \(i=m,\cdots,N-1\) in the step 2, \(a_{im}^{(m)}\) and \(l_{ik}^{(m)}\) in the right hand side may change after the calculation. However, if the \(i\)th row is swapped with the \(m\)th row at the immediately after step 3, then the values will be fixed at that stage. Therefore, \(a_{im}^{(m)}\) and \(l_{ik}^{(m)}\) are “the value of final \(a_{mm}^{(N-1)}\) and \(l_{mk}\), respectively, which would be obtained if the \(m\)th and \(i\)th rows were swapped”. In eq. (\ref{eq.c}), a change of \(\myvector{L}\) due to the row swappings was not taken into account. Eq. (\ref{eq.c.use}) is a modified version taking into account this aspect for evaluating \(u_{mm}\) after the row swapping, so that choosing \(I\) that maximizes \(|\hat{c}_i^{(m)}|\) based on this formula would maximize \(u_{mm}\).

3-6. ステップ4における\(\myvector{U}\)の計算式(\ref{eq.umm})が 本来の計算式(\ref{eq.u.pivot})と等価であることの確認 (Confirmation of the equivalence between eq. (\ref{eq.umm}) in the step 4 and the original formula (\ref{eq.u.pivot}) for the computation of \(\myvector{U}\))
(\ref{eq.uim})式の右辺に登場するのは \(\myvector{A}\)と\(\myvector{L}\)の第\(m\)行成分のみであるので その後の行の入れ替え操作によって影響を受けない最終的な値である。 すなわち(\ref{eq.umm})式において \(a_{mm}^{(m+1)}=a_{mm}^{(N-1)}\), \(l_{mk}^{(m+1)}=l_{mk}^{(N-1)}=l_{mk}\)であり、 ステップ4で用いられる(\ref{eq.umm})式は 本来の\(u_{ij}\)の計算式(\ref{eq.u.pivot})で\(i=j=m\)とした式になっている。
Only the \(m\)th row of \(\myvector{A}\) and \(\myvector{L}\) appears in the right hand side of eq. (\ref{eq.umm}). Therefore all the quantities in the right hand side of this equation are the final values which will not change by later row swapping operations. Namely, relations \(a_{mm}^{(m+1)}=a_{mm}^{(N-1)}\) and \(l_{mk}^{(m+1)}=l_{mk}^{(N-1)}=l_{mk}\) hold in eq. (\ref{eq.umm}), so that this equation is equivalent to the original formula (\ref{eq.u.pivot}) for \(u_{ij}\) with \(i=j=m\).

3-7. ステップ5における\(\myvector{L}\)の計算式(\ref{eq.lim})が 本来の計算式(\ref{eq.l.pivot})と等価であることの確認 (Confirmation of the equivalence between eq. (\ref{eq.lim}) in the step 5 and the original formula (\ref{eq.l.pivot}) for the computation of \(\myvector{L}\))
(\ref{eq.lim})式では\(\myvector{L}\)の第\(i\)行を計算するために \(\myvector{A}\)と\(\myvector{L}\)の第\(i\)行成分のみを用いている。 したがって、後に(より大きな\(m\)に対する)ステップ3で行を入れ替える際に、 (\ref{eq.lim})式の左辺と右辺に登場する \(\myvector{A}\)や\(\myvector{L}\)の成分間の関係は維持される。 すなわち、計算時点で第\(i\)行だった行が最終的に第\(i’\)行に移動する場合、 \(l_{im}^{(m+1)}=l_{i’m}\), \(l_{ik}^{(m+1)}=l_{i’k}\), \(a_{im}^{(m+1)}=a_{i’m}^{(N-1)}\)であり、 これらを(\ref{eq.lim})に代入すると \[\begin{equation} l_{i’m} =\frac{1}{u_{mm}}\left(a_{i’m}^{(N-1)}-\sum_{k=0}^{m-1} l_{i’k}u_{km}\right) \label{eq.lim.arranged} \end{equation}\] となって(\ref{eq.l.pivot})式で\(i=i’\), \(j=m\)とした式と一致する。
In eq. (\ref{eq.lim}), only the \(i\)th row components of \(\myvector{A}\) and \(\myvector{L}\) are used to compute the \(i\)th row of \(\myvector{L}\). Therefore the relation between the left and right hand sides of eq. (\ref{eq.lim}) is kept by the later row swappings at the step 3 for greater \(m\). Namely, if the \(i\)th row at the computation will finally move to the \(i’\)th column, then relations \(l_{im}^{(m+1)}=l_{i’m}\), \(l_{ik}^{(m+1)}=l_{i’k}\), and \(a_{im}^{(m+1)}=a_{i’m}^{(N-1)}\) hold. Inserting these relations into (\ref{eq.lim}) results in (\ref{eq.lim.arranged}), which is equivalent to eq. (\ref{eq.l.pivot}) with \(i=i’\) and \(j=m\).

4. 非正則行列とゼロ割りの問題の扱い (Treatments of irregular matrices and zero-division problems)

LU分解の計算において割り算は2ヶ所に登場する。 ステップ2(ピボット選択)において\(c_i^{(m)}\)を\(\myvector{A^{(m)}}\)の第\(i\)行の 絶対値最大値で割る操作と、 ステップ5((\ref{eq.lim})式)における\(u_{mm}\)による割り算である。 これらの割り算において分母が0となる場合の扱いについて述べる。
Division operations appear at two places in the calculation of the LU decomposition: division of \(c_i^{(m)}\) by the absolute maximum of the \(i\)th row of \(\myvector{A^{(m)}}\) in the step 2 (pivot selection), and division by \(u_{mm}\) in the step 5 (eq. (\ref{eq.lim})). We describe the treatment in cases where the denominators are zero in these division operations.

4-1. \(c_i^{(m)}\)を\(\myvector{A^{(m)}}\)の第\(i\)行の絶対値最大値で 割る操作(ステップ2) (Division of \(c_i^{(m)}\) by the absolute maximum of the \(i\)th row of \(\myvector{A^{(m)}}\) (step 2))
この操作で分母(\(\myvector{A^{(m)}}\)の第\(i\)行の絶対値最大値)がゼロになるのは \(\myvector{A^{(m)}}\)の第\(i\)行成分が全て0の場合である。 このとき分子(\(c_i^{(m)}\))もまた0になることが次のようにして確かめられる。
The denominator of this operation (the absolute maximum of the \(i\)th row of \(\myvector{A^{(m)}}\)) becomes zero if and only if all the components of the \(i\)th row of \(\myvector{A^{(m)}}\) are zero. In this case, the numerator of this operation (\(c_i^{(m)}\)) is also zero, as is shown below.

すなわち、\(\myvector{A^{(m)}}\)の第\(i\)行成分が全て0の場合には \(\hat{c}_i^{(m)}\)は0/0の不定値となる。
Therefore, in cases where all the \(i\)th row components of \(\myvector{A^{(m)}}\) are zero, \(\hat{c}_i^{(m)}\) becomes 0/0 (indefinite).
実際の処理としてはこの場合には\(\hat{c}_i^{(m)}=0\)とおく。 というのも\(\hat{c}_i^{(m)}\)はピボット選択において 最適な行を見つけるために計算する量であるので 大小関係のみが意味を持ち、値自体は利用されない。 \(|\hat{c}_i^{(m)}|\)が最大になる行が選ばれるのであるから、 \(\hat{c}_i^{(m)}=0\)とおくことでこの行が選ばれるのをなるべく避けられる。
In the program, \(\hat{c}_i^{(m)}=0\) is used in this case. This is because \(\hat{c}_i^{(m)}\) is used only for finding the best row in the pivot selection; only the magnitude correlation is important, and the value of \(\hat{c}_i^{(m)}\) itself is not used. Since the row which gives the maximum of \(|\hat{c}_i^{(m)}|\) is selected in the pivot selection, setting \(\hat{c}_i^{(m)}=0\) results in avoiding this row to be selected as long as possible.

4-2. (\ref{eq.lim})式における\(u_{mm}\)による割り算(ステップ5) (Division by \(u_{mm}\) in eq. (\ref{eq.lim}) (step 5))
(\ref{eq.c.use})-(\ref{eq.lim})式より \[\begin{eqnarray} l_{im}^{(m+1)} &=& \frac{1}{u_{mm}} \left(a_{im}^{(m+1)}-\sum_{k=0}^{m-1} l_{ik}^{(m+1)}u_{km}\right) \nonumber \\ &=& \frac{1}{u_{mm}} \left(a_{im}^{(m)}-\sum_{k=0}^{m-1} l_{ik}^{(m)}u_{km}\right) \nonumber \\ &=& \frac{1}{u_{mm}}c_i^{(m)} \hspace{1em} (i\neq I) \label{eq.lim.use_c} \end{eqnarray}\] \[\begin{eqnarray} l_{Im}^{(m+1)} &=& \frac{1}{u_{mm}} \left(a_{Im}^{(m+1)}-\sum_{k=0}^{m-1} l_{Ik}^{(m+1)}u_{km}\right) \nonumber \\ &=& \frac{1}{u_{mm}} \left(a_{mm}^{(m)}-\sum_{k=0}^{m-1} l_{mk}^{(m)}u_{km}\right) \nonumber \\ &=& \frac{1}{u_{mm}}c_m^{(m)} \label{eq.lIm.use_c} \end{eqnarray}\] \[\begin{eqnarray} u_{mm} &=& a_{mm}^{(m+1)}-\sum_{k=0}^{m-1} l_{mk}^{(m+1)}u_{km} \nonumber \\ &=& a_{Im}^{(m+1)}-\sum_{k=0}^{m-1} l_{Ik}^{(m+1)}u_{km} \nonumber \\ &=& c_I^{(m)} \label{eq.umm.use_c} \end{eqnarray}\] である。 このことと\(|\hat{c}_i^{(m)}|\)の最大値を与える行を第\(I\)行とおいていることから、 \(u_{mm}=0\)のときは全ての\(i\)について \(l_{im}^{(m+1)}\)の計算式の分子が0となることが分かる (\(u_{mm}=0\) \(\rightarrow\) \(c_I^{(m)}=0\) \(\rightarrow\) \(\hat{c}_I^{(m)}=0\) \(\rightarrow\) \(\hat{c}_i^{(m)}=0\) \(\rightarrow\) \(c_i^{(m)}=0\))。 すなわち\(u_{mm}=0\)のときは\(l_{im}^{(m+1)}\)は0/0の不定値となる。
From eqs. (\ref{eq.c.use})-(\ref{eq.lim}), we can derive eqs. (\ref{eq.lim.use_c})-(\ref{eq.umm.use_c}). In addition, remember that the row which gives the maximum of \(|\hat{c}_i^{(m)}|\) was defined as the \(I\)th row. These facts lead to the following conclusion: when \(u_{mm}=0\), the numerator of the formula of \(l_{im}^{(m+1)}\) is zero for every \(i\). This is shown as: \(u_{mm}=0\) \(\rightarrow\) \(c_I^{(m)}=0\) \(\rightarrow\) \(\hat{c}_I^{(m)}=0\) \(\rightarrow\) \(\hat{c}_i^{(m)}=0\) \(\rightarrow\) \(c_i^{(m)}=0\). Therefore, when \(u_{mm}=0\), \(l_{im}^{(m+1)}\) becomes 0/0 (indefinite).
このようなことが起きるのは\(\myvector{A}\)が非正則行列の場合であり、 そもそもLU分解ができないということである。 この場合、プログラムはエラーとして終了する。 しかし計算が止まってしまうと不都合なときもあるので、 マクロFORCE_LU_DECOMPOSITIONの値がTRUEのときは\(u_{mm}=0\)であってもエラーとせずに \(l_{im}^{(m+1)}=0\)とおいて計算を続行する。
This occurs when the matrix \(\myvector{A}\) is a irregular matrix, for which the LU decomposition is impossible. In this case the program finishes as an error. However, you can avoid the termination of the program by setting a macro FORCE_LU_DECOMPOSITION to TRUE; in this case, the program does not issue an error even when \(u_{mm}=0\), and instead the program continues the calculation assuming \(l_{im}^{(m+1)}=0\).

4-3. 0と見なす条件 (Condition to regard a quantity as zero)
\(\myvector{A^{(m)}}\)の第\(i\)行成分の絶対値最大値が0であると見なすのは 本当に0.0の場合のみである。 というのも、LU分解を連立方程式の解の計算に用いる場合を考えると \(\myvector{A^{(m)}}\)の各行には定数倍の不確定があり、 値そのものが大きいか小さいかを行をまたいで比較しても意味がない。 またLU分解では\(\myvector{A}\)の行の入れ替えのみで 足し算や引き算等は行わないので、 丸め誤差によって本当は0なのに見かけ上ノンゼロになるという心配は要らない。
The absolute maximum of the \(i\)th row of \(\myvector{A^{(m)}}\) is regarded to be zero only when the value is exactly 0.0. This is because an ambiguity of a constant factor exists in each row of \(\myvector{A^{(m)}}\) when the LU decomposition is applied to solving a simultaneous equation, so that comparison of the values among different rows makes no sense. In addition, only the row swappings of \(\myvector{A}\) appear in the LU decomposition, with no addition and subtraction operations, so that an apparent non-zero value of an indeed zero quantity due to rounding errors does not occur.
\(u_{mm}=0\)と見なすのは \(|u_{mm}|\)がそれまでに出くわした \(|u_{m’m’}|\)(\(m’<m\))の最大値の ZERO_THRESHOLD倍よりも小さい場合とする。 但し\(u_{00}\)に限り、本当に0.0の場合のみ0であると見なす。
The value of \(u_{mm}\) is regarded to be zero when \(|u_{mm}|\) is less than ZERO_THRESHOLD times the maximum of \(|u_{m’m’}|\), where \(m’<m\), for \(m\geq 1\). The value of \(u_{00}\) is regarded to be zero when it is exactly 0.0.

5. 計算式とプログラムでの変数の関係 (Relation between the formulas and variables in the program)

計算式中の変数とプログラム中で用いている変数との対応関係を 表2に示す。 以下の2点を理解すると分かりやすい。
Relation between the variables in the formulas and those in the program code are shown in Table 2. Understanding the following two points would help reading the code.

表2. 計算式とプログラムの変数の対応関係
Table 2. Relation between the variables in the formulas and program code
計算式
Formula 		プログラム
Program code
\(\myvector{A}\)A
\(\myvector{L}\)re.L
\(\myvector{U}\)re.U
\(c_i^{(m)}\)c.main[i]
\(\hat{c}_i^{(m)}\)c_normalized.main[i]
\(a_{ij}^{(m)}\)A.main[re.pivot[i]][j]
\(u_{ij}^{(m)}\)re.U.main[i][j]
\(l_{ij}^{(m)}\)re.L.main[i][j]


◆検証(Validation)

プログラムの動作チェックのため、以下のテストを行った。 まず行列\(\myvector{A}\)を以下のように与えた。 \[\begin{equation} \myvector{A}= \begin{pmatrix} 24 & 27 & 35 & 12 & 14 \\ -15 & -25 & 13 & -26 & -22 \\ -18 & 16 & -31 & -23 & 21 \\ 28 & 11 & 17 & 33 & 20 \\ -29 & -34 & -19 & 30 & 32 \end{pmatrix} \label{eq.A.test} \end{equation}\] この行列を関数LU_decompositionに渡すことにより、 以下のLU分解が得られた。 \[\begin{equation} \myvector{L}= \begin{pmatrix} 1 & 0 & 0 & 0 & 0 \\ 0.62069 & 1 & 0 & 0 & 0 \\ 0.517241 & -0.199814 & 1 & 0 & 0 \\ -0.827586 & -0.0306691 & 0.984045 & 1 & 0 \\ -0.965517 & -0.58829 & -0.665835 & 0.0508279 & 1 \end{pmatrix} \label{eq.L.test_result} \end{equation}\] \[\begin{equation} \myvector{U}= \begin{pmatrix} -29 & -34 & -19 & 30 & 32 \\ 0 & 37.1034 & -19.2069 & -41.6207 & 1.13793 \\ 0 & 0 & 18.9898 & -49.8336 & -38.3243 \\ 0 & 0 & 0 & 84.5897 & 78.2306 \\ 0 & 0 & 0 & 0 & 22.072 \end{pmatrix} \label{eq.U.test_result} \end{equation}\] またpivot選択は以下のようになった。
したがって、LU分解が正しく行われていれば 上記の\(\myvector{L}\)と\(\myvector{U}\)を掛けると 以下の行列が得られるはずである。
\[\begin{equation} \myvector{A^{(N-1)}}= \begin{pmatrix} -29 & -34 & -19 & 30 & 32 \\ -18 & 16 & -31 & -23 & 21 \\ -15 & -25 & 13 & -26 & -22 \\ 24 & 27 & 35 & 12 & 14 \\ 28 & 11 & 17 & 33 & 20 \end{pmatrix} \label{eq.AN1.test_expected_result} \end{equation}\] 実際に(\ref{eq.L.test_result})(\ref{eq.U.test_result})式の \(\myvector{L}\)と\(\myvector{U}\)を 関数multiply_matrix (matrix/operation.h) に渡して積を計算したところ、 (\ref{eq.AN1.test_expected_result})式の\(\myvector{A^{(N-1)}}\)が得られた。 したがってLU分解は正しく動いたと判断できる。
To examine the program, we conducted the following test. First, a matrix \(\myvector{A}\) was given by eq. (\ref{eq.A.test}). The result of the LU decomposition of this matrix, calculated by function LU_decomposition, was those shown by eqs. (\ref{eq.L.test_result}) and (\ref{eq.U.test_result}). The pivot selection was conducted as follows:
Thus, if the LU decomposition was correct, then multiplying \(\myvector{L}\) and \(\myvector{U}\) would result in \(\myvector{A^{(N-1)}}\) shown by eq. (\ref{eq.AN1.test_expected_result}). Indeed, the product of \(\myvector{L}\) and \(\myvector{U}\), calculated by a function multiply_matrix (matrix/operation.h), was consistent with eq. (\ref{eq.AN1.test_expected_result}). Therefore it is concluded that the LU decomposition was conducted correctly.