Mathematics for AI

Vectors & Vector Spaces

What is a Vector?

A vector is an ordered list of numbers representing a point or direction in space. In AI, every data point is a vector. An image is a vector of pixel values. A sentence is a vector of word embeddings. A customer profile is a vector of features.

A vector in ℝ² has 2 components (2D plane), ℝ³ has 3 (3D space), and ℝⁿ has n components (n-dimensional hyperspace). In ML, n can be thousands or millions — a 224×224 color image is a vector in ℝ¹⁵⁰,⁵²⁸.

Vector Operations

Addition and Scalar Multiplication

Python
import numpy as np

u = np.array([3, 4, 1])
v = np.array([1, -2, 5])

print("Addition:", u + v)          # [4, 2, 6]
print("Scalar mult:", 3 * u)       # [9, 12, 3]
print("Magnitude:", np.linalg.norm(u))  # √(9+16+1) = √26 ≈ 5.1
print("Unit vector:", u / np.linalg.norm(u))  # Normalize to length 1

Dot Product

The dot product measures how aligned two vectors are. It is the single most important operation in ML — used in every neural network layer, every linear model, every attention mechanism.

If the dot product is positive, vectors point in similar directions. If zero, they are perpendicular (orthogonal). If negative, they point in opposite directions.

Python
u = np.array([1, 2, 3])
v = np.array([4, 5, 6])

dot = np.dot(u, v)  # 1*4 + 2*5 + 3*6 = 32

# Cosine similarity — measures angle between vectors
cos_sim = np.dot(u, v) / (np.linalg.norm(u) * np.linalg.norm(v))
print(f"Cosine similarity: {cos_sim:.4f}")  # 0.9746 — very similar

# This is how search engines find similar documents!

Cross Product (3D only)

The cross product produces a vector perpendicular to both input vectors. Used in computer graphics, 3D physics, and robotics.

Python
u = np.array([1, 0, 0])
v = np.array([0, 1, 0])
cross = np.cross(u, v)
print("Cross product:", cross)  # [0, 0, 1] — perpendicular to both

Vector Spaces, Basis & Dimension

A vector space is a collection of vectors that is closed under addition and scalar multiplication. The basis is the smallest set of vectors that can represent every vector in the space through linear combinations. The dimension is the number of basis vectors.

In ℝ³, the standard basis is {[1,0,0], [0,1,0], [0,0,1]}. Any 3D vector can be written as a combination of these three. In ML, finding a good basis (like PCA does) means finding the most informative directions in your data.

Linear Independence

Vectors are linearly independent if no vector can be written as a combination of the others. If vectors are dependent, they carry redundant information. In ML, redundant features waste computation and can cause numerical instability.

Python
# Check linear independence via matrix rank
vectors = np.array([
    [1, 0, 0],
    [0, 1, 0],
    [1, 1, 0]  # This = vec1 + vec2 → dependent!
])
print("Rank:", np.linalg.matrix_rank(vectors))  # 2, not 3 → dependent

Matrices & Matrix Operations

What is a Matrix?

A matrix is a 2D array of numbers with m rows and n columns (m × n). In AI, matrices represent: datasets (rows = samples, columns = features), weight matrices in neural networks, images (pixel grids), and transformation operations.

Python
import numpy as np

# Dataset: 4 samples, 3 features
X = np.array([
    [5.1, 3.5, 1.4],
    [4.9, 3.0, 1.4],
    [7.0, 3.2, 4.7],
    [6.3, 3.3, 6.0]
])
print(f"Shape: {X.shape}")  # (4, 3) — 4×3 matrix

Matrix Multiplication

Matrix multiplication is not element-wise. For A (m×n) and B (n×p), the result C = AB is (m×p) where each element is a dot product of a row of A with a column of B.

Python
# Neural network forward pass: y = Wx + b
W = np.random.randn(64, 128)    # 64 neurons, 128 inputs
x = np.random.randn(128, 1)     # input vector
b = np.random.randn(64, 1)      # bias vector

y = W @ x + b                    # (64×128) @ (128×1) = (64×1)
print(f"Output shape: {y.shape}")  # (64, 1) — 64 neuron outputs

