Preface

This is the shared mathematical companion for Eigen Times and Eigen Hacks. It accompanies Eigen Times: A Newspaper in the Eigenbasis of the News, which states the methods compactly and reports their results. Here we slow down. A reader needs arithmetic, fractions, and a willingness to work through small lists of numbers. We explain the other tools when they become necessary: vectors, matrices, logarithms, covariance, eigenvectors, derivatives, and the singular value decomposition. The aim is to understand what each calculation measures, why we want it, and how to reproduce it.

The same mathematics can organize newspaper articles and technology discussions. Their data distributions need not be the same. A similarity threshold measured on one archive does not automatically work on another. Historical corpus counts and threshold observations in this book retain their original dates and source context; they are not current site counters or new measurements of Eigen Hacks. Synthetic examples, by contrast, are small teaching inputs that the reader can change.

A chapter first explains the problem in words, introduces the necessary terms, and then works through a calculation. Read the intermediate arithmetic before the compact formula. When a formula still looks dense, identify what goes in, what comes out, and the size of each object. The upfront guide and final glossary provide a second route through the terminology; the subject index links back to explanations.

Every matrix in the original figures is drawn as a grid of coloured squares: blue for positive entries, rust for negative, darker for larger. The figure inputs and dated aggregate observations remain frozen. The complete Python and native OCaml notebooks contain this text, all nineteen figures, and executable lessons interleaved with the chapters. Their cells calculate the worked examples from declared inputs, test identities and reference answers, and finish with matching numerical receipts. Neither notebook needs access to the production database or a private archive. The public teaching repository is querygraph/eigenmath; a delivered study bundle also contains the exact notebook editions and reading copies.

Version 1.1.2. This patch expands compressed teaching steps throughout the shared companion. It introduces the meaning of a residual before measuring its length, develops differentiation from a small change in a scalar function before using it for minimization, and unrolls normalization, covariance, decomposition, rotation, statistical thresholds, and streaming sums. It adds an upfront notation guide, a glossary, and a linked index, and repairs clipped or overlapping figure annotations. One clustering sentence is corrected to agree with its existing frozen fixture: a threshold of 0.9 leaves one linked pair on 24 February 2022. No archive is refitted and no historical threshold or measurement is replaced.

Reading and notation guide

0.1 How to work through the book

Chapters 1–3 build the numerical language. Chapters 4–6 find recurring directions and give them interpretable names. Chapter 7 measures a new article or story. Chapters 8–10 explain grouping, unusual days, and how directions retain their identity over time. Chapters 11–12 show how the same calculations can be accumulated in pieces rather than loading an archive into memory.

In a notebook, run the setup cell first and then run cells in reading order. A cell may use a result computed earlier. An assertion is an executable check: it stops the lesson if a mathematical identity or reference result fails. If you edit a teaching input, a reference answer may intentionally fail while an identity should still hold. Restarting the kernel—the running language process—and running all cells checks that the lesson does not depend on hidden state. Python and OCaml are alternative implementations of the same examples; neither language is a prerequisite for reading the prose.

The later volume Eigen Times Math History places these methods in a historical sequence. Its companion adds techniques beyond the present book; its complete teaching text is available in the History Math notebooks. The present chapters are the direct prerequisites for the calculations here; the historical volume is a route to the foundational literature, not a substitute for a missing teaching step.

0.2 Reading symbols before calculating

A scalar is a single number. A vector is an ordered list of numbers. A matrix is a rectangular table. An entry is one number inside a vector or matrix; a coordinate describes an amount along a chosen direction. Coordinates depend on which directions we choose.

Lowercase subscripts identify entries: $x_i$ is observation number $i$, while $x_{ij}$ is its coordinate number $j$. Context distinguishes an observation index from a coordinate index. An index is a label, not a quantity to multiply. Superscripts usually indicate powers, so $a^2=a\,a$, except that $\top$ means transpose: exchange rows and columns. An accent is also part of the name: $\hat{x}$ means an estimate or reconstruction of $x$, not a new multiplication.

For numbers $a_1,\ldots,a_n$, the dots mean the intervening entries. The expression $\sum_{i=1}^{n}a_i$ means $a_1+\cdots+a_n$. The symbol $\prod$ means multiplication over an indicated set of entries. Parentheses group operations: first calculate what is inside. The square root $\sqrt{a}$ is the nonnegative number whose square is $a$, for $a\geq0$. The relations $\leq$ and $\geq$ allow equality; $<$ and $>$ do not. The relation $\approx$ means approximate equality, usually because a displayed decimal is rounded.

The real numbers, denoted $\mathbb{R}$, include positive and negative numbers and fractions. Writing $x\in\mathbb{R}^{d}$ says that $x$ has $d$ real entries. A matrix with $n$ rows and $d$ columns has shape $n\times d$, also written $\mathbb{R}^{n\times d}$. Shapes count entries; they do not indicate physical units. The zero vector has every entry zero. The identity matrix $I$ is square, with ones on its main diagonal and zeros elsewhere; multiplying by it leaves a compatible vector unchanged.

Single bars $|a|$ mean the absolute value of a scalar, its size without its sign. Double bars $\lVert x\rVert$ mean the Euclidean length of a vector: square its entries, add, and take a square root. The dot product multiplies corresponding entries of two equally long vectors and adds the products. Both operations receive explicit calculations in Chapter 2. Adjacent matrices or vectors mean a matrix product only when their inner dimensions agree; Chapter 3 works through that rule.

0.3 The book’s main objects

This is a map of symbols, not a formula to memorize. The chapters repeat the definitions at the point of use.

The symbol $V$ appears in both the eigenvector and singular value chapters because the directions are closely related; the surrounding data table determines its shape. The symbol $\Sigma$ always means covariance here, except that an original SVD diagram uses it for the singular-value matrix and its caption explains that local convention. This prose writes that matrix as $\Sigma_{\mathrm{svd}}$.

Across the series, the mathematical objects stay the same while storage conventions can differ. Here observations are columns in projection formulas and their transposes form data-table rows; in row-based formulas the same multiplication is transposed. Here $\Sigma$ means covariance and $Z$ means total fitting weight. The History companion uses $\Gamma$ for covariance and $Z$ for a standardized news table in some lessons. Read the local definition before transferring a formula. The present book reserves $\Gamma$ for the covariance after a naming rotation and keeps these conventions consistent throughout.

1 Why a newspaper needs linear algebra

Eigen Times starts from a claim about news: that most of it recurs. Wars, scandals, market shocks, elections, disasters happen again and again with different names attached. If that is true, then the enormous variety of news articles should be mostly explained by a small number of recurring directions, and each day’s news should be describable as “so much of this recurring thing, so much of that one, plus a remainder those directions do not capture.”

Linear algebra gives that sentence a precise meaning. We will build the meaning before using the machinery. A scalar is one number. A vector is an ordered list of numbers: in a two-coordinate illustration, $(3,1)$ says three units in the first coordinate and one in the second. Order matters; exchanging the entries changes the vector. A nonzero vector points from the origin, where every coordinate is zero, in a particular direction.

For a first example, suppose our representation records only the horizontal direction $(1,0)$. Taking three times that direction produces $(3,0)$. This is a reconstruction: the part of the observation that our chosen direction can rebuild. Subtract reconstruction from observation, coordinate by coordinate:

$$(3,1)-(3,0)=(0,1).$$

That remainder is the residual. It records what the reconstruction missed. A residual is a vector before it is a length or a score; here its first entry is zero and its second entry is one. The later measurement chapter will explain exactly how to measure its size. For now, the decomposition is visible: observation equals reconstructed part plus residual.

The number three is a coordinate relative to the chosen direction: it tells us how much of that direction to use. A weighted combination adds several directions after multiplying each by a scalar. All combinations obtainable from a collection of directions form their span, a subspace. Our single horizontal direction spans a line; adding a genuinely different direction can span a plane. A low-dimensional representation keeps a limited collection of such directions. The question for Eigen Times is which directions let us reconstruct recurring patterns in news while making the remainder useful to inspect.

The notebook begins with three synthetic observations, $(1,0)$, $(2,0)$, and $(3,1)$. Their horizontal coordinates are 1, 2, and 3. The first two are rebuilt completely; only the third has a nonzero remainder. These are arithmetic examples, not measured levels of war or banking. Actual text representations acquire their coordinates in the next chapter.

The chapters follow the pipeline. Text becomes vectors (§2). Vectors are stacked into matrices, and matrices are multiplied (§3). The covariance of a million article vectors is computed and its eigenvectors found (§4). The singular value decomposition works directly on the data matrix and gives us latent semantic analysis (§5). The eigenvectors are oriented and rotated into axes with names (§6). A new story is projected, and what is left over is measured (§7). Then come the thresholds (§8–§9), matching axes across bases and nights (§10), streaming the calculation (§11), and the shapes of everything (§12).

2 Text as vectors

2.1 Counting words

The oldest way to turn a document into numbers is to count its words. Let $p$ count the terms in a fixed vocabulary; a document becomes a row of $p$ counts. A corpus is the collection of documents being analysed. Six tiny “articles”, labelled d1 through d6, and seven retained terms are enough to see the construction:

d1  troops shell border town
d2  ceasefire troops border
d3  bank raises rate
d4  rate rise hits mortgage
d5  troops ceasefire talks border
d6  bank mortgage rate rise

Terms that occur in only one document (shell, town, raises, hits, talks) are dropped in this example—they provide no shared term for connecting documents—leaving troops, border, ceasefire, bank, rate, mortgage, rise. Fix that order before writing any numbers. Document d1 then has row $(1,1,0,0,0,0,0)$: one occurrence of troops, one of border, and none of the other retained terms. Document d3 has row $(0,0,0,1,1,0,0)$. Their first entries refer to the same term even though the documents concern different subjects.

Stacking the six rows creates the count matrix, the left grid of Figure 1. Reading across answers “which retained words does this document use?” Reading down the troops column finds three documents containing that word. A count inside one row and a count of documents down a column answer different questions. The weighting rule below uses both.

2.2 TF‑IDF

TF‑IDF means term frequency–inverse document frequency. Its two factors answer two questions: how often does this document use a term, and how widely does that term occur across documents? Raw counts favour long documents and give common and rare words the same weight per occurrence. We soften the first effect and adjust the second.

2.2.1 Logarithms before the weighting rule

A power repeats multiplication: $2^3=2\times2\times2=8$. A logarithm reverses this question: to what power must a chosen base be raised to produce a given positive number? The natural logarithm, written $\ln$, uses the fixed positive base $e\approx2.71828$. Thus $\ln(e)=1$ and $\ln(1)=0$, because the first and zeroth powers of $e$ are $e$ and 1. The operator $\approx$ means approximately equal. Real-number powers extend this relationship between whole-number powers; we do not need to construct that extension to use a calculator’s log function.

The useful property is that doubling a positive input adds the same amount, $\ln(2)$, to its logarithm. Doubling a count repeatedly therefore adds rather than doubles its contribution. For positive counts 1, 2, 4, and 8, the values of $1+\ln(\text{count})$ are approximately 1.000, 1.693, 2.386, and 3.079. Eight mentions still count more than one, but they do not receive eight times the weight. This is damping.

Let $t$ identify a vocabulary term, $N$ count documents, and $\mathrm{df}(t)$ count the documents containing $t$. For a positive within-document count, define damped term frequency by $\mathrm{tf}=1+\ln(\text{count})$. An absent term gets zero directly; we never ask for $\ln(0)$, which is not defined as a real number. Define inverse document frequency, denoted $\mathrm{idf}(t)$, by

