[인사이드 머신러닝] 다중회귀모델 (Multiple Linear Regression Model)

cleansky·2021년 8월 15일

인사이드 머신러닝

목록 보기

5/12

앞의 두 포스트
1. [인사이드 머신러닝] 단순회귀모델: 회귀계수의 추정
2. [인사이드 머신러닝] 단순회귀모델: 회귀선의 적합도 평가
에서는 독립변수가 하나인 단순회귀모델에 대하여 알아보았다. 그러나 우리가 마주하게 되는 거의 대부분의 문제에서는 종속변수 $y$ 의 변화는 2개 이상의 독립변수에 영향을 받는 경우가 많다. 이와 같이 종속 변수의 변화를 설명하기 위해 2개 이상의 독립 변수가 사용되는 선형회귀모델을 다중선형회귀모델(multiple linear regression model) 혹은 간단하게 다중회귀모델(multiple regression model)이라고 한다.

이번 포스트에서는 다중회귀모델을 다루려고 한다. 사실 독립변수가 2개 이상이라는 점을 제외하면 다중회귀 모델의 기본적인 개념과 분석방법은 단순회귀모델과 같다. 그럼에도 불구하고 단순회귀모델과 다중회귀모델을 나누어 살펴보는 이유는 보다 단순화를 통해 통찰력을 얻고 이를 일반화하여 살펴보기 위함이다.

❗ 다중회귀 vs. 다항회귀
다중회귀와 다항회귀(polynomial regression model) 완전히 다른 용어이다.

다중회귀는 $y=a_1x_1 + a_2x_2 + \dots$ 와 같이 2개 이상의 독립변수가 있는 경우를 말하며 모든 독립변수에 대하여 선형이다.

다항회귀는 $y=a_1x + a_2x^2 + \dots$ 와 같이 하나의 독립변수와 종속변수의 관계가 고차 다항식으로 표현되는 경우를 말한다. 이때, 당연히 $x$ 와 $y$ 는 비선형 관계이다. 다항회귀는 다음 포스트에서 살펴보기로 한다.

다중회귀모델

모델의 정의

앞에서 설명한대로 다중회귀모델은 2개 이상의 독립변수가 있는 다음과 같은 모델을 의미한다,

\tag{1} y_i = \beta_0 + \beta_{1}x_{i1} + \beta_{2} x_{i2} + \dots + \beta_{k}x_{ik} + \epsilon_{i}

여기서 아래첨자 $y_i$ 는 $i$ 번째 관측치를 의미하고, $\epsilon_i$ 는 이때의 오차항을 나타낸다. 우리가 추정하고자하는 값, 즉 회귀계수는 $\beta_{j}$ ( $0 \leq j \leq k$ )이고, 독립변수 $x_{ij}$ 는 known value이다. 식 (1)을 $N$ 개의 샘플에 대하여 확장한 후, vector-matrix 형태로 표기하면 다음과 같다.

\tag{2} \underbrace{ \begin{bmatrix} y_1\\ y_2\\ \vdots\\ y_N \end{bmatrix}}_{\mathbf{y}}= \underbrace{ \begin{bmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1k}\\ 1 & x_{21} & x_{22} & \cdots & x_{2k}\\ \vdots & \vdots & \vdots & \ddots & \vdots &\\ 1 & x_{N1} & x_{N2} & \cdots & x_{Nk}\\ \end{bmatrix}}_{\mathbf{X}} \underbrace{ \begin{bmatrix} \beta_0 \\ \beta_1 \\ \vdots \\ \beta_k \end{bmatrix}}_{\bm{\beta}}+ \underbrace{ \begin{bmatrix} \epsilon_1\\ \epsilon_2\\ \vdots\\ \epsilon_N \end{bmatrix}}_{\bm{\epsilon}}

여기서 $\mathbf{e} \sim \mathcal{N}(0,\,\sigma^{2}\mathbf{I}_N)$ .

회귀선의 추정

단순회귀모델에서와 마찬가지로 다중회귀모델에서도 회귀계수들로 이루어진 가중치 벡터 $\bm{\beta}$ 는 오차제곱합(Sum of Squared Error: SSE)을 최소화하는 다음과 같은 최적화 문제의 해로 얻어진다.

