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

[chapter8] 직교 행렬과 QR 분해: 선형대수학의 핵심 분해법 1

8-1

난수 행렬 Q를 생성하고 $Q^T$와 $Q^{-1}$을 계산하는 코드를 작성해서 아래 수식 구현하기

$$ Q^TQ = QQ^T = Q^{-1}Q = QQ^{-1} = I $$

정답 코드

A = np.random.randn(5,5)
Q = np.linalg.qr(A)[0]
Qt = Q.T
Qin = np.linalg.inv(Q)

QtQ = Qt @ Q
print(np.round(QtQ, 8)), print(' ')
QQt = Q @ Qt
print(np.round(QQt, 8)), print(' ')
QinQ = Qin @ Q
print(np.round(QinQ,8)), print(' ')
QQin = Q @ Qin
print(np.round(QQin, 8)), print(' ')

# 모두 단위행렬이 나온다.

8-2

그람-슈미트 과정 구현하기. 4x4 난수 행렬을 사용하고, 결과를 np.linalg.qr의 Q와 대조한다.

하우스홀더 변환과 같은 변환에는 근본적으로 부호 불확실성이 존재해서, 사소산 차이에 따라 벡터가 뒤집힐 수 있다.

파이썬의 Q에서 만든 Q를 빼고 만든 Q와 파이썬의 Q를 더하면 한 쪽의 0이아닌 열이 다른 쪽에서는 0이 된다.

정답 코드

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

Q = np.zeros((n,n))

for i in range(n):
    
    # 원본 벡터 복사
    Q[:,i] = A[:,i]
    
    a = A[:,i] 
    for j in range(i):
        q = Q[:,j] 
        Q[:,i]=Q[:,i]-np.dot(a,q)/np.dot(q,q)*q # q방향 성분을 계산한 뒤에 뺀다.
        # np.dot(a,q)/np.dot(q,q) = a가 q방향으로 얼마나 기울어있는지 구함(스칼라)
        # *q = 벡터를 빼야하기 때문에 q를 곱해 q방향의 벡터를 만듦
    
    # 정규화
    Q[:,i] = Q[:,i] / np.linalg.norm(Q[:,i])

Q2,R = np.linalg.qr(A)

print( np.round( Q-Q2 ,10) ), print(' ')
print( np.round( Q+Q2 ,10) )
#[[ 1.82397971 -0.          0.          0.30413865]
# [-0.20295416  0.          0.          0.79317595]
# [-0.7945823  -0.         -0.          0.44391235]
# [ 0.02337943 -0.         -0.         -1.75536703]]
 
#[[ 0.          0.6975557  -0.30661661  0.        ]
# [-0.         -0.57618663 -1.73138398  0.        ]
# [ 0.          1.75739636 -0.28833806  0.        ]
# [ 0.          0.30493122 -0.90838149  0.        ]]

# 둘 중 하나는 항상 0열이 나오는 열이 생김
# → 결국 Q와 Q2는 같은 직교공간을 표현한다는 뜻

8-3

거의 직교에 가깝지만 직교는 아닌 행렬에 QR 분해를 적용하면 어떤 일이 발생하는지 알아보기.

6x6 난수 행렬의 QR 분해로부터 U라는 직교행렬을 만들기

1. U의 QR분해를 계산하고 R = I임을 확인하기

2. U의 각 열의 노름을 수정하기. 1~6열의 노름을 10~15값의 노름으로 설정한다.(U의 첫 번째 열의 노름은 10, 두번째 열의 노름은 11)

변조된 U 행렬을 QR 분해하여 그 R의 대각선 원소가 10에서 15인 대각 행렬인지 살펴보기.

이 행렬의 $Q^TQ$는 무엇인지 보기

3. 원소 $u_{1.4} = 0$으로 설정하여 U의 직교성을 깨뜨리기. R은 어떻게 되며 그 이유는 무엇인지 살펴보기.

정답 코드

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

U = np.linalg.qr(A)[0]

#1
Q,R = np.linalg.qr(U)
print(np.round(R,8))
#[[ 1. -0.  0.  0.  0.  0.]
# [ 0.  1.  0. -0.  0. -0.]
# [ 0.  0.  1.  0.  0.  0.]
# [ 0.  0.  0.  1.  0. -0.]
# [ 0.  0.  0.  0.  1.  0.]
# [ 0.  0.  0.  0.  0.  1.]]

