Bài 31: Vectorization — phép tính element-wise trên array

• Hiểu vectorization và lý do nhanh hơn vòng for Python.
• Dùng toán tử số học, so sánh, logic element-wise trên array.
• Gọi đúng ufunc (np.exp, np.log, ...) thay vì hàm math.
• Aggregate theo axis=0 / axis=1 / axis=None.
• Implement ReLU, sigmoid, softmax, cosine similarity vectorized.
• Phân biệt np.dot, np.matmul, toán tử @.

Vectorization là viết phép tính cho cả array một lần, thay vì for loop chạy từng phần tử ở mức Python. Code Python chỉ gọi 1 lệnh; loop thực sự nằm trong C bên trong NumPy.

import numpy as np

a = np.array([1, 2, 3, 4])
b = np.array([10, 20, 30, 40])

a + b          # array([11, 22, 33, 44]) — 1 lệnh, loop trong C
a * b          # array([10, 40, 90, 160])

Vì sao nhanh: (1) loop trong C bỏ qua chi phí interpret Python từng vòng, (2) dữ liệu lưu liên tiếp trong bộ nhớ giúp CPU cache hiệu quả, (3) ufunc gọi instruction SIMD và một số phép như matmul còn dùng BLAS / LAPACK đã tối ưu sẵn.

Yêu cầu: 2 array tham gia phép tính phải cùng shape, hoặc shape tương thích theo quy tắc broadcasting (xem Bài 33).

Toán tử số học áp dụng cho từng cặp phần tử cùng vị trí:

a = np.array([6, 8, 10, 12])
b = np.array([2, 4,  5,  3])

a + b          # array([ 8, 12, 15, 15])
a - b          # array([ 4,  4,  5,  9])
a * b          # array([12, 32, 50, 36])
a / b          # array([3., 2., 2., 4.])     — luôn ra float
a // b         # array([3, 2, 2, 4])         — chia lấy nguyên
a % b          # array([0, 0, 0, 0])
a ** 2         # array([ 36,  64, 100, 144]) — luỹ thừa từng phần tử

Toán tử so sánh trả boolean array cùng shape:

a == b         # array([False, False, False, False])
a >  b         # array([ True,  True,  True,  True])
a <= 10        # array([ True,  True,  True, False])

Toán tử logic bitwise dùng cho boolean array — phải dùng &, |, ~, KHÔNG dùng and, or, not (xem Bài 30):

mask1 = a > 5
mask2 = a < 11
mask1 & mask2  # array([ True,  True,  True, False])
~mask1         # array([False, False, False, False])

Với 2 array khác shape sẽ raise ValueError trừ khi shape tương thích broadcasting.

ufunc là hàm element-wise đã được biên dịch sẵn cho ndarray. Hầu hết phép toán phổ biến đều có ufunc tương ứng.

x = np.array([0.0, 0.5, 1.0, 2.0])

np.sin(x)      # array([0.   , 0.479, 0.841, 0.909])
np.cos(x)      # array([1.   , 0.878, 0.540, -0.416])
np.tan(x)      # array([0.   , 0.546, 1.557, -2.185])
np.exp(x)      # array([1.   , 1.649, 2.718, 7.389])
np.log(x + 1)  # array([0.   , 0.405, 0.693, 1.099])  — ln, tránh log(0)
np.log2(x + 1) # log cơ số 2
np.log10(x+1)  # log cơ số 10
np.sqrt(x)     # array([0.   , 0.707, 1.   , 1.414])
np.abs([-3, -1, 2])   # array([3, 1, 2])

Rounding:

y = np.array([1.2, 1.5, 1.8, -1.5, -1.7])

np.round(y)    # array([ 1.,  2.,  2., -2., -2.]) — banker's rounding (round half to even)
np.floor(y)    # array([ 1.,  1.,  1., -2., -2.])
np.ceil(y)     # array([ 2.,  2.,  2., -1., -1.])

Lưu ý quan trọng: module math của Python chỉ nhận scalar; truyền array vào math.exp(a) sẽ TypeError. Luôn dùng np.exp, np.log, np.sqrt khi làm việc với array.

import math
math.exp(np.array([1, 2]))   # TypeError
np.exp(np.array([1, 2]))     # array([2.718, 7.389]) — OK

Aggregation thu nhiều phần tử về 1 con số (hoặc 1 array nhỏ hơn). Phổ biến nhất:

a = np.array([3, 1, 4, 1, 5, 9, 2, 6])

np.sum(a)      # 31
np.mean(a)     # 3.875
np.std(a)      # 2.667... — độ lệch chuẩn (ddof=0)
np.var(a)      # 7.109... — phương sai
np.min(a)      # 1
np.max(a)      # 9
np.median(a)   # 3.5

argmin, argmax trả về index của phần tử min/max — rất hay dùng để chọn class có xác suất cao nhất trong output classifier:

np.argmin(a)   # 1  — index của giá trị nhỏ nhất
np.argmax(a)   # 5  — index của giá trị lớn nhất

probs = np.array([0.1, 0.6, 0.3])
predicted_class = np.argmax(probs)   # 1

cumsum và cumprod tính tích luỹ — output cùng shape với input:

np.cumsum([1, 2, 3, 4])     # array([ 1,  3,  6, 10])
np.cumprod([1, 2, 3, 4])    # array([ 1,  2,  6, 24])

Cũng có thể gọi dưới dạng method: a.sum(), a.mean(), a.argmax().

Với array nhiều chiều, axis chỉ định trục bị collapse (thu lại). Quy tắc nhớ: aggregate theo axis nào, axis đó biến mất.

A = np.array([[ 1,  2,  3,  4],
              [ 5,  6,  7,  8],
              [ 9, 10, 11, 12]])           # shape (3, 4)

A.sum(axis=0)    # array([15, 18, 21, 24])      shape (4,) — cộng theo hàng, kết quả per cột
A.sum(axis=1)    # array([10, 26, 42])          shape (3,) — cộng theo cột, kết quả per hàng
A.sum(axis=None) # 78                           collapse toàn bộ → scalar
A.sum()          # 78                           mặc định axis=None

Diễn giải bằng shape: (3, 4) với axis=0 → bỏ axis 0 → kết quả (4,). Với axis=1 → bỏ axis 1 → kết quả (3,).

Áp dụng cho mọi aggregation, không chỉ sum:

A.mean(axis=0)     # mean của mỗi cột:  array([5., 6., 7., 8.])
A.max(axis=1)      # max của mỗi hàng:  array([ 4,  8, 12])
A.argmax(axis=1)   # index max theo hàng: array([3, 3, 3])

Trong ML, hay gặp ma trận shape (n_samples, n_features):

• axis=0 → thống kê theo từng feature (mean/std từng cột) — dùng để standardize.
• axis=1 → thống kê theo từng sample — ví dụ tổng feature của 1 sample.

Nếu muốn giữ chiều bị collapse (size 1) để tiện broadcasting, thêm keepdims=True:

A.sum(axis=1, keepdims=True)
# array([[10],
#        [26],
#        [42]])    shape (3, 1) — vẫn 2D

Anti-pattern hay gặp: lặp Python để tính trên array. Code sai dưới đây chạy được nhưng chậm, khó đọc:

# SAI — vòng for Python trên ndarray
n = len(a)
result = np.zeros(n)
for i in range(n):
    result[i] = a[i] * b[i] + c[i]

Viết lại vectorized — 1 dòng, nhanh hơn vài chục đến vài trăm lần:

# ĐÚNG — vectorized
result = a * b + c

Vài pattern thường gặp cần biết để khỏi viết loop:

• Tính min/max/mean toàn array: dùng np.min, np.max, np.mean thay for-min.
• Đếm phần tử thoả điều kiện: (a > 0).sum() thay vì loop đếm.
• If-else từng phần tử: np.where(cond, x, y) (xem Bài 30).
• Apply hàm cho từng phần tử: nếu là hàm phổ biến (sin, exp, log) thì dùng ufunc; nếu là hàm tự định nghĩa, cân nhắc np.vectorize (nhưng chỉ là wrapper for-loop, KHÔNG nhanh — chỉ giúp code gọn).

Đo thử với 1 triệu phần tử trên máy laptop điển hình:

import numpy as np
import time

n = 1_000_000
a = np.random.rand(n)
b = np.random.rand(n)

# Vectorized
t0 = time.perf_counter()
c = a + b
t_vec = time.perf_counter() - t0

# Python loop
t0 = time.perf_counter()
c = np.zeros(n)
for i in range(n):
    c[i] = a[i] + b[i]
t_loop = time.perf_counter() - t0

print(f"vectorized: {t_vec*1000:.2f} ms")
print(f"loop:       {t_loop*1000:.2f} ms")
print(f"speedup:    {t_loop / t_vec:.1f}x")