Transpose, Inverse & Determinant

Python
A = np.array([[2, 1], [5, 3]])

# Transpose: flip rows and columns
print("Aᵀ =\n", A.T)

# Determinant: scalar that measures how A scales area/volume
det = np.linalg.det(A)
print(f"det(A) = {det:.1f}")  # 2*3 - 1*5 = 1

# Inverse: A⁻¹ such that A·A⁻¹ = I (only if det ≠ 0)
A_inv = np.linalg.inv(A)
print("A⁻¹ =\n", A_inv)
print("A·A⁻¹ =\n", (A @ A_inv).round(2))  # Identity matrix

Special Matrices

A = np.array([[1,2],[3,4]])
B = np.array([[5,6],[7,8]])
print(A @ B)  # [[19,22],[43,48]]

Eigenvalues, Eigenvectors & Matrix Decompositions

Eigenvalues & Eigenvectors

An eigenvector of matrix A is a special vector that, when transformed by A, only gets scaled (not rotated). The scaling factor is the eigenvalue.

Think of it this way: most vectors change direction when multiplied by A. Eigenvectors are the stubborn ones that refuse to change direction — they only stretch or shrink. The eigenvalue tells you by how much.

Python
import numpy as np

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

eigenvalues, eigenvectors = np.linalg.eig(A)
print("Eigenvalues:", eigenvalues)    # [5, 2]
print("Eigenvectors:\n", eigenvectors)

# Verify: A·v = λ·v
v = eigenvectors[:, 0]   # First eigenvector
lam = eigenvalues[0]     # First eigenvalue
print("A·v =", A @ v)
print("λ·v =", lam * v)  # Same! ✓

Singular Value Decomposition (SVD)

SVD decomposes any matrix (even non-square) into three matrices:

Where U contains left singular vectors, Σ is diagonal with singular values (importance scores), and Vᵀ contains right singular vectors. The singular values are sorted from largest to smallest.

Python
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

U, sigma, Vt = np.linalg.svd(A)
print("Singular values:", sigma.round(3))

# Low-rank approximation: keep only top-k singular values
# This is how image compression and recommendation engines work!
k = 2
A_approx = U[:, :k] @ np.diag(sigma[:k]) @ Vt[:k, :]
print("Rank-2 approximation:\n", A_approx.round(2))

PCA via Eigendecomposition

Python
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris

iris = load_iris()
X = iris.data  # 150 samples, 4 features

# Reduce 4D → 2D using PCA
pca = PCA(n_components=2)
X_2d = pca.fit_transform(X)

print(f"Original: {X.shape}")            # (150, 4)
print(f"Reduced: {X_2d.shape}")            # (150, 2)
print(f"Variance retained: {pca.explained_variance_ratio_.sum():.2%}")

# Under the hood, PCA computes eigenvectors of the covariance matrix
cov_matrix = np.cov(X.T)
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
print("Top 2 eigenvalues:", sorted(eigenvalues, reverse=True)[:2])

Tensors (Multi-dimensional Arrays)

The Tensor Hierarchy

Python
import numpy as np

# Batch of 32 color images, 224×224 pixels, 3 channels (RGB)
batch = np.random.randn(32, 3, 224, 224)
print(f"Shape: {batch.shape}")  # (32, 3, 224, 224)
print(f"Rank: {batch.ndim}")    # 4
print(f"Total elements: {batch.size:,}")  # 4,816,896

# Reshaping — critical for connecting conv layers to dense layers
flat = batch.reshape(32, -1)   # (32, 150528)
print(f"Flattened: {flat.shape}")

# Broadcasting — NumPy auto-expands dimensions
image = np.random.randn(3, 224, 224)
mean = np.array([0.485, 0.456, 0.406]).reshape(3, 1, 1)
normalized = image - mean  # mean broadcasts across H×W

Einstein Summation (einsum)

Einsum is a powerful notation that expresses any tensor operation in a compact string:

Python
A = np.random.randn(3, 4)
B = np.random.randn(4, 5)

