1. Giới thiệu

Việc xây dựng các mô hình dự đoán chính xác đóng vai trò then chốt trong nhiều lĩnh vực. Hồi quy tuyến tính không chỉ là một công cụ thống kê cổ điển mà còn là nền tảng cho nhiều thuật toán học máy hiện đại. Phương pháp này cho phép xác định mức độ ảnh hưởng của từng biến độc lập đến biến phụ thuộc, từ đó đưa ra dự đoán và phân tích xu hướng.

2. Cơ sở lý thuyết

2.1. Khái niệm

Hồi quy tuyến tính (Linear Regression) là một thuật toán học có giám sát ước lượng chính xác mức độ thay đổi của biến phụ thuộc khi biến độc lập thay đổi. Mô hình này nhằm xác định mức độ ảnh hưởng cụ thể của biến độc lập đến biến phụ thuộc thông qua một mối quan hệ tuyến tính. Mô hình toán học của hồi quy tuyến tính đơn giản được biểu diễn bằng phương trình đường thẳng:

$$\hat{y} = w x + b$$

Trong đó:

• ŷ: Giá trị dự đoán.
• x: Giá trị của đặc trưng đầu vào.
• w: Trọng số (weight) hay hệ số góc (slope), cho biết mức độ ảnh hưởng của x lên giá trị dự đoán ŷ. Thay đổi w sẽ làm thay đổi độ dốc của đường thẳng.
• b: Độ lệch (bias) cho phép đường thẳng dịch chuyển lên hoặc xuống để khớp với dữ liệu tốt hơn

Để đánh giá độ phù hợp của mô hình hồi quy tuyến tính, chúng ta sử dụng hàm mất mát nhằm đo lường mức độ sai lệch giữa dự đoán và thực tế (ground truth). Mô hình được xem là hiệu quả khi hàm mất mát đạt giá trị tối thiểu.

$$L = (\hat{y}-y)^2 = (w x + b - y)^2$$

Trong đó:

• ŷ: Giá trị dự đoán.
• y: Giá trị thực tế
• (ŷ − y): Sai số.

Việc bình phương sai số đảm bảo giá trị mất mát luôn không âm, luôn khả vi tại mỗi điểm và phạt các sai số lớn nặng hơn.

2.2. Cơ chế hoạt động

Để tìm ra các giá trị w và b tối ưu nhằm cực tiểu hóa hàm mất mát, chúng ta sử dụng thuật toán tối ưu hóa Gradient Descent. Ý tưởng cốt lõi là sử dụng đạo hàm để xác định hướng thay đổi của hàm mất mát và cập nhật tham số theo hướng giảm dần giá trị hàm này.

• Nếu đạo hàm tại điểm hiện tại là dương f'(x) > 0, điều đó cho thấy hàm đang tăng, ta cần giảm giá trị của x để tiến về cực tiểu.
• Ngược lại, nếu đạo hàm là âm f'(x) < 0, hàm đang giảm, ta cần tăng x để tiếp tục đi xuống dốc.

Cơ chế hoạt động của Gradient Descent:

1. Khởi tạo: Bắt đầu với các giá trị ngẫu nhiên cho w và b
2. Tính toán gradient: Lấy đạo hàm riêng của hàm mất mát theo từng tham số:
   • $\dfrac{\partial L}{\partial w} = 2x(\hat{y}-y)$
   • $\dfrac{\partial L}{\partial b} = 2(\hat{y}-y)$
3. Cập nhật tham số: Điều chỉnh w và b bằng cách di chuyển một bước nhỏ theo hướng ngược lại của gradient. Điều này đảm bảo hàm mất mát sẽ giảm sau mỗi lần cập nhật
   • $w_{new} = w_{old} - \eta \dfrac{\partial L}{\partial w}$
   • $b_{new} = b_{old} - \eta \dfrac{\partial L}{\partial b}$
4. Lặp lại: Lặp lại bước 2 và 3 cho đến khi hàm mất mát hội tụ (giảm không đáng kể).

Trong đó, η là tốc độ học (learning rate), có vai trò kiểm soát kích thước của bước cập nhật:

• Nếu η quá lớn, mô hình có thể vượt qua điểm tối ưu và không thể hội tụ.
• Nếu η quá nhỏ, quá trình huấn luyện sẽ diễn ra rất chậm

