Chapter 3: Mathematical Foundations

A matrix transforms vectors by mapping each input to an output. The columns of a matrix show where the basis vectors (x and y axes) end up after transformation.

Key insight: Any linear transformation—rotation, scaling, shearing—can be represented as a matrix. This is why matrices are fundamental to computer graphics and vision.

What makes a transformation linear?

• $T(\vec{u} + \vec{v}) = T(\vec{u}) + T(\vec{v})$ (additivity)
• $T(c\vec{v}) = cT(\vec{v})$ (scaling)

Common Transformations:

• Scaling: Diagonal entries stretch/compress axes
• Rotation: Preserves lengths and angles
• Shear: Skews one axis relative to another
• Reflection: Flips across an axis (determinant = -1)

Matrix multiplication composes transformations. If $A$ rotates and $B$ scales, then $AB$ does both (scale first, then rotate).

Important: Matrix multiplication is not commutative—$AB \neq BA$ in general. The order of transformations matters!

Dimensions: For $A$ (m×n) and $B$ (n×p), the result $AB$ is (m×p). Inner dimensions must match.

Properties:

• Associative: $(AB)C = A(BC)$ — compute in any order
• Distributive: $A(B + C) = AB + AC$
• Transpose: $(AB)^T = B^T A^T$ — order reverses

Computational Cost: Multiplying two $n \times n$ matrices naively costs $O(n^3)$. Optimized algorithms like Strassen's achieve $O(n^{2.37})$, which matters for large neural network layers.

The rotation matrix rotates points counterclockwise by angle $\theta$. Key properties:

• Orthogonal: $R^T R = I$ (columns are perpendicular unit vectors)
• Determinant = 1: Preserves area/volume (proper rotation)
• Inverse = Transpose: $R^{-1} = R^T$ (rotate by $-\theta$)
• Composition: $R(\alpha)R(\beta) = R(\alpha + \beta)$

In 3D: Rotation matrices are 3×3. Common representations:

• Euler angles: Three successive rotations (but suffer from gimbal lock)
• Axis-angle: Rotate by $\theta$ around unit vector $\hat{n}$
• Quaternions: 4D representation avoiding gimbal lock, used in 3D vision

Rodrigues' Formula converts axis-angle to matrix: $R = I + (\sin\theta)K + (1-\cos\theta)K^2$ where $K$ is the skew-symmetric matrix of the axis.

The inverse $A^{-1}$ undoes transformation $A$. If $A$ rotates 30°, $A^{-1}$ rotates -30°.

When it exists: Only square matrices with non-zero determinant are invertible. Singular matrices (det = 0) have no inverse—they collapse space into lower dimensions.

For 2×2: $A^{-1} = \frac{1}{ad-bc}\begin{pmatrix} d & -b \\ -c & a \end{pmatrix}$

Pseudoinverse: For non-square or singular matrices, the Moore-Penrose pseudoinverse $A^+$ provides a best-fit solution: $A^+ = (A^T A)^{-1} A^T$ Used in least squares problems and solving overdetermined systems.

Numerical Considerations: Direct inversion is often unstable. In practice, use:

• LU decomposition for general systems
• Cholesky for symmetric positive definite matrices
• SVD for ill-conditioned matrices

Linearity is the property that makes a transformation predictable from very little data. The formula says $T$ respects linear combinations: transforming $\alpha\mathbf{u} + \beta\mathbf{v}$ gives the same result as transforming $\mathbf{u}$ and $\mathbf{v}$ separately and combining them with the same scalars $\alpha$ and $\beta$. Consequently $T$ is fully determined by where it sends the basis vectors, and those images become the columns of its matrix. Two quick tests catch most non-examples. A linear map must send $\mathbf{0}$ to $\mathbf{0}$, so translation by a fixed offset fails, and it cannot contain squares or products of inputs. A layer $W\mathbf{x}$ is linear; adding a bias makes it affine.

A scaling matrix is diagonal: $s_x$ multiplies every x-coordinate and $s_y$ every y-coordinate, and the two axes never mix. Factors above 1 stretch, factors between 0 and 1 compress, and a negative factor also reflects across the other axis. The determinant $s_x s_y$ is the factor by which areas change. Equal factors give uniform scaling, which preserves shapes; unequal factors distort circles into ellipses. Standardizing features, dividing each column by its standard deviation, is exactly a diagonal scaling of the data. A zero factor collapses an axis entirely, which makes the matrix singular and destroys the information along that axis.