$$\mathrm{idf}(t) = \ln\frac{N + 1}{\mathrm{df}(t) + 1} + 1,$$

Read the formula from inside outward. Add one to both document counts, divide, take the natural logarithm, then add one. The additions inside the fraction are a smoothing convention; the addition outside ensures that a term present in every document still has idf 1. Fewer containing documents make the denominator smaller, the ratio larger, and the idf larger.

There are six documents. Troops appears in three, so its calculation is $1+\ln(7/4)\approx1.56$. Ceasefire appears in two, so its calculation is $1+\ln(7/3)\approx1.85$. Multiply each term’s damped frequency by its idf. In these tiny documents each retained term occurs at most once, so each nonzero damped frequency is 1. Weighting changes how much each present word contributes without inventing a contribution for an absent word.

2.2.2 Turning a weighted row into a unit row

We next remove the row’s overall scale. For a row $a$ with entries $a_t$, the squared Euclidean length is the sum of the squared entries, $\sum_t a_t^2$. The summation sign means add over all vocabulary positions $t$. Squaring means multiplying an entry by itself. The Euclidean length, also called the norm and written $\lVert a\rVert$, is the nonnegative square root of that sum. A square root undoes squaring for a nonnegative result. In symbols,

$$\lVert a\rVert=\sqrt{\sum_t a_t^2}.$$

The double bars denote a vector’s length. Single bars around one scalar will mean its absolute value, or magnitude without sign: $|-2|=2$. A length is nonnegative even when some vector entries are negative.

For d2, only the first three weighted entries are nonzero. Using the unrounded idf values, its length is approximately 2.8770. Dividing each entry by that same length gives approximately $(0.5421,0.5421,0.6421,0,0,0,0)$. Squaring and adding the unrounded normalized entries gives 1. Such a vector is a unit vector. Dividing every entry by one common positive number changes length without changing direction. This is normalization; it leaves the relative proportions of words in the row intact. An all-zero row has length zero and cannot be normalized by division.

Eigen Times builds this representation over 1,093,166 articles and a 50,000‑term vocabulary, keeping terms that occur in at least 20 documents and at most 30% of them. The result has 149 million non‑zero entries out of 55 billion cells—it is 99.7% zeros. Sparse storage keeps only nonzero entries and their positions; dense storage allocates every cell. The algorithms of §5 avoid forming this full matrix densely.

2.3 Cosine similarity

Let $a$ and $b$ be two nonzero rows with the same vocabulary order, and $a_t,b_t$ their entries for term $t$. Their dot product, also called the Euclidean inner product and written $a\cdot b$, multiplies corresponding entries and adds them. For three entries this means $a_1b_1+a_2b_2+a_3b_3$; for a vocabulary it means $\sum_t a_tb_t$. Shared nonzero coordinates can contribute; a product involving an absent coordinate is zero.

The dot product alone grows if we multiply either vector by a large positive number. Cosine similarity, written $\cos(a,b)$, removes this dependence on overall length by dividing by both norms:

$$\cos(a, b) = \frac{a \cdot b}{\lVert a\rVert \, \lVert b\rVert } = \sum_t a_t b_t \quad \text{when } \lVert a\rVert = \lVert b\rVert = 1.$$

A cosine of 1 means the same direction, 0 means perpendicular directions, and $-1$ means opposite directions. In geometric language cosine is the horizontal coordinate of a unit direction measured relative to the first direction; this connects the formula to angles. Its value always lies between −1 and 1. For nonnegative word-count or TF‑IDF rows, negative cosines cannot occur; later centered or embedded vectors can have signed entries. The zero row has no defined cosine because the denominator would be zero.

For d2 compared with itself, the dot product is its first normalized entry squared, plus the second squared, plus the third squared; the remaining products are zero. Normalization made that sum 1. Document d5 has exactly the same retained row, so the same calculation gives 1 for d2 versus d5. For d1 versus d3, every coordinate product is zero: wherever one row has a nonzero retained term, the other has zero. Their dot product, and therefore their cosine, is zero.

In Figure 1, d2 and d5 share troops, border, ceasefire in the same proportions and have cosine 1.0; d1 and d3 share nothing and have cosine 0. Cosine is the similarity used throughout Eigen Times—for clustering a day’s articles into stories, for threading stories into episodes, for finding precedents—so §8 is largely about what values of it mean.

2.4 Embeddings

A TF‑IDF vector knows only which words appear. Two headlines using no shared retained terms have perpendicular, or orthogonal, TF‑IDF rows even if their meanings agree. A sentence embedding is a vector produced by a trained neural network—a numerical model that learns from examples—in which cosine similarity can track meaning beyond vocabulary. Eigen Times uses bge-small-en-v1.5, trained on hundreds of millions of sentence pairs, to turn a headline and lede (opening text: the first 300 characters, at most 96 tokens, or model text units) into 384 numbers. We treat the network as a black box that gives us a vector; everything downstream is linear algebra.

An embedding coordinate need not name a recognizable topic. Meaning is represented by the pattern across coordinates and relationships among vectors. This is why a direction in embedding space will later need evidence from words before we give it a readable name.

These embedding vectors are dense and not centred: their mean—the coordinate-by-coordinate average—has not been subtracted. To see why this matters, consider the synthetic vectors $(1,1,0)$ and $(1,0,1)$. Their dot product is 1 and each has length $\sqrt2$, so their cosine is $1/(\sqrt2\sqrt2)=0.5$. Remove the shared component $(1,0,0)$ from both. The remaining vectors are $(0,1,0)$ and $(0,0,1)$, whose cosine is zero. Normalizing the original vectors would not remove that shared component; subtracting it changes their directions.

News embeddings also share a large common component, and the observed cosine between unrelated articles is around 0.5, not 0. The synthetic calculation demonstrates a mechanism; it does not reproduce that empirical distribution. The observation drives two of the thresholds in §8 and §9.

3 Matrices, pictured

3.1 A matrix times a vector

Let $n$ count document rows and $p$ count their coordinates. Stacking them gives an $n\times p$ matrix $X$: the shape always lists rows before columns. Individual observation vectors will be written as columns when used in projection or outer-product formulas; their transposes supply the rows of a data matrix. A transpose exchanges rows and columns and is marked $\top$. Thus a column $x$ has row form $x^{\top}$.

For Figure 2, let $A$ be a $4\times3$ example matrix and $v$ a $3\times1$ column of weights, with entries $v_1$ through $v_3$. Their product $Av$ is a $4\times1$ column. Adjacent mathematical symbols mean multiplication when their shapes permit it.

The actual inputs and result are

$$A=\begin{pmatrix}2&0&1\\1&3&0\\0&1&2\\1&1&1\end{pmatrix},\qquad v=\begin{pmatrix}1\\2\\-1\end{pmatrix},\qquad Av=\begin{pmatrix}1\\7\\0\\2\end{pmatrix}.$$

The second output entry is $1\times1+3\times2+0\times(-1)=7$. Each row of $A$ has three entries and meets all three entries of $v$. This is why the number of columns in the matrix must equal the number of entries in the input column. There are four rows to evaluate, so the output has four entries. A negative weight is allowed: it subtracts a contribution rather than adding one.

The definition of $Av$ is “dot each row of $A$ with $v$”, but the picture to keep is the other one (Figure 2): $Av$ is a weighted sum of the columns of $A$, with the entries of $v$ as the weights. If the columns of $A$ are directions, $Av$ is the point you reach by going $v_1$ of the way along the first, $v_2$ along the second, and so on. That is what “coordinates” will mean in §7: the coordinates of a story are the weights that best rebuild it from the basis directions.

In this example, take the first column once, the second column twice, and subtract the third column once. The second entries combine as $1+2\times3-0=7$, just as in the row calculation. Both descriptions perform the same multiplications; they group them differently. We will use row products to measure coordinates and column combinations to rebuild observations.

3.2 A matrix times a matrix

For Figure 3, take a new $3\times4$ example $A$, let $B$ be $4\times2$, and call their product $C=AB$. Indices $i$ and $j$ select a row and column: entry $C_{ij}$ is the dot product of row $i$ of $A$ with column $j$ of $B$. Shapes must agree in the middle: $(3\times4)(4\times2)$ gives $3\times2$. Each column of $B$ is a separate input vector; multiplying by a matrix carries out both matrix–vector calculations together.

Here are the inputs, with the second row and second column providing a short calculation:

$$A=\begin{pmatrix}1&2&0&1\\0&1&1&2\\2&0&1&0\end{pmatrix},\qquad B=\begin{pmatrix}1&0\\2&1\\0&3\\1&1\end{pmatrix}.$$

The entry $C_{22}$ is $0\times0+1\times1+1\times3+2\times1=6$. The four intermediate products correspond to the common inner dimension, 4. The two outer dimensions say that there will be three output rows and two output columns. Matrix multiplication is not entry-by-entry multiplication, and reversing the order need not give the same shape or result.

Two special cases recur. If $X$ has $n$ rows and $p$ columns, $X^{\top}X$ has shape $p\times p$. A row of $X^{\top}$ is a column of $X$, so each output entry compares two columns of $X$ by their dot product. This is called a Gram matrix, and it is the raw material of covariance. In contrast, $XX^{\top}$ compares pairs of observation rows and has shape $n\times n$.

For the other special case, let $k$ count selected, mutually perpendicular unit directions and let $V$ contain them as $p\times1$ columns. Then $XV$ has shape $n\times k$. Entry $(i,j)$ asks how much observation $i$ points along direction $j$. This table contains every document’s coordinates on those directions: the archive projected onto them. We will derive the reconstruction corresponding to those coordinates rather than assuming that taking a projection preserves everything.

3.3 Transpose, symmetry, orthogonality

$A^{\top}$ flips rows and columns. A square matrix has the same number of rows and columns. Call such a matrix $S$; it is symmetric when $S=S^{\top}$. Real symmetric matrices, including covariances, admit a complete set of real, mutually perpendicular eigenvector directions, defined in the next chapter.

3.3.1 Perpendicular directions and what their products say

Two vectors are perpendicular when their dot product is zero. Unit-length perpendicular columns are orthonormal: “ortho” refers to perpendicularity and “normal” to unit length. The identity matrix $I$ has ones on its diagonal and zeros elsewhere; multiplying by it changes nothing. A $p\times k$ matrix $V$ with orthonormal columns satisfies $V^{\top}V=I$, with a $k\times k$ identity. Each diagonal entry in this product is a column’s squared length, so it is 1. Each off-diagonal entry compares two different columns, so it is 0.

Directions are independent when none can be formed as a weighted combination of the others. A basis is an independent collection spanning the space under discussion: its combinations reach every point in that space. The dimension counts the number of basis directions. A matrix’s rank counts the independent directions among its columns, equivalently among its rows. Listing a direction twice does not give us an extra independent direction.

3.3.2 A rotation keeps a whole space; a projection keeps part

If $k=p$, the matrix $V$ is square and orthogonal. It has enough orthonormal columns to describe the whole space. Multiplication by $V$ or $V^{\top}$ rotates or reflects coordinates and preserves lengths and angles. It also satisfies $VV^{\top}=I$, so projecting and rebuilding recovers every input exactly.

If $k<p$, the map $V^{\top}x$ keeps only selected coordinates of an observation column $x$. Rebuilding gives $VV^{\top}x$, its orthogonal projection onto the retained subspace. Now $VV^{\top}$ need not be the identity: the missing directions cannot be recovered merely by changing back to the original number of coordinates.