2.3. Phân loại Gradient Descent

2.3.1. Stochastic Gradient Descent

Cập nhật tham số sau mỗi mẫu (Stochastic Gradient Descent - SGD) là phương pháp xử lý từng mẫu dữ liệu một cách độc lập và tuần tự. Quá trình này diễn ra theo các bước sau:

1. Khởi tạo các tham số a và b.
2. Lấy mẫu đầu tiên $(x_1, y_1)$.
3. Tính toán $\frac{\partial L_1}{\partial a}$ và $\frac{\partial L_1}{\partial b}$ dựa trên mẫu đầu tiên.
4. Cập nhật a và b mới
5. Lấy mẫu thứ hai $(x_2, y_2)$.
6. Tính toán $\frac{\partial L_2}{\partial a}$ and $\frac{\partial L_2}{\partial b}$ dựa vào mẫu thứ hai (dừng những giá trị a và b đã cập nhật ở mục 4).
7. Cập nhật lại a và b.
8. Lặp lại quy trình cho đến khi tất cả các điểm dữ liệu đã được xử lý.

Hiệu quả tính toán cao nhờ cập nhật tham số thường xuyên mà không cần xử lý toàn bộ tập dữ liệu, hiệu quả với tập dữ liệu quy mô lớn. Tuy nhiên, gradient có độ nhiễu cao và biến động mạnh do chỉ dựa trên một mẫu đơn lẻ mà không đại diện cho hướng tối ưu của toàn bộ tập dữ liệu và có thể dao động liên tục quanh điểm cực tiểu mà không hội tụ hoàn toàn.

2.3.2. Batch Gradient Descent

Phương pháp này có cách tiếp cận hệ thống hơn bằng cách xem xét tất cả các mẫu dữ liệu cùng một lúc trước khi thực hiện bất kỳ cập nhật tham số nào.

1. Khởi tạo tham số a và b.
2. Lấy tất cả mẫu của bộ dữ liệu
3. Đối với mỗi mẫu, hãy tính gradient riêng của nó
4. Tính trung bình tất cả các gradient đã tính để thu được một gradient tổng hợp duy nhất:
   • $\frac{\partial L_{total}}{\partial a} = \frac{1}{N} \sum_{i=1}^{N} \frac{\partial L_i}{\partial a}$
   • $\frac{\partial L_{total}}{\partial b} = \frac{1}{N} \sum_{i=1}^{N} \frac{\partial L_i}{\partial b}$
5. Thực hiện một lần cập nhật duy nhất cho a và b bằng cách sử dụng gradient trung bình này.

Gradient được ước lượng chính xác thông qua trung bình mẫu, đảm bảo hội tụ ổn định và trực tiếp đến điểm cực tiểu toàn cục đối với hàm convex. Tuy nhiên, chi phí tính toán cao do phải xử lý toàn bộ tập dữ liệu trước mỗi lần cập nhật tham số.

2.3.3. Mini-batch Gradient Descent

Mini-batch Gradient Descent (MBGD) là phương pháp kết hợp ưu điểm của cả SGD và Batch GD bằng cách chia tập dữ liệu thành các batch nhỏ có kích thước cố định. Mỗi lần cập nhật tham số được thực hiện dựa trên gradient trung bình của một mini-batch thay vì một mẫu đơn lẻ hoặc toàn bộ tập dữ liệu. Trong đó, kích thước thường được chọn theo lũy thừa của 2 (32, 64, 128, 256) để tối ưu hóa hiệu suất tính toán trên GPU. Tuy nhiên, cần điều chỉnh thêm siêu tham số batch size, ảnh hưởng trực tiếp đến độ ổn định gradient và hiệu suất tính toán.

Thực nghiệm

import numpy as np
import matplotlib.pyplot as plt
import random

# Tải dữ liệu từ file csv, bỏ qua dòng tiêu đề
data = np.genfromtxt('advertising.csv', delimiter=',', skip_header=1)
# Lấy tổng số mẫu dữ liệu
N = data.shape[0]

# Tách dữ liệu thành các cột tương ứng
tv_data = data[:, 0]
radio_data = data[:, 1]
newspaper_data = data[:, 2]
sales_data = data[:, 3]

