ML Coding (from scratch)

The normal equations solution minimizes the mean squared error: w = (X^TX)^{-1}X^Ty. Steps: form the Gram matrix X^TX (O(nd²)), compute X^Ty (O(nd)), solve the linear system with Cholesky or LU (O(d³)). Total: O(nd² + d³). For n >> d, the dominant term is O(nd²). For large d, use gradient descent or QR decomposition (more stable than directly inverting X^TX). The pseudoinverse solution X^+y is preferred for rank-deficient X. Adding regularization (ridge): w = (X^TX + λI)^{-1}X^Ty—guaranteed invertible for λ>0.

import numpy as np
def linear_regression(X, y, lam=0.0):
    n, d = X.shape
    A = X.T @ X + lam * np.eye(d)
    b = X.T @ y
    return np.linalg.solve(A, b)  # more stable than inv

X = np.random.randn(100, 5)
y = X @ np.array([1.,2.,3.,4.,5.]) + np.random.randn(100)
w = linear_regression(X, y, lam=1e-4)
print(w)

Logistic regression models P(y=1|x) = σ(w^Tx+b), σ(z)=1/(1+e^{-z}). Binary cross-entropy loss: L = -[y·log(p) + (1-y)·log(1-p)] averaged over samples. Gradients: ∂L/∂w = X^T(p-y)/n, ∂L/∂b = mean(p-y). Update with gradient descent. Numerically stable sigmoid: for z<0 use e^z/(1+e^z), else 1/(1+e^{-z}). Use np.clip on p to prevent log(0). Converges to maximum likelihood estimate under logistic model. L2 regularization adds λ||w||² to loss, gradient contribution +2λw.

import numpy as np
def sigmoid(z): return 1/(1+np.exp(-np.clip(z,-500,500)))
class LogReg:
    def __init__(self, lr=0.1, iters=1000, lam=0.):
        self.lr,self.iters,self.lam=lr,iters,lam
    def fit(self,X,y):
        n,d=X.shape; self.w=np.zeros(d); self.b=0.
        for _ in range(self.iters):
            p=sigmoid(X@self.w+self.b)
            dw=(X.T@(p-y))/n+2*self.lam*self.w
            db=(p-y).mean()
            self.w-=self.lr*dw; self.b-=self.lr*db
    def predict_proba(self,X): return sigmoid(X@self.w+self.b)

K-means minimizes the within-cluster sum of squared distances. Algorithm: (1) Initialize k centroids (random data points or k-means++). (2) Assignment step: assign each point to the nearest centroid—O(nkd) per iteration. (3) Update step: recompute centroids as mean of assigned points—O(nd). (4) Repeat until centroids don't change (convergence) or max iterations reached. Convergence is guaranteed because the total inertia (sum of squared distances) decreases monotonically and there are finitely many assignments. k-means++ initialization chooses each centroid with probability proportional to distance² from nearest chosen centroid, improving solution quality by O(log k) approximation guarantee.

import numpy as np
def kmeans(X, k, max_iter=300, seed=42):
    rng=np.random.default_rng(seed)
    centroids=X[rng.choice(len(X),k,replace=False)]
    for _ in range(max_iter):
        dists=np.linalg.norm(X[:,None]-centroids[None,:],axis=2)
        labels=dists.argmin(axis=1)
        new_c=np.array([X[labels==j].mean(0) for j in range(k)])
        if np.allclose(new_c,centroids): break
        centroids=new_c
    return labels, centroids

k-NN stores all training data and classifies a query point by the majority vote of its k nearest neighbors. Prediction: compute distances from query to all training points (O(nd)), find k smallest (O(n) with partial sort), take majority class. No training phase—k-NN is a lazy learner. Time: O(nd) per query; O(1) 'training'. Space: O(nd) for storing data. For large n, KD-trees (O(log n) query for low d) or Ball trees reduce prediction time. Euclidean distance is standard but any valid metric works. k controls the bias-variance tradeoff: small k = low bias, high variance; large k = high bias, low variance.

import numpy as np
from scipy.stats import mode
class KNN:
    def __init__(self,k=5): self.k=k
    def fit(self,X,y): self.X_tr=X; self.y_tr=y
    def predict(self,X):
        preds=[]
        for x in X:
            dists=np.linalg.norm(self.X_tr-x,axis=1)
            idx=np.argpartition(dists,self.k)[:self.k]
            preds.append(mode(self.y_tr[idx],keepdims=False).mode)
        return np.array(preds)

Gaussian Naive Bayes assumes features are conditionally independent given the class and follow Gaussian distributions per class. Training: compute class priors P(y=c), per-class feature means μ_c and variances σ_c² from training data. Prediction: compute log P(y=c|x) = log P(y=c) + Σlog N(xⱼ; μ_cj, σ_cj²), choose class with highest log-posterior. Using log space prevents underflow from multiplying many small probabilities. Log Gaussian density: -0.5·log(2πσ²) - (x-μ)²/(2σ²). Var smoothing (add ε to all variances) prevents division by zero for zero-variance features. Despite strong independence assumption, Naive Bayes often performs well for text classification.

import numpy as np
class GaussianNB:
    def fit(self,X,y):
        self.classes=np.unique(y)
        self.priors,self.means,self.vars={},{},{}
        for c in self.classes:
            Xc=X[y==c]; self.priors[c]=len(Xc)/len(X)
            self.means[c]=Xc.mean(0); self.vars[c]=Xc.var(0)+1e-9
    def predict(self,X):
        log_posts=[]
        for c in self.classes:
            lp=np.log(self.priors[c])
            lp+=(-0.5*np.log(2*np.pi*self.vars[c])-(X-self.means[c])**2/(2*self.vars[c])).sum(1)
            log_posts.append(lp)
        return self.classes[np.argmax(np.array(log_posts),axis=0)]

A decision tree recursively splits data to maximize information gain (IG = H(parent) - weighted average H(children)), where H is entropy. Training: at each node, try all features and split values; choose the split maximizing IG. Stop when max_depth reached, min_samples_split not met, or all samples have same class. Prediction: traverse the tree following the feature conditions until a leaf, return majority class. Time: O(n·d·log n) for training (n splits per level, log n levels), O(d·log n) per prediction. Memory: O(2^depth) nodes in the worst case. Gini impurity is an alternative to entropy: G = 1 - Σpᵢ².

import numpy as np
def entropy(y):
    ps=np.bincount(y)/len(y)
    return -np.sum(ps[ps>0]*np.log2(ps[ps>0]))
def info_gain(parent,left,right):
    n=len(parent)
    return entropy(parent)-(len(left)/n*entropy(left)+len(right)/n*entropy(right))
def best_split(X,y):
    best_ig,best_feat,best_thr=-1,None,None
    for f in range(X.shape[1]):
        for thr in np.unique(X[:,f]):
            l,r=y[X[:,f]<=thr],y[X[:,f]>thr]
            if len(l)>0 and len(r)>0:
                ig=info_gain(y,l,r)
                if ig>best_ig: best_ig,best_feat,best_thr=ig,f,thr
    return best_feat,best_thr

PCA steps: (1) Center data X: X -= X.mean(0). (2) Compute SVD: U, s, Vt = SVD(X). (3) Principal components are rows of Vt; scores are U * s (or X @ Vt.T). (4) Variance explained by component k: s[k]² / Σs². Choose number of components via: (a) explained variance threshold (e.g., keep components summing to ≥ 95% variance); (b) scree plot elbow; (c) cross-validation on downstream task performance. Reconstruct with X_approx = scores[:,:k] @ Vt[:k]. SVD is preferred over eigendecomposition of X^TX for numerical stability when n < d. Whitening: divide scores by s (or singular values) to equalize variance.

import numpy as np
def pca(X, n_components):
    X_c = X - X.mean(0)
    U, s, Vt = np.linalg.svd(X_c, full_matrices=False)
    explained = s**2 / (s**2).sum()
    print('Explained:', explained[:n_components].sum())
    return X_c @ Vt[:n_components].T, Vt[:n_components]

X = np.random.randn(200, 20)
scores, components = pca(X, 5)

