2. Mathematical Foundations

This chapter covers the mathematics that every later technique quietly depends on: linear algebra for representing and transforming data, calculus and optimization for fitting models, probability for reasoning about uncertainty, and information theory for measuring how much a signal tells you. You do not need to derive these results from scratch, but you do need to recognize them, because almost every modelling failure later in the book is a violation of an assumption stated here.

2.1 Linear Algebra

Linear algebra is the language of data representation: a dataset is a matrix, a model is often a transformation, and most efficient algorithms are matrix operations in disguise.

Vectors and Vector Spaces

A dataset row is a list of numbers: age 34, income 52000, tenure 19. Treating that list as a point in space, rather than three unrelated facts, is the move that makes the whole subject useful. Once rows are points, "similar customers" becomes "nearby points", "a trend" becomes "a direction", and "redundant features" becomes "points that lie on a lower-dimensional sheet".

A vector is an ordered list of numbers. Geometrically it is an arrow from the origin to a point, which gives it two readings: a location, and a direction with a magnitude. Both readings are used constantly. A data row is a location. A gradient is a direction.

A vector space is the set of all vectors you can build from a starting set using two operations: adding vectors, and scaling them by a number. Being closed under these operations is what makes the space well-behaved. The span of a set of vectors is everything reachable this way.

Two related ideas matter for data work. Vectors are linearly independent if none of them can be built from the others; if one can, it adds no new direction and is redundant. A basis is a minimal set that spans the space, and its size is the dimension.

An analogy: a basis is a set of ingredients. Flour, water, and salt span a range of doughs. Adding "flour again, but pre-measured" adds no new dough you could not already make, which is exactly what linear dependence means.

  VECTORS AS ARROWS AND AS DATA ROWS

        y
        ^
      4 |          * b = (3, 4)
        |         /|
      3 |        / |
        |       /  |          a + b = (5, 5)
      2 |      /   |            . - - - - - *
        |     /    |          /
      1 |    /     |  * a = (2, 1)
        |   /   ___/__/
      0 +--/---/-------------------> x
        0  1  2   3   4   5

  Same object, two readings:
    location  -> "a customer with 2 years tenure and 1 support ticket"
    direction -> "move 2 right and 1 up"

  SPAN of {a} is the infinite line through a.
  SPAN of {a, b} is the whole plane, because a and b are independent.
  SPAN of {a, 2a} is still just the line: 2a adds nothing.

import numpy as np

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

print(a + b)        # [5 5]  -- tip-to-tail addition
print(3 * a)        # [6 3]  -- scaling changes magnitude, not direction
print(a @ b)        # 10     -- dot product: 2*3 + 1*4

# Independence check: stack as rows and ask for the rank.
# rank = number of genuinely distinct directions present.
print(np.linalg.matrix_rank(np.vstack([a, b])))      # 2 -> independent
print(np.linalg.matrix_rank(np.vstack([a, 2 * a])))  # 1 -> dependent

# Why this matters on real data: a redundant feature is a dependent column.
X = np.array([
    [1.0, 2.0, 3.0],     # feature 3 = feature 1 + feature 2, exactly
    [2.0, 1.0, 3.0],
    [4.0, 5.0, 9.0],
])
print(np.linalg.matrix_rank(X))   # 2, not 3 -- one column is redundant
# Linear regression on this X has no unique solution. The redundancy is
# not a nuisance to clean up later; it breaks the model outright.

Key Takeaways

Treating rows as points in space turns similarity, trend, and redundancy into geometric questions.

A vector reads both as a location (a data row) and as a direction (a gradient or a trend).

Linear dependence means a vector adds no new direction, which on a data matrix means a redundant feature.

The dimension of a space is the size of its basis, and matrix rank counts the distinct directions actually present.

🧪 Practice

By hand, determine whether (1, 2), (2, 4), and (0, 1) span the plane. Verify with matrix_rank.

Construct a 4-by-4 data matrix whose rank is 3, and say which column is redundant and why.

Interview: A colleague one-hot encodes a three-category feature into three columns and adds an intercept, and the regression fails. What happened? (Hint: sum the three dummy columns and compare to the intercept column.)

Matrix Operations

A matrix is a grid of numbers, but that description hides why it matters. A matrix does two jobs in data science, and keeping them separate prevents most confusion. First, it stores a dataset: rows are observations, columns are features. Second, it acts as a transformation: multiplying a vector by a matrix rotates, scales, projects, or otherwise moves it.

Matrix multiplication is the operation that makes the second reading work. The product AB is defined so that applying AB to a vector is the same as applying B first and then A. The mechanical rule follows from that: the entry at row i, column j of AB is the dot product of row i of A with column j of B. Inner dimensions must match, which is the single most common source of runtime errors in numerical code.

Two properties surprise newcomers. Multiplication is not commutative: AB and BA are usually different and may not both exist. It is associative, which is what allows the optimization of reordering a chain of products.

An analogy: matrices are verbs, not nouns. "Rotate then stretch" is a different instruction from "stretch then rotate", which is exactly why order matters.

  SHAPE ALGEBRA -- the rule that prevents most errors

     A          B          AB
  (m x n) @ (n x p)  =  (m x p)
        \      /
         must match

  Example: predictions = X @ w
     X: (1000 x 5)   1000 rows, 5 features
     w: (5 x 1)      one weight per feature
     ->  (1000 x 1)  one prediction per row

  A MATRIX AS A TRANSFORMATION

   unit square        after A = [[2, 1],      sheared and
                                 [0, 1]]      stretched
   (0,1) +----+ (1,1)             (1,1) +------+ (3,1)
         |    |          --->           /      /
         |    |                        /      /
   (0,0) +----+ (1,0)           (0,0) +------+ (2,0)

  The columns of A are where the basis vectors land:
  column 1 = (2,0) is the image of (1,0); column 2 = (1,1) is the image
  of (0,1). Reading a matrix by its columns tells you what it does.

import numpy as np

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

print(A @ B)        # [[5 1]  matrix product: apply B, then A
                    #  [3 1]]
print(B @ A)        # [[2 1]  DIFFERENT: order changes the result
                    #  [6 4]]
print(A * B)        # [[2 0]  element-wise, NOT matrix multiplication
                    #  [0 1]]  -- a frequent and silent bug

# The three multiplications above are all legal and all mean different
# things. Only `@` composes transformations.

X = np.random.default_rng(0).normal(size=(1000, 5))
w = np.array([0.5, -1.2, 0.0, 2.0, 0.3])

preds = X @ w                     # (1000,5) @ (5,) -> (1000,)
print(preds.shape)                # (1000,)

# Associativity is a real performance lever, not a curiosity:
# (X @ W1) @ W2  vs  X @ (W1 @ W2) give identical results at very
# different costs when the inner matrices are small.
W1 = np.random.default_rng(1).normal(size=(5, 5))
W2 = np.random.default_rng(2).normal(size=(5, 5))
np.allclose((X @ W1) @ W2, X @ (W1 @ W2))   # True -- but time them

Operation	Notation	Shape rule	Meaning
Addition	`A + B`	identical shapes	Element-wise sum
Scalar multiple	`c * A`	any	Scale every entry
Element-wise	`A * B`	identical (or broadcastable)	Hadamard product
Matrix product	`A @ B`	`(m,n) @ (n,p)`	Compose transformations
Transpose	`A.T`	`(m,n)` becomes `(n,m)`	Swap rows and columns
Matrix-vector	`A @ v`	`(m,n) @ (n,)`	Transform a vector

Key Takeaways

A matrix both stores a dataset and acts as a transformation; keep the two readings distinct.

Matrix multiplication composes transformations, so it is associative but not commutative.

The columns of a matrix are the images of the basis vectors, which is the fastest way to read what it does.

A * B and A @ B are both legal and mean different things, which makes the confusion silent rather than loud.

🧪 Practice

Multiply a (3,4) by a (4,2) matrix by hand for one entry, then check the whole product with NumPy.

Find a 2-by-2 matrix that rotates the plane 90 degrees, by deciding where (1,0) and (0,1) must land.

Interview: You must compute A @ B @ C where shapes are (1000,2), (2,1000), and (1000,2). Which grouping do you choose and why? (Hint: count the multiplications each order performs.)

Matrix Inverse and Rank

Solving a linear system means undoing a transformation: given Ax = b, find the x that A turned into b. The inverse A^-1 is the transformation that undoes A, so that A^-1 A x = x. When it exists, the solution is x = A^-1 b.

The catch is that most matrices have no inverse, and understanding when is more useful than knowing the formula. A transformation can be undone only if it loses no information. If A squashes two distinct inputs onto the same output, there is no way to recover which one you started from. Such a matrix is singular.

Rank is what measures this. The rank of a matrix is the number of linearly independent columns, which is the dimension of the space its outputs live in. A square n-by-n matrix is invertible exactly when its rank is n, called full rank. Lower rank means the transformation collapsed the space.

The practically important case is the near-singular matrix: technically invertible but numerically hopeless. The condition number measures this. A high condition number means tiny changes in b cause huge changes in x, so the solution is meaningless in floating point.

  RANK AS INFORMATION LOSS

  Full rank (rank 2)              Rank deficient (rank 1)
  A = [[2, 0],                    A = [[2, 4],
       [0, 3]]                         [1, 2]]

    plane -> plane                   plane -> line
                                      \
   +----+      +------+                \  every input lands on
   |    |  ->  |      |                 \ one line; infinitely many
   +----+      +------+                  \ inputs share an output
                                          \
  invertible: each output              singular: cannot invert,
  came from exactly one input          information was destroyed

  IN DATA TERMS
  rank deficiency in X'X means two features are (nearly) redundant,
  so the regression cannot decide how to split the credit between them.

import numpy as np

A = np.array([[2.0, 0.0],
              [0.0, 3.0]])
print(np.linalg.matrix_rank(A))      # 2 -> full rank, invertible
print(np.linalg.inv(A))              # [[0.5, 0], [0, 0.333]]

S = np.array([[2.0, 4.0],
              [1.0, 2.0]])           # column 2 = 2 * column 1
print(np.linalg.matrix_rank(S))      # 1 -> singular
# np.linalg.inv(S)  -> LinAlgError: Singular matrix

# The dangerous case is NEARLY singular, which raises no error at all.
N = np.array([[1.0, 1.0],
              [1.0, 1.0000001]])
print(f"condition number = {np.linalg.cond(N):.3e}")   # ~4e+07
b = np.array([2.0, 2.0])
x1 = np.linalg.solve(N, b)
x2 = np.linalg.solve(N, b + np.array([0.0, 1e-6]))     # tiny perturbation
print(x1, x2)     # wildly different answers from nearly identical inputs

# Practical rule: never invert to solve a system. `solve` is faster and
# numerically better than forming the inverse and multiplying.
x = np.linalg.solve(A, b)            # preferred
x_bad = np.linalg.inv(A) @ b         # works, but slower and less stable

# For rank-deficient or ill-conditioned systems, use the pseudo-inverse,
# which returns the minimum-norm least-squares solution instead of failing.
print(np.linalg.pinv(S) @ b)

Key Takeaways

A matrix is invertible exactly when it destroys no information, which is equivalent to having full rank.

Rank counts independent columns and therefore the true dimensionality of the transformation's output.

Near-singular matrices raise no error but produce meaningless solutions; check the condition number.

Use solve rather than forming an explicit inverse, and pinv when the system is rank deficient.

🧪 Practice

Build a 3-by-3 matrix of rank 2 and confirm that inv fails while pinv returns an answer.

Compute the condition number of a matrix whose two columns are nearly parallel, then increase the angle and watch it fall.

Interview: A regression on 40 features and 30 rows returns wildly large coefficients. What is going on structurally? (Hint: compare the number of rows to the number of columns and ask what the rank can possibly be.)

Norms and Distance Metrics

Almost every data science method needs to answer "how far apart are these two things?" Clustering groups nearby points, nearest-neighbour methods retrieve them, regularization penalizes large coefficients, and loss functions measure the distance between prediction and truth. All of these need a definition of size and distance, and crucially, the definition is a choice.

A norm measures the size of a single vector. A distance metric measures separation between two, and the most common ones are just the norm of the difference.

The main family is the p-norms. The L2 (Euclidean) norm is straight-line distance and is the default almost everywhere. The L1 (Manhattan) norm sums absolute values and behaves as if you must travel along a grid. The L-infinity norm takes the largest single component.

The choice has real consequences. L2 penalizes one large error far more than several small ones, because it squares. L1 treats them proportionally, which makes it robust to outliers and, when used as a regularizer, drives coefficients to exactly zero.

An analogy: L2 is the distance a bird flies, L1 is the distance a taxi drives, and L-infinity is the number of moves a chess king needs.

  UNIT BALLS: the set of points at distance 1 from the origin

     L1 (Manhattan)        L2 (Euclidean)       L-infinity (max)
                                    
         /\                    _--_                +------+
        /  \                  /    \               |      |
       <    >                |      |              |      |
        \  /                  \ _ _/               |      |
         \/                                        +------+

   diamond: corners on       circle: no            square: corners
   the axes                  preferred direction   on the diagonals

  WHY L1 PRODUCES SPARSITY
  The diamond's corners lie ON the axes. When a loss contour first
  touches the L1 ball, it very often touches at a corner, where one
  coordinate is exactly zero. The circle has no corners, so L2 shrinks
  coefficients toward zero without ever setting them to zero.

import numpy as np
from scipy.spatial import distance

v = np.array([3.0, -4.0])

print(np.linalg.norm(v, ord=1))        # 7.0   = |3| + |-4|
print(np.linalg.norm(v, ord=2))        # 5.0   = sqrt(9 + 16)
print(np.linalg.norm(v, ord=np.inf))   # 4.0   = max(|3|, |-4|)

# Different norms rank the same errors differently.
errors_spread   = np.array([1.0, 1.0, 1.0, 1.0])   # four small errors
errors_spiked   = np.array([0.0, 0.0, 0.0, 2.0])   # one larger error
for name, p in [("L1", 1), ("L2", 2)]:
    s, k = np.linalg.norm(errors_spread, p), np.linalg.norm(errors_spiked, p)
    print(f"{name}: spread={s:.2f} spiked={k:.2f} -> "
          f"{'spread' if s > k else 'spiked'} is penalized more")
