関数solution_LU_givenLUdecomposition マニュアル

(The documentation of function solution_LU_givenLUdecomposition)

Last Update: 2021/12/6

◆機能・用途(Purpose)

LU分解を用いて連立1次方程式 \[\begin{equation} \myvector{A}\myvector{x}=\myvector{b} \label{eq.problem} \end{equation}\] (\(\myvector{A}\):既知の正方行列、\(\myvector{b}\):既知の列ベクトル、 \(\myvector{x}\):未知の列ベクトル)を解く。係数行列\(\myvector{A}\)そのものではなく \(\myvector{A}\)のLU分解を与えるバージョン。
Solve a simultaneous linear equation (\ref{eq.problem}) using an LU decomposition, where where \(\myvector{A}\) is a known square matrix, \(\myvector{b}\) is a known column vector, and \(\myvector{x}\) is an unknown column vector. In this function, the LU decomposition of the coefficient matrix \(\myvector{A}\) is given instead of \(\myvector{A}\) itself.

◆形式(Format)

#include <matrix/inverse.h>
inline struct columnvector solution_LU_givenLUdecomposition
(const struct LU AA_LU,const struct columnvector b)

◆引数(Arguments)

AA_LU

連立方程式の左辺の係数行列\(\myvector{A}\)のLU分解。 \(\myvector{A}\)を表す行列を関数LU_decomposition (matrix/inverse.h)に渡せば計算できる。厳密には\(\myvector{A}\)自体のLU分解ではなく、 \(\myvector{A}\)の行を交換した行列\(\myvector{A’}\) のLU分解である。
The LU decomposition of the coefficient matrix \(\myvector{A}\) in the left hand side of the simultaneous equation, which can be generated by passing a matrix representing \(\myvector{A}\) into function LU_decomposition (matrix/inverse.h). Strictly, the function LU_decomposition does not yield the LU decomposition of \(\myvector{A}\) itself; instead, the LU decomposition of \(\myvector{A’}\) in which several rows of \(\myvector{A}\) are swapped is given.

連立方程式の右辺の定数ベクトル (\ref{eq.problem})式の\(\myvector{b}\)。
A vector composed of the constants in the right hand side of the simultaneous equations; the column vector \(\myvector{b}\) of eq. (\ref{eq.problem}).

◆戻り値(Return value)

方程式(\ref{eq.problem})式の解\(\myvector{x}\)を表す列ベクトル。
A column vector representing the solution \(\myvector{x}\) of eq. (\ref{eq.problem}).

◆使用例(Example)

struct matrix A;
struct columnvector b;
struct LU AA_LU=LU_decomposition(A);
struct columnvector x=solution_LU_givenLUdecomposition(AA_LU,b);

◆計算式(Formula)

連立方程式(\ref{eq.problem})をLU分解を用いて解くことを考える。 LU分解ではピボット選択を行う関係上、 \(\myvector{A}\)そのもののLU分解はできない。得られるのは\(\myvector{A}\)の行を交換した行列\(\myvector{A’}\)のLU分解 \[\begin{equation} \myvector{A’}=\myvector{L}\myvector{U} \label{eq.LUdecomposition} \end{equation}\] である(ここで\(\myvector{L}\)は下三角行列、\(\myvector{U}\)は上三角行列)。そこで、\(\myvector{A}\)を\(\myvector{A’}\)に変換する際に行ったのと同一の行の交換を \(\myvector{b}\)に対して行って得られる列ベクトルを \(\myvector{b’}\)としよう。このとき、 \[\begin{equation} \myvector{A’}\myvector{x}=\myvector{b’} \label{eq.problem.corrected} \end{equation}\] が成り立つ。(\ref{eq.problem.corrected})に (\ref{eq.LUdecomposition})を代入すると \[\begin{equation} \myvector{L}\myvector{U}\myvector{x}=\myvector{b’} \label{eq.problem.corrected.LU} \end{equation}\] となるが、これは2つの連立方程式 \[\begin{equation} \myvector{L}\myvector{y}=\myvector{b’} \label{eq.problem.L} \end{equation}\] \[\begin{equation} \myvector{U}\myvector{x}=\myvector{y} \label{eq.problem.U} \end{equation}\] に分割できる。まず(\ref{eq.problem.L})式を前進代入で解いて\(\myvector{y}\)を求め、次にそれを用いて(\ref{eq.problem.U})式を後退代入で解けば \(\myvector{x}\)が得られる。
Let us consider to solve eq. (\ref{eq.problem}) using an LU decomposition. Due to a pivot selection used during the LU decomposition, the result of the decomposition is not the LU decomposition of \(\myvector{A}\) but the decomposition of \(\myvector{A’}\) (eq. \ref{eq.LUdecomposition}), where \(\myvector{A’}\) is a matrix obtained by swapping the rows of \(\myvector{A}\), and \(\myvector{L}\) and \(\myvector{U}\) are lower and upper triangular matrices, respectively. Now, let \(\myvector{b’}\) be a column vector obtained by swapping the rows of \(\myvector{b}\), where the swapping is same as that used to constitute \(\myvector{A’}\) from \(\myvector{A}\). Then eq. (\ref{eq.problem.corrected}) holds. Substituting this equation into (\ref{eq.LUdecomposition}) results in (\ref{eq.problem.corrected.LU}), which can be divided to (\ref{eq.problem.L}) and (\ref{eq.problem.U}). The solution \(\myvector{x}\) is obtained by first solving for \(\myvector{y}\) using eq. (\ref{eq.problem.L}) with a forward substitution, and then solving for \(\myvector{x}\) using eq. (\ref{eq.problem.U}) with a backward substitution.