Forward: compute batch mean μ = mean(X, axis=0), variance σ² = var(X, axis=0), normalize X_hat = (X-μ)/√(σ²+ε), scale/shift Y = γX_hat + β. Cache {X, μ, σ², X_hat, γ, ε} for backward. Backward: dγ = sum(dY * X_hat, axis=0); dβ = sum(dY, axis=0); dX_hat = dY * γ; dσ² = sum(dX_hat * (X-μ) * -0.5*(σ²+ε)^{-1.5}, axis=0); dμ = sum(-dX_hat/√(σ²+ε), axis=0) + dσ²·sum(-2(X-μ),axis=0)/n; dX = dX_hat/√(σ²+ε) + dσ²·2(X-μ)/n + dμ/n. During inference, use running statistics (exponential moving averages computed during training).

import numpy as np
def bn_forward(X, gamma, beta, eps=1e-5):
    mu=X.mean(0); var=X.var(0)
    Xh=(X-mu)/np.sqrt(var+eps)
    return gamma*Xh+beta, (X,mu,var,Xh,gamma,eps)
def bn_backward(dY, cache):
    X,mu,var,Xh,gamma,eps=cache; n=X.shape[0]
    dg=(dY*Xh).sum(0); db=dY.sum(0)
    dXh=dY*gamma
    dvar=(dXh*(X-mu)*(-0.5)*(var+eps)**(-1.5)).sum(0)
    dmu=(-dXh/np.sqrt(var+eps)).sum(0)+dvar*(-2*(X-mu)).sum(0)/n
    dX=dXh/np.sqrt(var+eps)+dvar*2*(X-mu)/n+dmu/n
    return dX,dg,db

Dropout randomly zeros activations with probability p during training, scaling remaining activations by 1/(1-p) (inverted dropout). This ensures expected output magnitude matches inference, where no dropout is applied. Inverted dropout implementation: generate mask = (random uniform < (1-p)) of same shape; multiply input by mask, then divide by (1-p). Store mask for backward pass: gradient multiplied by same mask, divided by (1-p). During inference: simply return input unchanged (no masking, no scaling). Non-inverted dropout (scale during inference by 1-p instead) is equivalent but less efficient. Dropout rate p=0.5 for hidden layers, p=0.1-0.2 for input layers; rarely used after batch normalization.

import numpy as np
def dropout_forward(X, p, training=True, seed=None):
    if not training: return X, None
    rng=np.random.default_rng(seed)
    mask=(rng.random(X.shape)<(1-p))/(1-p)
    return X*mask, mask
def dropout_backward(dout, mask):
    return dout * mask  # mask already includes 1/(1-p) scaling

Two-layer net: z1 = XW1+b1; a1 = relu(z1); z2 = a1W2+b2; loss = cross-entropy(softmax(z2), y). Backward: dz2 = softmax(z2) - one_hot(y); dW2 = a1.T @ dz2 / n; db2 = dz2.mean(0); da1 = dz2 @ W2.T; dz1 = da1 * (z1>0) [ReLU derivative]; dW1 = X.T @ dz1 / n; db1 = dz1.mean(0). This is the chain rule applied layer by layer, propagating gradients from output to input. All O(n·d·h) per batch for input dimension d and hidden dimension h. Store intermediate activations in the forward pass; use them in the backward pass.

import numpy as np
def softmax(z): e=np.exp(z-z.max(1,keepdims=True)); return e/e.sum(1,keepdims=True)
class TwoLayerNN:
    def __init__(self,d,h,c,lr=1e-2):
        self.W1=np.random.randn(d,h)*0.01; self.b1=np.zeros(h)
        self.W2=np.random.randn(h,c)*0.01; self.b2=np.zeros(c); self.lr=lr
    def forward(self,X,y):
        self.X=X; self.z1=X@self.W1+self.b1; self.a1=np.maximum(0,self.z1)
        self.z2=self.a1@self.W2+self.b2; p=softmax(self.z2)
        n=len(y); oh=np.eye(self.W2.shape[1])[y]
        return -np.sum(oh*np.log(p+1e-15))/n, p
    def backward(self,y):
        n=len(y); oh=np.eye(self.W2.shape[1])[y]; p=softmax(self.z2)
        dz2=(p-oh)/n; self.W2-=self.lr*self.a1.T@dz2; self.b2-=self.lr*dz2.sum(0)
        da1=dz2@self.W2.T; dz1=da1*(self.z1>0)
        self.W1-=self.lr*self.X.T@dz1; self.b1-=self.lr*dz1.sum(0)

Scaled dot-product attention computes: Attention(Q,K,V) = softmax(QK^T/√d_k)V, where Q (queries), K (keys), V (values) are matrices of shape (n, d_k) or (n, d_v). Steps: (1) Compute raw scores S = QK^T / √d_k (scale prevents softmax saturation in high dimensions). (2) Apply mask (optional, for causal/padding masking): add -inf to masked positions before softmax. (3) Apply softmax row-wise to get attention weights A. (4) Compute output = A @ V. The output is a weighted sum of values, where weights reflect query-key similarity. Computational complexity O(n²d_k) in sequence length n—the bottleneck for long sequences, motivating linear attention variants.

import numpy as np
def scaled_dot_product_attention(Q, K, V, mask=None):
    d_k = Q.shape[-1]
    scores = Q @ K.swapaxes(-2,-1) / np.sqrt(d_k)  # (...,n,n)
    if mask is not None:
        scores = np.where(mask, scores, -1e9)
    weights = np.exp(scores - scores.max(-1,keepdims=True))
    weights /= weights.sum(-1,keepdims=True)
    return weights @ V, weights

Adam maintains per-parameter first moment (m, exponential MA of gradients) and second moment (v, exponential MA of squared gradients). Initialize m=0, v=0, t=0. At each step: t+=1; m=β₁m+(1-β₁)g; v=β₂v+(1-β₂)g²; m̂=m/(1-β₁^t); v̂=v/(1-β₂^t); θ -= η*m̂/(√v̂+ε). Defaults: β₁=0.9, β₂=0.999, ε=1e-8. The bias correction (1-β^t) is critical in early steps when m and v are near-zero (initialized as zero). Implementing correctly requires separate m, v buffers per parameter. Decoupled weight decay (AdamW): add -η*λ*θ after the adaptive step, not inside the gradient.

import numpy as np
class Adam:
    def __init__(self,params,lr=1e-3,b1=0.9,b2=0.999,eps=1e-8):
        self.params=params; self.lr=lr; self.b1=b1; self.b2=b2; self.eps=eps
        self.m=[np.zeros_like(p) for p in params]
        self.v=[np.zeros_like(p) for p in params]
        self.t=0
    def step(self,grads):
        self.t+=1
        for i,(p,g) in enumerate(zip(self.params,grads)):
            self.m[i]=self.b1*self.m[i]+(1-self.b1)*g
            self.v[i]=self.b2*self.v[i]+(1-self.b2)*g**2
            mh=self.m[i]/(1-self.b1**self.t)
            vh=self.v[i]/(1-self.b2**self.t)
            p-=self.lr*mh/(np.sqrt(vh)+self.eps)

Gradient checking numerically approximates gradients using finite differences: grad_approx[i] = (f(θ+ε·eᵢ) - f(θ-ε·eᵢ)) / (2ε) for each parameter θᵢ. Compare with analytical gradient: relative error = ||grad_analytic - grad_approx|| / (||grad_analytic|| + ||grad_approx||). Error < 1e-7 indicates correct implementation; > 1e-3 indicates a bug. Use ε=1e-5. Caveats: use only during debugging (2p forward passes for p parameters is expensive); disable dropout/batch norm (stochastic/batch-dependent); check a random subset of parameters for large networks. This is the gold standard for validating custom backprop implementations.

import numpy as np
def gradient_check(f, params, analytic_grads, eps=1e-5):
    errs=[]
    for p,ag in zip(params,analytic_grads):
        num_g=np.zeros_like(p)
        it=np.nditer(p,flags=['multi_index'])
        while not it.finished:
            idx=it.multi_index; orig=p[idx]
            p[idx]=orig+eps; fp=f()
            p[idx]=orig-eps; fm=f()
            p[idx]=orig; num_g[idx]=(fp-fm)/(2*eps); it.iternext()
        err=np.linalg.norm(ag-num_g)/(np.linalg.norm(ag)+np.linalg.norm(num_g)+1e-10)
        errs.append(err)
    return errs

