[개발자를 위한 실전 선형대수학] Chapter 9 연습문제

[chapter9] 행 축소와 LU 분해: 선형대수학의 핵심 분해법 2

9-1

LU 분해 실행하는 데 걸리는 시간 측정해보기

정답 코드

import time

tic = time.time()

for i in range(1000):
  A = np.random.randn(100,100)
  P,L,U = scipy.linalg.lu(A)

toc = time.time() - tic
toc # 0.6329073905944824

LU 분해와 QR 분해는 각각 어떨때 쓸까? with ChatGPT

LU 분해가 QR 분해보다 더 빠르다면, QR 분해는 안쓰는걸까? 가 궁금해져서 ChatGPT에게 물어봤다.

LU분해는 L: 아래 삼각 / U: 위 삼각으로 분해하는 방법으로, 계산량이 적어 빠르다. 하지만 수치적으로 불안정할 수 있고, A가 특이하면 결과가 망가질 수 있다는 단점이 있다.

따라서 LU 분해는 속도가 중요한 선형방정식 풀이용으로 사용한다.

QR분해는 Q: 직교 행렬 / R : 위 삼각행렬로 분해하는 방법으로, LU보다 계산량이 많아 느리다. 하지만 수치적으로 매우 안정적이고, 열 벡터의 독립성, 투영 구조가 깔끔하게 드러난다.

따라서 QR 분해는 정확성과 안정성이 중요한 문제용으로 사용한다.

QR분해는 최소제곱 문제, 데이터가 많고 잡음이 있을 때, 그람-슈미트, 직교기저, 차원 분석 등에서 주로 사용된다.

9-2

행렬 곱셈을 사용해 계수-3의 6x8 행렬을 만든다.

LU 분해 결과의 세 행렬을 행렬의 계수와 함께 제목에 표기한다. L의 대각선이 모두 1에 존재함에 주목한다.

정답 코드

M = 5
N = 7
r = 5

A = np.random.randn(M,r) @ np.random.randn(r,N)

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

_,axs = plt.subplots(1,3,figsize=(12,7))

axs[0].imshow(A,vmin=-1,vmax=1,cmap='gray')
axs[0].set_title(f'A, rank={np.linalg.matrix_rank(A)}')

axs[1].imshow(L,vmin=-1,vmax=1,cmap='gray')
axs[1].set_title(f'L, rank={np.linalg.matrix_rank(L)}')

axs[2].imshow(U,vmin=-1,vmax=1,cmap='gray')
axs[2].set_title(f'U, rank={np.linalg.matrix_rank(U)}')

plt.tight_layout()
plt.savefig('Figure_09_02.png',dpi=300)
plt.show()

여기서 L의 대각선에 모두 1인 것이 왜인지 궁금했다.

LU 분해는 유일하지 않기 때문에, 대각 성분을 L과 U 사이에서 마음대로 나눌 수 있다.

따라서 관례적으로 L의 대각선 = 1, U = 실제 피벗 크기를 모두 가지도록 설정해 자유도를 하나 제거해서 분해를 유일하게 만드는 것이다.

9-3

LU 분해를 활용하면 행렬식을 계산할 수 있다. 아래는 행렬식의 두 가지 속성이다.

• 삼각 행렬의 행렬식은 대각선의 곱이다
• 행렬 곱셈의 행렬식은 각 행렬식의 곱과 같다($det(AB) = det(A)det(B)$)

이 두 가지를 종합하면 행렬식은 L의 대각선 곱에 U의 대각선 곱을 곱한 값으로 계산할 수 있다.

이 때, L의 대각선의 곱은 모두 1이기 때문에 행렬 A의 행렬식은 단순히 U의 대각선의 곱이다.

정답 코드

M = 6
A = np.random.randn(M,M)

P,L,U = scipy.linalg.lu(A) 
# scipy에서는 PA = LU
# P = 치환행렬(행을 바꾼 기록을 모아둔 행렬)
# 따라서 계산 안정성을 위해 P를 곱한다.
# P의 행렬식 : 행을 한 번 교환(-1), 두 번 교환(+1) ...

detLU = np.prod(np.diag(U)) * np.linalg.det(P)
detNp = np.linalg.det(A)

print(detLU, detNp) # -14.420330528637672 -14.420330528637665

 

9-4

4x4 난수 행렬에 대해 scipy.linalg.lu를 호출한 결과를 사용해 아래 방정식을 직접 구현하기

$$ A^{-1} = U^{-1}L^{-1}P $$

$AA^{-1}$은 단위행렬일 수도 있고 아닐수도 있는데, P에 따라 달라진다. 이러한 불일치는 scipy.linalg.lu의 출력으로 인해 발생하는데, 수학 규칙이 아닌 Scipy의 규칙을 따르도록 수정하기

정답 코드

m = 4
A = np.random.randn(m,m)

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

right = np.linalg.inv(U) @ np.linalg.inv(L) @ P.T
left = np.linalg.inv(A)

np.round(A@left, 10)
#array([[ 1.,  0.,  0.,  0.],
#       [-0.,  1., -0.,  0.],
#       [-0., -0.,  1., -0.],
#       [-0., -0.,  0.,  1.]])

9-5

행렬 A = PLU에서, 치환 행렬을 사용하지 않고도 $A^TA$를 $U^TL^TLU$로 계산할 수 있다.

치환 행렬을 삭제할 수 있는 이유가 무엇인지 살펴보기.

그리고 파이썬에서 무작위 행렬을 사용해 $P\neq I$일 때에도 $A^TA = U^TL^TLU$임을 확인하기.

이유 살펴보기

LU 분해는 보통 아래처럼 쓸 수 있다.

$$ PA = LU $$

이 때, 치환행렬 P는 직교행렬이므로 아래와 같다.

$$ P^T = P^{-1}, P^TP = I $$

수식으로 바꿔보면

$ A = P^TLU $ 이고, 전치하면 $A^T = U^TL^TP$ 이다.

두 개를 곱하면 $A^TA = (U^TL^TP)(P^TLU)$가 되는데, $PP^T= I$이다.

따라서 아래처럼 된다.

$$ A^TA = U^TL^TLU$$

파이썬으로 확인하기

A = np.random.randn(4,4)

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

AtA_lu = U.T @ L.T @ L @ U
AtA_direct = A.T @ A

np.round(AtA_lu - AtA_direct, 10)
#array([[ 0.,  0.,  0., -0.],
#       [ 0., -0.,  0., -0.],
#       [-0.,  0.,  0., -0.],
#       [-0., -0.,  0., -0.]]) 같다!