1 How to use this companion

Eigen Times History Math teaches techniques that the earlier Mathematics of Eigen Times and Mathematics of Anthropology do not yet develop as complete executable lessons. Its parent, Eigen Times Math History, explains the historical arc and original publications. Here the emphasis is on learning how to perform a calculation, understanding why it works, and discovering when its answer should not be trusted.

The five editions share one text: PDF, EPUB, HTML, a Python notebook, and a native OCaml notebook. The notebooks contain the explanations, definitions, examples, exercises, glossary, references, and index. Code cells follow the relevant explanations. Each numerical example starts from declared inputs and is computed in both languages. Saved outputs are a record of a checked run; changing an input and running the cells again is how to investigate a new question.

This is a complement to the other two companions. A covariance calculation is not repeated merely because it occurs in another historical chapter. Conversely, calling a library routine is not the same as explaining the algorithm inside it. For example, the earlier OCaml setup already diagonalizes matrices and orthogonalizes columns. The new lessons expose the rotations, projections, reflections, and conditioning questions that a learner could otherwise miss.

1.1 A route through the three companions

Begin with the Eigen Times Math reader if vectors, matrix multiplication, covariance, or singular values are unfamiliar. Principal component analysis (PCA) finds orthogonal directions of variation in centered observations; its construction is taught there. Use the Anthropology Math reader for person profiles, people directions, reciprocal navigation, and proposed ontology predictors. The prerequisite catalog names exact lessons and the Python and OCaml cells that compute them. The delivered study bundle includes unchanged copies and reading versions of those four notebooks from their frozen 1.0.0 baselines, preserving the exact lesson and cell references used by this prerequisite catalog. This companion is the expanded 1.1.0 edition; separately delivered expanded companions provide their current teaching. The baseline citations are not silently repinned to changed cells. The public study repository provides a study landing page; the local bundle does not require access to a private repository.

The new companion has two connected routes. The numerical route moves from fitting observations to reliable matrix calculations, graphs, and retrieval. The information and learning route moves from uncertainty and forecast quality to formal concepts, derivatives, predictive embeddings, attention, and calibration. Read each route in order on a first pass. Later, the subject index can take you directly to a method.

A useful teaching rhythm is to read the definitions, predict the output, run the cell, inspect the invariant being tested, and change one input. An invariant is a property that the stated operation should preserve, such as an orthogonal transformation preserving squared length. A reference assertion about a fixed toy answer may need changing after a deliberate edit; an invariant should still hold whenever its assumptions hold. Cross-language agreement is an additional check, not a replacement for such mathematical reasoning.

1.2 What a teaching computation establishes

All new arrays and graphs in this book are synthetic. They are small enough to inspect and require no credentials, newspaper archive, model download, or live service. A two-layer derivative exercise does not train BERT. A scaled-attention calculation does not reproduce the BGE embedding model. A limited triple-rule closure is not an implementation of an OWL reasoner. Quantizing a handful of numbers illustrates rounding and error; it does not reproduce a production model’s complete quantization procedure.

These distinctions matter because the three sites combine mathematics with choices about data, identities, timestamps, and evidence. A cosine or learned representation can suggest related coverage. It cannot establish an employment relationship, an investment, or the reliability of a source. Likewise, a score between zero and one is not automatically a calibrated probability. The new lessons explain the computations while keeping those evidential boundaries.

Numerical examples are computed with ordinary finite-precision arithmetic. Rounding error is the discrepancy introduced when a mathematical number or operation cannot be represented exactly by the machine. Approximation error is the loss intentionally accepted by a model or truncation. A stable algorithm can reduce the first while leaving the second unchanged. Neither is the same as an incorrect identity match or an unsupported factual claim.

The notebooks execute without network access after their runtimes and packages are installed. Optional reading links and the browser’s mathematics renderer may need a connection. Use a fresh kernel and run all cells from the beginning; later lessons can reuse helpers defined earlier. The setup code is a toolbox, not a prerequisite vocabulary test: its operations are introduced in the corresponding lessons. Python uses NumPy where a library operation helps make a comparison; OCaml performs the calculations natively.

2 Notation and prerequisites before calculation

An observation is one measured item. A scalar is one number. A vector is an ordered list of numbers, and a matrix is a rectangular table of them. Observation vectors are rows in this companion; projection directions are columns. A matrix with $u$ rows and $v$ columns has shape $u\times v$, where $u$ and $v$ here are positive integer dimensions. Multiplying a $u\times v$ matrix by a $v\times k$ matrix gives a $u\times k$ matrix, with $k$ another positive dimension. Each product entry is a sum of coordinatewise products.

The symbol $\mathbb R$ denotes the real numbers. Superscript $\mathsf T$ means transpose, which exchanges rows and columns. For a vector $x$, its Euclidean norm $\lVert x\rVert$ is the square root of the sum of squared coordinates. For a matrix $M$, the Frobenius norm $\lVert M\rVert_F$ is the square root of the sum of squared entries. Its operator norm $\lVert M\rVert_2$ is the largest length of $Mx$ over column vectors $x$ of length one. The identity matrix $I_k$ leaves every $k$-coordinate column unchanged; $\mathbf1_n$ is a column of $n$ ones. The symbol $\operatorname{diag}$ forms a diagonal matrix from the listed numbers, and $\operatorname{tr}$ sums a square matrix’s diagonal entries.

A linear combination multiplies vectors by scalars and adds the results. A span is the set of all such combinations. A list is linearly independent when none of its vectors is a combination of the others. A basis is an independent list spanning a space. Its number of members is the space’s dimension. A matrix’s rank is the dimension of its column span, equal to the dimension of its row span. A nullspace consists of inputs the matrix maps to zero. These ideas explain both uniqueness and ambiguity in fitting.

The dot product of two equally sized vectors is the sum of their coordinatewise products. Two vectors are orthogonal when their dot product is zero. A set is orthonormal when its members are orthogonal and have length one. An orthogonal matrix is a square matrix whose transpose is its inverse. An eigenvector of a square matrix is a nonzero column whose image is a scalar multiple of that column; the multiplier is its eigenvalue. The earlier Eigen Times companion calculates these objects; the new lessons explain more of the numerical machinery for constructing them.

The natural logarithm is written $\ln$ or $\log$; $\log_2$ has base two. The exponential function is $\exp$, with $e=\exp(1)$. A sum sign $\sum$ adds terms over its stated index. A hat denotes an estimate or reconstruction, and a bar denotes a mean. The notation $\arg\min$ asks which argument makes an objective smallest. A derivative measures local change of a function with respect to one input; the neural lesson develops its use from an explicit calculation. A gradient is the vector of derivatives with respect to several inputs. Local auxiliary symbols are defined before the formula using them.

Reading a calculation without skipping its operations

A subscript names an entry rather than multiplying: $e_a$ is entry $a$ of a residual column $e$. If the column has $n$ entries, the notation $\sum_{a=1}^n e_a^2$ instructs us to square each entry and add all $n$ results. The squared Euclidean length is this sum; the length itself requires one more operation, the square root:

$$\lVert e\rVert^2=e^{\mathsf T}e=\sum_{a=1}^n e_a^2, \qquad \lVert e\rVert=\sqrt{\sum_{a=1}^n e_a^2}.$$

A residual can have negative entries even though its length cannot be negative. The sign records whether an individual prediction is above or below an observation; squaring removes that sign when measuring total error. Squared length and length have different units. If residual entries are measured in a unit of response, length has that unit and squared length has its square. Dividing squared length by $n$ instead gives mean squared error. Taking its square root gives root mean squared error. These quantities answer related questions but must not be exchanged under the single word “error.”

Products are read from their shapes. If $E$ has $n$ observation rows and $k$ predictor columns and $\vartheta$ is a $k\times1$ coefficient column, then $E\vartheta$ is an $n\times1$ prediction column. Its $a$th entry is $\sum_{j=1}^k E_{aj}\vartheta_j$, where $j$ temporarily indexes predictors. The transpose $E^{\mathsf T}$ is $k\times n$, so $E^{\mathsf T}e$ contains one dot product with each predictor column. Reading those entries explicitly will explain the least-squares equations, rather than asking us to memorize a matrix formula.

An equality sign says two expressions have the same value. The sign $\approx$ marks an approximation, often a rounded display. An inequality such as $x\geq0$ includes equality. Brackets $[a,b]$ describe an interval including its endpoints; set braces describe a collection of elements. Set inclusion $E\subseteq F$ means every element of $E$ belongs to $F$. A prime in the formal-concept lesson is a named set operation, while a prime on a scalar function in the fitting lesson means a derivative; each lesson states its local use.

When we differentiate a scalar function, we ask how much its output changes per small change of input. For a function $f$ of a scalar $t$, its derivative $f'(t)$ is the limit of $[f(t+\epsilon)-f(t)]/\epsilon$ as the nonzero displacement $\epsilon$ approaches zero, when that limit exists. The symbol $\partial$ denotes a partial derivative: change one input while holding the others fixed. The fitting lesson develops this idea first for a square, and the network lesson builds a chain of such local changes.

2.1 The shared series convention

The following symbols retain the meanings used by the history paper. New local algorithms declare their own temporary arrays explicitly. A locally subscripted neural quantity must not silently become a news or people matrix.

A formula copied from a source using column observations may require a transpose before it can be compared with this book. The notation is a teaching convention, not a claim that every cited paper or production source file used these letters.

3 Making the numerical steps visible

The earlier companions already calculate covariance, principal components, singular-value decompositions, and people profiles. The lessons here open the operations between those results. They use small disclosed inputs so that an answer can be checked by a second identity, not merely accepted because a library printed it. Observations remain rows. A coefficient vector or a column extracted while factoring a matrix is explicitly a column; that does not change the observation convention.

3.1 Fit an error before trusting a coefficient

Least squares starts with a choice about error. Adrien-Marie Legendre’s 1805 formulation supplies the historical reference; the small examples below are modern teaching fixtures, not reconstructions of his astronomical observations. Legendre 1805

Let the three observed responses be $y_1=1$, $y_2=2$, and $y_3=3$. Let $a$ index these observations, let $\vartheta$ be one fitted scalar, and let $\mathcal J(\vartheta)$ denote the sum of squared discrepancies. A completed-square form writes a quadratic as a nonnegative squared term plus a constant, making its minimum visible. Here it is

