NumPy

Broadcasting allows NumPy to operate on arrays of different shapes without copying data. The three rules: (1) If arrays have different number of dimensions, prepend 1s to the shape of the smaller-rank array. (2) Arrays with size 1 along a dimension are stretched to match the other array's size in that dimension. (3) Arrays whose sizes differ along a dimension and neither is 1 raise a ValueError. Broadcasting never actually copies data—it uses stride tricks internally, making it memory-efficient. Example: adding shape (3,1) to (1,4) yields (3,4) output without allocating intermediate arrays.

import numpy as np
a = np.array([[1],[2],[3]])   # shape (3,1)
b = np.array([10,20,30,40])  # shape (4,)
print((a + b).shape)          # (3,4) - broadcasted

A NumPy view shares the underlying data buffer with the original array; modifying one modifies the other. Slicing creates views by default. A copy has its own data buffer—changes are independent. Check with arr.base: if base is not None, it is a view. np.copy() or arr.copy() forces a copy. Fancy indexing (integer arrays, boolean masks) always returns a copy, not a view. Understanding this distinction is critical for performance: views avoid memory allocation for large arrays, but accidental in-place modification through a view can introduce subtle bugs in ML pipelines where data should remain immutable.

import numpy as np
a = np.arange(6)
b = a[2:5]       # view
b[0] = 99
print(a)         # [0 1 99 3 4 5] - original modified
c = a[2:5].copy()
c[0] = 0
print(a)         # unchanged

Strides are the number of bytes to step in each dimension to move to the next element. A C-contiguous array of shape (m, n) with dtype float64 has strides (8n, 8). Transposing simply swaps the strides tuple: shape becomes (n, m), strides become (8, 8n)—no data is copied. Similarly, reshaping only changes shape and strides when the array is contiguous. np.lib.stride_tricks.as_strided creates arbitrary views by specifying custom shapes and strides, enabling sliding window views, Toeplitz matrices, and im2col operations without data copies. Incorrect strides can cause memory access violations, so caution is needed.

import numpy as np
a = np.zeros((3,4), dtype=np.float64)
print(a.strides)       # (32, 8)  -> 4 floats * 8 bytes per row
print(a.T.strides)     # (8, 32)  -> transposed, no copy

In C-order (default in NumPy), the last index varies fastest in memory—elements of the same row are contiguous. In Fortran-order, the first index varies fastest—elements of the same column are contiguous. For matrix operations: row-wise access is fast in C-order, column-wise in F-order. np.ascontiguousarray() and np.asfortranarray() convert between layouts (copying if needed). BLAS/LAPACK routines expect F-order for some operations. When calling Fortran or C extensions, passing the wrong layout causes implicit copies. Check with arr.flags['C_CONTIGUOUS'] or arr.flags['F_CONTIGUOUS'].

import numpy as np
a = np.array([[1,2,3],[4,5,6]], order='C')
b = np.array([[1,2,3],[4,5,6]], order='F')
print(a.strides)  # (24, 8)
print(b.strides)  # (8, 16)

np.einsum uses Einstein summation notation to express tensor contractions, outer products, and transposes concisely. The subscript string maps input axes to output axes; repeated indices imply summation over that axis. ML use cases: (1) Batched matrix multiply: np.einsum('bij,bjk->bik', A, B) handles a batch of matrix multiplications efficiently. (2) Attention score computation: np.einsum('bik,bjk->bij', Q, K) computes pairwise dot products. (3) Outer product: np.einsum('i,j->ij', a, b). Compared to np.matmul or np.dot, einsum is more expressive for non-standard contractions and can be optimized with optimize='optimal' to choose the best contraction order.

import numpy as np
A = np.random.randn(32, 10, 16)  # batch, seq, dim
B = np.random.randn(32, 16, 8)   # batch, dim, out
C = np.einsum('bij,bjk->bik', A, B)  # (32,10,8)
print(C.shape)

A ufunc (universal function) operates element-wise on ndarrays, supporting broadcasting, type casting, and output argument specification. It is implemented in compiled C, making it orders of magnitude faster than Python loops. Key ufunc features: output argument (out=) writes results to a pre-allocated array avoiding allocation overhead; reduce method (np.add.reduce) performs cumulative reduction; accumulate gives running totals; outer produces an outer product. Creating custom ufuncs uses np.frompyfunc (wraps Python functions) or np.vectorize (convenience wrapper, not truly compiled). For true performance, use numba.vectorize to JIT-compile custom ufuncs to C speed.

import numpy as np
a = np.arange(1e6)
out = np.empty_like(a)
np.sqrt(a, out=out)   # in-place, avoids allocation
print(np.add.reduce(a))   # sum via ufunc reduce

Basic slicing uses integer ranges, slices, and newaxis—it returns views with no data copy, O(1) time. Fancy indexing uses integer arrays or boolean masks—it always returns copies, requiring O(k) time and memory for k selected elements. Fancy indexing is more flexible: it can select non-contiguous rows, reorder elements, and select duplicates. Boolean fancy indexing (masking) is equivalent to np.where and is the recommended pattern for conditional selection. Performance-wise, avoid fancy indexing in tight loops; instead batch all fancy index operations at once. In ML, boolean masking is ubiquitous for filtering samples satisfying constraints.

import numpy as np
a = np.arange(10)
print(a[[0, 3, 7]])          # fancy: copy, [0 3 7]
print(a[a > 5])              # boolean fancy: copy, [6 7 8 9]
print(a[2:6])                # basic slice: view

