[Nhập môn Machine Learning] Bài 10: Cài đặt Linear Regression

Ở bài này, chúng ta sẽ sử dụng thư viện Numpy để tiến hành cài đặt thuật toán Linear Regression. Không giống những lần trước, khi chương trình của tôi chỉ xử lý được dữ liệu một chiều. Lần này, chúng ta có thể mở rộng ra cho dữ liệu nhiều chiều bằng cách vector hóa chúng. (Nếu bạn không biết vector hóa, bạn có thể xem lại bài viết này).

Trước khi bắt đầu, ta cần phải quy ước trước kích thước của các tham số trong thuật toán Linear Regression. Để cho thuận tiện, tôi sẽ sử dụng quy ước dưới đây:

• $\boldsymbol{X}$ sẽ là ma trận chứa các mẫu dữ liệu. Ma trận này có kích thước $m\times n$ tức là mỗi dòng sẽ là một điểm dữ liệu, còn mỗi dòng sẽ là một feature.
• $\boldsymbol{y}$ sẽ là một vector chứa các nhãn có độ dài $m$.
• $\boldsymbol{\theta}$ sẽ là vector chứa các tham số của ta có độ dài $n$. Tới đây bạn có thể thắc mắc tại sao lại là $n$ mà không phải $n+1$ vì ta cần một chỗ cho $x_0$ và $\theta_0$. Lý do sẽ được giải thích ngay sau đây.
• Khác với lần trước là ta phải thêm một cột các giá trị $1$ vào trước $\boldsymbol{X}$ mà ta hay gọi là $x_0$ để nhân được với tham số $\theta_0$ thì bây giờ ta không cần làm vậy nữa. Cơ chế broadcasting sẽ giúp chúng ta trong việc này nên tham số $\theta_0$ ta có thể để riêng, tôi sẽ gọi tham số này là $b$. Việc này có lợi trên thực nghiệm vì ta không phải thay đổi ma trận $\boldsymbol{X}$.

Bây giờ ta có thể bắt đầu. Đầu tiên, ta cần import thư viện Numpy vào đầu file để có thể sử dụng.

import numpy as np

Hypothesis Function

Tiếp theo, ta sẽ tiến hành cài đặt hàm tính hypothesis function của Linear Regression. Như đã biết ở phần vector hóa, ta có công thức
$$h_\theta (x) = \boldsymbol{X\theta} = \begin{bmatrix} - & (x^{(1)})^T & -\\ - & (x^{(2)})^T & -\\ & \vdots & \\ - & (x^{(m)})^T & -\\ \end{bmatrix} \begin{bmatrix} \theta_1\\ \theta_2\\ \vdots\\ \theta_n \end{bmatrix}$$
Trong code, ta làm như sau

def hypothesis_function(X, theta, b):
    """ Hàm tính dự đoán của mô hình"""
    return np.matmul(X, theta) + b # shape (m,)

Ở đây, sau khi nhân $\boldsymbol{X}$ với $\theta$, kết quả ta có được là một ma trận kích thước $m\times 1$ là kết quả của phép toán $\theta_1 x_1 + \theta_2 x_2 + \dots + \theta_n x_n$. Tuy nhiên ta vẫn cần phải cộng thêm $\theta_0$ (là biến $b$) cho $m$ kết quả này. Quả thật, ta chỉ cần cộng thẳng b vô thôi, các ndarray sẽ tự động cộng $b$ cho tất cả các phần tử thông qua cơ chế broadcasting.

Cost Function

Đây là hàm quan trọng thứ hai trong Linear Regression, công thức mà lúc trước ta đã biết được là
$$J(\boldsymbol{\theta}) = \frac{1}{2m}(\boldsymbol{X\theta} - \boldsymbol{y})^T(\boldsymbol{X\theta} - \boldsymbol{y})$$
hay
$$J(\boldsymbol{\theta}) = \frac{1}{2m}(h_\theta(\boldsymbol{X}) - \boldsymbol{y})^T(h_\theta(\boldsymbol{X}) - \boldsymbol{y})$$
Với code, ta làm như sau

def cost_function(X, y, theta, b):
    """Hàm tính toán chi phí"""
    m = X.shape[0]
    # Lấy dự đoán của mô hình
    y_hat = hypothesis_function(X, theta, b) # shape (m,)
    # Tính độ chênh lệch giữa dự đoán và nhãn
    z = y_hat - y # shape (m,)
    
    return 1/(2*m) * np.matmul(z.T, z).squeeze()

Có một hàm các bạn có thể thấy không quen đó là np.squeeze. np.squeeze giúp giảm số chiều của ndarray đi. Có thể nói là nó sẽ bỏ hết tất cả các số $1$ trong shape đi. Ví dụ một ndarray có shape (m,1,1) hay (1,m) sẽ chỉ còn là (m) sau khi squeeze. Trong trường hợp của ta trong đoạn code trên, sau khi nhân, ta sẽ có một ndarray shape (1,1), tức là nó vẫn là một ma trận hai chiều và điều này rất thiếu tự nhiên. squeeze giúp bỏ hết những chiều thừa đi, khiến giá trị chỉ còn là một số thực.

Tính toán Gradient