# Hàm dự đoán giá trị đầu ra
def predict(x1, x2, x3, w1, w2, w3, b):
    return w1*x1 + w2*x2 + w3*x3 + b

# Hàm tính giá trị mất mát (loss)
def compute_loss(y_hat, y):
    return (y_hat - y)**2

# Hàm tính gradient
def compute_gradient(x1, x2, x3, y, y_hat):
    dl_dw1 = 2 * x1 * (y_hat - y)
    dl_dw2 = 2 * x2 * (y_hat - y)
    dl_dw3 = 2 * x3 * (y_hat - y)
    dl_db = 2 * (y_hat - y)
    return dl_dw1, dl_dw2, dl_dw3, dl_db

# Hàm cập nhật các trọng số
def update_weights(w1, w2, w3, b, dw1, dw2, dw3, db, lr):
    w1 = w1 - lr * dw1
    w2 = w2 - lr * dw2
    w3 = w3 - lr * dw3
    b = b - lr * db
    return w1, w2, w3, b

# Hàm khởi tạo các tham số ban đầu
def initialize_params():
    w1 = random.gauss(mu=0.0, sigma=0.01)
    w2 = random.gauss(mu=0.0, sigma=0.01)
    w3 = random.gauss(mu=0.0, sigma=0.01)
    b = 0
    return w1, w2, w3, b

# --- Phương pháp 1: Stochastic Gradient Descent (SGD) ---
def train_sgd(epochs, lr):
    w1, w2, w3, b = initialize_params()
    losses = []

    for epoch in range(epochs):
        for i in range(N):
            x1, x2, x3, y = tv_data[i], radio_data[i], newspaper_data[i], sales_data[i]

            y_hat = predict(x1, x2, x3, w1, w2, w3, b)
            loss = compute_loss(y, y_hat)
            losses.append(loss)

            dw1, dw2, dw3, db = compute_gradient(x1, x2, x3, y, y_hat)
            # Cập nhật trọng số ngay sau mỗi mẫu.
            w1, w2, w3, b = update_weights(w1, w2, w3, b, dw1, dw2, dw3, db, lr)

    return losses, (w1, w2, w3, b)

# --- Phương pháp 2: Mini-batch Gradient Descent ---
def train_mini_batch(epochs, lr, batch_size):
    w1, w2, w3, b = initialize_params()
    losses = []

    for epoch in range(epochs):
        # Lặp qua dữ liệu theo từng bước có kích thước là batch_size.
        for i in range(0, N, batch_size):
            # Khởi tạo các biến để tích lũy gradient cho mỗi lô (batch).
            acc_dw1, acc_dw2, acc_dw3, acc_db = 0, 0, 0, 0

            # Lấy dữ liệu cho lô hiện tại.
            end = i + batch_size
            x1_batch, x2_batch, x3_batch, y_batch = tv_data[i:end], radio_data[i:end], newspaper_data[i:end], sales_data[i:end]
            current_batch_size = len(x1_batch)

            # Tính tổng gradient cho lô.
            for j in range(current_batch_size):
                x1, x2, x3, y = x1_batch[j], x2_batch[j], x3_batch[j], y_batch[j]

                y_hat = predict(x1, x2, x3, w1, w2, w3, b)
                loss = compute_loss(y, y_hat)

                dw1, dw2, dw3, db = compute_gradient(x1, x2, x3, y, y_hat)

                acc_dw1 += dw1
                acc_dw2 += dw2
                acc_dw3 += dw3
                acc_db += db

            # Cập nhật trọng số một lần cho mỗi lô, sử dụng gradient trung bình.
            w1, w2, w3, b = update_weights(w1, w2, w3, b, acc_dw1/current_batch_size, acc_dw2/current_batch_size, acc_dw3/current_batch_size, acc_db/current_batch_size, lr)
            # Ghi lại giá trị loss của mẫu cuối cùng trong lô.
            losses.append(loss)

    return losses, (w1, w2, w3, b)

