Một vòng for thuần Python duyệt một triệu phần tử mất cỡ trăm mili giây; cùng phép tính viết bằng mảng NumPy chỉ mất dưới một mili giây. Chênh lệch hai bậc độ lớn đó là lý do gần như mọi thư viện machine learning trong Python, từ pandas đến scikit-learn, được xây trên NumPy, còn PyTorch có kiểu tensor riêng nhưng thiết kế API theo tư duy mảng của NumPy và chuyển qua lại với ndarray được (torch.from_numpy dùng chung bộ nhớ). Bài này giải thích hai khái niệm quyết định khả năng viết NumPy hiệu quả: vector hoá (thay vòng lặp bằng phép toán trên cả mảng) và broadcasting (quy tắc cho phép mảng khác shape cùng tính toán). Toàn bộ code chạy thật với NumPy 1.23.1 trên Python 3.10, máy Apple M4, 16 GB RAM; output dán nguyên văn.

NumPy đứng ngay dưới pandas. Bài trước Làm sạch dữ liệu crawl bằng pandas: chuẩn hoá, khử trùng lặp, kiểu dữ liệu cho ra bảng có kiểu đúng; bài này đi xuống tầng bên dưới, nơi các cột đó thật sự là mảng số. Cách dựng môi trường sạch ở Môi trường Python cho dự án dữ liệu và AI: venv, Poetry và quản lý phiên bản vẫn áp dụng: pip install numpy==1.23.1 trong venv riêng.

Vì sao NumPy nhanh hơn vòng lặp Python

Để thấy vì sao NumPy nhanh: một danh sách Python là mảng con trỏ tới các đối tượng float rải rác trong bộ nhớ; mỗi phép cộng phải kiểm tra kiểu, tra phương thức và cấp phát đối tượng kết quả. Một mảng ndarray của NumPy là một khối bộ nhớ liền mạch chứa các số cùng kiểu (dtype) cùng thông tin shape và strides. Vòng lặp của NumPy trên khối đó chạy bằng mã C đã biên dịch, không qua trình thông dịch Python cho từng phần tử. Bài báo Array Programming with NumPy (Harris và cộng sự, Nature 2020) mô tả đây chính là ý tưởng lõi: một vài khái niệm mảng đơn giản tạo nên cách tổ chức và phân tích dữ liệu khoa học gọn và mạnh, đồng thời NumPy trở thành nền tảng của cả hệ sinh thái Python khoa học.

Ta đo thử bằng hàm tính sai số bình phương trung bình (MSE) giữa giá trị thật và dự đoán trên một triệu phần tử. Phiên bản vòng lặp và phiên bản vector hoá cho cùng một kết quả:

import time
import numpy as np

rng = np.random.default_rng(42)

# --- A1. Vòng lặp Python và vector hoá: tổng bình phương sai số ---
n = 1_000_000
y_true = rng.normal(100, 15, n)
y_pred = y_true + rng.normal(0, 5, n)

def mse_loop(a, b):
    s = 0.0
    for i in range(len(a)):
        d = a[i] - b[i]
        s += d * d
    return s / len(a)

def mse_vec(a, b):
    d = a - b
    return float(np.mean(d * d))

def best_of(fn, *args, repeat=5):
    ts = []
    for _ in range(repeat):
        t0 = time.perf_counter(); r = fn(*args); ts.append(time.perf_counter() - t0)
    return r, min(ts)

r1, t1 = best_of(mse_loop, y_true, y_pred, repeat=3)
r2, t2 = best_of(mse_vec, y_true, y_pred)
print(f"MSE vòng lặp : {r1:.6f}  ({t1*1000:8.1f} ms)")
print(f"MSE vector   : {r2:.6f}  ({t2*1000:8.1f} ms)")
print(f"Nhanh hơn    : {t1/t2:.0f} lần")
print("Kết quả khớp :", np.isclose(r1, r2))
MSE vòng lặp : 24.962240  (   117.0 ms)
MSE vector   : 24.962240  (     0.6 ms)
Nhanh hơn    : 206 lần
Kết quả khớp : True

Cả hai cho ra cùng giá trị MSE, còn phiên bản vector hoá nhanh hơn 206 lần trong lần đo này (best of 5 với bản vector, best of 3 với bản vòng lặp). Con số cụ thể đổi theo máy và kích thước mảng, thứ tự độ lớn thì ổn định. Với mảng chỉ vài chục phần tử, chi phí gọi hàm NumPy có thể lấn át lợi ích.

Shape, dtype và bộ nhớ của mảng NumPy

