行列のLU分解（LU decompotision）とは、正方行列を下三角行列（lower triangular matrix）と上三角行列（upper triangular matrix）の積に分解することです。
正方行列がLU分解可能であるならば、下三角行列の全ての対角成分が $1$ であるようなLU分解は一意に定まります。すなわち、$n$ 次正方行列 $A$ がLU分解可能であるならば、
$$ L= \begin{pmatrix} 1 & 0 & 0 & \cdots & 0\\ l_{21} & 1 & 0 & \cdots & 0\\ l_{31} & l_{32} & 1 & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & 0\\ l_{n1} & l_{n2} & \cdots & l_{n(n-1)} & 1 \end{pmatrix} ,\quad U= \begin{pmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n}\\ 0 & u_{22} & u_{23} & \cdots & u_{2n}\\ 0 & 0 & u_{33} & \cdots & u_{3n} \\ \vdots & \vdots & \ddots & \ddots & \vdots\\ 0 & 0 & \cdots & 0 & u_{nn} \end{pmatrix} $$の形の下三角行列 $L$ と上三角行列 $U$ を用いて
$$ A=LU $$と一意的に表現できます。
正方行列がLU分解可能であるための必要十分条件は、その全ての首座小行列の行列式が $0$ でないことです。すなわち、$n$ 次正方行列
$$ A= \begin{pmatrix} a_{11} & a_{12} & a_{13} & \cdots & a_{1n}\\ a_{21} & a_{22} & a_{23} & \cdots & a_{2n}\\ a_{31} & a_{32} & a_{33} & \cdots & a_{3n} \\ \vdots & \vdots & \vdots & \ddots & \vdots\\ a_{n1} & a_{n2} & a_{n3} & \cdots & a_{nn} \end{pmatrix} $$に対して
$$ A_1= \begin{pmatrix} a_{11} \end{pmatrix} ,\quad A_2= \begin{pmatrix} a_{11} & a_{12}\\ a_{21} & a_{22} \end{pmatrix} ,\quad A_3= \begin{pmatrix} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{23}\\ a_{31} & a_{32} & a_{33} \end{pmatrix} ,\quad \ldots,\quad A_n=A $$とおくとき、$A$ がLU分解可能であるための必要十分条件は
$$ \forall i\in\{1,\ldots,n\}, \ \det A_i\neq 0 $$です。これは、ガウスの消去法の前進消去の過程でピボット選択を行う必要がない（注目する対角成分に0が現れない）場合に相当します。
後述するように、$A$ のLU分解が分かってさえいれば、連立一次方程式 $Ax=b$ の解を比較的簡単に求めることができます。そこで、LU分解をどのように求めればよいかを次に扱います。
行列のLU分解を求める方法は複数知られていますが、ここではガウスの消去法の前進消去を用いる方法を扱います。以下では、考える行列 $A$ がLU分解可能であるという前提の下で話を進めることにします。
ピボット選択を行わないガウスの消去法の前進消去においては、$n$ 次正方行列 $A$ に対して
という行基本変形を繰り返して行うことで、上三角行列 $U$ まで変形します。
ここで、
という行基本変形は、単位行列の $(i,j)$ 成分を $\alpha_{ij}$ に変えた行列
$$ E_{ij}= \begin{pmatrix} 1 & 0 & \cdots & \cdots & 0\\ 0 & 1 & \ddots & & \vdots\\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \vdots & \alpha_{ij} & \ddots & 1 & 0\\ 0 & \cdots & \cdots & 0 & 1 \end{pmatrix} $$を左から掛けることに対応します。この行列 $E_{ij}$ は正則であり、逆行列
$$ E_{ij}^{-1}= \begin{pmatrix} 1 & 0 & \cdots & \cdots & 0\\ 0 & 1 & \ddots & & \vdots\\ \vdots & \ddots & \ddots & \ddots & \vdots \\ \vdots & -\alpha_{ij} & \ddots & 1 & 0\\ 0 & \cdots & \cdots & 0 & 1 \end{pmatrix} $$をもちます。
したがって、ピボット選択を行わないガウスの消去法の前進消去は
$$ E_{n(n-1)}\left(E_{n(n-2)}E_{(n-1)(n-2)}\right)\cdots\left(E_{n2}\cdots E_{32}\right)\left(E_{n1}\cdots E_{21}\right)A=U $$と表現することができ、
$$ L=\left(E_{21}^{-1}\cdots E_{n1}^{-1}\right)\left(E_{32}^{-1}\cdots E_{n2}^{-1}\right)\cdots\left(E_{(n-1)(n-2)}^{-1}E_{n(n-2)}^{-1}\right)E_{n(n-1)}^{-1} $$とおけば、$L$ は全ての対角成分が $1$ であるような下三角行列
$$ L= \begin{pmatrix} 1 & 0 & 0 & \cdots & 0\\ -\alpha_{21} & 1 & 0 & \cdots & 0\\ -\alpha_{31} & -\alpha_{32} & 1 & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & 0\\ -\alpha_{n1} & -\alpha_{n2} & \cdots & -\alpha_{n(n-1)} & 1 \end{pmatrix} $$となり、LU分解 $A=LU$ が得られます。
上の理論的内容をまとめると、行列 $A$ に対してピボット選択を行わないガウスの消去法の前進消去を適用した変形結果が上三角行列 $U$ になり、途中の行基本変形で考えるスカラー $\alpha_{ij}$ の符号を変えた $-\alpha_{ij}$ を集めることで下三角行列 $L$ が得られることが分かりました。
前回作成したプログラムでは、前進消去が終わった段階の変数 $A$ に上三角行列 $U$ の要素が格納されていましたが、「値が $0$ で記憶する必要のない対角より下の成分」に 「$L$ において記憶する必要のある対角より下の成分」を格納することで、メモリの節約ができます。
つまり、前進消去が終わった段階の変数 $A$ が
$$ \begin{pmatrix} u_{11} & u_{12} & u_{13} & \cdots & u_{1n}\\ l_{21} & u_{22} & u_{23} & \cdots & u_{2n}\\ l_{31} & l_{32} & u_{33} & \cdots & u_{3n} \\ \vdots & \vdots & \ddots & \ddots & \vdots\\ l_{n1} & l_{n2} & \cdots & l_{n(n-1)} & u_{nn} \end{pmatrix} $$となるようにします。
なお、LU分解においては、前回と違って行列 $A$ の行基本変形に合わせたベクトル $b$ の変形は不要です。
例として、$4$ 次正方行列
$$ A= \begin{pmatrix} 2 & 1 & -2 & 3\\ 2 & 0 & 4 & 1\\ 3 & 0 & 5 & 2\\ 1 & -1 & 2 & 1 \end{pmatrix} $$のLU分解を考えます。
A = [2,1,-2,3; 2,0,4,1; 3,0,5,2; 1,-1,2,1]
n = size(A,1)
A = 2 1 -2 3 2 0 4 1 3 0 5 2 1 -1 2 1 n = 4
第1行を用いた変形は、次のようになります。
alpha = -A(2,1)/A(1,1) %Aの(2,1)成分を0にするためのスカラー
A(2,2:4) = A(2,2:4)+alpha*A(1,2:4); %Aの第2行に第1行のalpha倍を加える（第1列の計算はしない）
A(2,1) = -alpha %下三角行列の成分を格納する
alpha = -A(3,1)/A(1,1) %Aの(3,1)成分を0にするためのスカラー
A(3,2:4) = A(3,2:4)+alpha*A(1,2:4); %Aの第3行に第1行のalpha倍を加える（第1列の計算はしない）
A(3,1) = -alpha %下三角行列の成分を格納する
alpha = -A(4,1)/A(1,1) %Aの(4,1)成分を0にするためのスカラー
A(4,2:4) = A(4,2:4)+alpha*A(1,2:4); %Aの第4行に第1行のalpha倍を加える（第1列の計算はしない）
A(4,1) = -alpha %下三角行列の成分を格納する
alpha = -1 A = 2 1 -2 3 1 -1 6 -2 3 0 5 2 1 -1 2 1 alpha = -1.5000 A = 2.0000 1.0000 -2.0000 3.0000 1.0000 -1.0000 6.0000 -2.0000 1.5000 -1.5000 8.0000 -2.5000 1.0000 -1.0000 2.0000 1.0000 alpha = -0.50000 A = 2.00000 1.00000 -2.00000 3.00000 1.00000 -1.00000 6.00000 -2.00000 1.50000 -1.50000 8.00000 -2.50000 0.50000 -1.50000 3.00000 -0.50000
第2行を用いた変形は、次のようになります。
alpha = -A(3,2)/A(2,2) %Aの(3,2)成分を0にするためのスカラー
A(3,3:4) = A(3,3:4)+alpha*A(2,3:4); %Aの第3行に第2行のalpha倍を加える（第1,2列の計算はしない）
A(3,2) = -alpha %下三角行列の成分を格納する
alpha = -A(4,2)/A(2,2) %Aの(4,2)成分を0にするためのスカラー
A(4,3:4) = A(4,3:4)+alpha*A(2,3:4); %Aの第4行に第2行のalpha倍を加える（第1,2列の計算はしない）
A(4,2) = -alpha %下三角行列の成分を格納する
alpha = -1.5000 A = 2.00000 1.00000 -2.00000 3.00000 1.00000 -1.00000 6.00000 -2.00000 1.50000 1.50000 -1.00000 0.50000 0.50000 -1.50000 3.00000 -0.50000 alpha = -1.5000 A = 2.00000 1.00000 -2.00000 3.00000 1.00000 -1.00000 6.00000 -2.00000 1.50000 1.50000 -1.00000 0.50000 0.50000 1.50000 -6.00000 2.50000
第3行を用いた変形は、次のようになります。
alpha = -A(4,3)/A(3,3) %Aの(4,3)成分を0にするためのスカラー
A(4,4:4) = A(4,4:4)+alpha*A(3,4:4); %Aの第4行に第3行のalpha倍を加える（第1,2,3列の計算はしない）
A(4,3) = -alpha %下三角行列の成分を格納する
alpha = -6 A = 2.00000 1.00000 -2.00000 3.00000 1.00000 -1.00000 6.00000 -2.00000 1.50000 1.50000 -1.00000 0.50000 0.50000 1.50000 6.00000 -0.50000
ガウスの消去法の前進消去によるLU分解のアルゴリズムを実現するコードを書き、
$$ A= \begin{pmatrix} 2 & 1 & -2 & 3\\ 2 & 0 & 4 & 1\\ 3 & 0 & 5 & 2\\ 1 & -1 & 2 & 1 \end{pmatrix} $$に対する変形の結果を表示してください。ただし、$A$ の次数が $4$ でない場合にも正しく動作するコードにすること。
%演習1のコード
A = [2,1,-2,3; 2,0,4,1; 3,0,5,2; 1,-1,2,1];
n = size(A,1);
%ここに、不足しているコードを書く
%for文のネスト（for文の中にfor文を入れる）を使うべき
%4ではなく変数nを使うことが重要
連立一次方程式 $Ax=b$ の解を求めることを考えます。LU分解 $A=LU$ が得られているとすると、
$$ LUx=b $$において $y=Ux$ とおくことにより、二つの連立一次方程式
$$ Ly=b,\quad Ux=y $$が得られます。
$L$ と $b$ を用いて $Ly=b$ の解 $y$ を求めた後に、$U$ と $y$ を用いて $Ux=y$ の解 $x$ を求めれば、元々の $Ax=b$ の解が求まったことになります。
ここで、$L$ は全ての対角成分が $1$ であるような下三角行列なので、
$$ L= \begin{pmatrix} 1 & 0 & 0 & \cdots & 0\\ l_{21} & 1 & 0 & \cdots & 0\\ l_{31} & l_{32} & 1 & \ddots & \vdots \\ \vdots & \vdots & \ddots & \ddots & 0\\ l_{n1} & l_{n2} & \cdots & l_{n(n-1)} & 1 \end{pmatrix} ,\quad b= \begin{pmatrix} b_1\\ b_2\\ b_3\\ \vdots\\ b_n \end{pmatrix} $$とすれば、次のようにして $Ly=b$ の解 $y$ の成分 $y_i$ を前から順に求めることができます。これを前進代入（forward substitution）と言います。
$y_1=b_1$
$y_2=b_2-l_{21}y_1$
$y_3=b_3-l_{31}y_1-l_{32}y_2$
$\vdots$
$y_n=b_n-l_{n1}y_1-\cdots-l_{n(n-1)}y_{n-1}$
その後の $Ux=y$ については、$U$ が上三角行列なので前回扱った後退代入を用いて解 $x$ を求めることができます。
なお、プログラム上では、前進消去が終わった段階で下三角行列 $L$ と 上三角行列 $U$ の成分が変数 $A$ に一緒に格納されていることに注意してください。
演習1のコードに前進代入と後退代入のアルゴリズムを実現するコードを追加し、
$$ A= \begin{pmatrix} 2 & 1 & -2 & 3\\ 2 & 0 & 4 & 1\\ 3 & 0 & 5 & 2\\ 1 & -1 & 2 & 1 \end{pmatrix} ,\quad b= \begin{pmatrix} 11\\ -5\\ -5\\ 4 \end{pmatrix} $$に対して $Ax=b$ の解を求めてください。ただし、$A$ と $b$ の次数が $4$ でない場合にも正しく動作するコードにすること。
%演習2のコード
前進代入の過程は、ガウスの消去法の前進消去において $A$ と同じ行基本変形を $b$ に対して行うことに相当します。
このことから、ガウスの消去法の前進消去による $LU$ 分解を利用して前進代入と後退代入で $Ax=b$ の解を求める上記の方法は、計算量の観点では前回扱ったピボット選択を行わないガウスの消去法と変わらないと言えます。
上記の方法が威力を発揮するのは、一つの行列 $A$ と複数のベクトル $b_1,\ldots,b_m$ に対して、複数の連立一次方程式
$$ Ax=b_1,\quad\ldots,\quad Ax=b_m $$の解をそれぞれ求めるような場合です。
演習1～演習2に取り組んでください。