# L1: spread=4.00 spiked=2.00 -> spread is penalized more
# L2: spread=2.00 spiked=2.00 -> a tie; push the spike higher and L2 loses

# Cosine distance ignores magnitude entirely, which is why it dominates
# text and embedding retrieval where document length should not matter.
short_doc = np.array([1.0, 2.0, 1.0])
long_doc  = np.array([10.0, 20.0, 10.0])      # same topic mix, 10x longer
print(np.linalg.norm(short_doc - long_doc))   # 22.05 -- "very different"
print(distance.cosine(short_doc, long_doc))   # 0.0  -- "identical direction"

Metric	Formula sketch	Robust to outliers	Typical use
L1 / Manhattan	sum of `abs(xi - yi)`	Yes	Lasso, robust loss, grids
L2 / Euclidean	sqrt of sum of squares	No	Default distance, ridge, MSE
L-infinity	max of `abs(xi - yi)`	No (max-driven)	Worst-case bounds
Cosine	1 minus normalized dot	Magnitude-blind	Text, embeddings, retrieval
Hamming	count of differing positions	N/A	Categorical, binary strings
Mahalanobis	L2 after decorrelating	No	Correlated features, outliers

Key Takeaways

Distance is a modelling choice, not a fixed fact, and it changes which points count as similar.

L2 squares errors and so punishes single large errors; L1 is proportional and therefore robust.

The L1 ball's corners lie on the axes, which is the geometric reason Lasso produces exact zeros.

Cosine distance discards magnitude, which is desirable for text and embeddings and wrong for physical quantities.

🧪 Practice

Compute L1, L2, and cosine distance between two feature vectors and explain why the rankings differ.

Show empirically that scaling one feature by 100 changes k-nearest neighbours under L2 but not under cosine.

Interview: You cluster customers by spend and visit count and every cluster is driven by spend alone. Why, and what do you change? (Hint: compare the numeric ranges of the two features.)

Eigenvalues and Eigenvectors

Most vectors, when hit with a matrix, change direction. A few special ones do not: they come out pointing the same way, merely stretched or shrunk. Those are the eigenvectors, and the factor by which each is scaled is its eigenvalue. The defining equation is A v = lambda v.

Why care? Because eigenvectors reveal the intrinsic axes of a transformation. Instead of describing what a matrix does to every possible input, you can describe what it does along a handful of special directions, and everything else follows. This is the engine behind principal component analysis, spectral clustering, PageRank, and stability analysis.

For a covariance matrix specifically, the eigenvectors are the directions of greatest variance in the data, and the eigenvalues are how much variance lies along each. That is precisely the statement that PCA rests on.

An analogy: spinning a slightly lopsided object, most ways you spin it, it wobbles. There are a few axes about which it spins cleanly. Those axes are the eigenvectors, intrinsic to the object rather than to how you happened to set up coordinates.

  A v = lambda v  -- direction preserved, length scaled

   generic vector u                eigenvector v
   +                               +
   | u                             | v
   |                               |
   +------->            A          +------->
                       ===>
   A u points somewhere else       A v = 2v, same line, twice as long
   +                               +
   |        \  Au                  |
   |         \                     |          Av
   +          v                    +------------------->

  EIGENVECTORS OF A COVARIANCE MATRIX = PRINCIPAL COMPONENTS

        x2
         ^        . ' .  . '
         |    . ' . '. ' .        <- first eigenvector: direction of
         |  ' . ' . ' .              maximum variance (long axis)
         | . ' . '
         +--------------------> x1
                                  <- second eigenvector: perpendicular,
                                     the remaining variance (short axis)

  eigenvalue 1 >> eigenvalue 2 means the cloud is nearly one-dimensional,
  so one component captures almost everything.

import numpy as np

A = np.array([[2.0, 1.0],
              [1.0, 2.0]])

values, vectors = np.linalg.eig(A)
print(values)          # [3. 1.]
print(vectors)         # columns are the eigenvectors, unit length
# column 0 ~ [0.707, 0.707]  -> the (1,1) direction, scaled by 3
# column 1 ~ [-0.707, 0.707] -> the (-1,1) direction, scaled by 1

v = vectors[:, 0]
print(np.allclose(A @ v, values[0] * v))    # True -- the defining property

# PCA in five lines, to show the connection is not metaphorical.
rng = np.random.default_rng(0)
data = rng.normal(size=(500, 2)) @ np.array([[3.0, 1.0],
                                             [0.0, 0.5]])   # correlated
centered = data - data.mean(axis=0)          # eigen-analysis assumes centred
cov = np.cov(centered, rowvar=False)         # (2,2) covariance matrix
vals, vecs = np.linalg.eigh(cov)             # eigh: symmetric, sorted, stable
order = np.argsort(vals)[::-1]               # largest variance first
vals, vecs = vals[order], vecs[:, order]

print(f"variance explained: {vals / vals.sum()}")   # e.g. [0.93 0.07]
# 93% of the spread lies along one direction, so the second axis is
# nearly discardable. That sentence IS dimensionality reduction.

Key Takeaways

Eigenvectors are the directions a transformation leaves pointing the same way, scaled by their eigenvalues.

They expose a matrix's intrinsic axes, which is why so many methods reduce to an eigen-decomposition.

For a covariance matrix, eigenvectors are directions of variance and eigenvalues are the amount along each.

Use eigh for symmetric matrices such as covariance; it is faster and numerically better than eig.

🧪 Practice

Compute eigenvalues and eigenvectors of [[4,0],[0,1]] by hand and explain geometrically what the matrix does.

Generate correlated 2-D data, find the covariance eigenvectors, and plot them over the scatter to confirm they follow the cloud's axes.

Interview: A covariance matrix has one eigenvalue near zero. What does that say about the features? (Hint: near-zero variance along a direction means the data barely moves that way.)

Singular Value Decomposition

Eigen-decomposition has a serious limitation: it requires a square matrix, and works cleanly only for some of those. Data matrices are almost never square. The singular value decomposition is the generalization that always exists, for any matrix of any shape, real or complex.

SVD factors any matrix A into three pieces: A = U S V'. Here V' describes a rotation in the input space, S is a diagonal matrix of non-negative singular values that stretch along each axis, and U is a rotation in the output space. The claim is remarkable and true: every linear transformation, however complicated, is a rotation, then a scaling, then another rotation.

The singular values are ordered from largest to smallest, and they measure how much of the matrix's action happens along each direction. This ordering is what makes SVD the workhorse of compression and noise reduction: keep the largest few, discard the rest, and you have the best possible low-rank approximation of the original.

An analogy: SVD is a spectral analysis for matrices. Just as a chord decomposes into notes ranked by loudness, a matrix decomposes into rank-one pieces ranked by importance. Dropping the quietest notes is barely audible.

  A  =  U  @  S  @  V'
 (m,n)  (m,m) (m,n) (n,n)        full SVD

  ANY TRANSFORMATION = ROTATE, SCALE, ROTATE

   input circle      V' rotates      S scales         U rotates
       _--_             _--_          _------_          /\
      /    \    ->     /    \   ->   |        |   ->   |  |
      \ _ _/           \ _ _/         \______/          \/

  TRUNCATION = COMPRESSION

  A  ~=  U[:, :k] @ S[:k, :k] @ V'[:k, :]
         \_____________________________/
          best rank-k approximation of A (Eckart-Young theorem)

  storage: m*n   ->   k*(m + n + 1)
  for a 1000x1000 image at k=50: 1,000,000 -> 100,050 numbers

import numpy as np

rng = np.random.default_rng(0)
A = rng.normal(size=(200, 50))

U, s, Vt = np.linalg.svd(A, full_matrices=False)   # economy SVD
print(U.shape, s.shape, Vt.shape)    # (200,50) (50,) (50,50)
print(np.allclose(A, U @ np.diag(s) @ Vt))         # True -- exact

# s is a 1-D array of singular values, already sorted descending.
print(s[:5].round(2))

# Low-rank approximation: keep the k largest singular values.
def truncate(U, s, Vt, k):
    return U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]

for k in (1, 5, 20, 50):
    approx = truncate(U, s, Vt, k)
    # Relative error in the Frobenius norm, i.e. total squared difference.
    err = np.linalg.norm(A - approx) / np.linalg.norm(A)
    kept = s[:k].sum() / s.sum()
    print(f"k={k:2d}  relative error={err:.3f}  singular mass kept={kept:.3f}")

# Energy captured is often the criterion for choosing k:
energy = np.cumsum(s**2) / np.sum(s**2)
k90 = int(np.searchsorted(energy, 0.90) + 1)
print(f"rank needed for 90% of squared variation: {k90}")

# Connection to PCA: on centred data, the right singular vectors (rows of
# Vt) ARE the principal components, and s**2/(n-1) are the eigenvalues
# of the covariance matrix. PCA is usually computed via SVD for stability.

Key Takeaways

SVD exists for every matrix, unlike eigen-decomposition, which needs a square and well-behaved one.

Every linear map is a rotation, then an axis-aligned scaling, then another rotation.

Truncating to the largest k singular values gives the provably best rank-k approximation.

PCA is normally computed through SVD of the centred data rather than by eigen-decomposing the covariance matrix.

🧪 Practice

Run SVD on a grayscale image and view reconstructions at k = 5, 20, and 100. Note where quality stops improving visibly.

Plot the cumulative energy curve for a dataset and choose k for 95%.

Interview: Why is PCA usually implemented with SVD rather than by eigen-decomposing X'X? (Hint: think about what squaring does to the condition number.)

Matrix Factorization Basics

The recurring idea behind the previous two topics deserves naming on its own. Matrix factorization means writing one matrix as a product of simpler ones. This is not an aesthetic preference. Factorizations expose structure, make computation cheaper, and let you fill in missing entries.

The last use is the one that built the recommender systems industry. Consider a matrix of users by items where the entry is a rating. It is almost entirely missing, since no user rates most items. Factorizing it as R ~= P Q', with P a small matrix of user factors and Q a small matrix of item factors, forces every rating to be explained by a handful of latent traits. Once fitted, the product gives a value for every cell, including the ones that were never observed. Those are the predictions.

An analogy: factoring 143 into 11 times 13 does not change the number but reveals structure invisible in the original. Matrix factorization does the same for data, except the factors have interpretable meaning: latent genres, topics, or customer archetypes.

  RATINGS MATRIX             USER FACTORS      ITEM FACTORS
  (mostly missing)

        i1  i2  i3  i4              f1   f2         i1   i2   i3   i4
   u1 [ 5   ?   3   ? ]        u1 [ 1.2 -0.3 ]  f1 [ 2.1  0.4  1.8  0.2 ]
   u2 [ ?   4   ?   1 ]   ~=   u2 [ 0.1  1.5 ] @ f2 [ 0.3  2.2  0.1  0.9 ]
   u3 [ 4   ?   4   ? ]        u3 [ 1.1  0.0 ]
   u4 [ ?   5   ?   2 ]        u4 [-0.2  1.7 ]

     (4 x 4), 8 known            (4 x 2)            (2 x 4)
     8 unknown                   \________ 16 parameters ________/

  The product fills EVERY cell. Cell (u1, i2) becomes a prediction.
  The low rank is the point: it forces generalization instead of
  memorization, exactly as a small model does.

  COMMON FACTORIZATIONS
   LU    A = LU        solving systems, determinants
   QR    A = QR        least squares, stable regression
   SVD   A = U S V'    compression, PCA, pseudo-inverse
   NMF   A ~= W H      parts-based, non-negative, topics
   Chol  A = L L'      symmetric positive definite, sampling, speed

import numpy as np
from sklearn.decomposition import NMF, TruncatedSVD

# NMF constrains both factors to be non-negative, which makes the
# factors additive parts rather than signed mixtures -- much easier
# to interpret as "topics" or "genres".
counts = np.array([                 # documents x words
    [3, 0, 1, 0, 0],
    [2, 0, 2, 0, 1],
    [0, 4, 0, 3, 0],
    [0, 3, 0, 4, 1],
], dtype=float)

nmf = NMF(n_components=2, init="nndsvd", random_state=0, max_iter=500)
W = nmf.fit_transform(counts)       # documents x topics
H = nmf.components_                 # topics x words

print("document-topic weights:\n", W.round(2))
print("topic-word weights:\n", H.round(2))
print("reconstruction:\n", (W @ H).round(2))
# Docs 1-2 load on one topic, docs 3-4 on the other, and H shows which
# words define each. No labels were supplied; the structure was latent.

# Truncated SVD does the same job without the non-negativity constraint,
# and handles sparse matrices, which is why it is standard for text.
svd = TruncatedSVD(n_components=2, random_state=0)
embedded = svd.fit_transform(counts)
print("explained variance:", svd.explained_variance_ratio_.round(3))

Key Takeaways

Factorization rewrites a matrix as a product of simpler matrices, exposing structure and reducing cost.

Low rank forces generalization, which is what lets a factorization predict entries it never observed.

Different factorizations serve different goals: QR for least squares, SVD for compression, NMF for interpretable parts.

Non-negativity constraints make factors readable as additive components rather than signed mixtures.

🧪 Practice

Build a small ratings matrix with missing entries and use NMF to predict two held-out ratings.

Compare NMF and truncated SVD factors on the same document-term matrix and describe which is easier to interpret.

Interview: In a recommender, what does increasing the number of latent factors do to training error and to test error? (Hint: the rank is a capacity parameter like any other.)

2.2 Calculus and Optimization

Fitting a model means searching for parameters that minimize a loss, and calculus provides both the compass that says which way to move and the guarantees that say when moving will work.

Derivatives and Partial Derivatives

Training a model is a search problem: adjust the parameters until the error is as small as possible. To search efficiently you need to know, from where you currently stand, whether nudging a parameter up or down makes the error better or worse. That is exactly what a derivative reports.

A derivative measures the instantaneous rate of change of a function. If f(x) is the error and x is a weight, then f'(x) tells you how fast the error changes per unit change in the weight, right here. Its sign says which direction is uphill; its magnitude says how steep.

