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

[chapter10] 일반 선형 모델 및 최소제곱법 : 우주를 이해하기 위한 방법

10-1

잔차는 예측 데이터와 직교한다.

앞에서 예로 든 데이터 집합으로 오차에 따른 예측 데이터의 산점도를 만든다. 그 다음 잔차와 모델 예측 데이터 사이의 내적과 상관계수를 계산한다.

정답 코드

numcourses = [13,4,12,3,14,13,12,9,11,7,13,11,9,2,5,7,10,0,9,7]
happiness  = [70,25,54,21,80,68,84,62,57,40,60,64,45,38,51,52,58,21,75,70]

# X는 Nx1 열벡터
X = np.array(numcourses,ndmin=2).T
X_leftinv = np.linalg.inv(X.T@X) @ X.T
# 최소제곱 해
beta = X_leftinv @ happiness
# 예측값
pred_happiness = X@beta
# 잔차
res = happiness-pred_happiness

# 내적
print('Dot product: ' + str(np.dot(pred_happiness,res)) )
# Dot product: 5.4569682106375694e-12
# 0과 가깝다. 즉 직교한다.

# 상관계수
print('Correlation: ' + str(np.corrcoef(pred_happiness,res)[0,1]))
# Correlation: -0.6444852258446189
# 0과 가깝다.

print(' ')


# show in a plot
plt.figure(figsize=(6,6))
plt.plot(res,pred_happiness,'ko',markersize=12)
plt.xlabel('Residual error')
plt.ylabel('Model-predicted values')
plt.title(f'r = {np.corrcoef(pred_happiness,res)[0,1]:.20f}')
plt.savefig('Figure_10_06.png',dpi=300)
plt.show()

# 그래프를 보면 점이 산개해있다.
# 즉 예측값이 커져도 잔차는 커지지도 작아지지도 않는다.
# = 선형 관계 없다 = 직교

출력하면 아래처럼 나온다.

x축은 잔차 오류, y축은 모델 예측값이다.

10-2

잔차벡터는 그 하나의 선형 가중 결합에만 직교하지 않는다.

대신 잔차벡터는 설계 행렬을 생성하는 전체 부분공간과 직교한다. 파이썬에서 이를 구현하기.

= 잔차는 X의 열공간과 직교한다

= X의 모든 열벡터와 내적이 0이다

따라서 X의 모든 열벡터와 직교인 벡터들의 집합을 구한다.

정답 코드

# X의 열공간 외 나머지 공간
# X의 열공간에 직교한 모든 방향들의 집합이다.
nullspace = null_space(X.T)


# 잔차 res가 left-null space 안에 있는지 검사
# 어떤 벡터 r이 어떤 공간의 span 안에 있으면 기저에 r을 추가해도 랭크가 안 늘어난다
# = 기존 기저 벡터들 뒤에 벡터 r을 하나 더 붙여서 새 행렬을 만든다
nullspaceAugment = np.hstack( (nullspace,res.reshape(-1,1)) )


# 랭크가 같은지 검사
print(f'dim(  N(X)    ) = {np.linalg.matrix_rank(nullspace)}')
print(f'dim( [N(X)|r] ) = {np.linalg.matrix_rank(nullspaceAugment)}')

# dim(  N(X)    ) = 19
# dim( [N(X)|r] ) = 19

10-3

QR 분래를 통해 최소제곱법 계산하기. 아래 식들을 계산하고 서로 비교한다.

• 왼쪽 역인 $(X^TX)^{-1}X^Ty$
• 역을 $R^{-1}Q^Ty$로 하는 QR
• $Q^Ty$로 증강된 행렬에서 가우스-조던 소거법 계산

세 가지 방법의 베타 변수를 출력하고, 마지막으로 QR 분해의 결과 행렬을 출력하기

정답 코드

numcourses = [13,4,12,3,14,13,12,9,11,7,13,11,9,2,5,7,10,0,9,7]
happiness  = [70,25,54,21,80,68,84,62,57,40,60,64,45,38,51,52,58,21,75,70]