# Matrix multiplication
C = np.einsum('ik,kj->ij', A, B)  # Same as A @ B

# Batch matrix multiplication (used in attention)
Q = np.random.randn(8, 10, 64)  # 8 heads, 10 tokens, 64 dim
K = np.random.randn(8, 10, 64)
attn = np.einsum('bij,bkj->bik', Q, K)  # (8, 10, 10) attention scores

Derivatives & Partial Derivatives

The Derivative: Rate of Change

The derivative f'(x) tells you how fast f(x) changes as x changes. In ML, the loss function L(w) measures error, and ∂L/∂w tells you how the error changes as you adjust each weight — this is how neural networks learn.

Essential Differentiation Rules

Partial Derivatives

When a function has multiple variables, a partial derivative measures the rate of change with respect to one variable while holding all others constant. In a neural network with millions of weights, we compute partial derivatives with respect to each weight separately.

Python
import numpy as np

# Numerical derivative approximation
def numerical_derivative(f, x, h=1e-7):
    return (f(x + h) - f(x - h)) / (2 * h)

# Example: f(x) = x³ + 2x, f'(x) = 3x² + 2
f = lambda x: x**3 + 2*x

x = 3.0
numerical = numerical_derivative(f, x)
analytical = 3 * x**2 + 2
print(f"Numerical: {numerical:.6f}")   # 29.000000
print(f"Analytical: {analytical:.6f}")  # 29.000000 ✓

# Partial derivatives of f(x,y) = x²y + 3xy²
def f_xy(x, y):
    return x**2 * y + 3 * x * y**2

# ∂f/∂x at (2, 3)
df_dx = (f_xy(2 + 1e-7, 3) - f_xy(2 - 1e-7, 3)) / (2e-7)
print(f"∂f/∂x at (2,3) = {df_dx:.2f}")  # 2*2*3 + 3*9 = 39

Chain Rule, Gradients & Multivariate Calculus

The Chain Rule: Core of Backpropagation

If y = f(g(x)), then dy/dx = f'(g(x)) · g'(x). This lets us differentiate compositions of functions, which is exactly what a neural network is — layers of composed functions.

A 3-layer neural network computes: y = f₃(f₂(f₁(x))). The chain rule gives us:

This is backpropagation — computing gradients layer by layer from output back to input.

Python
import numpy as np

# Example: y = sin(x²)
# dy/dx = cos(x²) · 2x  (chain rule: outer' × inner')
x = 2.0
analytical = np.cos(x**2) * 2 * x

# Numerical verification
h = 1e-7
numerical = (np.sin((x+h)**2) - np.sin((x-h)**2)) / (2*h)
print(f"Analytical: {analytical:.6f}")
print(f"Numerical:  {numerical:.6f}")  # Match! ✓

# Backprop through a simple 2-layer network
def forward_and_backward(x, w1, w2):
    # Forward pass
    z1 = w1 * x           # Layer 1
    a1 = np.maximum(0, z1) # ReLU
    z2 = w2 * a1           # Layer 2
    loss = (z2 - 1.0)**2   # MSE loss (target = 1.0)

    # Backward pass (chain rule!)
    dL_dz2 = 2 * (z2 - 1.0)       # ∂L/∂z2
    dL_dw2 = dL_dz2 * a1          # ∂L/∂w2 = ∂L/∂z2 · ∂z2/∂w2
    dL_da1 = dL_dz2 * w2          # ∂L/∂a1 = ∂L/∂z2 · ∂z2/∂a1
    dL_dz1 = dL_da1 * (z1 > 0)   # ∂L/∂z1 (ReLU derivative)
    dL_dw1 = dL_dz1 * x           # ∂L/∂w1

    return loss, dL_dw1, dL_dw2

loss, grad_w1, grad_w2 = forward_and_backward(2.0, 0.5, 0.3)
print(f"Loss: {loss:.4f}, ∂L/∂w1: {grad_w1:.4f}, ∂L/∂w2: {grad_w2:.4f}")

The Gradient Vector

The gradient ∇f collects all partial derivatives into a vector. It points in the direction of steepest ascent. Moving in the negative gradient direction = steepest descent = fastest way to reduce the loss.