$$\mathcal J(\vartheta)=\sum_{a=1}^{3}(y_a-\vartheta)^2=3(\vartheta-2)^2+2.$$

The minimum is therefore at two, and the remaining squared error is two. The notebook checks the identity at several candidate values. Squaring prevents opposite discrepancies from canceling. It also makes a large discrepancy more influential, so the criterion is a modeling decision rather than a universal meaning of accuracy.

3.1.1 Expand the error before minimizing it

Here the model predicts the same number $\vartheta$ for every response. At the candidate value two, subtract prediction from observation to obtain residuals $(-1,0,1)^{\mathsf T}$. Their sum is zero, although two observations are missed. Their squared length is $(-1)^2+0^2+1^2=2$ and their length is $\sqrt2$. A zero sum of signed errors therefore does not mean a perfect fit.

To obtain the completed square, expand each term:

$$\begin{aligned} \mathcal J(\vartheta)&=(1-\vartheta)^2+(2-\vartheta)^2+(3-\vartheta)^2\\ &=3\vartheta^2-2(1+2+3)\vartheta\\ &\quad+(1+4+9)\\ &=3\vartheta^2-12\vartheta+14\\ &=3(\vartheta^2-4\vartheta+4)+2\\ &=3(\vartheta-2)^2+2. \end{aligned}$$

The square is nonnegative and becomes zero exactly at two. This establishes the minimum without calculus. Differentiation reaches the same conclusion, but it is worth seeing what that operation does. At a candidate $\vartheta$, change the parameter by a nonzero displacement $\epsilon$. Subtract the old loss from the new one and divide by the displacement:

$$\frac{\mathcal J(\vartheta+\epsilon)-\mathcal J(\vartheta)}{\epsilon} =6\vartheta-12+3\epsilon.$$

As $\epsilon$ tends to zero, the last term vanishes, leaving derivative $\mathcal J'(\vartheta)=6\vartheta-12$. Below two this rate is negative, so increasing the prediction a little decreases loss. Above two it is positive. The derivative is zero at the bottom of the bowl. A zero derivative alone does not generally guarantee a minimum; the completed square supplies that guarantee here. The lab checks both the expanded polynomial and the derivative’s finite change formula against the original sum of residual squares.

3.1.2 Turn scalar changes into the normal equations

For a linear fit, define a predictor matrix $E$ with $n$ observation rows and $k$ columns, a response column $y$ with $n$ entries, and a coefficient column $\vartheta$ with $k$ entries. Define prediction column $\hat y=E\vartheta$ and residual column $e=y-\hat y$. Let $j$ index a coefficient. Its residual entry is $e_a=y_a-\sum_{j=1}^k E_{aj}\vartheta_j$. The squared residual length is $\mathcal J(\vartheta)=\sum_a e_a^2$.

Hold all coefficients except $\vartheta_j$ fixed. Increasing that coefficient by $\epsilon$ changes residual $a$ by $-E_{aj}\epsilon$. The minus sign comes from subtracting the prediction. Expanding the changed square gives $e_a^2-2e_aE_{aj}\epsilon+E_{aj}^2\epsilon^2$. Summing, subtracting the old loss, dividing by $\epsilon$, and letting it approach zero gives

$$\frac{\partial\mathcal J}{\partial\vartheta_j} =-2\sum_{a=1}^n E_{aj}e_a.$$

At a minimum, every one of these partial derivatives is zero. Each sum is an entry of $E^{\mathsf T}e$. Substituting $e=y-E\vartheta$ and moving terms gives the normal equations

$$E^{\mathsf T}E\vartheta=E^{\mathsf T}y,\qquad E^{\mathsf T}e=0.$$

The second equality says that the residual is perpendicular to every predictor column. It is a useful independent check. In the fixture, $E$ has rows $(1,0),(1,1),(1,2)$ and the response column is $(1,2,2)^{\mathsf T}$. The two fitted coefficients are $7/6$ and $1/2$, with squared residual $1/6$. Independent predictor columns make this minimizer unique. Dependent columns need a different discussion, which follows in the next lesson.

The first column of this design is all ones; its coefficient is the intercept, the fitted response when the second predictor is zero. The second coefficient is the slope, the prediction change per unit change of that predictor. The products can be calculated before solving:

$$E^{\mathsf T}E=\begin{pmatrix}3&3\\3&5\end{pmatrix}, \qquad E^{\mathsf T}y=\begin{pmatrix}5\\6\end{pmatrix}.$$

Thus $3\vartheta_1+3\vartheta_2=5$ and $3\vartheta_1+5\vartheta_2=6$. Subtracting the first equation from the second gives $2\vartheta_2=1$. Substitution gives $\vartheta_1=(5-3/2)/3=7/6$. The fitted responses are $(7/6,5/3,13/6)^{\mathsf T}$, so the residual is $(-1/6,1/3,-1/6)^{\mathsf T}$. Squaring and adding gives $1/36+1/9+1/36=1/6$. Its length is $\sqrt{1/6}$; its mean squared error is $(1/6)/3=1/18$. The residual’s dot product with the all-ones column is $-1/6+1/3-1/6=0$, and with the second predictor it is $0+1/3-2/6=0$.

Why does perpendicularity minimize rather than merely stop the loss? Let $\hat\vartheta$ be the fitted column and let $c$ be any $k$-entry coefficient change. The new residual is $e-Ec$. Expanding its squared length gives $\lVert e\rVert^2-2c^{\mathsf T}E^{\mathsf T}e+\lVert Ec\rVert^2$. The middle term is zero and the last term cannot be negative. Independent columns make it positive for every nonzero $c$, establishing uniqueness.

3.1.3 Separate a stable prediction from a stable coefficient

Uniqueness does not guarantee stability. Let $t=(-1,0,1)^{\mathsf T}$ and $s=(1,-2,1)^{\mathsf T}$ be two centered columns, and set $\varepsilon=0.001$. Assemble $t$ and $t+\varepsilon s$ as the two columns of a new $3\times2$ predictor matrix $E_\varepsilon$. Its columns are nearly the same. Fitting response $t$ gives coefficients $(1,0)^{\mathsf T}$; changing the response to $t+\varepsilon s$ gives $(0,1)^{\mathsf T}$. A small response change causes a large coefficient change. Their fitted responses can nevertheless remain close.

The Anthropology companion’s ridge lesson already defines a squared-coefficient penalty. Let $I_2$ be the two-coordinate identity matrix. Here a fixed penalty $\rho=0.1$ demonstrates its motivation: repeat both fits with $E_\varepsilon^{\mathsf T}E_\varepsilon+\rho I_2$. These centered predictors have no fitted intercept in this comparison. The coefficient change is much smaller, although penalization changes the objective and introduces bias. Stability is not a substitute for held-out evaluation. Hoerl and Kennard 1970

Exercise. Why does the tiny perturbation not imply tiny coefficient changes? The two predictor columns nearly describe the same direction. Large, opposing coefficient adjustments can nearly cancel in the fitted response. Inspect both the response-change norm and the coefficient-change norm; the notebook computes them separately.

3.2 Recover the shortest solution when a fit is not unique

The word inverse can hide two different tasks: solving for a response that the model can reproduce and choosing among several coefficient vectors that reproduce it. A rank-deficient matrix can leave the first task possible while making the second ambiguous. A pseudoinverse supplies a declared choice. Gene Golub and William Kahan’s 1965 paper connects singular-value computation with this generalized inverse and least-squares applications. Golub and Kahan 1965

Let $M$ be a real matrix with $u$ rows, $v$ columns, and rank $r$. Rank counts independent directions. Its compact singular-value decomposition uses orthonormal-column matrices $U$ of size $u\times r$ and $V_M$ of size $v\times r$, together with a diagonal $r\times r$ matrix $\Sigma$ of positive singular values. Define $\Sigma^{-1}$ by taking the reciprocal of each positive diagonal entry. The pseudoinverse, denoted $M^+$ and sized $v\times u$, is

$$M=U\Sigma V_M^{\mathsf T},\qquad M^+=V_M\Sigma^{-1}U^{\mathsf T}.$$

Only positive singular values are inverted. A thin factorization may contain additional zero singular values; those entries contribute zero to the pseudoinverse rather than an attempted division by zero. In floating-point arithmetic, deciding which small values count as zero requires a scale-aware tolerance. The example uses an exact rank deficiency and clearly separated positive scale, avoiding an ambiguous threshold experiment.

Take $M$ with rows $(1,1)$ and $(2,2)$, and response column $y=(1,2)^{\mathsf T}$. Both $(1,0)^{\mathsf T}$ and $(0,1)^{\mathsf T}$ fit perfectly. So does every coefficient column whose entries sum to one. The nullspace consists of coefficient changes proportional to $(1,-1)^{\mathsf T}$: multiplying such a change by $M$ gives zero. Denote the chosen coefficient column by $\vartheta$. The pseudoinverse gives $\vartheta=M^+y=(1/2,1/2)^{\mathsf T}$. Its squared length is $1/2$, smaller than the squared length one of either endpoint solution.

Now change the response to $(1,3)^{\mathsf T}$. It no longer lies exactly along the single output direction. The shortest least-squares coefficient column is $(0.7,0.7)^{\mathsf T}$, giving fitted response $(1.4,2.8)^{\mathsf T}$ and squared residual $0.2$. The residual remains perpendicular to the columns of $M$. “Shortest” describes the coefficient norm among all least-squares minimizers; it does not mean the smallest prediction error has become zero.

3.2.1 Separate the best prediction from the shortest coefficients

Write the sum of the two coefficients as $s=\vartheta_1+\vartheta_2$. Multiplication by $M$ gives $(s,2s)^{\mathsf T}$: the matrix can see the sum but cannot distinguish how it was divided between coefficients. For the exact response, $s=1$. Parameterize all its solutions as $(1/2+t,1/2-t)^{\mathsf T}$, where $t$ is any real scalar. Their squared length is

$$(1/2+t)^2+(1/2-t)^2=1/2+2t^2.$$

It is smallest at $t=0$. This is the geometric reason for equal coefficients; the pseudoinverse formula performs the same choice in general dimensions.

