오류에 대한 지적이나 질문, 토의 환영합니다. 자유롭게 댓글 남겨주세요!.!

---

선형대수 (Linear Algebra) (2): 가우스 소거법, LU분해

선형시스템의 해

$Ax = b$

• 해가 하나인 경우 $3x = 6$
• 해가 여러개인 경우 $0 \times x = 0$
• 해가 없는 경우 $0 \times x = 6$

• 해가 없으면($A=0$)이라면, $A$의 역행렬아 존재하지 않는 경우 $A$가 특이(singular)하다고 합니다.
• 해가 없으면, 선형시스템이 inconsistent 하다고 합니다.
• 해가 있으면, 선형시스템이 consistent 하다고 합니다.

---

$A$의 역행렬 구하기 - 가우스 소거법 (Gauss Elimination)

• 임의의 $m \times n$ 선형 시스템의 해를 구하는 가장 대표적인 방법입니다.
  1. 전방 소거법 (Forward Elimination)
     • 주어진 선형시스템을 아래로 갈수록 더 단순한 형태의 선형 방정식을 가지도록 변형하는 방법입니다.
  2. 후방 대입법 (Back-substitution)
     • 아래에서부터 위로 미지수를 실제 값으로 대체하는 방법입니다.

전방소거법 (Forward Elimination)

$Ax=b$
$\left[\begin{matrix}*\ * \ * \\*\ * \ * \\*\ * \ * \end{matrix}\right] \left[\begin{matrix}x_{1} \\x_{2} \\x_{3} \\\end{matrix}\right] =\left[\begin{matrix}* \\* \\* \\\end{matrix}\right]$

• 위와 같은 A의 꼴을 아래와 같이 변형하는 방법입니다.

$\left[\begin{matrix}*\ * \ * \\0\ * \ * \\0\ \ 0 \ * \end{matrix}\right] \left[\begin{matrix}x_{1} \\x_{2} \\x_{3} \\\end{matrix}\right] =\left[\begin{matrix}* \\* \\* \\\end{matrix}\right]$

• 이후 후방 대입법을 이용하여 $x_{1}, x_{2}, x_{3}$의 값을 찾아갑니다.
• 이렇게 가장 아래 식에 의해서 $x_{3} = *$이 구해지고, 두번 째 식에 의해서 $x_{2} = *$이 구해지고, 가장 위의 식에 의해서 $x_{1} = *$이 구해집니다.

• 전방 소거법에 사용되는 기본행 연산
  1. 치환 (replacement)
  2. 교환 (interchange)
  3. 스케일링 (scailing)

• 전방소거법의 가치

  • 주어진 선형시스템을 가장 풀기 쉬운 꼴로 변형해줍니다.
  • 주어진 선형시스템의 랭크(Rank: 복잡한지 아닌지)를 알려줍니다.
  • 선형시스템이 해가 있는지 아니면 없는지 알려줍니다.

---

행렬분해: 인수분해의 행렬 버전

행렬분해 (1): LU 분해

• 가우스 소거법의 전방소거법을 행렬로 코드화한 것입니다.

• 행렬 $A$를 $PLU$꼴로 변형시키는 분해입니다.

  • $A = P$(interchange) $\times$ $L$(lower triangular matrix) $\times$ $U$(Upper triangular matrix)

• LU 분해의 장점

  $Ax = b$ <=> $(LU)x = b$ <=> $L(Ux) = b$
  => $L(Ux) = b$ => $Ux$: 전방 대치법으로 쉽게 계산 가능
  => $Ly = b \ (단, Ux = y)$ => $Ly$: 후방 대치법으로 쉽게 계산 가능

• LU 분해의 활용

  • 수치적 안정성: 역행렬을 구하는 것보다 수치적으로 안정적입니다.
  • b가 자주 업데이트되는 경우 사용합니다.

LU 분해 실습

import numpy as np
import scipy        
import scipy.linalg   # LU 분해를 사용하기 위한 import

행렬 A, 벡터 b 코딩

A = np.array([[3, 1, 1], [1, -2, -1], [1, 1, 1]])
b = np.array([4, 1, 2])

print("A:", A)
print(np.shape(A))

print("b:", b)
print(np.shape(b))

A: [[ 3  1  1]
 [ 1 -2 -1]
 [ 1  1  1]]
(3, 3)
b: [4 1 2]
(3,)

LU 분해 결과 확인하기

P, L, U = scipy.linalg.lu(A)

print("P:", P)
print("L:", L)
print("U:", U)

AA = P @ L @ U
print("AA:", AA)

P: [[1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]]
L: [[ 1.          0.          0.        ]
 [ 0.33333333  1.          0.        ]
 [ 0.33333333 -0.28571429  1.        ]]
