Data

선형 회귀(Linear Regression)

Bongho, Lee

2022년 7월 13일 — 5 min read

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))