Kết quả tham khảo (con số tuyệt đối tuỳ máy, tỉ lệ thường giữ nguyên độ lớn):

• Vectorized a + b: vài ms.
• Python for-loop: vài trăm ms đến vài giây.
• Speedup: vài chục đến vài trăm lần.

Khoảng cách càng lớn khi (1) array càng lớn, (2) phép tính bên trong càng đơn giản — vì khi đó overhead của vòng for Python chiếm tỷ trọng càng cao.

Ngoài element-wise, NumPy có phép dot product / matrix multiplication — không phải element-wise. Đây là loại phép toán xuất hiện ở mọi layer dense / attention trong neural network.

# 1D × 1D → dot product (scalar)
a = np.array([1, 2, 3])
b = np.array([4, 5, 6])
np.dot(a, b)        # 1*4 + 2*5 + 3*6 = 32

# 2D × 2D → matrix multiplication
A = np.array([[1, 2],
              [3, 4]])         # shape (2, 2)
B = np.array([[5, 6],
              [7, 8]])         # shape (2, 2)
np.matmul(A, B)
# array([[19, 22],
#        [43, 50]])
A @ B                # tương đương np.matmul(A, B) — Python 3.5+

Quy tắc shape: (m, k) @ (k, n) → (m, n). Số cột của ma trận trái phải bằng số hàng của ma trận phải; nếu lệch sẽ ValueError.

Phân biệt np.dot và np.matmul:

• np.dot: 1D × 1D ra scalar, 2D × 2D giống matmul, nhưng với array >2D thì xử lý batch khác matmul — dễ gây bug.
• np.matmul / @: chuẩn cho ma trận, hỗ trợ batch ma trận đúng cách (batch, m, k) @ (batch, k, n) → (batch, m, n).

Khuyến nghị: với code có array nhiều chiều, dùng @ / np.matmul mặc định, chỉ dùng np.dot khi cần dot product 1D-1D.

Đừng nhầm A @ B (matmul) với A * B (element-wise, Hadamard product). Hai phép này khác hẳn:

A * B
# array([[ 5, 12],
#        [21, 32]])         — element-wise, KHÁC matmul

1. Mean per feature của batch. Batch X shape (batch_size, n_features), muốn mean từng feature trên cả batch:

batch_mean = X.mean(axis=0)   # shape (n_features,)

2. Standardization (z-score). Biến đổi mỗi feature về mean 0, std 1:

\[ z = \frac{x - \mu}{\sigma} \]

mu    = X.mean(axis=0)        # shape (n_features,)
sigma = X.std(axis=0)         # shape (n_features,)
X_std = (X - mu) / sigma      # broadcasting (B33)

3. Activation function vectorized.

ReLU: \( \mathrm{relu}(x) = \max(0, x) \).

def relu(x):
    return np.maximum(0, x)

Sigmoid: \( \sigma(x) = \dfrac{1}{1 + e^{-x}} \).

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

4. Cosine similarity. Đo độ tương đồng giữa 2 vector embedding:

\[ \mathrm{cos}(u, v) = \frac{u \cdot v}{\|u\| \cdot \|v\|} \]

def cosine_similarity(u, v):
    return np.dot(u, v) / (np.linalg.norm(u) * np.linalg.norm(v))

Giá trị trong khoảng [-1, 1]; gần 1 nghĩa là 2 vector cùng hướng. Đây là phép tính cốt lõi khi search trong vector database / RAG.

import numpy as np
import time

# ----- Activation functions vectorized -----
def relu(x):
    return np.maximum(0, x)

def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))

x = np.array([-2.0, -0.5, 0.0, 0.5, 2.0])
print(relu(x))      # [0.  0.  0.  0.5 2. ]
print(sigmoid(x))   # [0.119 0.378 0.5   0.622 0.881]

# ----- Aggregation theo axis: ma trận điểm thi -----
# 3 môn (toán, văn, anh) x 5 học sinh
scores = np.array([[ 8.5,  7.0,  9.0,  6.5,  8.0],   # toán
                   [ 7.0,  8.5,  6.5,  9.0,  7.5],   # văn
                   [ 9.0,  7.5,  8.0,  8.0,  6.0]])  # anh
# shape (3, 5): axis=0 → môn, axis=1 → học sinh

