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

[chapter7] 역행렬 : 행렬 방정식의 만능 키

7-1

역행렬의 역행렬은 원래 행렬이 된다. 파이썬을 사용해 이를 증명하기

내 코드

A = np.random.randn(3,3)
Ainv = np.linalg.inv(A)
Ainvinv = np.linalg.inv(Ainv)
print(Ainvinv)
print(A) # 두 개가 같다

7-2

7-3-3 절에서 설명한 전체 알고리즘을 구현하고 그림 7-3을 재현하기.

정답 코드

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

# 소행렬
M = np.zeros((n,n))
# 격자행렬
G = np.zeros((n,n))

# 여인수행렬 구현
for i in range(n):
  for j in range(n):
    rows = [True]*n
    rows[i] = False # i행 제거

    cols = [True]*n
    cols[j] = False # j열 제거

    M[i,j] = np.linalg.det(A[rows,:][:,cols]) # i행,j열 제거 후 행렬식
    G[i,j] = (-1)**(i+j) # 격자행렬 완성

# 여인수행렬
C = M * G

# 역행렬 = 여인수행렬의 전치에서 행렬식을 나눈다
Ainv = C.T / np.linalg.det(A)

# 역행렬2
Ainv2 = np.linalg.inv(A)

# 두 개가 같은지 확인
np.round(Ainv2 - Ainv, 8)

# array([[ 0., -0., -0., -0.],
#       [ 0., -0.,  0.,  0.],
#       [-0.,  0.,  0., -0.],
#       [ 0.,  0.,  0., -0.]])

7-4

넓은 행렬에 대한 오른쪽 역행렬을 유도하기. 넓은 행렬에 대해 그림 7-4 재현하기

정답 코드

W = np.random.randint(-10,11,size=(4,40))

print( f'This matrix has rank={np.linalg.matrix_rank(W)}\n\n' )

WWt = W@W.T
WWt_inv = np.linalg.inv(WWt)
print( np.round(WWt_inv@WWt,4) )

R = W.T @ WWt_inv

7-5

파이썬에서 의사역행렬(np.linalg.pinv)이 가역  행렬의 완전 역행렬(np.linalg.inv)과 같다는 것을 구현하기.

그리고 의사역행렬이 높은 최대열계수 행렬의 경우 왼쪽 역행렬과 같고, 넓은 최대계수 행렬의 경우 오른쪽 역행렬과 같다는 것을 나타내기.

내 코드

A = np.random.randint(-10,11,size=(4,4))

Ainv = np.linalg.inv(A)
Apinv = np.linalg.pinv(A)

np.round(Apinv - Ainv, 8)
# array([[-0., -0., -0., -0.],
#       [-0.,  0.,  0.,  0.],
#       [ 0.,  0.,  0.,  0.],
#       [ 0., -0.,  0., -0.]])

m,n = 14,4
# 높은 행렬
B = np.random.randint(-10,11,size=(m,n))
Bleft = np.linalg.inv(B.T @ B) @ B.T
Bpinv = np.linalg.pinv(B)
np.round(Bpinv - Bleft, 8)
#array([[ 0., -0., -0., -0., -0., -0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,
#         0.],
#       [ 0.,  0.,  0., -0., -0., -0.,  0., -0.,  0.,  0., -0.,  0.,  0.,
#         0.],
#       [ 0., -0., -0., -0., -0., -0.,  0.,  0.,  0.,  0., -0.,  0.,  0.,
#         0.],
#       [ 0., -0.,  0., -0.,  0., -0.,  0.,  0.,  0., -0., -0.,  0., -0.,
#        -0.]])

# 넓은 행렬
C = np.random.randint(-10, 11, size=(n,m))
Cright = C.T @ np.linalg.inv(C@C.T)
Cpinv = np.linalg.pinv(C)

np.round(Cpinv - Cright, 8)
# array([[ 0.,  0.,  0., -0.],
#       [ 0., -0.,  0., -0.],
#       [ 0., -0.,  0., -0.],
#       [-0., -0., -0.,  0.],
#       [ 0.,  0.,  0.,  0.],
#       [ 0., -0.,  0.,  0.],
#       [ 0.,  0.,  0., -0.],
#       [ 0.,  0.,  0., -0.],
#       [-0., -0., -0., -0.],
#       [-0., -0., -0.,  0.],
#       [ 0.,  0.,  0., -0.],
#       [-0., -0.,  0., -0.],
#       [-0., -0., -0., -0.],
#       [ 0.,  0.,  0., -0.]])