LASSO minimizes ||y - Xw||² + λ||w||₁. Coordinate descent cycles over each weight wⱼ: compute the partial residual r_j = y - X_{-j}w_{-j} = y - Xw + X[:,j]*w[j]; then rho_j = X[:,j].T @ r_j / n; update w[j] = soft_threshold(rho_j, λ/2) / (||X[:,j]||²/n). Soft thresholding: sign(rho)·max(|rho|-threshold, 0). Each coordinate update is O(n) and the solution for other coordinates is implicitly updated via the residual. Converges because the LASSO objective is coordinate-wise strictly convex. Cycling coordinate order: cyclic (each feature once per pass) or randomized (sample uniformly). With feature standardization, convergence is faster.

import numpy as np
def soft_thresh(x, t): return np.sign(x)*np.maximum(np.abs(x)-t,0)
def lasso_cd(X,y,lam,max_iter=1000):
    n,d=X.shape; w=np.zeros(d)
    Xnorms=(X**2).sum(0)/n
    for _ in range(max_iter):
        for j in range(d):
            r=y-X@w+X[:,j]*w[j]
            rho=X[:,j]@r/n
            w[j]=soft_thresh(rho,lam/2)/Xnorms[j]
    return w

RNN processes a sequence step by step: h_t = tanh(W_h @ h_{t-1} + W_x @ x_t + b). Forward: initialize h_0=0, compute h_1,...,h_T, compute output o_t = W_y @ h_t + b_y at each step. Store all h_t for BPTT. Backward (full BPTT): backpropagate dL/dh_t through time, accumulating gradients. Gradient w.r.t. W_h: Σ_t dh_t^T @ h_{t-1}. Truncated BPTT limits the backward pass to k steps to avoid O(T) memory and exploding gradients. The derivative through tanh: dh_t/dz_t = (1-h_t²) (element-wise). The vanishing gradient problem arises because gradient is multiplied by W_h^T at each step backward, and (W_h^T)^T → 0 for spectral radius < 1.

import numpy as np
class SimpleRNN:
    def __init__(self,input_dim,hidden_dim):
        self.Wh=np.random.randn(hidden_dim,hidden_dim)*0.01
        self.Wx=np.random.randn(hidden_dim,input_dim)*0.01
        self.b=np.zeros(hidden_dim)
    def forward(self,xs):
        h=np.zeros(self.Wh.shape[0])
        self.hs=[h]; self.xs=xs
        for x in xs:
            h=np.tanh(self.Wh@h+self.Wx@x+self.b)
            self.hs.append(h)
        return self.hs[1:]

Gini impurity G(node) = 1 - Σpᵢ² (where pᵢ is the fraction of class i samples). Range [0, 0.5] for binary; [0, 1-1/K] for K classes. Gini split: choose feature/threshold maximizing ΔG = G(parent) - (n_L/n)G(left) - (n_R/n)G(right). Entropy H = -Σpᵢ log pᵢ; range [0, log K]. Computationally: Gini avoids the log computation, making it slightly faster. Quality: both produce similar trees in practice. Key difference at pure node: Gini=0, entropy=0 (both identify purity correctly). CART (sklearn DecisionTreeClassifier) uses Gini by default. For regression trees, use MSE (mean squared error) as the impurity measure.

import numpy as np
def gini(y):
    n=len(y)
    if n==0: return 0.0
    ps=np.bincount(y)/n
    return 1-(ps**2).sum()
def gini_gain(parent,left,right):
    n=len(parent)
    return gini(parent)-(len(left)/n*gini(left)+len(right)/n*gini(right))

Gradient Boosting builds an additive model F_M(x) = Σf_m(x) by fitting each new tree f_m to the negative gradient (pseudo-residuals) of the loss w.r.t. the current prediction. For MSE loss: pseudo-residuals = y - F_{m-1}(x) (just the residuals). Steps: (1) Initialize F_0 = mean(y). (2) For m=1 to M: compute r = y - F_{m-1}(X); fit a shallow decision tree to (X, r); compute optimal leaf values (mean of residuals in each leaf for MSE); add η * tree to F. Learning rate η controls regularization. Each tree corrects remaining errors. For classification, pseudo-residuals come from log-loss gradient: r = y - sigmoid(F).

import numpy as np
from sklearn.tree import DecisionTreeRegressor
class SimpleGBM:
    def __init__(self,n_trees=100,lr=0.1,max_depth=3):
        self.n_trees=n_trees; self.lr=lr; self.max_depth=max_depth
        self.trees=[]; self.F0=None
    def fit(self,X,y):
        self.F0=y.mean(); F=np.full(len(y),self.F0)
        for _ in range(self.n_trees):
            r=y-F
            t=DecisionTreeRegressor(max_depth=self.max_depth)
            t.fit(X,r); self.trees.append(t)
            F+=self.lr*t.predict(X)
    def predict(self,X):
        return self.F0+self.lr*sum(t.predict(X) for t in self.trees)

SVM primal: minimize (1/2)||w||² + C·Σmax(0, 1-yᵢ(w·xᵢ+b)). The hinge loss max(0,1-yᵢzᵢ) is convex but non-smooth. Subgradient: ∂L/∂w = w + C·Σ(-yᵢxᵢ·1[yᵢzᵢ<1]); ∂L/∂b = C·Σ(-yᵢ·1[yᵢzᵢ<1]). Update with step η_t = 1/(λt) (decaying). Misclassified points (margin < 1) contribute to gradient; correctly classified points with margin > 1 do not. This is the Pegasos algorithm. Each step is O(d) for sparse updates on support vectors. Kernel SVM is solved more efficiently via the dual problem with quadratic programming.

import numpy as np
class SVMPrimal:
    def __init__(self,C=1.0,lr=0.01,n_iter=1000):
        self.C=C; self.lr=lr; self.n_iter=n_iter
    def fit(self,X,y):
        n,d=X.shape; self.w=np.zeros(d); self.b=0.
        for t in range(1,self.n_iter+1):
            eta=self.lr/t
            margins=y*(X@self.w+self.b)
            mask=margins<1
            dw=self.w-self.C*(y[mask,None]*X[mask]).sum(0)
            db=-self.C*y[mask].sum()
            self.w-=eta*dw; self.b-=eta*db
    def predict(self,X): return np.sign(X@self.w+self.b)

GMM assumes data from K Gaussian components with parameters {πₖ, μₖ, Σₖ}. E-step: compute responsibilities rᵢₖ = πₖN(xᵢ;μₖ,Σₖ)/Σⱼπⱼ N(xᵢ;μⱼ,Σⱼ) (soft cluster assignments). M-step: Nₖ=Σrᵢₖ; μₖ=Σrᵢₖxᵢ/Nₖ; Σₖ=Σrᵢₖ(xᵢ-μₖ)(xᵢ-μₖ)^T/Nₖ; πₖ=Nₖ/n. Iterate until log-likelihood converges. Log-likelihood: Σlog(Σₖπₖ N(xᵢ;μₖ,Σₖ)) should increase monotonically (EM guarantees this). Compute in log space for stability. Singular covariance matrices (when Nₖ is small) are regularized by adding ε·I.

import numpy as np
from scipy.stats import multivariate_normal
class GMM:
    def __init__(self,K,n_iter=100):
        self.K=K; self.n_iter=n_iter
    def fit(self,X):
        n,d=X.shape; K=self.K
        idx=np.random.choice(n,K,replace=False)
        self.mu=X[idx].copy(); self.Sigma=[np.eye(d)]*K
        self.pi=np.ones(K)/K
        for _ in range(self.n_iter):
            R=np.array([self.pi[k]*multivariate_normal.pdf(X,self.mu[k],self.Sigma[k]) for k in range(K)]).T
            R/=R.sum(1,keepdims=True)
            Nk=R.sum(0)
            self.mu=(R.T@X)/Nk[:,None]
            self.Sigma=[((R[:,k:k+1]*(X-self.mu[k])).T@(X-self.mu[k]))/Nk[k]+1e-6*np.eye(d) for k in range(K)]
            self.pi=Nk/n