np.where has two modes. (1) np.where(condition): equivalent to np.nonzero(condition), returns a tuple of arrays of indices where condition is True—useful for finding positions satisfying a constraint. (2) np.where(condition, x, y): element-wise selects from x where condition is True, else from y. Both x and y are fully evaluated before the selection (unlike Python's conditional expression which short-circuits), so avoid side effects. Broadcasting applies to condition, x, and y. This is the NumPy vectorized equivalent of the Python ternary operator and is heavily used in ML for numerical stability guards like np.where(x > 0, x, 0) for ReLU.

import numpy as np
x = np.array([-2, -1, 0, 1, 2])
relu = np.where(x > 0, x, 0)
print(relu)  # [0 0 0 1 2]
idx = np.where(x < 0)  # indices of negatives
print(idx)   # (array([0, 1]),)

np.linalg.svd decomposes an m×n matrix A into U @ np.diag(s) @ Vh, where U is m×m orthogonal, s is min(m,n) singular values in descending order, and Vh is n×n orthogonal (conjugate transpose of V). With full_matrices=False (economy SVD), U is m×k and Vh is k×n where k=min(m,n), saving memory. Singular values are non-negative and sorted descending. The rank of A equals the number of non-zero singular values. SVD is numerically stable: condition number = s[0]/s[-1]. Applications in ML: PCA (truncate to top-k singular values), pseudoinverse, low-rank approximation, and recommendation systems.

import numpy as np
A = np.random.randn(100, 20)
U, s, Vh = np.linalg.svd(A, full_matrices=False)
print(U.shape, s.shape, Vh.shape)  # (100,20) (20,) (20,20)
A_approx = U[:,:5] @ np.diag(s[:5]) @ Vh[:5,:]  # rank-5 approx

np.linalg.solve(A, b) computes the exact solution to the square system Ax=b, requiring A to be square and invertible. np.linalg.lstsq(A, b) finds the minimum-norm least-squares solution minimizing ||Ax-b||², working for any m×n A (overdetermined, underdetermined, or rank-deficient). It uses SVD internally and returns the solution, residuals, rank, and singular values. Use lstsq when: A is not square, A is rank-deficient, the system is overdetermined (more equations than unknowns, e.g., linear regression with n_samples >> n_features), or numerical stability is a concern. lstsq is the correct approach for fitting linear models.

import numpy as np
# Linear regression: find weights minimizing ||Xw - y||^2
X = np.random.randn(100, 5)  # 100 samples, 5 features
y = np.random.randn(100)
w, res, rank, sv = np.linalg.lstsq(X, y, rcond=None)
print(w.shape)  # (5,)

np.random.default_rng() creates a Generator object using the PCG64 algorithm (cryptographically stronger than the legacy MT19937 Mersenne Twister). The new API is stateless by design—each generator is independent—making it safe for parallel use. Legacy functions like np.random.randn modify a global state, causing non-reproducible behavior when multiple threads call them. The new Generator supports more distributions, better statistical properties, and seeded reproducibility: rng = np.random.default_rng(42). Methods: rng.random(), rng.integers(), rng.standard_normal(), rng.choice(). Always use default_rng for ML experiments requiring reproducibility and portability.

import numpy as np
rng = np.random.default_rng(42)
X = rng.standard_normal((100, 10))
idx = rng.choice(100, size=20, replace=False)
print(X.shape, idx.shape)

sliding_window_view creates a read-only view of an array with an extra trailing dimension for the window elements, using stride tricks to avoid any data copy. For a 1D array of length n and window size k, the output shape is (n-k+1, k) with strides (element_stride, element_stride)—each 'row' is a window, adjacent rows overlap. Modifying the input modifies all overlapping windows. Applications in ML: efficient feature extraction (rolling statistics), 1D convolution without explicit loops, and batch processing of time series segments. For 2D (image patches), it produces shape (h-kh+1, w-kw+1, kh, kw).

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
arr = np.arange(10)
wins = sliding_window_view(arr, window_shape=3)
print(wins.shape)   # (8, 3)
print(wins)         # [[0 1 2],[1 2 3],...,[7 8 9]]

Naive softmax exp(x)/sum(exp(x)) overflows for large x and underflows for very negative x. The stable version subtracts the maximum: shift = x - x.max(axis=-1, keepdims=True); e = np.exp(shift); return e / e.sum(axis=-1, keepdims=True). Subtracting the maximum doesn't change the result mathematically (the shift cancels in numerator and denominator) but ensures all exponents are ≤ 0, preventing overflow. The sum is always ≥ 1 (the maximum term contributes exp(0)=1), preventing division by zero. This pattern is universal in ML: log-softmax uses log(e) - log(sum(e)) = shift - log(sum(exp(shift))).

import numpy as np
def softmax(x):
    x = x - x.max(axis=-1, keepdims=True)
    e = np.exp(x)
    return e / e.sum(axis=-1, keepdims=True)

scores = np.array([[1000., 1001., 1002.]])
print(softmax(scores))  # [0.09, 0.24, 0.67] - stable

np.newaxis is an alias for None used in indexing to insert a new axis of size 1, reshaping without copying data. a[:, np.newaxis] adds an axis at position 1. np.expand_dims(a, axis) is functionally equivalent but more readable for programmatic axis selection, especially when the axis is a variable. Both return views. Use np.newaxis for inline broadcasting adjustments (e.g., column_vector @ row_vector = outer product via (n,1) @ (1,m) → (n,m)). np.squeeze is the reverse—removes size-1 axes. In ML, expand_dims is preferred in model code for clarity; np.newaxis is compact for one-liner manipulations.

import numpy as np
a = np.array([1,2,3])      # shape (3,)
b = a[:, np.newaxis]       # shape (3,1)
c = np.expand_dims(a, 0)   # shape (1,3)
print(b.shape, c.shape)
# outer product via broadcasting:
print((b * a).shape)  # (3,3)

np.tensordot(A, B, axes) contracts specified axes of A and B. With axes=1, it is equivalent to np.dot for 2D arrays. With axes=2, it contracts the last two axes of A with the first two of B. With axes=([0,1],[1,2]), it contracts custom axis pairs. np.matmul is a special case restricted to 2D matrices (or last two axes of batched tensors). np.tensordot is more flexible for high-dimensional tensor operations common in deep learning (e.g., weight tensors in multi-head attention, batch outer products). It uses BLAS internally for contiguous data, achieving near-theoretical FLOPS. Equivalent to np.einsum but sometimes faster.

import numpy as np
A = np.random.randn(3,4,5)
B = np.random.randn(5,6)
C = np.tensordot(A, B, axes=1)  # contracts last axis of A with first of B
print(C.shape)  # (3,4,6)

A structured array is an ndarray with a compound dtype, where each element contains multiple named fields of potentially different types. Defined with a list of (name, type) tuples: dt = np.dtype([('name','U20'),('age','i4'),('score','f8')]). Access fields by name: arr['age']. Useful for processing tabular data with mixed types without converting everything to float64 (which would lose string fields or waste memory). Structured arrays are used to interface with C structs via ctypes, read binary file formats, and process CSV-like data before pandas. They support vectorized operations field-by-field and can be sorted by a specific field.

import numpy as np
dt = np.dtype([('name','U10'),('score','f4'),('label','i2')])
data = np.array([('Alice',0.92,1),('Bob',0.78,0)], dtype=dt)
print(data['score'])    # [0.92 0.78]
print(data[data['label']==1]['name'])  # ['Alice']

np.memmap maps a binary file on disk to a NumPy array using the OS's virtual memory system. Only the accessed pages are loaded into RAM on demand; unaccessed pages remain on disk. This enables processing arrays larger than available RAM with standard NumPy syntax. Usage: np.memmap('data.npy', dtype='float32', mode='r+', shape=(1000000, 512)). Modes: 'r' (read-only), 'r+' (read-write, modifies file), 'w+' (create/overwrite), 'c' (copy-on-write, changes not written to disk). In ML training, memmap is used for large embedding tables and pre-processed feature matrices stored on SSD.

import numpy as np
# Create a 1GB float32 array on disk
mm = np.memmap('big.dat', dtype='float32', mode='w+', shape=(1000,250000))
mm[0] = np.arange(250000, dtype='float32')
del mm  # flushes to disk
mm2 = np.memmap('big.dat', dtype='float32', mode='r', shape=(1000,250000))
print(mm2[0,:5])  # [0. 1. 2. 3. 4.]

np.fft.fft computes the 1D Discrete Fourier Transform of an array in O(n log n) using the Cooley-Tukey FFT algorithm. Output is a complex array of n frequency-domain components; the first element is the DC component (sum of signal), elements 1 to n//2-1 are positive frequencies, and n//2 to n-1 are negative frequencies (mirrored for real input). Frequencies in Hz = np.fft.fftfreq(n, d=1/sample_rate). Power spectrum = np.abs(fft)**2. Applications in ML: audio feature extraction (spectrum analysis before MFCC), periodic pattern detection in time series, and computing convolution via pointwise multiplication in frequency domain (O(n log n) vs O(n²) direct convolution).

import numpy as np
fs = 1000; t = np.linspace(0,1,fs)
signal = np.sin(2*np.pi*50*t) + np.sin(2*np.pi*120*t)
spectrum = np.fft.rfft(signal)
freqs = np.fft.rfftfreq(len(signal), 1/fs)
print(freqs[np.argsort(np.abs(spectrum))[-2:]])  # ~[50, 120]

A masked array (numpy.ma module) pairs a regular ndarray with a boolean mask of the same shape; True means the value is masked (invalid/missing). Arithmetic operations skip masked elements, and aggregations (mean, sum) ignore them. Create with np.ma.masked_array(data, mask) or np.ma.masked_where(condition, data) or np.ma.masked_invalid(data) for NaN/Inf. Masked arrays are useful for sensor data with dropout, astronomical imaging (masking bright stars), and climatological data. However, they add overhead; for ML pipelines, prefer explicit NaN handling with np.nanmean/np.nansum for performance-critical paths.

import numpy as np
data = np.array([1.,2.,np.nan,4.,np.nan])
ma = np.ma.masked_invalid(data)
print(ma.mean())   # 2.333... ignores NaN
print(ma.count())  # 3

Pairwise distances between sets A (n×d) and B (m×d) using the identity ||a-b||² = ||a||² - 2a·b + ||b||². Compute: sq_a = (A**2).sum(axis=1, keepdims=True), sq_b = (B**2).sum(axis=1), cross = A @ B.T. Then dist² = sq_a + sq_b - 2*cross (broadcasting handles shapes). Take np.sqrt for actual distances. This uses BLAS matrix multiplication (A @ B.T) for the dominant term, giving O(n·m·d) with excellent cache performance. Avoid explicit loops which are O(n·m) Python iterations plus O(d) per pair. Use np.clip(dist2, 0, None) before sqrt to handle numerical negatives.

import numpy as np
def pairwise_sq_dist(A, B):
    sq_A = (A**2).sum(1, keepdims=True)  # (n,1)
    sq_B = (B**2).sum(1)                 # (m,)
    return sq_A + sq_B - 2 * (A @ B.T)  # (n,m)
A = np.random.randn(100,32)
B = np.random.randn(200,32)
D = np.sqrt(np.clip(pairwise_sq_dist(A,B), 0, None))

np.vectorize wraps a scalar Python function to accept arrays, applying it element-wise with broadcasting support. It handles type inference and is convenient for prototyping. However, it is not a true vectorization—it internally uses a Python loop, so performance is similar to an explicit Python for-loop (no SIMD, no compiled C). Avoid it when performance matters; instead rewrite using NumPy ufuncs, boolean indexing, or numba.vectorize. Appropriate uses: wrapping a function that calls external libraries per element, handling irregular types, or applying functions with complex Python logic during experimentation. np.frompyfunc provides similar functionality with more control over output dtype.

import numpy as np
# SLOW: np.vectorize is a Python loop
f_vec = np.vectorize(lambda x: x**2 + 1)
a = np.arange(1_000_000)
# FAST: use NumPy directly
result = a**2 + 1

np.einsum(subscripts, *operands, optimize=False) by default evaluates the contraction left to right, which may be suboptimal for chains of three or more tensors. Setting optimize='optimal' uses dynamic programming to find the best contraction order, minimizing total FLOPs. optimize='greedy' uses a heuristic that is faster to compute but may not find the global optimum. For two-tensor contractions, optimize has no effect. Matters most for contractions with three+ tensors of high rank: a naive ordering can be exponentially worse. Pre-compute the optimal path with np.einsum_path and cache it, then pass path to repeated calls of the same contraction.

import numpy as np
A=np.random.randn(100,100)
B=np.random.randn(100,100)
C=np.random.randn(100,100)
path, info = np.einsum_path('ij,jk,kl->il', A, B, C, optimize='optimal')
print(info)
result = np.einsum('ij,jk,kl->il', A, B, C, optimize=path)

np.concatenate joins arrays along an existing axis; all arrays must have the same shape except along the concatenation axis. np.stack creates a new axis and stacks arrays along it; all arrays must have identical shapes. np.hstack is concatenate along axis 1 (columns) for 2D, or axis 0 for 1D. np.vstack is concatenate along axis 0 (rows), adding a new axis for 1D arrays. np.dstack stacks along axis 2 (depth). Key difference: concatenate never adds a new axis; stack always adds one. In ML, stack is used to batch samples (each sample shape (d,) → batch shape (n,d)), while concatenate combines feature matrices.

import numpy as np
a = np.ones((3,4)); b = np.ones((3,4))
print(np.concatenate([a,b], axis=0).shape)  # (6,4)
print(np.stack([a,b], axis=0).shape)        # (2,3,4)
print(np.hstack([a,b]).shape)               # (3,8)
print(np.vstack([a,b]).shape)               # (6,4)

Batch normalization standardizes each feature across the batch: mean = X.mean(axis=0, keepdims=True), std = X.std(axis=0, keepdims=True) (or use var for numerical control). Normalized: X_norm = (X - mean) / (std + epsilon). Then scale and shift: Y = gamma * X_norm + beta, where gamma and beta are learnable parameters. Epsilon (e.g., 1e-5) prevents division by zero. During inference, use running mean and variance instead of batch statistics. All operations are vectorized O(n·d) with broadcasting handling the (1,d) statistics against (n,d) batch. The keepdims=True is essential for correct broadcasting.

import numpy as np
def batch_norm(X, gamma, beta, eps=1e-5):
    mean = X.mean(axis=0, keepdims=True)
    var  = X.var(axis=0, keepdims=True)
    X_norm = (X - mean) / np.sqrt(var + eps)
    return gamma * X_norm + beta

X = np.random.randn(32, 128)  # batch_size=32, features=128
gamma = np.ones(128); beta = np.zeros(128)
out = batch_norm(X, gamma, beta)

np.argsort returns the indices that would sort the array. arr[np.argsort(arr)] == np.sort(arr). For ranking, np.argsort(np.argsort(arr)) gives the rank of each element (0-indexed). For descending sort, use np.argsort(-arr) or np.argsort(arr)[::-1]. For partial sorting (top-k), np.argpartition(arr, k) places the k smallest elements in the first k positions (unordered) in O(n) time—faster than full sort O(n log n). Then argsort on the first k gives the sorted top-k. In ML: top-k retrieval in nearest neighbor search, ranking model predictions by confidence.

import numpy as np
arr = np.array([3.,1.,4.,1.,5.,9.,2.,6.])
print(np.argsort(arr))           # indices for ascending sort
print(np.argsort(np.argsort(arr)))  # ranks
k=3
idx=np.argpartition(arr,k)[:k]  # bottom-3 in O(n)
print(arr[idx])                  # 3 smallest (unordered)

np.linalg.eig computes eigenvalues and right eigenvectors of a general square matrix. Eigenvalues may be complex. No assumption about symmetry; uses general algorithms (QZ or QR iteration). np.linalg.eigh is specialized for real symmetric or complex Hermitian matrices. It guarantees real eigenvalues, orthogonal eigenvectors, and sorts eigenvalues in ascending order. Internally uses divide-and-conquer DSYEVD (faster and more accurate than the general path). For PCA, covariance matrices are symmetric positive semi-definite—always use eigh. It is twice as fast as eig for symmetric matrices and numerically more stable. eigh also provides better accuracy for ill-conditioned symmetric systems.

import numpy as np
X = np.random.randn(100, 5)
C = X.T @ X  # symmetric positive semi-definite (5x5)
vals, vecs = np.linalg.eigh(C)  # sorted ascending, real
print(vals)      # all >= 0
print(vecs.shape)  # (5,5), orthogonal columns

Python's Global Interpreter Lock (GIL) prevents multiple threads from executing Python bytecode simultaneously, effectively making CPU-bound Python code single-threaded. However, NumPy operations release the GIL during the underlying C/BLAS computations. This means multiple threads can run BLAS matrix multiplications simultaneously, achieving true parallelism. Large np.dot, np.matmul, and np.linalg operations are already multi-threaded via OpenBLAS/MKL (controlled by OMP_NUM_THREADS). For embarrassingly parallel element-wise operations, use numpy + multiprocessing (not threading) since those operations are GIL-bound Python calls. Numba @njit(parallel=True) or CuPy (GPU) avoids the GIL entirely.

One-hot encoding converts integer class labels to binary row vectors. Efficient approach: allocate a zeros matrix of shape (n, num_classes) and use advanced indexing to set the appropriate column to 1 for each row. np.eye(num_classes)[labels] achieves this in one line: it indexes rows of the identity matrix by label value, directly returning the one-hot matrix without any explicit loop. Both are O(n·C) time and space. For sparse datasets with many classes, use scipy.sparse.csr_matrix instead. In ML training, one-hot encoding is used for categorical targets in cross-entropy loss and as input for certain embedding layers.

import numpy as np
labels = np.array([0, 2, 1, 2, 0])
num_classes = 3
one_hot = np.eye(num_classes)[labels]
print(one_hot)
# [[1,0,0],[0,0,1],[0,1,0],[0,0,1],[1,0,0]]

np.clip(a, a_min, a_max) clamps all values in array a to the range [a_min, a_max] element-wise, equivalent to min(max(a, a_min), a_max). In gradient computation: gradient clipping np.clip(grads, -threshold, threshold) prevents exploding gradients in RNNs and deep networks. In numerical stability: sigmoid inputs are clipped before log to prevent log(0): np.clip(p, 1e-7, 1-1e-7). Cross-entropy loss: -np.sum(y * np.log(np.clip(p, 1e-15, 1))). Also used in reward clipping in reinforcement learning (clip to [-1,1]). np.clip supports scalar and array bounds and works in-place with out= parameter.

import numpy as np
# Stable cross-entropy
def cross_entropy(y_true, y_pred):
    y_pred = np.clip(y_pred, 1e-15, 1-1e-15)
    return -np.mean(np.sum(y_true * np.log(y_pred), axis=1))

# Gradient clipping
grads = np.random.randn(100)
clipped = np.clip(grads, -1.0, 1.0)

np.linalg.norm(x, ord, axis) computes vector or matrix norms. For vectors: ord=None (default) = L2 norm (Euclidean); ord=1 = L1 norm (sum of abs values); ord=np.inf = max abs value; ord=-np.inf = min abs value; ord=p = Lp norm. For matrices: ord='fro' or None = Frobenius norm (sqrt of sum of squares); ord=1 = max column sum; ord=np.inf = max row sum; ord=2 = largest singular value (spectral norm). The axis parameter computes norms along a specific axis for batched inputs. In ML: L2 weight decay uses Frobenius norm, L1 regularization uses entry-wise L1 norm, and gradient norm monitoring uses L2 vector norm.

import numpy as np
W = np.random.randn(10, 20)
print(np.linalg.norm(W, 'fro'))  # Frobenius
print(np.linalg.norm(W, ord=2))  # spectral norm
vecs = np.random.randn(100, 32)
norms = np.linalg.norm(vecs, axis=1)  # (100,) row norms
normed = vecs / norms[:, np.newaxis]

np.unique(ar) returns the sorted unique elements of an array. Optional return values: return_index (indices in original where unique values first appear), return_inverse (indices to reconstruct original from unique), return_counts (occurrence counts of each unique value). Set operations between arrays: np.intersect1d(a,b), np.union1d(a,b), np.setdiff1d(a,b) (elements in a not in b), np.in1d(a,b) (boolean mask of which elements of a are in b, now prefer np.isin). All internally use sorting, O(n log n). In ML: computing class frequencies with return_counts, finding overlapping vocabularies, and filtering seen/unseen examples.

import numpy as np
labels = np.array([2,0,1,2,0,0,1,2])
uniq, inv, cnt = np.unique(labels, return_inverse=True, return_counts=True)
print(uniq)   # [0 1 2]
print(cnt)    # [3 2 3] - class frequencies
print(inv)    # indices to reconstruct labels from uniq

Im2col ('image to column') reshapes a 2D input image into a matrix where each column (or row) is a flattened receptive field patch. A convolution with kernel K then becomes a single matrix multiply: output = K_matrix @ im2col(input). This converts O(C_out × C_in × H_out × W_out × Kh × Kw) element-wise convolution into a GEMM (general matrix multiply) call, which benefits from highly optimized BLAS routines. The trade-off is memory: im2col duplicates overlapping patches, expanding memory by a factor of Kh×Kw. Deep learning frameworks (PyTorch, TensorFlow) use im2col or FFT-based convolution internally. NumPy stride tricks implement a zero-copy version using sliding_window_view.

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
def im2col(x, kh, kw):
    H,W=x.shape
    patches=sliding_window_view(x,(kh,kw))
    return patches.reshape(-1, kh*kw)  # (H-kh+1)*(W-kw+1) x kh*kw
img=np.random.randn(6,6); K=np.random.randn(3,3)
col=im2col(img,3,3)
out=(col @ K.reshape(-1)).reshape(4,4)

NumPy follows dtype promotion: when operating on two arrays of different dtypes, the result uses the 'higher' dtype sufficient to represent both. int32 + float32 → float64 (not float32, to preserve precision). float32 + float64 → float64. int64 + float32 → float64. uint8 + int8 → int16. This can cause unexpected memory expansion: a float32 neural network accidentally mixed with float64 constants creates float64 intermediates, doubling memory usage and slowing computation. Always explicitly set dtypes: np.array([1.0], dtype=np.float32). Use astype() to control precision. NumPy 2.0 introduced stricter value-based promotion to reduce surprise upcasting.

import numpy as np
a = np.array([1], dtype=np.float32)
b = np.array([2.0])  # float64 by default
print((a+b).dtype)   # float64 - upcasted!
# Fix:
c = np.array([2.0], dtype=np.float32)
print((a+c).dtype)   # float32

The Gram matrix K[i,j] = k(x_i, x_j) for a kernel function k. For the linear kernel k(x,y) = x·y, the Gram matrix is simply X @ X.T (n×d times d×n = n×n), computed in one BLAS SGEMM call. For the RBF kernel k(x,y) = exp(-γ||x-y||²), compute pairwise squared distances using the identity ||x-y||² = ||x||²-2x·y+||y||², then apply exp. For the polynomial kernel k(x,y) = (x·y + c)^d, compute (X @ X.T + c)**d element-wise. All O(n²·d) time. The Gram matrix is symmetric positive semi-definite, enabling Cholesky decomposition and efficient linear solves via np.linalg.cholesky.

import numpy as np
def rbf_gram(X, gamma=1.0):
    sq = (X**2).sum(1)
    dist2 = sq[:,None] + sq[None,:] - 2*(X@X.T)
    return np.exp(-gamma * np.clip(dist2, 0, None))
X = np.random.randn(200, 50)
K = rbf_gram(X)
print(K.shape)  # (200, 200)

np.testing provides assertion functions for comparing arrays with tolerance, essential for testing numerical code. Key functions: np.testing.assert_allclose(actual, desired, rtol=1e-7, atol=0) checks relative and absolute tolerances simultaneously—use rtol for values of varying magnitude, atol for values near zero. assert_array_equal checks exact equality (use for integer arrays). assert_array_less checks element-wise ordering. assert_almost_equal checks to a specified decimal place. In ML: use assert_allclose to verify forward pass outputs match reference implementations, check gradient computations against finite differences (numerical gradient checking), and validate that refactored code produces identical results within floating-point precision.

import numpy as np
def my_sigmoid(x):
    return 1 / (1 + np.exp(-np.clip(x,-500,500)))
def ref_sigmoid(x):
    return np.where(x>=0, 1/(1+np.exp(-x)), np.exp(x)/(1+np.exp(x)))
x = np.linspace(-10,10,1000)
np.testing.assert_allclose(my_sigmoid(x), ref_sigmoid(x), rtol=1e-6)

np.polynomial (specifically np.polynomial.polynomial) uses power series representation with coefficients from lowest to highest degree—opposite to np.polyfit which uses highest to lowest. np.polynomial.polynomial.polyfit(x, y, deg) fits a polynomial of given degree to data via least squares. polyval(c, x) evaluates the polynomial. The Polynomial class provides object-oriented access: p = Polynomial.fit(x, y, deg); p(new_x) evaluates. For numerical stability, polynomial fitting in the numpy.polynomial subpackage uses Chebyshev or Legendre basis internally. np.polyfit uses the standard basis, which can be ill-conditioned for high degrees. Prefer the polynomial subpackage for degrees above 5.

import numpy as np
from numpy.polynomial import polynomial as P
x = np.linspace(-1, 1, 50)
y = np.sin(np.pi*x) + 0.1*np.random.randn(50)
c, stats = P.polyfit(x, y, 5, full=True)
print('coeffs:', c)
pred = P.polyval(x, c)

np.sum(a, axis) and np.add.reduce(a, axis) produce identical results—np.sum internally calls np.add.reduce. The distinction is semantic: reduce is lower-level, closer to the ufunc API. The axis parameter is critical: axis=0 reduces over rows (result shape loses axis 0), axis=1 reduces over columns, axis=-1 reduces over the last axis. axis=None (default) reduces all elements. keepdims=True preserves the reduced axis as size 1, enabling broadcasting of the result back against the original. Common mistake: computing batch statistics without keepdims, causing broadcast shape errors. np.sum is numerically stable for float64; for float32, use np.float64 accumulation with keepdims or use Kahan summation for extreme precision.

import numpy as np
X = np.random.randn(32, 10, 5)  # batch, seq, features
mu = X.mean(axis=(0,1), keepdims=True)  # (1,1,5)
X_norm = X - mu  # broadcasts correctly
print(X_norm.shape)  # (32,10,5)

For a multivariate Gaussian N(μ, Σ), sampling requires drawing z ~ N(0, I) and computing x = μ + L @ z, where L is the lower Cholesky factor of Σ (Σ = L @ L.T). Cholesky decomposition is computed in O(n³) and is twice as fast as LU for positive definite matrices. L @ z adds the correct covariance since Cov(Lz) = L @ I @ L.T = Σ. np.linalg.cholesky(Sigma) returns L. This avoids computing the matrix square root directly. Also used for: solving Σx=b via forward/back substitution (scipy.linalg.cho_solve), computing log-det Σ = 2·sum(log(diag(L))), and Gaussian process inference.

import numpy as np
mu = np.zeros(5)
A = np.random.randn(5,5); Sigma = A.T @ A + np.eye(5)*0.1
L = np.linalg.cholesky(Sigma)
z = np.random.randn(1000, 5)
samples = mu + z @ L.T  # shape (1000, 5)
print(np.cov(samples.T))  # should approximate Sigma

np.diff(a, n=1, axis=-1) computes the n-th order discrete difference along the given axis. First-order difference (default): diff[i] = a[i+1] - a[i], reducing array length by 1 per order. Used in time series: first difference removes linear trends (stationarization), second difference removes quadratic trends. np.cumsum is the inverse of np.diff. For financial data, np.diff(prices)/prices[:-1] gives percentage returns. np.gradient is more accurate for non-uniform spacing and computes derivatives using central differences at interior points. np.diff with prepend/append arguments adds boundary values, preserving array length for in-place replacements.

import numpy as np
prices = np.array([100., 105., 102., 110., 108.])
returns = np.diff(prices) / prices[:-1]  # % returns
print(returns)  # [0.05, -0.028, 0.078, -0.018]
print(np.diff(prices, n=2))  # second differences

np.searchsorted(a, v, side='left') returns the index where v should be inserted in the sorted array a to maintain order. side='left' gives the leftmost valid position; side='right' gives the rightmost (one past the last equal element). Time is O(log n) per query, O(m log n) for m queries. For m >> 1, all queries are batched in a single NumPy call. Uses: (1) quantile computation: np.searchsorted(sorted_data, quantile_value). (2) Histogram computation: bucketize data into bins. (3) Inverse CDF sampling. (4) Online lookup: maintain a sorted array and use searchsorted for O(log n) membership test. More cache-friendly than Python bisect for large arrays.

import numpy as np
a = np.array([1,3,5,7,9,11])
print(np.searchsorted(a, 6))         # 3 (insert before 7)
print(np.searchsorted(a, 5,'right')) # 3 (after 5)
# Bucketize: which bin does each value fall into?
bins = np.array([0,25,50,75,100])
vals = np.array([10,40,60,80])
print(np.searchsorted(bins, vals))   # [1 2 3 4]

NumPy links against BLAS (Basic Linear Algebra Subprograms) and LAPACK (Linear Algebra PACKage) for dense linear algebra. BLAS provides three levels: Level 1 (vector ops, O(n)), Level 2 (matrix-vector, O(n²)), Level 3 (matrix-matrix, O(n³)). np.matmul and np.dot call BLAS DGEMM (Level 3), which exploits cache blocking, SIMD, and multi-threading to achieve near-theoretical FLOPS. LAPACK provides higher-level routines: eigendecomposition, SVD, Cholesky, LU. OpenBLAS (default open-source) and Intel MKL (commercial, often 2-3× faster) are the common backends. Check with np.show_config(). For float32 computations, use MKL on Intel or cuBLAS on GPU for maximum throughput.

import numpy as np
# Check BLAS/LAPACK backend
np.show_config()  # prints BLAS/LAPACK info
# Large matmul uses BLAS DGEMM automatically
A = np.random.randn(1000, 1000).astype(np.float32)
B = np.random.randn(1000, 1000).astype(np.float32)
C = A @ B  # dispatched to SGEMM

np.pad adds values around the border of an array. Signature: np.pad(array, pad_width, mode). pad_width specifies (before, after) padding per axis. Modes: 'constant' (default, value=0), 'reflect' (mirror), 'edge' (replicate border), 'wrap' (periodic). In CNNs, zero-padding (mode='constant') before convolution preserves spatial dimensions: with kernel k×k, pad each side by (k-1)//2 to maintain the same output shape as input. 'same' padding in frameworks corresponds to np.pad with mode='constant'. Reflect padding is used in image stylization to avoid edge artifacts. Batched images use pad_width=((0,0),(p,p),(p,p),(0,0)) for (N,H,W,C) tensors.

import numpy as np
img = np.random.randn(28, 28)
padded = np.pad(img, pad_width=2, mode='constant', constant_values=0)
print(padded.shape)  # (32, 32)
# Batch padding (N,H,W,C):
batch = np.random.randn(32, 28, 28, 3)
padded_batch = np.pad(batch, ((0,0),(2,2),(2,2),(0,0)))

np.gradient computes derivatives using central differences at interior points: df/dx[i] = (f[i+1]-f[i-1])/(x[i+1]-x[i-1]). At boundaries, it uses forward/backward differences to maintain full-length output. For non-uniform spacing, pass the coordinate array as a second argument: np.gradient(f, x). For 2D arrays, np.gradient returns a list of gradients per axis. Unlike np.diff (which reduces array length by 1), np.gradient preserves the original shape. Accuracy: central differences are O(h²) accurate vs O(h) for one-sided. Used in image processing (Sobel filter approximation) and finite-difference optimization of discrete functions.

import numpy as np
x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x)
dydx = np.gradient(y, x)  # numerical derivative
print(np.allclose(dydx, np.cos(x), atol=1e-3))  # True