# U는 이미 직교행렬이기 때문에 U를 QR분해하면 Q는 U 자기자신이 되고, R는 단위행렬이 된다.

#2
for i in range(U.shape[0]):
  U[:,i] = U[:,i]*(10+i)
Q2, R2 = np.linalg.qr(U)
print(np.round(R2,8))
#[[10. -0.  0.  0.  0. -0.]
# [ 0. 11. -0.  0.  0. -0.]
# [ 0.  0. 12.  0.  0. -0.]
# [ 0.  0.  0. 13.  0.  0.]
# [ 0.  0.  0.  0. 14.  0.]
# [ 0.  0.  0.  0.  0. 15.]]

print(np.round(Q2.T @ Q,8))
#[[ 1. -0.  0. -0. -0.  0.]
# [-0.  1. -0. -0. -0. -0.]
# [-0. -0.  1. -0. -0. -0.]
# [ 0. -0.  0.  1.  0. -0.]
# [ 0. -0. -0.  0.  1. -0.]
# [-0.  0. -0.  0.  0.  1.]]

# QR 분해 = 벡터의 길이를 R의 대각 성분으로 뽑아내고 방향은 Q에 집어넣음.

#3
U[0,3] = 0
print(np.round(np.linalg.qr(U)[1]))
#[[10.  0. -0.  0.  0. -0.]
# [ 0. 11.  0.  1. -0.  0.]
# [ 0.  0. 12.  0.  0.  0.]
# [ 0.  0.  0. 13.  0. -1.] 
# [ 0.  0.  0.  0. 14.  0.]
# [ 0.  0.  0.  0.  0. 15.]]

# 직교가 깨지면, QR 분해는 이를 복원하기 위해 앞 열의 방향을 제거하는 과정(투영)을 수행하고, 
# 그 투영 크기가 R의 비대각 성분에 나타난다.

8-4

전통적인 역행렬 계산 방식의 수치 오차와 QR 기반의 오차를 비교하는 문제

연습문제 7-2의 코드를 복사한 다음, 행렬을 입력으로 받고 그 역행렬을 출력하는 파이썬 함수에 붙여넣는다.(oldSchoolInv)

다음으로 5x5 난수 행렬을 생성한다.

전통적인 방법과 QR 분해 방법을 사용해 그 역행렬을 계산한다.

역 추정 오차는 행렬과 계산된 역행렬의 곱으로부터 np.eye로 구한 단위행렬의 유클리드 거리로 계산한다.

결과를 막대 그래프로 만들어 두 가지 방법을 x축에, 오차에 대한 유클리드 거리를 y축에 표시한다.

코드를 여러 번 실행하고 그래프를 살펴보고, 30x30 행렬을 사용해 다시 시도한다.

정답 코드

def oldSchoolInv(A):
  n = A.shape[0]

  if not np.diff(A.shape)[0]==0:
    raise Exception('Matrix must be square.')

  # 소행렬
  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)
  return Ainv

m = 30
B = np.random.randn(m,m)

Q,R = np.linalg.qr(B)
Binv_qr = oldSchoolInv(R) @ Q.T
Binv_old = oldSchoolInv(B)

BBi_qr = B @ Binv_qr
BBi_old = B @ Binv_old

trueI = np.eye(m)
sse = [0,0]
# 각각 유클리드 거리를 구함
sse[0] = np.sqrt(np.sum((BBi_old-trueI)**2))
sse[1] = np.sqrt(np.sum((BBi_qr-trueI )**2))

plt.figure(figsize=(6,6))

plt.bar(range(2),sse,color=[.7,.7,.7])
plt.xticks(range(2),labels=['OldSchool','QR'])
plt.ylim([0,np.max(sse)*1.1])
plt.ylabel('Eucl. distance to identity')
plt.title(f'Inverse error ({m}x{m} matrix)',ha='center')
plt.savefig('Figure_08_03.png',dpi=300)
plt.show()

출력하면 아래와 같이 나온다(30x30 일 때)

5x5 행렬일 때는 oldSchool이 더 작을 때도 있었지만, 30x30 행렬이 되니 QR이 OldSchool 보다 작을 때가 훨씬 많아졌다.

즉, 이 문제는 QR 기반 역행렬이 고전 adjoint 방식보다 훨씬 정확하고 안정적이라는 것을 보여준다.

