[파이썬과 선형대수] 일반선형모델(General Linear Model)과 최소제곱해(Least Square Problem)

일반선형모델(General Linear Model)

우리가 흔히 말하는 일반 선형 모델은 독립변수와 종속변수를 연관시킨 방정식의 일종이다. 간단히 예를 들자면, 학생들이 다니는 학원의 개수($n$)과 일주일 동안의 공부 시간($t$)에 따른 전체 과목 평균 성적($y$)를 예측하는 모델이 있다고 가정하자.
방정식은 아래와 같을 것이다.
$$y= \beta_{1}n+ \beta_{2}t+\beta_{3}$$
학원의 개수와 공부 시간은 독립변수이고, 앞에 계수가 붙어서 종속변수 $y$를 이룬다.
이때, 새로 측정한 데이터를 입력해서 결과를 예측했을 때, 실제 $y$와 같도록 계수 $\beta$를 정해주는 것이 중요한 과제이다.

최소제곱해의 정의

앞에서 알아본 일반선형모델의 해를 구하는 방법을 최소제곱해라고 책에서는 소개하고 있다. 독립변수(위의 예에서는 $n,t$)의 쌍으로 이루어진 행렬을 $A$, 계수 $\beta$들의 쌍을 변수로 설정하면 $Ax=b$라는 행렬 방정식로 표현할 수 있다.

예를 들어 작성해 보면,

 

y	n	t
40	3	20
50	2	55
90	4	50

 

로 데이터 세트를 만들었다면,

$$90=4\beta_{1}+50\beta_{2}+\beta_{3}$$
$$40=3\beta_{1}+20\beta_{2}+\beta_{3}$$
$$50=2\beta_{1}+55\beta_{2}+\beta_{3}$$
로 표현할 수 있다.

$$\begin{bmatrix} 4&50&1 \\ 3&20&1 \\ 2&55&1\end{bmatrix} \begin{bmatrix} \beta_{1} \\ \beta_{2} \\ \beta_{3}\end{bmatrix}=\begin{bmatrix} 90 \\ 40 \\ 50\end{bmatrix}$$
로 행렬 방정식으로 변환할 수 있다.

데이터가 많아질수록 행렬은 직사각형 행렬이 된다.
그리고 지난번 [[9. 역행렬 2]]에서 알아보았듯이 직사각형 행렬도 최대 계수라면 역행렬을 구할 수 있었다.
$$Ax=b$$
$$(A^TA)^{-1}A^TAx=(A^TA)^{-1}A^T b$$
$$x=(A^TA)^{-1}A^T b$$

최소제곱해의 배경

이 책에서는 최소제곱해를 설명하는 순서를 이 글과 같게 설명해 놓았다. 물론 최소제곱해가 일반선형모델의 해를 구하는 방법이라는 것에는 동의하지만 공부하면서 실질적인 이유가 따로 있다는 생각이 들었다.

모든 입력, 데이터가 행렬 $X$의 열공간에 있다면 정확한 해($\beta$)를 찾을 수 있겠지만, 수많은 데이터의 선형 결합이 열공간에 정확히 위치하기란 불가능하다. 그렇다면 차선책은 무엇일까? 위의 선형 모델을 예로 들어보면, 얻은 선형모델($\beta$들의 세트)을 가지고 입력 데이터($n,t$)를 대입했을 때, 실제 $y$와 그 오차가 가장 적게 나오게 하는 것이 차선책이라 할 수 있다.

최소제곱해의 기하학적 관점

 

위 그림을 한번 살펴보자.
$Col; A$는 입력 데이터 행렬 $A$가 이루는 열공간이다. 실제 결과 $b$는 열공간 밖에 있다. 즉, 어떤 해($\beta$)를 집어넣더라도 $b$가 되는 것은 불가능하다.
배경 설명에서 그러면 최소 오차가 되는 해를 찾자고 했다. 이때, 최소 오차가 되는 해는 무엇일까?
그림에서 $Ax_{1},;Ax_{2},;A \hat{x}$ 세 가지의 결괏값을 볼 수 있다. 한눈에 보더라도 열공간에 수직이 되는 $A \hat{x}$이 $b$와의 거리가 가장 작다.
오차는 $b-A\hat{x}$ 로 표현할 수 있다.