Trước khi viết vector hoá bằng NumPy cần nắm ba thuộc tính. shape là bộ số chiều, ví dụ (1000, 20) là 1000 mẫu, 20 đặc trưng. dtype là kiểu phần tử, mặc định float64 (8 byte). nbytes là dung lượng khối dữ liệu. Quy ước chung của machine learning: ma trận đặc trưng có shape (n_mẫu, n_đặc_trưng), đúng như tài liệu scikit-learn mô tả.

# --- A2. Shape, dtype, bộ nhớ ---
x = rng.normal(size=(1000, 20))
print("\nshape", x.shape, "dtype", x.dtype, "ndim", x.ndim, "bytes", x.nbytes)
print("float32:", x.astype(np.float32).nbytes, "bytes")
shape (1000, 20) dtype float64 ndim 2 bytes 160000
float32: 80000 bytes

Đổi sang float32 giảm một nửa bộ nhớ, đổi lại độ chính xác thấp hơn (khoảng 7 chữ số thập phân so với khoảng 16 của float64).

Broadcasting trong NumPy: quy tắc so shape từ phải sang trái

Khi hai mảng khác shape, NumPy không báo lỗi ngay mà thử broadcasting. Theo tài liệu chính thức về broadcasting, NumPy so sánh shape của hai mảng từng phần tử, bắt đầu từ chiều cuối (ngoài cùng bên phải) và đi sang trái. Hai chiều tương thích khi chúng bằng nhau hoặc một trong hai bằng 1. Hai mảng không cần cùng số chiều: chiều thiếu được xem như 1 ở phía trái. Chiều có kích thước 1 được "kéo giãn" để khớp chiều kia mà không tạo bản sao thật.

Shape AShape BKết quảLý do
(4, 3)(3,)(4, 3)B được xem là (1, 3)
(4, 1)(3,)(4, 3)cả hai chiều bị kéo giãn
(8, 1, 6, 1)(7, 1, 5)(8, 7, 6, 5)so từ phải: 1 và 5, 6 và 1, 1 và 7, thiếu và 8
(4, 3)(4,)lỗichiều cuối 3 khác 4, không bằng và không phải 1

Chạy thật bốn trường hợp trên, cộng cách sửa lỗi bằng cách thêm trục mới. Các đoạn code nhãn B1 đến B6 (từ đây đến hết phần View và copy, trừ đoạn đo hiệu năng perf.py ở giữa) là các phần của cùng một script broadcasting.py: chúng dùng chung biến X, Z, rng, nên hãy chép nối tiếp theo thứ tự, bắt đầu từ đoạn dưới đây (đã có sẵn phần import và khởi tạo), thì sẽ ra đúng output dán nguyên văn:

import numpy as np

np.set_printoptions(suppress=True, precision=3)
rng = np.random.default_rng(0)

# --- B1. Quy tắc broadcasting: so shape từ phải sang trái ---
a = np.zeros((4, 3)); b = np.array([1., 2., 3.])
print("(4,3) + (3,)  ->", (a + b).shape)
print("(4,1) + (3,)  ->", (np.arange(4)[:, None] + np.arange(3)).shape)
print("(8,1,6,1) + (7,1,5) ->", (np.zeros((8, 1, 6, 1)) + np.zeros((7, 1, 5))).shape)
try:
    np.zeros((4, 3)) + np.zeros(4)
except ValueError as e:
    print("Lỗi:", e)
print("Sửa bằng cột:", (np.zeros((4, 3)) + np.zeros(4)[:, None]).shape)
(4,3) + (3,)  -> (4, 3)
(4,1) + (3,)  -> (4, 3)
(8,1,6,1) + (7,1,5) -> (8, 7, 6, 5)
Lỗi: operands could not be broadcast together with shapes (4,3) (4,) 
Sửa bằng cột: (4, 3)

Trường hợp lỗi là điều bạn sẽ gặp thường nhất. Ở đây (4, 3) + (4,) thất bại vì NumPy căn phải: chiều cuối của hàng là 3, của vector là 4. Muốn cộng một giá trị cho mỗi hàng, ta biến vector thành cột bằng [:, None] (tương đương np.newaxis), shape thành (4, 1) và broadcasting kéo nó thành (4, 3). Kỹ thuật đáng thuộc: thêm trục bằng None để chỉ định rõ chiều nào được kéo giãn.

Ví dụ 1: chuẩn hoá đặc trưng (z-score) bằng NumPy, không vòng lặp