Softmax converts logits z (n×C) to probabilities: pᵢ = exp(zᵢ)/Σexp(zⱼ) (computed stably by subtracting max). Cross-entropy loss: L = -mean(Σ yᵢⱼ log pᵢⱼ) where y is one-hot. Combined gradient (softmax + CE): ∂L/∂zᵢ = (pᵢ - yᵢ)/n. This elegant result arises because the Jacobian of softmax and the CE gradient cancel beautifully. Numerically: log(pᵢ) = zᵢ - max(z) - log(Σexp(z-max(z))) to prevent overflow. The gradient with respect to the logits is the error between predicted probabilities and true labels, normalized by batch size.

import numpy as np
def softmax_cross_entropy(logits, y_true):
    """logits: (n, C), y_true: (n,) integer labels"""
    n = len(y_true)
    shift = logits - logits.max(1, keepdims=True)
    log_sum_exp = np.log(np.exp(shift).sum(1, keepdims=True))
    log_probs = shift - log_sum_exp
    loss = -log_probs[np.arange(n), y_true].mean()
    probs = np.exp(log_probs)
    dlogits = probs.copy()
    dlogits[np.arange(n), y_true] -= 1
    return loss, dlogits / n

Ridge regression in primal: w = (X^TX + λI)^{-1}X^Ty, prediction X_test @ w. Kernel trick: via Woodbury, solution can be written as α = (K + λI)^{-1}y where K_{ij}=k(xᵢ,xⱼ) (n×n kernel matrix). Prediction: f(x_test) = Σαᵢk(xᵢ,x_test). This enables nonlinear regression by replacing the dot product with any PSD kernel. Solving is O(n³) for the kernel matrix inversion. For RBF kernel k(x,y) = exp(-γ||x-y||²). The dual representation expresses the solution as a weighted sum of kernel evaluations, working in possibly infinite-dimensional feature space without explicitly computing features.

import numpy as np
def rbf_kernel(X1, X2, gamma=1.0):
    sq=((X1[:,None,:]-X2[None,:,:])**2).sum(-1)
    return np.exp(-gamma*sq)
class KernelRidgeRegression:
    def __init__(self,lam=1.0,gamma=1.0):
        self.lam=lam; self.gamma=gamma
    def fit(self,X,y):
        self.X_tr=X; K=rbf_kernel(X,X,self.gamma)
        self.alpha=np.linalg.solve(K+self.lam*np.eye(len(X)),y)
    def predict(self,X):
        K_test=rbf_kernel(X,self.X_tr,self.gamma)
        return K_test@self.alpha

Softmax regression extends binary logistic regression to C classes: P(y=c|x) = exp(w_c·x+b_c)/Σexp(w_k·x+b_k). Parameters: W (C×d) and b (C,). Loss: cross-entropy L = -mean(Σ yᵢc log pᵢc). Gradient: dW = (P-Y)^T @ X / n; db = mean(P-Y, axis=0). P is the n×C softmax output; Y is one-hot. Update W and b with gradient descent. With C=2, equivalent to binary logistic regression (redundant parameter). Regularization: add λ||W||²_F to loss, gradient +2λW. Decision boundary is linear: argmax W@x+b is a linear function of x, partitioning space into convex Voronoi regions.

import numpy as np
class SoftmaxRegression:
    def __init__(self,lr=0.1,iters=500,lam=0.):
        self.lr=lr; self.iters=iters; self.lam=lam
    def fit(self,X,y):
        n,d=X.shape; C=len(np.unique(y))
        self.W=np.zeros((C,d)); self.b=np.zeros(C)
        Y=np.eye(C)[y]
        for _ in range(self.iters):
            Z=X@self.W.T+self.b
            Z-=Z.max(1,keepdims=True)
            P=np.exp(Z); P/=P.sum(1,keepdims=True)
            dW=(P-Y).T@X/n+2*self.lam*self.W
            db=(P-Y).mean(0)
            self.W-=self.lr*dW; self.b-=self.lr*db
    def predict(self,X):
        Z=X@self.W.T+self.b; return Z.argmax(1)

Log-sum-exp(x) = log(Σexp(xᵢ)) is computed stably by subtracting the maximum: LSE(x) = max(x) + log(Σexp(xᵢ-max(x))). This prevents overflow (all shifted terms ≤ 1) and underflow (at least one term is exp(0)=1, so sum ≥ 1). LSE appears ubiquitously: log-likelihood of mixture models (log Σπₖ p(x|θₖ)), computing log-partition function in CRFs and Boltzmann machines, numerically stable softmax (subtract LSE from logits), and forward algorithm in HMMs. The log-sum-exp trick is fundamental: any computation involving log(sum(exp(...))) should use it. In probability, LSE is related to free energy in thermodynamics and the log-normalizer of exponential families.

import numpy as np
def logsumexp(x, axis=None, keepdims=False):
    m = x.max(axis=axis, keepdims=True)
    result = np.log(np.exp(x - m).sum(axis=axis, keepdims=True)) + m
    if not keepdims:
        result = result.squeeze(axis=axis)
    return result

# Stable log-softmax:
def log_softmax(x):
    return x - logsumexp(x, axis=-1, keepdims=True)

An embedding layer maps integer indices to dense vectors via a lookup table (matrix E of shape vocab_size × embedding_dim). Forward: output = E[indices], which selects rows. Backward: gradient flows to only the rows corresponding to input indices—a sparse update. For a batch of indices, accumulate gradients: dE = sparse matrix where dE[indices[i]] += dout[i]. Implementing efficiently: use np.add.at(dE, indices, dout) for accumulation (safe for repeated indices). This is fundamentally a matrix multiplication: one-hot(indices) @ E, but using indexing avoids materializing the one-hot matrix. In practice, embedding gradients are sparse; sparse optimizers or gradient accumulation strategies exploit this.

import numpy as np
class Embedding:
    def __init__(self,vocab_size,embed_dim):
        self.E=np.random.randn(vocab_size,embed_dim)*0.01
        self.dE=None
    def forward(self,indices):
        self.indices=indices
        return self.E[indices]  # (batch, embed_dim)
    def backward(self,dout):
        self.dE=np.zeros_like(self.E)
        np.add.at(self.dE,self.indices,dout)  # sparse accumulation
        return self.dE

DBSCAN finds clusters as dense regions separated by low-density areas. Parameters: ε (neighborhood radius) and minPts (minimum neighbors to be a core point). Algorithm: (1) For each unvisited point, find all ε-neighbors. (2) If neighbors ≥ minPts, the point is a core point—start a new cluster. (3) Expand the cluster: add all ε-neighbors and recursively expand from core points among them. (4) Non-core points on the boundary of a cluster are assigned to it; isolated points are noise (-1). Complexity: O(n²) naively; O(n log n) with spatial indexing (KD-tree). Unlike k-means, DBSCAN finds arbitrary-shaped clusters and automatically detects outliers without specifying k.

import numpy as np
from sklearn.metrics import pairwise_distances
def dbscan(X, eps, min_pts):
    n=len(X); labels=np.full(n,-1); visited=np.zeros(n,bool)
    D=pairwise_distances(X); cluster_id=-1
    def expand(i,cid):
        nb=np.where(D[i]<=eps)[0]
        if len(nb)<min_pts: labels[i]=-1; return False
        labels[nb]=cid
        for j in nb:
            if not visited[j]: visited[j]=True; expand(j,cid)
        return True
    for i in range(n):
        if visited[i]: continue
        visited[i]=True
        if expand(i,cluster_id+1): cluster_id+=1
    return labels