그림에서 보는 것과 같이 $b-A\hat{x}$ 가 $Col ;A$와 수직이다. 평면은 행렬의 열공간이기 때문에
행렬의 모든 열과 수직이다.
$$a_{j} \cdot (b-A \hat{x})=0$$
행렬 계산으로 변환해 보면,
$$A^T(b-A\hat{x})=0$$
가 성립한다.

위의 식을 풀어서 정리해 보면
$$A^TA\hat{x}=A^Tb$$
$$\hat{x}=(A^TA)^{-1}A^Tb$$
라는 같은 결과를 얻을 수 있다.

QR분해로 구하는 최소제곱해

위에서 구한 방법은 역행렬을 계산해야 하기 때문에 수치적으로 불안정하다는 단점이 있다. 반면에 이전에 알아보았던 QR분해를 사용한다면 조금 더 간단한 방식으로 구할 수 있다.

$$Ax=b$$
$$QRx=b$$
이때, Q는 직교행렬이므로 역행렬과 전치행렬이 동일하다는 것을 알고 있다.
$$\hat{x}=R^{-1}Q^Tb$$
로 해를 구할 수 있다.

하지만 $R$의 역행렬을 구해야 한다는 점에서 위 방법과 별 차이가 없다고 느낄 수도 있다.

반대로 $R^{-1}$를 구하지 않고도 해를 구할 수 있다.
$R$이 상삼각행렬이라는 점을 생각해 보면 역치환을 통해서 빠르게 해를 구할 수 있다.
$$Rx=Q^Tb$$
방정식을 풀어버리는 것이다. 계수행렬을 $Q^{T}b$와 증강시킨 후 RREF(Reduced Row Echelon Form)을 만들어서 바로 해를 구할 수도 있다.

역행렬을 구해야 한다는 불편함을 차치하고 QR분해를 할 때, 그람-슈미트 과정을 사용하면 수치적으로 불안정해진다는 단점이 있었다. 따라서 하우스홀더 변환을 통해서 QR분해를 하고 증강행렬(Augmented Matrix)를 만든 후에 RREF를 계산, 해를 구하면 수치적으로 안정되게 구할 수 있다.

---

파이썬으로 구현하기

책에서 학생들의 수강한 수업 수에 따른 행복도 데이터를 가지고 최소제곱해를 구하는 문제를 다루고 있다.
데이터는 아래 numcourses와 happiness 변수와 동일하다.

코드

import numpy as np
import matplotlib.pyplot as plt
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 = np.hstack((np.ones((20,1)),np.array(numcourses,ndmin=2).T))

# 최소제곱해를 계산
X_leftinv = np.linalg.inv(X.T@X) @ X.T
beta = X_leftinv @ happiness

# 예측한 데이터
pred_happiness = X@beta

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

# 데이터 산점도 plotting
plt.plot(numcourses,happiness,'ks',markersize=15)
plt.xlabel('Number of courses taken')
plt.ylabel('General life happiness')
plt.xlim([-1,15])
plt.ylim([0,100])
plt.grid()
plt.xticks(range(0,15,2))

# 예측 데이터 plotting
plt.plot(numcourses,pred_happiness,'o-',color=[.6,.6,.6],linewidth=3,markersize=8)

# 데이터와 예측 사이 오차 표시(점선)
for n,y,yHat in zip(numcourses,happiness,pred_happiness):
  plt.plot([n,n],[y,yHat],'--',color=[.8,.8,.8],zorder=-10)

# 총 오차의 합 표시
plt.title(f'SSE = {np.sum((pred_happiness-happiness)**2):.2f}')

plt.show()

실행결과