\tag{3} \hat{\bm{\beta}}=\argmin_{\bm{\beta}} (\mathbf{y}-\mathbf{X}\bm{\beta})^T(\mathbf{y}-\mathbf{X}\bm{\beta})

위 식의 목적함수는 SSE이며, 이를 $S$ 로 표기하여 전개하면 다음과 같다.

\tag{4} \begin{aligned} S &= (\mathbf{y}-\mathbf{X}\bm{\beta})^T(\mathbf{y}-\mathbf{X}\bm{\beta})\\ &= \mathbf{y}^T\mathbf{y} - \mathbf{y}^T\mathbf{X}\bm{\beta}-\bf{\beta}^T\mathbf{X}^T\mathbf{y} + \bm{\beta}^T\mathbf{X}^T\mathbf{X}\bm{\beta}\\ &= \mathbf{y}^T\mathbf{y} -2\bf{\beta}^T\mathbf{X}^T\mathbf{y} + \bm{\beta}^T\mathbf{X}^T\mathbf{X}\bm{\beta} \end{aligned}

식 (4)의 두 번째 줄에서 $\mathbf{y}^T\mathbf{X}\bm{\beta}$ 와 $\bf{\beta}^T\mathbf{X}^T\mathbf{y}$ 는 서로 transpose 관계이면서 스칼라이므로 둘은 같다. 따라서 식 (4)의 세번째 줄이 성립한다. $S$ 를 최소화하는 $\bm{\beta}$ 를 찾기 위해 $S$ 를 $\bm{\beta}$ 로 미분한다. 이때, $S$ 는 스칼라이고 $\bm{\beta}$ 는 벡터이기 때문에 미분 결과는 벡터가 된다. 따라서 미분 결과가 영벡터가 되는 $\bm{\beta}$ 를 찾으면 그것이 식 (3)의 해가 된다. $S$ 를 $\bm{\beta}$ 로 미분하면 다음과 같다.

\tag{5} \frac{\partial S}{\partial \bm{\beta}} = -2\mathbf{X}^T\mathbf{y} + 2\mathbf{X}^T\mathbf{X}\bm{\beta}

따라서 식 (3)의 해는 다음과 같다.

\tag{6} \hat{\bm{\beta}}=\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}

식 (2)에 의해 $E[\mathbf{y}]=\mathbf{X}\bm{\beta}$ 이므로, $E[\hat{\bm{\beta}}]=\bm{\beta}$ 이므로 $\hat{\bm{\beta}}$ 은 $\bm{\beta}$ 의 불편추정량(unbiased estimate)이다 1.식 (6)을 이용하여 종속변수 $\mathbf{y}$ 의 추정량은 다음과 같이 표현할 수 있다.

\tag{7} \begin{aligned} \hat{\mathbf{y}} &= \mathbf{X}\hat{\bm{\beta}} \\ &= \mathbf{X}\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}\\ &= \mathbf{H} \mathbf{y} \end{aligned}

여기서 $\mathbf{H}$ 는 projection matrix, influence matrix 또는 hat matrix라 부르며 idempotent matrix이다. 즉, $\mathbf{H}^2=\mathbf{H}\mathbf{H}=\mathbf{H}$ 이면서 $\mathbf{H}^T\mathbf{H}$ 이다. 식 (7)의 결과를 이용하여 오차 벡터는 다음과 같이 나타낼 수 있다.

\tag{8} \begin{aligned} \mathbf{e} &= \mathbf{y}-\hat{\mathbf{y}} \\ &= \mathbf{y} - \mathbf{H} \mathbf{y} \\ &= (\mathbf{I}_N-\mathbf{H}) \mathbf{y} \end{aligned}

Residual의 성질

다중회귀모델은 단순회귀모델의 확장이므로 단순회귀모델: 회귀계수의 추정 포스트에서 살펴본 residual의 성질을 그대로 따른다. 즉,