7-6

LIVE EVIL 규칙은 곱 행렬의 역행렬에 적용된다. 코드에서 이를 테스트해보기.

두 개의 정방 최대계수 행렬 A와 B를 생성한 다음, 유클리드 거리를 사용해 $(AB)^{-1}$, $A^{-1}B^{-1}$, $B^{-1}A^{-1}$을 비교하기.

정답 코드

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

op1 = np.linalg.inv(A@B)
op2 = np.linalg.inv(A) @ np.linalg.inv(B)
op3 = np.linalg.inv(B) @ np.linalg.inv(A)

# 유클리드 거리
dist12 = np.sqrt(np.sum( (op1-op2)**2 ))
dist13 = np.sqrt(np.sum((op1-op3) ** 2))

print(f'Distance between (AB)^-1 and (A^-1)(B^-1) is {dist12:.8f}')
print(f'Distance between (AB)^-1 and (B^-1)(A^-1) is {dist13:.8f}')

# Distance between (AB)^-1 and (A^-1)(B^-1) is 2.08813520
# Distance between (AB)^-1 and (B^-1)(A^-1) is 0.00000000
# 즉 (AB)^-1은 (B^-1)(A^-1)과 같다.(LIVE EVIL) 법칙

7-7

LIVE EVIL 규칙이 단방향 역행렬에도 적용되는지 테스트하기.

M,N = 14,4
A = np.random.randn(M,N)

op1 = np.linalg.inv(A@A.T)
op2 = np.linalg.inv(A) @ np.linalg.inv(A.T)
# 정방이 아니기 때문에 역행렬 정의 불가.

7-8

그림 7-6을 재현하는 코드 작성하기.

먼저 연습문제 6-3의 코드를 복사하고, 그림을 재현한 후 왼쪽 아래 원소를 1로 설정해 변환행렬을 비가역으로 만든다.

정답 코드

T = np.array([
    [1, .5],
    [0, .5]
])

Tinv = np.linalg.inv(T)

# np.linspace(start, end, num) = start부터 end까지 포함해서 num개의 숫자를 균등 간격으로 생성
# 원을 20등분한 각도를 하나씩 만들되, 시작점(0)과 끝점(2pi)을 중복하지 않게 함
theta = np.linspace(0, 2*np.pi - 2*np.pi/20, 20)
origPoints = np.vstack((np.cos(theta), np.sin(theta))) # 단위 원 위의 좌표들

# 변형
transformedPoints = T @ origPoints

# 다시 되돌리기
backTransformed = Tinv @ transformedPoints

# 이미지화
plt.figure(figsize=(6,6))
plt.plot(origPoints[0,:],origPoints[1,:],'ko',label='Original')
plt.plot(transformedPoints[0,:],transformedPoints[1,:],'s',
         color=[.7,.7,.7],label='Transformed')
plt.plot(backTransformed[0,:],backTransformed[1,:],'rx',markersize=15,
         color=[.7,.7,.7],label='Inverse-transformed')

plt.axis('square')
plt.xlim([-2,2])
plt.ylim([-2,2])
plt.legend()
plt.savefig('Figure_07_06.png',dpi=300)
plt.show()

실행하면 아래와 같이 나온다.

7-9

힐버트 행렬을 생성하고, 식 7-1에 맞춰 정수를 입력으로 받아 힐버트 행렬을 생성하는 파이썬 함수 작성하기.

그리고 그림 7-5를 재현하기.

이중 for문 사용한 버전, 사용하지 않은 버전(외적) 두 가지의 정확성 살펴보기.

def hilbmat(k):
  H = np.zeros((k,k))
  for i in range(k):
    for j in range(k):
      H[i,j] = 1/((i+1)+(j+1)-1) # 인덱스 0부터니까 +1 해준다.
  return H

def hilbmat2(k):
  # np.arange로 1부터 k까지 일정 간격으로 숫자를 만든다.
  # reshape로 열벡터로 바꿔준다.(-1 = 알아서 개수만큼 = 40)
  k = np.arange(1,k+1).reshape(1,-1)

  # 브로드캐스팅을 활용한다.
  # k가 4라면 k.T는 4x1 행렬, k는 1x4 행렬
  # k.T + k의 결과는 4x4
  # 따라서 i+j-1이 된다.
  return 1/(k.T+k-1)

print(hilbmat(5)), print(' ')
print(hilbmat2(5)), print(' ')
from scipy.linalg import hilbert
print( hilbert(5) )
# 3개가 모두 같다.