Real models have many parameters, so the function has many inputs. A partial derivative is the derivative with respect to one input while holding all others fixed. Written df/dx, it answers: if I change only this weight and touch nothing else, how does the loss respond?

An analogy: standing on a hillside, the derivative in the east-west direction is how steeply the ground rises if you step east. The north-south partial is a different number. Neither alone tells you the steepest way up; together they do, which is the next topic.

  DERIVATIVE = SLOPE OF THE TANGENT

   f(x)
    ^                        .
    |                      .'
    |                    .'  <- steep: large derivative, big step helps
    |                 ..'
    |            ...''
    |      ..''''            <- flat: derivative near 0, you are near
    |  .''                      a minimum (or a plateau)
    +-------------------------> x
         x0        x1

  PARTIAL DERIVATIVES ON A SURFACE

   loss(w1, w2), a bowl seen from above as contour lines

        w2
         ^   ____
         |  /    \        df/dw1 : slope moving purely right
         | |  .   |       df/dw2 : slope moving purely up
         |  \____/        At the centre both are 0 -> minimum
         +-----------> w1

import numpy as np

# Mean squared error for a one-parameter model: f(w) = mean((w*x - y)^2)
x = np.array([1.0, 2.0, 3.0, 4.0])
y = np.array([2.0, 4.0, 6.0, 8.0])        # true relationship: y = 2x

def loss(w):
    return np.mean((w * x - y) ** 2)

def dloss(w):
    # d/dw mean((wx - y)^2) = mean(2x(wx - y))  -- chain rule applied
    return np.mean(2 * x * (w * x - y))

for w in [0.0, 1.0, 2.0, 3.0]:
    print(f"w={w}: loss={loss(w):6.2f}  derivative={dloss(w):7.2f}")
# w=0.0: loss= 30.00  derivative= -30.00   <- negative: increase w
# w=1.0: loss=  7.50  derivative= -15.00   <- still increase, less urgently
# w=2.0: loss=  0.00  derivative=   0.00   <- zero: this is the minimum
# w=3.0: loss=  7.50  derivative=  15.00   <- positive: decrease w

# Numerical check by finite differences. Useful for verifying hand-derived
# gradients, which is where most training bugs actually live.
def numeric_derivative(f, w, h=1e-6):
    return (f(w + h) - f(w - h)) / (2 * h)      # central difference

print(np.isclose(dloss(1.5), numeric_derivative(loss, 1.5)))   # True

Key Takeaways

A derivative gives the local rate of change: its sign is the direction and its magnitude is the steepness.

A partial derivative varies one parameter while holding the rest fixed.

A derivative of zero marks a flat point, which is a candidate minimum, maximum, or saddle.

Verify analytic derivatives against finite differences; mismatched gradients are a leading cause of silent training failure.

🧪 Practice

Differentiate f(w) = (w - 3)^2 + 5 by hand, find where the derivative is zero, and confirm it is the minimum.

Write a gradient-check function and use it to catch a deliberately wrong analytic gradient.

Interview: Your loss stops decreasing and the gradient is near zero, but the loss is still high. What are two possible explanations? (Hint: zero gradient does not by itself mean a good minimum.)

Gradients and the Jacobian

One partial derivative tells you about one direction. With a hundred parameters you have a hundred partials, and you need them assembled into something you can act on. That assembly is the gradient.

The gradient, written grad f, is the vector of all partial derivatives of a scalar-valued function. Its crucial property is geometric: it points in the direction of steepest increase, and its length is how steep that is. This single fact is the basis of essentially all model training. To decrease the loss, step in the direction of the negative gradient.

When the output is a vector rather than a scalar, the generalization is the Jacobian: a matrix whose row i is the gradient of output i. Neural network layers are vector-to-vector functions, so their local behaviour is described by Jacobians, and backpropagation is a sequence of Jacobian products arranged to avoid ever forming them explicitly.

An analogy: on a hillside in fog, the gradient is what a marble tells you. Set it down and it rolls along the negative gradient, the steepest way down, without needing a map of the whole terrain.

  THE GRADIENT POINTS STEEPEST UPHILL

   contour map of loss(w1, w2)

     w2 ^
        |    .-''''''-.
        |   /   .--.   \        grad f  ->  perpendicular to the contour,
        |  |   / .. \   |                   pointing toward higher loss
        |  |  | (  ) |  |
        |  |   \ '' /   |       -grad f ->  the descent direction we take
        |   \   '--'   /
        |    '-......-'
        +------------------> w1

  Key property: the gradient is always perpendicular to the contour line
  it sits on. Moving along a contour changes nothing; moving across it
  changes the loss fastest.

  SHAPES

   scalar output:  f: R^n -> R      grad is (n,)        one gradient
   vector output:  f: R^n -> R^m    Jacobian is (m, n)  m gradients stacked

import numpy as np

# f(w1, w2) = w1^2 + 3*w2^2  -- an elongated bowl
def f(w):
    return w[0]**2 + 3 * w[1]**2

def grad_f(w):
    # partials: df/dw1 = 2*w1,  df/dw2 = 6*w2
    return np.array([2 * w[0], 6 * w[1]])

w = np.array([3.0, 1.0])
g = grad_f(w)
print(g)                         # [6. 6.] -- equal partials here
print(np.linalg.norm(g))         # 8.49 -- steepness at this point

# The descent step goes AGAINST the gradient.
step = 0.1
print(f(w), f(w - step * g))     # 12.0 -> 6.24, the loss fell

# Jacobian for a vector-valued function g: R^2 -> R^3
def vec_fn(w):
    return np.array([w[0] * w[1], w[0]**2, np.sin(w[1])])

def jacobian(w):
    # row i = gradient of output i  -> shape (3, 2)
    return np.array([
        [w[1],        w[0]      ],   # d(w1*w2)/dw1, d(w1*w2)/dw2
        [2 * w[0],    0.0       ],   # d(w1^2)/dw1,  d(w1^2)/dw2
        [0.0,         np.cos(w[1])], # d(sin w2)/dw1, d(sin w2)/dw2
    ])

print(jacobian(np.array([2.0, 0.5])).shape)   # (3, 2)

# Note the anisotropy in f: the w2 direction is 3x steeper than w1.
# A single shared learning rate must be small enough for the steep
# direction, so it crawls along the shallow one. That mismatch is the
# reason adaptive optimizers such as Adam exist.

Key Takeaways

The gradient collects all partial derivatives and points in the direction of steepest increase.

Training descends by stepping along the negative gradient.

The gradient is perpendicular to the contour line through the current point.

The Jacobian generalizes the gradient to vector-valued functions and underlies backpropagation.

🧪 Practice

Compute the gradient of f(w1,w2) = w1^2 + w1*w2 + w2^2 at (1, 2) by hand and verify numerically.

Write the Jacobian of the softmax function for a 3-element input and check its shape.

Interview: Why does a loss surface with very different curvature along two axes slow down plain gradient descent? (Hint: one learning rate must serve both directions.)

Chain Rule

Models are compositions. A neural network applies a linear map, then an activation, then another linear map, then a loss. To train it you need the derivative of the final loss with respect to a weight buried several layers back. The chain rule is what makes that computable, and it is the single most important calculus result in machine learning.

The rule states that if y = f(g(x)), then dy/dx = f'(g(x)) * g'(x). Derivatives of composed functions multiply. Sensitivity propagates: if a small change in x moves g a lot, and a small change in g moves f a lot, then x moves f by the product.

Backpropagation is the chain rule applied to a computation graph, evaluated from the output backwards. Going backwards is not arbitrary: with one scalar loss and many parameters, the reverse direction reuses shared work and costs roughly the same as one forward pass, while the forward direction would cost one pass per parameter.

An analogy: a chain of gears. If gear A turns gear B three times as fast, and B turns C twice as fast, then A turns C six times as fast. Ratios multiply along the chain.

  FORWARD PASS: values flow right

   x ---[ g ]---> h ---[ f ]---> y ---[ loss ]---> L

  BACKWARD PASS: derivatives flow left, multiplying as they go

   dL/dx  <--- dL/dh <--- dL/dy <--- dL/dL = 1
    |          |           |
    |          |           +-- local derivative of loss w.r.t. y
    |          +-- multiply by f'(h)
    +-- multiply by g'(x)

   dL/dx = dL/dy * dy/dh * dh/dx
           \_______ a product of local derivatives _______/

  WHY BACKWARD IS CHEAPER
   forward mode:  cost proportional to NUMBER OF INPUTS
   reverse mode:  cost proportional to NUMBER OF OUTPUTS
   deep learning: millions of inputs (weights), one output (loss)
   => reverse mode wins by a factor of millions

import numpy as np

# A two-layer chain: z = w2 * tanh(w1 * x), loss = (z - target)^2
x, target = 2.0, 1.0
w1, w2 = 0.5, 1.5

# ---- forward pass: save intermediate values, they are needed later ----
a = w1 * x                 # pre-activation
h = np.tanh(a)             # hidden activation
z = w2 * h                 # output
L = (z - target) ** 2      # loss
print(f"forward: a={a:.4f} h={h:.4f} z={z:.4f} L={L:.4f}")

# ---- backward pass: multiply local derivatives, right to left ----
dL_dz = 2 * (z - target)          # d/dz (z - t)^2
dL_dw2 = dL_dz * h                # z = w2*h, so dz/dw2 = h
dL_dh = dL_dz * w2                # dz/dh = w2
dL_da = dL_dh * (1 - h ** 2)      # d/da tanh(a) = 1 - tanh(a)^2
dL_dw1 = dL_da * x                # a = w1*x, so da/dw1 = x

print(f"dL/dw1 = {dL_dw1:.6f}   dL/dw2 = {dL_dw2:.6f}")

# Verify against finite differences -- always do this when hand-deriving.
def loss_at(w1_, w2_):
    return (w2_ * np.tanh(w1_ * x) - target) ** 2

eps = 1e-6
num_w1 = (loss_at(w1 + eps, w2) - loss_at(w1 - eps, w2)) / (2 * eps)
print(f"numeric dL/dw1 = {num_w1:.6f}")     # matches to ~6 decimals

# Note dL_da multiplies by (1 - h^2), which is at most 1 and near 0 when
# |h| is close to 1. Stack fifty such layers and the product collapses.
# That is the vanishing gradient problem, and it is pure chain rule.

Key Takeaways

Derivatives of composed functions multiply, which is what makes deep models trainable.

Backpropagation is the chain rule evaluated backwards over a computation graph.

Reverse mode costs scale with the number of outputs, so one scalar loss with millions of weights is the ideal case.

Repeated multiplication by small local derivatives causes vanishing gradients, and by large ones, exploding gradients.

🧪 Practice

Differentiate f(x) = sin(3x^2 + 1) by hand, naming each local derivative in the chain.

Extend the code above to three layers and verify every gradient numerically.

Interview: Why do sigmoid activations cause vanishing gradients in deep networks? (Hint: find the maximum value of the sigmoid's derivative and raise it to the power of the depth.)

Convexity

Gradient descent follows the slope downhill. Whether that strategy actually finds the best answer, or merely a locally acceptable one, depends on the shape of the loss surface. Convexity is the property that makes the difference, and knowing whether your problem has it tells you how much to trust your optimizer.

A function is convex if the line segment joining any two points on its graph lies on or above the graph. Visually, it is bowl-shaped with no dents. The consequence is decisive: a convex function has no local minima other than the global one. Any downhill path reaches the best answer.

Non-convex functions have multiple valleys, plateaus, and saddle points. Gradient descent can settle in a poor valley, and the result depends on where you started. Neural network losses are thoroughly non-convex, which is why random initialization and restarts matter there and not in linear regression.

A practical test in one dimension: a twice-differentiable function is convex when its second derivative is non-negative everywhere. In many dimensions, the equivalent condition is that the Hessian matrix is positive semi-definite, meaning all its eigenvalues are non-negative.

  CONVEX                          NON-CONVEX

   \                    /          \      /\        /
    \                  /            \    /  \      /
     \                /              \  /    \    /
      \              /                \/      \  /
       \____________/                          \/
                                   local min   global min
   any chord lies ABOVE           a chord can cut BELOW the curve
   one minimum, always found      where you start decides where you end

  THE CHORD TEST

   f(x)                             f(x)
    |    .                           |   .        .
    |   . \___  chord above          |  . \  .-'/   chord dips below
    |  .      '-.                    | .   \/  /    the curve here
    | .          '.                  |.    /\ /
    +----------------                +----------------

Property	Convex problem	Non-convex problem
Local minima	All are global	Can be far from global
Effect of initialization	None on final answer	Decides which basin you reach
Stopping guarantee	Zero gradient is optimal	Zero gradient means little
Typical models	Linear, logistic, SVM, Lasso	Neural nets, k-means, GMM
Practical tactic	Just run the optimizer	Restarts, schedules, momentum

import numpy as np

# Second-derivative test in one dimension.
def is_convex_1d(second_derivative, grid):
    return np.all(second_derivative(grid) >= -1e-12)   # tolerance for float

grid = np.linspace(-5, 5, 1001)
print(is_convex_1d(lambda x: np.full_like(x, 2.0), grid))        # x^2  -> True
print(is_convex_1d(lambda x: 12 * x**2 - 12, grid))              # x^4-6x^2 -> False

# Multivariate test: all Hessian eigenvalues non-negative.
def hessian_convex(H):
    return np.all(np.linalg.eigvalsh(H) >= -1e-12)

print(hessian_convex(np.array([[2.0, 0.0], [0.0, 6.0]])))    # bowl -> True
print(hessian_convex(np.array([[2.0, 0.0], [0.0, -6.0]])))   # saddle -> False

# Why it matters in practice: the same optimizer, two problems.
def descend(grad, start, lr=0.05, steps=300):
    x = float(start)
    for _ in range(steps):
        x -= lr * grad(x)
    return x

convex_grad     = lambda x: 2 * x                     # from x^2
nonconvex_grad  = lambda x: 4 * x**3 - 12 * x         # from x^4 - 6x^2

print([round(descend(convex_grad, s), 4) for s in (-4, 0.5, 7)])
# [0.0, 0.0, 0.0]        -- same answer from every start
print([round(descend(nonconvex_grad, s), 4) for s in (-4, 0.5, 7)])
# [-1.7321, 1.7321, 1.7321]  -- the starting point decided the outcome

Key Takeaways

A convex function is bowl-shaped, so every local minimum is the global minimum.

On convex problems the starting point is irrelevant; on non-convex ones it decides the basin you land in.

Convexity is checked by a non-negative second derivative, or a positive semi-definite Hessian in higher dimensions.

Linear and logistic regression are convex; neural networks and k-means are not, which is why the latter need restarts.

🧪 Practice

Show that f(x) = x^4 - 3x^2 is non-convex and find both local minima.

Verify that the logistic regression loss is convex in the weights for a small dataset by checking the Hessian numerically.

Interview: k-means gives different clusters on different runs. Is that a bug? (Hint: consider the shape of the objective it minimizes.)

Gradient Descent

With a gradient available and convexity understood, the algorithm follows directly: to reduce the loss, repeatedly step in the direction of the negative gradient. Gradient descent is the most important optimization algorithm in machine learning, not because it is the best, but because it is the only one that scales to millions of parameters.

The update rule is w := w - lr * grad(w), where lr is the learning rate. That one hyperparameter controls everything about whether training works. Too small and progress is glacial. Too large and steps overshoot the valley, bouncing up the far wall and diverging. There is a critical value above which divergence is guaranteed, determined by the curvature.

Why not solve for the minimum directly? For linear regression you can, via the normal equations. But that requires inverting a matrix whose size is the number of features, costing roughly cubic time, and it does not exist for most models. Gradient descent needs only the ability to compute a gradient, which is why it generalizes to everything.

An analogy: descending a mountain in fog by always stepping downhill. Learning rate is stride length. Tiny steps take forever; huge steps step over the valley and land higher up the opposite slope.

  THE LEARNING RATE DECIDES EVERYTHING

  lr too small          lr well chosen        lr too large
   \        /            \        /            \        /
    \  ....            .  \      /              \  1   /
     \ ....           2 \  \    /                \ / \/ 3
      \...           .   \1 \  /                  X   X
       \..          3 \   \  \/                  / \ / \
        \_/            \___\_/                  /   V   2
                                               /    4    \
  crawls, may stop      converges              diverges: each step
  before converging     in a few steps         overshoots further

  CONVERGENCE CRITERION
   stop when any of:
     - gradient norm < tolerance          (flat: likely a minimum)
     - loss change between steps < tol    (no progress)
     - max iterations reached             (give up, investigate)

import numpy as np

# Linear regression by gradient descent, with the closed form for comparison.
rng = np.random.default_rng(0)
X = np.c_[np.ones(200), rng.normal(size=(200, 2))]      # bias column + 2 feats
true_w = np.array([1.0, 2.0, -3.0])
y = X @ true_w + rng.normal(scale=0.5, size=200)

def mse_grad(w):
    residual = X @ w - y
    # d/dw mean(residual^2) = (2/n) X' residual
    return (2 / len(y)) * X.T @ residual

def gradient_descent(lr, steps=500, tol=1e-8):
    w = np.zeros(X.shape[1])
    history = []
    for i in range(steps):
        g = mse_grad(w)
        w = w - lr * g                       # THE update rule
        history.append(np.mean((X @ w - y) ** 2))
        if np.linalg.norm(g) < tol:          # converged: gradient is flat
            break
    return w, history

for lr in (0.001, 0.1, 1.05):
    w, hist = gradient_descent(lr)
    status = "diverged" if not np.isfinite(hist[-1]) or hist[-1] > 1e3 else "ok"
    print(f"lr={lr:<6} final loss={hist[-1]:>12.4f}  {status}")
# lr=0.001  loss still falling   ok        (far too slow, not converged)
# lr=0.1    loss near the noise floor of 0.25  ok   (converged)
# lr=1.05   loss overflows to inf  diverged  (overshoots every step)
#
# The divergence threshold is not arbitrary: gradient descent on a
# quadratic diverges once lr exceeds 2 / (largest Hessian eigenvalue).
# Here the Hessian is (2/n) X'X, whose largest eigenvalue is near 2,
# so the critical learning rate is near 1.0.

# Closed form, available for linear regression only:
w_exact = np.linalg.solve(X.T @ X, X.T @ y)   # normal equations
print("exact:   ", w_exact.round(3))
print("descent: ", gradient_descent(0.1, 5000)[0].round(3))   # same answer

# Feature scaling is not optional. Unscaled features create elongated
# contours, and the largest safe learning rate is set by the steepest
# direction, so the shallow ones barely move.

Key Takeaways

The update w := w - lr * grad(w) is the core of nearly all model training.

The learning rate is the decisive hyperparameter: too small stalls, too large diverges.

Gradient descent needs only a gradient, which is why it works where closed-form solutions do not exist.

Unscaled features distort the loss surface and force a smaller usable learning rate.

🧪 Practice

Implement gradient descent for f(w) = (w - 4)^2 and plot the loss for three learning rates including one that diverges.

Fit a linear regression by gradient descent and confirm it matches the normal equations to three decimals.

Interview: Your loss decreases then suddenly becomes NaN. What are the likely causes? (Hint: think about the step size and the magnitude of the inputs.)

Stochastic and Mini-Batch Gradient Descent

Plain gradient descent computes the gradient over the entire dataset before taking a single step. With ten million rows, that is ten million evaluations to move once. Training would take days to perform a few hundred updates. The fix is to accept a noisy estimate of the gradient in exchange for taking far more steps.

Stochastic gradient descent uses one randomly chosen example per update. The gradient from one example is a poor estimate of the true gradient, but it is unbiased, and cheap enough that you can take thousands of steps in the time full-batch takes one. Mini-batch gradient descent is the compromise used in practice: estimate the gradient on a small batch, typically 32 to 512 examples, which reduces noise and maps efficiently onto hardware that performs matrix operations in parallel.

The noise turns out to be a feature, not just a tolerated cost. On non-convex surfaces, the jitter helps escape shallow local minima and saddle points that would trap a deterministic optimizer.

An analogy: full-batch is surveying the whole mountain before each step. Stochastic is glancing at the ground and stepping immediately. You wobble, but you cover far more ground per hour, and the wobble sometimes shakes you out of a ditch.

  PATHS TO THE MINIMUM

   FULL BATCH              MINI-BATCH            STOCHASTIC
    smooth, few steps       slight wobble         very noisy

      \                       \                     \  ,
       \                       \.                    '\ '
        \                       \'                   ,\'
         \                       \.                 '  \,
          *                       *                    '*

   1 step per epoch        ~n/B steps per epoch   n steps per epoch
   exact gradient          low-variance estimate  high-variance estimate
   slow on big data        the practical default  rarely used alone

  COST PER EPOCH IS ROUGHLY EQUAL; NUMBER OF UPDATES IS NOT.
   n = 1,000,000 rows, batch size B = 256
   full batch:  1 update per epoch
   mini-batch:  3,906 updates per epoch

import numpy as np

rng = np.random.default_rng(0)
n = 10_000
X = np.c_[np.ones(n), rng.normal(size=(n, 3))]
true_w = np.array([0.5, 2.0, -1.0, 3.0])
y = X @ true_w + rng.normal(scale=0.5, size=n)

def batch_grad(w, Xb, yb):
    return (2 / len(yb)) * Xb.T @ (Xb @ w - yb)

def train(batch_size, epochs=10, lr=0.05, seed=0):
    rng_local = np.random.default_rng(seed)
    w = np.zeros(X.shape[1])
    updates = 0
    for _ in range(epochs):
        order = rng_local.permutation(n)        # reshuffle every epoch:
        for start in range(0, n, batch_size):   # without it, batch order
            idx = order[start:start + batch_size]   # becomes a fixed bias
            w -= lr * batch_grad(w, X[idx], y[idx])
            updates += 1
    return w, updates

for bs in (n, 256, 1):
    w, updates = train(bs)
    label = {n: "full batch", 256: "mini-batch", 1: "stochastic"}[bs]
    err = np.linalg.norm(w - true_w)
    print(f"{label:<12} batch={bs:<6} updates={updates:<6} |w-w*|={err:.4f}")
# full batch   batch=10000  updates=10     |w-w*|~1.3     barely moved
# mini-batch   batch=256    updates=400    |w-w*|~0.01    converged
# stochastic   batch=1      updates=100000 |w-w*|~0.02    converged, noisy
#
# Full batch is not worse per update; it simply took 10 updates while
# mini-batch took 400 for the same number of passes over the data.

Variant	Batch size	Updates per epoch	Gradient quality	Notes
Full batch	All rows	1	Exact	Smooth, impractical at scale
Mini-batch	32 to 512	n / B	Low variance	The practical default
Stochastic	1	n	High variance	Noisy, poor hardware use

Key Takeaways

Mini-batching trades gradient accuracy for far more updates in the same compute budget.

The batch gradient is an unbiased estimate of the full gradient, so the noise averages out over steps.

Gradient noise helps escape saddle points and shallow minima on non-convex surfaces.

Reshuffle the data each epoch, or the fixed batch composition becomes a systematic bias.

🧪 Practice

Train the same model at batch sizes 1, 32, and full, and plot loss against wall-clock time rather than epochs.

Remove the per-epoch shuffle and describe how the loss curve changes.

Interview: Doubling the batch size, how should you change the learning rate and why? (Hint: consider what happens to the variance of the gradient estimate.)

Constrained Optimization and Lagrange Multipliers

So far the search has been unconstrained: any parameter values are allowed. Real problems often forbid that. Portfolio weights must sum to one. Probabilities must be non-negative. A regularized model must keep its coefficients within a budget. These are constraints, and they change the answer, because the best point overall may be outside the permitted region.

The key insight is geometric. At a constrained optimum, you cannot improve the objective without leaving the feasible set. For an equality constraint g(x) = 0, that happens exactly when the objective's gradient is parallel to the constraint's gradient: every direction that keeps you on the constraint surface is perpendicular to the objective's gradient, so no legal move improves anything.

Writing "parallel" as "equal up to a scalar" gives grad f = lambda grad g. That scalar is the Lagrange multiplier, and it has a useful interpretation of its own: it is the rate at which the optimal objective value would improve if the constraint were relaxed by one unit, which economists call the shadow price.

This is not an abstraction in data science. Ridge and Lasso regression are literally constrained least squares: minimize error subject to a budget on coefficient size. The regularization strength people tune is a Lagrange multiplier.

  TANGENCY AT THE CONSTRAINED OPTIMUM

        contours of f (objective)
             ___
           /     \            * unconstrained minimum (infeasible)
          | . ' . |
          |'  *  '|
           \ ___ /
        ______________  g(x) = 0   (the feasible line)
              @
              ^
              constrained optimum: the LOWEST contour that still
              touches the line. At that point the contour is TANGENT
              to the constraint, so their gradients are parallel:
                   grad f = lambda * grad g

  REGULARIZATION IS THIS PICTURE

   minimize ||y - Xw||^2   subject to   ||w||^2 <= t     (ridge)
   minimize ||y - Xw||^2   subject to   ||w||_1 <= t     (lasso)

   The L1 ball is a diamond with corners on the axes, so tangency
   commonly occurs AT a corner, where a coefficient is exactly zero.
   That is why Lasso selects features and Ridge does not.

import numpy as np
from scipy.optimize import minimize

# Minimize f(x,y) = x^2 + y^2 subject to x + y = 1.
# By symmetry the answer should be (0.5, 0.5).

f = lambda v: v[0]**2 + v[1]**2
constraint = {"type": "eq", "fun": lambda v: v[0] + v[1] - 1}   # g(v) = 0

result = minimize(f, x0=[2.0, 0.0], constraints=[constraint])
print(result.x.round(4))       # [0.5 0.5]
print(result.fun.round(4))     # 0.5

# Analytic solution via the Lagrangian:
#   L(x,y,lam) = x^2 + y^2 - lam*(x + y - 1)
#   dL/dx = 2x - lam = 0   -> x = lam/2
#   dL/dy = 2y - lam = 0   -> y = lam/2
#   dL/dlam = x + y - 1 = 0 -> lam = 1, so x = y = 0.5
# The multiplier lam = 1 says: relaxing the constraint to x+y = 1.01
# would change the optimal objective by about 0.01 * 1.

# Ridge regression as constrained least squares, two equivalent views.
from sklearn.linear_model import Ridge, Lasso
rng = np.random.default_rng(0)
X = rng.normal(size=(100, 8))
y = X @ np.array([3.0, -2.0, 0, 0, 0, 0, 0, 1.5]) + rng.normal(scale=.5, size=100)

ridge = Ridge(alpha=5.0).fit(X, y)     # alpha IS the Lagrange multiplier
lasso = Lasso(alpha=0.2).fit(X, y)
print("ridge:", ridge.coef_.round(2))  # all shrunk, none exactly zero
print("lasso:", lasso.coef_.round(2))  # several are exactly 0.0
# The difference is entirely the shape of the constraint region.

Key Takeaways

At a constrained optimum the objective's gradient is parallel to the constraint's gradient, which is the Lagrange condition.

The multiplier is a shadow price: the improvement gained per unit of constraint relaxation.

Ridge and Lasso are constrained least squares, and their regularization strength is a Lagrange multiplier.

The geometry of the constraint region decides the behaviour, which is why L1 zeroes coefficients and L2 only shrinks them.

🧪 Practice

Minimize x^2 + 4y^2 subject to x + y = 3 by hand using a Lagrangian, then verify with scipy.optimize.minimize.

Fit Ridge and Lasso across a range of alpha values and plot how each coefficient path behaves.

Interview: Why does Lasso produce exactly zero coefficients while Ridge does not? (Hint: draw both constraint regions and ask where a contour first touches each.)

2.3 Probability Theory

Probability is the mathematics of uncertainty, and since every dataset is a sample rather than the whole truth, it is the language in which every honest data science conclusion is stated.

Sample Spaces and Events

Before computing any probability you must be precise about what could happen. Vagueness here produces paradoxes that look profound and are merely sloppy. The formal setup exists to force that precision.

The sample space, written S or Omega, is the set of all possible outcomes of an experiment, where outcomes are mutually exclusive and collectively exhaustive. Exactly one occurs. An event is any subset of the sample space: a collection of outcomes you care about. "The die shows an even number" is the event {2, 4, 6}.

Because events are sets, set operations have probabilistic meaning. Union is "A or B", intersection is "A and B", complement is "not A". Three axioms, due to Kolmogorov, generate everything else: probabilities are non-negative, the whole sample space has probability 1, and probabilities of mutually exclusive events add.

The one rule worth memorizing is inclusion-exclusion, because forgetting it is the most common arithmetic error in the subject. P(A or B) is not P(A) + P(B) unless the events cannot co-occur; otherwise the overlap is counted twice and must be subtracted.

  SAMPLE SPACE AND EVENTS AS SETS

   S = all outcomes
   +-----------------------------------+
   |                                   |
   |    +---------+                    |
   |    |    A    |                    |     A          = event A
   |    |     +---+-----+              |     A and B    = overlap
   |    |     |///|     |              |     A or B     = both regions
   |    +-----+---+     |              |     not A      = everything else
   |          |    B    |              |
   |          +---------+              |
   |                                   |
   +-----------------------------------+

  P(A or B) = P(A) + P(B) - P(A and B)
                              ^
                              the shaded overlap, otherwise double-counted

  MUTUALLY EXCLUSIVE means the overlap is empty:
   +-------+   +-------+     P(A or B) = P(A) + P(B), no correction
   |   A   |   |   B   |
   +-------+   +-------+

from itertools import product

# Sample space for two fair dice: 36 equally likely outcomes.
S = list(product(range(1, 7), repeat=2))
print(len(S))                                  # 36

def prob(event):
    """Probability of an event given as a predicate over outcomes."""
    return sum(1 for outcome in S if event(outcome)) / len(S)

A = lambda o: sum(o) == 7           # the sum is seven
B = lambda o: o[0] == 6             # the first die shows six

print(f"P(A)        = {prob(A):.4f}")                  # 0.1667 = 6/36
print(f"P(B)        = {prob(B):.4f}")                  # 0.1667 = 6/36
print(f"P(A and B)  = {prob(lambda o: A(o) and B(o)):.4f}")   # 0.0278 = 1/36

# The naive sum double-counts the (6,1) outcome that satisfies both.
naive = prob(A) + prob(B)
correct = prob(lambda o: A(o) or B(o))
print(f"naive       = {naive:.4f}")     # 0.3333  WRONG
print(f"correct     = {correct:.4f}")   # 0.3056  = 11/36
print(f"difference  = {naive - correct:.4f}")   # 0.0278, exactly the overlap

Key Takeaways

A sample space lists mutually exclusive, collectively exhaustive outcomes, and exactly one occurs.

Events are subsets, so set operations translate directly into "and", "or", and "not".

P(A or B) = P(A) + P(B) - P(A and B); dropping the overlap term is the most common error.

Most probability paradoxes dissolve once the sample space is written down explicitly.

🧪 Practice

Write the sample space for three coin flips and compute the probability of at least two heads.

Verify inclusion-exclusion empirically for two overlapping events by enumerating a sample space in code.

Interview: A family has two children and at least one is a girl. What is the probability both are girls? (Hint: enumerate the sample space and remove only the outcomes the information excludes.)

Conditional Probability

Probabilities are not fixed properties of events; they depend on what you know. Learning that a patient is a smoker changes the probability of lung disease. Conditional probability is the machinery for updating in light of information, and it is the foundation of every predictive model, since a prediction is a conditional probability given features.

The definition is P(A | B) = P(A and B) / P(B), read as "the probability of A given B". The intuition is that learning B occurred shrinks the sample space to B alone. You then ask what fraction of that reduced world also has A. The division by P(B) renormalizes so probabilities in the new world still sum to one.

Rearranging gives the multiplication rule, P(A and B) = P(A | B) P(B), which is how joint probabilities are built up.

Two events are independent when conditioning changes nothing: P(A | B) = P(A), equivalently P(A and B) = P(A) P(B). Independence is an assumption that must be justified, not a default. Assuming it wrongly is one of the most consequential errors in applied statistics, because it makes compound risks look far smaller than they are.

  CONDITIONING SHRINKS THE WORLD

  BEFORE: full sample space S      AFTER learning B: only B remains
  +---------------------------+       +-------------+
  |          +--------+       |       |    +----+   |
  |   +------+--------+---+   |       |    |////|   |   the new universe
  |   |  A   |////////| B |   |  -->  |    |////|   |   is B; P(A|B) is
  |   +------+--------+---+   |       |    +----+   |   the shaded share
  |          +--------+       |       +-------------+
  +---------------------------+
                                      P(A|B) = shaded / total(B)
                                             = P(A and B) / P(B)

  INDEPENDENCE
   P(A|B) = P(A)  -- knowing B tells you nothing about A
   Careful: independent is NOT the same as mutually exclusive.
   Mutually exclusive events are maximally DEPENDENT: if one happens,
   the other definitely did not.

import pandas as pd
import numpy as np

rng = np.random.default_rng(0)
n = 10_000
smoker = rng.random(n) < 0.20
# Disease risk depends on smoking: 15% for smokers, 2% otherwise.
disease = rng.random(n) < np.where(smoker, 0.15, 0.02)
df = pd.DataFrame({"smoker": smoker, "disease": disease})

p_disease = df["disease"].mean()
p_smoker = df["smoker"].mean()
p_both = (df["smoker"] & df["disease"]).mean()

print(f"P(disease)            = {p_disease:.4f}")           # ~0.0458
print(f"P(disease | smoker)   = {p_both / p_smoker:.4f}")   # ~0.1500
print(f"P(disease | nonsmoker)= "
      f"{df.loc[~df['smoker'], 'disease'].mean():.4f}")     # ~0.0200

# Independence test: does P(A and B) equal P(A)*P(B)?
print(f"P(A and B)   = {p_both:.4f}")
print(f"P(A)*P(B)    = {p_disease * p_smoker:.4f}")
# 0.0300 vs 0.0092 -- far apart, so the events are strongly dependent.

# Why assuming independence is dangerous. Suppose three subsystems each
# fail with probability 0.01 and share a power supply.
print(f"if independent: {0.01 ** 3:.8f}")       # 1 in a million
# If a power failure takes out all three at once, the true joint failure
# rate is closer to the probability of THAT single common cause, perhaps
# 0.005 -- five thousand times larger than the independent estimate.

Key Takeaways

Conditioning restricts the sample space and renormalizes within it.

Every predictive model outputs a conditional probability given the features.

Independence means conditioning changes nothing, and it must be argued for rather than assumed.

Mutually exclusive events are strongly dependent, not independent; the two terms are frequently confused.

🧪 Practice

From a contingency table of 100 customers, compute P(churn | premium) and P(premium | churn) and explain why they differ.

Simulate two independent and two dependent binary variables and verify the multiplication rule for each pair.

Interview: An insurer prices flood policies assuming each house floods independently. What goes wrong? (Hint: ask what causes floods and whether neighbouring houses share that cause.)

Bayes' Theorem

Often you know the probability of evidence given a cause, but you need the probability of the cause given the evidence. A test's manufacturer reports P(positive | disease), while the patient wants P(disease | positive). These are different numbers and confusing them is called the base rate fallacy. Bayes' theorem is the conversion.

The theorem follows immediately from the definition of conditional probability written two ways:

                P(E | H) * P(H)
  P(H | E)  =  -----------------
                     P(E)

  posterior  =  likelihood * prior / evidence

  P(H)     prior       what you believed before seeing evidence
  P(E|H)   likelihood  how well the hypothesis explains the evidence
  P(E)     evidence    total probability of the evidence, a normalizer
  P(H|E)   posterior   the updated belief

The decisive term is the prior. When a condition is rare, even an accurate test produces mostly false positives, because there are so many more healthy people to generate them. This is not a flaw in the test. It is arithmetic, and it is counterintuitive to nearly everyone on first encounter.

An analogy: a car alarm that is 99% accurate still mostly signals wind, because thefts are rare and wind is not. The alarm is not broken; the base rate is doing the work.

def bayes(prior, sensitivity, specificity):
    """P(disease | positive test).

    sensitivity = P(positive | disease)      -- true positive rate
    specificity = P(negative | no disease)   -- true negative rate
    """
    true_pos = sensitivity * prior
    false_pos = (1 - specificity) * (1 - prior)
    return true_pos / (true_pos + false_pos)     # the denominator is P(E)

# A very good test: 99% sensitive, 99% specific.
for prevalence in (0.5, 0.01, 0.001):
    p = bayes(prevalence, 0.99, 0.99)
    print(f"prevalence {prevalence:<6} -> P(disease | positive) = {p:.4f}")
# prevalence 0.5    -> 0.9900   as intuition expects
# prevalence 0.01   -> 0.5000   a coin flip, despite a 99% accurate test
# prevalence 0.001  -> 0.0902   nine out of ten positives are FALSE

# Why: count the people out of 100,000 at 0.1% prevalence.
#   100 have the disease  -> 99 test positive
#   99,900 are healthy    -> 999 test positive (1% of them)
#   total positives = 1,098, of which only 99 are real => 9%.

# Sequential updating: today's posterior is tomorrow's prior.
posterior = bayes(0.001, 0.99, 0.99)         # 0.0902 after one test
posterior_2 = bayes(posterior, 0.99, 0.99)   # retest, using the update
print(f"after a second independent positive: {posterior_2:.4f}")   # 0.9075
# Two independent positives move the belief decisively; one does not.

Key Takeaways

Bayes' theorem converts P(evidence | hypothesis) into P(hypothesis | evidence), which is usually what you actually want.

When the base rate is low, most positives from an accurate test are false; this is arithmetic, not a defect.

The posterior of one update becomes the prior of the next, which is how evidence accumulates.

Reporting a classifier's precision without stating the class prevalence is close to meaningless.

🧪 Practice

A test is 95% sensitive and 90% specific for a condition affecting 2% of people. Compute the probability of disease given a positive result.

Plot posterior probability against prevalence for a fixed test and find where it crosses 50%.

Interview: Your fraud model has 90% precision in testing but flags mostly legitimate transactions in production. What changed? (Hint: compare the fraud rate in your test set with the live rate.)

Random Variables

Sample spaces of raw outcomes are awkward to compute with. A random variable is a function that assigns a number to each outcome, which turns questions about events into questions about numbers, where arithmetic works.

Rolling two dice has a 36-outcome sample space, but you usually care about the sum. Define X = sum of the dice, and you now have a numeric quantity with a distribution, an average, and a spread.

Random variables come in two kinds. Discrete variables take countably many values and are described by a probability mass function, P(X = x), which gives the probability of each exact value. Continuous variables take values in a range and are described by a probability density function, where P(X = x) is zero for every single point and only intervals carry probability. Density is probability per unit length, not probability.

Both are summarized by the cumulative distribution function, F(x) = P(X <= x), which is defined the same way for either kind and is what quantiles and p-values are computed from.

  DISCRETE: mass at points          CONTINUOUS: density over a range

   P(X=x)                             f(x)
    |                                  |      ____
    |      |                           |    .'    '.
    |    | | |                         |   /        \
    |  | | | | |                       |  /          \
    |  | | | | | |                     | /            \
    +--+-+-+-+-+-+--> x                +----------------> x
       2 3 4 5 6 7                        a        b

   heights ARE probabilities          AREA is probability
   sum of all heights = 1             total area under curve = 1
   P(X = 4) is meaningful             P(X = 4) = 0 exactly
                                      P(a < X < b) = shaded area

  CDF (both kinds): F(x) = P(X <= x), non-decreasing from 0 to 1

    1 |          .----------
      |        .'
      |      .'
      |   .-'
    0 +-'------------------> x

import numpy as np
from itertools import product
from scipy import stats

# Discrete: X = sum of two dice, derived from the raw sample space.
outcomes = list(product(range(1, 7), repeat=2))
sums = [a + b for a, b in outcomes]
pmf = {s: sums.count(s) / len(sums) for s in sorted(set(sums))}
print({k: round(v, 4) for k, v in pmf.items()})
# {2: 0.0278, 3: 0.0556, ..., 7: 0.1667, ..., 12: 0.0278}
print(f"total mass = {sum(pmf.values()):.4f}")        # 1.0000

# The CDF answers range questions directly.
cdf = {s: sum(p for k, p in pmf.items() if k <= s) for s in pmf}
print(f"P(X <= 5) = {cdf[5]:.4f}")                    # 0.2778

# Continuous: density is NOT probability. It can exceed 1.
narrow = stats.norm(loc=0, scale=0.1)
print(f"density at 0 = {narrow.pdf(0):.4f}")          # 3.9894, above 1
print(f"P(X = 0)     = 0 exactly")
print(f"P(-0.1<X<0.1)= {narrow.cdf(0.1) - narrow.cdf(-0.1):.4f}")  # 0.6827

# Only intervals have probability, and they come from the CDF.
# This is why a probability density plot's y-axis is not interpretable
# as a probability, a routine misreading of histograms with density=True.

Key Takeaways

A random variable maps outcomes to numbers so probability questions become arithmetic.

Discrete variables have a mass function where heights are probabilities; continuous ones have a density where areas are.

For a continuous variable the probability of any exact value is zero, and density can exceed one.

The cumulative distribution function works for both kinds and is what quantiles and p-values come from.

🧪 Practice

Build the probability mass function for the maximum of two dice and confirm it sums to one.

Show numerically that a normal density with small scale exceeds 1 at its peak, and explain why that is not a contradiction.

Interview: What does the y-axis mean on a density histogram, and why can it exceed 1? (Hint: check what has to integrate to one.)

Expectation, Variance, Covariance

A distribution can be complicated, so we summarize it with a few numbers. The three that matter most describe where it sits, how spread out it is, and how two variables move together.

Expectation, E[X], is the probability-weighted average: the long-run mean over many repetitions. It is the balance point of the distribution, not necessarily a value the variable can take. The expected value of a die roll is 3.5.

Variance, Var(X) = E[(X - E[X])^2], is the expected squared deviation from the mean. Squaring makes deviations positive and penalizes large ones disproportionately. Its square root, the standard deviation, returns to the original units and is what you should usually report.

Covariance measures joint variation: positive when the two tend to exceed their means together, negative when one rises as the other falls. Its size is uninterpretable because it carries the product of both units, which is why correlation, covariance divided by both standard deviations, exists. It is bounded in [-1, 1] and is unitless.

The critical limitation: correlation measures linear association only. A perfect deterministic relationship can have zero correlation.

  EXPECTATION IS A BALANCE POINT       VARIANCE IS SPREAD

    |  |                                 low variance    high variance
    |  |  |                                  ||            |  |  |
   _|__|__|__|___                          __||||__      __|__|__|__
       ^                                     ^               ^
     E[X]                                 same mean, different spread

  CORRELATION IS LINEAR ASSOCIATION ONLY

   r = +0.9        r = -0.9        r = 0.0         r = 0.0 (!)
     .  .'           '.  .            . . .          .  .  .
    . .'              '. .           . . . .        .       .
   . .'                '. .           . . .        .         .
  ..'                    '..         . . . .        '.     .'
                                                       '''''
   strong linear   strong inverse   no pattern     PERFECT parabola,
                                                   zero correlation

  Always plot. Correlation collapses a scatter into one number and
  discards exactly the structure that would have warned you.

