Why ML uses NumPy

• Represents data as arrays: vectors, matrices, tensors.
• Fast elementwise ops without Python loops (vectorization).
• Broadcasting makes shape math ergonomic.
• Most ML frameworks (PyTorch, TF, JAX, scikit-learn) accept or return NumPy arrays.
• The universal interchange format for numeric data in Python.

Key ideas

Array creation

NumPy provides many constructors. Use np.random.default_rng(seed) (the modern API) rather than the legacy np.random.seed() for reproducible random arrays.

import numpy as np

# From Python sequences
a = np.array([1, 2, 3])                    # 1-D, dtype=int64
b = np.array([[1.0, 2.0], [3.0, 4.0]])    # 2-D, dtype=float64

# Built-in constructors
np.zeros((3, 4))          # all zeros, shape (3,4)
np.ones((2, 3))           # all ones
np.eye(4)                 # 4x4 identity matrix
np.arange(0, 10, 2)      # [0 2 4 6 8]
np.linspace(0, 1, 5)     # [0.   0.25 0.5  0.75 1.  ]
np.full((2, 2), 7)        # [[7 7] [7 7]]

# Random
rng = np.random.default_rng(seed=42)   # preferred modern API
rng.standard_normal((3, 3))            # Gaussian N(0,1)
rng.integers(0, 10, size=(4,))         # random ints in [0,10)
rng.uniform(0, 1, size=(100,))

dtypes

NumPy infers dtype from the input. ML code almost always wants float32 to match GPU memory layout; the NumPy default of float64 doubles memory usage and slows data transfer to PyTorch/TF.

import numpy as np

# Specify dtype explicitly
a = np.array([1, 2, 3], dtype=np.float32)
b = np.zeros((4,), dtype=np.int16)

# Check and cast
print(a.dtype)              # float32
c = a.astype(np.float64)   # cast to float64

# Common dtypes in ML
# float32  – standard GPU/CPU training dtype
# float64  – NumPy default; use when precision matters
# int32    – class labels, indices
# bool_    – masks

Shapes, indexing, and slicing

Shape bugs are the most common NumPy mistake. Print .shape liberally during development and distinguish between rank-1 arrays (n,) and column/row vectors (n,1) / (1,n).

import numpy as np

a = np.arange(24).reshape(2, 3, 4)   # shape (2,3,4)

# Indexing
a[0]            # first slice along axis 0 → shape (3,4)
a[0, 1, 2]     # scalar element
a[:, 1, :]     # all along axis 0 and 2, index 1 on axis 1

# Slicing
a[0, :2, 1:]   # rows 0-1, cols 1-3 of first block

# Boolean indexing
x = np.array([10, 20, 30, 40])
mask = x > 15
x[mask]        # [20 30 40]
x[x % 20 == 0] # [20 40]

# Fancy indexing
idx = np.array([0, 2])
x[idx]         # [10 30]

# Shape manipulation
a.reshape(6, 4)        # same data, new shape
a.flatten()            # always returns a copy
a.ravel()              # returns view when possible
a.T                    # transpose
np.expand_dims(x, 0)  # (4,) -> (1,4)
x[:, np.newaxis]       # (4,) -> (4,1)

Vectorized operations and aggregations

Every arithmetic operator and universal function (ufunc) works elementwise on arrays. Aggregations accept an axis argument — axis 0 collapses rows, axis 1 collapses columns.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0])

# Elementwise arithmetic
x + 10          # [11. 12. 13. 14.]
x * x           # [ 1.  4.  9. 16.]
np.sqrt(x)      # [1.   1.41 1.73 2.  ]
np.exp(x)
np.log(x)

# Aggregations
x.sum()         # 10.0
x.mean()        # 2.5
x.std()         # 1.118...
x.min(), x.max()
x.argmin(), x.argmax()   # index of min/max

# Axis-wise aggregation (2-D example)
m = np.array([[1, 2, 3], [4, 5, 6]])
m.sum(axis=0)   # sum each column → [5 7 9]
m.sum(axis=1)   # sum each row    → [ 6 15]
m.mean(axis=0)  # [2.5 3.5 4.5]

Broadcasting — rules and pitfalls

import numpy as np

# Rule: align shapes from the right; size-1 dims are stretched.

a = np.ones((3, 4))     # shape (3,4)
b = np.ones((4,))       # shape   (4,) -> broadcast to (3,4)
(a + b).shape           # (3,4)  ✓

# Outer product via broadcasting
col = np.array([[1], [2], [3]])   # (3,1)
row = np.array([10, 20, 30])      # (3,)  treated as (1,3)
col * row   # (3,3) — outer product without np.outer