The forward algorithm computes P(observations|HMM parameters) in O(TN²) where T is sequence length and N is number of states. Define α_t(i) = P(o₁,...,o_t, s_t=i|λ). Initialization: α_1(i) = π_i · B(i, o₁). Recursion: α_t(i) = [Σⱼ α_{t-1}(j) A_{ji}] · B(i, o_t). Final: P(O|λ) = Σᵢ α_T(i). Compute in log space to prevent underflow: use log-sum-exp for the summation. Scaling: normalize α at each step by the sum (store scale factors for log-likelihood computation). Compared to naive O(N^T) enumeration of all state sequences, the forward algorithm uses DP to reuse computations.

import numpy as np
def hmm_forward(obs, pi, A, B):
    """pi: initial probs (N,), A: transition (N,N), B: emission (N,V)"""
    T=len(obs); N=len(pi)
    log_alpha=np.zeros((T,N))
    log_alpha[0]=np.log(pi+1e-300)+np.log(B[:,obs[0]]+1e-300)
    for t in range(1,T):
        for j in range(N):
            log_alpha[t,j]=np.logaddexp.reduce(log_alpha[t-1]+np.log(A[:,j]+1e-300))
            log_alpha[t,j]+=np.log(B[j,obs[t]]+1e-300)
    return np.logaddexp.reduce(log_alpha[-1])

t-SNE (van der Maaten 2008) maps high-dimensional points to 2D preserving local structure. The high-dimensional affinities use Gaussian kernels: p_{j|i} = exp(-||xᵢ-xⱼ||²/(2σᵢ²)) / Σₖ≠ᵢ exp(-||xᵢ-xₖ||²/(2σᵢ²)). Symmetrize: p_{ij} = (p_{j|i}+p_{i|j})/(2n). Each σᵢ is chosen by binary search to achieve a target perplexity 2^H(pᵢ) where H is entropy. Low-dimensional affinities use Student-t distribution (heavier tails prevent crowding): q_{ij} ∝ (1+||yᵢ-yⱼ||²)^{-1}. Minimize KL(P||Q) via gradient descent. The t-distribution in low-d allows moderately dissimilar points to be pushed far apart, separating clusters.

import numpy as np
def compute_pij(X, perplexity=30.0):
    n=len(X)
    dists=((X[:,None]-X[None,:])**2).sum(-1)
    P=np.zeros((n,n))
    for i in range(n):
        lo,hi=1e-30,1e30; sigma=1.0
        for _ in range(50):
            p=np.exp(-dists[i]/(2*sigma**2)); p[i]=0
            p/=p.sum()+1e-15
            H=-np.sum(p[p>0]*np.log2(p[p>0]))
            if abs(2**H-perplexity)<1e-5: break
            if 2**H>perplexity: hi=sigma
            else: lo=sigma
            sigma=(lo+hi)/2
        P[i]=p
    return (P+P.T)/(2*n)

Log-likelihood of x ~ N(μ, Σ): log p(x) = -d/2·log(2π) - 1/2·log|Σ| - 1/2·(x-μ)^TΣ^{-1}(x-μ). Stable implementation via Cholesky: L = chol(Σ); log|Σ| = 2·sum(log(diag(L))); solve Lv = (x-μ) (forward substitution); Mahalanobis = v·v. For n data points: solve L @ V.T = (X-μ).T; mahal = (V**2).sum(1). Total: -d/2·log(2π) - sum(log(diag(L))) - 0.5·mahal. Using Cholesky avoids computing Σ^{-1} explicitly (numerically unstable for near-singular Σ) and log|Σ| via determinant (numerically unstable for high-d). Add ε·I to Σ for regularization.

import numpy as np
def mvn_log_prob(X, mu, Sigma, eps=1e-6):
    d=X.shape[1]
    Sigma_reg=Sigma+eps*np.eye(d)
    L=np.linalg.cholesky(Sigma_reg)
    log_det=2*np.sum(np.log(np.diag(L)))
    diff=X-mu  # (n,d)
    V=np.linalg.solve(L,diff.T).T  # (n,d)
    mahal=(V**2).sum(1)
    return -0.5*(d*np.log(2*np.pi)+log_det+mahal)

Viterbi finds the maximum likelihood state sequence in O(TN²). Define δ_t(i) = max probability of any path ending in state i at time t. Initialization: δ_1(i) = π_i·B(i,o₁); Ψ_1(i)=0. Recursion: δ_t(j) = max_i[δ_{t-1}(i)·A_{ij}]·B(j,o_t); Ψ_t(j) = argmax_i[δ_{t-1}(i)·A_{ij}]. Backtrace: s*_T = argmax δ_T; s*_t = Ψ_{t+1}(s*_{t+1}). Work in log space to prevent underflow: replace product with sum of logs, and argmax remains valid. Unlike the forward algorithm (sum over paths), Viterbi uses max over paths. Used in speech recognition, POS tagging, and protein sequence labeling.

import numpy as np
def viterbi(obs, pi, A, B):
    T=len(obs); N=len(pi)
    log_delta=np.log(pi+1e-300)+np.log(B[:,obs[0]]+1e-300)
    psi=np.zeros((T,N),int)
    for t in range(1,T):
        trans=log_delta[:,None]+np.log(A+1e-300)
        psi[t]=trans.argmax(0)
        log_delta=trans.max(0)+np.log(B[:,obs[t]]+1e-300)
    path=np.zeros(T,int); path[-1]=log_delta.argmax()
    for t in range(T-2,-1,-1): path[t]=psi[t+1,path[t+1]]
    return path

REINFORCE estimates the policy gradient ∇_θJ(θ) = E[G_t·∇_θ log π_θ(a|s)] using Monte Carlo rollouts. For a k-armed bandit (no state): sample action a from π_θ(a), receive reward r, compute log probability of action, update θ += α·r·∇_θ log π_θ(a). Using softmax policy: π(a) = softmax(θ)[a]. Gradient: ∂log π(a)/∂θ = 1[i=a] - π(a) for each logit θᵢ. This is the score function (REINFORCE) gradient. Baseline b subtracted from reward reduces variance without biasing the gradient: θ += α·(r-b)·∇log π(a). Common baseline: running average of rewards. This is the foundation of PPO and actor-critic methods.

import numpy as np
class REINFORCEBandit:
    def __init__(self,k,lr=0.1): self.theta=np.zeros(k); self.lr=lr; self.baseline=0.
    def policy(self):
        e=np.exp(self.theta-self.theta.max())
        return e/e.sum()
    def update(self,action,reward,t):
        self.baseline+=(reward-self.baseline)/t
        pi=self.policy()
        grad=np.zeros_like(self.theta); grad[action]=1-pi[action]
        for i in range(len(self.theta)):
            if i!=action: grad[i]=-pi[i]
        self.theta+=self.lr*(reward-self.baseline)*grad
    def select_action(self):
        return np.random.choice(len(self.theta), p=self.policy())

Multi-head attention projects Q, K, V into h parallel subspaces, computes attention in each, and concatenates: MultiHead(Q,K,V) = Concat(head₁,...,headₕ)W_O, where headᵢ = Attention(QW_Qi, KW_Ki, VW_Vi). Each head has dimension d_k = d_model/h. Implementation: stack all projections as a single tensor operation by reshaping: reshape Q to (batch, seq, h, d_k) then transpose to (batch, h, seq, d_k). Compute attention scores for all heads simultaneously: (batch, h, seq, seq). This is memory-efficient—all heads are computed in parallel. Each head can attend to different aspects of the input (positional, syntactic, semantic), and the concatenated output captures all aspects. Causal masking adds -inf to upper-triangular positions.

import numpy as np
class MultiHeadAttention:
    def __init__(self,d_model,h):
        self.h=h; self.dk=d_model//h
        scale=np.sqrt(d_model)
        self.WQ=np.random.randn(d_model,d_model)/scale
        self.WK=np.random.randn(d_model,d_model)/scale
        self.WV=np.random.randn(d_model,d_model)/scale
        self.WO=np.random.randn(d_model,d_model)/scale
    def forward(self,Q,K,V,mask=None):
        B,S,D=Q.shape; h=self.h; dk=self.dk
        def proj_split(X,W):
            return (X@W).reshape(B,S,h,dk).transpose(0,2,1,3)
        Qp,Kp,Vp=proj_split(Q,self.WQ),proj_split(K,self.WK),proj_split(V,self.WV)
        sc=Qp@Kp.swapaxes(-2,-1)/np.sqrt(dk)
        if mask is not None: sc=np.where(mask,-1e9,sc)
        A=np.exp(sc-sc.max(-1,keepdims=True)); A/=A.sum(-1,keepdims=True)
        out=(A@Vp).transpose(0,2,1,3).reshape(B,S,D)
        return out@self.WO