For the notebook’s example, take $x=(1,2,3)^{\top}$ and let $V$ contain the first two coordinate directions in three dimensions. The superscript $\top$ makes the displayed list a column. Projection gives $(1,2)^{\top}$ and rebuilding gives $(1,2,0)^{\top}$. The squared lengths are $1+4+9=14$ before projection and $1+4=5$ afterward. Nine units of squared length lie in the discarded coordinate. The rectangular matrix obeys $V^{\top}V=I$, but the reconstruction operator is the diagonal matrix with entries 1, 1, 0, not a three-dimensional identity.

A quarter-turn in the retained plane maps $(1,2)^{\top}$ to $(-2,1)^{\top}$. The squared length remains $4+1=5$. This separates two operations that the later pipeline uses for different purposes: retaining a subspace loses information; rotating coordinates within it does not.

4 The covariance and its eigenvectors

4.1 Centering and covariance

4.1.1 From an average to a measure of spread

The mean is an average: add the observations and divide by their count. For the scalar values 1, 2, and 3, the mean is $(1+2+3)/3=2$. Their deviations are the observed values minus that mean: −1, 0, and 1. Subtracting a common reference mean is called centering. The centered values say where each observation lies relative to that reference rather than relative to zero on the original scale.

Adding signed deviations gives zero here. That does not mean there is no spread; the negative and positive deviations cancel. Squaring avoids this cancellation. The squared deviations are 1, 0, and 1, whose mean is $2/3$. This average squared deviation is the variance under the total-count convention used for the fitted covariance. Its nonnegative square root, $\sqrt{2/3}\approx0.8165$, is the standard deviation. If observations are measured in some unit, variance is measured in that unit squared; standard deviation returns to the original unit.

Some statistical estimates divide the sum of squared deviations by one less than the count, to correct a bias when estimating a population variance from a sample. That is a different convention. The fitted covariance below divides by total weight. We will identify the reference and divisor when using a standard deviation elsewhere.

4.1.2 How two coordinates vary together

For more than one coordinate, calculate a mean for each coordinate. Consider the two observation rows $(1,2)$ and $(3,4)$. Their mean row is $(2,3)$, and their centered rows are $(-1,-1)$ and $(1,1)$. The first-coordinate squared deviations average to 1; the second-coordinate squared deviations also average to 1. The products between the first and second deviations are $(-1)(-1)=1$ and $(1)(1)=1$, so their average is also 1.

That last average is the covariance between the two coordinates. It is positive when deviations tend to have matching signs, negative when they tend to have opposite signs, and zero when their average product is zero. Statistical independence means that knowing one variable supplies no information about the other. Zero covariance does not generally establish independence: two variables may still have a relationship that this signed product does not detect.

Arrange these averages in a table. Diagonal entries compare a coordinate with itself and are variances; off-diagonal entries are covariances. For the two-row example the table is

$$\begin{pmatrix}1&1\\1&1\end{pmatrix}.$$

To compute every pairwise product at once, place a centered observation in a column and multiply it by its transpose. This outer product produces a square table. Its entry in row $j$, column $l$ is the product of centered coordinates $j$ and $l$. The two indices select coordinates, not articles. Averaging these tables produces the covariance matrix.

4.1.3 Fit weights change the average, not the operation

Let $i$ index fitted articles, from 1 to $n$, and let $d$ count embedding coordinates. Each $x_i$ is a $d\times1$ column; row $i$ of the $n\times d$ data matrix $X$ is $x_i^{\top}$. Let $w_i$ be its nonnegative fit weight. A weighted mean gives each observation the stated weight, adds the weighted contributions, then divides by the total weight. Let $Z>0$ be that total and $\mu$ the $d\times1$ weighted mean. A sum over $i$ includes every fitted article, and $Z^{-1}$ means $1/Z$:

$$Z = \sum_i w_i, \qquad \mu = Z^{-1} \sum_i w_i x_i,$$

For example, give the two observation rows above weights 1 and 3. The total is 4, the first mean coordinate becomes $(1\times1+3\times3)/4=2.5$, and the second becomes 3.5. Centering must use this new weighted mean. The centered rows are $(-1.5,-1.5)$ and $(0.5,0.5)$. The weighted mean of either squared coordinate is $(1\times2.25+3\times0.25)/4=0.75$. Every entry of this example’s weighted covariance is 0.75.

Let $\Sigma$ name the $d\times d$ weighted covariance. Apply precisely the same weights to the outer products of deviations from the weighted mean:

$$\Sigma = Z^{-1} \sum_i w_i (x_i - \mu)(x_i - \mu)^{\top}.$$

The first edition uses uniform article weights; the later archive fit uses inverse daily article counts, explained in §11. These weights are distinct from the story energy used to rank news. In Eigen Times $d=384$, so $\Sigma$ is a $384\times384$ matrix—147,456 numbers summarising the fitted article cloud. Figure 4 pictures the operations in their working order: center each observation, take its outer product, and average the products.

4.1.4 Why three sums suffice to update the covariance

A practical point matters for updates. Let $m$ be the $d\times1$ weighted coordinate sum and $M$ the $d\times d$ weighted outer-product sum:

$$m=\sum_i w_i x_i,\qquad M=\sum_i w_i x_i x_i^{\top}.$$

The mean is immediately $m/Z$. To recover covariance, first expand a single centered outer product:

$$(x_i-\mu)(x_i-\mu)^{\top} =x_ix_i^{\top}-x_i\mu^{\top}-\mu x_i^{\top}+\mu\mu^{\top}.$$

Now take the weighted average of all four terms. The first becomes $M/Z$. In the second and third terms, the weighted average of $x_i$ is $\mu$, so each becomes $-\mu\mu^{\top}$. The last term is the same for every observation and averages to $+\mu\mu^{\top}$. Two subtractions and one addition leave one subtraction:

$$\mu=m/Z,\qquad \Sigma=M/Z-\mu\mu^{\top}.$$

The three running sums $Z,m,M$ therefore retain everything needed for these mean and covariance calculations. They do not retain the full article collection. New observations with fixed weights can be added without revisiting old observations. If a day’s final article count changes its existing weights, the old contributions must also be corrected. Although the identity is exact, subtracting large nearly equal floating-point numbers can lose accuracy; the streaming chapter distinguishes algebraic identities from numerical reliability.

For compact matrix notation, let $W=\mathrm{diag}(w_1,\ldots,w_n)$ be the $n\times n$ diagonal weight matrix: $\mathrm{diag}$ puts the listed values on the diagonal and zeros elsewhere. Let $\mathbf1$ be a column of $n$ ones. Multiplying $\mathbf1\mu^{\top}$ repeats the mean row $n$ times. The centered data matrix is therefore $X_c=X-\mathbf1\mu^{\top}$. Multiplying $WX_c$ scales each centered row by its weight; multiplying on the left by $X_c^{\top}$ then sums the weighted products. This gives the same covariance as $\Sigma=X_c^{\top}WX_c/Z$.

4.2 Eigenvectors: the axes of the cloud

A cloud of points in $d$ dimensions has a shape. It can be broad in one direction and narrow in another, even when neither direction follows an original coordinate axis. We want to measure that shape using perpendicular directions. The special directions of the covariance matrix will supply them.

A nonzero $d\times1$ column $v$ is an eigenvector of $\Sigma$ with eigenvalue $\lambda$ when $\Sigma v=\lambda v$. The left side multiplies the vector by a matrix; the right side multiplies it by a scalar. Equality says that this input stays on its own line under the matrix operation. For a positive eigenvalue it points the same way and changes length by that factor; a zero eigenvalue sends it to zero.

4.2.1 Check an eigenpair by multiplication

Consider the separate exact example

$$\Sigma=\begin{pmatrix}2&1\\1&2\end{pmatrix}.$$

Multiplying by $(1,1)^{\top}$ gives $(3,3)^{\top}=3(1,1)^{\top}$. Multiplying by $(1,-1)^{\top}$ gives $(1,-1)^{\top}=1(1,-1)^{\top}$. We have found two eigendirections and eigenvalues 3 and 1. Their dot product is zero. Their lengths are both $\sqrt2$, so divide each by $\sqrt2$ to obtain unit columns, called $v_1$ and $v_2$. Scaling an eigenvector to unit length does not change its eigenvalue.

This is a way to verify an answer, not an assumption that every covariance’s eigenvectors are easy to guess. Numerical algorithms find them for larger matrices. The purpose of the example is to make the equation an operation we can inspect.

For a covariance, eigenvalues are nonnegative. Here is why. Let $u$ be any unit column direction and let $a_i=u^{\top}(x_i-\mu)$ be observation $i$’s scalar projection score along it. Its weighted mean is zero because the observations were centered. Its weighted variance is the average squared score. Substituting the covariance definition gives

$$Z^{-1}\sum_i w_i a_i^2=u^{\top}\Sigma u.$$

The left side is an average of nonnegative squares. If $u$ is a unit eigenvector with eigenvalue $\lambda$, the right side becomes $u^{\top}(\lambda u)=\lambda(u^{\top}u)=\lambda$. Thus the eigenvalue is exactly the variance along that unit direction and cannot be negative.

4.2.2 Write the whole covariance in its own directions

Choose $d$ unit, mutually perpendicular eigenvectors $v_j$, ordered by decreasing eigenvalue $\lambda_j$, with $j=1,\ldots,d$. Real symmetric matrices permit such a choice, although repeated eigenvalues need not specify unique individual directions. Let $V$ contain those columns and $\Lambda$ the diagonal matrix of their eigenvalues, both $d\times d$. They diagonalise the covariance, meaning

$$\Sigma = V\Lambda V^{\top},\qquad V^{\top}V=I,\qquad \Lambda=\mathrm{diag}(\lambda_1,\ldots,\lambda_d),$$

where $\lambda_1\ge\lambda_2\ge\cdots\ge\lambda_d\ge0$ and $\ge$ means greater than or equal to. The diagonal representation has zero covariances between different eigenvector coordinates. The principal components are these directions, ranked by variance; using them to summarize the cloud is principal component analysis (PCA).

Read the factorization from right to left when applying it to a column: $V^{\top}$ expresses the input in eigenvector coordinates, $\Lambda$ scales each coordinate separately, and $V$ expresses the result back in the original coordinates. The covariance is simple in this basis because no off-diagonal terms mix the coordinates.

Why does the first direction capture the most variance? Express an arbitrary unit direction as a combination $u=\sum_j\alpha_jv_j$, where each scalar $\alpha_j$ is its coordinate in the full orthonormal eigenbasis. Its unit length means $\sum_j\alpha_j^2=1$. Its projected variance becomes

$$u^{\top}\Sigma u=\sum_j\lambda_j\alpha_j^2.$$

The squared coordinates are nonnegative and add to one. The variance is therefore a weighted average of the eigenvalues and cannot exceed the largest, $\lambda_1$. Choosing $u=v_1$ reaches that largest value. Requiring the next direction to be perpendicular to $v_1$ sets its first coordinate to zero; the same argument then selects $v_2$. This explains the ordering without assuming calculus.

Figure 5 shows 300 points in two dimensions: the long axis has eigenvalue 4.06, the short one 0.44, so about 90% of the displayed variation lies along one direction. The synthetic distribution has known mean zero. This illustration forms second moments about that origin, rather than subtracting the finite sample’s empirical mean; the displayed eigenvalues estimate population variances. The notebook computes both choices so that centering is an explicit decision.

This is the content of “the eigen‑news are the directions that explain the most variance”. Variance describes spread in the chosen numerical representation; the argument alone does not establish significance, truth, or editorial importance. Computing eigenvectors for a $384\times384$ symmetric matrix is a standard numerical operation (eigh in common linear-algebra libraries; faer in the Rust implementation) and takes milliseconds in the reported fit.