$\mathbf{e}^T\mathbf{e}=0$ (잔차의 합은 0)
$\hat{\mathbf{y}}^T\mathbf{e}=0$ (종속변수와 오차는 직교)
$\mathbf{X}^T\mathbf{e}=\mathbf{0}$ (독립변수와 오차는 직교)
$\mathbf{X}^T\mathbf{e}= \begin{bmatrix} 1 & 1 & \dots & 1 \\ x_{11} & x_{21} & \dots & x_{N1} \\ x_{11} & x_{22} & \dots & x_{N2} \\ \vdots & \vdots & \ddots & \vdots \\ x_{1k} & x_{2k} & \dots & x_{Nk} \\ \end{bmatrix} \begin{bmatrix} e_1\\ e_2\\ \vdots\\ e_N \end{bmatrix}= \begin{bmatrix} \sum e_i \\ \sum x_{i1}e_i\\ \sum x_{i2}e_i\\ \vdots \\ \sum x_{ik} e_i \end{bmatrix}= \begin{bmatrix} 0 \\ 0 \\ 0 \\ \vdots \\ 0 \end{bmatrix}$
오차항 $\epsilon_i$ 는 서로 상관관계가 없더라도 보통 residual $e_i$ 들 간에는 상관관계가 존재한다. 즉,
$\tag{9} \begin{aligned} Cov[\mathbf{e}\mathbf{e}^T] &= (\mathbf{I}_N-\mathbf{H}) Cov[\mathbf{y}] (\mathbf{I}_N-\mathbf{H})^T \\ &= (\mathbf{I}_N-\mathbf{H})(\mathbf{I}_N-\mathbf{H})^T\sigma^2\\ &= (\mathbf{I}_N-\mathbf{H})\sigma^2 \end{aligned}$
여기서 다음의 성질을 이용하였다.
$\begin{aligned} (\mathbf{I}_N-\mathbf{H})(\mathbf{I}_N-\mathbf{H})^T &= \mathbf{I}_N-\mathbf{H}^T-\mathbf{H}+\mathbf{H}\mathbf{H}^T=\mathbf{I}_N-\mathbf{H} \\ E[\mathbf{y}] &= \mathbf{X}\bm{\beta} \\ Cov[\mathbf{y}] &= E[(\mathbf{y}-\mathbf{X}\bm{\beta})(\mathbf{y}-\mathbf{X}\bm{\beta})^T] = E[\mathbf{e}\mathbf{e}^T] = \sigma^2\mathbf{I}_N \end{aligned}$
식 (9)에서와 같이 $\mathbf{I}_N-\mathbf{H}$ 이 대각행렬이 아니므로 residual 간의 공분산이 존재한다2.

회귀직선의 적합도 평가

다중회귀모델에서도 단순회귀모델에서와 같이 분산분석에 의한 $F$ 검정, 결정계수, MSE 등이 회귀직선의 적합도 평가에 사용된다. 따라서 기본적인 개념은 이전 포스트를 참고하면 된다.

본 절에서는 이전 포스트의 내용을 다중회귀모델에 보다 손쉽게 적용하기 위해 분산분석에 사용되는 $SST$ , $SSR$ , $SSE$ 의 표현식을 살펴보고 수정결정계수(adjusted R-squared)에 대하여 알아본다.

분산분석

총 제곱합, $SST$ 는 다음과 같다.

\tag{10} \begin{aligned} SST &= \sum_{i=1}^N (y_i-\bar{y})^2 \\ &=\sum_{i=1}^N y_i^2 - N\bar{y}^2 \\ &= \mathbf{y}^T\mathbf{y} - N \bar{y}^2 \\ &= \mathbf{y}^T\mathbf{y} - \frac{(\sum_{i=1}^N y_i)^2}{N} \\ &= \mathbf{y}^T\mathbf{y} - \mathbf{y}^T\mathbf{1}\mathbf{1}^T\mathbf{y}/N \\ &= \mathbf{y}^T\mathbf{y} - \mathbf{y}^T(\mathbf{J}/N)\mathbf{y} \\ &= \mathbf{y}^T\left(\mathbf{I}_N-\mathbf{J}/N\right)\mathbf{y} \end{aligned}

여기서 $\mathbf{1}$ 은 모든 요소가 $1$ 인 열벡터를 의미하고, $\mathbf{J}$ 는 모든 요소가 $1$ 인 $N \times N$ 행렬을 의미한다. $SST$ 의 자유도는 $N-1$ 이다.

오차에 의한 제곱합, $SSE$ 는 다음과 같다.

\tag{11} \begin{aligned} SSE &= \sum_{i=1}^N (y_i - \hat{y}_i)^2 \\ &= \mathbf{e}^T\mathbf{e} \\ &= \mathbf{y}^T (\mathbf{I}_N-\mathbf{H}) \mathbf{y} \end{aligned}