# --- Phương pháp 3: Batch Gradient Descent ---
def train_batch(epochs, lr):
    w1, w2, w3, b = initialize_params()
    losses = []

    for epoch in range(epochs):
        # Khởi tạo các biến để tích lũy gradient cho mỗi epoch.
        total_dw1, total_dw2, total_dw3, total_db = 0, 0, 0, 0
        total_loss = 0

        # Lặp qua toàn bộ tập dữ liệu để tính tổng gradient và loss.
        for i in range(N):
            x1, x2, x3, y = tv_data[i], radio_data[i], newspaper_data[i], sales_data[i]

            y_hat = predict(x1, x2, x3, w1, w2, w3, b)
            loss = compute_loss(y, y_hat)
            total_loss += loss

            dw1, dw2, dw3, db = compute_gradient(x1, x2, x3, y, y_hat)

            total_dw1 += dw1
            total_dw2 += dw2
            total_dw3 += dw3
            total_db += db

        # Cập nhật trọng số MỘT LẦN cho mỗi epoch, sử dụng gradient trung bình của toàn bộ tập dữ liệu.
        w1, w2, w3, b = update_weights(w1, w2, w3, b, total_dw1/N, total_dw2/N, total_dw3/N, total_db/N, lr)
        losses.append(total_loss/N)

    return losses, (w1, w2, w3, b)

if __name__ == "__main__":
    # Thiết lập các siêu tham số (hyperparameters)
    epochs = 5
    learning_rate = 1e-5
    batch_size = 32

    # Huấn luyện mô hình với ba phương pháp khác nhau
    print("Đang huấn luyện với SGD...")
    losses_sgd, params_sgd = train_sgd(epochs, learning_rate)

    print("Đang huấn luyện với Mini-batch GD...")
    losses_mini_batch, params_mini_batch = train_mini_batch(epochs, learning_rate, batch_size)

    print("Đang huấn luyện với Batch GD...")
    losses_batch, params_batch = train_batch(epochs, learning_rate)

    # Vẽ đồ thị so sánh loss
    plt.figure(figsize=(12, 6))
    plt.plot(losses_sgd, label='Loss của SGD', alpha=0.5, linestyle='--')
    plt.plot(losses_mini_batch, label=f'Loss của Mini-batch (size={batch_size})', linewidth=2)
    # Loss của Batch có ít điểm dữ liệu hơn, nên chúng ta cần xử lý việc vẽ đồ thị cho nó
    batch_iterations = np.arange(0, len(losses_sgd), N)
    if len(batch_iterations) > len(losses_batch):
        batch_iterations = batch_iterations[:len(losses_batch)]
    plt.plot(batch_iterations, losses_batch, label='Loss của Batch (mỗi epoch)', marker='o', linestyle='-')

    plt.title('So sánh sự hội tụ của Loss giữa các phương pháp GD')
    plt.xlabel('Lần lặp (cho SGD) / Bước xử lý lô (cho Mini-batch)')
    plt.ylabel('Loss (MSE)')
    plt.legend()
    plt.ylim(0, max(losses_batch) * 1.1)
    plt.show()

2.4. Đồ thị tính toán

Từ mô hình tuyến tính: ŷ = w x + b. Đồ thị tính toán (Computational graph) trực quan hóa các phép toán của mô hình. Dữ liệu và tham số là các nút đầu vào, đi qua các nút phép toán để tạo ra đầu ra và giá trị mất mát.

Cách thức hoạt động:

• Bước lan truyền xuôi (Forward pass): Các giá trị đầu vào và tham số đi từ trái sang phải qua đồ thị. Các giá trị trung gian được tính toán tuần tự cho đến khi ra được giá trị mất mát cuối cùng.
• Bước lan truyền ngược (Backward pass): Bắt đầu từ giá trị mất mát, tính toán và lan truyền các gradient từ phải sang trái bằng cách sử dụng quy tắc chuỗi (Chain Rule). Quá trình này tạo ra đạo hàm của hàm mất mát theo w và b.
• Cập nhật tham số: Sử dụng các gradient đã tính để cập nhật w và b thông qua thuật toán Gradient Descent.

Minh họa cách thực hiện đồ thị tính toán:

• Lan truyền xuôi: Với x = 6.7, y = 9.1, w = −0.34, b = 0.04, ta thực hiện lần lượt việc tính toán t, ŷ và L.

  • t = w · x = −0.34 · 6.7 = −2.278.

  • ŷ = t + b = −2.278 + 0.04 = −2.238

  • L = (ŷ − y)² = (−2.238 − 9.1)² ≈ 128.55