4.3 How many axes to keep

Let $k\le d$ count retained directions; $\le$ means less than or equal to. The explained-variance fraction is $\sum_{j=1}^{k}\lambda_j/\sum_{j=1}^{d}\lambda_j$, provided total variance is positive. It reports how much of the fitted cloud’s variation those directions retain. For the Eigen Times embedding basis, 60 components carry 55% of the variance, and the curve is flat beyond; a small eigenvalue alone does not prove that a direction is meaningless.

In the exact two-direction example, total variance is $3+1=4$. Keeping the first direction retains $3/4=75\%$ and discards $1/4=25\%$. Keeping both retains all the variance. With a longer list of eigenvalues, make the same running calculation: add the retained values, divide by the sum of all values, and inspect what each extra direction contributes. A flat part of that curve means that each new direction adds little to the numerator. It does not mean that the discarded articles or subjects are unimportant.

The Marchenko–Pastur edge asks a different question: how large might covariance eigenvalues be in a specified model of featureless variation? Here noise means variation generated by that reference model, not a judgment about an article. The spectrum of a covariance is its collection of eigenvalues. In this illustration, let $N$ count independent observations in $d$ coordinates with equal noise variance $\sigma^2$, unrelated to the singular-value notation used later. For a concrete reference, suppose all entries are independent draws from the same zero-mean Gaussian distribution, a symmetric bell-shaped probability model, with variance $\sigma^2$. Independence applies across both observations and coordinates; repeated draws use the same distribution. These assumptions specify an illustrative noise model, not a description of news embeddings. Let $\lambda_+$ denote the approximate upper edge of the model’s covariance spectrum. Under that model, in the large-sample limit,

$$\lambda_{+} = \sigma^2 \left(1 + (d/N)^{1/2}\right)^2,$$

Here raising to $1/2$ means taking a square root. First form the dimension-to-observation ratio $d/N$, take its square root, add one, square, and multiply by the stipulated variance. The notebook uses $\sigma^2=0.01$, $N=1{,}000{,}000$, and $d=384$, giving approximately 0.010396. That is slightly larger than 0.01: finite sampling can spread estimated variances even when the model gives every coordinate the same underlying variance. This teaching input is not an estimate fitted to news.

Eigenvalues near or below such a reference may be noise; this is not a universal significance test for correlated, weighted news data. The median eigenvalue supplies a rough noise-variance estimate in the reported comparison. With $N\approx10^6$ and $d=384$, the ratio $d/N$ is tiny and the edge is close to that estimate. The deployed bases use fixed $k=60$, with the variance curve as a check. From here onward $V$ usually means only the retained $d\times k$ columns and $\Lambda$ their $k\times k$ positive eigenvalue matrix; full decompositions will be identified explicitly.

5 The singular value decomposition and LSA

5.1 What the SVD is

The eigendecomposition works on a square covariance. The singular value decomposition (SVD) works directly on a rectangular data matrix. It separates a matrix operation into input directions, nonnegative scale factors, and output directions. The input and output can have different dimensions, so the directions on the two sides need not live in the same space.

5.1.1 Follow one input direction through a rectangular table

For this chapter let $X$ be any $n\times d$ matrix. A column input has $d$ entries and the output $Xv$ has $n$ entries. In the small example

$$X=\begin{pmatrix}3&0\\0&2\\0&0\end{pmatrix},$$

the input $(1,0)^{\top}$ maps to $(3,0,0)^{\top}$, which is three times the unit output direction $(1,0,0)^{\top}$. The input $(0,1)^{\top}$ maps to $(0,2,0)^{\top}$, which is twice $(0,1,0)^{\top}$. These perpendicular input directions map to perpendicular output directions, with scale factors 3 and 2. The output has three entries even though the input has two.

For a general matrix, SVD finds suitable directions even when the input table is not already diagonal. Let $v_j$ be a unit right singular direction in the input space, let $u_j$ be its unit left singular direction in the output space, and let $\sigma_j$ be their nonnegative scale. The key relation is $Xv_j=\sigma_ju_j$. “Right” and “left” describe where the directions appear in the factorization; neither means an orientation on a page. A zero scale contributes no output in that direction.

5.1.2 Assemble the separate directions into matrices

Let $U$ be an $n\times n$ orthogonal matrix of left directions, and $V$ a $d\times d$ orthogonal matrix of right directions. Let $\Sigma_{\mathrm{svd}}$ be an $n\times d$ rectangular diagonal matrix: its only potentially nonzero entries lie where row and column indices agree. Those entries are the singular values $\sigma_j$ in decreasing order. Then

$$X=U\Sigma_{\mathrm{svd}}V^{\top}.$$

Figure 6 labels this singular-value matrix with the traditional bare $\Sigma$; it is a different object from the covariance $\Sigma$ in §4. The columns of $V$ are directions indexed by features, and those of $U$ are directions indexed by observations. Actual coordinate scores are $XV=U\Sigma_{\mathrm{svd}}$, not $U$ alone: each left direction must receive its singular-value scale.

Read the product from right to left, as with eigendecomposition. $V^{\top}$ resolves an input into right-direction coordinates; $\Sigma_{\mathrm{svd}}$ scales those coordinates and accounts for the change in dimension; $U$ writes the result in the original output coordinates. For the $3\times2$ example, full factors can use a $3\times3$ identity for $U$, a $2\times2$ identity for $V$, and the displayed rectangular table for $\Sigma_{\mathrm{svd}}$. Nothing needs turning because the table already uses suitable directions.

Let $q=\min(n,d)$, where $\min$ selects the smaller number. The thin SVD removes the unused full-matrix directions: it has $U$ of shape $n\times q$, $\Sigma_{\mathrm{svd}}$ of shape $q\times q$, and $V$ of shape $d\times q$, with the same product $X$. Keeping only $k<q$ directions is the further operation of truncation.

In the example, $q=2$. The thin $U$ keeps the first two columns of the three-dimensional identity, the scale matrix becomes $\mathrm{diag}(3,2)$, and $V$ remains the two-dimensional identity. Multiplying these factors still reconstructs every entry. Thin factorization changes storage without losing any information; truncation will deliberately lose contributions.

5.1.3 Why covariance eigenvalues contain squared singular values

Transposing a product reverses its order: $(AB)^{\top}=B^{\top}A^{\top}$. Apply that rule to the SVD, multiply, and use $U^{\top}U=I$. For the full rectangular decomposition the steps are

$$X^{\top}X =V\Sigma_{\mathrm{svd}}^{\top}U^{\top}U\Sigma_{\mathrm{svd}}V^{\top} =V(\Sigma_{\mathrm{svd}}^{\top}\Sigma_{\mathrm{svd}})V^{\top}.$$

Each diagonal entry in the middle is a singular value multiplied by itself. For the thin decomposition the middle product is the square diagonal matrix $\Sigma_{\mathrm{svd}}^2$. Thus the right singular directions are also eigenvectors of the Gram matrix, with squared singular values as eigenvalues. In the exact example, the Gram matrix is $\mathrm{diag}(9,4)$.

If $X$ is centered and equally weighted, covariance is that Gram matrix divided by $n$. Its right directions therefore stay the same, and its eigenvalues are $\lambda_j=\sigma_j^2/n$. The diagonal example above is not centered, so its Gram matrix should not silently be called an empirical centered covariance. The notebook also applies the identity to the centered six-document table.

For the weighted covariance of §4, instead decompose $W^{1/2}X_c$. Here $W^{1/2}$ is diagonal with entries $\sqrt{w_i}$, so it scales each centered row by the square root of its weight. When we form a Gram product, one square-root weight comes from each factor; their product is the full weight $w_i$. Divide the resulting squared singular values by total weight $Z$. Scores for the original, unweighted observations are still $X_cV$, not the weighted rows’ scores. A direction’s canon is its representative collection of high-scoring stories; a left singular vector alone is not generally the original article-score column in a weighted fit.

5.2 Truncation

Let $U_k$ and $V_k$ contain the first $k$ columns of the thin factors, and let $\Sigma_k$ contain their $k$ singular values on its diagonal. Their shapes are $n\times k$, $d\times k$, and $k\times k$, respectively. The approximation $X_k=U_k\Sigma_kV_k^{\top}$ has rank at most $k$. It keeps only $k$ independent contributions with which to rebuild the table.

5.2.1 Measure what a reconstruction leaves behind

The residual table is original minus reconstructed, $X-X_k$. A positive residual entry says that the reconstruction was smaller than the original there; a negative entry says it was larger. Adding those signed errors would allow overestimates and underestimates to cancel. Instead square every entry and add. This is the sum of squared errors. Least squares means choosing a reconstruction that makes such a sum as small as the allowed family permits.

The Frobenius norm, written $\lVert A\rVert_F$ for a matrix $A$, is the square root of the sum of all its squared entries. It extends vector length to a whole table by treating every cell as a coordinate. Minimizing this norm or its square gives the same minimizing table, because the square-root function preserves order among nonnegative numbers.

In the rectangular diagonal example, keeping the scale 3 and dropping the scale 2 reconstructs only the entry 3. The residual has a single nonzero entry, 2. Its squared Frobenius norm is therefore 4, and its norm is 2. The original table’s squared norm is $9+4=13$. Define relative error as residual norm divided by the positive original norm; here it is $2/\sqrt{13}\approx0.5547$. An error ratio is not a proportion of articles reconstructed correctly; it compares geometric lengths in this representation.

5.2.2 A small detour: what minimizing and differentiating mean

Before asking a matrix theorem to minimize anything, solve a scalar fitting problem. Suppose we must replace the observations 1, 2, and 3 by one common fitted value, called $\theta$ (Greek theta). Their residuals are $1-\theta$, $2-\theta$, and $3-\theta$. Let $J(\theta)$ name their sum of squared residuals; $J$ is a function, a rule returning a number for each proposed input $\theta$. Expanding the three squares gives

$$\begin{aligned} J(\theta)&=(1-\theta)^2+(2-\theta)^2+(3-\theta)^2\\ &=3\theta^2-12\theta+14\\ &=3(\theta-2)^2+2. \end{aligned}$$

The last form makes the minimum visible. A square cannot be negative, and $(\theta-2)^2$ is zero exactly when $\theta=2$. Thus the smallest squared error is 2. Trying $\theta=1$ or $\theta=3$ gives squared error 5. At the best fit the signed residuals are −1, 0, and 1; their sum is zero, while their sum of squares is two. These are different quantities.

Differentiation provides another way to inspect how changing the fit changes the error. Let $h$ be a small nonzero change in $\theta$. The difference $J(\theta+h)-J(\theta)$ is the change in error. Divide by $h$ to obtain change in error per unit change in the fitted value, an average slope over that step. Expanding the quadratic and canceling equal terms gives

$$\frac{J(\theta+h)-J(\theta)}{h}=6\theta-12+3h.$$

As the step becomes arbitrarily small, the last term approaches zero. The resulting local slope is the derivative, written $J'(\theta)$ or $\mathrm{d}J/\mathrm{d}\theta$; the prime here means derivative, not the rotated coordinate notation used later:

$$J'(\theta)=6\theta-12.$$

We differentiated to learn which way the error changes. Below 2 the slope is negative, so increasing the fit locally lowers error; above 2 the slope is positive, so increasing it locally raises error. At 2 the slope is zero. A zero derivative is a stationary point, not automatically a minimum for every function. Here the completed-square formula proves it is the global minimum. This example supplies the meaning of derivative, the calculation from a small change, and the reason for setting it to zero. No matrix differentiation is required to use the SVD that follows.

