Linear Algebra with Sage

<LU-Decomposition>


Made by SKKU Linear Algebra Lab (2011)



Solve the following LSE(@@A \textbf{x} = \textbf{v}@@) using LU-decomposition (@@A \textbf{x} = LU \textbf{x} = \textbf{v}@@).
 

  @@\begin{bmatrix} 1 & 1 & 2\\ 2 & 4 & -3\\ 3 & 6 & -5 \end{bmatrix} \begin{bmatrix} x\\y\\z\end{bmatrix} =\begin{bmatrix} 9\\1\\0\end{bmatrix}@@



Step 1 : Find LU-decomposition of @@A@@

(주어진 행렬 @@A@@는 full rank matrix 이므로 permutation matrix @@P@@를 곱하지 않아도 LU분해를 할 수 있다. 하지만 주어진 행렬 @@A@@가 full rank matrix가 아닌 경우에는 적당한 permutation matrix @@P@@를 먼저 곱하고 LU분해를 시행한다.)


Define matrices (행렬 생성 및 확인)

A=matrix(QQ,[[1,1,2],[2,4,-3],[3,6,-5]]);

I=matrix(QQ,[[1,0,0],[0,1,0],[0,0,1]]);

v=matrix(3,1,[9,1,0]);

A

[ 1  1  2]

[ 2  4 -3]

[ 3  6 -5]


<Sage를 이용하여 간결한 계산(Find @@ L = L^{'-1} @@, @@U@@)을 하기 위해, 다음과 같이 단위행렬과 행렬 @@A@@로 이루어진 첨가행렬을 이용한다.>


@@(I | PA)@@
@@\Rightarrow E_1 (I | PA) = (E_1 I | E_1 PA)@@
@@\Rightarrow E_2 E_1 (I | PA) = (E_2 E_1 I | E_2 E_1 PA)@@
@@\Rightarrow \cdots@@
@@\Rightarrow E_n \cdots E_2 E_1 (I | PA) = ( E_n \cdots E_2 E_1 I | E_n \cdots E_2 E_1 PA)@@


Let  @@( E_n \cdots E_2 E_1 I | E_n \cdots E_2 E_1 PA)=(L^{'} | U)@@, then

   @@E_n \cdots E_2 E_1 = L^{'}@@, @@E_n \cdots E_2 E_1 PA=U @@.

@@\Rightarrow PA = (E_n \cdots E_2 E_1)^{-1} U @@
@@\Rightarrow PA = (L^{'-1}) U @@


Since @@E_n \cdots E_2 E_1@@ is lower triangular, @@L^{'-1}@@is also lower triangular.

@@\Rightarrow PA=(L^{'-1})U = LU@@


Make a augmented matrix @@[I@@ @@|@@ @@A]@@ (수반행렬 생성)

IA=I.augment(A);IA

[ 1  0  0  1  1  2]

[ 0  1  0  2  4 -3]

[ 0  0  1  3  6 -5]


1st ERO

IA.add_multiple_of_row(1, 0, -2);IA

[ 1  0  0  1  1  2]

[-2  1  0  0  2 -7]

[ 0  0  1  3  6 -5]


2nd ERO

IA.add_multiple_of_row(2, 0, -3);IA

[  1   0   0   1   1   2]

[ -2   1   0   0   2  -7]

[ -3   0   1   0   3 -11]


3rd ERO (첨가행렬 안의 행렬 @@A@@가 상삼각행렬이 될 때까지 기본행연산을 한다.
Upper triangular matrix: @@E_3 E_2 E_1 A = U @@ and Lower triangular matrix: @@E_3 E_2 E_1 = L^{'}@@)

IA.add_multiple_of_row(2, 1, -3/2);IA

[   1    0    0    1    1    2]

[  -2    1    0    0    2   -7]

[   0 -3/2    1    0    0 -1/2]


Define (lower triangular) matrix @@L@@ (하삼각행렬 @@L = L^{'-1} = (E_3 E_2 E_1)^{-1}@@)

L=IA.submatrix(0,0,3,3).inverse();L # submatrix(0,0,3,3) = 1행 1열의 성분으로 부터 3행, 3열로 이루어진 부분행렬

[  1   0   0]

[  2   1   0]

[  3 3/2   1]


Define (upper triangular) matrix @@U@@ of @@[I@@ @@|@@ @@A]@@ (상삼각행렬 @@U=E_3 E_2 E_1 A@@)

U=IA.submatrix(0,3,3,3);U # submatrix(0,3,3,3) = 1행 4열의 성분으로 부터 3행, 3열로 이루어진 부분행렬

[   1    1    2]

[   0    2   -7]

[   0    0 -1/2]


@@A=LU@@ (LU분해 확인)

L*U

[ 1  1  2]

[ 2  4 -3]

[ 3  6 -5]



Step 2 : Solve @@L \textbf{y} = \textbf{v}@@ (@@A \textbf{x} = (LU)\textbf{x} = L(U\textbf{x}) = L \textbf{y} = \textbf{v}@@)

y=L\v; y

[   9]

[ -17]

[-3/2]



Step 3 : Solve @@ U \textbf{x} = \textbf{y}@@ (@@  L (U \textbf{x}) = L \textbf{y} = \textbf{v}@@)

x=U\y; x

[1]

[2]

[3]



Theorem. Let @@A@@ be an invertible matrix. Then for fixed suitable permutation matrix @@P@@ the matrix @@PA@@ has unique @@LDU@@ factorization.



* Command for LU-Decomposition of @@B@@ : (LU분해에 관한 자동 명령어)


Redefine (RDF type) matrix (RDF형식의 행렬 생성 및 확인)

B=matrix(RDF,[[1,1,2],[2,4,-3],[3,6,-5]]);

B

[ 1.0  1.0  2.0]

[ 2.0  4.0 -3.0]

[ 3.0  6.0 -5.0]


@@PB=LU@@ (LU분해)

P,L,U=B.LU()

print P

print

print L

print

print U

[0.0 0.0 1.0]

[1.0 0.0 0.0]

[0.0 1.0 0.0]


[           1.0            0.0            0.0]

[0.333333333333            1.0            0.0]

[0.666666666667            0.0            1.0]


[           3.0            6.0           -5.0]

[           0.0           -1.0  3.66666666667]

[           0.0            0.0 0.333333333333]