Projection keeps the part of a vector that lies along a chosen direction and discards the rest. For a unit vector $\mathbf{u}$, the matrix $P = \mathbf{u}\mathbf{u}^T$ sends $\mathbf{x}$ to $(\mathbf{u}^T\mathbf{x})\mathbf{u}$: the dot product measures how far $\mathbf{x}$ reaches along $\mathbf{u}$, and multiplying by $\mathbf{u}$ points the result in that direction. Two properties identify it: $P^2 = P$, since projecting twice changes nothing, and $P^T = P$. Its rank is 1 and its determinant is 0, so it cannot be undone. The formula requires unit length; for any other $\mathbf{u}$, divide by $\mathbf{u}^T\mathbf{u}$. PCA projects data onto its principal directions this way.

The determinant condenses a square matrix into one number that says how the transformation scales space. Its absolute value is the factor by which areas, in 2D, or volumes, in 3D, are multiplied, and its sign says whether orientation is preserved or flipped. For a 2×2 matrix [[a, b], [c, d]] it equals $ad - bc$. A determinant of zero means some direction is crushed to nothing, so the matrix is singular and has no inverse. It is multiplicative, $\det(AB) = \det(A)\det(B)$, matching the composition of area scalings, and it equals the product of the eigenvalues. Its size is a poor measure of near-singularity, though; the condition number measures that properly.

Orthogonal matrices move space rigidly, without stretching it. The formula $Q^TQ = I$ says each column has length 1 and every pair of columns is perpendicular, and it gives the cheapest possible inverse, $Q^{-1} = Q^T$. Because $\|Q\mathbf{x}\| = \|\mathbf{x}\|$ and dot products are preserved, lengths and angles survive, and rounding errors are not amplified: the condition number is 1. The determinant is +1 for rotations and -1 for reflections, and every eigenvalue has absolute value 1. Orthogonal weight initialization uses these properties to keep signal norms stable across layers. A matrix whose columns are orthogonal but not unit length is not an orthogonal matrix.

The normal, or Gaussian, distribution is the bell curve, and it earns its central place through the Central Limit Theorem: sums and averages of many independent effects tend toward it, whatever their individual distributions. In the formula, $\mu$ is the mean, which sets the center, $\sigma$ is the standard deviation, which sets the width, and the factor $1/(\sigma\sqrt{2\pi})$ makes the total area 1. The exponent penalizes squared distance from the mean, so values beyond three standard deviations are rare. Assuming Gaussian noise is what turns maximum likelihood into least squares. Real data often has heavier tails, however, and a Gaussian model then badly underestimates extreme events.

Bayes' theorem reverses a conditional probability. We usually know how likely the evidence $B$ is under a hypothesis $A$, the likelihood $P(B|A)$, but we want the probability of the hypothesis after seeing the evidence, the posterior $P(A|B)$. The theorem multiplies the likelihood by the prior $P(A)$, our belief before the evidence, and divides by $P(B)$, the total probability of the evidence, found by summing over all hypotheses. The classic mistake, base-rate neglect, is to treat $P(B|A)$ as if it were $P(A|B)$. Naive Bayes classifiers, Bayesian inference, and probabilistic graphical models are all built on this single equation.

The mean $\mu$ is the balance point of the data, the sum of the values $x_i$ divided by their count $n$. The variance $\sigma^2$ is the average squared distance from that mean, and the standard deviation $\sigma$, its square root, returns the spread to the original units. Squaring makes every deviation positive and weights large deviations heavily, so a single outlier can inflate the variance a lot. Dividing by $n$ gives the population variance; dividing by $n - 1$ gives the unbiased sample estimate, and NumPy's var divides by $n$ unless you pass ddof=1. Standardization subtracts the mean and divides by the standard deviation, and batch normalization applies it inside networks.

Covariance measures whether two variables move together. In the formula, each pair contributes the product of its deviations from the means $\mu_X$ and $\mu_Y$: positive when both sit on the same side of their means, negative when they sit on opposite sides. Averaging those products gives $\text{Cov}(X, Y)$. Its size depends on units, so the correlation $\rho$ divides by both standard deviations and always lands in [-1, 1]. Values near 1 or -1 mean a nearly perfect linear relationship. A value of 0 means no linear relationship, not independence: $Y = X^2$ with symmetric $X$ has zero correlation. Collecting all pairwise covariances gives the covariance matrix used by PCA.