import numpy as np

rng = np.random.default_rng(0)

# Expectation of a die: a value the die can never show.
faces = np.arange(1, 7)
print(f"E[X] = {(faces * (1/6)).sum():.2f}")          # 3.50

x = rng.normal(loc=10, scale=3, size=100_000)
print(f"mean={x.mean():.3f} var={x.var(ddof=1):.3f} sd={x.std(ddof=1):.3f}")
# ddof=1 gives the unbiased sample variance; ddof=0 (the NumPy default)
# divides by n and underestimates the population variance.

# Covariance and correlation on correlated data.
a = rng.normal(size=1000)
b = 2 * a + rng.normal(scale=0.5, size=1000)
cov = np.cov(a, b)[0, 1]
corr = np.corrcoef(a, b)[0, 1]
print(f"cov={cov:.3f}  corr={corr:.3f}")              # cov=2.01 corr=0.970

# Scaling changes covariance but not correlation -- which is the point.
print(f"cov after scaling b by 100: {np.cov(a, 100 * b)[0,1]:.1f}")  # 201.0
print(f"corr after scaling:         {np.corrcoef(a, 100*b)[0,1]:.3f}")  # 0.970

# The trap: perfect dependence, zero correlation.
u = rng.uniform(-3, 3, size=10_000)
v = u ** 2                                   # v is fully determined by u
print(f"corr(u, u^2) = {np.corrcoef(u, v)[0,1]:.4f}")   # ~0.0
# Reporting "no correlation, therefore no relationship" here would be
# exactly wrong. The relationship is perfect; it is just not linear.