For the changed response the residual is $(1-s,3-2s)^{\mathsf T}$. Its squared length expands to $5s^2-14s+10=5(s-7/5)^2+1/5$. The best sum is therefore $s=7/5$ and splitting it equally gives each coefficient $7/10$. The actual residual is $(-0.4,0.2)^{\mathsf T}$, with squared length $0.16+0.04=0.2$. No choice of split can reduce that residual; adding a nullspace vector changes coefficients while leaving the fitted response fixed. Conversely, “shortest” does not mean choosing the zero coefficients, since their prediction fails the first requirement of having minimum residual error.

The laboratory constructs the inverse from the factors and checks four identities: $MM^+M=M$, $M^+MM^+=M^+$, and symmetry of both $MM^+$ and $M^+M$. The two products project into different spaces and have different roles. Anthropology’s transpose-based return is a special case with orthonormal loading columns, not a general replacement for this construction.

Exercise. Double the first response column. Linearity doubles the chosen coefficients and quadruples their squared length. The notebook checks both changes. The result illustrates a mathematical convention for underdetermined fitting; it does not identify missing social evidence or reconstruct a discarded person attribute.

3.3 Build an orthogonal frame one subtraction at a time

The existing notebooks use orthonormalization inside larger calculations. This lesson makes its construction visible. Jørgen Pedersen Gram and Erhard Schmidt supply the historical references; Schmidt explicitly credited Gram. We use finite matrices rather than reproducing their function-space settings. Gram 1883; Schmidt 1907

Let $M$ be a matrix with $n$ rows and $k$ linearly independent columns, where $k\leq n$. Its observations are rows, but the algorithm constructs an orthogonal frame for its column space. Call its $j$th column $c_j$. Let $q_i$ denote an already constructed unit column, and let $u_j$ be the part of $c_j$ left after removing its projections onto earlier columns. In exact arithmetic the construction is

$$u_j=c_j-\sum_{i<j}q_i(q_i^{\mathsf T}c_j),\qquad q_j=\frac{u_j}{\lVert u_j\rVert}.$$

Independence ensures a nonzero denominator. Assemble the columns $q_j$ into $Q_{\mathrm{qr}}$, an $n\times k$ orthonormal-column matrix. Define the $k\times k$ matrix $R_{\mathrm{qr}}=Q_{\mathrm{qr}}^{\mathsf T}M$. Its entries below the diagonal are zero because each new input column has no component along a future constructed direction. Thus $M=Q_{\mathrm{qr}}R_{\mathrm{qr}}$. The name QR describes these two factors; it does not identify this factorization with the residual statistic called $Q$ elsewhere.

Take the columns $c_1=(1,1,0)^{\mathsf T}$ and $c_2=(1,0,1)^{\mathsf T}$. The first unit column is $q_1=c_1/\sqrt2$. The coefficient removed from the second is $1/\sqrt2$, leaving $u_2=(1/2,-1/2,1)^{\mathsf T}$, whose squared length is $3/2$. The notebook prints these intermediate quantities, the final factors, and checks both orthonormality and reconstruction. Checking only one of those identities would be insufficient: a set of unit columns might reconstruct the wrong matrix, and a reconstruction could use nonorthogonal columns.

3.3.1 Read the projection as a coefficient and a vector

A dot product with a unit column returns a scalar coordinate along that column. Multiplying the coordinate by the unit column puts that component back in the original space. Here $q_1^{\mathsf T}c_2=(1+0+0)/\sqrt2=1/\sqrt2$, so the vector removed is $(1/2,1/2,0)^{\mathsf T}$. Subtracting it from $(1,0,1)^{\mathsf T}$ gives the stated remainder. Its squared length is $1/4+1/4+1=3/2$, and its length is $\sqrt{3/2}$. Dividing the remainder by this positive length gives $q_2=(1,-1,2)^{\mathsf T}/\sqrt6$.

The dot product $q_1^{\mathsf T}u_2=(1/2-1/2)/\sqrt2=0$ checks that the removed component has not been left behind. The decomposition $c_2=q_1/\sqrt2+q_2\sqrt{3/2}$ explains the second column of $R_{\mathrm{qr}}$: it stores reconstruction coefficients rather than another set of observation rows. The original squared length two splits into projection energy $1/2$ plus remainder energy $3/2$. This is the Pythagorean identity for perpendicular components. It explains why residual energy measures what the selected direction failed to represent.

The executable procedure updates each remainder immediately after a projection. This is the modified Gram-Schmidt ordering. It repeats the projections once to reduce lost orthogonality from rounding. Exact formulas that agree algebraically need not have the same numerical behavior, so the order and the second pass are stated openly. The small tolerance detects a numerically empty remainder relative to the input matrix; it is not a general-purpose rank estimator for arbitrary data scales.

QR also gives a route to least squares without first forming a Gram matrix. For a response column $y$ with $n$ entries, compute the response coordinates $Q_{\mathrm{qr}}^{\mathsf T}y$, then solve the upper-triangular system with $R_{\mathrm{qr}}$. The lab compares this answer with normal equations on the same well-conditioned fixture. Agreement here does not prove equal stability on a difficult matrix.

Exercise. Replace the second column by the first. The remainder vanishes, so there is no second direction to normalize. The code rejects that full-column QR request. A rank-revealing algorithm could return a smaller factorization instead; silently dividing by a tiny remainder would conceal the change in the problem.

3.4 Reflect a column and obtain a QR factorization

Subtracting one projection at a time is not the only way to create zeros. A Householder transformation reflects an entire vector onto a coordinate direction while preserving lengths. Golub and Christian Reinsch’s published SVD procedures use such transformations as part of a larger numerical algorithm. The fixture here teaches the reflection and QR construction, not a complete replacement for their implementation. Golub and Reinsch 1970

Let $x$ be a nonzero column with $n$ entries, and let $e_1$ be the first coordinate column of length $n$. Define $\operatorname{sign}_+(x_1)$ to be one for a nonnegative first entry and minus one otherwise. Set the scalar $\alpha=-\operatorname{sign}_+(x_1)\lVert x\rVert$ and define the column $v=x-\alpha e_1$. Let $I_n$ be the identity and define the reflection matrix $P_{\mathrm{ref}}$ by

$$P_{\mathrm{ref}}=I_n-2\frac{vv^{\mathsf T}}{v^{\mathsf T}v}.$$

The denominator is the squared length of a nonzero column. The sign choice makes the first component of $v$ add magnitudes rather than subtract nearly equal quantities. Because $P_{\mathrm{ref}}^{\mathsf T}=P_{\mathrm{ref}}$ and $P_{\mathrm{ref}}^2=I_n$, the matrix is symmetric and orthogonal. It sends $x$ to $\alpha e_1$. A zero input needs no change, and the implementation returns the identity in that case.

For $x=(4,3,0)^{\mathsf T}$, the norm is five and the reflected column is $(-5,0,0)^{\mathsf T}$. The notebook computes the matrix from the formula and checks the mapped vector, symmetry, and orthogonality. The minus sign is useful numerical bookkeeping, not a negative length. Applying the reflection twice returns the original vector.

The length calculation is $\sqrt{4^2+3^2+0^2}=5$. The sign rule gives $\alpha=-5$, so $v=x+5e_1=(9,3,0)^{\mathsf T}$ and $v^{\mathsf T}v=81+9=90$. Multiplying out the reflection gives

$$P_{\mathrm{ref}}=\begin{pmatrix}-4/5&-3/5&0\\-3/5&4/5&0\\0&0&1\end{pmatrix}.$$

The first output entry is $(-4/5)4+(-3/5)3=-5$; the second is $(-3/5)4+(4/5)3=0$. More generally, split a column into its component along $v$ and a component perpendicular to $v$. The fraction $vv^{\mathsf T}/(v^{\mathsf T}v)$ extracts the first component. Subtracting twice that component reverses it while preserving the perpendicular part. Reversing a component preserves its squared length, which explains the word reflection and the orthogonality check.

To factor a matrix, embed each reflection in a larger identity so that already processed rows remain fixed. Let $M$ be the $4\times3$ matrix whose rows are $(4,1,2),(3,2,0),(0,2,1),(0,0,2)$. Repeated left reflections make entries below the diagonal zero. Their accumulated product is a $4\times4$ orthogonal matrix $Q_{\mathrm{qr}}$, and the resulting $4\times3$ upper-trapezoidal matrix is $R_{\mathrm{qr}}$. They satisfy $M=Q_{\mathrm{qr}}R_{\mathrm{qr}}$. The full orthogonal factor has four columns, unlike the thin factor in the preceding lesson; the last row of the rectangular triangular factor is zero.

The checks verify the zeros, the factorization, and preservation of the Frobenius norm, which is the square root of the sum of squared matrix entries. Multiplication order matters: the sequence used to reduce the matrix is reversed by transposition when reconstructing it. Printing only a triangular-looking matrix would not test this bookkeeping.

Exercise. Use $x=(-4,3,0)^{\mathsf T}$. The chosen first output becomes positive five, with the other entries zero. The lab checks this separately. Both choices represent the same length-preserving construction; neither assigns a stable semantic sign to a news direction.

3.5 Reduce both sides without squaring the difficult scales

The Gram identity is exact: squared singular values of a matrix are eigenvalues of its column-product matrix. That identity does not imply that forming the product is always a good numerical route. Golub and Kahan’s 1965 work and Golub and Reinsch’s 1970 procedures explain the importance of reducing the original matrix directly. Golub and Kahan 1965; Golub and Reinsch 1970

Let $M$ now be a matrix with $u$ rows and $v$ columns, with $u\geq v$. Define $L_{\mathrm{bd}}$ as an orthogonal $u\times u$ matrix and $R_{\mathrm{bd}}$ as an orthogonal $v\times v$ matrix. We seek a $u\times v$ matrix $B_{\mathrm{bd}}$ whose only potentially nonzero entries lie on its main diagonal and the diagonal immediately above it. Such a matrix is called upper bidiagonal. The reduction satisfies

$$B_{\mathrm{bd}}=L_{\mathrm{bd}}^{\mathsf T}MR_{\mathrm{bd}},\qquad M=L_{\mathrm{bd}}B_{\mathrm{bd}}R_{\mathrm{bd}}^{\mathsf T}.$$