Maximum likelihood chooses the parameters $\theta$ under which the observed data is most probable. For independent samples $x_i$, the likelihood is the product of the individual probabilities $p(x_i | \theta)$, and the estimate $\hat{\theta}$ is the value that maximizes it. Products of many small numbers underflow, so in practice we maximize the log-likelihood, a sum, which has the same maximizer because log is increasing. Minimizing cross-entropy, and minimizing squared error under Gaussian noise, are both MLE in disguise. With little data MLE overfits: after three heads in three flips it estimates heads with probability 1. Adding a prior gives MAP estimation, which acts as regularization.

The expectation is the long-run average of a random variable: each possible value $x$ is weighted by its probability $P(x)$, summed for discrete variables and integrated against the density $p(x)$ for continuous ones. It need not be a value the variable can actually take. Expectation is linear, $E[aX + b] = aE[X] + b$ and $E[X + Y] = E[X] + E[Y]$, even for dependent variables, which makes it easy to compute. Training minimizes expected loss over the data distribution, estimated by an average over a mini-batch. The trap is nonlinearity: $E[f(X)]$ generally differs from $f(E[X])$, and the gap between $E[X^2]$ and $E[X]^2$ is exactly the variance.

Entropy measures how unpredictable a distribution is. Each outcome $x$ carries surprise $-\log p(x)$: a certain event carries none, and a one-in-a-million event carries a lot. Entropy $H(X)$ averages this surprise, weighting each outcome by its probability $p(x)$. With base-2 logs it is the minimum average number of bits per symbol that any lossless code can achieve, Shannon's source coding theorem. It is zero for a deterministic outcome and largest, $\log_2 n$ bits for $n$ outcomes, for the uniform distribution. Outcomes with zero probability contribute nothing, since $p \log p$ tends to 0, so code must guard against evaluating log 0 directly.

Cross-entropy scores a model distribution $q$ against the true distribution $p$. It averages the surprise $-\log q(x)$ that the model assigns to each outcome, weighted by how often that outcome really occurs, $p(x)$. It is never smaller than the entropy $H(p)$, and it equals $H(p)$ only when $q = p$. For classification with a one-hot label, every term vanishes except the true class, so the loss is simply the negative log-probability the model gives to the correct answer. Minimizing it over a dataset is maximum likelihood. If the model gives the true class probability 0, the loss is infinite, which is why implementations clip probabilities or work with logits.

KL divergence measures the extra cost of describing data from $p$ with a model $q$. The formula shows it equals the cross-entropy $H(p, q)$ minus the entropy $H(p)$, which no model can remove. It is always at least 0 and equals 0 only when $p = q$, so it behaves like a distance, but it is not one, because it is asymmetric. The forward direction heavily penalizes a $q$ that puts low probability where $p$ has mass, while the reverse direction penalizes $q$ for mass where $p$ has none. When $p$ is fixed, minimizing KL and minimizing cross-entropy are equivalent. VAEs, knowledge distillation, and RLHF penalties all use it.

Mutual information measures how much learning one variable tells you about another. The first form says $I(X; Y)$ is the entropy of $X$ minus its conditional entropy $H(X|Y)$, the uncertainty that remains once $Y$ is known. The second says it is the KL divergence between the joint distribution $p(x, y)$ and the product of marginals $p(x)p(y)$, which is what the joint would be under independence. So it is symmetric, never negative, and zero exactly when $X$ and $Y$ are independent. Unlike correlation, it detects nonlinear dependence. It is used for feature selection and in contrastive learning, where the InfoNCE loss bounds it from below.

The unit of entropy depends only on the base of the logarithm. Base 2 gives bits, the natural fit for binary codes; base $e$ gives nats. The conversion follows from $\ln x = \ln 2 \cdot \log_2 x$: one bit is $\ln 2 \approx 0.693$ nats, and one nat is $1/\ln 2 \approx 1.443$ bits, as the formula states. Machine-learning libraries use natural logs, because the derivative of $\ln x$ is simply $1/x$ and softmax is defined with the exponential. So a cross-entropy loss printed by PyTorch is in nats. A common mistake is to compare a loss in nats with a compression rate or benchmark figure reported in bits.