Ridge regression minimizes ||y - Xw||² + λ||w||². Closed-form: w* = (X^TX + λI)^{-1}X^Ty, O(d³) for the solve. Gradient descent: ∇L = -2X^T(y-Xw) + 2λw; update w -= η·∇L. Convergence: O(κ/ε) iterations where κ = condition number of (X^TX/n + λI). Regularization λ reduces condition number (adds λ to all eigenvalues), speeding GD. Closed-form is exact in one step but O(d³); GD is iterative but cheap per step O(nd). For d > 10⁴, GD is preferred. The optimal learning rate is 1/(max eigenvalue of (X^TX + λI)/n) = 1/(||X||²/n + λ). Stochastic gradient descent further reduces per-step cost to O(d) per sample.

import numpy as np
def ridge_gd(X,y,lam,lr=1e-3,n_iter=1000):
    n,d=X.shape; w=np.zeros(d)
    for _ in range(n_iter):
        grad=-2*X.T@(y-X@w)/n+2*lam*w
        w-=lr*grad
    return w
def ridge_closed(X,y,lam):
    d=X.shape[1]
    return np.linalg.solve(X.T@X+lam*np.eye(d), X.T@y)

K-means++ chooses centroids sequentially: (1) Select first centroid uniformly at random. (2) For each subsequent centroid, select point x with probability proportional to D(x)² (squared distance to nearest already-chosen centroid). (3) Repeat until k centroids selected. This probabilistic seeding spreads centroids across the data, avoiding the degenerate case where random initialization places multiple centroids in the same dense cluster while leaving others uncovered. Theoretical guarantee: expected cost of k-means++ initialization is O(log k) times the optimal k-means cost (vs arbitrarily bad random init). Runtime of initialization: O(nkd). sklearn.KMeans uses k-means++ by default (init='k-means++').

import numpy as np
def kmeans_plus_plus(X, k, seed=42):
    rng=np.random.default_rng(seed)
    centroids=[X[rng.integers(len(X))]]
    for _ in range(1,k):
        dists=np.array([min(np.linalg.norm(x-c)**2 for c in centroids) for x in X])
        probs=dists/dists.sum()
        centroids.append(X[rng.choice(len(X),p=probs)])
    return np.array(centroids)

An autoencoder compresses input to a lower-dimensional latent code and reconstructs it. Architecture: Encoder: z = relu(X @ W_e + b_e) (input_dim → latent_dim). Decoder: X_hat = sigmoid(z @ W_d + b_d) (latent_dim → input_dim, sigmoid for [0,1] pixel values). Loss: MSE ||X - X_hat||² or binary cross-entropy -[X·log(X_hat) + (1-X)·log(1-X_hat)]. Backpropagation proceeds through decoder then encoder. Applications: dimensionality reduction (PCA alternative, nonlinear), denoising (train on noisy inputs, reconstruct clean), anomaly detection (high reconstruction error = anomaly), and feature learning for downstream tasks. Variational autoencoders add a KL regularization term to enforce a structured latent space.

import numpy as np
class Autoencoder:
    def __init__(self,input_dim,latent_dim,lr=1e-3):
        scale=0.01; self.lr=lr
        self.We=np.random.randn(input_dim,latent_dim)*scale; self.be=np.zeros(latent_dim)
        self.Wd=np.random.randn(latent_dim,input_dim)*scale; self.bd=np.zeros(input_dim)
    def encode(self,X): return np.maximum(0,X@self.We+self.be)
    def decode(self,z):
        s=z@self.Wd+self.bd; return 1/(1+np.exp(-s))
    def loss(self,X):
        z=self.encode(X); Xh=self.decode(z)
        return 0.5*((X-Xh)**2).mean(), Xh, z

Numerical gradient checking verifies analytical gradients by computing central finite differences: (f(x+ε)-f(x-ε))/(2ε) for each parameter. Compare with analytical gradient using relative error = ||grad_num - grad_anal|| / (||grad_num|| + ||grad_anal|| + 1e-10). Relative error < 1e-7 is excellent; > 1e-3 indicates a bug. Use ε=1e-5. The central difference is O(ε²) accurate, much better than forward difference O(ε). Test on a small batch with a fixed seed; disable stochastic operations. Check all parameters: weights, biases, embeddings. This technique is essential for implementing custom backward passes and verifying custom CUDA kernels.

import numpy as np
def check_grad(loss_fn, params, analytic_grads, eps=1e-5):
    for i, (p, ag) in enumerate(zip(params, analytic_grads)):
        ng = np.zeros_like(p)
        for idx in np.ndindex(*p.shape):
            p[idx] += eps; fp = loss_fn()
            p[idx] -= 2*eps; fm = loss_fn()
            p[idx] += eps
            ng[idx] = (fp - fm) / (2*eps)
        err = np.linalg.norm(ag-ng)/(np.linalg.norm(ag)+np.linalg.norm(ng)+1e-10)
        print(f'Param {i}: rel_err={err:.2e}', 'OK' if err<1e-5 else 'BUG')

ReLU (Rectified Linear Unit) f(x) = max(0, x) sets negative values to zero. Derivative: f'(x) = 1 if x > 0, 0 if x < 0, undefined at x=0 (set to 0 by convention). In NumPy: np.maximum(0, x) for forward, (x > 0).astype(float) for derivative. Leaky ReLU: f(x) = x if x>0 else αx (typically α=0.01); gradient = 1 if x>0 else α. Avoids dead neurons (where all inputs negative, gradient permanently zero). ELU: f(x) = x if x>0 else α(e^x-1); smooth at zero, negative saturation. GELU: f(x) = x·Φ(x) where Φ is standard normal CDF; smoother than ReLU, used in transformers. All computed element-wise with no parameters (except Leaky ReLU's α).

import numpy as np
def relu(x): return np.maximum(0, x)
def relu_grad(x): return (x > 0).astype(float)
def leaky_relu(x, alpha=0.01): return np.where(x>0, x, alpha*x)
def leaky_relu_grad(x, alpha=0.01): return np.where(x>0, 1., alpha)
def elu(x, alpha=1.0): return np.where(x>0, x, alpha*(np.exp(x)-1))

L2 regularization (weight decay) adds λ||w||² to the loss: gradient becomes ∇L + 2λw; update: w -= lr*(grad + 2λw). Equivalently, w *= (1-2λ·lr) then w -= lr*grad. L1 regularization adds λ||w||₁: gradient becomes ∇L + λ·sign(w); update: w -= lr*(grad + λ·sign(w)). L1 promotes sparsity (zeros out small weights via soft-thresholding effect); L2 shrinks all weights proportionally without zeroing. Combined elastic net: λ₁||w||₁ + λ₂||w||². In PyTorch, L2 is implemented via weight_decay in the optimizer (directly decays weights each step); L1 requires manual implementation. For neural networks, L2 is standard; L1 is used in sparse NLP applications.

import numpy as np
def update_l2(w, grad, lr, lam):
    return w - lr * (grad + 2*lam*w)  # L2
def update_l1(w, grad, lr, lam):
    return w - lr * (grad + lam*np.sign(w))  # L1
def update_elastic(w, grad, lr, l1, l2):
    return w - lr*(grad + l1*np.sign(w) + 2*l2*w)

Multiclass SVM (Crammer-Singer) hinge loss for scores s=Xw: for each sample i, loss = Σⱼ≠y max(0, sⱼ - sy + 1) (sum of margin violations for all incorrect classes). Gradient w.r.t. scores: for incorrect class j, if sⱼ-sy+1>0, d_sⱼ=+1 and d_sy-=1; else 0. Gradient w.r.t. W: dW_j += x for violating j; dW_y -= (number of violations)·x for true class y. Vectorized implementation: compute scores S=XW, margins = max(0, S - S[y]+1), set S[y]=0 per sample, sum margins for loss, count violations per sample for gradient. Total O(nkd) per batch. More interpretable than cross-entropy: encourages a specific margin between correct and incorrect classes.