Các thuật toán dựa trên khoảng cách hay gradient nhạy với thang đo: thu nhập hàng chục triệu đồng sẽ át tuổi vài chục. Cách chuẩn hoá chuẩn là trừ trung bình và chia độ lệch chuẩn theo từng cột. Với X shape (500, 3), X.mean(axis=0) cho shape (3,), và phép (X - mu) / sd là broadcasting (500, 3) với (3,):

# --- B2. Chuẩn hoá cột (z-score): X là (n_mẫu, n_đặc_trưng) ---
X = np.column_stack([rng.normal(30, 8, 500),               # tuổi
                     rng.normal(15_000_000, 4_000_000, 500),  # thu nhập (VND)
                     rng.integers(0, 12, 500)])             # số đơn hàng
mu = X.mean(axis=0); sd = X.std(axis=0)                    # shape (3,)
Z = (X - mu) / sd                                          # (500,3) - (3,) -> (500,3)
print("\nmean trước:", X.mean(axis=0))
print("mean sau  :", np.abs(Z.mean(axis=0)).round(6), "| std sau:", Z.std(axis=0).round(6))
mean trước: [      29.785 14723333.38         5.64 ]
mean sau  : [0. 0. 0.] | std sau: [1. 1. 1.]

Trung bình ba cột trước khi chuẩn hoá chênh nhau từ khoảng 30 đến hơn 14 triệu; sau khi chuẩn hoá, trị tuyệt đối trung bình bằng 0 và độ lệch chuẩn bằng 1 ở cả ba cột (in dạng làm tròn). Trong dự án thật, mu và sd phải tính từ tập huấn luyện rồi mới áp cho tập kiểm thử để tránh rò rỉ dữ liệu.

Ví dụ 2: khoảng cách giữa mọi cặp điểm

Bài toán lặp đi lặp lại trong phân cụm, k-NN và tìm kiếm tương đồng: cho P gồm 5 điểm và Q gồm 3 điểm trong không gian 2 chiều, tính ma trận khoảng cách 5 × 3. Thêm trục để P[:, None, :] có shape (5, 1, 2) và Q[None, :, :] có shape (1, 3, 2), hiệu của chúng broadcast thành (5, 3, 2); bình phương, cộng theo trục cuối rồi lấy căn:

# --- B3. Khoảng cách giữa mọi cặp điểm, không dùng vòng lặp ---
P = rng.normal(size=(5, 2)); Q = rng.normal(size=(3, 2))
diff = P[:, None, :] - Q[None, :, :]                       # (5,1,2)-(1,3,2) -> (5,3,2)
D = np.sqrt((diff ** 2).sum(axis=-1))                      # (5,3)
D_loop = np.array([[np.linalg.norm(p - q) for q in Q] for p in P])
print("\ndiff shape:", diff.shape, "| D shape:", D.shape, "| khớp vòng lặp:", np.allclose(D, D_loop))
diff shape: (5, 3, 2) | D shape: (5, 3) | khớp vòng lặp: True

Kết quả khớp với phiên bản hai vòng lặp lồng nhau dùng np.linalg.norm. Cái giá của NumPy ở đây: biểu thức này tạo mảng trung gian (5, 3, 2), tức tăng bộ nhớ theo tích số điểm và số chiều. Tài liệu NumPy cũng cảnh báo mảng tạm có thể lớn. Phần hiệu năng bên dưới đo đúng chuyện này.

Ví dụ 3: softmax ổn định số học theo từng hàng

Softmax biến vector điểm số thành xác suất và xuất hiện ở cuối hầu hết bộ phân loại. Viết thẳng exp(z) / sum(exp(z)) dễ tràn số khi điểm số lớn. Cách chuẩn là trừ giá trị lớn nhất của từng hàng trước khi lấy exp; kết quả toán học không đổi vì tử và mẫu cùng nhân một hệ số. Tham số keepdims=True giữ shape (2, 1) thay vì (2,) để phép trừ và chia broadcast đúng theo hàng:

# --- B4. Softmax ổn định số học theo từng hàng ---
logits = np.array([[2.0, 1.0, 0.1], [1000.0, 1001.0, 1002.0]])
def softmax(z):
    z = z - z.max(axis=1, keepdims=True)                   # trừ max: tránh tràn exp
    e = np.exp(z)
    return e / e.sum(axis=1, keepdims=True)
print("\nsoftmax:\n", softmax(logits).round(4))
with np.errstate(over="ignore", invalid="ignore"):
    naive = np.exp(logits) / np.exp(logits).sum(axis=1, keepdims=True)
print("softmax ngây thơ (hàng 2):", naive[1])
softmax:
 [[0.659 0.242 0.099]
 [0.09  0.245 0.665]]
softmax ngây thơ (hàng 2): [nan nan nan]

