[ Tensorflow 2 / sklearn ] Linear Regression과 Gradient Descent Algorithm의 종류

이기는 것이 중요한 것은 아니지만,

이기기를 원하는 것은 중요하다

- 에릭 슈미트(Google 전 Ceo) -

 이번 포스팅에서는 머신러닝의 가장 간단한 모델 중 하나인 회귀(Regression)에 속하는 선형 회귀(Linear Regression)에 대해 알아보고, 이를 최적화하는 훈련 방식인 경사 하강 알고리즘과 Normal Equation에 대해 공부해보겠습니다.

---

보통 우리가 Regression에서 모델을 훈련시킨다고 하면 경사 하강 알고리즘(Gradient Descent algorithm)을 사용합니다. 좀 더 세부적으로는 Batch GD, Mini-Batch GD, SGD가 있는데, 이는 잠시 후에 다시 다뤄보겠습니다.

또 다른 훈련 방식에는 Normal Equation가 있습니다. 이 또한 Linear Regression을 다룬 후 다시 언급드리겠습니다.

---

4.1 Linear Regression

$$ 삶의 만족도 = θ_0+θ_1*(1인당 GDP) $$

 여기 1인당 GDP를 이용하여 삶의 만족도를 예측하고자 하는 Regression 모델이 있습니다. 이 모델은 Linear model이고, $θ_0$, $θ_1$이 모델 파라미터입니다.

$$y = θ_0 + θ_1*x_1 + θ_2*x_2 + θ_3*x_3 + ... + θ_n*x_n$$

 앞에서 본 것이 예시였다면 위 식은 좀 더 일반적인 Linear model을 표현한 것입니다. Linear model은 Input의 feature의 Weight 합과 bias 혹은 intercept이라 부르는 상수를 더해 최종 예측값을 만들어 냅니다.

$$y = h_θ(x) = θ*x$$

앞에서 본 일반적인 식을 벡터를 이용하여 위처럼 더욱 간단히 표현할 수 있습니다.

- $θ$는 $θ_0 ... θ_n$까지 Weight를 담은 Paraeter 벡터입니다.

- $x$는 $x_0 ... x_n$까지 Input의 feature를 담은 feature 벡터입니다. 이때 $x_0$는 항상 0입니다.

- $h_θ(x)$는 가설 함수입니다.

---

 이제 Linear Regression 모델의 구조에 대해 알아보았으니 본격적으로 훈련을 시켜보겠습니다. 모델을 훈련시킨다는 것은 Training set에 최적화되게끔 parameter를 설정한다는 의미입니다.

 그리고 이를 위해 필요한 것은 예측값과 실제값의 차이를 측정할 loss function을 정의하는 것입니다. 보통 Linear Regression 모델에서 사용하는 loss function은 MSE입니다.

---

4.1.1 Normal Equation

 Normal Equation은 앞서 언급했듯이 model parameter를 최적화하는 하나의 방식입니다. 우리가 보통 또 하나의 방법인 Gradient Descent Algorithm을 떠올리며 Gradient를 반복적으로 계산하며 조정하는 원리를 떠올립니다.

Normal Equation은 이러한 반복적인 연산없이 공식을 이용해 한번의 최적의 parameter를 도출해냅니다.

- $θ$ : loss function을 최소화하는 즉 모델에 최적화된 parameter입니다.

- $X$ : Input입니다. 이떄 $X^T$는 $X$의 전치 행렬입니다.

- $y$는 Input data에 대한 Label을 담고 있는 벡터입니다.

우선 Normal Equation을 적용하기 위해 랜덤으로 Linear한 데이터를 생성해보겠습니다.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

X = 2 * np.random.rand(100,1)
y = 4+ 3*X + np.random.randn(100,1)

plt.plot(X, y, "b.")
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$y$", rotation=0, fontsize=18)
plt.axis([0, 2, 0, 15])

plt.show()

 이제 여기에서 생성된 $x_1$과 $y$를 Normal Equation에 대입하여 최적화된 $θ$를 구해보겠습니다.

 아래 코드는 첫번째 줄에서 모든 데이터에 $x_0 = 1$을 추가한 후 np.linalg의 inv() 함수를 이용하여 역행렬(Inverse Matrix)를 만들고, dot() 메서드를 이용하여 Matrix 곱셈을 진행하였습니다.