# Common pitfall: (n,) vs (n,1) vs (1,n)
x = np.arange(4)        # shape (4,)
x + x[:, np.newaxis]    # (4,4) — outer sum, intended?
x + x                   # (4,)  — elementwise, different

Linear algebra (np.linalg)

Use the @ operator for matrix multiplication — it is cleaner than np.dot and works on stacks of matrices. Prefer np.linalg.solve(A, b) over computing inv(A) @ b — it is faster and numerically more stable.

import numpy as np

A = np.array([[1., 2.], [3., 4.]])
b = np.array([5., 6.])

# Matrix multiply
A @ A               # matrix product (preferred over np.dot)
np.dot(A, b)        # also works

# Solve linear system Ax = b
x = np.linalg.solve(A, b)

# Decompositions
U, S, Vh = np.linalg.svd(A)     # singular value decomposition
Q, R     = np.linalg.qr(A)      # QR decomposition
vals, vecs = np.linalg.eig(A)   # eigenvalues / eigenvectors

# Norms
np.linalg.norm(b)           # L2 norm (default)
np.linalg.norm(b, ord=1)    # L1 norm
np.linalg.norm(A, ord='fro') # Frobenius norm

# Other useful ops
np.linalg.det(A)    # determinant
np.linalg.inv(A)    # inverse (prefer solve over inv @ b)
np.linalg.matrix_rank(A)

Common ML patterns in NumPy

Implementing these from scratch is the best way to understand what frameworks do internally — and to catch shape/dtype bugs before they hide behind autograd.

import numpy as np

# ── Sigmoid ──────────────────────────────────────────────────
def sigmoid(z):
    return 1 / (1 + np.exp(-z))

# ── Softmax (numerically stable) ─────────────────────────────
def softmax(z):
    e = np.exp(z - z.max(axis=-1, keepdims=True))
    return e / e.sum(axis=-1, keepdims=True)

# ── Binary cross-entropy loss ────────────────────────────────
def bce_loss(y_pred, y_true, eps=1e-7):
    y_pred = np.clip(y_pred, eps, 1 - eps)
    return -np.mean(y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred))

# ── Cosine similarity ────────────────────────────────────────
def cosine_sim(a, b):
    return np.dot(a, b) / (np.linalg.norm(a) * np.linalg.norm(b))

# ── Z-score normalisation ────────────────────────────────────
def z_norm(X):
    return (X - X.mean(axis=0)) / (X.std(axis=0) + 1e-8)

# ── Min-max normalisation ────────────────────────────────────
def minmax(X):
    mn, mx = X.min(axis=0), X.max(axis=0)
    return (X - mn) / (mx - mn + 1e-8)

# ── Train / val / test split (no sklearn) ───────────────────
def split(X, y, val=0.1, test=0.1, seed=0):
    rng = np.random.default_rng(seed)
    idx = rng.permutation(len(X))
    n_test = int(len(X) * test)
    n_val  = int(len(X) * val)
    return (X[idx[n_test+n_val:]], y[idx[n_test+n_val:]],
            X[idx[n_test:n_test+n_val]], y[idx[n_test:n_test+n_val]],
            X[idx[:n_test]], y[idx[:n_test]])

Useful functions quick-reference

Saving and loading arrays

Prefer .npy / .npz over CSV — they preserve dtype and shape and are an order of magnitude faster to read and write.

import numpy as np

arr = np.arange(100).reshape(10, 10)

# Single array
np.save("arr.npy", arr)
loaded = np.load("arr.npy")

# Multiple arrays in one file
np.savez("data.npz", X=arr, y=arr[:, 0])
data = np.load("data.npz")
X, y = data["X"], data["y"]

# Compressed (good for large arrays)
np.savez_compressed("data_c.npz", X=arr)

# Text (CSV-like, slower, loses dtype)
np.savetxt("arr.csv", arr, delimiter=",")
loaded_txt = np.loadtxt("arr.csv", delimiter=",")

Interoperability

NumPy arrays move between PyTorch, Pandas, and Pillow with near-zero overhead thanks to the buffer protocol. PyTorch tensors and NumPy arrays share memory when the tensor is on CPU — mutating one mutates the other.

import numpy as np

# NumPy <-> PyTorch
import torch
t = torch.from_numpy(arr)   # zero-copy when possible
a = t.numpy()               # back to NumPy (CPU tensors only)

# NumPy <-> Pandas
import pandas as pd
df = pd.DataFrame(arr, columns=[f"f{i}" for i in range(arr.shape[1])])
back = df.to_numpy()        # or df.values

# NumPy <-> PIL image
from PIL import Image
img = Image.fromarray(np.uint8(np.random.randint(0,255,(64,64,3))))
arr_img = np.array(img)     # HxWxC uint8

Common pitfalls — in depth