Python
# Gradient of f(x,y) = x² + y² (a bowl)
def f(xy):
    return xy[0]**2 + xy[1]**2

def gradient(xy):
    return np.array([2*xy[0], 2*xy[1]])

# Gradient descent: start at (5, 3), walk downhill
pos = np.array([5.0, 3.0])
lr = 0.1

for i in range(20):
    grad = gradient(pos)
    pos = pos - lr * grad  # Step opposite to gradient
    if i % 5 == 0:
        print(f"Step {i:2d}: pos=({pos[0]:.3f}, {pos[1]:.3f}), f={f(pos):.4f}")
# Converges to (0, 0) — the minimum!

Jacobian, Hessian Matrices & Optimization via Calculus

The Jacobian Matrix

When a function maps multiple inputs to multiple outputs, the Jacobian is the matrix of all first-order partial derivatives. It generalizes the gradient to vector-valued functions.

Python
import numpy as np

# f(x,y) = [x²+y, xy+y²]
# Jacobian at (1, 2):
# J = [[∂f₁/∂x, ∂f₁/∂y],   = [[2x, 1 ],   = [[2, 1],
#      [∂f₂/∂x, ∂f₂/∂y]]     [ y, x+2y]]     [2, 5]]

def jacobian_numerical(f, x, h=1e-7):
    n = len(x)
    f0 = f(x)
    m = len(f0)
    J = np.zeros((m, n))
    for j in range(n):
        x_plus = x.copy(); x_plus[j] += h
        x_minus = x.copy(); x_minus[j] -= h
        J[:, j] = (f(x_plus) - f(x_minus)) / (2 * h)
    return J

f = lambda x: np.array([x[0]**2 + x[1], x[0]*x[1] + x[1]**2])
J = jacobian_numerical(f, np.array([1.0, 2.0]))
print("Jacobian:\n", J.round(4))

The Hessian Matrix

The Hessian is the matrix of second-order partial derivatives. It tells you about the curvature of the loss surface — is it a bowl (convex, minimum), a saddle, or a ridge?

If all eigenvalues of H are positive → local minimum (bowl). If all negative → local maximum. If mixed → saddle point. This information is used by second-order optimizers like Newton's method.

Python
# Hessian of f(x,y) = x⁴ + y⁴ - 2x²y² at (1,1)
# ∂²f/∂x² = 12x² - 4y² = 8
# ∂²f/∂y² = 12y² - 4x² = 8
# ∂²f/∂x∂y = -4xy = -4 (mixed partial, symmetric by Schwarz)

H = np.array([[8, -4],
              [-4, 8]])
eigenvalues = np.linalg.eigvalsh(H)
print("Eigenvalues:", eigenvalues)  # [4, 12] — both positive → minimum

Probability Distributions & Bayes' Theorem

Probability Distributions

A probability distribution describes the likelihood of each possible outcome. In ML, data is assumed to come from some underlying distribution, and the model tries to learn it.

Key Discrete Distributions

The Gaussian (Normal) Distribution

The most important continuous distribution. Data in nature often follows a bell curve. Neural network weights are initialized from Gaussians.

Python
import numpy as np
import matplotlib.pyplot as plt

# Generate samples from different distributions
normal_data = np.random.normal(loc=0, scale=1, size=10000)
uniform_data = np.random.uniform(low=-3, high=3, size=10000)

# Neural network weight initialization
# Xavier/Glorot: Normal with std = √(2 / (fan_in + fan_out))
fan_in, fan_out = 256, 128
weights = np.random.normal(0, np.sqrt(2.0 / (fan_in + fan_out)), size=(fan_in, fan_out))
print(f"Weight mean: {weights.mean():.4f}, std: {weights.std():.4f}")

Bayes' Theorem

Bayes' theorem lets us update our beliefs when we observe new evidence. It reverses conditional probabilities — going from P(evidence|hypothesis) to P(hypothesis|evidence).

Posterior = (Likelihood × Prior) / Evidence