X_b = np.c_[np.ones((100,1)), X] # np.c_ : 배열 붙이기, np.ones : array 생성 함수
theta_best = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y) #dot() : 행렬 곱셈

 Normal Equation에 의해 최적화된 θ를 확인해보겠습니다.

theta_best

##결과값##
>>
array([[3.93328217],
       [3.08032243]])

사실 $θ_0 = 4$, $θ_1 = 3$을 기대했기 때문에 완전히 정확하다고 할 순 없겠지만, 매우 비슷합니다. 이제 최적화된 $θ$를 사용하여 예측을 해보겠습니다.

X_new = np.array([[0], [2]]) # 2차원 열 벡터
X_new_b = np.c_[np.ones((2,1)), X_new] # 모든 샘플에 x0 = 1을 추가
y_predict = X_new_b.dot(theta_best)
y_predict

##결과값##
>>
array([[ 3.93328217],
       [10.09392703]])

 모델의 예측을 그래프에 나타내보겠습니다.

plt.plot(X_new, y_predict, "r-")
plt.plot(X, y, "b.")
plt.axis([0,2,0,15])
plt.show()

---

 이번엔 잠시 Normal Equation은 접어두고, sklearn에서 Lienar Regression을 사용해보겠습니다. 보시다시피 매우 간단합니다.

from sklearn.linear_model import LinearRegression
lin_reg = LinearRegression()
lin_reg.fit(X, y)

lin_reg.intercept_, lin_reg.coef_ # .intercept_:추정된 상수항, .coef_ : 추정된 weight 벡터
>>
(array([4.27746389]), array([[2.91097464]]))

lin_reg.predict(X_new)
>>
array([[ 4.27746389],
       [10.09941316]])

 sklearn의 LinearRegression 클래스는 scipy.linalg.lstsq()함수를 기반으로 하고, 이를 직접 호출할 수도 있습니다.

theta_best_svd, residuals, rank, s = np.linalg.lstsq(X_b, y, rcond=1e-6)
theta_best_svd

##결과값##
array([[4.27746389],
       [2.91097464]])

위에서 사용한 np.linalg.lstsq는 $θ=X^+*y$을 계산하여 최적화 된 $θ$를 도출합니다. 이때 $X^+$는 $X$의 유사 역행렬(pseudoinverse matrix)입니다. 이를 pseudoinverse matrix을 직접 구할 수도 있긴 합니다.

np.linalg.pinv(X_b).dot(y)

##결과값##
>>
array([[4.27746389],
       [2.91097464]])

---

 여긴 수학적인 부분이므로 건너뛰어도 무방합니다.

 pseudoinverse matrix는 특이값 분해(SVD)라는 Matrix 분해 기법을 사용하여 계산됩니다. SVD는 아래와 같이 Training set를 담은 Matrix $X$를 분해합니다.

Matrix A가 위로 길쭉한 직사각형 Matrix라면 아래와 같을 것입니다.

pseudoinverse matrix는 이 Matrix A를 아래와 같이 변형하여 도출되며, 그 형태는 아래와 같을 것입니다.

---

그렇다면 왜 np.linalg.lstsq는Normal Equation이 아닌 pseudoinverse matrix를 이용하여 parameter를 최적화하는 것일까요?

 그 이유는 Normal Equation에서는 $X^T$를 필요로 하는데, Matrix의 행이 열보다 작거나(m < n), feature가 중복되는 Matrix는 Inverse Matrix를 가지지 않는 것이 가장 큰 이유입니다. Inverse Matrix와 달리 pseudoinverse matrix는 항상 구할 수 있습니다.

---