Key Takeaways

Expectation is the probability-weighted mean and need not be an attainable value.

Variance squares deviations, so it punishes large ones and is reported in squared units; prefer the standard deviation.

Correlation is covariance rescaled to [-1, 1] and is unitless, which makes it comparable across variables.

Zero correlation rules out a linear relationship only, so always plot the scatter as well.

🧪 Practice

Compute expectation and variance for a biased coin paying 10 on heads and -5 on tails with P(heads) = 0.3.

Generate data with correlation near zero but obvious visual structure, and plot it.

Interview: Two stocks each have 8% expected return. One has twice the variance. Why might an investor still prefer it? (Hint: variance penalizes upside and downside identically.)

Discrete Distributions

Rather than deriving a distribution for every new problem, it pays to recognize the handful of standard shapes that recur. Each arises from a specific generative story, and identifying the story tells you which distribution applies and what its parameters mean.

The Bernoulli distribution models a single trial with two outcomes and probability p of success. The binomial counts successes in n independent Bernoulli trials with the same p. The Poisson counts events in a fixed interval when they occur independently at a constant average rate lambda; it is the limit of the binomial when n is huge and p is tiny. The geometric counts trials until the first success. The multinomial generalizes the binomial to more than two categories.

The commonest modelling error is applying a Poisson to overdispersed count data. Poisson assumes variance equals the mean. Real counts, such as website visits or insurance claims, are usually more variable than that, and the negative binomial exists precisely for those.