# np.array(numcourses,ndmin=2).T = numcourses를 일반 파이썬 리스트에서 열벡터로 만듦
# np.ones = 행 20개
# 둘다 행이 20개니까 np.hstack하면 (20,2)모양(열 기준으로 붙임)
# 예) X =
#    [[ 1, 13],
#     [ 1,  4],
#     [ 1, 12], ...
X = np.hstack((np.ones((20,1)),np.array(numcourses,ndmin=2).T))

# 왼쪽 역
beta1 = np.linalg.inv(X.T @ X) @ X.T @ happiness

# QR 분해
Q, R = np.linalg.qr(X)
beta2 = np.linalg.inv(R) @ (Q.T @ happiness)

# 가우스-조던
tmp = (Q.T@happiness).reshape(-1,1)
Raug = np.hstack( (R,tmp) ) # augment the matrix
Raug_r = sym.Matrix(Raug).rref()[0] # this gets the matrix
# Raug_r은 아래 형태가 됨
#[1  0 | β₁]
#[0  1 | β₂]

# 마지막 열, 즉 beta를 가져오기
beta3 = np.array(Raug_r[:,-1]) # convert back to numpy

print('Betas from left-inverse: ')
print(np.round(beta1,3)), print(' ')

print('Betas from QR with inv(R): ')
print(np.round(beta2,3)), print(' ')

print('Betas from QR with back-substitution: ')
print(np.round(np.array(beta3.T).astype(float),3))

# Betas from left-inverse: 
# [23.13   3.698]
 
# Betas from QR with inv(R): 
# [23.13   3.698]
 
# Betas from QR with back-substitution: 
# [[23.13   3.698]]

# 같은 최소제곱(OLS) 회귀 문제를 세 가지 서로 다른 선형대수 방법으로 풀면 결과가 같아야 한다는 걸 확인하는 문제

# show the matrices
print('Matrix R:')
print(np.round(R,3)) # note that it's upper-triangular (as you know!)

print(' ')
print("Matrix R|Q'y:")
print(np.round(Raug,3))

print(' ')
print("Matrix RREF(R|Q'y):")
print(np.round(np.array(Raug_r).astype(float),3)) # convert to numpy floats

# Matrix R:
# [[ -4.472 -38.237]
#  [  0.     17.747]]
 
# Matrix R|Q'y:
# [[  -4.472  -38.237 -244.849]
#  [   0.      17.747   65.631]]
 
# Matrix RREF(R|Q'y):
# [[ 1.     0.    23.13 ]
# [ 0.     1.     3.698]]

10-4

이상치는 드물거나 대표적이지 않은 데이터값이다.

이 문제에서는 행복도 데이터에 이상치를 생성해 최소제곱법에 미치는 영향을 관찰한다.

데이터 벡터에서 첫 번째 관측된 데이터 점을 70에서 170으로 변경한다.(데이터 입력 오타 시뮬레이션)

그리고 최소제곱 적합도를 다시 계산하고 데이터로 그래프를 그린다.

다음으로는 마지막 점을 70에서 170으로 변경하고 똑같이 적합도 계산 및 그래프를 그린다.

그리고 시각화 자료를 만들어서 비교한다.

이상치는 결과 변수에서 동일하지만 해당 x축 값으로 인해 데이터에 대한 모델의 적합도에 미치는 영향은 상당히 다르다.

이상치의 이러한 차별적 영향을 지렛대라고 한다.

정답 코드

happiness_oops1 = [170,25,54,21,80,68,84,62,57,40,60,64,45,38,51,52,58,21,75,70]
happiness_oops2 = [70,25,54,21,80,68,84,62,57,40,60,64,45,38,51,52,58,21,75,170]

# X는 고정, y만 바뀐다.
# 즉 투영 공간 C(X)는 동일하나, y의 위치가 달라진다.
# 즉, y를 투영한 결과(beta)가 달라진다.

X = np.hstack((np.ones((20,1)),np.array(numcourses,ndmin=2).T))
X_leftinv = np.linalg.inv(X.T@X) @ X.T

_,axs = plt.subplots(1,3,figsize=(16,5))