$SSE$ 에 대한 자유도는 $N-k-1$ 이다.

회귀에 의한 제곱합, $SSR$ 은 다음과 같다.

\tag{12} \begin{aligned} SSE &= \sum_{i=1}^N (\hat{y}_i-\bar{y})^2 \\ &=\hat{\mathbf{y}}^T\hat{\mathbf{y}}-N\bar{y}^2 \\ &=\mathbf{y}^T\mathbf{H}\mathbf{y} - \mathbf{y}^T(\mathbf{J}/N)\mathbf{y} \\ &=\mathbf{y}^T\left(\mathbf{H}-\mathbf{J}/N\right)\mathbf{y} \end{aligned}

여기서 두 번째 줄의 유도하는데 있어 $E[\hat{\mathbf{y}}]=E[\mathbf{X}\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T\mathbf{y}]=E[\mathbf{X}\left(\mathbf{X}^T\mathbf{X}\right)^{-1}\mathbf{X}^T(\mathbf{X}\bm{\beta}+\bm{\epsilon})]=\mathbf{X}\bm{\beta}=E[\mathbf{y}]$ 의 성질을 사용하였고3, 세 번째 줄을 유도하는데 있어 식 (7)과 식 (10)의 결과를 이용하였다. $SSR$ 의 자유도는 독립변수의 수와 같은 $k$ 가 된다.

$SST$ , $SSE$ , $SSR$ 을 계산할 수 있으므로 이로부터 단순회귀모델에서와 같이 분산분석표와 $F$ 검정통계량을 유도하여 검정을 실시한다.

수정결정계수

다중회귀모델에서 독립변수가 추가되면 결정계수는 항상 증가한다. 왜냐하면 $SST$ 는 고정된 값이고, 독립변수가 추가될수록 오버피팅에 의해 $SSE$ 는 작아지기 때문이다. 따라서 단순회귀모델에서 사용했던 결정계수를 그대로 사용할 수 없으며 다음과 같이 수정결정계수를 이용한다.

\tag{13} R_{\text{adjusted}}^2 = 1-\frac{SSE/(n-k-1)}{SST(n-1)} = 1-\left(1-R^2\right)\frac{n-1}{n-k-1}

수정결정계수는 식 (1)의 정의에서와 같이 독립변수가 추가되더라도 $SST$ 와 $SSE$ 를 각각의 자유도로 정규화하기 때문에 설명력이 떨어지는 독립변수가 추가될 때는 감소하는 성질을 갖는다. 참고로 단순회귀모델에서는 $k=1$ 이므로 결정계수와 수정결정계수가 항상 같게 된다.

다중회귀와 표준화

두 개의 독립변수가 존재하는 간단한 다중회귀모델을 생각해보자. 여기서 종속변수는 행복지수라고 하고, 독립변수는 월평균소득과 월평균근로시간이라고 해보자. 아마도 소득은 많을수록 그리고 근로시간은 짧을수록 행복지수가 높게 나오는 경향이 있을 것이다. 그리고 실제 빅데이터를 기반으로 회귀계수를 추정해보니 다음과 같은 회귀모델을 찾았다고 해보자.

\text{행복지수} = 0.1\times\text{월평균소득} - 100\times\text{월평균근로시간} + 0.05

위 회귀모델을 보면 월평균소득에는 $0.1$ 라는 값이 곱해졌고 월평균근로시간에는 $-100$ 이라는 가중치가 곱해졌다. 월평균근로시간에 곱해진 가중치의 크기가 월평균소득에 곱해진 가중치의 1000배나 된다. 그렇다면 월평균근로시간에 곱해진 가중치의 절대값이 더 크기 때문에 행복지수는 월평균소득보다 월평균근로시간에 더 많은 영향을 받는다고 할 수 있을까? 그렇지 않다. 왜냐하면 월평균소득은 값의 범위가 보통 수백만 이상인데 반해 월평균근로시간은 많아야 수백이기 때문이다. 가령 월평균소득이 백만원이라면 0.1을 곱하더라도 10만인데 반해, 월평균 근로시간이 200시간이라고 하더라도 100을 곱해봐야 2만 밖에 되지 않는다. 따라서 다중회귀모델에서는 회귀계수의 크기만을 갖고 독립변수의 영향력을 설명할 수 없다.

