NumPy

Mảng và shape

Mảng NumPy (ndarray) giống mảng của C++ hơn là list của Python: mọi phần tử cùng một kiểu số, gọi là dtype, và mảng có kích thước cố định theo từng chiều, gọi là shape. Vì mọi phần tử cùng kiểu, NumPy chạy vòng lặp bên trong bằng code C, nhanh hơn nhiều so với vòng lặp Python.

Ảnh là ví dụ gặp nhiều nhất. Ảnh màu là mảng shape (cao, rộng, 3): ảnh $(480, 640, 3)$ cao 480 và rộng 640 điểm ảnh. Mỗi điểm ảnh có ba kênh đỏ, lục, lam; mỗi kênh lưu bằng 8 bit kiểu uint8, tức $2^8 = 256$ mức từ 0 tới 255. Số càng lớn càng sáng: $(0, 0, 0)$ là đen, $(255, 255, 255)$ là trắng.

shape trả lời “có bao nhiêu phần tử trên mỗi trục”, ndim là số trục, size là tích các độ dài trục. np.zeros((3, 4)) tạo mảng shape $(3, 4)$ toàn số 0, size là 12. np.arange(0, 10, 2) giống range: được [0 2 4 6 8], không lấy 10. Với bảng hai hàng và ba cột, phần tử a[1, 2] nằm ở hàng thứ hai, cột thứ ba:

import numpy as np

a = np.array([[1, 2, 3], [4, 5, 6]], dtype=np.int64)
print(a.shape, a.ndim, a.size, a.dtype)
print(np.zeros((2, 3)).dtype, np.arange(6).reshape(2, 3).tolist())
print(np.linspace(0, 1, 5))
print(np.arange(12).reshape(3, -1).shape)    # -1: NumPy tự tính trục còn lại

Kết quả(2, 3) 2 6 int64
float64 [[0, 1, 2], [3, 4, 5]]
[0.   0.25 0.5  0.75 1.  ]
(3, 4)

reshape giữ nguyên dãy giá trị và thứ tự, chỉ đổi cách chia hàng: 12 số có thể đọc thành $3 \times 4$ hay $2 \times 6$, miễn tích các trục bằng 12. Chỉ được viết -1 ở một trục; NumPy tự tính độ dài trục đó.

Kiểu dữ liệu là chỗ NumPy khác Python nhiều nhất. int Python có độ chính xác tùy ý, nhưng số nguyên NumPy có độ rộng cố định. Phép toán trên mảng nguyên có thể quay vòng khi vượt miền biểu diễn, không nhất thiết báo lỗi:

big = np.array([2**62], dtype=np.int64)
print(big * 4, big.dtype)

pixels = np.array([200, 100], dtype=np.uint8)   # ảnh thường lưu bằng uint8, giá trị 0..255
print(pixels + 100)

x = np.array([1, 2, 3])
x[0] = 2.7                                       # gán số thực vào mảng nguyên: bị cắt
print(x)

Kết quả[0] int64
[ 44 200]
[2 2 3]

$2^{62} \cdot 4 = 2^{64}$ vượt int64 và quay vòng về 0 trong phép toán này. Điểm ảnh $200 + 100 = 300$ vượt 255 và thành $300 - 256 = 44$. Đổi sang float32 trước khi xử lý rồi kẹp về miền hợp lệ bằng np.clip sẽ tránh phép cộng uint8 bị quay vòng. Dòng cuối cho thấy gán số thực vào mảng nguyên làm mất phần lẻ.

Với số thực, float64 dùng 8 byte mỗi số, float32 dùng 4 byte. float64 chính xác hơn, nhưng mạng neural thường dùng float32 vì tốn nửa bộ nhớ và chạy nhanh hơn. Nhãn lớp và chỉ số thì vẫn để kiểu nguyên.

Tính trên cả mảng

Phép toán số học, so sánh và các hàm như np.exp, np.sqrt, np.log áp dụng lên từng phần tử của mảng mà không cần vòng lặp. Cách viết này gọi là vector hóa (vectorization):