5.2.3 The best reconstruction within a rank budget

The Eckart–Young theorem says that the truncated SVD minimizes $\lVert X-X_k\rVert_F$ among all matrices of rank at most $k$. This comparison matters: we are not saying that an arbitrary reconstruction is best, or that the reconstruction is a true explanation of the subjects. We are choosing the best approximation under this particular squared-error measure and rank constraint.

The squared error equals the sum of squared discarded singular values. To see the accounting, each singular component has a unit left direction, a unit right direction, and a scale $\sigma_j$. Its contribution to squared table length is $\sigma_j^2$. Distinct components are perpendicular when their tables are compared entry by entry, so their squared lengths add without cross terms. Keeping the largest scales retains the most of this length budget. The theorem adds the stronger claim that no different rank-$k$ table does better.

Figure 7 applies this to the six-document matrix of Figure 1. Its singular values are 1.64, 1.51, 0.82, 0.55, and then almost nothing. Relative error means $\lVert X-X_k\rVert_F/\lVert X\rVert_F$ for nonzero $X$. Compute it by squaring and adding the discarded singular values, dividing by the sum of all squared singular values, and taking a square root. Rank 1 rebuilds only the war documents (relative error 0.74); rank 2 rebuilds both blocks, war and banking, and the error drops to 0.42; rank 3 reaches 0.24. Two components already separate the two kinds of document in this small example. The notebook uses unrounded singular values, so its checks do not accumulate rounding error from the printed labels.

5.3 Latent semantic analysis

Latent semantic analysis (LSA) uses a truncated SVD of a document–term matrix such as TF‑IDF; Eigen Times additionally centres and, in the archive edition, weights that matrix. Each retained column of $V$ is a direction over terms. Its entries are loadings: coefficients telling how strongly each input feature contributes to the direction. Its largest term loadings provide a word list (Figure 8). In the toy corpus the first component loads on troops, border, ceasefire, the second on bank, rate, mortgage, rise. In the real Eigen Times v0 basis the axes read bank · financial · credit · rate · loan, Russian · Russia · Ukraine · Putin · Moscow, flight · airline · passenger · airport · crash, and so on—sixty word lists that a reader recognises at once.

There are two distinct roles for a loading. When measuring an article’s coordinate, multiply each term entry by the direction’s loading for that term, then add the products. When reconstructing an article, multiply that whole direction by the article’s coordinate on it. A loading’s sign determines which pole of the direction it supports; its magnitude measures strength within this numerical representation. Neither alone is a probability that a document concerns a subject. A word list is a reading aid extracted from the direction, not a replacement for its full vector.

The small factorization can be computed directly. The real document–term table needs two additional ideas: applying a centered table without storing it, and finding its leading directions through a much smaller working space. We develop them separately.

5.3.1 Apply centering without storing a dense centered table

For this calculation, $X$ has $n$ document rows and $p$ term columns, $\mu$ is its $p\times1$ reference mean, and $\mathbf1$ has $n$ entries. Forming $X_c=X-\mathbf1\mu^{\top}$ would turn many absent-term zeros into negative mean values and destroy sparsity. We can still calculate a product involving $X_c$ without storing those changed cells.

Let $l$ count a small number of working columns, let $\Omega$ (capital omega) be any $p\times l$ matrix, and let $Y$ be $n\times l$. Each column of $\Omega$ is an input direction. Expand the centered product using the rule that multiplication distributes over subtraction:

$$X_c\Omega=(X-\mathbf1\mu^{\top})\Omega =X\Omega-(\mathbf1\mu^{\top})\Omega.$$

Associativity lets us change the grouping in the last product without changing its order. Compute the small $1\times l$ row $\mu^{\top}\Omega$ first, then repeat that row $n$ times using $\mathbf1$. For the transposed operation, apply the same reasoning. The resulting identities are

$$X_c\Omega=X\Omega-\mathbf1(\mu^{\top}\Omega),$$ $$X_c^{\top}Y=X^{\top}Y-\mu(\mathbf1^{\top}Y).$$

The results have shapes $n\times l$ and $p\times l$. In the second correction, $\mathbf1^{\top}Y$ adds each column of $Y$ over observations; multiplying those column totals by $\mu$ constructs the needed adjustment. Each correction is an outer product, of rank at most one, because all its columns are scalar multiples of the same column. The sparse $X$ is used as is; the centered dense matrix need never be stored. The notebook verifies both identities against a directly centered small table before using them in the larger algorithm.

5.3.2 Build a small working space

A $1{,}093{,}166\times50{,}000$ matrix is far too large for a full classical SVD here, but only the top 60 directions are wanted. The method of Halko, Martinsson and Tropp approximates them with repeated matrix products. Draw a random $50{,}000\times80$ matrix $\Omega$: $l=80$ working columns provide $k=60$ desired directions plus 20 extra directions for approximation. Each random column mixes many input features; applying $X_c$ gives one corresponding mixture of output directions. The resulting $n\times l$ table $Y=X_c\Omega$ is called a sketch: a smaller collection of vectors with which to explore the large matrix’s output space. The extra columns are oversampling, a margin for approximation rather than extra final news axes.

The range of a matrix is the set of outputs it can produce. The goal is for the sketch’s span to cover the most important part of that range. Random combinations can capture the leading directions well when the singular values decay sufficiently; the method offers an approximation whose quality must still be assessed.

5.3.3 Make the working directions perpendicular

To orthonormalise is to replace spanning vectors by perpendicular unit vectors spanning the same space. Here is one step of the Gram–Schmidt idea used to explain that operation. Take columns $a=(1,1)^{\top}$ and $b=(1,0)^{\top}$. First normalize $a$ to the unit column $q_1=a/\sqrt2$. The scalar $q_1^{\top}b=1/\sqrt2$ measures how much of $b$ already lies along $q_1$. The corresponding vector contribution is $q_1(q_1^{\top}b)=(1/2,1/2)^{\top}$.

Subtract that contribution from $b$. The residual is $(1/2,-1/2)^{\top}$ and has dot product zero with $q_1$. Normalize it to get $q_2=(1,-1)^{\top}/\sqrt2$. The columns $q_1,q_2$ now have unit lengths and zero mutual dot product, while their combinations still reach the same plane as $a,b$. If the residual were zero, the new input would add no independent direction and should be skipped. Practical routines also address floating-point accuracy, often with repeated orthogonalization or other stable factorizations. The notebook checks the small construction directly.

5.3.4 Strengthen leading directions, then factor the small table

Two rounds of power iteration amplify stronger singular directions relative to weaker ones. The update $Y\leftarrow X_c(X_c^{\top}Y)$ means replace $Y$ by that product. The left arrow denotes assignment: calculate the right side using the current $Y$, then store the result as the next $Y$.

Why does this help? A left singular direction with scale $\sigma_j$ receives a factor $\sigma_j$ when mapped back by $X_c^{\top}$ and another factor $\sigma_j$ when mapped forward by $X_c$. Its coefficient is multiplied by $\sigma_j^2$. In the rectangular diagonal example the two factors are 9 and 4, so the stronger contribution gains relative to the weaker. Orthonormalising between rounds prevents all the stored columns from merely growing huge and becoming numerically indistinguishable. The calculation applies two products in sequence; it does not form a huge $X_cX_c^{\top}$ matrix.

Call the final $n\times l$ orthonormal working matrix $Q_{\mathrm{range}}$. Its columns describe a smaller output space in which to seek the answer. The product $Q_{\mathrm{range}}^{\top}X_c$ expresses every original feature column in those working coordinates and has shape $l\times p$. Compute an SVD of that smaller table. Its left directions have $l$ entries; multiplying them by $Q_{\mathrm{range}}$ returns them to the original $n$-entry output space. Retain the desired $k$ directions. If a sketch has fewer than $l$ independent columns, the working dimension is reduced accordingly.

This gives an approximate truncated SVD of $X_c$. Weighted fitting applies the same method to $W^{1/2}X_c$. Each sparse product touches the nonzero entries once. On the reported laptop the first LSA fit, including tokenising 1.6 GB of text twice, takes about five minutes. A fixed random seed, input order, and numerical procedure make the computation reproducible. For teaching, the notebooks use the same explicitly specified four-column sketch in both languages, so differences in random-number generators cannot hide disagreements in the algebra.

6 From eigenvectors to named axes

6.1 Two ambiguities

Eigenvectors come with two kinds of freedom that the mathematics does not fix and a newspaper must.

6.1.1 Choose which end of an axis points positive

If $v$ is an eigenvector, so is $-v$: the same line pointing the other way. A score is measured by a dot product with the direction, so flipping the direction also flips its scores. The contribution to reconstruction is unchanged because $(-v)(-c)=vc$ for a scalar score $c$. Its positive and negative poles, the two ends of the line, exchange names. Choosing a sign gives the display a consistent convention without changing the represented observations.

Eigen Times fixes orientation using skewness, a measure of asymmetry in the score distribution. Here a distribution describes how frequently different score values occur, and a tail is the part extending toward unusually large or small scores. A central moment averages powers of deviations from the mean. We already know the second central moment: averaging squared deviations gives variance. The third central moment instead averages cubed deviations. Cubing retains signs, so a large negative deviation contributes negatively and a large positive deviation contributes positively.

Let $m_2$ and $m_3$ denote the second and third central moments on one axis. For the notebook’s scores −5, −1, 0, 1, 1, the mean is −0.8. The centered scores are −4.2, −0.2, 0.8, 1.8, 1.8. Squaring and averaging gives $m_2=4.96$; cubing and averaging gives $m_3=-12.384$. The negative third moment says that, by this signed cubic measure, the negative side is more pronounced.

When $m_2>0$, standardized skewness is $m_3/m_2^{3/2}$. The power $3/2$ means the standard deviation cubed: $m_2^{3/2}=(\sqrt{m_2})^3$. This division removes the units of the scores, but it does not change the sign. The orientation rule therefore needs only the sign of $m_3$: flip if $m_3<0$, making the third moment positive. Reversing all centered scores preserves their squares and reverses their cubes. This is a convention based on the measured third moment, not a claim that every asymmetric distribution has one unambiguous longer tail. A zero third moment supplies no preferred sign. The streaming chapter computes these moments without retaining every score.

6.1.2 Choose directions inside an already selected subspace

Equal eigenvalues are called degenerate: any perpendicular unit pair in their plane is an eigenbasis. For example, a covariance equal to twice the identity multiplies every direction by two. There is no uniquely preferred diagonal in its plane because every unit direction has variance two. The notebook rotates this covariance and verifies that its entries remain unchanged.

Close eigenvalues can make individual directions sensitive to small data changes, even when their shared subspace is stable (Figure 9). A perturbation is such a change to the covariance; the eigenvalue gap is the separation from neighbouring eigenvalues. The tracking chapter explains the bound relating them. Even distinct principal axes can be difficult to name: variance-maximising directions may blend subjects such as markets and party politics. Naming keeps the selected subspace and chooses different axes inside it; the rotated axes need no longer be covariance eigenvectors. The purpose of the new coordinates is readability, while the retained geometric information stays the same.

6.2 Varimax

6.2.1 Turn the axes and counter-turn the coordinates

Recall the retained $d\times k$ orthonormal basis $V$. Let $R$ be a $k\times k$ orthogonal rotation and let $V'=VR$ be the named basis; the prime here means rotated. If $c$ is a $k\times1$ raw coordinate column, its named column is $y=R^{\top}c$ (also sometimes written $c'$). Why does the coordinate transformation use the transpose? The coordinates must change oppositely to the basis so that the represented point stays fixed. Substituting $y$ and using $RR^{\top}=I$ gives