np.repeat(a, repeats, axis) repeats each element or row/column individually: np.repeat([1,2,3], 2) = [1,1,2,2,3,3]. np.tile(a, reps) repeats the whole array as tiles: np.tile([1,2,3], 2) = [1,2,3,1,2,3]. For 2D arrays: np.repeat(A, 2, axis=0) doubles each row; np.tile(A, (2,1)) stacks A twice vertically. Both return copies. In ML: np.tile creates repeated feature matrices for batch operations; np.repeat is used for element-level upsampling (as in nearest-neighbor upsampling in image super-resolution). Neither produces a view—use broadcasting instead when possible for memory efficiency.

import numpy as np
print(np.repeat([1,2,3], 2))   # [1 1 2 2 3 3]
print(np.tile([1,2,3], 2))     # [1 2 3 1 2 3]
A = np.array([[1,2],[3,4]])
print(np.repeat(A, 2, axis=0)) # doubles each row

np.linalg.inv supports batched inputs: for an array of shape (..., n, n), it inverts each n×n matrix. For example, a (batch, n, n) array is inverted to (batch, n, n) by calling np.linalg.inv(A). Internally uses LAPACK DGETRI for each matrix in the batch. Better alternative for solving Ax=b (multiple b vectors): np.linalg.solve(A, B) avoids forming the inverse and is more numerically stable. For banded, triangular, or SPD structure, use scipy.linalg.solve_triangular or np.linalg.cholesky + solve. Avoid explicit inversion for numerical stability—direct solves are preferable whenever possible.