Distribution	Generative story	Parameters	Mean	Variance
Bernoulli	One trial, success or failure	`p`	`p`	`p(1-p)`
Binomial	Successes in `n` independent trials	`n`, `p`	`np`	`np(1-p)`
Poisson	Events per interval at constant rate	`lambda`	`lambda`	`lambda`
Geometric	Trials until the first success	`p`	`1/p`	`(1-p)/p^2`
Negative binom.	Overdispersed counts	`r`, `p`	varies	exceeds mean
Multinomial	Counts across `k` categories in `n` trials	`n`, `p_k`	`n*p_k`	`n*p_k(1-p_k)`

  BINOMIAL n=20, p=0.3             POISSON lambda=3

   P(X=k)                            P(X=k)
    |     ||                          |  ||
    |    ||||                         |  |||
    |   ||||||                        | ||||
    |  ||||||||                       | |||||
    +--+-+-+-+-+-+---> k              +-+-+-+-+-+-+---> k
       2 4 6 8 10                      0 2 4 6 8

   bounded above by n                 unbounded, right-skewed
   "how many of 20 converted?"        "how many arrive per hour?"

import numpy as np
from scipy import stats

# Binomial: 20 visitors, each converts with probability 0.3.
binom = stats.binom(n=20, p=0.3)
print(f"P(exactly 6)   = {binom.pmf(6):.4f}")          # 0.1916
print(f"P(at least 10) = {1 - binom.cdf(9):.4f}")      # 0.0480
print(f"mean={binom.mean():.2f} var={binom.var():.2f}")  # 6.00, 4.20

# Poisson: support tickets arriving at 3 per hour.
pois = stats.poisson(mu=3)
print(f"P(zero tickets) = {pois.pmf(0):.4f}")          # 0.0498
print(f"P(more than 6)  = {1 - pois.cdf(6):.4f}")      # 0.0335

# The Poisson assumption, stated and testable: variance == mean.
rng = np.random.default_rng(0)
real_counts = np.r_[rng.poisson(1, 800), rng.poisson(20, 200)]  # mixed rates
mean, var = real_counts.mean(), real_counts.var(ddof=1)
print(f"mean={mean:.2f} variance={var:.2f} ratio={var/mean:.2f}")
# mean=4.80 variance=53.84 ratio=11.22 -- massively overdispersed.
# A Poisson model here will produce far too narrow prediction intervals
# and will look badly calibrated. Use a negative binomial instead.

# Binomial converges to Poisson when n is large and p is small.
print(stats.binom(n=10_000, p=0.0003).pmf(3).round(6))   # 0.224042
print(stats.poisson(mu=3).pmf(3).round(6))               # 0.224042

Key Takeaways

Each standard distribution encodes a generative story; identify the story and the distribution follows.

Binomial counts successes in a fixed number of trials; Poisson counts events in an interval with no fixed maximum.

Poisson assumes variance equals the mean, which real count data frequently violates.

Overdispersed counts need a negative binomial, or prediction intervals will be far too narrow.

🧪 Practice

A page converts at 4%. Using the binomial, find the probability of at least 10 conversions from 150 visitors.

Test whether a count column in a real dataset is overdispersed by comparing its variance to its mean.

Interview: When is a Poisson model inappropriate for daily website visits? (Hint: consider weekends, marketing campaigns, and whether the rate is really constant.)

Continuous Distributions

Continuous distributions describe quantities measured on a scale rather than counted: durations, revenues, measurement errors, latencies. As with discrete distributions, a handful of shapes cover most situations, and recognizing them saves work.

The normal (Gaussian) distribution is the default for quantities formed by adding many small independent effects, which the central limit theorem explains. It is symmetric, and fully described by its mean and standard deviation.

The uniform distribution assigns equal density across a range and encodes complete ignorance within bounds. The exponential models waiting times between Poisson events and is memoryless. The log-normal applies when a quantity is the product of many small effects rather than their sum, giving the right-skewed shape typical of incomes, file sizes, and session lengths. The beta distribution lives on [0, 1] and is the natural choice for modelling a probability or a rate.

The most consequential mistake is assuming normality for skewed data. Revenue per customer, request latency, and time on site are almost never normal, and treating them as such produces confidence intervals that are wrong in both directions.

  NORMAL              LOG-NORMAL            EXPONENTIAL
     _--_                |                  |
    /    \               |\                 |\
   /      \              | \___              \
  /        \             |     ----__         \___
 -----------          ---+-----------       ------------
  symmetric           right-skewed,          decays from 0,
  mean = median       mean > median          memoryless

  UNIFORM              BETA(2,5)             BETA(5,5)
  +---------+            _-_                    _--_
  |         |           /   \_                 /    \
  |         |          /      ---_            /      \
 -+---------+---     -+----------+--        -+--------+--
  0         1         0          1           0        1
  ignorance           rate near 0.25         rate near 0.5

  WHY LOG-NORMAL MATTERS: mean != median for skewed data.
  Reporting "average revenue per user" on a log-normal spend
  distribution describes almost nobody in the dataset.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# Normal: the 68-95-99.7 rule, worth knowing by heart.
z = stats.norm()
for k in (1, 2, 3):
    print(f"within {k} sd: {z.cdf(k) - z.cdf(-k):.4f}")
# 0.6827, 0.9545, 0.9973

# Log-normal revenue: the shape almost all spend data actually has.
revenue = rng.lognormal(mean=3.0, sigma=1.0, size=100_000)
print(f"mean   = {revenue.mean():.2f}")       # 33.13
print(f"median = {np.median(revenue):.2f}")   # 20.07
print(f"90th   = {np.percentile(revenue, 90):.2f}")   # 72.29
# The mean sits at roughly the 65th percentile. "Average revenue" is
# above what two thirds of customers actually spend.

# Test normality before assuming it.
stat, p = stats.shapiro(revenue[:5000])       # Shapiro-Wilk, capped sample
print(f"shapiro p = {p:.2e}")                 # ~0 -> reject normality
log_rev = np.log(revenue)
print(f"after log, p = {stats.shapiro(log_rev[:5000])[1]:.3f}")  # ~0.5

# Exponential is memoryless: having waited does not change the future.
exp = stats.expon(scale=10)                   # mean wait 10 minutes
p_gt_5 = 1 - exp.cdf(5)
p_gt_15_given_10 = (1 - exp.cdf(15)) / (1 - exp.cdf(10))
print(f"{p_gt_5:.4f} vs {p_gt_15_given_10:.4f}")   # identical: 0.6065

Key Takeaways

Normal arises from sums of small effects; log-normal from products, which is why skewed positive quantities are usually log-normal.

For skewed data the mean and median differ substantially, and the mean describes almost no individual.

The exponential distribution is memoryless: elapsed waiting time carries no information about remaining wait.

Test the normality assumption rather than inheriting it, and consider a log transform when it fails.

🧪 Practice

Generate log-normal data, compare mean and median, and show how a log transform makes it approximately normal.

Model support wait times as exponential and compute the probability of waiting more than twice the mean.

Interview: Why is reporting average revenue per user often misleading? (Hint: sketch the distribution and mark where the mean falls.)

Joint and Marginal Distributions

Single-variable distributions are insufficient because data has multiple columns that interact. Joint distributions describe several random variables at once, and the operations of marginalizing and conditioning let you move between the full picture and any slice of it.

The joint distribution P(X, Y) gives the probability of each combination. The marginal distribution of X is obtained by summing or integrating the joint over all values of Y, which discards information about Y entirely. The name comes from the practice of writing these sums in the margins of a table.

Conditioning fixes Y at a value and renormalizes the remaining row or column, giving P(X | Y = y).

The essential warning is that marginal and conditional relationships can point in opposite directions. Simpson's paradox is the name for this, and it is not an exotic curiosity. It appears routinely in real analyses whenever a lurking variable differs across groups, and it is the reason a result that holds overall can reverse in every subgroup.

  JOINT, MARGINAL, CONDITIONAL IN ONE TABLE

                 Y = rainy   Y = dry   |  MARGINAL P(X)
   X = commute      0.15      0.25     |     0.40
   X = wfh          0.10      0.50     |     0.60
   ---------------------------------- +
   MARGINAL P(Y)    0.25      0.75     |     1.00
                     ^
                sum down the column -> marginal of Y
                sum across the row  -> marginal of X

   CONDITIONAL: fix Y = rainy, then renormalize that column by 0.25
     P(commute | rainy) = 0.15 / 0.25 = 0.60
     P(wfh     | rainy) = 0.10 / 0.25 = 0.40   (sums to 1, as required)

import pandas as pd
import numpy as np

# Simpson's paradox: admissions across two departments.
# Department A admits generously; department B admits almost nobody.
applicants = pd.DataFrame({
    "dept":     ["A", "A", "B", "B"],
    "gender":   ["M", "F", "M", "F"],
    "applied":  [800, 100, 200, 900],
    "admitted": [496,  82,  12,  63],
})
applicants["rate"] = applicants["admitted"] / applicants["applied"]

# Within EVERY department, women are admitted at a higher rate.
print(applicants.pivot(index="dept", columns="gender", values="rate").round(3))
#        F      M
# A  0.820  0.620
# B  0.070  0.060

# Marginalizing over department REVERSES the conclusion.
overall = applicants.groupby("gender")[["applied", "admitted"]].sum()
overall["rate"] = overall["admitted"] / overall["applied"]
print(overall["rate"].round(3))
# F  0.145     <- lower overall
# M  0.508     <- higher overall

# The lurking variable: women applied overwhelmingly to department B,
# which admits almost nobody regardless of gender. Marginalizing threw
# away the department, and with it the explanation.

Key Takeaways

Marginalizing sums the joint over a variable and discards its information entirely.

Conditioning fixes a variable and renormalizes the remaining slice.

Simpson's paradox means a marginal relationship can reverse within every subgroup, and it occurs routinely in practice.