$$V'y=VR(R^{\top}c)=Vc.$$

Both bases therefore reconstruct the same projection. Subtracting that same projection from an observation leaves the same residual, so its length is unchanged. The covariance-adjusted distance called $T^2$ in §7 is also unchanged when its covariance is transformed consistently; it is not determined by the subspace alone.

6.2.2 Read the criterion as a variance calculation

Kaiser’s varimax criterion selects a rotation making each direction load strongly on some features and weakly on others. A criterion, or objective, is a numerical rule for comparing candidate choices. A larger value of this criterion indicates greater contrast in the sizes of loadings within the columns. It favors concentration but does not force entries to become exactly zero.

Let $p$ count the features, let $i=1,\ldots,p$ index them, and let $j=1,\ldots,k$ index axes. Write $\ell_{ij}$ for row $i$, column $j$ of the rotated loading table. Hold one column $j$ fixed. Square each loading, average those squared loadings, then measure how much the squared loadings vary around their mean. Squaring treats equally strong positive and negative loadings equally; it is concentration of magnitude that matters here.

We already derived the identity “variance equals mean square minus square of mean.” Apply it to the values $\ell_{ij}^2$. Their squares are $\ell_{ij}^4$, so this column’s variance is the average fourth power minus the square of the average second power. Add those variances over columns. That sequence of operations gives the criterion

$$\sum_{j=1}^{k}\left[\frac1p\sum_{i=1}^{p}\ell_{ij}^{4}-\left(\frac1p\sum_{i=1}^{p}\ell_{ij}^{2}\right)^2\right].$$

The fourth power appears because it is the square of a squared loading. It is not a new assumption that extreme feature values should be counted four times. For a two-feature column with loadings 1 and 0, squared loadings are 1 and 0, their mean is 0.5, and the variance of those squared loadings is $((1-0.5)^2+(0-0.5)^2)/2=0.25$. Spread the same total squared loading equally by using loadings $1/\sqrt2$ and $1/\sqrt2$. The squared loadings are now 0.5 and 0.5, so their variance is zero. The first column has greater concentration even though both columns have squared length one. In a rotation we cannot choose each column independently: the whole loading table must turn together. The notebooks compute the variance of squared loadings by directly subtracting each column’s mean and also by the displayed formula; the two methods agree.

Figure 10 shows six variables in two groups mixed by 30°. Their small secondary loadings make the computed maximizing rotation approximately $-33.3°$, rather than an exact reversal of the mixing angle. The rotation concentrates each variable on one axis without making its other loading exactly zero.

6.2.3 Rotate two columns at a time

Eigen Times sweeps over every pair of axes, rotating just that pair by its best angle. This two-coordinate operation is a Givens rotation. It leaves the other columns unchanged, reducing one step of a many-axis search to choosing one angle.

Let $\phi$ (Greek phi) denote that angle. Angles may be measured in degrees, with 360° in a full turn, or in radians, where angle is arc length divided by circle radius. One full turn is $2\pi$ radians, with $\pi$ the circle constant. On the unit circle, $\cos\phi$ is the horizontal coordinate and $\sin\phi$ the vertical coordinate after turning through $\phi$ from the positive horizontal direction. Their squared sum is 1 because the point lies on that circle. A two-dimensional rotation matrix is

$$R(\phi)=\begin{pmatrix}\cos\phi&-\sin\phi\\\sin\phi&\cos\phi\end{pmatrix}.$$

For one pair only, let $x,y$ be its $p\times1$ loading columns, distinct from article and score vectors elsewhere. Multiplying the loading table by $R(\phi)$ changes row $i$’s pair to

$$x_i\cos\phi+y_i\sin\phi,\qquad -x_i\sin\phi+y_i\cos\phi.$$

These two entries’ squared sum remains $x_i^2+y_i^2$. Rotation redistributes a row’s contribution between axes while preserving its total squared size.

6.2.4 Why the best angle contains a quarter of an angle

We can evaluate the criterion for many trial angles, but the two-column case also has a direct formula. Elementwise arithmetic treats each row separately. Define $u_i=x_i^2-y_i^2$ as the difference of the squared loadings and $v_i=2x_iy_i$ as twice their product. Define four scalar sums

$$A=\sum_i u_i,\quad B=\sum_i v_i,\quad C=\sum_i(u_i^2-v_i^2),\quad D=\sum_i2u_iv_i.$$

These local capital letters are scalars, not earlier matrices. The first two sums measure the column totals of the new quantities; the last two collect their squared and cross-product terms. Define centered combinations $C_*=C-(A^2-B^2)/p$ and $D_*=D-2AB/p$. The subscript star marks these adjusted scalar values.

Here is the route from the criterion to the angle. After rotation, the difference of the two squared loadings is $u_i\cos(2\phi)+v_i\sin(2\phi)$, obtained by expanding the two squared rotated entries and using the double-angle identities. In those identities, $\cos(2\phi)=\cos^2\phi-\sin^2\phi$ and $\sin(2\phi)=2\sin\phi\cos\phi$. Substituting this difference into the variance criterion and expanding once more leaves an angle-independent part plus

$$\frac{C_*\cos(4\phi)+D_*\sin(4\phi)}{4p}.$$

The appearance of $4\phi$ comes from squaring expressions that already involve $2\phi$. We can maximize the displayed part geometrically: it is a dot product between the fixed vector $(C_*,D_*)$ and the unit vector $(\cos(4\phi),\sin(4\phi))$, divided by a positive constant. The dot product is largest when the unit vector points along the fixed vector. The function $\operatorname{atan2}(a,b)$ returns the angle of the point with horizontal coordinate $b$ and vertical coordinate $a$, preserving the quadrant that a ratio alone would lose. It supplies $4\phi$; divide by four to get

$$\phi = \tfrac{1}{4}\,\operatorname{atan2}\left(D - 2AB/p,\; C - (A^2 - B^2)/p\right).$$

This is a formula for choosing a maximizing angle, not a requirement to memorize trigonometric manipulation. The notebooks check its objective against a dense grid of trial angles and verify the angle-dependent expression. If both arguments of $\operatorname{atan2}$ are zero, that expression is constant and the pair supplies no preferred rotation.

A sweep visits every axis pair. Sweeps stop when no pair rotates by more than $10^{-6}$ radians (one millionth of a radian), or after fifty sweeps. Exact pairwise maximization cannot decrease the criterion, so its bounded objective improves monotonically, meaning it never moves downward. This does not guarantee the global best orientation: the choices across many column pairs interact. The reported alternative, an SVD-based fixed-point iteration (repeatedly applying one update rule), oscillated between two mirror-image states on a symmetric configuration, motivating the pairwise method.

6.2.5 Give an embedding direction evidence from words

For the embedding basis, naming rotates associations with words rather than unexplained embedding coordinates. The articles supply the common rows connecting the two descriptions: each has term measurements and raw eigen-coordinates.

Let $T$ be the $N\times p$ article TF‑IDF matrix, $\mu_T$ its $p\times1$ article mean, and $\mathbf1$ a column of $N$ ones. Let $S$ be an $N\times k$ table of raw eigen-coordinates divided by their positive raw-axis standard deviations. Dividing a coordinate by its standard deviation makes one unit mean one reference standard deviation on every axis. Because the unrotated eigen-coordinates already have zero reference mean and zero reference cross-covariance, this scaling gives unit reference variance and zero reference cross-covariance. Together those properties are called whitening. Scaling arbitrary correlated coordinates to unit variance would not by itself remove their cross-covariances.

Let $\beta$ (Greek beta) be the $p\times k$ term-loading matrix. Its entry for one term and one axis sums, over articles, the centered term value multiplied by that article’s standardized raw score. Positive products arise when the term and score are both above their respective references, or both below. The sum measures association between term use and the direction; it is not a probability and has not been divided by the number of articles. Matrix multiplication performs all these sums together:

$$\beta=\left(T-\mathbf1\mu_T^{\top}\right)^{\top}S.$$

Here $p=50{,}000$ and $k=60$. Lexical whitening floors its raw variance denominators at $10^{-12}$ before taking square roots: if a variance is smaller, the implementation uses that positive minimum to avoid division by zero or an excessively small scale. Exact unit-variance whitening applies to positive axes whose denominators have not been changed by this safeguard, under the same reference weights used to fit them.

The embedding naming fit rotates $\beta$; the LSA fit rotates its right singular directions, already indexed by terms. Applying the resulting $R$ to the corresponding raw coordinate map gives named directions. Because scaling and rotation need not give the same result in opposite orders, word associations fitted from whitened scores are labels for investigation, not exact term-by-term decompositions of every stored named score. The block-sum chapter expands the centered loading formula. The older shorthand $T^{\top}U$ assumes compatible centered, scaled score columns; the SVD factor $U$ must not generally be substituted for $S$. The distinction is practical: a unit left singular direction, an article coordinate, and a variance-standardized coordinate carry different scales even when they point along related patterns.

7 Measuring a story

7.1 Projection, reconstruction, residual

The question is now about one observation: how much of it can the retained directions reconstruct, and what do they leave out? Here an observation is an article, or a representative vector chosen for a story. It is already a numerical vector; we are not yet assigning it a claim about historical novelty.

Let $d$ count its original coordinates and $k$ count retained directions. Write $x$ for the $d\times1$ observation column and $\mu$ for the $d\times1$ fitted article mean. The $d\times k$ matrix $V$ contains the retained raw directions as columns. They are orthonormal: each has length one and different columns have dot product zero. The following steps have different jobs.

7.1.1 First find coordinates, then reconstruct

Subtract the reference point. The centred column $x-\mu$ is the displacement from the fitted mean. Projection concerns this displacement, so the mean must be removed before measuring coordinates.

Measure along each retained direction. Define the $k\times1$ coordinate column $c$. Entry $c_j$, with $j=1,\ldots,k$, is the dot product with retained direction $j$. Stacking those dot products is matrix multiplication:

$$c=V^{\top}(x-\mu).$$

Rebuild what those coordinates describe. Multiplying by $V$ takes the weighted sum of the retained directions. Adding the mean puts the result back in the original coordinate system. Call this $d\times1$ reconstruction $\hat x$; the hat marks an approximation:

$$\hat x=\mu+Vc.$$

Subtract the reconstruction from the observation. The $d\times1$ difference $r$ is the residual: what remains after the attempted reconstruction. The word names a difference vector, not its length and not its probability:

$$r=x-\hat x.$$

For a complete synthetic example, retain the first two coordinate directions in three dimensions. Take $x=(2,1,3)^{\top}$ and $\mu=(0,0,0)^{\top}$. Then

$$V=\begin{pmatrix}1&0\\0&1\\0&0\end{pmatrix},\qquad c=\begin{pmatrix}2\\1\end{pmatrix},\qquad \hat x=\begin{pmatrix}2\\1\\0\end{pmatrix},\qquad r=\begin{pmatrix}0\\0\\3\end{pmatrix}.$$

The first coordinate survives unchanged, the second survives unchanged, and the third is absent from the retained plane. Subtraction reveals that absence explicitly: $3-0=3$ in the last residual entry.

7.1.2 Why this is the nearest reconstruction

An orthogonal projection keeps the component in the chosen subspace and leaves a residual perpendicular to it. Here this can be checked by taking its coordinates again. Since $V^{\top}V=I$, where $I$ is the $k\times k$ identity,