import numpy as np
# Batched inversion
A = np.random.randn(10, 4, 4)  # batch of 10 4x4 matrices
A = A @ A.swapaxes(-2,-1) + np.eye(4)  # make PD
A_inv = np.linalg.inv(A)  # (10, 4, 4)
# Better: solve instead of invert
b = np.random.randn(10, 4)
x = np.linalg.solve(A, b)  # (10, 4)

np.roll(a, shift, axis) shifts array elements cyclically: elements shifted off one end reappear at the other. np.roll([1,2,3,4,5], 2) = [4,5,1,2,3]. For 2D images, np.roll(img, shift, axis=0) shifts rows (vertical), axis=1 shifts columns (horizontal). In circular convolution: convolve(a, b) = IFFT(FFT(a)*FFT(b)), where multiplication corresponds to circular (periodic) convolution. np.roll implements the circular shift operator directly. Data augmentation use: randomly shift images with np.roll to create invariant representations. np.roll creates a copy for non-zero shifts; use stride tricks for zero-copy rolling if needed.

import numpy as np
arr = np.array([1,2,3,4,5])
print(np.roll(arr, 2))   # [4 5 1 2 3]
print(np.roll(arr, -1))  # [2 3 4 5 1]
# Image shift:
img = np.random.randn(100, 100)
shifted = np.roll(np.roll(img, 5, axis=0), -3, axis=1)