Python
# Medical diagnosis example
# Disease prevalence: 1% of population has disease
# Test sensitivity: 99% (positive if sick)
# Test specificity: 95% (negative if healthy)
# Question: You test positive. What's the probability you're actually sick?

p_disease = 0.01         # Prior: P(D)
p_pos_given_disease = 0.99 # Likelihood: P(+|D)
p_pos_given_healthy = 0.05 # False positive: P(+|¬D)

# P(+) = P(+|D)P(D) + P(+|¬D)P(¬D)
p_positive = p_pos_given_disease * p_disease + p_pos_given_healthy * (1 - p_disease)

# P(D|+) = P(+|D)P(D) / P(+)
p_disease_given_pos = (p_pos_given_disease * p_disease) / p_positive

print(f"P(disease | positive test) = {p_disease_given_pos:.2%}")
# Only ~16.7%! Despite 99% sensitive test!
# This is because the disease is rare (low prior)

Expectation, Variance & Maximum Likelihood Estimation

Expectation & Variance

Expectation is the weighted average — the "center" of the distribution. Variance measures spread — how far values typically deviate from the mean. Standard deviation σ = √Var is in the same units as the data.

Python
import numpy as np

data = np.array([72, 85, 90, 78, 95, 88, 76, 92])

print(f"Mean (E[X]):     {np.mean(data):.2f}")
print(f"Variance:        {np.var(data):.2f}")
print(f"Std deviation:   {np.std(data):.2f}")

# Covariance — how two variables move together
height = np.array([165, 170, 175, 180, 185])
weight = np.array([55, 62, 70, 78, 85])
cov_matrix = np.cov(height, weight)
print(f"\nCovariance matrix:\n{cov_matrix}")
print(f"Correlation: {np.corrcoef(height, weight)[0,1]:.4f}")

Maximum Likelihood Estimation (MLE)

MLE finds the parameter values that make the observed data most probable. It asks: "Given this data, which parameters θ maximize the probability of seeing this data?"

For computational convenience (products → sums), we maximize the log-likelihood:

Python
import numpy as np

# MLE for Gaussian: find μ and σ that maximize P(data)
data = np.array([2.1, 1.8, 2.5, 2.3, 1.9, 2.7, 2.0, 2.4])

# MLE solutions (derivable analytically):
mu_mle = np.mean(data)     # MLE for μ = sample mean
sigma_mle = np.std(data)   # MLE for σ = sample std (biased)

print(f"MLE μ = {mu_mle:.3f}")
print(f"MLE σ = {sigma_mle:.3f}")

# Connection to ML: minimizing cross-entropy loss
# = maximizing log-likelihood of correct labels!
# So training a classifier with cross-entropy IS doing MLE.

Information Theory & Sampling Methods (MCMC)

Entropy: Measuring Uncertainty

Entropy measures the average "surprise" or uncertainty in a probability distribution. A fair coin (p=0.5) has maximum entropy — you can't predict the outcome. A biased coin (p=0.99) has low entropy — almost always heads.

Python
import numpy as np

def entropy(probs):
    probs = np.array(probs)
    probs = probs[probs > 0]  # Avoid log(0)
    return -np.sum(probs * np.log2(probs))

print(f"Fair coin:   {entropy([0.5, 0.5]):.3f} bits")    # 1.000
print(f"Biased coin: {entropy([0.99, 0.01]):.3f} bits")  # 0.081
print(f"Fair die:    {entropy([1/6]*6):.3f} bits")       # 2.585
print(f"Certain:     {entropy([1.0, 0.0]):.3f} bits")    # 0.000

KL Divergence

KL divergence measures how different two distributions P and Q are. It's the "extra bits" needed to encode data from P using a code optimized for Q. Used in VAEs, GANs, knowledge distillation.

Python
def kl_divergence(p, q):
    p, q = np.array(p), np.array(q)
    mask = p > 0
    return np.sum(p[mask] * np.log(p[mask] / q[mask]))

P = [0.4, 0.3, 0.3]   # True distribution
Q = [0.33, 0.33, 0.34] # Approximate distribution