지금까지 우리는 회귀모델을 적합하는데 있어 각 독립변수를 날 것 그대로 사용하였다. 이러한 회귀모델을 비표준화 회귀모델이라고 한다. 반면, 각 독립변수를 평균이 $0$ 이고 표준편차가 $1$ 이 되도록 표준화한 후 회귀계수를 추정하는 모델을 표준화된 회귀모델(standardized regression model)이라 한다. 즉, 데이터를 이용하여 각 독립변수 $x_k$ 의 표본평균 $\mu_k$ 와 표본표준편차 $\sigma_k$ 를 구한 후, 다음과 같이 $x_{ik}$ 에 z-score 변환을 적용한다.

\tag{14} x_{ij} = \frac{x - \mu_k}{\sigma_k}

참고로 종속변수도 반드시 표준화를 하여야 한다. 각 독립변수와 종속변수를 표준화한 후 회귀계수를 추정하면, 회귀계수의 절대값 비교를 통해 독립변수의 영향력을 판단할 수 있다.

계산 복잡도 문제

식 (6)과 같이 정규방적식(normal equation)을 통해 $\hat{\bm{\beta}}$ 의 해를 구하면 최적의 해를 구할 수 있지만 $\mathbf{X}^T\mathbf{X}$ 의 역행렬을 계산하여야 한다. 행렬 $\mathbf{X}$ 의 크기는 $N \times k+1$ 이므로 $\mathbf{X}^T\mathbf{X}$ 를 계산하는데 걸리는 시간은 샘플 수 $N$ 에 선형적으로 증가한다. 반면, 크기는 $(k+1) \times (k+1)$ 인 $\mathbf{X}^T\mathbf{X}$ 의 역행렬을 구하기 위한 계산 복잡도는 알고리즘에 따라 $O(k^{2.4})$ ~ $O(k^{3})$ 이므로 독립변수의 수가 증가할 수록 역행렬을 계산하는데 어려움이 발생한다. 가령 독립변수의 수가 2배 늘어난다면 역행렬을 구하는 시간은 대략 8배 정도 늘어난다.

따라서 특성의 수가 많거나( $k$ 가 크거나), 훈련 샘플이 너무 많을 경우( $N$ 이 클 경우)에는 식 (6)을 그대로 사용할 수 없다. 게다가 훈련 샘플이 추가된 경우에는 매번 역행렬을 다시 구하여야 하므로 굉장히 비효율적이라고 할 수 있다. 이 경우, 우리는 온라인 학습4을 위해 경사하강법(gradient descent, GD) 알고리즘을 사용할 수 있다.

경사하강법을 이용한 회귀계수의 추정

text [그림 1] 경사하강법이 최적의 해를 찾아가는 과정

경사하강법은 다양한 최적화 문제의 해를 구하는데 널리 사용되는 알고리즘이다. 경사하강법의 기본적인 아이디어는 간단하다. 우리가 원하는 것은 오차를 최소화하는 답을 구하는 것이다. 오차를 최소화한다는 것은 식 (5)와 같이 목적함수의 미분이 0이 되는 지점을 찾는 것과 같다. 식 (6)은 미분값이 0이 되는 해를 한 번에 찾아주지만 앞서 설명한 바와 같이 경우에 따라 복잡도 문제로 인해 식 (6)과 같은 방법으로 해를 찾는데는 제약이 따른다. 반면, 경사하강법은 식 (6)의 해와 같거나 유사한 값의 해를 점진적으로 찾아준다5. 경사하강법은 아래 수식과 같이 동작하며, [그림 1]과 같이 시각화할 수 있다.

\tag{15} \hat{\bm{\beta}}_{i+1} = \hat{\bm{\beta}}_{i} -\eta \nabla_{\bm{\beta}}S