np.linalg.matrix_power(A, n) computes A^n using repeated squaring in O(d³ log n) time for d×d matrix A. For n=0, returns identity; n<0 returns (A^{-1})^|n|; n=1 returns A. Applications: (1) Markov chain long-run distributions: A^T converges to stationary distribution for large T. (2) Computing Fibonacci numbers: [[1,1],[1,0]]^n gives F(n) in O(log n). (3) Graph adjacency powers: A^k[i][j] = number of paths of length k from i to j. (4) PageRank power iteration is equivalent to repeated matrix-vector multiplication. For large sparse matrices, scipy.sparse.linalg.expm_multiply is more efficient than dense matrix powers.

import numpy as np
# Fibonacci via matrix power
F = np.array([[1,1],[1,0]])
n = 10
result = np.linalg.matrix_power(F, n)
print(result[0,1])  # F(10) = 55
# Markov chain steady state
P = np.array([[0.7,0.3],[0.4,0.6]])
print(np.linalg.matrix_power(P, 100))  # converges

For 1D arrays a (length m) and b (length n): np.dot(a,b) = scalar inner product (requires m==n); np.inner(a,b) = same as dot for 1D; np.outer(a,b) = m×n matrix (outer product, no shape requirement); np.matmul(a,b) = same as dot for 1D. For 2D arrays, differences emerge: np.dot(A,B) is matrix multiply; np.inner(A,B) computes inner products of rows (result[i,j]=A[i]·B[j], equivalent to A@B.T); np.matmul treats the last two dims as matrix multiply and broadcasts over batch dims; np.dot does not broadcast over batch dims. Always prefer np.matmul (@) for matrix multiplication—it has explicit broadcasting semantics and is unambiguous.