First use a left reflection to remove entries below the current diagonal position. Then use a right reflection to remove entries farther right than the next position in that row. Embed each reflection so earlier zeros remain untouched. Repeating this sequence on the preceding lesson’s $4\times3$ matrix yields a bidiagonal result. The notebook checks the zero pattern, both orthogonal factors, exact reconstruction within rounding tolerance, and preserved total squared size. It does not claim that reaching bidiagonal form completes a singular-value decomposition. Diagonalizing that reduced matrix requires an additional algorithm.

Why avoid an unnecessary squaring? Set $\varepsilon=10^{-8}$, and let $M_\varepsilon$ have rows $(1,1)$ and $(0,\varepsilon)$. Define $s_{\max}$ and $s_{\min}$ as its larger and smaller singular values. Define $t_\varepsilon=2+\varepsilon^2$, the sum of their squares. The matrix maps a unit square to a parallelogram of area $\varepsilon$, also the product of its two singular stretches. A stable analytic calculation for this particular fixture is

$$s_{\max}=\sqrt{\frac{t_\varepsilon+\sqrt{t_\varepsilon^2-4\varepsilon^2}}2},\qquad s_{\min}=\frac{\varepsilon}{s_{\max}}.$$

The small value is approximately $7.0711\times10^{-9}$. Computing it from the product relation avoids subtracting nearly equal quantities. This formula is a special two-dimensional reference, not a general direct-SVD implementation. Python additionally compares it with its numerical SVD; OCaml computes the reference natively without using its Gram-based helper as the authority.

The spectral condition number, denoted $\kappa_2$, is the ratio $s_{\max}/s_{\min}$ for this nonsingular matrix. The exact Gram matrix has condition number $\kappa_2^2$. In binary64, the usual 64-bit floating-point format, its entry $1+\varepsilon^2$ rounds to one for this fixture, erasing the small distinction before an eigensolver sees it. A sophisticated eigensolver cannot recover information already lost in its input.

3.5.1 Locate the loss of information

The column-product matrix in this example is

$$M_\varepsilon^{\mathsf T}M_\varepsilon =\begin{pmatrix}1&1\\1&1+\varepsilon^2\end{pmatrix}.$$

Its trace is $2+\varepsilon^2=t_\varepsilon$ and its determinant is $\varepsilon^2$. The determinant of a two-by-two matrix with entries $a,b,c,d$ is $ad-bc$; here it equals the product of the two squared singular values. Solving the resulting quadratic gives the large squared singular value in the displayed formula. Instead of obtaining the small squared value by subtracting that large root from the trace, use the product relation to recover the small positive singular value.

Conditioning is a property of the problem: some input directions are stretched far more than others, so reconstructing their inputs can amplify small errors. Stability is a property of a computational procedure: how much extra error its rounding introduces. The Gram route both squares the ratio of scales and rounds the product entries before solving. The positive input entry $\varepsilon$ still exists in the original matrix even when $1+\varepsilon^2$ has rounded to one in the product. Working from the original matrix gives a direct algorithm a chance to retain that information; it does not promise exact arithmetic.

Exercise. Is the native Gram-based teaching SVD evidence that a direct SVD also loses this direction? No. It is evidence about that particular computational route. The notebook displays the positive analytic small value beside the rounded Gram eigenvalue and keeps those two claims separate. Production solver selection also requires documented precision, scale, and rank tolerances.

3.6 Watch Jacobi rotations remove off-diagonal entries

An eigensolver can be an inspectable sequence of geometric changes rather than one opaque call. Carl Gustav Jacob Jacobi’s 1846 paper is the historical point of reference. The following modern plane-rotation algorithm exposes the mechanism already used inside the native notebook helper. It is not a transcription of the original paper’s notation. Jacobi 1846

3.6.1 Read an angle on the unit circle

An angle measures a turn from a reference direction. Draw a circle of radius one centered at the origin of two perpendicular coordinates; this is the unit circle. Start at its point $(1,0)$ and turn through angle $\theta$, with positive angles measured counterclockwise. The horizontal coordinate of the resulting point is called $\cos\theta$, the cosine. Its vertical coordinate is $\sin\theta$, the sine. Because the point remains one unit from the origin, their squared coordinates satisfy

$$\cos^2\theta+\sin^2\theta=1.$$

A radian measures angle as the traveled arc length divided by the radius. Here $\pi$ is the circle constant, distinct from a probability: a full turn is $2\pi$ radians, a half-turn is $\pi$, and a right angle is $\pi/2$. The numerical functions in both notebooks take angles in radians. The rotation below uses these circle coordinates as matrix entries. The identity above makes its columns unit length, and their opposite signed cross-products cancel to make the columns perpendicular. Its displayed sign convention describes a clockwise turn when applied to a column; the paired basis change below uses that convention consistently.

Let $d$ count coordinates and let $\mathcal C$ be a real symmetric $d\times d$ working matrix. Select two different coordinate indices $p$ and $q$. Write $a=\mathcal C_{pp}$, $b=\mathcal C_{pq}$, and $c=\mathcal C_{qq}$ for the entries of their symmetric two-coordinate block. Let $\theta$ be a rotation angle and define its $2\times2$ block, denoted $R_{\mathrm{plane}}$, as

$$R_{\mathrm{plane}}=\begin{pmatrix}\cos\theta&\sin\theta\\-\sin\theta&\cos\theta\end{pmatrix}.$$

Embed this block in a $d\times d$ identity matrix and call the resulting rotation $R_{\mathrm J}$. Define the next working matrix by $\mathcal C_{\mathrm{next}}=R_{\mathrm J}^{\mathsf T}\mathcal C R_{\mathrm J}$. The new off-diagonal entry in the selected block is $b\cos(2\theta)+(a-c)\sin(2\theta)/2$. The two-argument inverse tangent $\operatorname{atan2}(y,x)$ returns the direction of the point $(x,y)$ in the appropriate quadrant. Choose

$$\theta=\tfrac12\operatorname{atan2}(2b,c-a).$$

This avoids an ambiguous division when $c-a$ is zero. If $b=0$, skip that block. The choice removes the selected off-diagonal entry up to rounding. The trace is the sum of the diagonal entries. An orthogonal similarity transformation preserves eigenvalues, trace, and Frobenius norm.

For $\mathcal C$ with rows $(2,1),(1,2)$, one rotation yields diagonal values one and three. The accumulated rotation columns are eigenvectors, so checking the diagonal alone is not enough: the laboratory also checks the original matrix multiplied by each recovered column against that column multiplied by its eigenvalue.

In this block $a=c=2$ and $b=1$, so the rotation angle is one half of a right angle. Its columns are $(1,-1)^{\mathsf T}/\sqrt2$ and $(1,1)^{\mathsf T}/\sqrt2$. Multiplying the original matrix by the first column gives $(1,-1)^{\mathsf T}/\sqrt2$, while the second gives $3(1,1)^{\mathsf T}/\sqrt2$. This verifies eigenvalues one and three without relying on a printed diagonal.

The two multiplications in $R_{\mathrm J}^{\mathsf T}\mathcal C R_{\mathrm J}$ have different roles: the right multiplication applies the original map to the new basis columns, and the left multiplication expresses those outputs in the new basis. A diagonal result means each new basis direction maps into itself. For the larger matrix a rotation removes only the selected pair; other entries can move. That is why the lab tracks total off-diagonal energy over successive steps and does not mistake one vanished entry for a completed diagonalization.

Next use the symmetric matrix with rows $(4,1,2),(1,3,0),(2,0,2)$. At each step select the largest remaining off-diagonal magnitude. Define off-diagonal energy as the sum of squared entries outside the diagonal. The code records its decrease, accumulates the rotations in the correct order, and stops when the largest off-diagonal entry is small relative to the matrix norm. A maximum iteration count prevents an unnoticed infinite loop. The final trace remains nine. Both languages compare the reconstructed matrix and eigenpair residuals with independent identities.

Repeated eigenvalues explain a limit on interpretation. Let $I_2$ denote the two-coordinate identity. The matrix $2I_2$ remains unchanged under any orthogonal rotation. It has no preferred pair of directions inside that plane. A deterministic numerical sign or ordering convention makes outputs reproducible, but does not create unique mathematical directions.

Exercise. Run the repeated-eigenvalue fixture and rotate its axes. The notebook verifies that the matrix remains the same. This is why an axis can move while the subspace remains stable; the next lesson measures subspaces directly rather than demanding identical columns.

3.7 Compare whole subspaces with principal angles

Two coordinate frames can disagree column by column while describing exactly the same plane. A successful Procrustes fit demonstrates one way to align paired coordinates; principal angles answer a related but different question about the spaces themselves. Åke Björck and Gene Golub studied their numerical computation in 1973. Björck and Golub 1973

Let $d$ count input coordinates and $K$ count retained directions. Let $V_A$ and $V_B$ be two $d\times K$ matrices with orthonormal columns. Their rows must refer to the same input coordinates, in the same units. Let $j$ index the singular values $c_j$ of the $K\times K$ cross-product $V_A^{\mathsf T}V_B$, arranged from largest to smallest. The inverse cosine, written $\arccos$, returns the angle associated with a cosine. Define the principal angles $\theta_j$ as numbers between zero and a right angle satisfying

$$\theta_j=\arccos c_j,\qquad j=1,\ldots,K.$$

Each singular value lies between zero and one in exact arithmetic. The implementation clips tiny floating-point excursions to that interval before taking the inverse cosine, after checking the orthonormal-input contract. Clipping is not permission to feed arbitrary unnormalized matrices into the formula.

Use three input coordinates and retain two directions. Here $\pi$ denotes the circle constant, not a predicted probability. Let $e_1,e_2,e_3$ be the three coordinate unit columns, set the tilt angle $\alpha=\pi/6$, and choose $V_A=[e_1,e_2]$ and $V_B=[e_1,\cos\alpha\,e_2+\sin\alpha\,e_3]$. The planes share their first direction, giving principal angle zero. Their second directions meet at $\alpha$, giving a second angle of thirty degrees. The code computes the cross-product, its singular values, and both angles rather than inserting the expected answer.

Now rotate the two columns inside each plane, using angles $\pi/4$ and $\pi/7$ respectively. The individual columns change, but the principal angles remain the same. Their projectors, defined as $P_A=V_AV_A^{\mathsf T}$ and $P_B=V_BV_B^{\mathsf T}$, also remain unchanged by internal rotations. An independent check relates their squared Frobenius separation to the principal angles:

$$\lVert P_A-P_B\rVert_F^2=2\sum_{j=1}^{K}\sin^2\theta_j.$$

For this fixture the squared separation is $1/2$. The factor two reflects that both directions of disagreement contribute to the projector difference. This check links the angle calculation to an independently constructed matrix quantity.

Here the cross-product is diagonal with entries $1$ and $\cos\alpha$: the shared direction has dot product one with itself and zero with the tilted direction; the second original direction has dot product $\cos\alpha$ with that tilt. At $\alpha=\pi/6$, the latter is $\sqrt3/2$. Taking inverse cosines recovers angles zero and $\pi/6$. A radian measures angle as arc length divided by radius; $\pi$ radians is a half-turn, making $\pi/6$ thirty degrees.

The projector maps any input column to its nearest column in the specified plane. Internal rotations change how coordinates are written inside that plane, but cannot change this projection. In the final two coordinates, the difference $P_A-P_B$ has entries $1/4,-\sqrt3/4,-\sqrt3/4,-1/4$. Squaring all four and adding gives $1/16+3/16+3/16+1/16=1/2$, agreeing with $2\sin^2(\pi/6)=1/2$. The projector calculation is independent of how an SVD routine chooses signs for its columns.

Exercise. Change $\alpha$ first to zero and then to a right angle. The angle pairs become zero/zero and zero/right-angle. The notebook runs both cases. Identical spaces need not have identical displayed axes, while one shared direction does not imply that all retained directions agree.

This calculation compares geometry under a shared feature description. It cannot directly compare a term-loading matrix with an embedding-loading matrix whose row meanings differ. It also says nothing about whether the articles, people, or source populations used in the two fits were comparable. Those prerequisites need evidence before a small angle can support a continuity claim.

3.8 Merge centered summaries without losing between-block variation

An archive can be summarized in blocks, but its variance is not simply the sum of the variances inside those blocks. The distances between their means matter too. Welford’s corrected-sum update and Tony Chan, Golub, and Randall LeVeque’s later analysis provide historical context for algorithms that keep deviations visible. Welford 1962; Chan et al. 1983

For a block of scalar observations, let $n$ be its positive count, $\mu$ its scalar mean, and $\mathcal M$ its corrected sum of squares, the sum of $(x_a-\mu)^2$ over its values $x_a$. This lesson uses unit weights. For two blocks labeled A and B, write their summaries as $(n_A,\mu_A,\mathcal M_A)$ and $(n_B,\mu_B,\mathcal M_B)$. Define $n=n_A+n_B$ and the scalar mean difference $\delta=\mu_B-\mu_A$. Their merged summary is

$$\mu=\mu_A+\delta\frac{n_B}{n},\qquad \mathcal M=\mathcal M_A+\mathcal M_B+\delta^2\frac{n_An_B}{n}.$$

The final term is the variation between the block means. It is nonnegative and vanishes when those means agree. It is not an optional correction for a particular division of the data: omitting it changes the quantity being calculated. Empty blocks are handled by returning the other summary before applying the formula, because an empty block has no defined sample mean.

To see where that term comes from, take an observation in block A and write $x_a-\mu=(x_a-\mu_A)+(\mu_A-\mu)$. Squaring produces the squared within-block deviation, twice the product of the two terms, and the squared mean shift. When summed over block A, the product term vanishes because $\sum_{a\in A}(x_a-\mu_A)=0$. The block therefore contributes $\mathcal M_A+n_A(\mu_A-\mu)^2$. Repeat for block B. Substituting $\mu_A-\mu=-\delta n_B/n$ and $\mu_B-\mu=\delta n_A/n$ reduces the two mean-shift contributions to $\delta^2n_An_B/n$.

This derivation explains what the stored summary must retain: count tells us how heavily a block’s mean matters, the mean locates the block, and the corrected sum records spread around that location. Keeping only a variance would lose the count and location needed to combine blocks.

For the first block $[2,3]$, the count is two, mean $2.5$, and corrected sum $0.5$. For the second block $[7]$, the count is one, mean seven, and corrected sum zero. The mean difference is $4.5$, and the between-block contribution is $13.5$. The merged mean is four and corrected sum fourteen. Dividing by three gives the population-normalized variance; dividing by two gives sample variance seven. The choice of divisor is a reporting convention after the same centered accumulation.

The within-block arithmetic is $(2-2.5)^2+(3-2.5)^2=0.5$ and $(7-7)^2=0$. The additional term is $4.5^2(2\cdot1)/3=13.5$, making the merged sum $0.5+0+13.5=14$. Check directly around the merged mean four: $(2-4)^2+(3-4)^2+(7-4)^2=4+1+9=14$. Thus almost all the variation in this fixture is between blocks, exactly the part a sum of within-block variances would omit. The paired lab checks this direct route against summary merging.

Next add $10^{12}$ to every observation. In exact arithmetic this translation cannot change any deviation from the mean. Binary64 arithmetic can nevertheless lose the answer when it computes a raw sum of squared values and then subtracts count times squared mean. Both terms are enormous compared with fourteen. The lab displays that naive result, then checks a shifted calculation, Welford updates, and centered block merging against the original corrected sum. The exact rounded failure value is implementation-dependent, so it is not treated as a cross-language reference constant.

The existing companion’s raw weighted-moment accumulation teaches a useful algebraic identity and memory strategy. This lesson adds its numerical boundary and a centered alternative. A method that is more stable for a fixed population does not repair changing weights: if adding articles changes a day’s normalized fitting weights, the summaries must reflect the revised weights too.

Exercise. Merge $[1,3]$ and $[0,4]$. Both means are two, so the between-block term is zero; their corrected sums add to ten. The code verifies this case and the empty-block identity. The exercise isolates why the correction is sometimes invisible in a convenient fixture but necessary in the general formula.

3.9 Find the same clusters in a tree and a complete graph

The news companion already shows why a chain of article links can join endpoints that are not themselves similar. A minimum spanning tree supplies a smaller representation of the same threshold connectivity. Florek and colleagues’ early connecting-tree work and Gower and Ross’s later single-linkage result provide the historical setting. The inspected story code uses pairwise comparisons and union-find, which maintains component identities as edges join groups; it does not construct a tree. Florek et al. 1951; Gower and Ross 1969

A graph has vertices and edges; a complete graph contains an edge for every pair of distinct vertices. A path is a sequence of edges meeting at successive vertices. A cycle is a closed path with no repeated intermediate vertex. A connected component consists of all vertices joined by paths. A spanning tree connects all vertices without a cycle, and a minimum spanning tree has the smallest total edge cost among such trees. Let $n$ count vertices, let $a$ and $b$ index them, and let $d_{ab}$ denote a nonnegative symmetric edge cost between distinct vertices. This subscripted cost is unrelated to the input dimension $d$. Let $\eta$ be a distance threshold. The full threshold graph keeps every edge whose cost is at most $\eta$.

The equivalence says that removing from a minimum spanning tree every edge with cost greater than $\eta$ gives exactly the same connected components. It does not say the two graphs have the same edges. A complete graph can retain many redundant routes. If a tree path needed an edge more expensive than an available route joining its two sides, substituting a cheaper crossing edge would contradict the tree’s minimal total cost. This exchange argument explains why the tree preserves the relevant threshold connections.

Use four points at positions $0,1,3,7$ on a line. Costs are absolute differences. The six possible edge costs are $1,2,3,4,6,7$, all distinct. A simple greedy construction sorts the edges and accepts one only if its endpoints are currently in different components. The accepted costs are $1,2,4$, with total seven. The program computes these choices and verifies that exactly three edges connect all four points.

Initially each vertex is its own component. The edge between positions zero and one costs one and joins those two components. The edge between one and three costs two and joins position three to them. The edge between zero and three costs three, but both endpoints already have a path through position one; accepting it would create a cycle. Skip it. The edge between three and seven costs four and brings in the last vertex. All vertices now connect with three edges. The larger remaining costs cannot improve the greedy construction.

A path’s existence, rather than a bound on the distance between the endpoints, determines a component. Rejecting edges inside existing components prevents cycles; accepting edges between components grows connectivity. A tree on $n$ vertices has $n-1$ edges, so its compactness is structural, not an arbitrary stopping threshold.

Then cut both representations at thresholds $0.5,1,2,3,4$. Their component counts are respectively four, three, two, two, and one. Equality includes an edge whose cost is exactly the threshold. At threshold two the first three points form one component, although the endpoints at zero and three are farther apart than two. The tree preserves chaining; it does not repair or conceal it.

For cosine-based article grouping, one can use cost one minus cosine, so a large similarity becomes a small cost. The input must first obey the nonzero-vector cosine convention, and the inequality must be translated consistently. No claim about event identity follows from this equivalence. It is a statement about a declared numerical graph.

Exercise. Why does moving the threshold from two to three leave the components unchanged? The newly admitted direct edge joins vertices already connected through the middle point. The notebook compares the full partitions, not merely their counts, so two different arrangements with the same number of groups could not pass by accident. Tied costs require a deterministic tree-selection policy, but need not make the threshold components ambiguous.

3.10 Fold a new query into a frozen lexical model

Factorizing a document table is not yet a complete retrieval example. A new query must use the training vocabulary, weighting rule, and coordinate map. The vector-space and latent-semantic-analysis papers provide the historical bridge; the precise TF-IDF recipe here is the one disclosed in the earlier companion, not a quotation of an original historical formula. Salton et al. 1975; Deerwester et al. 1990

Use six document rows and seven vocabulary columns, so $n=6$ and $v=7$. The retained terms are troops, border, ceasefire, bank, rate, mortgage, and rise. The six documents, labeled d1 through d6 in order, are “troops shell border town,” “ceasefire troops border,” “bank raises rate,” “rate rise hits mortgage,” “troops ceasefire talks border,” and “bank mortgage rate rise.” Tokens outside the retained vocabulary do not receive columns. These are synthetic examples, not current news.

Let $h$ index retained terms, $F$ be the $n\times v$ count matrix, and $\mathrm{df}_h$ count documents containing term $h$. Write $\ln$ for the natural logarithm. Define the training inverse-document-frequency (IDF) weight $\mathrm{idf}_h=1+\ln((n+1)/(\mathrm{df}_h+1))$. A positive count gets damped frequency one plus its natural logarithm; an absent term gets zero. Multiply these values by the training IDF weights and normalize each nonzero row to obtain the document matrix $A$. The earlier TF-IDF laboratory explains this recipe in detail; here its frozen application is the focus.