Before trusting an aggregate comparison, check whether the groups differ in composition on some third variable.

🧪 Practice

Build a two-way joint probability table and compute both marginals and both sets of conditionals.

Construct your own Simpson's paradox example with two groups and verify the reversal numerically.

Interview: Treatment A beats B overall, but B wins for both mild and severe cases. Which do you recommend? (Hint: ask how patients were assigned to each treatment.)

Law of Large Numbers

Why should a sample tell you anything about a population? The law of large numbers is the guarantee, and it is the reason statistics works at all.

The law states that as the sample size grows, the sample mean converges to the true expected value. Averaging many independent observations cancels their random deviations, because they are equally likely to fall above and below the mean.

Two points are routinely misunderstood. First, convergence is slow: the error shrinks in proportion to one over the square root of n, so quadrupling the sample halves the error. Getting one more decimal place of precision costs a hundred times more data.

Second, the law says nothing about short runs. The gambler's fallacy is the belief that a run of heads makes tails "due", as though the process compensates. It does not. Convergence happens because later observations dilute the early deviation, not because they offset it.

  CONVERGENCE OF A RUNNING MEAN (fair coin, true value 0.5)

  1.0 |*
      | *
      |  **                                  error shrinks like 1/sqrt(n):
  0.5 |----**---*--**----*---------------     n=100   -> about 0.05
      |      *    *    **     ***========     n=10000 -> about 0.005
      |                                       100x the data for 10x
  0.0 +------------------------------------>     the precision
      1    10    100   1000   10000   n

  The early wild swings are never "corrected". They are DILUTED:
  after 10,000 flips, an excess of 5 heads moves the mean by 0.0005.

import numpy as np

rng = np.random.default_rng(0)
flips = rng.integers(0, 2, size=100_000)
running = np.cumsum(flips) / np.arange(1, len(flips) + 1)

for n in (10, 100, 1_000, 10_000, 100_000):
    print(f"n={n:>7}: mean={running[n-1]:.5f}  error={abs(running[n-1]-0.5):.5f}")
# n=     10: mean=0.30000  error=0.20000
# n=    100: mean=0.47000  error=0.03000
# n=   1000: mean=0.49300  error=0.00700
# n=  10000: mean=0.50130  error=0.00130
# n= 100000: mean=0.49934  error=0.00066

# The 1/sqrt(n) rate, made explicit.
for n in (100, 400, 1600):
    errors = [abs(rng.integers(0, 2, n).mean() - 0.5) for _ in range(2000)]
    print(f"n={n:>5}: average error = {np.mean(errors):.4f}")
# n=  100: 0.0397     quadrupling n
# n=  400: 0.0200     roughly halves the error
# n= 1600: 0.0100

# Gambler's fallacy, tested directly.
after_three_heads = [flips[i+3] for i in range(len(flips)-3)
                     if flips[i] == flips[i+1] == flips[i+2] == 1]
print(f"P(heads | three heads) = {np.mean(after_three_heads):.4f}")   # ~0.50
# Exactly the base rate. The coin has no memory of the run.

Key Takeaways

The sample mean converges to the population mean as the sample grows, which is why sampling works.

Error shrinks like one over the square root of n, so precision is expensive.

Deviations are diluted by later data, never corrected by it.

The law describes long-run behaviour and says nothing useful about any short sequence.

🧪 Practice

Simulate a running mean for a fair die and plot the error against n on log axes to see the rate.

Determine empirically how many samples are needed to estimate a proportion within 1 percentage point.

Interview: An A/B test shows a 30% lift after 40 users. Do you ship? (Hint: estimate the error you would expect at that sample size.)

Central Limit Theorem

The law of large numbers says the sample mean converges. The central limit theorem says something much stronger and far more useful: it tells you the shape of the distribution of that sample mean, and the shape barely depends on the population it came from.

The theorem states that the sample mean of n independent observations is approximately normally distributed for large n, centred at the true mean, with standard deviation equal to the population standard deviation divided by the square root of n. That last quantity is the standard error.

The remarkable part is the universality. The underlying data can be skewed, bimodal, or discrete; the distribution of its mean still approaches a normal. This is why confidence intervals and t-tests work on data that is visibly non-normal, a point that confuses many beginners into over-testing for normality of the raw data.

The caveats matter. Convergence needs a finite variance, so heavy-tailed distributions can fail. Convergence is slower for more skewed populations, and the theorem describes the sample mean, not the raw data, which stays exactly as skewed as it ever was.

  POPULATION (heavily skewed)      DISTRIBUTION OF ITS SAMPLE MEAN

   |\                               n = 2      n = 30
   | \____                            /\         _-_
   |      -----___                   /  \_      /   \
   +---------------                 /     \_   /     \
                                  -+--------  -+-------+-
   nothing like a normal           still skew  nearly normal

  STANDARD ERROR SHRINKS WITH SAMPLE SIZE

   SE = sigma / sqrt(n)

   n:      10      40     160     640
   SE:  0.316   0.158   0.079   0.040       (for sigma = 1)
        \_____/  \_____/ \_____/
         4x the n halves the SE each time

  This single formula is why sample size planning is possible at all.

import numpy as np
from scipy import stats

rng = np.random.default_rng(0)

# A deliberately non-normal population: exponential, heavily right-skewed.
population = rng.exponential(scale=10, size=1_000_000)
print(f"population skew = {stats.skew(population):.3f}")     # ~2.0

# Take many samples of size n and look at the distribution of their means.
for n in (2, 5, 30, 100):
    means = population[rng.integers(0, len(population), size=(5000, n))].mean(axis=1)
    observed_se = means.std(ddof=1)
    predicted_se = population.std() / np.sqrt(n)      # sigma / sqrt(n)
    print(f"n={n:>3}: skew of means={stats.skew(means):>6.3f}  "
          f"SE observed={observed_se:.3f} predicted={predicted_se:.3f}")
# n=  2: skew of means= 1.410  SE observed=7.072 predicted=7.070
# n=  5: skew of means= 0.895  SE observed=4.475 predicted=4.472
# n= 30: skew of means= 0.365  SE observed=1.826 predicted=1.826
# n=100: skew of means= 0.199  SE observed=1.001 predicted=1.000
# The skew of the means falls toward zero; the SE formula is exact.

# This is what licenses a confidence interval on skewed data.
sample = rng.choice(population, size=200, replace=False)
se = sample.std(ddof=1) / np.sqrt(len(sample))
ci = (sample.mean() - 1.96 * se, sample.mean() + 1.96 * se)
print(f"95% CI for the mean: ({ci[0]:.2f}, {ci[1]:.2f}), truth = "
      f"{population.mean():.2f}")
# The raw data is wildly non-normal; the interval is still valid, because
# the theorem applies to the MEAN, not to the observations.

Key Takeaways

The sample mean is approximately normal for large n regardless of the population's shape, given finite variance.

The standard error is the population standard deviation divided by the square root of the sample size.

This universality is what makes confidence intervals and t-tests valid on non-normal data.

The theorem applies to the mean, not to the raw observations, and converges more slowly for more skewed populations.

🧪 Practice

Simulate sample means from a uniform, an exponential, and a bimodal population and plot each at n = 2, 10, and 50.

Verify the standard error formula empirically for three sample sizes.

Interview: Your revenue data is heavily skewed. Can you still build a confidence interval for the mean? (Hint: ask which quantity the theorem actually describes.)

2.4 Discrete Mathematics and Information Theory

Counting and set reasoning underpin probability itself, while information theory supplies the measures of surprise and shared signal that decide how decision trees split, how classifiers are trained, and how features are selected.

Combinatorics and Counting

Probability on equally likely outcomes reduces to counting: how many outcomes satisfy the event, divided by how many exist. Doing that reliably requires knowing which counting rule applies, and the choice hinges on two questions. Does order matter, and can items repeat?

The product rule is the foundation: if one choice has m options and an independent second choice has n, the pair has m * n possibilities. Everything else follows from it.

A permutation counts arrangements where order matters. The number of ways to order k items chosen from n is n! / (n-k)!. A combination counts selections where order does not matter, written C(n, k) and equal to n! / (k! (n-k)!). The division by k! removes the orderings you just decided to ignore.

Combinatorics matters in practice for computing the size of a search space. Before launching a grid search, count its points. The reason hyperparameter grids explode is multiplicative growth, and the reason random search often beats grid search is that the grid spends its budget on redundant combinations along dimensions that do not matter.

  THE DECISION TREE FOR COUNTING

   Does order matter?
        |
        +-- YES --> Can items repeat?
        |            +-- YES --> n^k           (passwords, sequences)
        |            +-- NO  --> n!/(n-k)!     (rankings, podium places)
        |
        +-- NO  --> Can items repeat?
                     +-- YES --> C(n+k-1, k)   (multiset selections)
                     +-- NO  --> C(n, k)       (committees, hands, subsets)

  CHOOSING 3 FROM {A, B, C, D}

   PERMUTATIONS (order matters): 4!/(4-3)! = 24
     ABC ABD ACB ACD ADB ADC BAC BAD ... (all 24 arrangements)

   COMBINATIONS (order ignored): C(4,3) = 4
     ABC  ABD  ACD  BCD
     each corresponds to 3! = 6 permutations, hence 24/6 = 4

from math import comb, perm, factorial
from itertools import combinations, product

print(perm(4, 3))      # 24  -- ordered selections
print(comb(4, 3))      # 4   -- unordered selections
print(comb(52, 5))     # 2,598,960 -- distinct five-card poker hands

# Probability of a flush, computed purely by counting.
flushes = 4 * comb(13, 5)          # pick a suit, then 5 of its 13 ranks
print(f"P(flush) = {flushes / comb(52, 5):.6f}")     # 0.001981

# Where this bites in practice: grid search size is a product.
grid = {
    "n_estimators":   [100, 300, 500],
    "max_depth":      [3, 5, 7, 10],
    "learning_rate":  [0.01, 0.05, 0.1],
    "subsample":      [0.7, 0.85, 1.0],
    "min_child_wt":   [1, 5, 10],
}
sizes = [len(v) for v in grid.values()]
total = 1
for s in sizes:
    total *= s                      # the product rule, literally
print(f"grid points = {total}")     # 324
print(f"with 5-fold CV: {total * 5} model fits")   # 1620

# Adding ONE more three-valued parameter triples the cost.
print(f"one more parameter: {total * 3 * 5} fits")  # 4860
# This multiplicative blowup is why random search over the same budget
# usually wins: it samples distinct values along every axis, whereas the
# grid repeats each value many times over the irrelevant dimensions.

Key Takeaways

Choose a counting rule by asking whether order matters and whether items may repeat.

Combinations divide permutations by k! to remove orderings you have decided to ignore.

Search spaces grow multiplicatively, so each added hyperparameter multiplies rather than adds cost.

Random search beats grid search at equal budget partly because the grid wastes evaluations repeating values along unimportant axes.

🧪 Practice

Compute the probability of being dealt a full house using counting arguments only.

Count the total fits implied by a grid you have used, including cross-validation folds.

Interview: Why does random search often outperform grid search at the same budget? (Hint: consider how many distinct values of the single most important hyperparameter each method actually tries.)

Set Theory Basics

Sets are the substrate underneath probability, database joins, and feature engineering. A set is an unordered collection of distinct elements, and the distinctness is what makes it useful: it models membership without duplication.

The operations have direct data-work meanings. Union combines, giving everything in either set. Intersection keeps only shared elements. Difference removes one set's elements from another. Symmetric difference keeps what is in exactly one. The cardinality is the count of elements.

Two results are worth remembering. De Morgan's laws state that the complement of a union is the intersection of complements, and vice versa, which is how a filter condition is correctly negated. The inclusion- exclusion principle gives the size of a union without double-counting, and is the counting version of the probability rule from earlier.

The immediate practical payoff is that set operations are how you reason about joins and about overlaps between cohorts, and they are dramatically faster than looping when implemented as hash sets.

  SET OPERATIONS AND THEIR SQL-JOIN COUNTERPARTS

   A union B          A intersect B       A minus B
   +-----+-----+      +-----+-----+       +-----+-----+
   |#####|#####|      |     |     |       |#####|     |
   |##+--+--+##|      |  +--+--+  |       |##+--+--+  |
   |##|#####|##|      |  |#####|  |       |##|  |  |  |
   +--+--+--+--+      +--+--+--+--+       +--+--+--+--+
      |     |            |     |             |     |
      +-----+            +-----+             +-----+
   FULL OUTER JOIN    INNER JOIN          LEFT ANTI JOIN

  DE MORGAN'S LAWS
   not (A or B)  ==  (not A) and (not B)
   not (A and B) ==  (not A) or  (not B)

   Applied to a filter: the opposite of
     "premium AND active"
   is NOT "not premium AND not active", it is
     "not premium OR not active".
   Getting this wrong silently drops rows from a cohort.

  INCLUSION-EXCLUSION FOR THREE SETS
   |A u B u C| = |A|+|B|+|C| - |AnB| - |AnC| - |BnC| + |AnBnC|
                              \___ remove doubles ___/  \_ add back _/

mobile = {"u1", "u2", "u3", "u4", "u7"}
web    = {"u3", "u4", "u5", "u6"}
email  = {"u1", "u4", "u6", "u8"}

print(mobile | web)            # union: any channel
print(mobile & web)            # intersection: both channels
print(mobile - web)            # difference: mobile only
print(mobile ^ web)            # symmetric difference: exactly one

# Inclusion-exclusion, verified against the direct count.
formula = (len(mobile) + len(web) + len(email)
           - len(mobile & web) - len(mobile & email) - len(web & email)
           + len(mobile & web & email))
print(formula, len(mobile | web | email))       # 8 8 -- they agree

# De Morgan in practice: negating a compound filter.
all_users = mobile | web | email
premium = {"u1", "u3", "u6"}
active  = {"u1", "u2", "u3", "u4"}

target = premium & active                       # "premium AND active"
wrong  = (all_users - premium) & (all_users - active)   # NOT the complement
right  = (all_users - premium) | (all_users - active)   # correct complement

print(sorted(target))                            # ['u1', 'u3']
print(sorted(wrong))                             # ['u5', 'u7', 'u8']  too few
print(sorted(right))                             # everything except u1, u3
print(target | right == all_users)               # True -- a valid partition

