Gần như mọi mô hình học máy hiện đại, từ hồi quy tuyến tính đến mạng nơ-ron hàng tỷ tham số, đều được huấn luyện bằng cùng một ý tưởng: đo sai số, tính độ dốc, bước một bước nhỏ ngược chiều độ dốc, lặp lại. Ý tưởng đó là gradient descent. Bài này tự cài đặt gradient descent bằng NumPy cho bài toán hồi quy tuyến tính (dự đoán giá thuê phòng), tự kiểm tra gradient bằng sai phân hữu hạn, chạy thật và đối chiếu với nghiệm đóng cùng LinearRegression của scikit-learn. Bạn sẽ thấy tận mắt vì sao learning rate quyết định hội tụ hay phân kỳ, và vì sao chuẩn hoá dữ liệu không phải là "mẹo" mà là điều kiện để thuật toán chạy ổn định.

Bài dựa trên NumPy cho machine learning: vector hoá và broadcasting (toàn bộ vòng lặp bên dưới được vector hoá) và đặt vào quy trình đã mô tả ở Machine learning là gì: bài toán, dữ liệu và quy trình một dự án AI. Môi trường chạy lại: Apple M4, 16 GB, Python 3.10, NumPy 1.23.1, scikit-learn 1.1.2 (chỉ để đối chiếu cuối bài), cùng bộ phiên bản với hai bài trước của series; mã chỉ dùng phép toán NumPy cơ bản (default_rng có từ NumPy 1.17) nên chạy được trên các bản cũ hơn.

Bài toán: hồi quy tuyến tính là gì

Hồi quy tuyến tính dự đoán một số thực từ các đặc trưng bằng một tổ hợp tuyến tính: y_hat = b + w1*x1 + w2*x2 + .... Với ma trận đặc trưng X kích thước (n mẫu, p đặc trưng), ta thêm một cột toàn số 1 để gộp b vào vector tham số theta, khi đó dự đoán là X @ theta. Mục tiêu là chọn theta để tổng bình phương sai số là nhỏ nhất. Tài liệu scikit-learn về Ordinary Least Squares mô tả đúng như vậy: mô hình cực tiểu hoá tổng bình phương phần dư giữa giá trị quan sát và giá trị dự đoán.

Hàm mất mát (loss) ta dùng là mean squared error:

L(theta) = (1/n) * sum( (X @ theta - y)^2 )
dL/dtheta = (2/n) * X.T @ (X @ theta - y)

Công thức thứ hai là gradient: vector chỉ hướng làm loss tăng nhanh nhất. Gradient descent đi ngược hướng đó: theta = theta - lr * gradient, trong đó lr là learning rate (bước học). Vì loss của hồi quy tuyến tính là hàm lồi (một cái "bát" duy nhất) nên gradient descent với lr đủ nhỏ luôn đi tới cùng một đáy với nghiệm đóng. Đó là lý do đây là bài tập lý tưởng để kiểm chứng cài đặt của chính mình.

Dữ liệu tổng hợp: giá thuê phòng

Ta sinh dữ liệu có đáp án biết trước để so sánh: giá thuê (triệu đồng/tháng) = 1,2 + 0,11 × diện tích (m²) − 0,25 × khoảng cách tới trung tâm (km) + nhiễu. Dữ liệu hoàn toàn giả lập, không phải dữ liệu thật của ai.

import numpy as np

rng = np.random.default_rng(42)
n = 200
dien_tich = rng.uniform(15, 60, n)
khoang_cach = rng.uniform(0.5, 12, n)
nhieu = rng.normal(0, 0.4, n)
y = 1.2 + 0.11 * dien_tich - 0.25 * khoang_cach + nhieu
X = np.column_stack([dien_tich, khoang_cach])
Xb = np.column_stack([np.ones(n), X])      # thêm cột bias

Nghiệm đóng làm mốc so sánh