print(f"KL(P||Q) = {kl_divergence(P, Q):.4f}")
print(f"KL(Q||P) = {kl_divergence(Q, P):.4f}")  # Different! KL is asymmetric

Cross-Entropy Loss

Cross-entropy = Entropy + KL divergence. Minimizing cross-entropy between true labels and predictions is equivalent to minimizing KL divergence (since entropy of true labels is constant).

Python
# Binary cross-entropy — the standard classification loss
def binary_cross_entropy(y_true, y_pred):
    epsilon = 1e-15
    y_pred = np.clip(y_pred, epsilon, 1 - epsilon)
    return -np.mean(y_true * np.log(y_pred) + (1-y_true) * np.log(1-y_pred))

y_true = np.array([1, 0, 1, 1, 0])
y_good = np.array([0.9, 0.1, 0.8, 0.95, 0.05])   # Good predictions
y_bad  = np.array([0.5, 0.5, 0.5, 0.5, 0.5])    # Random predictions

print(f"Good model loss: {binary_cross_entropy(y_true, y_good):.4f}")
print(f"Bad model loss:  {binary_cross_entropy(y_true, y_bad):.4f}")

Sampling Methods & MCMC

When we can't compute probabilities analytically (which is common in complex Bayesian models), we sample from the distribution instead. Monte Carlo methods use random sampling to approximate expectations.

Markov Chain Monte Carlo (MCMC) constructs a Markov chain whose stationary distribution equals the target distribution. The Metropolis-Hastings algorithm is the classic approach:

Python
import numpy as np

# Metropolis-Hastings to sample from a target distribution
# Target: mixture of two Gaussians
def target_pdf(x):
    return 0.3 * np.exp(-0.5 * (x + 2)**2) + 0.7 * np.exp(-0.5 * (x - 3)**2)

def metropolis_hastings(target, n_samples=10000, step_size=1.0):
    samples = []
    x = 0.0  # Starting point

    for _ in range(n_samples):
        # Propose new point
        x_new = x + np.random.normal(0, step_size)
        # Accept/reject based on probability ratio
        acceptance = target(x_new) / target(x)
        if np.random.random() < acceptance:
            x = x_new
        samples.append(x)

    return np.array(samples)

samples = metropolis_hastings(target_pdf, n_samples=50000)
print(f"Sample mean: {samples[1000:].mean():.3f}")  # Discard burn-in
print(f"Sample std:  {samples[1000:].std():.3f}")

Convex vs Non-Convex Optimization

What is Convexity?

A function f is convex if the line segment between any two points on the curve lies above the curve. Geometrically, it's a "bowl" shape. The critical property: any local minimum is also the global minimum. No getting stuck in suboptimal solutions.

Python
import numpy as np
import matplotlib.pyplot as plt

# Convex function: f(x) = x² — single global minimum
x = np.linspace(-5, 5, 200)
convex = x**2

# Non-convex function: multiple local minima
non_convex = x**4 - 8*x**2 + 5*np.sin(3*x)

# In deep learning, loss surfaces look more like the non-convex one
# with millions of dimensions — many local minima and saddle points

Convex vs Non-Convex in ML

Gradient Descent Variants & Loss Functions

Gradient Descent Variants

Python
import numpy as np

# Mini-batch gradient descent
def mini_batch_gd(X, y, lr=0.01, epochs=100, batch_size=32):
    n_samples, n_features = X.shape
    w = np.zeros(n_features)
    b = 0

    for epoch in range(epochs):
        # Shuffle data each epoch
        indices = np.random.permutation(n_samples)
        X_shuffled = X[indices]
        y_shuffled = y[indices]

        for i in range(0, n_samples, batch_size):
            X_batch = X_shuffled[i:i+batch_size]
            y_batch = y_shuffled[i:i+batch_size]

            y_pred = X_batch @ w + b
            error = y_pred - y_batch

            # Gradients averaged over batch
            dw = (2 / len(y_batch)) * X_batch.T @ error
            db = (2 / len(y_batch)) * np.sum(error)

            w -= lr * dw
            b -= lr * db

    return w, b

Advanced Optimizers

Momentum