여기서 $i$ 는 샘플 인덱스를 의미한다. 즉, 경사하강법은 $i$ 번째 샘플을 이용하여 추정된 $\hat{\bm{\beta}}_i$ 를 $i+1$ 번째 샘플을 이용하여 $\hat{\bm{\beta}}_{i+1}$ 로 갱신한다. 식 (15)에서 $\eta$ 는 학습률을 의미하며 학습 전에 설계자가 임의로 설정하여야 하는 값이다. $\eta$ 가 너무 크면 최소 값이 되는 지점을 스쳐지나가게 되어 발산하거나 최적의 해에 잘 수렴할 수 없다는 단점이 있지만 대신 학습 속도는 빠르다. 반면, $\eta$ 가 너무 작으면 학습 속도는 오래 걸리지만 충분한 iteration을 통해 최적의 해에 도달할 수 있다. 최적화 문제를 다루는데 있어 식 (15)의 경사하강법을 그대로 사용하기도 하지만 학습률을 적응적으로 변화시키거나 과거의 학습 결과와 현재의 학습 결과에 적당한 가중치를 두는 방법 등 다양한 방법의 변형이 존재하며, 이 문제는 나중에 별도의 포스트를 통해 다루기로 한다. 경사하강법은 최근 매우 핫한 딥러닝의 기본이 되는 알고리즘이기도 하다.

경사하강법은 독립변수(즉, feature)의 수가 많거나 훈련 데이터의 크기가 크더라도 사용하는데 무리가 없으며, 새로운 훈련 데이터가 주어진 경우에도 기존에 학습한 결과에 새로운 훈련 데이터에 의한 학습 결과를 반영해주면 된다는 장점이 있다. 그러나 경사하강법에 의해 찾은 해가 항상 최적해라는 보장은 없다. 예를 들어, [그림 2]와 같이 초기 값에서 경사하강법을 통해 미분이 0이 되는 지점을 찾았음에도 불구하고 최소 값에는 도달할 수 없게 되는 경우가 발생한다. 다행히도 선형회귀에서 비용함수로 다루는 SSE(혹은 MSE)는 convex 함수이기 때문에 항상 global minimum에 도달하는 것이 보장되지만 많은 문제에서 비용함수는 convex나 concave하지 않기 때문에 경사하강법은 local minimum(혹은 local maximum)만 찾을 수 있는 경우가 대부분이다.

text [그림 2] 경사하강법의 문제점

경사하강법을 사용할 때 주의하여야 할 점은 앞에서 살펴본 표준화 등을 통해 각 독립변수들이 반드시 같은 스케일을 갖도록 해주어야 한다는 점이다. 만약, 그렇지 않으면 gradient 값이 특정 독립변수에 지나치게 영향을 많이 받게 되어 수렴하는데 훨씬 오랜 시간이 걸리게 된다.

text [그림 3] 특성 스케일을 적용한 경우의 수렴성능(좌측)과 그렇지 않은 경우의 수렴성능(우측) [6]

그리고 경사하강법은 독립변수가 많은 경우에 식 (6)의 정규방정식 대비 복잡도 면에서 유리하기는 하지만 최적의 해에 수렴하는데는 역시 많은 시간이 걸리게 된다(최적의 해를 찾는데 비현실적으로 오랜 시간이 필요한 경우, 최적의 해에 도달하지 못한다.). 머신러닝 분야에서 이러한 문제를 차원의 저주라고 하는데 너무 잘 알려진 문제이기 때문에 별도의 언급은 생략한다. 마찬가지로 경사하강법을 사용하는데 있어 배치, 미니배치, 스토캐스틱 방법을 사용할 수 있으며, 이 내용 역시 별도의 언급은 생략한다. (나중에 별도로 포스팅할지는 모르겠지만 구글에 너무나도 많은 자료가 있다.)

마지막으로 선형회귀 뿐만 아니라 딥러닝 등에서도 경사하강법을 통해 $\epsilon$ 의 오차 범위 안에서 최적이 해를 찾겟다고 할 경우 $O(1/\epsilon)$ 의 epoch(=iteration)이 걸릴 수 있다는 점에 유의하여야 한다. 따라서 오차 범위를 절반으로 줄이기 위해서는 최소 2배 이상의 반복 학습이 필요하며, 10배 줄이고 싶을 경우 10배 이상의 반복 학습이 필요하다.

파이썬 예제

아래 예제는 특성(독립변수) 수가 2개인 다항회귀 모델의 최적해를 정규방정식과 경사하강법을 통해 구하는 예제코드이다. 추정된 회귀 계수를 출력해보면 두 방식 모두 거의 동일한 회귀 계수 값을 추정한 것을 확인할 수 있을 것이다. 특성수와 학습률에 따라 epoch 수는 더 늘려야 할 수도 있다.