import numpy as np
def svm_loss(X, y, W, delta=1.0):
    n = len(y)
    S = X @ W  # (n, C)
    correct = S[np.arange(n), y][:, None]  # (n,1)
    margins = np.maximum(0, S - correct + delta)  # (n,C)
    margins[np.arange(n), y] = 0
    loss = margins.sum() / n
    mask = (margins > 0).astype(float)
    mask[np.arange(n), y] = -mask.sum(axis=1)
    dW = X.T @ mask / n
    return loss, dW

Mini-batch SGD with momentum maintains a velocity buffer for each parameter. At each training step: (1) Sample a random mini-batch of size B indices. (2) Compute forward pass and loss on the mini-batch. (3) Backpropagate to get gradients. (4) Update velocity: v = β·v + g (or v = β·v + (1-β)·g for normalized momentum). (5) Update parameters: w -= lr·v. Momentum smooths the gradient estimates across batches, reducing oscillation in narrow valleys. With β=0.9, the effective learning rate is amplified by 1/(1-β)=10 in consistent gradient directions. Data shuffling before each epoch ensures each mini-batch is representative and prevents cyclic patterns in gradient estimates.

import numpy as np
def train_sgd_momentum(model, X, y, lr=1e-2, beta=0.9, batch_size=32, epochs=10):
    n = len(X)
    velocity = [np.zeros_like(p) for p in model.params]
    for ep in range(epochs):
        idx = np.random.permutation(n)
        for start in range(0, n, batch_size):
            batch = idx[start:start+batch_size]
            loss, grads = model.forward_backward(X[batch], y[batch])
            for i, (p, g) in enumerate(zip(model.params, grads)):
                velocity[i] = beta*velocity[i] + g
                p -= lr * velocity[i]

The Perceptron is the simplest binary classifier: predict sign(w·x+b). Update rule: if yᵢ(w·xᵢ+b) ≤ 0 (misclassified), update w += lr·yᵢxᵢ; b += lr·yᵢ. For linearly separable data, the Perceptron convergence theorem guarantees convergence in at most R²/γ² updates, where R = max||xᵢ|| and γ is the margin. For non-separable data, it cycles forever—early stopping or a fixed number of epochs is needed. The Perceptron with kernel trick (replacing x·x' with k(x,x')) becomes the kernel Perceptron, capable of learning nonlinear boundaries. It is the building block of neural networks: a single neuron is a Perceptron with a smooth activation function.

import numpy as np
class Perceptron:
    def __init__(self, lr=1.0, n_iter=1000):
        self.lr=lr; self.n_iter=n_iter
    def fit(self, X, y):
        n,d=X.shape; self.w=np.zeros(d); self.b=0.
        for _ in range(self.n_iter):
            for xi,yi in zip(X,y):
                if yi*(self.w@xi+self.b)<=0:
                    self.w+=self.lr*yi*xi
                    self.b+=self.lr*yi
    def predict(self,X): return np.sign(X@self.w+self.b)

PCA with whitening produces decorrelated, unit-variance features. Steps: (1) Center X. (2) Compute SVD of X: U, s, Vt. (3) Whitening matrix W = Vt.T @ diag(1/s * sqrt(n-1)). (4) Whitened data Z = X @ W has identity covariance. Verify: (1/(n-1)) Z.T @ Z ≈ I. ZCA whitening (zero-phase component analysis): W_ZCA = Vt.T @ diag(1/s) @ Vt, applied as Z = X @ W_ZCA—similar covariance but preserves the original feature directions as much as possible. ZCA whitening is used in image preprocessing (natural image statistics) and as a preprocessing step before PCA. Numerical stability: add ε to s before inverting to avoid division by near-zero singular values.

import numpy as np
def pca_whiten(X, eps=1e-6):
    n = len(X)
    Xc = X - X.mean(0)
    U, s, Vt = np.linalg.svd(Xc, full_matrices=False)
    W = Vt.T / (s / np.sqrt(n-1) + eps)  # whitening matrix
    Z = Xc @ W
    print('Cov check:', np.allclose((Z.T@Z)/(n-1), np.eye(W.shape[1]), atol=1e-5))
    return Z, W

A 1D conv layer applies F filters of length K to input X (batch×length×channels_in) producing output (batch×(L-K+1)×F). Naive implementation: slide each filter over the input, computing dot products. Efficient im2col approach: reshape input windows into columns, then matrix multiply with filter weights. For padding='same', add (K-1)//2 zeros on each side. Stride>1 reduces output length. The output at position j with filter f: y[b,j,f] = sum over k,c of x[b,j+k,c]*w[k,c,f] + bias[f]. In ML: used for sequence models, audio processing, and as an efficient approximation to time-series feature extraction. Dilation (atrous convolution) expands the receptive field without increasing parameters.

import numpy as np
from numpy.lib.stride_tricks import sliding_window_view
def conv1d_forward(X, W, b, stride=1, pad=0):
    """X:(B,L,Cin), W:(K,Cin,Cout), b:(Cout,)"""
    B,L,Cin=X.shape; K,_,Cout=W.shape
    if pad>0: X=np.pad(X,((0,0),(pad,pad),(0,0)))
    Lout=(X.shape[1]-K)//stride+1
    # im2col
    windows=sliding_window_view(X,window_shape=(K,Cin),axis=(1,2))
    windows=windows[:,::stride,0,:,:]  # (B,Lout,K,Cin)
    out=windows.reshape(B,Lout,-1)@W.reshape(-1,Cout)+b
    return out

Layer normalization (Ba et al. 2016) normalizes across the feature dimension (not batch dimension), making it batch-size independent and suitable for RNNs and transformers. For input x (batch×features), compute mean and variance along the feature axis: μ = x.mean(-1, keepdims=True); σ² = x.var(-1, keepdims=True). Normalize: x_hat = (x-μ)/√(σ²+ε). Scale and shift: y = γ·x_hat + β (learnable per-feature γ, β). Backward: similar to batch norm but reduced over the last axis. LayerNorm vs BatchNorm: LayerNorm is applied per sample (no statistics over batch), making inference identical to training and eliminating the need for running statistics. Standard in transformers, replacing BatchNorm.

import numpy as np
def layer_norm(x, gamma, beta, eps=1e-5):
    mu = x.mean(axis=-1, keepdims=True)
    var = x.var(axis=-1, keepdims=True)
    x_hat = (x - mu) / np.sqrt(var + eps)
    return gamma * x_hat + beta

X = np.random.randn(32, 512)  # batch, features
gamma = np.ones(512); beta = np.zeros(512)
out = layer_norm(X, gamma, beta)
print(out.mean(-1).max())  # ~0
print(out.std(-1).min())   # ~1

A transformer encoder block consists of: (1) Multi-head self-attention with residual connection and layer norm: x = LayerNorm(x + MultiHeadAttn(x, x, x)). (2) Feed-forward network with residual: x = LayerNorm(x + FFN(x)), where FFN(x) = W₂·ReLU(W₁x+b₁)+b₂ (with expansion factor 4). Both residual connections use pre-norm (LayerNorm before sublayer) or post-norm (after). Pre-norm is more stable for deep transformers. The attention mask handles padding tokens. Dropout is applied after each sublayer before the residual addition. Total parameters: d²·(4+4+8) for attention projections and FFN (approximately 12d² per layer for expansion=4).

import numpy as np
class TransformerEncoderBlock:
    def __init__(self,d,h,dff,eps=1e-5):
        self.attn=MultiHeadAttention(d,h) if False else None  # placeholder
        self.W1=np.random.randn(d,dff)*0.02; self.b1=np.zeros(dff)
        self.W2=np.random.randn(dff,d)*0.02; self.b2=np.zeros(d)
        self.gamma1=np.ones(d); self.beta1=np.zeros(d)
        self.gamma2=np.ones(d); self.beta2=np.zeros(d)
    def ffn(self,x):
        return np.maximum(0,x@self.W1+self.b1)@self.W2+self.b2
    def layer_norm(self,x,g,b,eps=1e-5):
        mu=x.mean(-1,keepdims=True); v=x.var(-1,keepdims=True)
        return g*((x-mu)/np.sqrt(v+eps))+b