• Lan truyền ngược:

  • Bắt đầu từ L: Đạo hàm của L theo chính nó là 1.
  • Ngược về phép bình phương: Ta cần tính ∂L/∂ŷ

  $$ \frac{\partial L}{\partial \hat{y}} = 2(\hat{y} - y) = 2(-2.238 - 9.1) \approx -22.676 $$

  • Ngược về phép cộng: Nút cộng sẽ truyền gradient đầu vào cho các nhánh của nó.

    • Đạo hàm của L theo b:

    $$ \frac{\partial L}{\partial b} = \frac{\partial L}{\partial \hat{y}} \cdot \frac{\partial \hat{y}}{\partial b} = (-22.676) \cdot 1 = -22.676 $$

    • Đạo hàm của L theo t:

    $$ \frac{\partial L}{\partial t} = \frac{\partial L}{\partial \hat{y}} \cdot \frac{\partial \hat{y}}{\partial t} = (-22.676) \cdot 1 = -22.676 $$

  • Ngược về x: Bây giờ tính ∂L/∂w.

    $$ \frac{\partial L}{\partial w} = \frac{\partial L}{\partial t} \cdot \frac{\partial t}{\partial w} = (-22.676) \cdot x = (-22.676) \cdot 6.7 \approx -151.92 $$

Sau bước cập nhật, ta thu được các tham số mới: w = 1.179, b = 0.266. Nếu tiếp tục quá trình này thêm một lần nữa, giá trị mất mát sẽ giảm từ 128.55 xuống còn 0.868.

3. Các hàm mất mát trong hồi quy tuyến tính

Khi tối ưu hóa hàm mất mát bằng Gradient Descent, hàm mất mát cần thỏa mãn đồng thời 2 tính chất: Tính liên tục và tính khả vi. Trong đó:

• Tính liên tục

Một hàm số $f(x)$ được cho là liên tục tại điểm $a$ khi và chỉ khi 3 điều kiện sau đây thỏa mãn:

1. $f(a)$ xác định
2. $\lim_{x \to a} f(x)$ tồn tại
3. $\lim_{x \to a} f(x) = f(a)$

Một hàm số gián đoạn tại điêm $a$ nếu nó không thỏa mãn các điều kiện liên tục trên tại $a$.

• Tính khả vi

Gọi $f$ là một hàm số, và $f'$ là đạo hàm của hàm đó, với miền xác định gồm các giá trị của $x$ sao cho giới hạn sau tồn tại:

$$ f'(x)= \lim_{h \to 0} \frac{f(x+h) - f(x)}{h} $$

Một hàm số $f(x)$ được gọi là khả vi tại điểm $a$ nếu đạo hàm $f'(a)$ tồn tại. Tổng quát hơn, hàm số được gọi là khả vi trên tập $S$ nếu nó khả vi tại mọi điểm trong tập mở $S$. Một hàm khả vi là hàm mà đạo hàm $f'(x)$ tồn tại trên toàn bộ miền xác định của nó.

3.1. Mean Squared Error (MSE)

MSE đo lường sai lệch trung bình bình phương giữa giá trị dự đoán $\hat{y}$ và giá trị thực tế $y$ trong tập dữ liệu, và phạt nặng hơn cho các sai số lớn.

Trong hồi quy tuyến tính, MSE là hàm bậc hai theo tham số w, b. Điều này đảm bảo tồn tại nghiệm tối ưu toàn cục.

$$L =(\hat y-y)^2$$

$$\frac{\delta L}{\delta \hat y} = 2(\hat y - y)$$

Tuy nhiên, MSE nhạy cảm với các giá trị ngoại lai. MSE chỉ hoạt động tốt khi nhiễu tuân theo phân phối Gaussian — tức là không có các điểm ngoại lai lớn. Vì MSE giả định rằng các phần dư $e$ được lấy từ phân phối chuẩn (normal distribution):

$$ e \sim N(0,\sigma ^2) $$

3.2. Mean Absolute Error (MAE)

Sai số tuyệt đối trung bình (MAE – Mean Absolute Error) là giá trị trung bình của các sai lệch tuyệt đối giữa dự đoán và giá trị thực tế, mà không xét đến hướng sai lệch (dự đoán cao hơn hay thấp hơn).