# 정상 데이터 : 직선들이 대충 중앙을 통과함
# oops1 : x가 큰쪽(오른쪽)에서 y가 튐 > 회귀선 전체가 끌려올라감
# oops2 : x가 작은쪽(왼쪽)에서 y가 튐 > 회귀선은 크게 흔들리지 않고 잔차만 엄청 큼.
# 잔차가 아주 크면 SSE를 압도함(SSE : 오차제곱합)
for axi,y in zip(axs,[happiness,happiness_oops1,happiness_oops2]):

  # compute the best-fit parameters
  beta = X_leftinv @ y

  # predicted data
  pred_happiness = X@beta


  # plot the data and predicted values
  axi.plot(numcourses,y,'ks',markersize=15)
  axi.plot(numcourses,pred_happiness,'o-',color=[.6,.6,.6],linewidth=3,markersize=8)

  # plot the residuals (errors)
  for n,y,yHat in zip(numcourses,y,pred_happiness):
    axi.plot([n,n],[y,yHat],'--',color=[.8,.8,.8],zorder=-10)

  # make the plot look nicer
  axi.set(xlabel='Number of courses taken',ylabel='General life happiness',
          xlim=[-1,15],ylim=[0,100],xticks=range(0,15,2))
  axi.legend(['Real data','Predicted data','Residual'])
  axi.set_title(f'SSE = {np.sum((pred_happiness-y)**2):.2f}')
  


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

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

10-5

최소제곱법을 사용하여 역행렬을 계산하는 문제

방정식 $XB = Y$에서 X는 역행렬을 구할 정방 최대계수 행렬, B는 미지의 계수 행렬(역행렬이 됨), Y는 관측 데이터(단위행렬)이다.

세 가지 방법으로 B를 계산한다.

먼저 왼쪽 역 최소제곱법을 사용해서 한 번에 한 열씩 행렬을 계산한다. 이것은 for문에서 행렬 X와 Y의 각 열 사이에 최소제곱 적합도를 계산해 구한다.

다음으로 왼쪽 역행렬을 방법을 사용해 한 줄의 코드로 전체 B 행렬을 계산한다.

마지막으로 np.linalg.inv() 함수를 사용해 $X^{-1}$을 계산한다.

각 B 행렬에 X를 곱해 그림으로 보여준다. 마지막으로 역을 계산하는 이 세 가지 다른 방법이 동등한지 시험한다.(역행렬은 고유하므로 동등해야 한다.)

정답 코드

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

Y = np.eye(n)

# np.zeros_like = X랑 똑같은 모양의 0으로만 된 배열을 만들기
Xinv1 = np.zeros_like(X)

for coli in range(n):
  # β=(X^TX)^{−1}X^Ty
  # OLS 해 공식
  # X가 정사각 & full rank이면 X^{-1} @ e_coli 가 X 역행렬의 coli번째 열
  # 즉, 각 열마다 Xbeta = e_i를 최소제곱으로 푸는 beta를 구해서 Xinv1의 i번째 열로 넣음
  Xinv1[:,coli] = np.linalg.inv(X.T @ X) @ X.T @ Y[:,coli]

Xinv2 = np.linalg.inv(X.T@X) @ X.T @ Y

Xinv3 = np.linalg.inv(X)

# visualize
_,axs = plt.subplots(1,3,figsize=(10,6))

# column-wise least-squares
axs[0].imshow( Xinv1@X ,cmap='gray')
axs[0].set_title('Via column-wise LS')

# matrix-wise least-squares
axs[1].imshow( Xinv2@X ,cmap='gray' )
axs[1].set_title('Via matrix-wise LS')

# inv()
axs[2].imshow( Xinv3@X ,cmap='gray' )
axs[2].set_title('Via inv()')


# don't need the tick marks
for a in axs: a.set(xticks=[],yticks=[])

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

출력하면 아래와 같다.

세 그림이 거의 동일 → 세 방법이 같은 결과

즉 역행렬은 표준기저들을 타깃으로 한 최소제곱 문제들의 해를 모아놓은 것이다.

이 때 X가 직사각 모양이면 최소제곱(OLS)만 살아남고, 역행렬은 죽는다.

X가 직사각이어도 각 열마다 OLS 문제는 정의된다. 즉 유사역행렬은 OLS(최소제곱) 문제를 풀어서 얻은 연산자다.