Triplet loss trains embeddings so similar (anchor-positive) pairs are closer than dissimilar (anchor-negative) pairs by at least a margin α. Loss = mean(max(0, d(a,p) - d(a,n) + α)) over triplets. With Euclidean distance d(x,y) = ||x-y||₂. Hard negative mining selects the hardest (closest) negatives to accelerate training. Semi-hard mining: select negatives farther than the positive but within the margin. Online triplet mining generates triplets within the batch for efficiency. Compared to contrastive loss (pairwise), triplet loss considers relative distances (ternary constraint), which is more informative. Used in face recognition (FaceNet), image retrieval, and few-shot learning.

import numpy as np
def triplet_loss(anchor, positive, negative, margin=0.2):
    d_ap = np.linalg.norm(anchor-positive, axis=1)  # (n,)
    d_an = np.linalg.norm(anchor-negative, axis=1)  # (n,)
    loss = np.maximum(0, d_ap - d_an + margin)
    return loss.mean()

a=np.random.randn(16,128); p=a+0.1*np.random.randn(16,128)
n=np.random.randn(16,128)
print(triplet_loss(a,p,n))

InfoNCE (used in SimCLR, MoCo) maximizes agreement between augmented views of the same image. For a batch of N pairs (zᵢ, z'ᵢ), treat each pair as a positive example and 2(N-1) others as negatives. NT-Xent loss: L = -mean(log(exp(sim(zᵢ,z'ᵢ)/τ) / Σⱼ≠ᵢ exp(sim(zᵢ,zⱼ)/τ))) where sim is cosine similarity and τ is temperature. Implementation: normalize all embeddings, compute (2N)×(2N) similarity matrix, mask the diagonal (self-similarity), and compute cross-entropy with positive pairs as targets. Temperature τ controls the sharpness: lower τ = harder negatives, larger loss gradients. Temperature τ=0.07 is typical for SimCLR.

import numpy as np
def nt_xent_loss(z1, z2, temperature=0.07):
    n = len(z1)
    z1 = z1 / np.linalg.norm(z1, axis=1, keepdims=True)
    z2 = z2 / np.linalg.norm(z2, axis=1, keepdims=True)
    Z = np.concatenate([z1, z2], axis=0)  # (2n, d)
    sim = Z @ Z.T / temperature  # (2n, 2n)
    np.fill_diagonal(sim, -1e9)  # mask diagonal
    labels = np.concatenate([np.arange(n,2*n), np.arange(n)])
    logp = sim - np.log(np.exp(sim).sum(1, keepdims=True))
    return -logp[np.arange(2*n), labels].mean()

Focal Loss (Lin et al. 2017, RetinaNet) addresses class imbalance by down-weighting easy negatives. Standard CE: -log(p_t). Focal Loss: -(1-p_t)^γ · log(p_t), where p_t = p for y=1 else 1-p. The modulating factor (1-p_t)^γ reduces loss for well-classified examples (high p_t, factor ≈ 0) and focuses training on hard misclassified examples. γ=2 is typical. Combined with α-balancing: -α_t(1-p_t)^γ log(p_t) where α_t=α for foreground and 1-α for background. The key insight: in object detection, background boxes vastly outnumber foreground, and standard CE leads to gradient domination by easy negatives. Focal loss makes the model focus on hard examples without explicit hard example mining.

import numpy as np
def focal_loss(y_true, y_pred, gamma=2.0, alpha=0.25, eps=1e-7):
    y_pred = np.clip(y_pred, eps, 1-eps)
    pt = np.where(y_true==1, y_pred, 1-y_pred)
    alpha_t = np.where(y_true==1, alpha, 1-alpha)
    focal = -alpha_t * (1-pt)**gamma * np.log(pt)
    return focal.mean()

y_true = np.array([1,0,1,0,1])
y_pred = np.array([0.9,0.1,0.4,0.8,0.6])
print(focal_loss(y_true, y_pred))

Huber loss combines MSE for small errors and MAE for large errors: L_δ(y,ŷ) = ½(y-ŷ)² if |y-ŷ|≤δ, else δ(|y-ŷ|-δ/2). Gradient: (y-ŷ) if |y-ŷ|≤δ, else δ·sign(y-ŷ). The transition threshold δ controls the boundary. Preferred over MSE: robust to outliers (linear instead of quadratic growth for large errors, preventing exploding gradients). Preferred over MAE: differentiable everywhere (MAE has undefined gradient at 0), enabling second-order methods. Used in regression with heavy-tailed noise, Q-learning in DQN (where reward differences can be large), and Faster R-CNN bounding box regression. δ=1 recovers smooth L1 loss in PyTorch (nn.SmoothL1Loss).

import numpy as np
def huber_loss(y_true, y_pred, delta=1.0):
    err = np.abs(y_true - y_pred)
    loss = np.where(err<=delta, 0.5*err**2, delta*(err-0.5*delta))
    return loss.mean()
def huber_grad(y_true, y_pred, delta=1.0):
    err = y_pred - y_true
    return np.where(np.abs(err)<=delta, err, delta*np.sign(err))

SGDR (Loshchilov & Hutter 2017) cycles the learning rate using a cosine decay: η_t = η_min + ½(η_max-η_min)(1+cos(πt_cur/T_cur)). When t_cur reaches T_cur, restart with t_cur=0 and optionally multiply T_cur by T_mult (e.g., 2) for increasingly long cycles. Warm restarts allow the optimizer to escape local minima by periodically increasing the learning rate. Ensemble by saving checkpoints at the end of each cycle—these diverse solutions can be averaged for better generalization (snapshot ensembles). The schedule requires only knowing the current step and cycle parameters, needing no validation feedback (unlike ReduceLROnPlateau). Commonly T_0=10 epochs, T_mult=2.

import math
class CosineAnnealingWR:
    def __init__(self, lr_max, lr_min=0., T0=10, T_mult=2):
        self.lr_max=lr_max; self.lr_min=lr_min
        self.T0=T0; self.T_mult=T_mult; self.T_cur=0; self.T_i=T0
    def step(self):
        lr=self.lr_min+0.5*(self.lr_max-self.lr_min)*(1+math.cos(math.pi*self.T_cur/self.T_i))
        self.T_cur+=1
        if self.T_cur>=self.T_i:
            self.T_cur=0; self.T_i*=self.T_mult
        return lr

mAP evaluates object detection by averaging precision-recall curves across IoU thresholds and classes. For each class: (1) Sort predictions by confidence. (2) Compute precision P(k) and recall R(k) at each rank k. (3) Compute area under P-R curve (AUC) as the Average Precision (AP) using the 11-point interpolation or all-points interpolation. (4) Average AP across all classes = mAP. For COCO mAP, average over IoU thresholds [0.5:0.05:0.95]. IoU between predicted and ground-truth box: IoU = intersection area / union area. A prediction is a true positive if IoU ≥ threshold (e.g., 0.5 for PASCAL VOC) and the class matches; otherwise a false positive. Each ground truth can be matched to at most one prediction.

import numpy as np
def iou(box1, box2):
    x1=max(box1[0],box2[0]); y1=max(box1[1],box2[1])
    x2=min(box1[2],box2[2]); y2=min(box1[3],box2[3])
    inter=max(0,x2-x1)*max(0,y2-y1)
    a1=(box1[2]-box1[0])*(box1[3]-box1[1])
    a2=(box2[2]-box2[0])*(box2[3]-box2[1])
    return inter/(a1+a2-inter+1e-6)
def compute_ap(precisions, recalls):
    recalls=np.concatenate([[0],recalls,[1]])
    precisions=np.concatenate([[1],precisions,[0]])
    for i in range(len(precisions)-2,-1,-1):
        precisions[i]=max(precisions[i],precisions[i+1])
    return np.sum(np.diff(recalls)*precisions[1:])