Với hàm tính gradient, ta có
$$\nabla J(\boldsymbol{\theta}) = \frac{1}{m}\boldsymbol{X}^T(h_\theta(\boldsymbol{X}) - \boldsymbol{y})$$
sẽ thỏa mãn phương trình $$\frac{\partial}{\partial \theta_j} = \frac{1}{m}\sum^{m}_{i=1}(h_\theta(\boldsymbol{X}) - y)x_j \ với\ 1\leq j \leq n$$
còn với $j=0$, ta chỉ cần cộng các chênh lệch lại rồi chia cho $m$. Ta sẽ code như sau

def calc_gradient(X, y, theta, b):
    """Hàm trả về tính toán gradient của biến theta và b"""
    m = X.shape[0]
    # Lấy dự đoán của mô hình
    y_hat = hypothesis_function(X, theta, b) # shape (m,)
    # Tính độ chênh lệch giữa dự đoán và nhãn
    z = y_hat - y # shape (m,)
    
    return 1/m*np.matmul(X.T, z).squeeze(), 1/m*np.sum(z)

Hàm này trả về $2$ giá trị là gradient của $\theta$ và đạo hàm của $b$. Đến đây, ta có thể bắt đầu cài đặt hàm cuối cùng là thuật toán gradient descent.

Gradient Descent

Hàm này sẽ nhận vào 2 biến X và y và trả về các tham số tối ưu cùng một mảng chứa các giá trị của Cost Function qua từng vòng lặp. Với các giá trị ban đầu của các tham số, ta có thể khởi tạo chúng với giá trị ngẫu nhiên vì gradient descent sẽ tìm được tham số tối ưu dù tham số ban đầu như thế nào.

def gradient_descent(X, y, alpha=0.01, epsilon=3e-6, max_iter=500):
    """Thực hiện thuật toán gradient descent trên tập dữ liệu X, y"""
    
    m, n = X.shape
    # Khởi tạo các giá trị
    theta = np.random.rand(n) # Khởi tạo mảng n phần tử ngẫu nhiên
    b  = 0 
    history = []
    
    while True:
        # Tính gradient
        grad_theta, grad_b = calc_gradient(X, y, theta, b)
        # Thực hiện update theo ngược chiều gradient
        theta -= alpha* grad_theta
        b -= alpha* grad_b
        #Tính giá trị của cost function
        j_cost = cost_function(X, y, theta, b)
        
        history.append(j_cost)
        # Ta dừng khi sự khác biệt trong cost đủ nhỏ hoặc số vòng lặp đủ lớn
        if abs(history[-2] - history[-1]) < epsilon:
            break
        elif len(history)> max_iter:
            break
    
    return theta, b, history

Bạn có thể tự hỏi tại sao lần trước tôi lấy điều kiện dừng là khi $\nabla \boldsymbol{\theta}$ đủ nhỏ nhưng lần này lại có điều kiện khác. Trên thực tế, có nhiều cách dừng thuật toán này lại tùy theo yêu cầu bài toán. Điều kiện lý tưởng nhất là khi mỗi đạo hàm có giá trị tuyệt đối cực nhỏ, nhưng cho đến khi đó sẽ tốn rất nhiều lần lặp của chúng ta nhưng hiệu quả (hàm Cost Function) của chúng ta chẳng giảm được bao nhiêu. Vậy nên đôi khi ta chỉ cần đáp số đủ tốt là ổn rồi và điều này phản ánh qua mức chênh lệch của Cost Function.

Và để chắc chắn hơn, tôi cũng giới hạn số lần lặp của Gradient Descent để không phải chờ quá lâu, bạn có quyền nâng giá trị max_iter nếu muốn nhiều vòng lặp hơn.

Thử nghiệm

Để chắc chắn là thuật toán chạy đúng, ta cần phải thử lại lần nữa. Ở đây, (vì lười) nên tôi sẽ lấy lại bộ dữ liệu lúc trước. Tuy nhiên cũng biến đổi đôi chút để phù hợp với quy ước ban đầu.

X = np.array([[3], 
              [5], 
              [3.25], 
              [1.5]]) # shape (4, 1)
y = np.array([1.5, 2.25, 1.625, 1.0]) # shape (4,)

Khi đã có dữ liệu, ta bắt đầu đưa chúng vào thuật toán và nhận lại các tham số sau khi chạy Gradient Descent cũng như lịch sử giá trị hàm Cost.

theta, b, history = gradient_descent(X, y, epsilon=0.0005)

Cuối cùng, để dễ thấy rằng hàm Cost của ta giảm theo thời gian, tốt nhất là nên đưa chúng lên một biểu đồ.

import matplotlib.pyplot as plt
plt.plot(range(len(history)), history)

Đây là kết biểu đồ tôi nhận được sau cùng.

Cột y là giá trị hàm Cost, cột x là số lần lặp

Tất nhiên các bạn không cần lấy dữ liệu một chiều này mà bất cứ dữ liệu nào, 2, 3 feature tùy ý.

Tổng kết

Vậy là ta đã có thể cài đặt thuật toán Linear Regression chỉ với Numpy. Đó là một bước tiến lớn. Tuy nhiên, với Linear Regression, vẫn có giới hạn các bài toán mà ta có thể giải. Nhỡ ta muốn làm một ứng dụng nhận diện khuôn mặt hay nhận diện chó mèo thì sao? Đó là lý do ở bài tiếp theo, ta sẽ tìm hiểu về bài toán phân loại và thuật toán cơ bản nhất để giải quyết bài toán này là Logistic Regression.