# Membership testing is O(1) in a set and O(n) in a list. On a million
# ids this is the difference between milliseconds and minutes.

Key Takeaways

Set operations map directly onto join types and cohort overlap questions.

De Morgan's laws govern how compound filters negate, and getting them wrong silently drops rows.

Inclusion-exclusion counts a union without double-counting the overlaps.

Hash-set membership is constant time, so replacing list scans with sets is often the single biggest easy speedup.

🧪 Practice

Given three campaign audiences, compute how many users saw at least one and how many saw all three.

Write the correct negation of "active AND (premium OR trial)" and verify it on a small set.

Interview: An anti-join returns more rows than the left table. What went wrong? (Hint: check whether the join key is unique on the right side.)

Entropy

How much uncertainty is in a random variable? "A lot" is not actionable. Entropy makes it a number, and that number turns out to govern decision tree splits, compression limits, and the loss functions used to train classifiers.

The starting intuition is surprise. An outcome that was almost certain carries little information when it occurs. A rare outcome carries a lot. Defining surprise as -log2(p) captures this: it is zero for a certain event and grows without bound as the probability shrinks. Using base 2 makes the unit the bit.

Entropy is the expected surprise: H(X) = -sum p(x) * log2(p(x)). It measures the average uncertainty per observation, and equivalently, the minimum average number of bits needed to encode outcomes from the distribution.

Entropy is maximized by the uniform distribution, where every outcome is equally likely and you learn the most from an observation. It is zero when one outcome has probability 1, since there is nothing left to learn. This is exactly the property decision trees exploit: a split is good when it reduces the entropy of the label within the resulting groups.

An analogy: entropy is the average number of yes-or-no questions needed to pin down the outcome, when you ask optimally. Guessing a fair coin needs one. Guessing one of eight equally likely items needs three.

  ENTROPY OF A BINARY VARIABLE

   H(p)
  1.0 |          _.-''-._
      |        .'        '.            maximum at p = 0.5:
      |      .'            '.          maximum uncertainty
  0.5 |    .'                '.
      |  .'                    '.
  0.0 +-'------------------------'--> p
      0        0.25  0.5  0.75      1
      ^                             ^
   certain (H=0)               certain (H=0)

  ENTROPY AS OPTIMAL ENCODING LENGTH

   Four symbols, UNIFORM: p = 1/4 each
     H = 4 * (1/4 * log2 4) = 2 bits    -> encode as 00, 01, 10, 11

   Four symbols, SKEWED: 1/2, 1/4, 1/8, 1/8
     H = 1.75 bits    -> encode as 0, 10, 110, 111
     the frequent symbol gets the short code, so the AVERAGE
     message is shorter than 2 bits. That is compression.

import numpy as np

def entropy(probs, base=2):
    p = np.asarray([x for x in probs if x > 0])     # 0 * log 0 is taken as 0
    return float(-(p * np.log(p) / np.log(base)).sum())

print(entropy([0.5, 0.5]))                  # 1.0    fair coin, one bit
print(entropy([0.9, 0.1]))                  # 0.469  biased, less uncertain
print(entropy([1.0]))                       # 0.0    certain, no information
print(entropy([0.25] * 4))                  # 2.0    four equally likely
print(entropy([0.5, 0.25, 0.125, 0.125]))   # 1.75   skewed, compressible

# Decision trees pick the split that reduces entropy the most.
def information_gain(parent_labels, left_labels, right_labels):
    def h(labels):
        _, counts = np.unique(labels, return_counts=True)
        return entropy(counts / counts.sum())
    n = len(parent_labels)
    weighted_child = (len(left_labels) / n * h(left_labels)
                      + len(right_labels) / n * h(right_labels))
    return h(parent_labels) - weighted_child      # the gain from splitting

parent = np.array([0]*50 + [1]*50)            # perfectly mixed: H = 1.0
clean_left, clean_right = np.zeros(50), np.ones(50)     # a perfect split

print(f"perfect split gain = "
      f"{information_gain(parent, clean_left, clean_right):.4f}")   # 1.0000
useless_l = np.array([0]*25 + [1]*25)
useless_r = np.array([0]*25 + [1]*25)
print(f"useless split gain = "
      f"{information_gain(parent, useless_l, useless_r):.4f}")      # 0.0000
# A split that leaves both children as mixed as the parent gains nothing,
# and the tree will not choose it. That comparison IS the algorithm.

Key Takeaways

Surprise is -log2(p), and entropy is its expected value, measured in bits.

Entropy peaks at the uniform distribution and is zero when an outcome is certain.

Entropy equals the minimum average code length, which is why skewed distributions compress.

Decision trees choose splits by maximizing the entropy reduction, called information gain.

🧪 Practice

Compute the entropy of a fair six-sided die and of a die loaded to show six half the time.

Implement information gain and use it to pick the best split on a small labelled dataset.

Interview: Why does a decision tree prefer a split producing pure child nodes? (Hint: express purity in terms of the entropy of each child.)

Cross-Entropy and KL Divergence

Entropy measures the uncertainty in one distribution. Training a classifier requires comparing two: the true labels and the model's predicted probabilities. Cross-entropy and Kullback-Leibler divergence are the tools for that, and cross-entropy is the loss function behind essentially every classifier and language model in use.

Cross-entropy H(p, q) is the average number of bits needed to encode outcomes from the true distribution p using a code optimized for the predicted distribution q. If the prediction is perfect, cross-entropy equals the entropy of p and is as small as it can be. If the prediction is wrong, you pay extra bits.

KL divergence is exactly that penalty: KL(p || q) = H(p, q) - H(p). It is the excess cost of using the wrong distribution, is always non-negative, and is zero only when the two distributions match.

KL divergence is not a distance, because it is asymmetric: KL(p || q) differs from KL(q || p). The asymmetry has real consequences in variational inference, where the choice of direction determines whether an approximation covers all the modes or concentrates on one.

The practical property that matters most: cross-entropy punishes confident wrong predictions brutally, because the loss involves the log of the predicted probability, and the log of a small number is a large negative.

  THE THREE QUANTITIES

   H(p)      = cost of encoding p with the BEST possible code
   H(p, q)   = cost of encoding p with a code built for q
   KL(p||q)  = H(p,q) - H(p)  = the WASTE from using the wrong code

         |<------------- H(p, q) ------------->|
         |<------ H(p) ------>|<-- KL(p||q) -->|
         |   irreducible      |   your model's |
         |   uncertainty      |   error        |

  CROSS-ENTROPY LOSS PER PREDICTION (true label = 1)

   predicted p:  0.99   0.9    0.7    0.5    0.1    0.01   0.001
   loss -ln(p):  0.01   0.11   0.36   0.69   2.30   4.61   6.91
                                                     ^
                          confident and wrong is catastrophically costly,
                          which is exactly the incentive you want

import numpy as np

def cross_entropy(p, q, eps=1e-12):
    p, q = np.asarray(p, float), np.clip(np.asarray(q, float), eps, 1)
    return float(-(p * np.log(q)).sum())        # natural log -> nats

def kl_divergence(p, q, eps=1e-12):
    p = np.clip(np.asarray(p, float), eps, 1)
    q = np.clip(np.asarray(q, float), eps, 1)
    return float((p * np.log(p / q)).sum())

p = np.array([1.0, 0.0, 0.0])          # true label: class 0 (one-hot)

for q in ([0.9, 0.05, 0.05], [0.5, 0.3, 0.2], [0.05, 0.9, 0.05]):
    print(f"pred={q}  CE={cross_entropy(p, q):.4f}")
# pred=[0.9, 0.05, 0.05]   CE=0.1054   confident and right: cheap
# pred=[0.5, 0.3, 0.2]     CE=0.6931   uncertain: moderate
# pred=[0.05, 0.9, 0.05]   CE=2.9957   confident and WRONG: expensive

# When p is one-hot, H(p) = 0, so cross-entropy and KL coincide. This is
# why classifier training can minimize either and call it the same thing.
print(kl_divergence(p, [0.9, 0.05, 0.05]))     # 0.1054, identical

# Asymmetry: KL is not a distance.
a = np.array([0.7, 0.2, 0.1])
b = np.array([0.3, 0.4, 0.3])
print(f"KL(a||b) = {kl_divergence(a, b):.4f}")   # 0.3446
print(f"KL(b||a) = {kl_divergence(b, a):.4f}")   # 0.3527
# Different numbers. If you need symmetry, use Jensen-Shannon divergence,
# which averages both directions and is bounded.

# Why clipping matters: an unclipped log(0) is -inf and poisons the loss.
# Every production implementation clips, and this is the reason.

Key Takeaways

Cross-entropy is the cost of encoding the truth with the model's distribution, and it is the standard classification loss.

KL divergence is the excess cost over the irreducible entropy, always non-negative and zero only for a perfect match.

KL is asymmetric and therefore not a distance; use Jensen-Shannon when symmetry is needed.

Cross-entropy punishes confident wrong predictions severely, which is why implementations clip probabilities away from zero.

🧪 Practice

Compute cross-entropy loss for three predicted distributions against a one-hot label and rank them.

Show numerically that KL(p||q) and KL(q||p) differ for two distributions of your choice.

Interview: Why is cross-entropy preferred over mean squared error for classification? (Hint: compare the gradient each supplies when the model is confidently wrong.)

Mutual Information

Correlation detects linear association and misses everything else. Mutual information detects any statistical dependence, linear or not, which makes it a far more general tool for deciding whether a feature carries signal about a target.

Mutual information I(X; Y) measures how much knowing one variable reduces uncertainty about the other. Formally it is the entropy of X minus the entropy of X given Y: the uncertainty you started with, minus what remains after learning Y. Equivalently, it is the KL divergence between the joint distribution and the product of the marginals, which is the precise statement of how far the variables are from independent.

Its key properties: it is symmetric, it is zero exactly when the variables are independent, and it is never negative. Unlike correlation, a zero value genuinely means no relationship of any kind.

The costs are real. Mutual information requires estimating a joint distribution, which means binning continuous variables, and the result is sensitive to how you bin. It is also measured in bits rather than on a fixed scale, so values are not comparable across variable pairs without normalization.

  MUTUAL INFORMATION AS OVERLAPPING UNCERTAINTY

     H(X)                      H(Y)
   .-------------------.  .-------------------.
   |                   |  |                   |
   |     H(X|Y)     .--+--+--.     H(Y|X)     |
   |                |  I(X;Y) |                |
   |                '--+--+--'                |
   '-------------------'  '-------------------'

   I(X;Y) = H(X) - H(X|Y) = H(Y) - H(Y|X)
          = the shared information, the overlap

   independent variables -> no overlap -> I = 0
   Y fully determined by X -> H(Y|X) = 0 -> I = H(Y), maximal

  WHERE IT BEATS CORRELATION

    y = x^2 relationship          correlation:        ~0.00
      .        .                  mutual info:        HIGH
       .      .
        .    .                    Correlation reports "no relationship"
         .  .                     and is wrong. MI reports the dependence
          ..                      that is plainly visible in the plot.

import numpy as np
from sklearn.feature_selection import mutual_info_regression
from sklearn.metrics import mutual_info_score

rng = np.random.default_rng(0)
x = rng.uniform(-3, 3, size=5000)

relationships = {
    "linear":     2 * x + rng.normal(scale=0.3, size=5000),
    "quadratic":  x**2 + rng.normal(scale=0.3, size=5000),
    "sinusoidal": np.sin(3 * x) + rng.normal(scale=0.3, size=5000),
    "noise":      rng.normal(size=5000),
}

for name, y in relationships.items():
    corr = abs(np.corrcoef(x, y)[0, 1])
    mi = mutual_info_regression(x.reshape(-1, 1), y, random_state=0)[0]
    print(f"{name:<11} |correlation|={corr:.3f}   mutual info={mi:.3f}")
# linear      |correlation|~0.99    mutual info HIGH
# quadratic   |correlation|~0.01    mutual info HIGH    <- corr blind, MI sees it
# sinusoidal  |correlation|~0.05    mutual info HIGH    <- same story
# noise       |correlation|~0.01    mutual info ~0      <- both agree: nothing
# Exact MI values depend on the estimator's neighbourhood parameter, so
# read the ranking rather than the magnitudes.

# For discrete variables, MI is computed directly from the joint counts.
a = rng.integers(0, 3, size=5000)
b = (a + rng.integers(0, 2, size=5000)) % 3       # dependent on a
c = rng.integers(0, 3, size=5000)                 # independent of a
print(f"I(a;b) = {mutual_info_score(a, b):.4f}")  # ~0.41 nats, dependent
print(f"I(a;c) = {mutual_info_score(a, c):.4f}")  # ~0.000, independent
# The 0.41 is not arbitrary: knowing `a` narrows `b` from three equally
# likely values to two, so I = ln(3) - ln(2) = 0.4055 nats.

# Caveat: MI on continuous data depends on the binning or the estimator's
# neighbourhood size. Compare MI values computed the same way, and treat
# the absolute number as far less trustworthy than the ranking.

Measure	Detects	Range	Symmetric	Needs binning
Pearson correlation	Linear only	`[-1, 1]`	Yes	No
Spearman correlation	Monotonic	`[-1, 1]`	Yes	No
Mutual information	Any dependence	`[0, inf)`	Yes	Yes, usually
KL divergence	Distribution mismatch	`[0, inf)`	No	Yes, usually

Key Takeaways

Mutual information measures how much one variable reduces uncertainty about another, capturing non-linear dependence.

It is zero exactly when variables are independent, unlike correlation, whose zero is much weaker evidence.

It is symmetric and non-negative, and equals the KL divergence between the joint and the product of marginals.

Estimates depend on binning, so trust the ranking of features more than the absolute values.

🧪 Practice

Generate a quadratic relationship and show that correlation is near zero while mutual information is clearly positive.

Rank the features of a dataset by mutual information with the target and compare the ordering to a correlation ranking.

Interview: When would you prefer correlation over mutual information for feature screening? (Hint: consider interpretability, sample size, and what binning does to a small dataset.)

← Data Science