[Incarnate the Algorithm] Linear Regression

Linear Regression problem은 주어진 (N + 1) - dimension space의 input에 대해 hyperplane을 fitting하는 것이다.

hyperplane이란 linear function으로 h(w) = a + b^Tw로 표현된다. (특별한 표기가 없으면 1D 행렬은 열벡터의 형태이므로 w와 행렬곱을 하기 위해 행벡터꼴인 b의 전치행렬(b^T)로 바꿔준 것이다.

구하고자 하는 것은 hyperplane이 정답을 의미하는 output y_p에 대해 최적의 parameter w vector를 구하는 것이다.

 

Linear Regression 문제에서 주로 쓰이는 cost function은 least square cost function이다.

least square cost function

자주 쓰이는 이유는 이 함수가 convex function이기 때문이다.

이는 매우 중요한데, local optimizatino method의 최대 단점은 cost function이 복잡할수록 global minimum을 구하지 못할 가능성이 높다는 것이다.

하지만 linear regression에서 least square를 cost function으로 적용하면 항상 global minimum을 구할 수 있다.

* 위 그림에는 표시가 되어있지 않지만, y의 개수인 p*2만큼을 시그마결과값에서 나눠주기도 한다. 실제로 밑의 코드에서는 abs_least_square함수에서 np.size(y), i.e., y의 개수로 나눠준다. 이는 사용자의 필요에 따라 유동적으로 변화가 가능한 부분이다.

 

 

이제 위의 정보를 바탕으로 코드를 구현해보자.

100개의 원소를 가진 주어진 data x와 y가 존재한다.

데이터는 다음의 분포를 가진다.

x, y 데이터 분포

 

Least square function을 cost function로 한 코드를 구현하였다.

# cost function 구현
def model(x, w):
    a = w[0] + np.dot(x.T, w[1:]) # w[0] : bias
    return a.T

def least_square(w, x, y):
    cost = np.sum((model(x, w) - y)**2)
    return cost/float(y.size)

def abs_least_square(w, x, y): # outlier가 존재할 때 유용한 function
    cost = np.sum(abs(model(x, w) - y))
    return cost/float(y.size)

 

먼저 least square function에 대해 계산하여 그래프를 그린 후에 outlier를 고려한 regression function을 구해줄 것이다.

이전의 포스트에서 사용한 gradient descent algorithm을 활용해주자.

https://yeppyteen2.tistory.com/19

[Incarnate the Algorithm] Gradient Descent

 

이때, g 함수의 인수가 달라졌으므로 이에 유의하여 function을 수정해주어야 한다.

from autograd import grad 
import autograd.numpy as np 

# gradient descent function 
def gradient_descent(g,alpha,max_its,w,x,y):
    # compute gradient module using autograd
    gradient = grad(g) # get the differentiated function of g

    # run the gradient descent loop
    weight_history = [w] # weight history container
    cost_history = [g(w,x,y)] # cost function history container
    for k in range(max_its):
        # evaluate the gradient
        grad_eval = gradient(w, x, y)

        # take gradient descent step
        w = w - alpha*grad_eval
        
        # record weight and cost
        weight_history.append(w)
        cost_history.append(g(w, x, y))
    return weight_history,cost_history

 

이제 gradient_descent 함수를 이용하여 계산할 수 있다.

x_data = x_data.reshape((1,100)) # model 함수 내에서 행렬 곱 연산을 해야 하므로 행렬의 모양을 (100,)에서 (1,100)으로 바꿔준다.
y_data = y_data.reshape((1,100))
g = least_square
w0 = 0.1*np.random.randn(2,1)
max_its = 1000
alpha_choice = 10**(-2)
w_history, cost_history = gradient_descent(g,alpha_choice,max_its,w0,x_data,y_data)

이때 배열의 형태를 곱 연산에 적절하게 변환시키는 과정이 필요하다.

 

 

구한 값들로 그려준 cost function과 w value의 history plot은 다믕과 같다.

초반 부분에 비해 일정 횟수 이상부터는 감쇄가 거의 일어나지 않는 것을 알 수 있다.

 

효과적으로 코드가 잘 작동하는 것을 알 수 있다.

마지막으로 구한 w값을 이용하여 regression function을 그려보자.

# regression function 그리기

x_data = x_data.reshape((1,100)) # 행렬 곱 연산을 위해 array의 형태를 바꿔준다
W = []
W.append(w1[-1]) # 최종적으로 구한 w value는 마지막 값이다
W.append(w2[-1])
W = np.array(W).reshape((2,1)) # 행렬 곱 연산을 위해 array의 형태를 바꿔준다
A = model(x_data, W)

x_data2 = x_data.reshape((100,)) # 그림을 그리기 위해 다시 일반적인 array의 형태로 바꿔준다
A = A.reshape((100,)) # 그림을 그리기 위해 다시 일반적인 array의 형태로 바꿔준다
idx = np.argsort(x_data2)  # get the indices that would sort x
x_sorted = x_data2[idx]  # sort x in ascending order
A_sorted = A[idx]

plt.figure()

plt.scatter(x_data, y_data, alpha=0.3) # x, y data는 각각 100개
plt.plot(x_sorted, A_sorted, c='r')

잘 그려진 것을 알 수 있다.

 

 

상당히 괜찮은 결과를 얻은 것을 알 수 있다.

하지만 지금까지 사용한 least square function은 데이터셋에 outlier가 있을 때 치명적인 단점을 보인다. 이 데이터는 outlier라고 할 만한 데이터가 존재하지 않지만 실제 데이터의 경우에는 outlier가 존재하는 경우가 훨씬 많다. 그럴 때 대체할 수 있는 것이 제곱 대신에 절댓값을 취하는 것이다. 위에서 구현한 abs_least_square function이 그것이다.

약간의 코드를 수정하여 같은 과정을 반복하면 다음의 history plot과 regression function을 얻을 수 있다.

 

least square function에서 제곱을 취하는 것이 절댓값을 취하는 것보다 더 높은 정확도를 보여야 할 것이다. 구한 두 function을 비교해보면 이를 쉽게 확인할 수 있다.

outlier가 없는 경우에는 least square가 더 효과적이다