Let $f_q$ be the count row of a new query in those same seven columns. Applying the same weighting and normalization gives query row $a_q$. Do not append the query to the training corpus and recompute document frequencies. That would change the fitted representation while pretending to perform only a lookup. The all-out-of-vocabulary query “radio quantum” has a zero row and returns no lexical direction, rather than dividing by zero.

Let $r$ be a positive retained rank no larger than the rank of $A$. Define $U_r$, of size $n\times r$, and $V_A$, of size $v\times r$, as its first $r$ left and right singular columns. Let the $r\times r$ diagonal matrix $\Sigma_r$ contain their positive singular values. Define the $n\times r$ document-score matrix $T_A$ and the $1\times r$ query-score row $t_q$ by

$$T_A=AV_A,\qquad t_q=a_qV_A.$$

Compare the query with each document by cosine in this reduced space, checking for a zero reduced row as well. The lab ranks “ceasefire border” and “bank mortgage” in the original lexical space and at ranks two and three. It displays all scores and resolves display ties by document identifier. Rank choice changes the comparison space; the example makes no claim that one synthetic ordering measures real retrieval quality.

The document scores are $U_r\Sigma_r$, not just $U_r$. Dividing both document and query coordinates by the corresponding singular values creates another metric. The notebook computes that alternative and exhibits changed similarities. It is a definable modeling choice, but should not be silently substituted for the stated projection.

3.10.1 Follow a query through the same operations as a document

For “ceasefire border,” set the ceasefire and border count entries to one and all others to zero. Damped frequency leaves each positive count at $1+\ln1=1$. Multiplication by each term’s stored IDF weight creates an unnormalized weighted row. Its length is the square root of the sum of those two squared weights. Divide every entry by that one length to obtain $a_q$. Normalization changes the row’s scale without changing which terms it contains.

The $j$th reduced coordinate is the dot product of $a_q$ with column $j$ of $V_A$. Each direction can involve many vocabulary entries, so this step is a projection rather than selecting $r$ individual words. For a nonzero document row $t_a$ of $T_A$, its reduced cosine with the query is $t_a t_q^{\mathsf T}/(\lVert t_a\rVert\lVert t_q\rVert)$. Both lengths belong to reduced coordinates: normalizing the original lexical row does not guarantee that its projection still has length one.

Why is the document table $U_r\Sigma_r$? Multiply the retained approximation $U_r\Sigma_rV_A^{\mathsf T}$ on the right by $V_A$ and use $V_A^{\mathsf T}V_A=I_r$. The right factor cancels, leaving $U_r\Sigma_r$. Dividing by singular values reverses those coordinate scales; in a cosine calculation this generally reweights directions unequally and can change rankings. The choice of ruler must be part of the retrieval specification.

Exercise. Repeat “bank” in the query “bank bank mortgage.” The positive count becomes two, so logarithmic damping is now exercised rather than receiving only zero-or-one counts. IDF stays fixed. The lab recomputes the row and confirms that the transformation changes while the training document matrix remains untouched. Frozen preprocessing makes a lookup interpretable; it does not make vocabulary omissions harmless.

4 Uncertainty, meaning, and learned representations

These lessons add a second kind of question to the geometry of observations. What does a probability promise? What follows from a recorded relation? How can a representation be learned? Small, fully specified examples let us inspect the arithmetic without treating a synthetic demonstration as evidence about the deployed newspapers. Every vector remains a row. Symbols for local neural computations carry explicit qualifiers so they do not replace the news and people notation.

4.1 Uncertainty and the question worth asking

Entropy describes uncertainty over specified possibilities. It does not measure the truth or social value of a news item. In Shannon’s communication theory, this separation lets us reason about distributions of messages without first resolving their meanings. The same mathematical distinction matters when an ontology interface proposes a question. Shannon 1948

Let $k$ count possible states, indexed locally by $i$. Let $\pi_i$ be the probability of state $i$, with nonnegative entries adding to one. Write $\boldsymbol\pi$ for their $1\times k$ row. Define entropy $\mathcal H$ in bits, where $\log_2$ is the base-two logarithm, and assign zero to the limiting expression $0\log_2 0$:

$$\mathcal H(\boldsymbol\pi)=-\sum_{i=1}^k\pi_i\log_2\pi_i.$$

The distributions $(0.5,0.5)$ and $(0.9,0.1)$ have entropies $1$ and approximately $0.468996$ bits. A known state, represented by $(1,0)$, has zero entropy. The outcome has become predictable, not necessarily desirable.

The logarithm asks what power produces a number. Since $2^{-1}=0.5$, $-\log_2(0.5)=1$: observing either equally likely outcome supplies one bit of information. Entropy averages this information over outcomes, weighting each by its probability. Thus the first calculation is $0.5\cdot1+0.5\cdot1=1$. For the unequal probabilities it is $-0.9\log_2(0.9)-0.1\log_2(0.1)$. The rare outcome is more informative when it occurs, but contributes only one tenth of the average. This distinction between one outcome and the expectation will also matter when comparing questions.

The zero convention follows a limit: as a positive probability becomes smaller, its negative logarithm grows, but probability times that logarithm tends to zero. An impossible outcome contributes no expected information. We do not ask the computer to take a logarithm of zero and multiply the resulting infinity by zero; the lab omits zero-probability terms explicitly.

Now define a candidate question with answers indexed by $b$. Let $\ell_{bi}$ be the probability of answer $b$ if state $i$ is true. Each state’s answer probabilities sum to one. The probability of answer $b$ is $\eta_b=\sum_i\pi_i\ell_{bi}$. If $\eta_b>0$, Bayes’ rule gives posterior state probabilities $\pi_{i\mid b}=\pi_i\ell_{bi}/\eta_b$. A posterior is the distribution after observing an answer under this specified likelihood model. Define the question’s expected information gain $\mathcal I$ by

$$\mathcal I=\mathcal H(\boldsymbol\pi)-\sum_{b:\eta_b>0}\eta_b\mathcal H(\boldsymbol\pi_{\mid b}).$$

Four equally likely concepts have entropy two bits. A deterministic question separating two concepts from two leaves one bit under either answer, gaining one bit. Separating one from three sometimes identifies the state completely, but its expected gain is only $2-\tfrac34\log_2 3\approx0.811278$ bits. The rare perfect answer cannot stand in for an average over answers.

For a noisy question, let the probabilities of “yes” in the four states be $(0.9,0.8,0.2,0.1)$. “Yes” has probability $0.5$ and gives posterior $(0.45,0.4,0.1,0.05)$. The opposite answer gives the reversed row. The notebook computes their entropy and the smaller expected gain. Answers with probability zero have no defined posterior and contribute no mass to the expectation.

To reproduce the update, multiply each prior probability $1/4$ by its conditional “yes” probability. The joint masses are $(0.225,0.2,0.05,0.025)$, summing to $0.5$. Dividing every mass by $0.5$ renormalizes the “yes” cases to a total of one and gives the stated posterior. This is Bayes’ rule as multiplication followed by normalization, rather than a change made because one answer sounds persuasive. For “no,” subtract each conditional “yes” probability from one before repeating the calculation.

Information gain is calculated before we know which answer will occur. Weight the entropy of each posterior by that answer’s probability, add, and subtract from prior entropy. For the deterministic one-versus-three split, the two answer probabilities are $1/4$ and $3/4$. Their posterior entropies are zero and $\log_2 3$, giving expected remaining uncertainty $(1/4)0+(3/4)\log_2 3$. A question that sometimes identifies the state can therefore be worse on average than a question that always narrows it evenly.

Exercise. Would a question whose answer is always “yes” reduce uncertainty? Explanation. Its posterior equals its prior, so the gain is zero. A question’s usefulness depends on its likelihood model, not its wording alone. The ontology dependency contains information-based selection machinery, while the audited Anthropology picker uses text matching. This exercise establishes neither deployment nor reliable answer probabilities for a future interface.

4.2 Probability forecasts and honest expected scores

A predicted probability becomes accountable when its event is specified and outcomes are observed. For example, the event might be that a reviewer, following a stated protocol, marks an article as substantially discussing databases. It is not an intrinsic probability that a person “is a database person.” A score summarizes forecasts against such outcomes; it does not repair an ambiguous target.

Let $n$ be the number of reviewed cases, $a$ index a case, $\pi_a\in[0,1]$ its forecast probability, and $o_a\in\{0,1\}$ its observed outcome. We use the modern one-category binary Brier score $\mathrm{BS}$, where lower is better:

$$\mathrm{BS}=\frac1n\sum_{a=1}^n(\pi_a-o_a)^2.$$

Forecast $0.8$ contributes $0.04$ when the event occurs and $0.64$ when it does not. For one of each outcome the mean is $0.34$. Brier’s original multicategory convention sums errors over mutually exclusive categories. In a binary problem, including the complementary category doubles this value to $0.68$. Neither arithmetic is wrong; comparisons require the same convention. Brier 1950

Another score is binary logarithmic loss. Define $\mathrm{LogLoss}$ using the natural logarithm $\ln$, with probabilities strictly between zero and one in this finite calculation:

$$\mathrm{LogLoss}=-\frac1n\sum_{a=1}^n\{o_a\ln\pi_a+(1-o_a)\ln(1-\pi_a)\}.$$

The corresponding contributions for forecast $0.8$ are approximately $0.223144$ and $1.609438$ nats. A nat is the information unit obtained with natural logarithms. A zero probability assigned to an event that occurs incurs infinite log loss. Silently clipping probabilities changes the evaluated rule; a numerical implementation must declare any such convention.

To explain strict propriety, imagine a single event whose stipulated true occurrence probability is $\pi_\star=0.7$. Let $o$ be the binary random outcome, $\pi$ the reported forecast, and $\mathbb E$ mean expectation over the two possible outcomes under that truth. Expanding the squared error gives

$$\mathbb E[(\pi-o)^2]=(\pi-\pi_\star)^2+\pi_\star(1-\pi_\star).$$