Hàng thứ hai có điểm số cỡ 1000: bản ngây thơ cho nan vì exp(1000) vượt giới hạn của float64 và thành vô cực, còn bản trừ max cho xác suất hợp lệ.

Ví dụ 4: phép nhân ma trận cho hồi quy tuyến tính

Với NumPy, dự đoán của hồi quy tuyến tính cho mọi mẫu cùng lúc chỉ là một phép nhân ma trận cộng độ lệch: y_hat = X @ w + b. Toán tử @ gọi vào thư viện đại số tuyến tính tối ưu (BLAS) nên nhanh hơn nhiều so với tự lặp. Đây cũng là nền tảng để tính cả gradient bằng vài phép nhân ma trận ở các bài về huấn luyện mô hình.

# --- B5. Dự đoán của hồi quy tuyến tính bằng phép nhân ma trận ---
w = np.array([0.5, -0.2, 0.1]); bias = 3.0
y_hat = Z @ w + bias                                       # (500,3)@(3,) -> (500,)
print("\ny_hat shape:", y_hat.shape, "| 3 giá trị đầu:", y_hat[:3].round(3))
y_hat shape: (500,) | 3 giá trị đầu: [2.882 2.962 3.8  ]

Hiệu năng: vòng lặp, broadcasting và phép nhân ma trận

Đo thật trên cùng bài toán khoảng cách, với 2000 điểm và 500 điểm tham chiếu trong 64 chiều. Ba cách cho cùng kết quả: hai vòng lặp lồng nhau; broadcasting (2000, 1, 64) - (1, 500, 64); và dùng đẳng thức ||p - q||² = ||p||² + ||q||² - 2·p·q để quy về một phép nhân ma trận P @ Q.T.

import time
import numpy as np

rng = np.random.default_rng(1)
N, M, D = 2000, 500, 64
P = rng.normal(size=(N, D)); Q = rng.normal(size=(M, D))

def t(fn, repeat=3):
    best = 1e9
    for _ in range(repeat):
        t0 = time.perf_counter(); r = fn(); best = min(best, time.perf_counter() - t0)
    return r, best * 1000

def loop():
    out = np.empty((N, M))
    for i in range(N):
        for j in range(M):
            out[i, j] = np.sqrt(((P[i] - Q[j]) ** 2).sum())
    return out

def bcast():
    return np.sqrt(((P[:, None, :] - Q[None, :, :]) ** 2).sum(-1))

def matmul():
    # ||p-q||^2 = ||p||^2 + ||q||^2 - 2 p.q
    d2 = (P ** 2).sum(1)[:, None] + (Q ** 2).sum(1)[None, :] - 2 * P @ Q.T
    return np.sqrt(np.maximum(d2, 0))

a, ta = t(loop, repeat=1)
b, tb = t(bcast)
c, tc = t(matmul)
print(f"Khoảng cách {N}x{M} điểm, {D} chiều")
print(f"vòng lặp Python : {ta:9.1f} ms")
print(f"broadcasting    : {tb:9.1f} ms  (tạm {N*M*D*8/1e6:.0f} MB)")
print(f"phép nhân ma trận: {tc:8.1f} ms")
print("khớp:", np.allclose(a, b), np.allclose(a, c))
Khoảng cách 2000x500 điểm, 64 chiều
vòng lặp Python :    1390.4 ms
broadcasting    :      65.3 ms  (tạm 512 MB)
phép nhân ma trận:      8.0 ms
khớp: True True

Thứ tự ngoài dự đoán với nhiều người mới: broadcasting nhanh hơn vòng lặp khoảng 21 lần, nhưng phép nhân ma trận nhanh hơn broadcasting khoảng 8 lần nữa trên máy này, vì không phải tạo mảng trung gian (2000, 500, 64). Con số 512 MB là kích thước của một mảng trung gian như vậy (2000 × 500 × 64 × 8 byte); biểu thức (P[:, None, :] - Q[None, :, :]) ** 2 còn tạo thêm một mảng cùng cỡ nên đỉnh bộ nhớ thực tế có thể gấp đôi, cỡ 1 GB.

Các tỉ lệ 21 lần và 8 lần phụ thuộc thư viện BLAS, phiên bản NumPy và máy (chạy lại trên máy khác có thể ra tỉ lệ rất khác, chẳng hạn BLAS nhanh hơn thì phép nhân ma trận càng vượt xa); hãy chạy lại để kiểm chứng trên môi trường của bạn. Bài học: chọn công thức ít bộ nhớ trung gian là bước tối ưu thứ hai sau vector hoá. Cái giá của cách thứ ba là sai số làm tròn: hiệu hai số gần nhau có thể ra âm rất nhỏ, vì vậy ta kẹp bằng np.maximum(d2, 0) trước khi lấy căn.