8-5

이전 연습문제의 코드를 매번 다른 난수행렬을 사용하여 실험을 100회 이상 반복하는 for문에 넣는다.

각 반복마다 오차(유클리드 거리)를 저장하고 모든 실험 결과의 평균과 모든 개별 오차를 보여주는 그래프를 그린다.

5x5 행렬과 30x30 행렬에 대해 실험을 진행한다.

전통적인 방법 대신 np.linalg.inv를 사용해 R을 반전시켜서 효과가 있는지 확인한다.

n = 30

numExprs = 100

sse = np.zeros((numExprs,2))

for expi in range(numExprs):

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

  # old school
  Ainv_old = oldSchoolInv(A)
  AAi_old  = Ainv_old@A

  # QR
  Q,R = np.linalg.qr(A)
  #Ainv_qr = oldSchoolInv(R)@Q.T # using the old-school method
  Ainv_qr = np.linalg.inv(R)@Q.T # using numpy's inv
  AAi_qr  = Ainv_qr@A

  # differences
  trueI = np.eye(n)
  sse[expi,0] = np.sqrt(np.sum((AAi_old-trueI)**2))
  sse[expi,1] = np.sqrt(np.sum((AAi_qr-trueI )**2))


# and plot
plt.figure(figsize=(6,6))

plt.plot(np.zeros(numExprs),sse[:,0],'ko')
plt.plot(np.ones(numExprs),sse[:,1],'ko')
plt.bar(range(2),np.mean(sse,axis=0),color=[.7,.7,.7])

plt.xticks(range(2),labels=['OldSchool','QR'])
plt.ylim([0,np.max(sse)*1.1])
plt.ylabel('Eucl. distance to identity')
plt.title(f'Inverse error ({n}x{n} matrix)',ha='center')
plt.savefig('Figure_08_04a.png',dpi=300)
plt.show()

30x30 행렬일 때 출력하면 아래와 같이 나온다.

QR에서 oldSchool로 R의 inverse를 구했을 때

QR에서 np.linalg.inv로 R의 inverse를 구했을 때

이 문제도 8-4처럼 QR 분해로 역행렬을 구하는 것이 훨씬 더 안정적이라는 것을 보여준다.

+ np.linalg.inv를 사용하면 더 안정적이다.

8-6

정방 직교 행렬은 모든 특잇값이 1이다. 즉, 유도된 2-노름(유도된 노름은 가장 큰 특잇값)이 1이고, 프로베니우스 노름이 M이다.

프로베니우스 노름이 M인 이유는 제곱된 특잇값을 더한 값의 제곱근이기 때문이다. 이러한 특성을 확인하는 문제이다.

난수 행렬의 QR 분해를 통해 MxN 직교 행렬을 만든다.

np.linalg.norm을 사용해 유도된 2-노름을 계산하고 5장에서 배운 방정식을 사용해 M의 제곱근으로 나누어 프로베니우스 노름을 계산한다. 두 값이 모두 1인지 확인한다.

다음으로 행렬-벡터 곱셈을 사용해 유도된 노름의 의미를 알아본다. 

무작위로 M개의 원소를 가진 열벡터 v를 생성한 다음 v와 Qv의 노름을 계산한다. 이 노름은 서로 같아야한다.

정답 코드

n = 13
Q,R = np.linalg.qr(np.random.randn(n,n))

print( np.round(np.linalg.norm(Q,2),8),               # induced 2-norm
       np.sqrt( np.sum(Q**2) )/np.sqrt(n) # 프로베니우스 노름을 n 제곱근으로 나눔
)
# 둘 다 1

# 직교행렬 Q는  길이를 절대 안 바꾸는 변환(회전/반사)이라서
# 공간에서 가장 많이 늘어난 벡터의 길이 증가량이 1이다.
# 따라서 유도된 2-노름 = 가장 큰 특잇값 = 1

# 직교행렬 Q의 각 열벡터의 길이는 1
# 프로베니우스 노름 = 각 열의 노름의 제곱을 모두 더해 루트 씌운 것
# 각 열의 노름이 모두 1이니까 프로베니우스 노름 = n의 제곱근
# 따라서 n의 제곱근으로 나누면 1