$$V^{\top}r=V^{\top}(x-\mu)-V^{\top}Vc=c-c=0.$$

Thus every retained direction has zero dot product with the residual. This is also why the reconstruction is a least-squares one. Any change to the reconstructed point within the retained subspace adds another perpendicular component to the error. Its squared length adds to the existing residual squared length, so it cannot improve the fit. No differentiation is needed for this geometric argument.

Recall the two distinct operations in a length calculation. Squaring entries and adding them gives squared length; taking its square root gives length. For our residual,

$$\lVert r\rVert^2=0^2+0^2+3^2=9,\qquad \lVert r\rVert=\sqrt9=3.$$

The centred observation has squared length $2^2+1^2+3^2=14$. The retained part has squared length $2^2+1^2=5$. The Pythagorean identity reads $14=5+9$: retained variation and residual variation account for the whole squared length. We use squared lengths below because perpendicular contributions add in this way. The ordinary lengths do not add: $\sqrt{14}$ is not $\sqrt5+3$.

7.2 T² and Q

We need two different comparisons. One asks how far the story lies along the retained directions, taking each direction’s usual spread into account. The other asks how much of the story lies outside those directions. A statistic is a numerical summary calculated from data; calling something a statistic does not by itself assign a probability to it.

7.2.1 Compare a coordinate with its usual spread

Recall that variance is an average squared deviation from a mean, and standard deviation is its nonnegative square root. Raw coordinate $j$ has fitted mean zero and variance $\lambda_j$. For this calculation retain only positive eigenvalues, so dividing by $\sqrt{\lambda_j}$ is possible. The quotient $c_j/\sqrt{\lambda_j}$ expresses the coordinate in units of its own standard deviation. A displacement of two units is ordinary on an axis whose usual spread is two, but larger relative to an axis whose usual spread is one.

A z-score subtracts a reference mean and divides by a positive reference standard deviation. The mean subtraction has already happened for the raw coordinates. Their variances can differ, but their fitted cross-covariances are zero. Dividing each by its own standard deviation therefore gives whitened coordinates: unit variance in every direction and zero cross-covariances under that same fitted reference.

Define $T^2$, Hotelling’s squared covariance-adjusted distance, by squaring these standardized coordinates and adding. Here $T^2$ is the statistic’s name, not the square of the TF‑IDF matrix $T$. With $j=1,\ldots,k$,

$$T^2=\sum_j\left(\frac{c_j}{\sqrt{\lambda_j}}\right)^2 =\sum_j\frac{c_j^2}{\lambda_j}.$$

Our synthetic story has $c=(2,1)^{\top}$. Suppose the fitted raw variances are $4$ and $1$, hence standard deviations $2$ and $1$. Its standardized coordinates are $(2/2,1/1)^{\top}=(1,1)^{\top}$. Squaring and adding gives

$$T^2=1^2+1^2=2.$$

Both retained coordinates are one reference standard deviation from their means, although their original sizes differ.

For compact matrix notation let $\Lambda$ be the $k\times k$ diagonal matrix of those positive raw eigenvalues. A superscript $-1$ denotes an inverse, the operation undoing multiplication by an invertible matrix. For a diagonal matrix, inversion replaces each nonzero diagonal entry by its reciprocal. Consequently $\Lambda^{-1}$ first divides coordinate $j$ by $\lambda_j$; the dot product with $c$ then multiplies by $c_j$ again and adds:

$$T^2=c^{\top}\Lambda^{-1}c.$$

7.2.2 Measure the part outside the retained space

Define $Q$ to be residual squared length:

$$Q=\lVert r\rVert^2.$$

For the same story, $Q=9$. This is a different measurement from $T^2=2$. A point can be far along a retained direction and have no residual at all; another can lie close to the mean within the retained directions but have a substantial perpendicular residual. A residual can contain discarded familiar variation as well as genuinely new material. Neither statistic alone proves that an event is unprecedented.

Figure 12 shows why reference spreads matter: equal Euclidean distances can have different covariance-adjusted distances. With 60 fitted positive directions, the fitted average of $T^2$ under the matching covariance weights is 60. To see why, average each term $c_j^2/\lambda_j$: the numerator’s fitted average is precisely $\lambda_j$, so every term contributes one. A value of 400 is much larger than that fitted average. This arithmetic does not yet provide the probability of reaching 400; that would require a justified probability model or a separately assessed reference distribution.

7.2.3 Keep the covariance when changing coordinates

The calculation uses raw eigen-coordinates before the naming rotation. The rotation changes the coordinate descriptions without changing the retained point. Recall the $k\times k$ orthogonal naming rotation $R$. In the named coordinate column $y=R^{\top}c$, define the covariance $\Gamma=R^{\top}\Lambda R$. Its off-diagonal entries need not be zero: named coordinates can vary together even though the raw eigen-coordinates did not.

To measure the same distance in this new coordinate system, use its full covariance:

$$T^2=y^{\top}\Gamma^{-1}y.$$

The inverse now adjusts for both different spreads and shared movement between coordinates. Reusing the old diagonal eigenvalues would attach old scales to new directions. Dividing by the named coordinates’ own variances is also insufficient when their covariances are nonzero. That operation is marginal standardization: it gives each coordinate unit variance separately but leaves correlations. Full whitening removes cross-covariances as well. The notebooks rotate the synthetic story and verify that the full calculation stays unchanged while the diagonal shortcut generally changes it.

The implementation sums $c_j^2/\lambda_j$ only for positive raw eigenvalues and contributes zero otherwise. This is an explicit handling rule for zero-variance directions, not permission to divide by zero. The inverse formulas above assume positive retained eigenvalues.

7.2.4 Turn residual size into a share

The two statistics originated in monitoring variation against a principal-component reference model. For browsing, it is also useful to ask what share of a centred observation remains outside the retained subspace. Let $\nu$ (Greek nu) denote this novelty ratio. When $x\ne\mu$, meaning the observation differs from the fitted mean, define

$$\nu=\frac{Q}{\lVert x-\mu\rVert^2}.$$

The denominator counts all centred squared length; the numerator counts only the residual part. In our example, $\nu=9/14\approx0.643$. About 64.3% of this observation’s squared length lies outside the retained plane. That is a statement about this representation, not a 64.3% probability of a historically new event.

Orthogonality makes the denominator the sum of retained and residual squared lengths, both nonnegative. The ratio therefore lies between zero and one. Zero means the centred vector is fully represented; one means it lies perpendicular to the retained subspace. At $x=\mu$, both numerator and denominator are zero, so the ratio is undefined. The implementation returns zero in that case as a software convention. This geometric use of “energy” means squared length; the attention score below uses the same ordinary word for a different quantity.

7.2.5 Average measurements in the order intended

These equations describe one vector. The site’s stored story statistics are means of the member articles’ $T^2$, $Q$, and $\nu$ values. The order matters: measure every member, then average the measurements. Measuring the mean embedding answers a different question.

For example, two member coordinate columns $(2,0)^{\top}$ and $(-2,0)^{\top}$ under variances $(4,1)$ each have $T^2=4/4=1$, so their mean statistic is one. Their mean coordinate column is zero, whose $T^2$ is zero. Averaging first cancels their opposite departures. Squared distances and ratios generally do not commute with averaging; that is why an article statistic, an average article statistic, and a centroid statistic must be kept distinct.

7.3 Dominant axis, spectrum, profile

These three objects answer questions at three different levels. The dominant axis selects a story’s label. The spectrum combines the day’s stories. A profile describes the stories usually carrying one label.

7.3.1 Select a label using comparable signed scores

For story $s$, let $y_s$ be its $k\times1$ named-coordinate column, with entry $y_{sj}$ on axis $j$. Here $s$ identifies a story. The separate notation $s_j>0$ denotes the empirical standard deviation of named coordinate $j$ across stories in the profile window: the latest ten years of loaded data. This is a reference scale, estimated separately from the raw eigenvalue $\lambda_j$. The implementation replaces scales smaller than $10^{-6}$ by that positive floor to avoid division by zero or an extremely small denominator.

Divide each named coordinate by its reference scale. The dominant archetype is the axis with the largest signed result. The notation $\arg\max_j$ means “return the index attaining the largest value”, whereas $\max_j$ would return the value itself. Thus the rule is $\arg\max_j y_{sj}/s_j$.

Take the notebook’s synthetic story $(2,1)^{\top}$ and scales $(2,0.5)$. The scores are $2/2=1$ and $1/0.5=2$, so the second axis wins despite its smaller raw coordinate. For a second story $(-1,-0.1)^{\top}$, the scores are $-0.5$ and $-0.2$. The second wins again because $-0.2$ is larger. Taking absolute values would change that decision. This particular rule does not subtract an empirical axis mean, so it is scale division rather than a general mean-subtracted z-score.

One possible “mixed” heuristic compares the positive top score $a_1$ with the runner-up $a_2$. Their relative gap is $(a_1-a_2)/a_1$; the denominator is the top score. A rule declaring a mixture when that gap is at most $0.20$ would accept scores $1$ and $0.85$, whose gap is $0.15$. This is an illustrative rule retained from the original explanation, not a rule implemented by the current site.

7.3.2 Combine stories into a daily spectrum

A simple average gives every story equal influence. A weighted mean first multiplies each value by its nonnegative weight, adds those products, and divides by the positive total weight. Division by the total makes the result an average rather than a total that grows automatically when more stories are present.

Let $\varepsilon_s$ denote story $s$’s attention weight, called story energy in the site:

$$\varepsilon_s=\text{distinct sources}\times\ln(1+\text{article count}).$$

The natural logarithm grows slowly; adding one makes the argument positive even for a zero count. This weight uses source and article counts. It is distinct from article fit weights, eigenvalues, and vector squared length.

For one chosen day, let $\bar y_{\mathrm{day}}$ be its $k\times1$ spectrum column. The bar indicates averaging, and the sum includes that day’s stories only:

$$\bar y_{\mathrm{day}}=\frac{\sum_s\varepsilon_s y_s}{\sum_s\varepsilon_s}.$$

The notebooks use three synthetic story rows $(2,1)$, $(-1,-0.1)$, and $(0.5,3)$, with source counts $2,1,3$ and article counts $3,1,5$. Their attention weights are $2\ln4\approx2.773$, $\ln2\approx0.693$, and $3\ln6\approx5.375$. For the first spectrum coordinate, multiply $2,-1,0.5$ by those weights, add, and divide by their sum. Repeating for the second coordinate gives approximately $(0.853,2.130)$. The more heavily weighted third story pulls the mixture towards its second coordinate; no new direction has been fitted in this averaging step.

7.3.3 Compare a day, or compare a story

The era statistics are each spectrum coordinate’s mean and standard deviation over all loaded day spectra. The algorithm named after Welford updates three running summaries: how many values have arrived, their current mean, and the sum of squared deviations from that current mean. It adjusts the deviation sum when the mean changes; simply squaring deviations from an outdated mean would give the wrong answer. The statistical-threshold chapter walks through its calculation.

Each spectrum z-score subtracts the coordinate’s era mean and divides by its positive era standard deviation. The front page highlights values beyond $\pm2$ (plus or minus two) and lists values below $-2$ under Silence. These are comparisons with all loaded data, including dates later than some selected historical dates. They are not historical forecasts computed using only the past available at the time.

An archetype’s profile uses a different reference set: the stories it dominates in the profile window. For each named coordinate, it records a mean and standard deviation among those stories. Same as always combines strong positive profile and story signals with agreement within one profile standard deviation; agreement alone is insufficient. Different this time lists departures exceeding 1.5 profile standard deviations. The symbol $\sigma$ in these descriptions is shorthand for the relevant reference standard deviation, not an SVD singular value. A zero reference deviation needs an explicit fallback before any z-score can be computed. Always identify the reference population before interpreting a standardized score.