ところでLU分解には丸め誤差が蓄積し易いという難点がある。これを解決するため、上記の方法で暫定的な解を求めた後で以下の方法で解の改良が行われる。 (\ref{eq.problem.corrected})式の真の解を\(\myvector{x_{true}}\)とし、上記の方法で得られた暫定的な数値解を\(\myvector{x_{numerical}}\)、その誤差を \[\begin{equation} \myvector{\delta x} \equiv \myvector{x_{numerical}}-\myvector{x_{true}} \label{eq.delta_x.def} \end{equation}\] とする。このとき \[\begin{eqnarray} \myvector{A’} \myvector{\delta x} &=& \myvector{A’} (\myvector{x_{numerical}}-\myvector{x_{true}}) \nonumber \\ &=& \myvector{A’}\myvector{x_{numerical}} - \myvector{A’}\myvector{x_{true}} \nonumber \\ &=& \myvector{A’}\myvector{x_{numerical}} - \myvector{b’} \label{eq.improvement} \end{eqnarray}\] が成り立つ。 (\ref{eq.improvement})式は2本の方程式 \[\begin{equation} \myvector{L} \myvector{\delta y} = \myvector{A’}\myvector{x_{numerical}} - \myvector{b’} \label{eq.improvement.L} \end{equation}\] \[\begin{equation} \myvector{U} \myvector{\delta x} = \myvector{\delta y} \label{eq.improvement.U} \end{equation}\] に分解できる。 (\ref{eq.improvement.L})を前進代入で解いて\(\myvector{\delta y}\)を求め、その結果を用いて(\ref{eq.improvement.U})を後退代入で解けば \(\myvector{\delta x}\)が求まる。このようにして求めた\(\myvector{\delta x}\)と最初の計算で得られた暫定解\(\myvector{x_{numerical}}\)を (\ref{eq.delta_x.def})に代入すれば \(\myvector{x_{true}}\)を推定できる。
In the LU decomposition, rounding errors accumulate. To suppress it, the solution is improved by the following way after obtaining the tantative solution by the procedure above. Let \(\myvector{x_{true}}\) be the true solution of eq. (\ref{eq.problem.corrected}), and \(\myvector{x_{numerical}}\) be the tantative solution obtained by the procedure above. The error of the tantitave solution is given by (\ref{eq.delta_x.def}\). Then eq. (\ref{eq.improvement}) holds, which can be divided to (\ref{eq.improvement.L}) and (\ref{eq.improvement.U}). The error \(\myvector{\delta x}\) is estimated by first solving for \(\myvector{\delta y}\) using eq. (\ref{eq.improvement.L}) with a forward substitution and then solving for \(\myvector{\delta x}\) using eq. (\ref{eq.improvement.U}) with a backward substitution. Substituting this \(\myvector{\delta x}\) and the tantative solution of \(\myvector{x_{numerical}}\) into eq. (\ref{eq.delta_x.def}), the value of \(\myvector{x_{true}}\) can be estimated.

以上まとめると方程式(\ref{eq.problem})を解く手順は以下のように整理できる。
Therefore the steps to solve eq. (\ref{eq.problem}) can be summarized as follows:

\(\myvector{A’}\)のLU分解(\ref{eq.LUdecomposition})を求める。
Compute the LU decomposition of \(\myvector{A’}\) (eq. \ref{eq.LUdecomposition}).
そのLU分解の際に\(\myvector{A}\)に対して行ったのと同じ行の交換を \(\myvector{b}\)に対して行い、\(\myvector{b’}\)を求める。
Compute \(\myvector{b’}\) by swapping the rows of \(\myvector{b}\), where the same swapping as that from \(\myvector{A}\) to \(\myvector{A’}\) is applied.
ステップ1で得た\(\myvector{A’}\)のLU分解とステップ2で得た\(\myvector{b’}\)を用いて (\ref{eq.problem.L})式の\(\myvector{y}\)を前進代入により求める。
Solve eq. (\ref{eq.problem.L}) for \(\myvector{y}\) with a forward substitution using the LU decomposition of \(\myvector{A’}\) obtained by step 1 and \(\myvector{b’}\) obtained by step 2.
ステップ3で得た\(\myvector{y}\)を用いて (\ref{eq.problem.U})式の\(\myvector{x}\)を後退代入により求める。
Solve eq. (\ref{eq.problem.U}) for \(\myvector{x}\) with a backward substitution using \(\myvector{y}\) obtained by step 3.
ステップ4で得た\(\myvector{x}\)を\(\myvector{x_{numerical}}\)として (\ref{eq.improvement.L})式の右辺を計算する。
Compute the right hand side of eq. (\ref{eq.improvement.L}), where the \(\myvector{x}\) value obtained by step 4 is used as \(\myvector{x_{numerical}}\).
ステップ5の結果を用いて (\ref{eq.improvement.L})式の\(\myvector{\delta y}\)を前進代入により求める。
Solve eq. (\ref{eq.improvement.L}) for \(\myvector{\delta y}\) with a forward substitution using the result of step 5.
ステップ6で得た\(\myvector{\delta y}\)を用いて (\ref{eq.improvement.U})式の\(\myvector{\delta x}\)を後退代入により求める。
Solve eq. (\ref{eq.improvement.U}) for \(\myvector{\delta x}\) with a backward substitution using \(\myvector{\delta y}\) obtained by step 6.
\(\myvector{x_{numerical}}\)から\(\myvector{\delta x}\)を引いて最終解\(\myvector{x_{true}}\)を得る(\ref{eq.delta_x.def}式)。
Calculate the final solution \(\myvector{x_{true}}\) by substituting \(\myvector{\delta x}\) from \(\myvector{x_{numerical}}\) (eq. \ref{eq.delta_x.def}).

このうちのステップ1(LU分解)は関数呼び出し前に予め行っておく必要がある。関数内ではステップ2-8が行われる。各ステップに対応する補助関数を用意してあるので (ステップ2はcolumnvector_swap_by_pivot、ステップ3,6はforward_substitution、ステップ4,7はbackword_substitution (いずれもmatrix/inverse.h))、この関数のコードはこれらの補助関数を組合せるだけの簡単なものになっている。上記の計算式に登場した変数と関数のソースコードで用いている変数名との対応関係を下表に示す。
The step 1 (LU decomposition) needs to be performed in advance to the function call. Within this function, steps 2-8 are conducted. The source code of this function is relatively simple, since individual steps can be executed by corresponding functions (step 2 by function columnvector_swap_by_pivot, steps 3 and 6 by function forward_substitution, and steps 4 and 7 by function backword_substitution; all of them are in matrix/inverse.h). Relations between the variables in the formula above and those in the source code of this function are listed in the table below.

計算式における変数名 Variable name in the formula	関数のソースコードにおける変数名 Variable name in the source code of this function
\(\myvector{A}\)	A
\(\myvector{b}\)	b
\(\myvector{A’}\)	使用しない Not used
\(\myvector{b’}\)	bb
\(\myvector{L}\)	AA_LU.L
\(\myvector{U}\)	AA_LU.U
\(\myvector{y}\)	y
\(\myvector{x_{numerical}}\)	x_numerical
\(\myvector{U}\myvector{x_{numerical}}\)	U_times_x
\(\myvector{A’}\myvector{x_{numetical}}\)	AA_times_x
\(\myvector{A’}\myvector{x_{numetical}}-\myvector{b’}\)	right
\(\myvector{\delta y}\)	dy
\(\myvector{\delta x}\)	dx
\(\myvector{x_{true}}\)	re