print("Mean mỗi môn   :", scores.mean(axis=1))       # axis=1 → collapse học sinh
print("Mean mỗi HS    :", scores.mean(axis=0))       # axis=0 → collapse môn
print("Max mỗi HS     :", scores.max(axis=0))
print("HS giỏi toán nhất:", scores[0].argmax())      # index trong hàng toán
print("Tổng toàn bộ   :", scores.sum())              # axis=None

# ----- Benchmark vectorized vs loop -----
n = 500_000
a = np.random.rand(n)
b = np.random.rand(n)
c = np.random.rand(n)

t0 = time.perf_counter()
result_vec = a * b + c
t_vec = time.perf_counter() - t0

t0 = time.perf_counter()
result_loop = np.zeros(n)
for i in range(n):
    result_loop[i] = a[i] * b[i] + c[i]
t_loop = time.perf_counter() - t0

print(f"vectorized: {t_vec*1000:.2f} ms")
print(f"loop:       {t_loop*1000:.2f} ms")
print(f"speedup:    {t_loop / t_vec:.0f}x")
assert np.allclose(result_vec, result_loop)

# ----- Cosine similarity -----
def cosine_similarity(u, v):
    return np.dot(u, v) / (np.linalg.norm(u) * np.linalg.norm(v))

u = np.array([1.0, 2.0, 3.0])
v = np.array([2.0, 4.0, 6.0])      # cùng hướng → 1.0
w = np.array([-1.0, -2.0, -3.0])   # ngược hướng → -1.0
print(cosine_similarity(u, v))     # 1.0
print(cosine_similarity(u, w))     # -1.0

Bài 1. Sinh ma trận M = np.random.rand(100, 5). Tính mean và std theo axis=0; xác nhận output có shape (5,). So sánh với mean / std theo axis=1 (shape (100,)) và giải thích ý nghĩa của từng output trong ngữ cảnh ma trận sample-feature.

Bài 2. Implement softmax vectorized cho vector 1D:

\[ \mathrm{softmax}(x_i) = \frac{e^{x_i}}{\sum_j e^{x_j}} \]

Test với x = np.array([1.0, 2.0, 3.0]); kết quả phải có tổng bằng 1. Gợi ý: để tránh tràn số khi x lớn, trừ x.max() trước khi np.exp (kết quả không đổi nhưng ổn định số học).

Bài 3. Tính khoảng cách Euclidean giữa 2 vector dùng vectorization: \( d(u, v) = \sqrt{\sum_i (u_i - v_i)^2} \). Test với u = np.array([1, 2, 3]), v = np.array([4, 6, 3]); kết quả mong đợi 5.0. Không dùng vòng for.

Bài 4. Cho X shape (n_samples, n_features). Viết hàm standardize(X) trả về ma trận đã chuẩn hoá (mỗi feature có mean 0, std 1). Sau khi chuẩn hoá, kiểm tra X_std.mean(axis=0) gần 0 và X_std.std(axis=0) gần 1.

Bài 5. Cho ma trận điểm scores shape (3, 50) (3 môn x 50 học sinh). Tìm: (a) học sinh có tổng điểm cao nhất, (b) môn có điểm trung bình thấp nhất, (c) số học sinh có cả 3 môn > 7.

Bài 6. Benchmark phép tính np.sqrt(a**2 + b**2) với n = 1_000_000: (a) vectorized, (b) Python for-loop với math.sqrt. In ra thời gian và tỉ số speedup.

• Vectorization = viết phép tính trên cả array; loop chạy trong C, nhanh hơn for Python vài chục đến vài trăm lần.
• Toán tử + - * / // % **, so sánh == != < > <= >=, logic bitwise & | ~ đều element-wise; cần shape khớp (hoặc broadcasting).
• Ufunc np.sin, np.cos, np.exp, np.log, np.sqrt, np.round, np.floor, np.ceil — dùng thay math.* khi làm việc với array.
• Aggregation sum, mean, std, var, min, max, median, argmin, argmax, cumsum, cumprod — gọi qua hàm np.* hoặc method a.*.
• axis=0 collapse hàng (per cột), axis=1 collapse cột (per hàng), axis=None collapse toàn bộ; keepdims=True để giữ chiều cho broadcasting.
• @ / np.matmul cho matrix multiplication (không phải element-wise); shape (m,k) @ (k,n) → (m,n). np.dot chỉ dùng cho 1D-1D dot product.
• Activation functions, standardization, cosine similarity, softmax — tất cả viết vectorized trong vài dòng.