Only the first term depends on the report. Its unique minimum is at the truthful report, with expected score $0.21$. Expected log loss also has its unique minimum there, approximately $0.610864$ nats. The executable sweep checks both using reports from $0.01$ to $0.99$. Strict propriety is a statement about expected optimization, not a guarantee that an estimated model is calibrated. Gneiting and Raftery 2007

The expectation is a weighted average over the two outcomes. With probability $\pi_\star$ the event occurs and loss is $(\pi-1)^2$; with probability $1-\pi_\star$ it does not and loss is $\pi^2$. Expanding gives

$$\begin{aligned} \pi_\star(\pi-1)^2+(1-\pi_\star)\pi^2 &=\pi^2-2\pi\pi_\star+\pi_\star\\ &=(\pi-\pi_\star)^2+\pi_\star-\pi_\star^2. \end{aligned}$$

The final term is unavoidable variability even under the correct probability. At $\pi_\star=0.7$ it is $0.7\cdot0.3=0.21$. A truthful probabilistic forecast does not predict each outcome with certainty; its optimum expected loss need not be zero.

Before applying a derivative to a logarithm, we need its local rate rule. For a positive scalar input $x$, the derivative of $\ln x$ is $1/x$. Two other rules used later are that the derivative of $\exp x$ is $\exp x$, and the derivative of $1/x$ is $-1/x^2$ when $x$ is nonzero. Here $\exp x=e^x$ and $\ln$ is its inverse: $\ln(\exp x)=x$. These rules describe limiting rates, in the same sense as the square derivative in the fitting lesson. Constants multiply derivatives, and the derivative of a sum is the sum of its terms’ derivatives.

For a composition, the chain rule multiplies the rate of the outer function by the rate of its inner input. In particular, the input $1-\pi$ changes at rate minus one when $\pi$ increases. Thus $\ln(1-\pi)$ has derivative $-1/(1-\pi)$, while $\ln\pi$ has derivative $1/\pi$. The minus signs in log loss must be applied after these inner changes are accounted for. The network lesson develops this same multiplication along a longer path.

For log loss, differentiating the expected expression $-\pi_\star\ln\pi-(1-\pi_\star)\ln(1-\pi)$ gives $-\pi_\star/\pi+(1-\pi_\star)/(1-\pi)$. Putting this over the positive denominator $\pi(1-\pi)$ leaves numerator $\pi-\pi_\star$. The derivative is negative below the truth and positive above it, explaining the unique interior minimum checked by the sweep.

Exercise. Does one surprising outcome prove that a forecast of $0.8$ was dishonest? Explanation. No. The remaining outcome still has probability $0.2$. Evaluate repeated forecasts under a defined sampling and review protocol. The proper-score principle explains a target for evaluation; it does not manufacture held-out data or justify displaying an arbitrary vector score as a percent confidence.

4.3 Calibration can be perfect while resolution is absent

Calibration asks whether outcomes occur at the frequencies forecasts describe. Resolution asks whether forecasts distinguish groups with different event frequencies. These properties can pull apart. A forecaster that always predicts the overall base rate can be calibrated while telling us nothing about which particular case is more likely to be positive. Murphy’s decomposition makes this distinction quantitative. Murphy 1973

Use the binary outcomes $o_a$ and probabilities $\pi_a$ defined in the preceding lesson. Partition the $n$ cases into groups indexed locally by $b$, each containing exactly the same forecast $f_b$. Let $n_b$ count the group, $\omega_b=n_b/n$ its fraction of all cases, $\bar o_b$ its observed positive fraction, and $\bar o$ the overall positive fraction. These group weights are evaluation frequencies, not the article-fitting weights used for news covariance.

Define reliability error $\mathrm{REL}$, resolution $\mathrm{RES}$, and uncertainty $\mathrm{UNC}$ by

$$\begin{aligned} \mathrm{REL}&=\sum_b\omega_b(f_b-\bar o_b)^2,\\ \mathrm{RES}&=\sum_b\omega_b(\bar o_b-\bar o)^2,\\ \mathrm{UNC}&=\bar o(1-\bar o). \end{aligned}$$

Reliability error is smaller when group forecasts match their outcome frequencies. Resolution is larger when group frequencies differ from the base rate. For groups of identical forecasts the binary Brier score satisfies the exact finite-sample identity

$$\mathrm{BS}=\mathrm{REL}-\mathrm{RES}+\mathrm{UNC}.$$

Consider ten cases in two groups of five. The first outcomes are $(0,0,0,0,1)$, with frequency $0.2$; the second are $(0,1,1,1,1)$, with frequency $0.8$. Overall prevalence is $0.5$, so uncertainty is $0.25$. A forecaster reporting $0.2$ and $0.8$ for these respective groups has reliability error zero and resolution $0.09$. Its Brier score is $0-0.09+0.25=0.16$.

Keep the outcomes fixed but change the reports to $0.4$ and $0.6$. Resolution stays $0.09$, because the two forecast groups are unchanged. Reliability error becomes $0.04$, increasing the Brier score to $0.20$. Finally, predict $0.5$ for every case. There is now one forecast group: reliability error is zero, resolution is zero, and the score is $0.25$. Perfect aggregate reliability has not supplied discrimination.

4.3.1 Derive the decomposition from within-group deviations

Inside group $b$, write $f_b-o_a=(f_b-\bar o_b)+(\bar o_b-o_a)$. When we square and average over that group, the cross term vanishes because the average of $\bar o_b-o_a$ is zero. The group’s mean squared error is therefore $(f_b-\bar o_b)^2+\bar o_b(1-\bar o_b)$. The second term follows because binary outcomes satisfy $o_a^2=o_a$.

Average over groups with weights $\omega_b$. The first terms give reliability error. For the remaining terms, expand the definition of resolution to obtain $\mathrm{RES}=\sum_b\omega_b\bar o_b^2-\bar o^2$, using $\sum_b\omega_b\bar o_b=\bar o$. Consequently, $\sum_b\omega_b\bar o_b(1-\bar o_b)=\mathrm{UNC}-\mathrm{RES}$. This proves the identity without assigning a separate empirical meaning to each case’s squared error.

For the well-calibrated fixture, each group contributes weight $1/2$. Resolution is $(1/2)(0.2-0.5)^2+(1/2)(0.8-0.5)^2=0.09$. Within the first group the direct score is $[4(0.2-0)^2+(0.2-1)^2]/5=(0.16+0.64)/5=0.16$; the second group has the same value. These separate arithmetic routes make the identity inspectable. Merging all reports to $0.5$ also merges the groups, which is why the resolution then disappears.

Exercise. Why may we not substitute arbitrary wide probability bins into this identity without qualification? Explanation. Forecasts that vary inside a bin are no longer identical. Replacing them with one bin average can discard terms. Our code groups equal forecasts exactly, then verifies the identity directly against every case’s squared error. Reliability diagrams remain useful, but their binning, sample counts, and uncertainty must be disclosed. These fictional ten cases are an arithmetic fixture, not evidence that an Anthropology concept classifier has achieved any particular quality.

4.4 An incidence table can generate a concept lattice

An incidence table records whether an object has an attribute. It can support a mathematics of exact shared properties, quite different from the covariance of numerical profiles. Formal concept analysis builds paired sets from such a table. Its word “concept” has a precise algebraic meaning here; it does not mean an eigenvector with an attractive label. Ganter and Wille 1999

Let $\mathcal O$ be a finite set of objects and $\mathcal A$ a finite set of attributes. These local sets are not article-background cells. Let $I_{\mathrm{FCA}}$ be their binary incidence table, with one row per object and one column per attribute: entry $I_{\mathrm{FCA},ij}$ is one when object $i$ has attribute $j$. We deliberately do not call this table $J$, which already denotes person–article association. For an object subset $E\subseteq\mathcal O$, define $E'$ to be the attributes shared by every object in $E$. For an attribute subset $F\subseteq\mathcal A$, define $F'$ to be all objects possessing every attribute in $F$. The prime denotes these two derivation operations, distinguished by their input sets:

$$E'=\{j\in\mathcal A: I_{\mathrm{FCA},ij}=1\text{ for every }i\in E\},$$ $$F'=\{i\in\mathcal O: I_{\mathrm{FCA},ij}=1\text{ for every }j\in F\}.$$

Applying the derivations twice gives the object closure $E''$. It contains the original objects and is unchanged by another closure. A formal concept is a pair $(E,F)$ satisfying both $E'=F$ and $F'=E$. Its object set $E$ is called its extent, and its attribute set $F$ its intent.

Use four fictional objects, numbered zero through three, with attributes respectively {database, distributed}, {database}, {distributed, open}, and {database, distributed, open}. Object zero’s shared attributes are database and distributed. Both objects zero and three have these attributes, so $\{0\}''=\{0,3\}$. Empty inputs require care: every attribute is shared by all objects in an empty object set, because no counterexample exists. Thus the closure of the empty extent here is $\{3\}$, not the empty set.

The code enumerates all sixteen object subsets and finds six formal concepts. Order these concepts by extent inclusion. For two closed extents $E_1,E_2$, their greatest lower bound is the largest closed extent contained in both, their intersection. Their least upper bound is the smallest closed extent containing both, the closure of their union. These operations are called meet and join. The notebook verifies the defining lower- and upper-bound conditions against every enumerated concept, rather than trusting the names.

4.4.1 Follow both directions of the set operation

Begin with the object set $\{0\}$. Read across its row to collect database and distributed. Then read down those two attribute columns and keep only objects having both: objects zero and three. Reading across both of these rows again returns exactly database and distributed. We have reached the formal concept $(\{0,3\},\{\text{database},\text{distributed}\})$. Closure adds objects indistinguishable under the shared attributes; it does not invent a new recorded attribute or turn a similarity score into a fact.

“Every” matters in the definition. Selecting more objects can only reduce or preserve their shared attribute set. Conversely, requiring more attributes can only reduce or preserve the matching object set. Applying both operations therefore gives an object closure that contains the starting set, respects set inclusion, and no longer changes after another closure. These properties are called extensiveness, monotonicity, and idempotence, respectively.

A lower bound under extent inclusion is a concept whose extent is contained in both selected extents; the greatest such bound contains every other lower bound. An upper bound contains both, and the least upper bound is contained in every other upper bound. Raw union can fail to be closed, so the join needs closure after union. A lattice is an ordered collection in which every pair has these two bounds. The lab checks this universal condition over the finite list, including the empty-input convention rather than silently omitting it.