v = np.random.randn(n,1) # 랜덤한 벡터 생성
norm_v  = np.linalg.norm(v) # 벡터의 노름을 구함
norm_Qv = np.linalg.norm(Q@v) # 직교행렬 Q로 회전/반사시킨 뒤 노름을 구함

print(norm_v)
print(norm_Qv)
# 두 값이 같다
# 즉, 직교행렬 Q는 벡터의 길이를 절대 바꾸지 않는다

# 유도된 2-노름은 어떤 벡터를 곱했을 때 가장 많이 늘어나는 비율을 의미함
# 직교행렬 Q의 2-노름은 1
# 즉, 직교행렬 Q는 직교행렬은 벡터를 절대 늘리지 않고 회전만 시킨다.

수식으로도 확인 가능하다.

벡터 노름 $ \left\| v\right\| $는 $v^Tv$로 계산할 수 있다.

따라서 $ \left\|Qv\right\| = \left ( Qv \right )^TQv = v^TQ^TQv = v^Tv = \left\| v\right\|$ 이렇게 계산할 수 있다.

즉, 직교 행렬을 통해 벡터를 회전시킬 수는 있지만 벡터의 크기는 조절할 수 없다.

8-7

A가 높고 최대열계수일 때, R 의 처음 N 행은 상삼각이고, 반면에 N+1부터 M까지의 행은 0이다.

파이썬에서 10x4 난수 행렬을 사용해 이를 확인한다. 경제형 분해가 아닌 완전형 QR 분해를 사용해야 한다.

R은 정방이 아니므로 비가역이다.

(1) 처음 N개의 행을 구성하는 행렬은 정방이고 최대계수이므로 완전 역행렬을 가진다.

(2) 높은 행렬은 의사역행렬을 가진다.

두 역행렬을 모두 계산하고, R의 처음 N행의 완전 역행렬이 R의 의사역행렬의 처음 N열과 같다는 것을 확인한다.

정답 코드

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

# get R
_,R = np.linalg.qr(A,'complete')

# examine R
np.round(R,3)
# array([[ 3.461, -0.716, -0.817, -0.628],
#       [ 0.   , -2.678, -0.591,  0.165],
#       [ 0.   ,  0.   ,  2.549, -0.259],
#       [ 0.   ,  0.   ,  0.   ,  3.317], < 4행까지 상삼각이다
#       [ 0.   ,  0.   ,  0.   ,  0.   ],
#       [ 0.   ,  0.   ,  0.   ,  0.   ],
#       [ 0.   ,  0.   ,  0.   ,  0.   ],
#       [ 0.   ,  0.   ,  0.   ,  0.   ],
#       [ 0.   ,  0.   ,  0.   ,  0.   ],
#       [ 0.   ,  0.   ,  0.   ,  0.   ]])

# 처음 4행까지만 자름
Rsub = R[:4,:]

# inverses
Rsub_inv = np.linalg.inv(Rsub) # 4행까지 자른 행렬의 역행렬
Rleftinv = np.linalg.pinv(R) # R 자체의 의사역행렬

# print out both
print('Full inverse of R submatrix:')
print(np.round(Rsub_inv,3)), print(f'\n\n')

# Full inverse of R submatrix:
# [[ 0.369 -0.134  0.242 -0.13 ]
# [-0.    -0.266  0.11  -0.045]
# [-0.    -0.    -0.33   0.301]
# [ 0.     0.     0.     0.405]]

print('Left inverse of R:')
print(np.round(Rleftinv,3))

# Left inverse of R:
# [[ 0.369 -0.134  0.242 -0.13   0.     0.     0.     0.     0.     0.   ]
# [-0.    -0.266  0.11  -0.045  0.     0.     0.     0.     0.     0.   ]
# [-0.     0.    -0.33   0.301  0.     0.     0.     0.     0.     0.   ]
# [-0.    -0.     0.     0.405  0.     0.     0.     0.     0.     0.   ]]

즉, 긴 R(10×4)의 의사역행렬 pinv(R)은, 위쪽 4×4 상삼각 부분 Rsub의 역행렬을 뒤에 0을 더해서 길게 만든 것과 같다.

A는 10×4이고 rank=4이면

• Q는 10차원 공간에서 4차원 부분공간을 만들어내고
• Rsub는 그 공간에서 좌표를 변환하는 진짜 역할을 하고
• 나머지 6개 행의 0은 그냥 쓸데없는 padding