z = np.array([-2.0, 0.0, 1.0, 3.0])
print(1 / (1 + np.exp(-z)))                      # sigmoid của bài 2.3, cho cả mảng
print(z > 0, (z > 0).sum(), (z > 0).mean())
print(np.where(z > 0, z, 0.0))                   # ReLU: giữ số dương, số âm thành 0
print(np.clip(z, -1, 1), z.max(), z.argmax())

Kết quả[0.11920292 0.5        0.73105858 0.95257413]
[False False  True  True] 2 0.5
[0. 0. 1. 3.]
[-1.  0.  1.  1.] 3.0 3

So sánh cho ra mảng bool; cộng mảng bool thì True tính là 1, nên (z > 0).sum() đếm số phần tử dương và (z > 0).mean() là tỉ lệ của chúng. Đây là cách tính accuracy: (pred == y).mean().

Dòng đầu tính hàm sigmoid $1/(1+e^{-z})$ cho cả bốn số cùng lúc; bài 0.4 sẽ học hàm này. Ở đây chỉ cần thấy một công thức được áp vào từng phần tử mà không cần vòng lặp. Thử so tốc độ với vòng lặp Python trên 200 000 số. Như ở bài 0.1, kiểm hai cách ra cùng kết quả (bằng np.allclose) rồi mới đo:

import math
from timeit import repeat

rng = np.random.default_rng(0)
z = rng.normal(size=200_000)
zs = z.tolist()

def sigmoid_loop(values):
    return [1 / (1 + math.exp(-value)) for value in values]

def sigmoid_numpy(values):
    return 1 / (1 + np.exp(-values))

p_loop = sigmoid_loop(zs)
p_vec = sigmoid_numpy(z)
assert np.allclose(p_loop, p_vec)
t_loop = min(repeat(lambda: sigmoid_loop(zs), repeat=5, number=3)) / 3
t_vec = min(repeat(lambda: sigmoid_numpy(z), repeat=5, number=3)) / 3
print(np.allclose(p_loop, p_vec), p_vec.shape)
print("NumPy nhanh hơn vòng lặp:", t_vec < t_loop)

Kết quảTrue (200000,)
NumPy nhanh hơn vòng lặp: True

Hai cách cho cùng kết quả, chỉ lệch ở chữ số cuối, và NumPy nhanh hơn. Số giây trên máy đang chạy nằm trong t_loop và t_vec; hình dưới vẽ từ một lần đo như vậy.

Chỉ số và lát cắt

Mảng nhiều chiều dùng một cặp ngoặc với các chỉ số ngăn bởi dấu phẩy: a[i, j], không phải a[i][j] (cách sau vẫn chạy nhưng chậm hơn và dễ sai với lát cắt). Mỗi trục nhận một chỉ số, một lát cắt start:stop:step, một list chỉ số, hoặc một mặt nạ bool:

a = np.arange(12).reshape(3, 4)
print(a)
print(a[1, 2], a[:, 0], a[-1, ::2])
print(a[[0, 2]].shape, a[a % 5 == 0])