4.1.2 계산 복잡도

 우리가 자료구조에서 빅오를 구하며, 메모리 사용의 효율성을 따지듯 머신러닝에서도 계산 복잡도를 이용하여 모델의 성능을 비교하는 것 또한 매우 중요합니다.

 우선 Normal Equation은 $X^T*X$를 계산하는데 이는 $(n+1) * (n+1)$의 크기를 가집니다. 이러한 Inverse Matrix를 구하는 계산 복잡도는 $O(n^2.4)$에서 $O(n^3)$ 사이입니다. 이를 다시 말하면 feature가 2배로 늘어나면 계산 시간이 5.3배에서 8배 정도 증가됨을 의미합니다.

 sklearn의 LinearRegression 클래스가 사용하는 SVD 방법은 $O(n^2)$입니다. 다시 말해서 feature가 두배가 되면 계산 시간은 4배가 됩니다.

 앞에서는 훈련 방식에 따른 계산 속도의 차이를 알아봤습니다. 이제 이와 별개로 Linear Regression의 예측 속도는 데이터와 feature 수에 선형적입니다. 이는 매우 빠른 축에 속한다고 할 수 있습니다.

---

4.2 Gradient Descent Algorithm

 우리는 현재 Linear Regression에서 θ를 최적화하는 첫번째 방식인 Normal Equation에 대해 알아보았습니다. 이제 가장 보편적으로 쓰이는 Gradient Descent Algorithm에 대해 알아보겠습니다.

Gradient Descent Algorithm의 아이디어는 loss function을 최소화하기 위해 반복적으로 θ를 조정하는 것입니다.

알고리즘을 단계화시켜 개념을 정립해보겠습니다.

1) 랜덤으로 θ를 정의하여 초기 위치를 정합니다.

2) θ에 대한 loss function의 기울기(gradient)를 계산하여 이 기울기가 감소하는 방향으로 θ를 조정합니다.

3) 기울기가 0이 되거나, 매우 가까울 때까지 2)를 반복합니다. 만약 0에 도달했다면 cost를 최소로 하는 최적화된 θ를 찾은 것입니다.

---

 Gradient Descent Algorithm에서 중요한 하이퍼 파라미터는 learning rate입니다. learning rate는 한번 이동할 때 얼마나 이동하냐를 결정짓습니다.

만약 learning rate가 너무 작다면 알고리즘이 수렴하기 위해 매우 많이 반복해야 하므로 시간이 오래 걸립니다.

그렇다고 너무 크다면 오히려 알고리즘이 더 큰 값으로 발산하게 만들 수도 있습니다.

---

 여기서 우리가 주목할 점은 모든 loss function이 위에서 본 그래프들처럼 무난한 포물선의 형태를 그리지 않는다는 것입니다. 만약 loss function이 훨씬 더 꼬불꼬불하여 복잡하다면 부분적으로 극소를 가지는 local minimum 그리고 전체에서 가장 최솟값을 가지는 지점인 global minimum 두가지가 존재하게 됩니다.

 근데 만약 알고리즘이 왼쪽부터 시작하여 오른쪽으로 이동하는 과정에서 global minimum에 도달하기 전에 local minimum을 만나면 이곳에서 수렴하게 된다는 문제점이 있습니다.

 다행스러운건 우리가 지금까지 다뤘던 Linear Regression에서 사용되는 loss function은 무난한 이차함수의 포물선 형태를 가지는 볼록 함수(convex function)입니다. 다시 말해서 local minimum 없이 global minimum만 존재한다는 뜻입니다.

---

 이제 feature의 범위(scale)에 따라 loss function의 형태는 어떤지, 어떤 결과를 불러오는지 알아보겠습니다.

보통 loss function은 위에서 내려다 봤을 때 원 모양이지만, feature 간 scale이 매우 다르면 길쭉한 모양이 됩니다.

 만약 일반적인 원 모양이라면 알고리즘이 global minimum으로 곧장 진행하고 있어 빠르게 수렴하게 됩니다. 그런데 오른쪽 그래프와 같이 길쭉할 경우 최솟값에는 수렴하겠으나, 시간이 오래 걸릴 것입니다.

---

4.2.1 Batch Gradient Descent

 이제 Batch GD를 시작으로 몇가지 GD 알고리즘의 종류에 대해 알아보겠습니다. 그전에 우리가 앞에서 기울기를 줄이는 방향으로 이동한다는 말을 도식화하여 Gradient Descent Algorithm식을 확인하겠습니다.