MAE có tính chất tuyến tính, mỗi sai số đều đóng góp như nhau vào trung bình. MAE ít nhạy cảm với ngoại lai, nhưng khó tối ưu hơn MSE do không khả vi tại 0. Trong gradient descent, ta dùng subgradient hay biến thể Huber loss để cập nhật tham số.

$$L = |\hat y-y|$$

$$\frac{\delta L}{\delta \hat y} = \frac{|\hat y - y|}{\hat y - y}$$

Trong hình học, đạo hàm của một hàm tại một điểm chính là độ dốc của tiếp tuyến với đường cong tại điểm đó. Độ dốc của một đường thẳng có liên hệ trực tiếp với góc mà nó tạo với trục hoành.

Giả sử tiếp tuyến tại điểm $x=a$ tạo với trục hoành một góc $\alpha$. Do đó, nếu đạo hàm của hàm số $y=f(x)$ tại $x=a$ là $f'(a)$ thì:

$$f'(a) = tan(\alpha)$$

Khi cố định learning rate η, quy luật cập nhật tham số như sau:

$$ w_{i+1} = w_i - \eta \frac{\delta L}{\delta w} $$

Đối với MSE: Giá trị đạo hàm giảm dần khi mô hình tiến gần đến điểm cực tiểu, vì góc nghiêng nhỏ dần ($\alpha_4 < \alpha_3 $) dẫn đến quá trình hội tụ dần ổn định.

Đối với MAE: Độ lớn của đạo hàm luôn bằng 1, nên $\alpha _1 = \alpha _2$ tại mọi điểm - mô hình liên tục thực hiện các bước cập nhật có kích thước như nhau, dẫn đến dao động quanh điểm tối ưu thay vì hội tụ hay hội tụ chậm, vì các sai số lớn lại chỉ tạo ra bước cập nhật nhỏ.

3.3. Huber Loss

Dựa trên tính liên tục và khả vi trên tập số thực $\mathbb{R}$ và hoạt động tốt với learning rate cố định của MSE và ít bị ảnh hưởng bởi giá trị ngoại lai của MAE, ta xây dựng Huber loss.

$$ k(x) = \begin{cases} \frac{1}{2}x^2, & \text{nếu } |x| \leq \delta \\[6pt] a|x|+b, & \text{ngược lại } \\[6pt] \end{cases} $$

Ý tưởng:

• Khi sai số $|x|$ nhỏ $(\leq \delta)$, ta sử dụng thành phần phạt bậc hai, giúp hàm mất mát trở nên trơn tru tương tự như MSE
• Khi sai số $|x|$ lớn $(> \delta)$, ta sử dụng thành phần phạt tuyến tính, mang lại khả năng kháng ngoại lai như MAE

Xây dựng huber loss tại điểm $k(x)$:

• Xây dựng $k(x)$ liên tục tại $x=\delta$

$$ \begin{aligned} &\lim_{x \to |\delta^{-}|}k(x) = \lim_{x \to |\delta^{+}|}k(x) \\ \rightarrow &\lim_{x \to |\delta^{-}|}\frac{1}{2}x^2 = \lim_{x \to |\delta^{+}|}a|x|+b \\ \rightarrow &\frac{1}{2}\delta^2 = a|\delta| + b \end{aligned} $$

• Xây dựng $k(x)$ khả vi tại $x=\delta$

$$ \begin{aligned} \lim_{x \to |\delta^{-}|}\frac{k(x)-k(\delta)}{x-\delta} &= \lim_{x \to |\delta^{+}|}\frac{k(x)-k(\delta)}{x-\delta} \\ \rightarrow \ \ \ \ \lim_{x \to |\delta^{-}|}\frac{\frac{1}{2}x^2-\frac{1}{2}\delta^2}{x-\delta} &= \lim_{x \to |\delta^{+}|}\frac{(a|x|+b)-(a\delta+b)}{x-\delta} \\ \rightarrow \ \ \ \ \lim_{x \to |\delta^{-}|}\frac{1}{2}(x+\delta) &= \lim_{x \to |\delta^{+}|}\frac{a(|x|-\delta)}{x-\delta} \\ \rightarrow \ \ \ \ a &= \delta \\ \end{aligned} $$

Cuối cùng:

$$ \begin{aligned} a &= \delta \\ b &= -\frac{1}{2}\delta^2 \end{aligned} $$

Với một điểm dữ liệu:

$$ L(\hat{y},y) = \begin{cases} \frac{1}{2}(\hat{y}-y)^2, & \text{nếu } |\hat{y}-y| \leq \delta \\[6pt] \delta|\hat{y}-y|-\frac{1}{2}\delta^2, & \text{ngược lại} \\[6pt] \end{cases} $$

$\delta$ là một ngưỡng điều khiển mức độ bền vững của mô hình; khi giá trị của nó nhỏ hơn, mô hình trở nên bền vững hơn trước các ngoại lai.

Với phương trình:

$$ \hat{y}=wx+b $$

Đạo hàm là:

$$ L'_w = \begin{cases} x(\hat{y}-y), & \text{nếu } |\hat{y}-y| \leq \delta \\[6pt] \delta x \cdot \frac{|\hat{y}-y|}{\hat{y}-y}, & \text{ngược lại} \end{cases} $$

$$ L'_b = \begin{cases} (\hat{y}-y), & \text{nếu } |\hat{y}-y| \leq \delta \\[6pt] \delta \cdot \frac{|\hat{y}-y|}{\hat{y}-y}, & \text{ngược lại} \end{cases} $$

Giá trị ngưỡng $\delta$ do người dùng định nghĩa là tham số quan trọng, quyết định khi nào sai số được xử lý bằng hình phạt bậc hai và khi nào chuyển sang hình phạt tuyến tính.

Tuy nhiên, việc xác định giá trị tối ưu của $\delta$ là không dễ dàng và thường cần đến:

• Thông tin sẵn có về đặc trưng của dữ liệu nhiễu
• Tinh chỉnh thực nghiệm thông qua cross-validation

Điều này làm cho mô hình hồi quy dựa trên Huber Loss có cấu hình phức tạp hơn đôi chút so với các mô hình dựa trên MSE hoặc MAE tiêu chuẩn, tuy nhiên nó thường cho kết quả chính xác và ổn định hơn khi được hiệu chỉnh phù hợp.

Thực nghiệm

import numpy as np
import matplotlib.pyplot as plt

# Create data with outliers
np.random.seed(42)
N = 80
X = np.linspace(0, 10, N)
y_true = 2.5 * X + 5
noise = np.random.laplace(0, 1.0, N)
y = y_true + noise

outliers_x_high = np.random.uniform(5, 10, 8)
outliers_y_high = 2.5 * outliers_x_high + 5 + np.random.uniform(15, 30, 8)
outliers_x_low = np.random.uniform(0, 3, 3)
outliers_y_low = 2.5 * outliers_x_low + 5 - np.random.uniform(10, 15, 3)

X = np.concatenate([X, outliers_x_high, outliers_x_low]).reshape(-1, 1)
y = np.concatenate([y, outliers_y_high, outliers_y_low]).reshape(-1, 1)

# Helper functions
def predict(X, w, b):
    return X.dot(w) + b

def mse_loss(y, y_hat):
    return np.mean((y - y_hat) ** 2)

def mae_loss(y, y_hat):
    return np.mean(np.abs(y - y_hat))

def huber_loss(y, y_hat, delta=1.0):
    err = y - y_hat
    cond = np.abs(err) <= delta
    return np.mean(np.where(cond, 0.5 * err**2, delta * (np.abs(err) - 0.5 * delta)))

# Trainer
def train(X, y, loss_type='mse', delta=1.0, lr=0.01, epochs=400):
    n = len(y)
    w, b = np.random.randn(1, 1), 0.0
    losses = []

    for _ in range(epochs):
        y_hat = predict(X, w, b)
        err = y_hat - y

        if loss_type == 'mse':
            grad = err
        elif loss_type == 'mae':
            grad = np.sign(err)
        elif loss_type == 'huber':
            cond = np.abs(err) <= delta
            grad = np.where(cond, err, delta * np.sign(err))
        else:
            raise ValueError("Unknown loss type")

        dw = (1 / n) * X.T.dot(grad)
        db = (1 / n) * np.sum(grad)

        # update weights
        w -= lr * dw
        b -= lr * db

        # compute loss
        if loss_type == 'mse':
            losses.append(mse_loss(y, y_hat))
        elif loss_type == 'mae':
            losses.append(mae_loss(y, y_hat))
        else:
            losses.append(huber_loss(y, y_hat, delta))

    return w, b, losses