Kết quả[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
6 [0 4 8] [ 8 10]
(2, 4) [ 0  5 10]

a[0] (hay a[0, :]) là hàng đầu, a[:, 0] là cột đầu, $(3,)$. a[[0, 2]] lấy hàng 0 và hàng 2 (fancy indexing). a[a % 5 == 0] lấy mọi phần tử chia hết cho 5 và trả về mảng một chiều.

Có một khác biệt quan trọng với list. Lát cắt của list Python tạo bản sao; lát cắt của mảng NumPy là một view: một cách nhìn khác vào cùng vùng nhớ. Sửa view thì mảng gốc đổi theo:

a = np.arange(6)
s = a[1:4]                  # view
s[0] = 100
print(a)

f = a[[1, 2, 3]]            # fancy indexing: luôn là bản sao
f[0] = -1
print(a, np.shares_memory(a, s), np.shares_memory(a, f))

Kết quả[  0 100   2   3   4   5]
[  0 100   2   3   4   5] True False

Quy tắc: lát cắt bằng start:stop:step cho view; chỉ số bằng list hay mặt nạ bool cho bản sao. View giúp NumPy không tốn bộ nhớ khi cắt mảng lớn, nhưng nếu định sửa lát cắt mà không muốn đụng tới gốc thì phải gọi .copy(). reshape cũng thường trả về view; riêng khi mảng đã chuyển vị hay cắt cách quãng, NumPy phải chép ra bản sao. Muốn biết hai mảng có dùng chung dữ liệu không thì gọi np.shares_memory(a, b).

Lưu ý: “mặt nạ cho bản sao” chỉ đúng khi lấy ra một mảng mới, như b = a[a > 5]; sửa b sau đó không đụng tới a. Còn gán thẳng a[a > 5] = 0 thì sửa đúng các ô đó của a.

axis

Giờ đến bảng điểm: 5 học sinh, 3 môn. Mỗi hàng là một học sinh, mỗi cột là một môn. Ý nghĩa của từng trục do ta quy ước; trong học máy, quy ước gần như luôn là trục 0 là các mẫu.

scores = np.array([[9, 8, 7],
                   [5, 6, 4],
                   [4, 4, 5],
                   [9, 8, 10],
                   [7, 8, 10]], dtype=float)       # (5 học sinh, 3 môn)
print(scores.mean(axis=0))                         # trung bình từng môn
print(scores.mean(axis=1).round(2))                # trung bình từng học sinh
print(scores.argmax(axis=1), scores.max())

Kết quả[6.8 6.8 7.2]
[8.   5.   4.33 9.   8.33]
[0 1 2 2 2] 10.0

Cách nhớ: trục nào được ghi trong axis thì trục đó biến mất. scores có shape $(5, 3)$; mean(axis=0) gộp theo trục 0 (gộp 5 học sinh lại) nên còn $(3,)$, một số cho mỗi môn. mean(axis=1) gộp 3 môn nên còn $(5,)$, một số cho mỗi học sinh. Không ghi axis thì gộp tất cả thành một số.

Quy tắc giữ nguyên với mảng nhiều chiều hơn. Một batch 100 ảnh xám $28 \times 28$ có shape $(100, 28, 28)$:

imgs = rng.random((100, 28, 28))
print(imgs.mean(axis=0).shape)                     # ảnh trung bình
print(imgs.mean(axis=(1, 2)).shape)                # độ sáng trung bình của từng ảnh
print(imgs.mean(axis=(1, 2), keepdims=True).shape)

Kết quả(28, 28)
(100,)
(100, 1, 1)

keepdims=True giữ lại trục bị gộp với độ dài 1 thay vì xóa nó. Mục sau sẽ cho thấy vì sao điều đó hữu ích.

Hình dưới làm hai phép gộp đó trên ảnh thật: 1797 ảnh chữ số viết tay $8 \times 8$ có sẵn trong scikit-learn. Gộp theo trục 0 được một ảnh trung bình; gộp theo trục 1 và 2 được độ sáng của từng ảnh.

Ghép hai mảng cũng phải chọn trục. concatenate nối hai mảng dọc theo một trục đã có, ví dụ ghép 2 học sinh đầu với 3 học sinh sau thành bảng 5 học sinh. stack thêm một trục mới, ví dụ xếp bảng điểm của hai buổi kiểm tra chồng lên nhau:

joined = np.concatenate([scores[:2], scores[2:]], axis=0)
two_sessions = np.stack([scores, scores], axis=0)
print(joined.shape, np.array_equal(joined, scores))
print(two_sessions.shape)    # (buổi, học sinh, môn)

Kết quả(5, 3) True
(2, 5, 3)

Shape (2, 5, 3) có thêm một trục mới ở đầu: trục buổi kiểm tra. Ghép sai trục thì code vẫn chạy, nhưng hàng và cột không còn mang ý nghĩa như ta nghĩ.

Broadcasting

Chuẩn hóa từng môn: lấy scores $(5, 3)$ trừ trung bình từng môn $(3,)$. Hai mảng khác shape, nhưng phép trừ vẫn chạy:

mu = scores.mean(axis=0)
sigma = scores.std(axis=0)
sigma_safe = np.where(sigma == 0, 1.0, sigma)
Z = (scores - mu) / sigma_safe
print(Z.shape, Z.mean(axis=0).round(6), Z.std(axis=0))

Kết quả(5, 3) [ 0.  0. -0.] [1. 1. 1.]

sigma là độ lệch chuẩn của từng môn: điểm của môn càng trải xa trung bình thì sigma càng lớn (bài 0.7 sẽ học kỹ). Mỗi ô của Z cho biết điểm đó cao hơn hay thấp hơn trung bình môn bao nhiêu lần độ lệch chuẩn, nên so được giữa các môn. Output cho thấy mỗi cột của Z có trung bình 0 và độ lệch chuẩn 1 (-0. là một số âm rất nhỏ làm tròn thành 0). sigma_safe thay độ lệch chuẩn bằng 1 ở môn nào mọi điểm giống nhau, để không chia cho 0.

Phép trừ scores - mu chạy được dù hai mảng khác shape: NumPy lấy mu $(3,)$ trừ vào từng hàng của bảng $(5, 3)$, như thể mu được chép thành 5 hàng. Cơ chế này gọi là broadcasting. NumPy không chép thật mu ra 5 lần, nhưng mảng kết quả vẫn đủ $5 \times 3$ ô. Có hai quy tắc:

1. Xếp shape của hai mảng thẳng hàng từ bên phải. Mảng nào ít trục hơn thì coi như có thêm các trục độ dài 1 ở bên trái.
2. Ở mỗi trục, hai độ dài phải bằng nhau, hoặc một trong hai bằng 1. Trục độ dài 1 được kéo dãn cho khớp trục kia.

Dòng cuối của bảng là lỗi thường gặp khi muốn trừ trung bình của từng học sinh: scores.mean(axis=1) có shape $(5,)$, xếp thẳng từ phải với $(5, 3)$ thì không khớp. Đây là lúc dùng keepdims:

row_mean = scores.mean(axis=1)
try:
    scores - row_mean
except ValueError as e:
    print("ValueError:", e)
print((scores - scores.mean(axis=1, keepdims=True)).shape)

Kết quảValueError: operands could not be broadcast together with shapes (5,3) (5,) 
(5, 3)

Với keepdims=True, trung bình có shape $(5, 1)$ và kéo dãn đúng theo trục môn học.

Broadcasting còn thay được cả vòng lặp lồng nhau. Tính khoảng cách giữa mọi cặp điểm của hai tập $A$ ($n$ điểm) và $B$ ($m$ điểm) trong $d$ chiều, việc mà bài 2.4 (k-nearest neighbors) cần, bằng cách chèn thêm trục độ dài 1:

A = np.array([[0.0, 0.0], [3.0, 4.0]])
B = np.array([[0.0, 4.0], [3.0, 0.0]])
diff = A[:, None, :] - B[None, :, :]  # (2, 1, 2) - (1, 2, 2) -> (2, 2, 2)
D = np.sqrt((diff ** 2).sum(axis=-1))
print(diff.shape, D.shape)
print(D)
print(diff.nbytes)

Kết quả(2, 2, 2) (2, 2)
[[4. 3.]
 [3. 4.]]
64

None (hoặc np.newaxis) chèn một trục độ dài 1. Trục cuối chứa tọa độ; gộp nó để lấy một khoảng cách cho mỗi cặp, shape (n, m). Cặp đầu là từ $(0,0)$ đến $(0,4)$, khoảng cách bằng 4; cặp hàng 1, cột 0 là từ $(3,4)$ đến $(0,4)$, khoảng cách bằng 3. Bốn ô đều tính tay được, nên ví dụ nhỏ này là phép thử cho code trước khi chạy trên dữ liệu thật.

Cách này tiện nhưng tốn bộ nhớ: mảng diff có $n \cdot m \cdot d$ số, mỗi số 8 byte. Ở đây chỉ 64 byte, nhưng với $n = m = 10\,000$ điểm, $d = 64$ chiều thì diff cần 51.2 GB, chưa kể diff ** 2. Vì vậy trước khi chạy trên dữ liệu lớn, hãy tính thử kích thước các mảng trung gian; bài 2.4 sẽ có cách tính khoảng cách mà không cần mảng diff.

$(n,)$ và $(n, 1)$

Hai shape này trông giống nhau và in ra gần giống nhau, nhưng là hai thứ khác hẳn. $(n,)$ là vector một trục. $(n, 1)$ là ma trận $n$ hàng, một cột. Trộn chúng với nhau là lỗi âm thầm nguy hiểm nhất của NumPy, vì broadcasting biến nó thành một phép tính hợp lệ nhưng sai:

y = np.array([1.0, 2.0, 3.0])                      # (3,): nhãn thật
y_hat = np.array([[1.0], [2.0], [3.0]])            # (3, 1): dự đoán, đúng hoàn toàn
print((y_hat - y).shape)
print("MSE sai:", np.mean((y_hat - y) ** 2), "  MSE đúng:", np.mean((y_hat.ravel() - y) ** 2))

Kết quả(3, 3)
MSE sai: 1.3333333333333333   MSE đúng: 0.0

Dự đoán đúng hoàn toàn mà MSE ra $1.33$ thay vì $0$. Phép trừ $(3, 1) - (3,)$ bị broadcast thành ma trận $3 \times 3$ chứa mọi hiệu $\hat y_i - y_j$, rồi np.mean lấy trung bình cả 9 số. NumPy không báo lỗi gì. Với $10^5$ mẫu, ma trận này có $10^{10}$ số, lúc đó notebook mới báo hết bộ nhớ.

Lỗi này hay gặp khi model trả dự đoán dạng $(n, 1)$, ví dụ mạng neural có một đầu ra. Cách sửa là đưa về $(n,)$ bằng .ravel() hoặc .reshape(-1). Cách phòng là kiểm shape ngay trong hàm tính loss:

def mse(y_hat, y):
    assert y_hat.shape == y.shape, (y_hat.shape, y.shape)
    return np.mean((y_hat - y) ** 2)

try:
    mse(y_hat, y)
except AssertionError as e:
    print("AssertionError:", e)

Kết quảAssertionError: ((3, 1), (3,))

Đại số tuyến tính và số ngẫu nhiên

Phép nhân ma trận dùng @. Phép * là nhân từng phần tử, không phải nhân ma trận:

M = np.array([[2.0, 1.0], [1.0, 3.0]])
v = np.array([1.0, 2.0])
print(M @ v, M * v)
print(M.T @ M)
print(np.linalg.solve(M, v), np.linalg.norm(v))

Kết quả[4. 7.] [[2. 2.]
 [1. 6.]]
[[ 5.  5.]
 [ 5. 10.]]
[0.2 0.6] 2.23606797749979

• M @ v: mỗi hàng của M nhân từng cặp với v rồi cộng lại. Hàng đầu: $2\cdot1+1\cdot2=4$; hàng sau: $1\cdot1+3\cdot2=7$. Kết quả có shape $(2,)$.
• M * v: nhân từng ô, theo broadcasting: cột đầu nhân 1, cột sau nhân 2.
• M.T: chuyển vị, đổi hàng thành cột; ma trận $(2, 3)$ thành $(3, 2)$.

Nhân hai ma trận theo quy tắc shape $(n, k) @ (k, m) \to (n, m)$: số cột bên trái phải bằng số hàng bên phải. Với hai vector cùng độ dài, a @ b là tích vô hướng: $(1, 2, 3)\cdot(4, 5, 6) = 4 + 10 + 18 = 32$. Bài 0.5 sẽ giải thích phép nhân này kỹ hơn.

np.linalg.solve(M, v) giải hệ phương trình $M\mathbf{x}=\mathbf{v}$, tức $2x_0 + x_1 = 1$ và $x_0 + 3x_1 = 2$. Kết quả $(0.2, 0.6)$; kiểm tay: $2\cdot0.2+0.6=1$ và $0.2+3\cdot0.6=2$. np.linalg.norm(v) là độ dài của vector: $\sqrt{1^2+2^2}\approx 2.236$.

Số ngẫu nhiên dùng một bộ sinh có seed, để chạy lại cho đúng kết quả cũ:

rng = np.random.default_rng(42)
print(rng.integers(0, 10, size=5), rng.normal(size=2).round(3))
print(rng.permutation(5), rng.choice(["A", "B", "C"], size=4))
rng2 = np.random.default_rng(42)
print(rng2.integers(0, 10, size=5))                # cùng seed: cùng dãy số

Kết quả[0 7 6 4 4] [ 0.941 -1.951]
[3 4 1 0 2] ['C' 'C' 'C' 'B']
[0 7 6 4 4]

Seed cố định cả một dãy số. Hai bộ sinh cùng seed, gọi theo cùng thứ tự, cho cùng các số. Nhưng mỗi lần gọi trên một rng lại lấy phần tiếp theo của dãy, nên hai lần gọi liên tiếp ra hai kết quả khác nhau. Trong code cũ hay gặp np.random.seed(...); cách này dùng một bộ sinh chung cho cả chương trình, dễ bị đoạn code khác gọi xen vào làm lệch dãy, nên nên dùng default_rng.

rng = np.random.default_rng(42)
a = rng.integers(0, 10, size=5)
b = rng.integers(0, 10, size=5)

Sắp xếp, đếm, cộng dồn và lưu mảng

NumPy có sẵn các thuật toán quen thuộc từ lập trình thi đấu: sắp xếp, đếm, cộng dồn, tìm nhị phân. Học máy dùng chúng để lấy $k$ lớp có điểm cao nhất, đếm số mẫu của mỗi nhãn, chia điểm thành các khoảng.

scores = np.array([7.5, 9.0, 4.0, 9.0, 6.5])
order = np.argsort(-scores, kind="stable")       # chỉ số theo điểm giảm dần
print(order, scores[order])
top2 = np.argpartition(-scores, 2)[:2]           # hai chỉ số điểm cao nhất, không sắp cả mảng
print(np.sort(top2))
labels = np.array([2, 0, 2, 1, 2, 0])
print(np.unique(labels, return_counts=True))
print(np.bincount(labels, minlength=4))
print(np.cumsum([3, 1, 4, 1, 5]))
edges = np.array([5.0, 6.5, 8.0])                # ngưỡng điểm: dưới 5, 5-6.5, 6.5-8, từ 8
print(np.searchsorted(edges, [4.0, 6.5, 9.0], side="right"))

Kết quả[1 3 0 4 2] [9.  9.  7.5 6.5 4. ]
[1 3]
(array([0, 1, 2]), array([2, 1, 3]))
[2 1 3 0]
[ 3  4  8  9 14]
[0 2 3]

argsort trả về chỉ số. Sắp -scores để được thứ tự giảm dần; kind="stable" giữ hai điểm 9.0 theo thứ tự xuất hiện. argpartition giống nth_element của C++: đưa $k$ chỉ số lớn nhất lên đầu trong thời gian tuyến tính nhưng không sắp chúng.

• unique: các giá trị khác nhau và số lần xuất hiện của mỗi giá trị.
• bincount: ô $i$ là số lần gặp số $i$; minlength=4 để có cả ô của nhãn 3 dù không xuất hiện.
• cumsum: mảng cộng dồn (prefix sum).
• searchsorted: tìm nhị phân trên mảng đã sắp; side="right" giống upper_bound của C++. Điểm 6.5 rơi vào khoảng 2, tức từ 6.5 đến dưới 8.

Đặc trưng hay embedding tính lâu thì lưu lại bằng định dạng của NumPy; lần sau nạp ngay, giữ nguyên shape và dtype:

from pathlib import Path
from tempfile import TemporaryDirectory

with TemporaryDirectory() as directory:
    path = Path(directory) / "features.npz"
    np.savez_compressed(path, X=np.arange(6, dtype=np.float32).reshape(2, 3), y=np.array([0, 1]))
    with np.load(path) as data:
        print(data.files, data["X"].shape, data["X"].dtype, data["y"])

Kết quả['X', 'y'] (2, 3) float32 [0 1]

np.save ghi một mảng ra tệp .npy; np.savez_compressed gom nhiều mảng có tên vào một tệp .npz. Lưu bằng CSV thì mất dtype và chậm với mảng lớn. Trên Kaggle, tệp ghi vào /kaggle/working được giữ lại sau khi lưu phiên bản notebook.