우리가 현재 θ 지점에서의 loss function의 기울기라는 것은 loss function의 편도 함수로 구합니다.

만약 feature가 여러개라면 벡터를 이용해 한번에 계산할 수도 있겠습니다.

---

 Batch GD는 θ를 조정하기 위한 기울기를 구할 때 Training data 전체를 사용합니다. 즉 우리는 매 스텝마다 기울기를 구해서 θ를 조정해나갈텐데, 만약 Training set이 매우 크다면 시간이 매우 오래 걸린다는 결론에 쉽게 도달할 수 있습니다.

 우리가 알고리즘에서 minimum에 도달하기 위해서는 기울기가 구해지면 이것이 작아지는 방향으로 가야합니다. 여기서 각 스텝의 이동 속도를 결정하는 learning rate가 사용됩니다.

이를 간단히 구현해보고 최적화된 θ까지 바로 확인해보겠습니다.

eta = 0.1
n_iterations = 1000
m = 100

theta = np.random.randn(2,1)

for iteration in range(n_iterations):
  gradients = 2/m * X_b.T.dot(X_b.dot(theta) - y)
  theta = theta - eta * gradients

theta

##결과값##
>>
array([[4.27746389],
       [2.91097464]])

---

 Normal Equation을 통해 도출했던 θ와 정확히 같습니다. 그렇다면 learnin rate가 달라짐에 따라 경사 하강 알고리즘의 스텝은 어떻게 나아가는지 그래프를 관찰해보겠습니다.

- learning rate가 작은 경우(왼쪽) : 최솟값엔 도달하겠지만 시간이 매우 오래걸릴 것입니다.

- learning rate가 적당한 경우(가운데) : 몇번의 반복만에 최솟값에 수렴하였습니다.

- learning rate가 클 경우(오른쪽) : 알고리즘이 이리저리 튀면서 스텝마다 최솟값에서 점점 멀어져 발산해갑니다.

 

그렇다면 적당한 learning rate는 어떻게 찾을 수 있을까요?

첫번재는 GridSearch를 사용합니다. 예를 들어, 0.000001부터 시작하여 10배로 불리며 각각의 결과를 관찰할 수 있겠습니다.

두번째는 반복 횟수를 크게 지정하고, 기울기 값이 어떤 ε라는 매우 작은 값보다 작으면 거의 minimum에 도달한 것으로 간주하고 중지하는 것입니다.

---

4.2.2 SGD

 Batch GD의 가장 큰 문제는 매 스텝마다 전체 Training set을 계산해야 하기에 데이터셋이 클 경우 매우 느리다는 점이었습니다.

 SGD는 매 스텝마다 전체 데이터셋이 아닌 하나의 데이터를 랜덤으로 선택하여 그 데이터에 대한 기울기를 계산합니다. 즉 매 스텝마다 사용하는 데이터가 작기 때문에 알고리즘이 매우 빠릅니다.

 그러나 빠른만큼 훨씬 불안정합니다. 다시 말해 최솟값으로 도달할 때 위아래로 튀면서 감소합니다. 그 결과 완벽하게 global minimum에 도달하진 못합니다.

하지만 불안정하다는 의미는 local minimum을 건너뛸 수 있다는 의미이기도 합니다. 따라서 SGD가 Batch GD보다 global minimum에 도달할 확률이 더 큽니다.

---

 그렇다면 SGD의 이러한 불안정성을 어떻게 해결할 수 있을까요?

 첫번째 방식은 learning rate를 점진적으로 감소기키는 것입니다. 다시 말해서 처음에는 큰 보폭으로 최솟값으로 다가가다가, 어느 순간부터 좁은 보폭으로 global minimum에 세심하게 다가가는 기법입니다.

 SGD를 구현해보겠습니다. 각 반복을 epoch라 하는데 아래 코드의 경우 1000번을 반복하도록 하는데 50번만 반복하고도 최솟값에 근접하게 도달한 케이스입니다.

n_epochs = 50
t0, t1 = 5, 50  # 학습 스케줄 하이퍼파라미터

def learning_schedule(t):
    return t0 / (t + t1)

theta = np.random.randn(2,1)  # 랜덤 초기화