Binary cross-entropy is cross-entropy for two classes. The label $y$ is 0 or 1, and the model outputs $\hat{y}$, the predicted probability of class 1, usually through a sigmoid. When $y = 1$ only the first term survives, giving a loss of $-\log \hat{y}$; when $y = 0$ only the second survives, giving $-\log(1 - \hat{y})$. The loss is near 0 for confident correct predictions and grows without bound for confident wrong ones. Combined with a sigmoid, its gradient with respect to the logit is simply $\hat{y} - y$. Implementations should use the logit form, such as BCEWithLogitsLoss, because computing $\log(1 - \hat{y})$ for $\hat{y}$ near 1 loses precision.

The derivative $f'(x)$ measures the rate of change of $f$ at point $x$. Geometrically, it's the slope of the tangent line.

• Positive derivative: Function is increasing
• Negative derivative: Function is decreasing
• Zero derivative: Local minimum or maximum (critical point)

Common Derivative Rules:

• Power rule: $\frac{d}{dx}x^n = nx^{n-1}$
• Product rule: $(fg)' = f'g + fg'$
• Quotient rule: $(f/g)' = (f'g - fg')/g^2$

Second Derivative: $f''(x)$ measures curvature. At critical points:

• $f'' > 0$: Local minimum (concave up)
• $f'' < 0$: Local maximum (concave down)
• $f'' = 0$: Inflection point or saddle point

The gradient generalizes the derivative to multiple variables. For $f(x, y)$:

• $\nabla f$ points in the direction of steepest increase
• $|\nabla f|$ gives the rate of steepest increase
• $-\nabla f$ points toward steepest decrease

Directional Derivative: Rate of change in direction $\vec{u}$: $D_{\vec{u}}f = \nabla f \cdot \vec{u}$

The Hessian Matrix: Matrix of second partial derivatives: $H = \begin{pmatrix} f_{xx} & f_{xy} \\ f_{yx} & f_{yy} \end{pmatrix}$

Eigenvalues of the Hessian reveal curvature: all positive = minimum, all negative = maximum, mixed = saddle point.

Image Gradients: Sobel/Prewitt operators compute discrete gradients. Edge magnitude $|\nabla I|$ and direction $\theta = \arctan(I_y/I_x)$ are fundamental to edge detection.

Gradient descent iteratively minimizes a function by taking steps opposite to the gradient:

• $\alpha$ (learning rate): Step size—too large causes overshooting, too small is slow
• Convergence: Stops when gradient is near zero (found a minimum)
• Local vs global: May find local minimum, not global

Variants Used in Deep Learning:

• SGD: Use mini-batch gradients for efficiency
• Momentum: Accumulate velocity to escape local minima
• Adam: Adaptive learning rates per parameter
• AdaGrad/RMSprop: Scale learning rate by historical gradients

Learning Rate Schedules: Start large for fast progress, decay over time for fine convergence. Common: step decay, cosine annealing, warmup.

Convergence Guarantees: For convex functions with Lipschitz gradients, GD converges at rate $O(1/k)$. Adding momentum can achieve $O(1/k^2)$.

The chain rule handles derivatives of composed functions. If $h(x) = f(g(x))$, then:

$h'(x) = f'(g(x)) \cdot g'(x)$

Multivariate Chain Rule: For $z = f(x, y)$ where $x = x(t)$, $y = y(t)$: $\frac{dz}{dt} = \frac{\partial f}{\partial x}\frac{dx}{dt} + \frac{\partial f}{\partial y}\frac{dy}{dt}$

Backpropagation: Neural networks are compositions of layers. For loss $L = f_n(f_{n-1}(...f_1(x)))$: $\frac{\partial L}{\partial W_i} = \frac{\partial L}{\partial a_n} \cdot \frac{\partial a_n}{\partial a_{n-1}} \cdots \frac{\partial a_{i+1}}{\partial a_i} \cdot \frac{\partial a_i}{\partial W_i}$

Automatic Differentiation: Modern frameworks (PyTorch, JAX) implement the chain rule automatically by building a computational graph and propagating gradients backward through it.

The Hessian collects every second partial derivative of $f$: entry $H_{ij}$ measures how the slope along $x_i$ changes as you move along $x_j$. For smooth functions the order of differentiation does not matter, so $H$ is symmetric and has real eigenvalues, one curvature per principal direction. At a critical point, where the gradient is zero, all positive eigenvalues mean a local minimum, all negative a maximum, and mixed signs a saddle. The ratio of largest to smallest eigenvalue, the condition number, predicts how badly gradient descent zig-zags. Newton's method steps by $-H^{-1}\nabla f$, but with $n$ parameters $H$ has $n^2$ entries, which is prohibitive for large networks.