Adds a "velocity" to gradient descent, accumulating past gradients like a ball rolling downhill. Helps escape shallow local minima and smooths noisy gradients.

Adam (Adaptive Moment Estimation)

Combines momentum (1st moment) with RMSProp (2nd moment) — tracks both the average gradient direction AND the average gradient magnitude. The most widely used optimizer in deep learning.

Python
class Adam:
    def __init__(self, lr=0.001, beta1=0.9, beta2=0.999, eps=1e-8):
        self.lr = lr
        self.beta1, self.beta2, self.eps = beta1, beta2, eps
        self.m = None   # 1st moment (mean of gradients)
        self.v = None   # 2nd moment (mean of squared gradients)
        self.t = 0

    def step(self, w, grad):
        self.t += 1
        if self.m is None:
            self.m = np.zeros_like(w)
            self.v = np.zeros_like(w)

        self.m = self.beta1 * self.m + (1 - self.beta1) * grad
        self.v = self.beta2 * self.v + (1 - self.beta2) * grad**2

        # Bias correction
        m_hat = self.m / (1 - self.beta1**self.t)
        v_hat = self.v / (1 - self.beta2**self.t)

        w -= self.lr * m_hat / (np.sqrt(v_hat) + self.eps)
        return w

Loss Functions

Lagrangian Optimization, Saddle Points & Local Minima

Lagrangian Optimization

When you need to optimize a function subject to constraints, you can't just use gradient descent. The Lagrangian method converts a constrained problem into an unconstrained one by introducing multiplier variables.

At the optimum, ∇ₓL = 0 and ∇λL = 0 (the constraint is satisfied).

Python
import numpy as np
from scipy.optimize import minimize

# Example: Minimize f(x,y) = x² + y² subject to x + y = 4
# Lagrangian: L = x² + y² + λ(x + y - 4)
# ∂L/∂x = 2x + λ = 0 → x = -λ/2
# ∂L/∂y = 2y + λ = 0 → y = -λ/2
# Constraint: x + y = 4 → -λ/2 - λ/2 = 4 → λ = -4
# Solution: x = 2, y = 2, f(2,2) = 8

# Numerical verification with scipy
result = minimize(
    lambda xy: xy[0]**2 + xy[1]**2,
    x0=[0, 0],
    constraints={'type': 'eq', 'fun': lambda xy: xy[0] + xy[1] - 4}
)
print(f"Optimal: x={result.x[0]:.2f}, y={result.x[1]:.2f}")
print(f"Minimum value: {result.fun:.2f}")

Saddle Points & Local Minima

A saddle point is where the gradient is zero but it's a minimum along some directions and a maximum along others — like a mountain pass. In a function of n variables, a saddle point has some positive and some negative Hessian eigenvalues.

Python
# Classic saddle point: f(x,y) = x² - y²
# ∇f = [2x, -2y] = [0,0] at origin
# Hessian = [[2, 0], [0, -2]]
# Eigenvalues: 2 and -2 → mixed signs → SADDLE POINT

# At the origin: minimum along x-axis, maximum along y-axis

H_saddle = np.array([[2, 0], [0, -2]])
eigenvalues = np.linalg.eigvalsh(H_saddle)
print(f"Eigenvalues: {eigenvalues}")
print("Mixed signs → Saddle point!")

Why Saddle Points Dominate in High Dimensions

At a critical point (∇f = 0), each Hessian eigenvalue is randomly positive or negative. For a saddle point, you need at least one of each sign. For a local minimum, ALL must be positive. In n dimensions:

For n = 1000 parameters: P ≈ 10⁻³⁰¹ — essentially zero! Saddle points are overwhelmingly more likely than true local minima in high-dimensional spaces. This is why neural networks can still find good solutions despite non-convexity.

Escaping Saddle Points

Several techniques help optimizers navigate past saddle points:

• SGD noise: Stochastic gradients provide natural perturbation that pushes away from saddle points
• Momentum: Accumulated velocity carries the optimizer through flat saddle regions
• Learning rate warmup: Starts with small steps to avoid early saddle-point traps
• Batch normalization: Smooths the loss landscape, reducing saddle point severity