선형 회귀(Linear Regression)
5 min read

선형 회귀(Linear Regression)

Base Class

  • init Function에서는 최적 $\beta$를 찾아가는 과정에서 Gradient Descent, 즉 을Loss Function을 최소화하는 방향으로 조금씩 값을 대입해보면서 최종 값을 찾아가는 형태이기 때문에 Learning Rate와 Iteration에 대한 변수를 설정하고 있다.
  • 그리고 Weight 값에 대해서는 Uniform Distribution을 초기화를 해주는데, 이 부분은 기본적으로 0으로 초기화해도 Training은 가능할 수 있으나, 이후에 Weight를 학습해나가는 과정에서 모든 Weight가 동일하게 학습될 가능성이 있고 이로 인해서 "Symmetry"를 Break하는데 실패하고 결과적으로 학습이 제대로 진행되지 않을 수 있다.
  • Fitting에서는 Bias항을 상수로 1을 제일 앞에 넣어주고, n_iteration만큼 돌면서 MSE(Mean Square Error)를 계산하고 예측과 실제값간의 차이, 즉 오차과와 입력값의 내적을 통해서  grad_w를 계산한다. 이 때 Ridge나 Lasso Regression의  경우 추가로 Panelty를 주기 위한 값을 추가한다.
  • Regularization이후 이 값만큼 가중치에서 빼주는 식으로 Weight 값을 찾는다. 이 때 빼주는 값은 Learning Rate를 통해서 조절한다.
class Regression(object):
    """ Base regression model. Models the relationship between a scalar dependent variable y and the independent 
    variables X. 
    Parameters:
    -----------
    n_iterations: float
        The number of training iterations the algorithm will tune the weights for.
    learning_rate: float
        The step length that will be used when updating the weights.
    """
    def __init__(self, n_iterations, learning_rate):
        self.n_iterations = n_iterations
        self.learning_rate = learning_rate

    def initialize_weights(self, n_features):
        """ Initialize weights randomly [-1/N, 1/N] """
        limit = 1 / math.sqrt(n_features)
        self.w = np.random.uniform(-limit, limit, (n_features, ))

    def fit(self, X, y):
        # Insert constant ones for bias weights
        X = np.insert(X, 0, 1, axis=1)
        self.training_errors = []
        self.initialize_weights(n_features=X.shape[1])

        # Do gradient descent for n_iterations
        for i in range(self.n_iterations):
            y_pred = X.dot(self.w)
            # Calculate l2 loss
            mse = np.mean(0.5 * (y - y_pred)**2 + self.regularization(self.w))
            self.training_errors.append(mse)
            # Gradient of l2 loss w.r.t w
            grad_w = -(y - y_pred).dot(X) + self.regularization.grad(self.w)
            # Update the weights
            self.w -= self.learning_rate * grad_w

    def predict(self, X):
        # Insert constant ones for bias weights
        X = np.insert(X, 0, 1, axis=1)
        y_pred = X.dot(self.w)
        return y_pred