import numpy as np
a = np.array([1.,2.,3.])
b = np.array([4.,5.,6.])
print(np.dot(a,b))    # 32.0 (inner product)
print(np.outer(a,b))  # 3x3 matrix
print(np.inner(a,b))  # 32.0 (same as dot for 1D)
A = np.random.randn(2,3,4)
B = np.random.randn(2,4,5)
print((A@B).shape)    # (2,3,5) - batch matmul

np.histogram(data, bins) computes the frequency distribution: returns (counts, bin_edges). For continuous features, histograms reveal distribution shape (skewness, multimodality). np.histogram2d(x, y, bins) creates a 2D joint histogram—useful for visualizing feature correlations and marginal distributions. The density=True option normalizes counts to form a probability density. In ML: histograms guide feature preprocessing decisions (log-transform for skewed features, clipping for heavy tails); 2D histograms detect spurious correlations; activation histograms during training reveal dead neurons (ReLU) or saturation (sigmoid). np.digitize bins new data points consistently with the same edges.

import numpy as np
data = np.random.randn(10000)
counts, edges = np.histogram(data, bins=50)
print(counts.shape, edges.shape)  # (50,) (51,)
# Normalize to PDF:
density, edges = np.histogram(data, bins=50, density=True)
print(np.sum(density * np.diff(edges)))  # ~1.0

np.shares_memory(a, b) returns True if arrays a and b share any memory (i.e., one is a view of the other or they overlap). This is stronger than checking a.base is b (which only checks the direct base). Important for debugging: if you modify a NumPy array and unexpectedly see changes elsewhere in your code, shares_memory reveals hidden aliasing. Frequently encountered with: slicing (a[2:5] shares memory with a), reshape (may share if contiguous), transpose (always shares), and astype if no copy needed. In ML pipelines, unintended in-place modification of cached activations is a common source of gradient bugs. Use a.copy() to guarantee independence.

import numpy as np
a = np.arange(10)
b = a[2:7]     # view
c = a.copy()   # copy
d = a.T        # view
print(np.shares_memory(a, b))  # True
print(np.shares_memory(a, c))  # False
print(np.shares_memory(a, d))  # True