Trước khi lặp, hãy lấy đáp án chuẩn bằng bình phương tối thiểu. Không nên tự nghịch đảo X.T @ X; dùng np.linalg.lstsq (giải bằng phân rã số ổn định hơn). scikit-learn cũng tính nghiệm bằng singular value decomposition của X, với chi phí cỡ O(n_samples × n_features²) theo tài liệu.

theta_ls, *_ = np.linalg.lstsq(Xb, y, rcond=None)
print("[lstsq] theta =", np.round(theta_ls, 4))
[lstsq] theta = [ 1.0692  0.1136 -0.2538]

Nghiệm không trùng 1,2 / 0,11 / −0,25 vì có nhiễu và chỉ 200 mẫu; đây là nghiệm tối ưu cho tập dữ liệu này và là đích mà gradient descent phải tới.

Cài đặt gradient descent từng bước

Hai hàm đủ cho toàn bộ thuật toán gradient descent: loss và gradient, cả hai là một dòng NumPy vector hoá. Vòng lặp chính cập nhật theta, ghi lại lịch sử loss, dừng khi loss gần như không đổi hoặc khi phát hiện phân kỳ.

def mse(theta, Xb, y):
    r = Xb @ theta - y
    return float(r @ r / len(y))

def grad(theta, Xb, y):
    r = Xb @ theta - y
    return 2.0 / len(y) * (Xb.T @ r)

def gradient_descent(Xb, y, lr, n_iter, tol=1e-10):
    theta = np.zeros(Xb.shape[1])
    hist = [mse(theta, Xb, y)]
    for i in range(n_iter):
        theta = theta - lr * grad(theta, Xb, y)
        hist.append(mse(theta, Xb, y))
        if not np.isfinite(hist[-1]) or hist[-1] > 1e12:
            return theta, hist, i + 1, "phân kỳ"
        if abs(hist[-2] - hist[-1]) < tol:
            return theta, hist, i + 1, "hội tụ"
    return theta, hist, n_iter, "hết vòng lặp"

Kiểm tra gradient bằng sai phân hữu hạn

Lỗi phổ biến nhất khi tự cài đặt là gradient sai dấu hoặc sai hệ số mà loss vẫn giảm "tạm được". Cách kiểm tra chuẩn là so gradient giải tích với gradient số: (L(theta+eps) - L(theta-eps)) / (2*eps) cho từng thành phần.

[gradient check] giải tích: [ -11.7752  -507.90364  -63.09334]
[gradient check] số học   : [ -11.7752  -507.90364  -63.09334]
[gradient check] sai lệch tối đa: 1.80e-08

Sai lệch cỡ 1e-9 cho thấy công thức đúng. Hãy tập thói quen chạy bước này mỗi khi tự viết gradient, đặc biệt khi chuyển sang mô hình phức tạp hơn.

Learning rate: nhỏ quá thì chậm, lớn quá thì phân kỳ

Chạy gradient descent trên dữ liệu thô (diện tích cỡ 15–60, khoảng cách cỡ 0,5–12) với ba learning rate:

--- Dữ liệu chưa chuẩn hoá ---
lr=0.01    phân kỳ      sau      4 vòng | MSE=1.227e+13 | theta=[ -2101.439 -87051.327 -13146.515]
lr=0.001   phân kỳ      sau     17 vòng | MSE=3.96e+12 | theta=[ 1193.885 49456.305  7468.715]
lr=0.0005  hết vòng lặp sau  50000 vòng | MSE=0.1562 | theta=[ 1.043  0.114 -0.253]

Chỉ có lr=0.0005 sống sót, và sau 50.000 vòng bias vẫn là 1,043 thay vì 1,069. Kết quả này có giải thích toán học rõ ràng. Với loss bình phương, ma trận Hessian là H = (2/n) * X.T @ X. Gradient descent ổn định khi lr < 2 / lambda_max, với lambda_max là trị riêng lớn nhất của H; còn tốc độ hội tụ theo hướng nhạy nhất bị chặn bởi trị riêng nhỏ nhất. Tính trực tiếp trên dữ liệu của ta (Xb là ma trận có thêm cột 1, Zb là ma trận sau z-score, cùng dữ liệu ở trên):