Linear Regression

  • Fitting을 할 때는 앞서 언급한 Gradient Descent를 써도 되고,  OLS(Ordinal Least Square)를 이용해서 구할 수 있는데 Feature 수에 따라 OLS는 행렬계산으로 인해서 Performance가 떨어지는 경향이 있고 반대로 Gradient Descent는 Feature 수에 따른 Performance 하락은 적은데 반해서 Hyparparameter 값을 초반에 정해야 하는 번거로움이 있다.
  • 여튼 OLS를 통해서 구할 때는 $\hat\beta$는 $ (\bm{X}'\bm{X})^{-1}\bm{X}'\bm{y}$이기 때문에 Moore-Penrose 유사 역행렬을 이용한다.
  • 기본적으로 역행렬은 $n×n$ 정방 행렬(square matrix)에서만 정의되는데 데이터포인트의 수와 Feature의 수가 항상 동일할 수가 없기 때문에 유사 역행렬(pseudo inverse matrix)을 활용해야 한다.
  • 그리고 Moore-Penrose 유사 역행렬은 SVD 특이값 분해($A = U \Sigma V^T$)를 이용한다. 중간에 U, S, V,X_sq_reg_inv 변수를 구하는 부분이 그 부분으로 Moore-Penrose는 행렬 $A$를 $A^+ = V \Sigma ^+ U^T$로 분해해준다.
class LinearRegression(Regression):
    """Linear model.
    Parameters:
    -----------
    n_iterations: float
        The number of training iterations the algorithm will tune the weights for.
    learning_rate: float
        The step length that will be used when updating the weights.
    gradient_descent: boolean
        True or false depending if gradient descent should be used when training. If 
        false then we use batch optimization by least squares.
    """
    def __init__(self, n_iterations=100, learning_rate=0.001, gradient_descent=True):
        self.gradient_descent = gradient_descent
        # No regularization
        self.regularization = lambda x: 0
        self.regularization.grad = lambda x: 0
        super(LinearRegression, self).__init__(n_iterations=n_iterations,
                                            learning_rate=learning_rate)
    def fit(self, X, y):
        # If not gradient descent => Least squares approximation of w
        if not self.gradient_descent:
            # Insert constant ones for bias weights
            X = np.insert(X, 0, 1, axis=1)
            # Calculate weights by least squares (using Moore-Penrose pseudoinverse)
            U, S, V = np.linalg.svd(X.T.dot(X))
            S = np.diag(S)
            X_sq_reg_inv = V.dot(np.linalg.pinv(S)).dot(U.T)
            self.w = X_sq_reg_inv.dot(X.T).dot(y)
        else:
            super(LinearRegression, self).fit(X, y)

Result

  • 실행을 하기 위해서는 데이터를 준비해야 하기 때문에 sklearn.datasets의 make_regression으로 X,y를 생성한다. 모두 array로 Return 된다.
  • 그리고 해당 데이터를 train_test_split해준다. sklearn Package를 불러올 수도 있지만, 간단하게 np.random.shuffle로 index를 shuffle한 이후 해당 index를 이용해서 Data를 섞고 Split 해주었다.
  • MSE 값을 가지고 와서 plot으로 보면 mse가 감소하다가 약 20회를 넘어서면 감소하지 앟는 pleatau 상태에 이르게 되는 것을 볼 수 있다.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.datasets import make_regression


def shuffle_data(X, y, seed=None):
    """ Random shuffle of the samples in X and y """
    if seed:
        np.random.seed(seed)
    idx = np.arange(X.shape[0])
    np.random.shuffle(idx)
    return X[idx], y[idx]

def train_test_split(X, y, test_size=0.5, shuffle=True, seed=None):
    """ Split the data into train and test sets """
    if shuffle:
        X, y = shuffle_data(X, y, seed)
    # Split the training data from test data in the ratio specified in
    # test_size
    split_i = len(y) - int(len(y) // (1 / test_size))
    X_train, X_test = X[:split_i], X[split_i:]
    y_train, y_test = y[:split_i], y[split_i:]

    return X_train, X_test, y_train, y_test

def mean_squared_error(y_true, y_pred):
    """ Returns the mean squared error between y_true and y_pred """
    mse = np.mean(np.power(y_true - y_pred, 2))
    return mse

X, y = make_regression(n_samples=100, n_features=1, noise=20)

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4)

n_samples, n_features = np.shape(X)

model = LinearRegression(n_iterations=100)

model.fit(X_train, y_train)

# Training error plot
n = len(model.training_errors)
training, = plt.plot(range(n), model.training_errors, label="Training Error")
plt.legend(handles=[training])
plt.title("Error Plot")
plt.ylabel('Mean Squared Error')
plt.xlabel('Iterations')
plt.show()

y_pred = model.predict(X_test)
mse = mean_squared_error(y_test, y_pred)
print ("Mean squared error: %s" % (mse))
image

References