#[[1.         0.5        0.33333333 0.25       0.2       ]
# [0.5        0.33333333 0.25       0.2        0.16666667]
# [0.33333333 0.25       0.2        0.16666667 0.14285714]
# [0.25       0.2        0.16666667 0.14285714 0.125     ]
# [0.2        0.16666667 0.14285714 0.125      0.11111111]]
 
#[[1.         0.5        0.33333333 0.25       0.2       ]
# [0.5        0.33333333 0.25       0.2        0.16666667]
# [0.33333333 0.25       0.2        0.16666667 0.14285714]
# [0.25       0.2        0.16666667 0.14285714 0.125     ]
# [0.2        0.16666667 0.14285714 0.125      0.11111111]]
 
#[[1.         0.5        0.33333333 0.25       0.2       ]
# [0.5        0.33333333 0.25       0.2        0.16666667]
# [0.33333333 0.25       0.2        0.16666667 0.14285714]
# [0.25       0.2        0.16666667 0.14285714 0.125     ]
# [0.2        0.16666667 0.14285714 0.125      0.11111111]]

7-10

힐버트 행렬 함수를 사용해서 힐버트 행렬을 만들고, np.linalg.inv를 사용해서 역을 계산하고 두 행렬의 곱을 계산하기.

이 곱은 단위행렬과 일치해야 한다.

즉, 이 곱과 np.eye가 생성한 실제 단위 행렬 사이의 유클리드 거리는 0이어야 한다.

이 코드를 3x3에서 12x12에 이르는 다양한 행렬 크기 범위에 대해 수행되도록 for문제 넣는다.

각 행렬 크기에 대해 유클리드 거리와 힐버트 행렬의 조건수를 저장한다.(np.linalg.cond 사용)

마지막으로 이전 코드를 되풀이하지만 힐버트 행렬 대신 가우스 난수 행렬 사용하기.

모든 결과를 그래프로 나타내기.

내 코드

from scipy.linalg import hilbert

sizes = np.arange(3, 13)

distances = np.zeros((len(sizes),2))
condNumbers = np.zeros((len(sizes),2))

for idx, val in enumerate(sizes):
  # 힐버트 행렬
  A = hilbert(val)
  Ai = np.linalg.inv(A)
  AAi = A @ Ai
  gap = AAi - np.eye(val)
  gap_dist = np.sqrt(np.sum(gap**2))
  distances[idx,0] = gap_dist
  condNumbers[idx,0] = np.linalg.cond(A)

  # 랜덤
  B = np.random.randn(val, val)
  Bi = np.linalg.inv(B)
  BBi = B @ Bi
  Bgap = BBi - np.eye(val)
  Bgap_dist = np.sqrt(np.sum(Bgap **2))
  distances[idx,1] = Bgap_dist
  condNumbers[idx,1] = np.linalg.cond(B)

# 그래프 그리기
fig,axs = plt.subplots(1,2,figsize=(14,5))

h = axs[0].plot(sizes, np.log(distances), 's-', markersize=12)
h[0].set_color('k')
h[0].set_marker('o')
h[1].set_color('gray')

axs[0].legend(['Hilbert','Random'])
axs[0].set_xlabel('Matrix size')
axs[0].set_ylabel('Log Euclidan distance')
axs[0].set_title('Distance to identity matrix')

h = axs[1].plot(sizes,np.log(condNumbers),'s-',markersize=12)
h[0].set_color('k') # adjust the individual line colors and shapes
h[0].set_marker('o')
h[1].set_color('gray')

axs[1].legend(['Hilbert','Random'])
axs[1].set_xlabel('Matrix size')
axs[1].set_ylabel('Log Kappa')
axs[1].set_title('Matrix condition number')

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

출력하면 아래와 같이 나온다.

힐버트 행렬은 수치적으로 매우 불안정한 행렬이다.

따라서 크기가 조금만 커져도 조건수가 폭발적으로 커지기 때문에, 역행렬 계산 오차도 커진다.

그래프에서는 오차가 직선모양으로 상승하는 것을 볼 수 있다.

랜덤 정규분포 행렬은 일반적으로 크기가 커져도 조건수가 비교적 안정적이다.

그래프에서도 오차가 불규칙하게 변하는 것을 알 수 있다.

즉, 같은 크기의 행렬이라도 구조적인 특성에 따라 수치 계산 안정성이 달라진다는 것을 보여주는 문제였다.