for epoch in range(n_epochs):
    for i in range(m):
        if epoch == 0 and i < 20:                    
            y_predict = X_new_b.dot(theta)           
            style = "b-" if i > 0 else "r--"         
            plt.plot(X_new, y_predict, style)        
        random_index = np.random.randint(m)
        xi = X_b[random_index:random_index+1]
        yi = y[random_index:random_index+1]
        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)
        eta = learning_schedule(epoch * m + i)
        theta = theta - eta * gradients
        theta_path_sgd.append(theta)              

plt.plot(X, y, "b.")                             
plt.xlabel("$x_1$", fontsize=18)                     
plt.ylabel("$y$", rotation=0, fontsize=18)           
plt.axis([0, 2, 0, 15])                              
save_fig("sgd_plot")                                 
plt.show()

theta

##결과값##
>>
array([[4.21076011],
       [2.74856079]])

 첫 20번의 훈련스텝을 보겠습니다. 매우 불안정하다는 것 그리고 learning rate가 점점 줄어든다는 것을 알 수 있습니다.

n_epochs = 50
t0, t1 = 5, 50  # 학습 스케줄 하이퍼파라미터

def learning_schedule(t):
    return t0 / (t + t1)

theta = np.random.randn(2,1)  # 랜덤 초기화

for epoch in range(n_epochs):
    for i in range(m):
        if epoch == 0 and i < 20:                    
            y_predict = X_new_b.dot(theta)          
            style = "b-" if i > 0 else "r--"         
            plt.plot(X_new, y_predict, style)        
        random_index = np.random.randint(m)
        xi = X_b[random_index:random_index+1]
        yi = y[random_index:random_index+1]
        gradients = 2 * xi.T.dot(xi.dot(theta) - yi)
        eta = learning_schedule(epoch * m + i)
        theta = theta - eta * gradients
        theta_path_sgd.append(theta)                

plt.plot(X, y, "b.")                                 
plt.xlabel("$x_1$", fontsize=18)                     
plt.ylabel("$y$", rotation=0, fontsize=18)           
plt.axis([0, 2, 0, 15])                              
plt.show()        

---

 SGD의 경우 매 스텝마다 데이터를 랜덤으로 선택하기 때문에 어떤 데이터는 여러 번 선택되는 반면, 어떤 데이터는 한번도 반영되지 않을 수도 있습니다. 이를 방지하려면 epoch마다 데이터를 섞어주면 됩니다.

sklearn에서 SGD를 구현해보았습니다. 반복 횟수는 1000번, 0.001보다 cost가 작을때까지 반복합니다.(tol=1e-3), learning rate는 0.1로 설정하였습니다.

from sklearn.linear_model import SGDRegressor

sgd_reg = SGDRegressor(max_iter = 1000, tol = 1e-3, penalty = None, eta0 = 0.1)
sgd_reg.fit(X, y.ravel())

sgd_reg.intercept_, sgd_reg.coef_
>>
(array([4.24365286]), array([2.8250878]))

---

4.2.3 Mini-Batch Gradient Algorithm

 이제 마지막으로 알아볼 알고리즘은 Mini-Batch Gradient Algorithm입니다. 이것은 Batch GD처럼 전체도 아니고, SGD처럼 하나도 아닌 mini-batch라 부르는 작은 데이터 세트에 대한 기울기를 계산합니다.

 이 mini-batch를 크게하면, Batch GD에 가까워지므로 SGD에서 관찰됐던 불안정성이 조금씩 해결됩니다. 그러나 local minimum에 도달한 확률은 더욱 커집니다.

다음은 세 가지 Gradient Descent Algorithm의 훈련 과정을 그래프로 표현한 것입니다.

- Batch GD는 실제로 minimum에서 멈췄습니다. 그러나 SGD와 Mini-Batch GD는 근처에서 맴돕니다.

- Batch GD가 학습에 가장 많은 시간이 소요됩니다.

---

 여기까지 Linear Regression의 원리와 구현, 그리고 Gradient Descent Algorithm의 종류와 원리 그리고 구현까지 알아보았습니다.

긴 글 읽어주셔서 감사합니다. 오늘도 행복한 하루 보내시길 바랍니다:)