U: [[ 3.          1.          1.        ]
 [ 0.         -2.33333333 -1.33333333]
 [ 0.          0.          0.28571429]]
AA: [[ 3.  1.  1.]
 [ 1. -2. -1.]
 [ 1.  1.  1.]]

LU 분해를 이용한 선형시스템 $Ax = b$ 풀기

lu, piv = scipy.linalg.lu_factor(A)
x = scipy.linalg.lu_solve((lu, piv), b)

print("x:", x)
print(np.shape(x))

bb = A@x
print("bb:", bb)

x: [ 1. -1.  2.]
(3,)
bb: [4. 1. 2.]

bbb = L@U@x

print('bbb:', bbb)

[4. 1. 2.]

행렬의 Rank 실습

1. $rank = 2$인 $2 \times 2$ 행렬 $A$

A = np.array([[1, 3], [-2, 1]])

print("A:", A)

A: [[ 1  3]
 [-2  1]]

print("rank:", np.linalg.matrix_rank(A))
A_inv = np.linalg.inv(A)  

print(A_inv)

rank: 2
[[ 0.14285714 -0.42857143]
 [ 0.28571429  0.14285714]]

LU 분해의 결과를 행렬로 확인하기

P, L, U = scipy.linalg.lu(A)

print("P:", P)
print("L:", L)
print("U:", U)

A_ = P @ L @ U
print("A_:", A_)

P: [[0. 1.]
 [1. 0.]]
L: [[ 1.   0. ]
 [-0.5  1. ]]
U: [[-2.   1. ]
 [ 0.   3.5]]
A_: [[ 1.  3.]
 [-2.  1.]]

b = np.array([2, 4])
# LU 분해
lu, piv = scipy.linalg.lu_factor(A)
x = scipy.linalg.lu_solve((lu, piv), b)

print("x:", x)
print(np.shape(x))

x: [-1.42857143  1.14285714]
(2,)

2. $rank = 1$인 $2 \times 2$ 행렬 $AA$

AA = np.array([[1, 3], [2, 6]])

print("AA:", AA)

AA: [[1 3]
 [2 6]]

• AA는 rank가 1이므로 역행렬이 존재하지 않습니다.
  그러므로 역행렬을 구하려고 하면 에러가 발생됩니다.

print("rank:", np.linalg.matrix_rank(AA))

AA_inv = np.linalg.inv(AA)
print(AA_inv)

rank: 1
---------------------------------------------------------------------------
LinAlgError                               Traceback (most recent call last)
<ipython-input-32-79b2f18d4fad> in <module>
      1 print("rank:", np.linalg.matrix_rank(AA))
----> 2 AA_inv = np.linalg.inv(AA)

<__array_function__ internals> in inv(*args, **kwargs)

c:\users\lmh-laptop\appdata\local\programs\python\python38\lib\site-packages\numpy\linalg\linalg.py in inv(a)
    543     signature = 'D->D' if isComplexType(t) else 'd->d'
    544     extobj = get_linalg_error_extobj(_raise_linalgerror_singular)
--> 545     ainv = _umath_linalg.inv(a, signature=signature, extobj=extobj)
    546     return wrap(ainv.astype(result_t, copy=False))
    547 

c:\users\lmh-laptop\appdata\local\programs\python\python38\lib\site-packages\numpy\linalg\linalg.py in _raise_linalgerror_singular(err, flag)
     86 
     87 def _raise_linalgerror_singular(err, flag):
---> 88     raise LinAlgError("Singular matrix")
     89 
     90 def _raise_linalgerror_nonposdef(err, flag):

LinAlgError: Singular matrix

LU 분해의 결과를 행렬로 확인하기

P, L, U = scipy.linalg.lu(AA)

print("P:", P)
print("L:", L)
print("U:", U)

AA_ = P @ L @ U
print("AA_:", AA_)

P: [[0. 1.]
 [1. 0.]]
L: [[1.  0. ]
 [0.5 1. ]]
U: [[2. 6.]
 [0. 0.]]
AA_: [[1. 3.]
 [2. 6.]]

• 위와 마찬가지로 AA는 역행렬이 존재하지 않기 때문에,
  LU 분해를 호출하는 과정에서 에라가 발생됩니다.

b = np.array([2, 4])
# LU 분해
lu, piv = scipy.linalg.lu_factor(AA)
x = scipy.linalg.lu_solve((lu, piv), b)

print("x:", x)
print(np.shape(x))

x: [nan nan]
(2,)
<ipython-input-35-3f70c98bbe5d>:3: LinAlgWarning: Diagonal number 2 is exactly zero. Singular matrix.
  lu, piv = scipy.linalg.lu_factor(AA)

---

이 글은 프로그래머스 스쿨 인공지능 데브코스 과정에서 공부한 내용을 바탕으로 정리한 글입니다.