for ten, M in [("chưa chuẩn hoá", Xb), ("z-score", Zb)]:
    H = 2.0 / n * M.T @ M              # Hessian của MSE
    ev = np.linalg.eigvalsh(H)         # trị riêng, sắp tăng dần
    print(f"{ten:<15} lambda_min={ev[0]:.4g} lambda_max={ev[-1]:.4g} "
          f"số điều kiện={ev[-1]/ev[0]:.4g} | lr tối đa ổn định = 2/lambda_max = {2/ev[-1]:.4g}")

Kết quả (đầy đủ trong stability.py):

chưa chuẩn hoá  lambda_min=0.148 lambda_max=3164 số điều kiện=2.138e+04 | lr tối đa ổn định = 2/lambda_max = 0.0006321
z-score         lambda_min=1.948 lambda_max=2.052 số điều kiện=1.053 | lr tối đa ổn định = 2/lambda_max = 0.9747

Ngưỡng 0,00063 giải thích vì sao 0,001 phân kỳ còn 0,0005 thì không. Số điều kiện (tỉ lệ trị riêng lớn trên nhỏ) là 21.380: mặt loss giống một thung lũng hẹp và dài, nên bước đủ nhỏ để không văng khỏi thành thung lũng lại quá nhỏ để đi dọc đáy. Đây là lý do thực tế để chuẩn hoá trước khi chạy gradient descent.

Chuẩn hoá z-score giúp gradient descent hội tụ nhanh

Chuẩn hoá z-score trừ trung bình và chia độ lệch chuẩn từng cột, đưa mọi đặc trưng về cùng thang đo. Quan trọng: trung bình và độ lệch chuẩn phải tính trên tập huấn luyện rồi áp lại cho dữ liệu mới. Tài liệu Stochastic Gradient Descent của scikit-learn viết rõ rằng SGD nhạy với thang đo đặc trưng nên rất nên chuẩn hoá, và cùng một phép chuẩn hoá phải áp cho dữ liệu kiểm tra.

mu, sd = X.mean(axis=0), X.std(axis=0)
Z = (X - mu) / sd
Zb = np.column_stack([np.ones(n), Z])
for lr in [1.2, 1.0, 0.3, 0.1, 0.01]:
    th, h, it, st = gradient_descent(Zb, y, lr, 50000)
    print(f"lr={lr:<5} {st:<12} sau {it:>6} vòng | MSE={h[-1]:.6f}")
--- Sau khi chuẩn hoá z-score ---
lr=1.2   phân kỳ      sau     35 vòng | MSE=1181440853141.061523
lr=1.0   phân kỳ      sau    264 vòng | MSE=1046219400655.340210
lr=0.3   hội tụ       sau     15 vòng | MSE=0.156107
lr=0.1   hội tụ       sau     57 vòng | MSE=0.156107
lr=0.01  hội tụ       sau    559 vòng | MSE=0.156107

So với 50.000 vòng chưa tới nơi, dữ liệu chuẩn hoá chỉ cần 57 vòng ở lr=0.1. Ngưỡng ổn định 0,9747 cũng khớp: lr=1.0 vượt ngưỡng và phân kỳ (chậm hơn lr=1.2 vì chỉ vượt chút ít), còn 0,3 hội tụ chỉ sau 15 vòng. Đường loss theo vòng lặp với lr=0.1:

MSE theo vòng lặp (lr=0.1): [16.7844, 10.7777, 1.9248, 0.3444, 0.1582, 0.1561]   (vòng 0, 1, 5, 10, 20, 50)

Tham số học được nằm trong không gian đã chuẩn hoá, nên cần đổi về thang gốc: w = w_z / sd và b = b_z - sum(w_z * mu / sd). Sau đó so với nghiệm đóng:

Tham số GD (đã quy về thang gốc): bias=1.0692, w=[ 0.1136 -0.2538]
Tham số lstsq                    : bias=1.0692, w=[ 0.1136 -0.2538]
Chênh lệch tối đa: 2.38e-06

Gradient descent tự cài đặt cho kết quả khớp nghiệm đóng tới cỡ 1e-6. Và LinearRegression của scikit-learn cho cùng nghiệm:

[sklearn LinearRegression] bias=1.0692, w=[ 0.1136 -0.2538]
MSE lstsq     : 0.156107
MSE sklearn   : 0.156107

Batch, mini-batch và stochastic gradient descent

Bản trên là batch gradient descent: mỗi lần cập nhật dùng toàn bộ n mẫu. Với dữ liệu hàng triệu dòng, mỗi bước quá đắt. Hai biến thể của gradient descent thường gặp, được tổng hợp trong bài khảo sát An overview of gradient descent optimization algorithms của Sebastian Ruder (arXiv 2016):

  • Stochastic gradient descent (SGD): cập nhật sau mỗi một mẫu; bước rẻ, nhiều bước, nhưng loss dao động.
  • Mini-batch gradient descent: cập nhật sau mỗi nhóm nhỏ (thường 32–256 mẫu); cân bằng giữa độ ổn định và tốc độ, tận dụng được phép nhân ma trận vector hoá, và là lựa chọn mặc định khi huấn luyện mạng nơ-ron.
def minibatch_gd(Zb, y, lr, epochs, batch, seed=0):
    r = np.random.default_rng(seed)
    theta = np.zeros(Zb.shape[1]); m = len(y)
    log = []
    for ep in range(epochs):
        idx = r.permutation(m)                 # xáo trộn mỗi epoch
        for s in range(0, m, batch):
            bi = idx[s:s + batch]
            theta -= lr * grad(theta, Zb[bi], y[bi])
        log.append(mse(theta, Zb, y))
    return theta, log
--- SGD / mini-batch (z-score, 30 epoch) ---
batch (full)   lr=0.1   MSE epoch 1=10.7777 | epoch 30=0.1561 | số lần cập nhật/epoch=1
mini-batch 32  lr=0.1   MSE epoch 1=0.8538 | epoch 30=0.1567 | số lần cập nhật/epoch=7
SGD (1 mẫu)    lr=0.02  MSE epoch 1=0.1651 | epoch 30=0.1641 | số lần cập nhật/epoch=200

Quan sát đáng giá: sau 1 epoch, batch đầy đủ mới đi được một bước (MSE 10,78), mini-batch đã xuống 0,85 và SGD đã gần đáy (0,165). Nhưng sau 30 epoch, batch đầy đủ chạm đáy chính xác (0,1561), mini-batch lệch nhẹ (0,1567) và SGD vẫn dao động quanh đáy (0,1641) vì mỗi bước chỉ nhìn một mẫu nhiễu. Muốn gradient descent kiểu SGD hội tụ sát đáy, người ta giảm learning rate dần; scikit-learn có sẵn các lịch giảm như invscaling và optimal cho SGDRegressor. Lưu ý lr=0.02 cho SGD là do tôi chọn nhỏ hơn để tránh dao động quá lớn; đây là một siêu tham số cần tinh chỉnh chứ không phải hằng số của thuật toán.

Gradient descent trong các thư viện lớn

  • scikit-learn: LinearRegression dùng nghiệm đóng (SVD), không lặp; SGDRegressor và SGDClassifier dùng SGD, hợp với dữ liệu rất lớn, và có partial_fit để học dần từng đợt (online/out-of-core). Tài liệu khuyên đặt max_iter khoảng 10^6 chia số mẫu làm giá trị khởi đầu và xáo trộn dữ liệu trước khi huấn luyện.
  • Deep learning: PyTorch và TensorFlow tính gradient tự động (automatic differentiation) thay vì ta tự đạo hàm, rồi dùng các biến thể như momentum hay Adam trên mini-batch. Phần lõi của gradient descent vẫn là vòng lặp theta -= lr * grad bạn vừa viết.
  • Hồi quy logistic, SVM tuyến tính: cùng khung, chỉ đổi hàm loss và gradient.