8 The clustering thresholds

A threshold turns a numerical measurement into a decision: link two articles if their similarity reaches a chosen value; highlight an axis if its standardized departure exceeds a chosen size. The measurement and the decision are separate. Cosine supplies a similarity, but mathematics alone does not say which cosine deserves a link.

The thresholds described here use cosines, z-scores, and counts. This chapter shows the empirical distributions used to choose them for the dated Eigen Times corpus. The methods also apply to Eigen Hacks, but these recorded newspaper measurements are not measurements of its technology-community corpus. A changed embedding model, input length, language, or corpus can change the distribution enough to require a fresh assessment.

8.1 Stories: single linkage on a day’s articles

Start with a graph: a collection of nodes and connections called edges. Here each article is a node. Compare each pair of that day’s article vectors; add an undirected edge when their cosine reaches the chosen threshold. Undirected means that a link from the first article to the second is also a link back.

Next follow links, including indirect links. A path is a sequence of edges joining successive nodes. A connected component is a maximal group joined by paths: no linked node has been left outside it. Each component becomes one story. An isolated article is a singleton, a component of size one.

This rule is called single linkage because one eligible link is enough to join groups. It does not require every pair in a story to meet the threshold. A union–find data structure implements that decision economically: each node begins in its own group; whenever an edge joins two groups, their group records are merged. It is a bookkeeping method for the same connected components, not an additional similarity rule.

The notebooks show the resulting chaining with four synthetic unit directions at 0, 30, 60, and 90 degrees. Adjacent directions have cosine about 0.866, so a threshold of 0.72 joins all four through a chain. The endpoints have cosine zero but still belong to one component. At 0.9, none of those synthetic pairs link. Figure 13 illustrates the danger: on a busy day, intermediate similarities can connect unrelated endpoints and let one chain absorb much of the news.

Figure 14 examines the cosines of article pairs on five busy days (12 September 2001, 16 September 2008, 24 June 2016, 12 March 2020, 24 February 2022), split by whether the two ended up in the same story. A distribution describes how values are spread across possible sizes. A histogram groups them into intervals, or bins, and counts values in each interval. The sample’s median is its middle value; its 99th percentile is a value at or below which about 99% of observations lie. The figure scales each histogram to its own highest bin, so equal displayed heights do not mean equal pair counts or equal probabilities. These are descriptive same-story labels from the grouping, not an independent human-labeled test set. In embedding space, pairs from different stories have median cosine 0.49 and a 99th percentile of 0.67—confirming that unrelated news embeddings sit around 0.5, not 0—while pairs from the same story have median 0.65 and a long tail toward 1. The two distributions overlap, as they must for anything less than a perfect representation, but the crossover is narrow.

Choosing a threshold balances two failures. Lowering it adds edges, which can join distinct events through chains. Raising it removes edges, which can split one event into several components; this is fragmentation. Figure 15 shows the search for a useful balance on 24 February 2022, the day Russia invaded Ukraine. At 0.6, one chain contains 143 of the day’s 152 articles—most of the day has become one “story”. At 0.66 the chain still holds 104. At 0.72 the largest cluster is the 62‑article invasion story, the intended grouping identified during inspection, with 67 singletons (features, comment, unrelated items). At 0.78 the invasion story has fragmented to 34 articles; at 0.82 no cluster is larger than 8; at 0.9 only one linked pair remains: 151 components, the largest of size 2, and 150 singletons among 152 articles. These are the counts in the frozen aggregate fixture; the earlier wording that nothing linked was incorrect. The deployed threshold is 0.72. The sweep shows the effect of this choice on that day’s graph; it does not prove an optimal threshold for every possible day. The original choice had been 0.82 (a number appropriate for a larger embedding model with longer inputs); the first look at the data—“Day of terror” on 12 September 2001 rendered as a three‑article story—showed it was wrong, and the sweep in Figure 15 replaced it.

The lower panel of Figure 14 explains why the LSA space needs a different threshold. LSA coordinates are projections onto the 60 leading right singular directions of the term matrix; articles about the same event share vocabulary but in a 60‑dimensional summary their cosines spread from 0 to 1 (median 0.35), and unrelated articles sit around 0 (median $-0.01$). The distributions are wider and cross lower, and the threshold that reproduces sensible stories is 0.6. The same numerical threshold therefore has a different meaning in the two representations. Their coordinate systems and similarity distributions differ. This is one reason the deployed newspaper clusters in embedding space and uses the LSA basis only for measurement. The notebooks can recompute the histogram scaling and proportions from frozen aggregate counts; they cannot recreate historical cluster assignments without the original article vectors.

8.2 Episodes: threading across days

An episode is a chain of related stories across days. First summarize each story by its centroid: average its member embedding vectors, then use that mean’s direction in the cosine comparison. Averaging is done coordinate by coordinate. A zero mean has no direction, so an implementation must handle that case explicitly rather than calculate an undefined cosine.

In this subsection $t$ indexes calendar days, distinct from the earlier term index. To thread a story on day $t$, compare its centroid with each story centroid on day $t-1$. Keep the greatest cosine, the best match. If it reaches 0.9, continue that previous story’s episode; otherwise start a new episode. This is a best-previous-story rule, not a one-to-one assignment between all stories on the two days. Several current stories can independently select the same previous episode.

The threshold acts after the maximum is found. In the notebook’s synthetic two-story example, one current centroid has best similarity 0.92 and continues; the other fails the 0.9 threshold and starts a new episode. A best match always exists when there are eligible nonzero previous vectors, but it need not be a good match. Figure 16 shows, for 678 stories over eight day-pairs, each best previous-day cosine. With raw embeddings the median is 0.66 because unrelated news shares a common component. Reported continuations above the strict threshold include “Battle for Kyiv” after “Putin invasion deepens” (0.92), the Hamas attack’s death toll two days running (0.90), and Brexit voters’ portraits (0.93). In late February 2022, 21 of 1,077 stories continued an episode, and the longest, the invasion, ran fourteen days from 18 February to 3 March.

Subtracting the mean first removes the common component and lowers the median to 0.29, with a sparser region around 0.55–0.75. A centred threshold near 0.75 is a candidate for separate evaluation, not a direct replacement proved by the figure. It would still miss the cited Crimea continuation at 0.69; lowering it enough to include that pair could also admit the unrelated New Zealand protest and abortion pair at 0.70. To assess this tradeoff, first obtain judgments of which candidate links really are continuations. Precision is the number of correct accepted links divided by all accepted links. Recall is the number of correct accepted links divided by all true continuations in the evaluated set. If no link is accepted, precision has a zero denominator; if the evaluated set has no true links, recall does too. Those cases need explicit reporting conventions.

The notebooks supply five labeled synthetic candidates: three true links at 0.92, 0.90, and 0.93, a false link at 0.70, and a true link at 0.69. A 0.9 threshold accepts three correct links: precision is $3/3=1$ and recall is $3/4=0.75$. Lowering it to 0.69 recovers all four true links but also accepts the false one: precision becomes $4/5=0.8$ and recall $4/4=1$. Those labels are teaching inputs, not a new evaluation of the historical stories. The strict raw rule aims to favor precision. An episode is a thread, not a day: “episode day 7” counts from the chain’s beginning, and an archetype’s canon counts episodes rather than repeated days.

8.3 Precedents

A precedent is an earlier candidate presented for comparison, not proof that history will repeat. The search has an eligibility stage and a ranking stage.

First retain only stories earlier than the query story and having the same dominant axis. This prevents future stories from entering a historical comparison and makes the archetype part of the search definition. Then compare named-coordinate vectors by cosine. Among stories from the same episode, keep the strongest match; finally take the five strongest remaining matches, or fewer if fewer are eligible. The distinct-episode restriction avoids filling the list with five days of the same event.

No similarity threshold is required to return the nearest eligible candidates. This is a ranking: the largest scores come first. A ranking can still return weak matches if all eligible candidates are weak. Printing the cosine beside each precedent lets the reader assess that distinction: recorded examples include 0.97 for the same week’s build-up briefing and 0.93 for Crimea in 2014. A shared dominant axis or a large cosine establishes resemblance in this representation, not a causal relationship between the events.

9 The statistical thresholds

9.1 ±2σ on the spectrum

The symbol $\sigma$ in this heading means a reference standard deviation. For each named axis, collect one spectrum value per loaded day. Subtract their mean to describe departures from the usual level, then divide by their positive standard deviation to put those departures on a common scale. The unit is now “reference standard deviations”, rather than the original coordinate unit. An axis is loud above $+2$ and silent below $-2$.

9.1.1 Calculate a reference spread before a z-score

Use the synthetic daily values $-1,0,1,2,3$. Their mean is $1$. Their deviations are $-2,-1,0,1,2$, whose squares add to $4+1+0+1+4=10$. The notebook uses the sample variance: divide by the count minus one, giving $10/(5-1)=2.5$. Its standard deviation is $\sqrt{2.5}\approx1.581$. This count-minus-one convention differs from the total-weight normalization of the fitted covariance earlier; name the convention instead of silently switching denominators.

Against this reference, a value $-3$ first gives deviation $-3-1=-4$, then z-score $-4/\sqrt{2.5}\approx-2.530$. It belongs to Silence. Neither $-3$ nor $-4$ alone could be compared directly with the standardized cutoff $-2$; dividing by the reference spread is essential.

To obtain the same summaries one value at a time, let $n$ be the count after a new scalar value $a$ arrives. Let $h_{\mathrm{old}}$ denote the previous running mean and $A_{\mathrm{old}}$ the previous sum of squared deviations from that mean. Define $\delta=a-h_{\mathrm{old}}$, the new value’s departure from the old mean. The updated mean $h_{\mathrm{new}}$ moves by that departure divided among all $n$ values:

$$h_{\mathrm{new}}=h_{\mathrm{old}}+\delta/n.$$

The updated squared-deviation sum is

$$A_{\mathrm{new}}=A_{\mathrm{old}}+\delta(a-h_{\mathrm{new}}).$$

These are Welford’s updates. Initialize the first mean at the first observation and its deviation sum at zero. Each later update uses both the old and new mean; this compensates for the movement of the reference point. With at least two values, divide the final sum by $n-1$ for the sample variance. With only one value, that sample variance is undefined. The notebooks verify these updates against the direct five-value calculation.

9.1.2 A standardized cutoff is not automatically a probability

A probability model specifies how frequently values are expected to fall in different ranges. A Gaussian, or normal bell-shaped distribution, is one such model. Under that model, about 4.55% of observations lie outside two standard deviations, often rounded to “one in twenty”. News series can be heavy-tailed, meaning that large departures occur more often than this model predicts. Standardization alone does not make them Gaussian.

There is also a many-comparisons issue. If 60 axes each really had the Gaussian 4.55% tail frequency, the expected number flagged would be about $60\times0.0455=2.73$. An expected number is a long-run average count, not the probability that any particular day has a flag. This average-count calculation does not require independent axes; independence would be an additional assumption when calculating the probability of at least one flag. The notebooks compute the Gaussian benchmark, not a calibrated newspaper alarm probability.

The empirical display intent is that a typical day lights one to three axes, not ten. On 24 February 2022 the Ukraine axis was at $+8$; on 12 September 2001 terrorism was at $+2.9$, Europe at $+2.6$, markets at $+2.5$. These recorded numbers are spectrum z-scores against the loaded reference, not probabilities.