Exercise. Why does the most specific concept contain object three? Explanation. It has all three recorded attributes. Under extent inclusion, that singleton is below every other closed extent in this fixture. Real incomplete evidence may not justify treating an unrecorded attribute as absent. Our binary table is stipulated complete. Anthropology’s supported numerical people profiles neither supply that completeness nor imply that formal-concept inference is deployed.

4.5 A recorded triple and a licensed inference are different

An identifier, a text label, and a rule have different jobs. RDF describes statements as subject–predicate–object triples. The subject identifies the thing discussed; the predicate identifies a relation; the object is another identified thing or a literal value. A literal can be a number or text. Recording a tuple is easy. Deciding what follows from it requires declared semantics. Cyganiak et al. 2014

Consider a finite set $\mathcal T$ of triples. A triple $(s,r,o)$ has subject $s$, relation $r$, and object $o$; these letters are local tuple positions, not story scores or matrix rank. We will use exactly two inference rules. The relation subclass means class inclusion in this toy rule system, and type assigns an individual to a class. For arbitrary identifiers $A,B,C$ and individual $x$, the rules are

$$(A,\mathrm{subclass},B),(B,\mathrm{subclass},C) \Longrightarrow(A,\mathrm{subclass},C),$$ $$(x,\mathrm{type},A),(A,\mathrm{subclass},B) \Longrightarrow(x,\mathrm{type},B).$$

A premise is a statement required by a rule; a conclusion is the statement that rule licenses when its premises match. Provenance records where an assertion came from or which premises supported a conclusion. The arrow means that both premises license the displayed conclusion. These capital letters name classes locally; they are not the term table or concept count elsewhere in the book. Starting with asserted triples, repeatedly add rule conclusions until no new tuple appears. This fixed point is the closure under these two rules. Because our set of identifiers is finite and the procedure only adds triples, it terminates.

The fixture asserts five statements: Ada has type CEO; CEO is a subclass of Person; Person is a subclass of Agent; Databases has the broader concept Technology; and Ada has topic Databases. Closure adds three statements: CEO is a subclass of Agent, Ada has type Person, and Ada has type Agent. Each derived record retains its two immediate supporting premises. Tracing these links reaches the initial assertions.

Read the first rule as a pattern with a shared middle identifier. The premises CEO–subclass–Person and Person–subclass–Agent match because Person occupies that middle position, allowing CEO–subclass–Agent. The second rule matches Ada–type–CEO with CEO–subclass–Person to license Ada–type–Person. Applying it again with Person–subclass–Agent licenses Ada–type–Agent. Alternatively, the newly derived CEO–subclass–Agent provides another route to the same type conclusion. Adding an already present triple does not add new knowledge, so the implementation uses set membership rather than accumulating duplicate rows forever.

Repeated rule application is necessary because one conclusion can become another rule’s premise. A fixed point means a further full pass adds nothing; it does not mean that all possible facts about the represented world are known.

Nothing in either rule turns broader into subclass, or hasTopic into type. The code therefore does not derive that Databases is a subclass of Technology or that Ada has type Databases. SKOS concept relations and OWL class axioms serve different purposes; their shared graph-shaped appearance is not permission to interchange them. This limited example is not an OWL implementation or a SKOS conformance test. Miles and Bechhofer 2009; W3C OWL Working Group 2012

Exercise. If no statement about Bo appears, should the query “Bo has type Person” return false? Explanation. Under our open-world interpretation, it returns unknown. Lack of an assertion or derivation is not evidence of its negation. “Open world” is an interpretation of missing knowledge, not a guarantee that supplied assertions are correct. Provenance still matters: a valid derivation from an erroneous premise remains erroneous about the world. The newspapers’ numerical similarities cannot substitute for the premises of a sourced employment or financing relationship.

4.6 A hidden representation learns through the chain rule

An embedding is useful partly because its coordinates can be learned rather than individually named in advance. Backpropagation efficiently applies the chain rule through a composed computation, allowing an error at the output to update parameters earlier in the network. The historical contribution of learned hidden representations should not be retold as the invention of differentiation itself. Rumelhart et al. 1986

Begin with one scalar input $x$, equivalent to a one-coordinate observation row, and two scalar parameters $\vartheta_1,\vartheta_2$. Define preactivation $u=\vartheta_1x$, hidden value $h=\tanh u$, and prediction $\hat o=\vartheta_2h$. Here $\tanh$ is the hyperbolic tangent, a smooth nonlinear function with derivative $1-\tanh^2u$. Let $o$ be a numerical target and let $\mathcal L_{\mathrm{net}}$ be half the squared prediction error:

$$\begin{aligned} u&=\vartheta_1x, & h&=\tanh u,\\ \hat o&=\vartheta_2h, & \mathcal L_{\mathrm{net}}&=\tfrac12(\hat o-o)^2. \end{aligned}$$

The hidden value is an intermediate representation. Calling it hidden means it lies between input and output; it does not imply secrecy or a named semantic property. Denote output error by $e_{\mathrm{net}}=\hat o-o$. The local symbol is distinct from the exponential constant.

4.6.1 Build the derivative from small local changes

The hyperbolic tangent can be written $\tanh u=(\exp u-\exp(-u))/(\exp u+\exp(-u))$. Its output lies between minus one and one. Near zero it responds strongly to small changes; near either extreme its derivative $1-h^2$ becomes small. The nonlinearity matters: without it, multiplying by two scalar parameters would just be one linear coefficient $\vartheta_2\vartheta_1$.

First consider loss as a function of prediction alone. Increase prediction by $\epsilon$ while keeping the target fixed. The loss change is $\tfrac12[(e_{\mathrm{net}}+\epsilon)^2-e_{\mathrm{net}}^2] =e_{\mathrm{net}}\epsilon+\epsilon^2/2$. Divide by $\epsilon$ and let it tend to zero. The local loss derivative is $e_{\mathrm{net}}$. The factor one half in the loss cancels the factor two from differentiating a square; it does not change the best-fitting parameters.

Next keep $h$ fixed and change $\vartheta_2$. Prediction changes at rate $h$, so the loss changes at rate $e_{\mathrm{net}}h$. Changing $\vartheta_1$ takes a longer route: preactivation changes at rate $x$, hidden value changes at rate $1-h^2$ with respect to preactivation, prediction changes at rate $\vartheta_2$ with respect to hidden value, and loss changes at rate $e_{\mathrm{net}}$ with respect to prediction. Multiplying these local rates is the chain rule. Intermediate changes shrink along with the parameter change, so it is the first-order rates that multiply. In a network with several paths from one parameter to the loss, the contributions of those paths add. The resulting parameter derivatives are

$$\frac{\partial\mathcal L_{\mathrm{net}}}{\partial\vartheta_2}=e_{\mathrm{net}}h,\qquad \frac{\partial\mathcal L_{\mathrm{net}}}{\partial\vartheta_1} =e_{\mathrm{net}}\vartheta_2(1-h^2)x.$$

Each factor in the second expression describes one dependency: output loss on prediction, prediction on hidden value, hidden value on preactivation, and preactivation on its parameter. Computing these local quantities once and reusing them is the central mechanism behind reverse accumulation in larger networks.

Set $x=2$, $o=0.3$, $\vartheta_1=0.4$, and $\vartheta_2=-0.7$. The forward pass gives $h\approx0.664037$, prediction $-0.464826$, and loss $0.292479$. The two derivatives are approximately $0.598612$ and $-0.507872$. A gradient-descent update subtracts a positive step size times each derivative. With step size $0.1$, the new loss is approximately $0.234023$.

The forward arithmetic begins with $u=0.4\cdot2=0.8$, then evaluates $h=\tanh(0.8)$, multiplies by $-0.7$ for the prediction, and subtracts target $0.3$ to obtain $e_{\mathrm{net}}\approx-0.764826$. The local rates along the first-parameter path are this error, $-0.7$, $1-h^2\approx0.559055$, and $2$. Their product is positive: both the output error and output weight are negative. Subtracting its step therefore decreases $\vartheta_1$. The derivative for $\vartheta_2$ is negative, so that parameter increases. The new parameters are approximately $(0.340139,-0.649213)$. Recompute the entire forward path at the new parameters; retaining the old hidden value would not evaluate the updated network.

A derivative is a local prediction of change, not the full change over a finite step. The finite-difference check perturbs only one parameter at a time, runs the forward computation at both nearby settings, and divides the loss difference by the total parameter separation. Agreement with the analytic formula checks that dependencies and signs were handled correctly. It does not turn a chosen step size into a proof of improvement for every possible input.

The notebook also approximates each derivative by central finite differences. For a parameter displacement $\epsilon>0$, this uses the loss difference at plus and minus $\epsilon$, divided by $2\epsilon$. Choosing $\epsilon=10^{-6}$ checks our formulas without treating finite differences as the intended training algorithm.

Exercise. Must an arbitrarily large descent step reduce the loss? Explanation. No. The derivative describes local behavior; a long step can overshoot or enter a differently curved region. Even a correctly decreasing training loss does not establish useful language representations. Our network has one invented example and no evaluation population. It teaches a differentiable dependency chain, not the training of the deployed embedding checkpoint.

4.7 A word vector learns by contrasting contexts

A representation can be trained by rewarding observed word–context pairs while contrasting them with sampled alternatives. Negative sampling made this idea computationally useful in the word2vec lineage. It avoids treating every vocabulary entry as an explicitly normalized candidate at every update. Its objective is a training construction; a score between two learned rows is not a documented social relationship. Mikolov et al. 2013a; Mikolov et al. 2013b

Let $e_{\mathrm{ctr}}$ be a $1\times d_{\mathrm{emb}}$ center-word row. The local dimension $d_{\mathrm{emb}}$ counts this toy embedding’s coordinates. Let $e_+$ be an observed context row and $e_-$ one sampled negative context row of the same shape. Define positive and negative pair scores $u_+=e_{\mathrm{ctr}}e_+^{\mathsf T}$ and $u_-=e_{\mathrm{ctr}}e_-^{\mathsf T}$. The logistic function is $\operatorname{sigmoid}(u)=1/(1+\exp(-u))$. For this one-positive, one-negative example, define loss