아래 예제에서는 경사하강법을 사용할 때, 훈련 데이터에서 랜덤하게 하나의 샘플을 선택하여 학습하는 확률적 경사하강법 (stochastic gradient descent) 방법을 사용하였다. 이 방법은 최적의 해를 찾아가는 과정에서 배치나 미니배치보다 learning curve가 더 요동치지만 무작위 선택에 의해 지역해를 건너뛸 수 있는 가능성이 증가하기 때문에 전역해를 찾을 가능성이 더 높다.

n_feats = 2

# generage regression dataset
from sklearn.datasets import make_regression
x, y = make_regression(n_samples=100, n_features=n_feats, n_targets=1, bias=2, noise=10, random_state=42)

# normal equation
from sklearn.linear_model import LinearRegression

# train medel (least sqaure estimation)
lr = LinearRegression()
lr.fit(x, y)

print(lr.coef_, lr.intercept_)

# stochastic gradient descent
import numpy as np

n_data = x.shape[0]

n_epochs = 1000
eta = 0.01

beta = np.random.randn(n_feats+1, 1) # 무작위 초기화, +1 for bias
for epoch in range(n_epochs):
    idx = np.random.randint(n_data)
    xi = np.reshape(np.append(x[idx,], np.reshape([1], (1, 1)) ), (1,-1))
    yi = y[idx]

    # 식 (5)에서 X행 하나와 그에 대응하는 y에 대하여 계산, 상수 2는 생략
    gradients = (xi.T@xi)@beta   - xi.T*yi     
    
    beta = beta - eta*gradients

print(beta)

Footnotes

1: $\theta$ 의 추정량 $\hat{\theta}$ 의 기대값이 $\theta$ 와 같도록 하는 estimator를 unbiased estimator이라 하고, 무수히 많은 unbiased estimator들 중에서 오차의 variance가 가장 작은 estimator를 efficient estimator라 한다. 이때, 이론적으로 유도된 minimum variance를 Cramer-Rao lower bound (CRLB)라 한다.

2: (넘어가도 되는 내용임. 개인적인 메모) 무선통신에서는 $\mathbf{X}$ 가 특별한 구조가 되도록 CAZAC 시퀀스 등을 사용하여 채널추정 혹은 동기 등에 활용한다. 이 경우, $\mathbf{H}$ 는 단위랭렬이 되며, residual은 항상 0이다. 혼동하지 말아야 할 것은 residual이 0이 perfect estimation을 의미하지는 않는다. 가령 채널추정에서 residual 0의 의미는 수신 신호를 CAZAC 시퀀스와 추정된 채널로 완벽히 복원이 가능하다는 의미이다. 채널추정은 항상 (SNR에 반비례하는) 추정오차가 존재한다.

3: 식 (12)에서는 sample mean을 사용하지만 샘플 수가 많으면 sample mean이 ensemble mean에 수렴하므로 식 (12)의두 번째 줄은 유효하다.

4: 모든 샘플을 한 번에 이용하여 학습하는 방식을 오프라인 학습이라 하고, 하나의 샘플(혹은 일부 샘플)에 대한 학습결과를 이용하여 가중치를 계속해서 갱신해나가는 학습방법을 온라인 학습이라 한다.

5: 단순회귀모델의 경우 훈련 데이터의 수가 많더라도 독립변수가 하나이기 때문에 경사하강법 대신 정규방정식을 이용하는 것이 더 유리하다.

참고자료

[1] 김성수, 강명욱, 강위창, "회귀모형," 한국방송통신대학교출판문화원, 2019.
[2] 장철원, "선형대수와 통계학으로 배우는 머신러닝 with 파이썬: 최적화 개념부터 텐서플로를 활용한 딥러닝까지," 비제이퍼블릭, 2021.
[3] A. Geron , "Hands-On Machine Learning with Scikit-Learn, Keras, and TensorFlow: Concepts, Tools, and Techniques to Build Intelligent Systems," 2nd Edition, O'Reilly Media, 2019.
[4] https://en.wikipedia.org/wiki/Projection_matrix
[5] https://en.wikipedia.org/wiki/Idempotent_matrix
[6] https://www.blog.nipunarora.net/ml_multi_variate_linear_regression/
[7] Kay, Steven M. Fundamentals of statistical signal processing: estimation theory. Prentice-Hall, Inc., 1993.