# Apply

w_mse, b_mse, loss_mse = train(X, y, 'mse')
w_mae, b_mae, loss_mae = train(X, y, 'mae')
w_hub, b_hub, loss_hub = train(X, y, 'huber', delta=1.5)

4. Regularization và Normalization

4.1. Regularization

Đa cộng tuyến (multicollinearity) xảy ra khi các đặc trưng đầu vào có tương quan cao, khiến ma trận thiết kế gần như suy biến. Khi đó, nhiều nghiệm của vector trọng số có thể cho kết quả dự đoán tương tự, nhưng một số nghiệm có hệ số rất lớn, làm mô hình nhạy cảm với nhiễu và dễ overfitting.

Từ trực giác hình học, khi các đặc trưng có mức độ tương quan cao, các cột trong ma trận thiết kế sẽ gần như hướng về cùng một phía trong không gian đặc trưng. Do đó, có rất nhiều tổ hợp trọng số (nhiều nghiệm khác nhau) có thể tạo ra các dự đoán gần như tương tự nhau.

Ví dụ trong dự đoán giá nhà:

• Area và Room có quan hệ chặt chẽ.
• Trọng số có thể phân bổ thành 100xArea + 0xRoom hoặc 10xArea + 5xRoom, đều cho kết quả gần giống nhau.
• Tuy nhiên, chỉ một nhiễu nhỏ cũng khiến mô hình chọn nghiệm rất khác nhau, dẫn đến thiếu ổn định.

Để khắc phục, Regularization bổ sung một thành phần phạt vào hàm mất mát nhằm giới hạn độ lớn trọng số. Cách tiếp cận này giúp mô hình ổn định hơn, giảm tác động của đa cộng tuyến và cải thiện khả năng tổng quát hóa.

Hàm mất mát L2 Regularization (Ridge) với hai đặc trưng:

$$L(w_1,w_2,b)=(\hat{y}-y)^2+\lambda(w_1^2+w_2^2)$$

• Nếu một trọng số $w_i$ quá lớn, hàm mất mát vẫn tăng lên do thành phần phạt, ngay cả khi sai số dự đoán $(\hat{y}-y)^2$ gần bằng 0.
• Khi trọng số nhỏ, giá trị hàm mất mát giảm, giúp mô hình ưu tiên nghiệm có hệ số nhỏ và ổn định hơn.

Đạo hàm của hàm mất mát sau khi thêm thành phần regularization được xác định như sau:

$$ \begin{aligned} L'_{w_i} &= 2x(\hat{y}-y)+2\lambda w_i \\ L'_b &= 2(\hat{y}-y) \end{aligned} $$

Ngoài L2, còn có:

• L1 Regularization (Lasso): thành phần phạt tỷ lệ với giá trị tuyệt đối của trọng số, khuyến khích nghiệm thưa (sparse solution).
• Elastic Net: kết hợp cả L1 và L2, vừa duy trì tính ổn định, vừa chọn lọc đặc trưng hiệu quả.

4.2. Normalization

Sự khác biệt lớn về thang đo giữa các đặc trưng có thể gây ra sự thiếu ổn định của Gradient Descent, trọng số khó diễn giải, hay regularization bị sai lệch do đặc trưng có giá trị lớn chiếm ưu thế.

Chuẩn hóa dữ liệu đưa các đặc trưng về cùng một thang đo, giúp ổn định và tăng tốc hội tụ, phản ánh đúng tầm quan trọng tương đối của đặc trưng, và duy trì tính tương thích với regularization. Chuẩn hóa không thay đổi giá trị dự đoán, nhưng ảnh hưởng đến cách diễn giải trọng số và hành vi số học của mô hình.

Trong quá trình tiền xử lý dữ liệu, hai phương pháp chuẩn hóa được sử dụng phổ biến nhất là:

• Standardization, hay Z-score normalization, chuẩn hóa dữ liệu sao cho trung bình của nó bằng 0 và độ lệch chuẩn bằng 1.

$$x'=\frac{x-\mu}{\sigma}$$

• Min-Max normalization chuẩn hóa dữ liệu sao cho mọi giá trị đều nằm trong phạm vi [0, 1].

$$x'=\frac{x-x_{min}}{x_{max}-x_{min}}$$