Khi nào không nên tự cài gradient descent

  • Hồi quy tuyến tính với vài chục nghìn mẫu và vài chục đặc trưng: dùng LinearRegression (nghiệm đóng) là nhanh, chính xác, không cần chọn learning rate.
  • Đặc trưng tương quan mạnh: nghiệm đóng có thể rất nhạy với nhiễu (tài liệu scikit-learn cảnh báo ma trận thiết kế gần suy biến làm phương sai lớn); nên thêm regularization như Ridge thay vì chỉ đổi thuật toán tối ưu.
  • Sản phẩm thật: tự cài gradient descent chỉ để hiểu và để dạy. Dùng thư viện đã được kiểm thử, có xử lý ổn định số học và các trường hợp biên.
  • Loss không lồi (mạng nơ-ron): gradient descent chỉ đảm bảo tới điểm dừng cục bộ; đây là chủ đề của các bài deep learning sau trong series.

Lỗi thường gặp

  1. Không chuẩn hoá đặc trưng: số điều kiện lớn làm gradient descent chậm, buộc dùng learning rate rất nhỏ như ví dụ ở trên.
  2. Gradient sai mà vẫn thấy loss giảm: luôn chạy gradient check với sai phân hữu hạn.
  3. Loss tăng hoặc thành inf/nan: learning rate vượt 2/lambda_max; giảm lr, thường chia 3 đến 10 lần.
  4. Chuẩn hoá bằng thống kê của toàn bộ dữ liệu (cả tập kiểm tra): rò rỉ dữ liệu; phần Pipeline ở bài sau xử lý tận gốc.
  5. Quên chuyển tham số về thang gốc khi báo cáo hệ số cho người dùng nghiệp vụ: hệ số "0,52 trên diện tích" ở không gian z-score không có nghĩa là 0,52 triệu mỗi m².
  6. Quên xáo trộn dữ liệu giữa các epoch khi dùng SGD hoặc mini-batch, nhất là khi dữ liệu sắp xếp theo thời gian hay theo nhãn.
  7. Dừng theo số vòng cố định không xét loss: dùng thêm tiêu chí dừng theo dung sai như tol ở trên.

Checklist

  • Viết loss và gradient dưới dạng vector hoá; kiểm tra gradient bằng sai phân hữu hạn.
  • Chuẩn hoá đặc trưng bằng thống kê của tập huấn luyện trước khi chạy gradient descent.
  • Chọn learning rate nhỏ hơn 2/lambda_max; theo dõi đường loss theo vòng lặp.
  • So kết quả với nghiệm đóng hoặc LinearRegression để chắc cài đặt đúng.
  • Dùng mini-batch khi dữ liệu lớn; giảm learning rate dần nếu cần hội tụ sát đáy.
  • Quy tham số về thang đo gốc trước khi diễn giải.

Mã nguồn đầy đủ nằm trong code/s7-10 (gd.py, stability.py, requirements.txt), chạy bằng python gd.py. Mã gradient descent này và mọi con số trong bài đều là output thật từ lần chạy trên Apple M4, 16 GB; dữ liệu là giả lập, hạt giống ngẫu nhiên 42.

Bài tiếp theo

Machine learning với scikit-learn: Pipeline, cross-validation và tuning. Sau khi hiểu cơ chế tối ưu bên dưới, ta dùng scikit-learn để ghép tiền xử lý và mô hình vào một Pipeline, đánh giá bằng cross-validation và tìm siêu tham số mà không rò rỉ dữ liệu.