View và copy: bẫy làm hỏng dữ liệu gốc

Theo tài liệu NumPy về copy và view, cắt lát cơ bản (X[:, 0], X[:2]) trả về view dùng chung bộ nhớ với mảng gốc, còn đánh chỉ mục nâng cao (danh sách chỉ số hoặc mảng boolean) trả về copy. Sửa một view là sửa luôn mảng gốc:

# --- B6. View và copy ---
v = X[:, 0]; v[0] = -1; print("\nSau khi sửa view, X[0,0] =", X[0, 0])
c = X[:, 1].copy(); c[0] = -1; print("Sau khi sửa copy, X[0,1] =", round(X[0, 1]))
print("Slice dùng chung bộ nhớ         :", np.shares_memory(X, X[:2, 0]))
print("Fancy indexing dùng chung bộ nhớ:", np.shares_memory(X, X[[0, 1], 0]))
Sau khi sửa view, X[0,0] = -1.0
Sau khi sửa copy, X[0,1] = 20171572
Slice dùng chung bộ nhớ         : True
Fancy indexing dùng chung bộ nhớ: False

Sau khi gán v[0] = -1 trên view của cột đầu, X[0, 0] đổi theo; còn sửa bản .copy() thì cột thu nhập gốc giữ nguyên giá trị. np.shares_memory cho câu trả lời dứt khoát: cắt lát dùng chung bộ nhớ, chỉ mục nâng cao thì không. Khi cần bản độc lập hãy gọi .copy() tường minh.

Khi nào không nên vector hoá bằng NumPy, và lỗi thường gặp

  • Mảng trung gian quá lớn: broadcasting của NumPy tạo mảng tạm theo tích các chiều. Với dữ liệu lớn, chia thành lô (chunk) hoặc đổi sang công thức dùng phép nhân ma trận như trên.
  • Logic rẽ nhánh phụ thuộc từng phần tử: nếu mỗi bước phụ thuộc kết quả bước trước (như bộ lọc đệ quy), vector hoá khó áp dụng; cân nhắc Numba hoặc viết lại thuật toán.
  • Nhầm shape (n,) và (n, 1): trong NumPy, vector một chiều không phải cột. (n, 1) - (n,) cho ma trận (n, n) thay vì (n,), lỗi lặng lẽ và rất khó thấy. Luôn print(x.shape) khi nghi ngờ.
  • Tràn và mất chính xác: mảng NumPy kiểu int32 tràn không báo lỗi; chọn dtype chủ động.
  • Lạm dụng vòng lặp np.vectorize: tài liệu NumPy nói rõ hàm này chỉ là tiện ích cú pháp, thực chất vẫn lặp trong Python và không đem lại tốc độ của vector hoá thật sự.

Checklist viết NumPy cho machine learning

  • Trong NumPy, quy ước shape (n_mẫu, n_đặc_trưng); in shape ở mỗi bước khi gỡ lỗi.
  • Thay vòng lặp qua phần tử bằng phép toán trên cả mảng, so kết quả với bản vòng lặp bằng np.allclose.
  • Dùng axis và keepdims=True để mean, sum, max broadcast ngược đúng chiều.
  • Thêm trục bằng [:, None] trong NumPy thay vì reshape thủ công.
  • Ước lượng bộ nhớ của mảng trung gian trước khi broadcast trên dữ liệu lớn.
  • Chọn float32 hoặc float64 có chủ đích; ghi vào code lý do.
  • Đo thời gian bằng time.perf_counter hoặc timeit, lấy giá trị nhỏ nhất của vài lần chạy, ghi rõ cấu hình máy.

Tóm lại, NumPy mạnh không nhờ cú pháp ngắn mà nhờ đẩy vòng lặp xuống mã C trên khối bộ nhớ liền mạch, còn broadcasting cho phép biểu diễn phép tính giữa các shape khác nhau mà không viết vòng lặp. Nắm quy tắc so shape từ phải sang trái, biết thêm trục và để ý mảng trung gian là đủ để đọc phần lớn code machine learning bằng NumPy.

Bài tiếp theo

Có công cụ mảng trong tay, ta quay lại câu hỏi nền tảng hơn: bài toán nào thì nên dùng học máy, dữ liệu cần ra sao và một dự án AI đi qua những bước nào. Mời bạn đọc tiếp Machine learning là gì: bài toán, dữ liệu và quy trình một dự án AI.