A Working Mathematics Reference

$\sin^2\theta+\cos^2\theta=1$ is the unit circle restated — the coordinates of a radius-1 point. Divide it through by $\cos^2\theta$ to get $1+\tan^2\theta=\sec^2\theta$; divide by $\sin^2\theta$ for $1+\cot^2\theta=\csc^2\theta$. You only need to remember the first.

The two angle-addition formulas are the parents of nearly every other identity:

$$\sin(\alpha\pm\beta)=\sin\alpha\cos\beta\pm\cos\alpha\sin\beta$$ $$\cos(\alpha\pm\beta)=\cos\alpha\cos\beta\mp\sin\alpha\sin\beta$$

Memorise these two and you can derive the rest rather than memorising a wall of formulas.

The double-angle formulas are just angle-addition with $\beta=\alpha$. Setting the two angles equal in the parent identities:

$$\sin 2\theta=2\sin\theta\cos\theta,\qquad \cos 2\theta=\cos^2\theta-\sin^2\theta.$$

Then $\sin^2\theta+\cos^2\theta=1$ rewrites the cosine version two more ways — $\cos 2\theta=1-2\sin^2\theta=2\cos^2\theta-1$ — which are exactly the forms you invert to integrate $\sin^2$ and $\cos^2$. Nothing was memorised beyond the two parents.

Inverse trig has restricted ranges. $\arcsin:[-1,1]\to[-\tfrac\pi2,\tfrac\pi2]$, $\ \arccos:[-1,1]\to[0,\pi]$, $\ \arctan:\mathbb R\to(-\tfrac\pi2,\tfrac\pi2)$. The restricted ranges are the price of inverting a periodic (hence not one-to-one) function — feed $\arcsin$ anything outside $[-1,1]$ and it is simply undefined. These are the functions whose derivatives appeared in Part I.

Radians, not degrees. Every calculus identity — derivatives of $\sin$/$\cos$, the small-angle approximations, Taylor series — assumes radians. The famous $\frac{d}{dx}\sin x=\cos x$ is false if $x$ is in degrees.

Euler's formula lets a single complex exponential carry both growth and rotation at once. Write an eigenvalue as $\lambda=\alpha+i\beta$; then

$$e^{\lambda t}=e^{\alpha t}\,e^{i\beta t}=e^{\alpha t}\big(\cos\beta t+i\sin\beta t\big).$$

Read it straight off: the real part $\alpha$ sets growth ($\alpha>0$) or decay ($\alpha<0$), and the imaginary part $\beta$ sets the oscillation frequency. A purely imaginary $\lambda$ ($\alpha=0$) gives undamped oscillation — a center. This is the entire content of "complex roots $\Rightarrow e^{\alpha t}(\cos\beta t,\sin\beta t)$" from Part VI, and of spirals, centers, and the Hopf bifurcation in Part VII.

Rearranging Euler's formula for $+\theta$ and $-\theta$ and adding/subtracting expresses the trig functions through exponentials:

$$\cos\theta=\frac{e^{i\theta}+e^{-i\theta}}{2},\qquad \sin\theta=\frac{e^{i\theta}-e^{-i\theta}}{2i}.$$

This is the bridge that turns Fourier analysis (Part VIII) from trigonometry into algebra: a sum of sinusoids becomes a sum of complex exponentials, and convolution becomes multiplication.

$\alpha$ and $\beta$ play different roles — don't swap them. In $e^{\lambda t}$ the real part is the rate (stability) and the imaginary part is the frequency. Reading them the wrong way round turns a stable spiral into a non-existent "fast decay."

A purely imaginary eigenvalue is the knife-edge. $\alpha=0$ is a center in the linearization, but the smallest nonlinear term can tip it into an inward or outward spiral — this is exactly where linear stability analysis stops being trustworthy (Part VII).

The definition is a challenge–response game. An adversary names a tolerance $\varepsilon$ — "get $f$ within this band of $L$." You must answer with a neighbourhood radius $\delta$ — "stay within this much of $a$ and I guarantee it." A limit exists precisely when you can always answer, no matter how small $\varepsilon$ gets. The value $f(a)$ never enters, which is why a function can have a clean limit at a point where it is undefined or even mis-defined.

L'Hôpital's rule handles the indeterminate cases $\tfrac00$ and $\tfrac{\infty}{\infty}$: $\displaystyle\lim_{x\to a}\frac{f(x)}{g(x)}=\lim_{x\to a}\frac{f'(x)}{g'(x)}$ when the right limit exists. It works because near $a$ each function is well-approximated by its tangent line, so the ratio of values becomes the ratio of slopes — a first-order Taylor argument in disguise.

At $x=0$ the expression is literally $\tfrac00$ — undefined. Yet the limit is exactly $1$, and you can see it three ways. Geometrically, on the unit circle the arc length $x$, the chord $\sin x$, and the tangent $\tan x$ are squeezed together as $x\to0$. Via Taylor, $\sin x = x - \tfrac{x^3}{6}+\cdots$, so $\frac{\sin x}{x}=1-\tfrac{x^2}{6}+\cdots\to1$. Via L'Hôpital, differentiate top and bottom: $\frac{\cos x}{1}\to\cos 0=1$. All three agree — this single limit is what makes $\frac{d}{dx}\sin x=\cos x$ true in the first place.

L'Hôpital only applies to genuine indeterminates. $\tfrac00$ and $\tfrac{\infty}{\infty}$ qualify; a determinate form like $\tfrac{1}{0}$ does not, and differentiating it gives nonsense. Check the form before reaching for the rule.

Continuity needs all three conditions. A limit existing is not enough — $f(a)$ must also be defined and equal to it. A removable discontinuity is exactly the case where the limit exists but the function value disagrees or is missing.

A one-sided mismatch kills the limit. If the left and right limits differ, the two-sided limit does not exist, no matter how nicely each side behaves on its own.

Why the chain rule works: rates multiply, because the outer function sees its input already scaled by the inner function's rate. Over a tiny step the inner map stretches the input by $g'(x)$, and the outer map then stretches that by $f'(g(x))$; the two stretches compose. That nesting is the seed of backpropagation, and it trips up beginners constantly.

Second derivative = concavity: $f''>0$ means the slope is increasing, so the graph is bowl-shaped (convex); $f''<0$ is dome-shaped (concave). This sign is concavity; "curvature" proper is the related but distinct quantity $\kappa=|f''|/(1+f'^2)^{3/2}$ — for the min-vs-max test only the sign matters. This is exactly what the Hessian tests in Part V to tell a minimum from a maximum.

This is a composition $f(g(x))$ with outer $f(u)=e^u$ and inner $g(x)=\sin x$. The chain rule says multiply the outer derivative (evaluated at the inner) by the inner derivative:

$$\frac{d}{dx}e^{\sin x}=\underbrace{e^{\sin x}}_{f'(g(x))}\cdot\underbrace{\cos x}_{g'(x)}.$$

The $\cos x$ is the part beginners drop. Nest one level deeper — $e^{\sin(x^2)}$ — and you simply multiply by one more inner rate: $e^{\sin(x^2)}\cos(x^2)\cdot 2x$. Backpropagation is nothing more than this bookkeeping carried through a whole computational graph.

Don't drop the inner derivative. $\frac{d}{dx}f(g(x))=f'(g(x))\,g'(x)$ — the $g'(x)$ factor is the single most common broken-backprop bug. Differentiating $\sin(2x)$ as $\cos(2x)$ instead of $2\cos(2x)$ is the canonical slip.

Concavity is the sign of $f''$, not its size. The second-derivative test tells a minimum ($f''>0$) from a maximum ($f''<0$) by sign alone; reaching for the full curvature formula $\kappa$ is unnecessary and invites arithmetic errors.

$a^x$ is not $x^n$. The power rule differentiates a variable base; for a variable exponent you need $\frac{d}{dx}a^x=a^x\ln a$. Mixing the two — writing $\frac{d}{dx}2^x=x2^{x-1}$ — is a frequent reflex error.

Demand a polynomial $P(x)=\sum_k c_k (x-a)^k$ that agrees with $f$ in value and in every derivative at $x=a$. Differentiate $P$ exactly $k$ times and set $x=a$: every term below degree $k$ has vanished, every term above still carries a factor of $(x-a)$ and dies, and the degree-$k$ term leaves $k!\,c_k$. Matching that to $f^{(k)}(a)$ forces $c_k=f^{(k)}(a)/k!$. The factorial is not decoration — it is exactly the derivative of $x^k$ taken $k$ times.

The remainder $R_n$ measures what the first $n$ derivatives cannot see; Taylor's theorem bounds it by the size of the next derivative on the interval, which is why the approximation is trustworthy precisely when higher derivatives stay tame.

Let $f(x)=e^x$ and $a=0$. Every derivative of $e^x$ is $e^x$, so $f^{(k)}(0)=e^0=1$ for all $k$. Drop that into $c_k=f^{(k)}(0)/k!$:

$$e^x=\sum_{k\ge0}\frac{1}{k!}x^k=1+x+\frac{x^2}{2}+\frac{x^3}{6}+\frac{x^4}{24}+\cdots$$

Check it at $x=1$: the partial sums $1,\,2,\,2.5,\,2.6\overline{6},\,2.708\ldots$ climb toward $e=2.71828\ldots$, and because the terms $1/k!$ shrink faster than any geometric ratio, the series converges for every $x$ — the radius of convergence is infinite. Sine and cosine are the same construction with the derivative cycle $\sin\to\cos\to-\sin\to-\cos$ stamping the alternating signs onto the odd and even powers.

A series can converge only near $a$. $\frac{1}{1-x}=\sum x^k$ holds only for $|x|<1$; push past the radius of convergence and the partial sums diverge even though $f$ itself is perfectly finite there. Always ask where a Taylor series is valid, not just what it is.

Convergence to the series $\neq$ convergence to the function. The function $f(x)=e^{-1/x^2}$ (with $f(0)=0$) is smooth, and every derivative at $0$ is zero. Its Maclaurin series is therefore identically $0$ — yet $f$ is positive everywhere else. The series converges; it just converges to the wrong thing. Smoothness alone does not buy you a usable Taylor expansion (that property is called analyticity).

Don't drop the inner derivative. Expanding $f(g(x))$ about a point, the linear term is $f'(g(a))\,g'(a)\,(x-a)$ — the $g'(a)$ is the chain-rule factor. Omitting it is the single most common broken-linearization (and broken-backprop) bug.

The FTC says the accumulation-so-far function $A(x)=\int_a^x f$ has derivative $f$. The reason is direct from the limit: $A(x+h)-A(x)$ is the area of a thin sliver of width $h$ and height $\approx f(x)$, so the difference quotient $\frac{A(x+h)-A(x)}{h}\to f(x)$. Integration and differentiation are therefore inverse operations — that single fact turns area problems into antiderivative lookups.

The two main techniques are the integral forms of the two main derivative rules, read backwards. Substitution is the chain rule reversed: it undoes the inner-derivative factor $g'(x)$ that the chain rule produced. Integration by parts is the product rule reversed: integrate $(uv)'=u'v+uv'$ and rearrange.

Pick $u$ by LIATE — the factor that gets simpler when differentiated. Here $u=x$ (algebraic, becomes $1$) and $dv=e^x\,dx$ (so $v=e^x$):

$$\int x e^x\,dx = uv-\int v\,du = x e^x-\int e^x\,dx = x e^x - e^x + C = (x-1)e^x+C.$$

The point of the LIATE choice is that $\int v\,du$ is easier than the original $\int u\,dv$ — picking $u=e^x$ instead would have made the remaining integral worse, not better.

The power rule has a hole at $n=-1$. $\int x^n\,dx=\frac{x^{n+1}}{n+1}$ would divide by zero when $n=-1$; that case is the separate fact $\int\frac1x\,dx=\ln|x|+C$. The absolute value matters — the antiderivative is valid on both sides of the origin.

Don't forget the constant (or the limits). An indefinite integral carries $+C$; a definite one requires evaluating $F(b)-F(a)$, and under a substitution the limits must be converted to the new variable (or back-substituted before plugging in).

Start with a basis of the kernel — the $\dim\ker(A)$ input directions that $A$ collapses to zero. Extend it to a full basis of $\mathbb R^n$; the vectors you added, $\operatorname{rank}(A)$ of them, map to an independent set spanning the column space. So the $n$ input dimensions split cleanly into the ones that die ($\ker$) and the ones that survive as output (rank). Nothing is created or destroyed — every input direction is accounted for exactly once, which is the whole content of the theorem.

Take rows $(1,2,3),(2,4,6),(0,1,1)$. Row 2 is exactly $2\times$ row 1, so it contributes no new direction; rows 1 and 3 are independent. The rank is $2$, not $3$ — the matrix maps $\mathbb R^3$ onto a 2D plane, and by rank–nullity its kernel is $3-2=1$ dimensional. Reducing to row-echelon form makes this mechanical: the rank is the number of non-zero pivot rows that remain.

Row rank from the original rows, not by eyeballing. Vectors can look independent yet hide a linear relation; only reduction (or a rank computation) settles it. Two "different-looking" rows that are scalar multiples add nothing.

Full column rank $\neq$ invertible unless square. A tall $m\times n$ matrix ($m>n$) can have rank $n$ (independent columns) yet have no inverse — invertibility is a square-matrix property. What it does have is a left inverse / least-squares solution.

Projection. The best approximation of $\mathbf x$ inside a subspace is the point whose error $\mathbf x-P\mathbf x$ is orthogonal to the subspace — any tilt away from perpendicular could be shortened. Writing that orthogonality condition against an orthonormal basis $U$ collapses to $P=UU^\top$, with the inner products $U^\top\mathbf x$ as the coordinates.

Why $\ell_1$ promotes sparsity. It penalises every non-zero entry equally regardless of size, so the optimiser prefers to zero components out entirely rather than keep many small ones. Geometrically the $\ell_1$ ball is a spiky diamond whose corners lie on the axes, and a shrinking constraint set first touches the loss at a corner — i.e. at a sparse point.

Project $\mathbf x=(3,4)$ onto the unit vector $\mathbf u=(1,0)$. The coordinate is $\mathbf u^\top\mathbf x=3$, so $P\mathbf x=(\mathbf u^\top\mathbf x)\,\mathbf u=(3,0)$ and the residual $(0,4)$ is orthogonal to $\mathbf u$ — as it must be. Generalise to a non-unit direction $\mathbf a$ and the same logic gives $\frac{\mathbf a^\top\mathbf x}{\mathbf a^\top\mathbf a}\,\mathbf a$, the least-squares projection that underlies linear regression.

$P=UU^\top$ needs orthonormal columns. For a general (non-orthonormal) basis $A$ the projector is $A(A^\top A)^{-1}A^\top$ — drop the middle factor and you get the wrong answer. Orthonormalise (Gram–Schmidt / QR) first, or use the full formula.

$UU^\top\neq U^\top U$. One is the (large) projection matrix onto the column space; the other is the identity when $U$ has orthonormal columns. Confusing the two is a frequent shape error.

The determinant is the signed factor by which $A$ scales volume: the unit cube maps to a parallelepiped of volume $|\det A|$, and the sign records whether orientation flipped. Everything else follows. $\det(AB)=\det A\det B$ because applying two maps multiplies their volume factors. $\det A=0$ exactly when the image collapses to lower dimension — no volume, no inverse. And since the eigenvalues are the stretch factors along the eigen-directions, their product is the total volume scaling: $\det A=\prod_i\lambda_i$. The trace is the matching sum, $\sum_i\lambda_i$, which is why it survives the change of basis in $\operatorname{tr}(AB)=\operatorname{tr}(BA)$.

For $A=\left(\begin{smallmatrix}2&1\\1&2\end{smallmatrix}\right)$: $\det A=2\cdot2-1\cdot1=3$ and $\operatorname{tr}A=4$. Cross-check against the eigenvalues $\lambda=3,1$ (from $\det(A-\lambda I)=0$): their product $3\cdot1=3=\det A$ and sum $3+1=4=\operatorname{tr}A$. The inverse is $\frac{1}{\det A}\left(\begin{smallmatrix}2&-1\\-1&2\end{smallmatrix}\right)$ — the $2\times2$ shortcut "swap the diagonal, negate the off-diagonal, divide by the determinant."

Order reverses under transpose and inverse. $(AB)^{-1}=B^{-1}A^{-1}$ and $(AB)^\top=B^\top A^\top$ — the factors flip. Writing $(AB)^{-1}=A^{-1}B^{-1}$ is one of the most common matrix-algebra errors; matrices don't commute, so the order is not optional.

$\det(A+B)\neq\det A+\det B$. The determinant is multiplicative, not additive. There is no useful formula for the determinant of a sum.

Never form $A^{-1}$ to solve $A\mathbf x=\mathbf b$. The textbook $\mathbf x=A^{-1}\mathbf b$ is correct mathematics but poor computation: explicitly inverting is slower, uses more memory, and amplifies rounding error. Use a solve instead — LU factorisation in general, or Cholesky when $A$ is symmetric positive-definite (twice as fast and more stable). In code this is solve(A, b), never inv(A) @ b.

The same applies to the determinant: computing it by cofactor expansion is exponential and numerically fragile; take it from the LU factorisation (product of the pivots) instead.

Stack the eigenvectors as the columns of $V$. By definition $AV=V\Lambda$ (each column is scaled by its own $\lambda$), and if the eigenvectors are independent $V$ is invertible, giving $A=V\Lambda V^{-1}$. The payoff is powers: $A^k=V\Lambda V^{-1}V\Lambda V^{-1}\cdots=V\Lambda^kV^{-1}$, because the inner $V^{-1}V$ pairs collapse. Raising a diagonal matrix to a power is just raising each $\lambda_i$ to that power — which is exactly why eigenvalues govern the long-term behaviour of iterations ($A^k$) and flows ($e^{At}$): the direction with the largest $|\lambda|$ eventually dominates.

For $A=\left(\begin{smallmatrix}2&1\\1&2\end{smallmatrix}\right)$, $\det(A-\lambda I)=(2-\lambda)^2-1=0$ gives $\lambda=3,1$. Solving $(A-3I)\mathbf v=0$ yields $\mathbf v_1=(1,1)$; $(A-I)\mathbf v=0$ yields $\mathbf v_2=(1,-1)$ — orthogonal, because $A$ is symmetric (the spectral theorem below). Then $A^{10}$, which would be tedious by repeated multiplication, is $V\operatorname{diag}(3^{10},1)V^{-1}$ — almost all of the action is the $3^{10}$ stretch along $(1,1)$.

Not every matrix diagonalizes. A defective matrix such as $\left(\begin{smallmatrix}1&1\\0&1\end{smallmatrix}\right)$ has a repeated eigenvalue but only one eigenvector direction, so no invertible $V$ exists and the factorisation fails. Repeated eigenvalues are the warning sign to check.

Don't assume real or orthogonal eigenstructure for a non-symmetric matrix. General real matrices can have complex eigenvalues (rotation — see Part 0½) and non-orthogonal eigenvectors. The clean "real eigenvalues, orthogonal eigenbasis" guarantee is special to symmetric matrices, not generic.

Take two eigenvectors with distinct eigenvalues, $A\mathbf u=\lambda\mathbf u$ and $A\mathbf w=\mu\mathbf w$. Compute $\mathbf u^\top A\mathbf w$ two ways: as $\mathbf u^\top(\mu\mathbf w)=\mu\,\mathbf u^\top\mathbf w$, and — using $A=A^\top$ — as $(A\mathbf u)^\top\mathbf w=\lambda\,\mathbf u^\top\mathbf w$. Subtract: $(\lambda-\mu)\,\mathbf u^\top\mathbf w=0$, and since $\lambda\neq\mu$ the eigenvectors must be orthogonal. Symmetry is exactly what forces the eigen-directions apart into a clean orthonormal frame; that is why $Q$ is orthogonal and $Q^{-1}=Q^\top$.

Test $A=\left(\begin{smallmatrix}2&1\\1&2\end{smallmatrix}\right)$. Its eigenvalues are $3$ and $1$, both positive, so $A\succ0$ — every $\mathbf x^\top A\mathbf x>0$. The quick $2\times2$ shortcut (Sylvester's criterion) agrees: leading minors $2>0$ and $\det=3>0$. Change the off-diagonal to $3$ and $\det=4-9<0$: now one eigenvalue is negative, the form is indefinite, and the corresponding surface is a saddle rather than a bowl.

Positive eigenvalues guarantee PD only when $A$ is symmetric. A non-symmetric matrix can have all-positive eigenvalues yet still give $\mathbf x^\top A\mathbf x<0$ for some $\mathbf x$. Always symmetrise (work with $\tfrac12(A+A^\top)$) before reading definiteness off the eigenvalues.

Semidefinite is not definite. A zero eigenvalue ($\lambda=0$) means $\mathbf x^\top A\mathbf x\ge0$ but can equal zero — the form has a flat direction. Covariance matrices of rank-deficient data are only positive semidefinite.

Jitter a covariance to keep it positive-definite. A covariance estimated from data is positive-definite in theory but, with rounding or near-collinear features, can come out with a tiny negative eigenvalue — which then breaks any Cholesky factorisation or Gaussian density evaluation. The standard fix is to add a small $\varepsilon I$ ($\varepsilon\sim10^{-6}$) to the diagonal, nudging every eigenvalue safely above zero.

Factor, don't invert. When you need $\Sigma^{-1}$ (e.g. a Gaussian exponent or a Mahalanobis distance), use the Cholesky factor $\Sigma=LL^\top$ and solve against $L$ — it is twice as fast as a general solve and numerically stable for symmetric positive-definite matrices.

Even when $A$ is not symmetric (or not even square), $A^\top A$ always is — and it is positive semidefinite. The spectral theorem hands it an orthonormal eigenbasis $V$ with non-negative eigenvalues; define $\sigma_i=\sqrt{\lambda_i}$. Setting $\mathbf u_i=A\mathbf v_i/\sigma_i$ produces an orthonormal output basis $U$, and assembling these gives $A=U\Sigma V^\top$. So SVD is the spectral theorem made universal: the symmetric guarantee, borrowed via $A^\top A$, applies to any matrix.

Suppose a $1000\times1000$ image has singular values that decay fast — say the first $20$ hold most of the magnitude. Eckart–Young says the rank-$20$ truncation $\sum_{i=1}^{20}\sigma_i\mathbf u_i\mathbf v_i^\top$ is the best possible rank-$20$ approximation in both Frobenius and spectral norm. Storing $20$ pairs $(\mathbf u_i,\mathbf v_i)$ plus $20$ scalars costs $\approx20\times2000$ numbers instead of $10^6$ — a 25× compression that keeps almost all the visible structure. This single fact is the engine behind PCA, image compression, and latent-factor models.

Singular values are not eigenvalues. They coincide only for symmetric positive-semidefinite matrices. In general $\sigma_i=\sqrt{\lambda_i(A^\top A)}$ — always real and non-negative even when $A$'s own eigenvalues are complex or negative.

$U$ and $V$ are different bases. For non-square $A$ they live in spaces of different dimension (output vs input); they are equal only in special symmetric cases. Don't conflate left and right singular vectors.

The condition number predicts trouble. $\kappa=\sigma_{\max}/\sigma_{\min}$ measures how badly a matrix amplifies error: solving $A\mathbf x=\mathbf b$ loses roughly $\log_{10}\kappa$ digits of accuracy. A $\kappa$ near $10^{16}$ means a double-precision solve retains essentially no correct digits — the matrix is numerically singular even if $\det A\neq0$ on paper.

This is why the determinant is a poor singularity test: a matrix can have a perfectly non-zero (even large) determinant while being catastrophically ill-conditioned. Look at the smallest singular value, not the determinant.

Write $X=U\Sigma V^\top$. Then $X^\top X=V\Sigma U^\top U\Sigma V^\top=V\Sigma^2 V^\top$, because $U^\top U=I$. That last expression is an eigendecomposition: $C=\tfrac1n X^\top X$ has eigenvectors $V$ (the right singular vectors of $X$) and eigenvalues $\sigma_i^2/n$. No separate covariance eigensolve is needed — the SVD of $X$ already contains it, which is also why computing PCA directly from $X$ is more numerically stable than forming $C$ first.

Centre first, always. PCA on un-centred data finds directions of maximum second moment, not variance — the first "component" just points at the mean. Subtract the column means before anything else.

Scaling changes the answer. Variance is unit-dependent, so features on larger numerical scales dominate the components. When units differ (e.g. kilograms vs millimetres) standardise to unit variance first, or PCA simply reports the biggest-numbered feature.

Don't compute variance as $E[X^2]-E[X]^2$. The textbook identity is exact in real arithmetic but suffers catastrophic cancellation in floating point: when the mean is large relative to the spread, the two big numbers $E[X^2]$ and $E[X]^2$ nearly cancel and the surviving digits are mostly rounding noise — you can even get a negative "variance." Use a two-pass formula (subtract the mean first, then sum squares) or Welford's online algorithm.

Prefer the SVD of $X$ over the eigendecomposition of $C$. Forming $C=X^\top X$ squares the condition number (it doubles the digit loss); running SVD on the centred $X$ directly avoids ever building $C$ and is the standard, stabler route to PCA.

The directional derivative is $D_{\mathbf u}f=\nabla f\cdot\mathbf u=\lVert\nabla f\rVert\cos\theta$, where $\theta$ is the angle between $\mathbf u$ and the gradient. This is largest when $\cos\theta=1$ — i.e. when $\mathbf u$ points along $\nabla f$ — so the gradient direction is steepest ascent, and its length is the steepest slope. Along a level set the function does not change, so $D_{\mathbf u}f=0$ there, forcing $\nabla f\perp\mathbf u$: the gradient is perpendicular to every contour. Both headline facts fall straight out of one dot product.

For $f(x,y)=x^2+3y^2$, $\nabla f=(2x,\,6y)$. At the point $(1,1)$ it is $(2,6)$ — steepest ascent points more in $y$ than $x$ because the surface climbs three times as fast in $y$. The contour through $(1,1)$ has tangent direction $(6,-2)$ (perpendicular to the gradient, as it must be). Step downhill by $-\eta\nabla f$ and you have one iteration of gradient descent on this bowl.

The gradient points uphill, not downhill. Minimisation moves along $-\nabla f$. Dropping the minus sign sends gradient descent climbing the loss — a classic sign bug.

The directional derivative needs a unit vector. $D_{\mathbf u}f=\nabla f\cdot\mathbf u$ only measures slope-per-unit-distance when $\lVert\mathbf u\rVert=1$; feeding an un-normalised direction silently scales the answer.

Each row of $J$ is the gradient of one output component, so $J\mathbf h$ collects all the first-order changes $\nabla f_i\cdot\mathbf h$ into a vector — exactly the multivariate version of $f(x+h)\approx f(x)+f'(x)h$. The Jacobian is therefore not a notational convenience but the derivative: the unique linear map that best approximates $\mathbf f$ near $\mathbf x$. The Hessian plays the same role for the gradient: it is the Jacobian of $\nabla f$, so it records how the slope itself changes — curvature — which is why its sign-definiteness sorts minima from maxima from saddles.

For $\mathbf f(x,y)=(x^2y,\ \sin x+y)$, differentiate each output by each input:

$$J=\begin{pmatrix}2xy & x^2\\ \cos x & 1\end{pmatrix}.$$

At $(0,1)$ this is $\left(\begin{smallmatrix}0&0\\1&1\end{smallmatrix}\right)$, whose determinant is $0$ — the map locally collapses a direction, so it is not invertible there. The Jacobian determinant is precisely the local volume-scaling factor that appears when you change variables in a multiple integral.

Jacobian shape and ordering. It is $m\times n$ — outputs index rows, inputs index columns. Transposing it (the numerator/denominator-layout confusion of the matrix-calculus section below) silently breaks any chain-rule product it feeds into.

The Hessian test needs a critical point. $H\succ0$ certifies a minimum only where $\nabla f=\mathbf 0$ already. Reading curvature at a non-stationary point tells you nothing about min vs max.

Symmetry assumes continuous second partials. $H$ is symmetric by Clairaut's theorem when the mixed partials are continuous; for pathological functions $\partial_x\partial_y\neq\partial_y\partial_x$ and the clean eigenvalue test no longer applies.

Near a point each map is its Jacobian: $\mathbf f(\mathbf x+\mathbf h)\approx\mathbf f(\mathbf x)+J_f\mathbf h$, then $\mathbf g$ acts on that with $J_g$. Composing, a small input nudge $\mathbf h$ produces output nudge $J_g J_f\mathbf h$ — so the Jacobian of the composition is the product $J_g J_f$. Backprop is this rule applied across many layers: the chain $J_L\cdots J_2 J_1$ is evaluated right-to-left because contracting a row-vector gradient against one matrix at a time is far cheaper than forming the full product. Each layer multiplies the incoming gradient by its own local Jacobian and passes it on.

Take $\mathbf y=W_1\mathbf x$, $\mathbf z=\sigma(\mathbf y)$ (elementwise), scalar loss $L=\tfrac12\lVert\mathbf z-\mathbf t\rVert^2$. The backward pass multiplies local Jacobians right-to-left:

$$\frac{\partial L}{\partial\mathbf z}=(\mathbf z-\mathbf t),\quad \frac{\partial L}{\partial\mathbf y}=(\mathbf z-\mathbf t)\odot\sigma'(\mathbf y),\quad \frac{\partial L}{\partial W_1}=\frac{\partial L}{\partial\mathbf y}\,\mathbf x^\top.$$

The $\odot\sigma'(\mathbf y)$ step is the elementwise Jacobian of the activation; forgetting it is the single most common backprop bug. Notice the gradient w.r.t. a weight matrix is an outer product of the incoming gradient with the layer input — the shape every framework computes.

Don't drop the inner Jacobian. Just as in 1D, every layer contributes its own factor; omitting an activation's $\sigma'$ or a layer's weight matrix is a silently-wrong gradient that still runs.

Order is fixed — matrices don't commute. $J_g J_f\neq J_f J_g$. The backward pass must multiply in reverse layer order; swapping factors gives a shape mismatch at best and a wrong-but-runnable gradient at worst.

Vanishing and exploding gradients are a conditioning problem. The backprop chain $J_L\cdots J_2 J_1$ is a product of many matrices; its overall gain is governed by the product of their singular values. If typical singular values sit below 1 the gradient shrinks geometrically with depth (vanishing); above 1 it blows up (exploding). This is the same repeated-multiplication instability that the condition number flags in Part II — depth simply compounds it.

The standard mitigations all target that product: careful weight initialisation (keep per-layer gain near 1), normalisation layers, residual/skip connections (which add an identity term so the product stays near 1), and gradient clipping for the exploding case.

Expand $f(\mathbf x)=\mathbf x^\top A\mathbf x=\sum_{i,j}x_i A_{ij} x_j$ and differentiate by one component $x_k$. The variable $x_k$ appears in two places — as the $i=k$ term and the $j=k$ term — so the product rule gives $\partial f/\partial x_k=\sum_j A_{kj}x_j+\sum_i x_i A_{ik}$, which is the $k$-th entry of $(A+A^\top)\mathbf x$. When $A$ is symmetric the two sums coincide and you get $2A\mathbf x$. The least-squares line follows by setting $\mathbf r=A\mathbf x-\mathbf b$ and chaining: $\nabla\lVert\mathbf r\rVert^2=2A^\top\mathbf r$.

Minimise $\lVert A\mathbf x-\mathbf b\rVert_2^2$. Set its gradient to zero using the last table row:

$$2A^\top(A\mathbf x-\mathbf b)=\mathbf 0\ \Longrightarrow\ A^\top A\,\mathbf x=A^\top\mathbf b\ \Longrightarrow\ \hat{\mathbf x}=(A^\top A)^{-1}A^\top\mathbf b.$$

This is the normal equation: the residual $A\hat{\mathbf x}-\mathbf b$ is orthogonal to the column space of $A$ (every column $A^\top$ kills it), which is exactly the projection condition from Part II's inner-product card.

Pick a layout and stay in it. Mixing denominator-layout gradients (columns) with numerator-layout ones (rows) inside a single derivation introduces stray transposes that don't cancel. Most "off by a transpose" bugs are a layout switch mid-calculation.

$\nabla(\mathbf x^\top A\mathbf x)=2A\mathbf x$ only when $A$ is symmetric. For general $A$ it is $(A+A^\top)\mathbf x$. Assuming symmetry on a non-symmetric quadratic form drops half the gradient.

Don't solve least squares by forming $(A^\top A)^{-1}$. The normal equation is correct mathematics but poor numerics: building $A^\top A$ squares the condition number, so you lose twice as many digits, and inverting compounds it. Solve the least-squares problem directly with a QR factorisation of $A$ (or SVD when $A$ is rank-deficient or ill-conditioned) — both work on $A$ itself and never form the squared, fragile $A^\top A$.

This is the same lesson as Part II's "never form $A^{-1}$," sharpened: here even assembling the matrix you would invert is the mistake.

Bayes is just the product rule read both ways. $P(A\cap B)$ can be factored as $P(A\mid B)P(B)$ or as $P(B\mid A)P(A)$; equate the two and divide by $P(B)$. The marginal rule (total probability) supplies the denominator $P(B)=\sum_a P(B\mid a)P(a)$, which is why the posterior automatically integrates to 1 — the normaliser is the sum of the numerator over every hypothesis.

$P(A\mid B)\neq P(B\mid A)$. The two are related by the base rates $P(A)$ and $P(B)$, not equal. This is the base-rate fallacy: a test that is 99% sensitive for a disease with 0.1% prevalence still yields mostly false positives, because the rare prior overwhelms the likelihood. Sensitivity is $P(\text{positive}\mid\text{disease})$; what you want is $P(\text{disease}\mid\text{positive})$ — Bayes converts one to the other.

Shifting by $b$ moves the mean by $b$ too, so the deviations $X-\mathbb EX$ are unchanged — additive constants never affect spread. Scaling by $a$ scales each deviation by $a$, and variance is built from squared deviations, so it scales by $a^2$. The sign of $a$ drops out, which is why variance is always non-negative and why standard deviation (its square root) carries the units of $X$.

Uncorrelated is not independent. $\operatorname{Cov}(X,Y)=0$ only kills linear dependence; $Y=X^2$ with symmetric $X$ has zero covariance yet is fully determined by $X$. The one exception: jointly Gaussian variables that are uncorrelated really are independent.

$\mathbb E[g(X)]\neq g(\mathbb E[X])$ for nonlinear $g$ (Jensen). For convex $g$ the expectation sits above the function of the mean. Plugging the mean into a nonlinearity — a tempting shortcut — systematically biases the answer.

Don't compute variance from $\mathbb E[X^2]-(\mathbb E X)^2$ in floating point. The identity is exact in real arithmetic but suffers catastrophic cancellation when the mean is large relative to the spread: two big nearly-equal numbers subtract to leave mostly rounding noise, and you can even get a negative "variance." Use a two-pass formula (subtract the mean first) or Welford's online algorithm. This is the same cancellation discussed in card-svd-pca.

The law of large numbers says a sample mean converges to the true mean as the sample grows — the formal licence for using averages at all. The central limit theorem says the fluctuations around that mean become Gaussian whatever the underlying distribution, which is the deeper reason Gaussians turn up everywhere. Together they underwrite most of what Dayan & Abbott assume when averaging spike counts or fitting distributions. ▶ 3Blue1Brown — Central Limit Theorem.

A density is not a probability. A pdf $p(x)$ can exceed 1 (a tight Gaussian peaks above 1); only its integral over an interval is a probability, bounded by 1. Reading a pdf value as a probability is a common slip.

Softmax / log-sum-exp overflow. Turning scores into a distribution via $e^{x_i}/\sum_j e^{x_j}$ overflows the moment any $x_i$ is large. Subtract the maximum first: $e^{x_i-m}/\sum_j e^{x_j-m}$ with $m=\max_j x_j$ is mathematically identical but keeps every exponent $\le 0$. The log-normaliser is then $m+\log\sum_j e^{x_j-m}$ (the log-sum-exp trick).

Sigmoid overflow. $1/(1+e^{-x})$ overflows for very negative $x$ (the $e^{-x}$ blows up). Compute it piecewise — use $e^{x}/(1+e^{x})$ when $x<0$ — so the exponent is always non-positive.

Marginals, conditionals, and linear transforms of Gaussians stay Gaussian; it is the maximum-entropy distribution for fixed mean and covariance; and it is self-conjugate, so Bayesian updates are closed-form. The quadratic in the exponent is why "completing the square" (Part 0) is the key skill.

Never form $\Sigma^{-1}$ explicitly. Evaluating the exponent needs $(\mathbf x-\boldsymbol\mu)^\top\Sigma^{-1}(\mathbf x-\boldsymbol\mu)$, but inverting $\Sigma$ squares its condition number. Cholesky-factor $\Sigma=LL^\top$ and solve $L\mathbf y=\mathbf x-\boldsymbol\mu$ instead; then the quadratic form is just $\|\mathbf y\|^2$ and $\log|\Sigma|=2\sum_i\log L_{ii}$. Add a small $\varepsilon I$ jitter to keep a near-singular covariance positive definite.

MAP maximises the posterior $p(\theta\mid D)\propto p(D\mid\theta)\,p(\theta)$. Taking logs turns the product into the sum $\ln p(D\mid\theta)+\ln p(\theta)$, and the evidence $p(D)$ drops out because it doesn't depend on $\theta$. MLE is the special case of a flat prior, where the second term is constant — so the only difference between the two is whether you carry prior knowledge into the objective.

Maximise the log-likelihood, never the raw likelihood. A likelihood is a product of many sub-1 probabilities; over a real dataset that product underflows to exactly 0 and all information is lost. Working in log-space turns the product into a sum $\sum_i\ln p(x_i\mid\theta)$ — which is also why the MLE objective is written as a sum in the first place.

Entropy is average surprise — bits to describe $X$, maximal when uniform. Mutual information is how many bits knowing $Y$ saves you about $X$ (zero iff independent) — the core measure of how much a neural response tells you about a stimulus. KL is a directed "distance" between distributions; minimising it is what variational inference and cross-entropy loss do.

KL is not symmetric and not a metric. $D_{\mathrm{KL}}(p\Vert q)\neq D_{\mathrm{KL}}(q\Vert p)$ in general, and it obeys no triangle inequality — calling it a "distance" is a convenient lie. The direction you minimise matters: $D_{\mathrm{KL}}(p\Vert q)$ (forward) is mean-seeking, $D_{\mathrm{KL}}(q\Vert p)$ (reverse) is mode-seeking.

Guard the $\log 0$ in every sum. Entropy, KL, and cross-entropy all contain $p\log p$ or $p\log(p/q)$ terms that go to $0\cdot(-\infty)$ when a probability is zero. Use the convention $0\log 0=0$ explicitly, and where probabilities come from a softmax, work from the log-probabilities (log-sum-exp, see #dists) rather than re-logging exponentiated values.

Suppose $\mathbf x^\star$ is a local min but some $\mathbf y$ has $f(\mathbf y)

"Local min is global" needs convexity. The guarantee is a property of convex functions, not a general fact. Deep-learning losses are non-convex, so the comforting theorem simply does not apply — there can be many minima of different quality.

A first-order Taylor step gives $f(\mathbf x-\eta\nabla f)\approx f(\mathbf x)-\eta\|\nabla f\|^2$. The correction is negative whenever $\nabla f\neq 0$, so for small enough $\eta$ each step lowers $f$ — the gradient is the direction of steepest ascent, so its negative is the steepest descent. The $L$-smoothness bound is exactly how large $\eta$ can be before the quadratic remainder overwhelms that gain.

The gradient points uphill. $\nabla f$ is the direction of steepest increase; you descend along $-\nabla f$. Dropping the minus sign climbs the loss.

$\nabla f=0$ is necessary, not sufficient. It also holds at maxima and saddles. Only the Hessian distinguishes them — and in high dimensions saddles vastly outnumber minima.

The learning rate sets the regime. Below $\eta=1/L$ the iterates converge monotonically; near $2/L$ they oscillate across the valley without diverging; above $2/L$ they diverge outright. Tuning $\eta$ is therefore a stability question, not just a speed one — too large doesn't slow learning, it breaks it.

Sigmoid/softmax inside the loss. When $f$ contains a sigmoid or softmax, evaluate them in the overflow-safe form (see Part IV #dists) before differentiating, or the gradients themselves return NaN.

Approximate $f$ near $\mathbf x_t$ by its second-order Taylor model $f(\mathbf x_t)+\nabla f^\top\mathbf p+\tfrac12\mathbf p^\top H\mathbf p$. Minimising this quadratic in the step $\mathbf p$ sets its gradient to zero: $\nabla f+H\mathbf p=0$, so $\mathbf p=-H^{-1}\nabla f$. Newton jumps exactly to the minimum of the local quadratic model — which is why it converges in one step on a true quadratic and quadratically fast near any smooth minimum.

Never invert the Hessian — solve the Newton system. Compute the step by solving $H\mathbf p=-\nabla f$ (Cholesky when $H\succ0$) rather than forming $H^{-1}$ and multiplying. It is cheaper, stabler, and avoids the squared condition number of an explicit inverse. Away from a minimum $H$ may be indefinite, in which case the raw step can point uphill — damp it or add $\tau I$ to restore positive definiteness.

On the constraint surface $g=0$, any feasible move is tangent to it — perpendicular to $\nabla g$. At an optimum no tangent move can decrease $f$, so $\nabla f$ has no component along the surface; it must lie entirely along $\nabla g$, i.e. $\nabla f=\lambda\nabla g$. KKT generalises this: an inactive inequality ($g_i<0$) leaves $f$ free, so its multiplier is forced to zero ($\mu_i g_i=0$), and only active constraints push back.

Complementary slackness is the whole game. Either a constraint is active ($g_i=0$) or its multiplier vanishes ($\mu_i=0$). Forgetting this and carrying nonzero multipliers on inactive constraints is the most common KKT error.

The characteristic roots of a linear ODE are exactly the eigenvalues of its companion matrix. This is no coincidence: solving an ODE and diagonalising a matrix are the same operation in different clothes. The same algebra — characteristic polynomial, real roots → exponentials, complex roots → oscillations — runs through Part II and Part VII because it is the same problem.

It's separable: $\frac{dV}{V-V_\infty}=-\frac{dt}{\tau}$. Integrating both sides gives $\ln|V-V_\infty|=-t/\tau+C$, so $V-V_\infty=(V_0-V_\infty)e^{-t/\tau}$. The right side of the original equation is "$-$(distance from rest)", so the system always moves toward $V_\infty$, fast when far, slow when near — exactly exponential approach.

Whenever you see $\tau\dot y=-(y-\text{target})$, the answer is "exponential relaxation to target, timescale $\tau$" — no integration needed. After one time constant $t=\tau$ the gap has shrunk to $e^{-1}\approx37\%$ of its start; after $3\tau$ to about $5\%$; after $5\tau$ it is essentially settled. You will see this shape dozens of times.

Insert $A=V\Lambda V^{-1}$ into the series: every term $\frac{(At)^k}{k!}=V\frac{(\Lambda t)^k}{k!}V^{-1}$, and the inner $V^{-1}V$ pairs collapse on summation, leaving $e^{At}=Ve^{\Lambda t}V^{-1}$. In eigen-coordinates $\mathbf z=V^{-1}\mathbf x$ the system becomes $\dot z_i=\lambda_i z_i$ — $n$ uncoupled scalar leaky-integrator equations, each solved by $e^{\lambda_i t}$. Diagonalisation literally decouples the modes.

Continuous-time stability is $\operatorname{Re}(\lambda)<0$. For a flow, modes decay when the real part is negative — not when $|\lambda|<1$. That magnitude test is for discrete maps (Part VII), and mixing the two criteria is the classic error when the same matrix is read as a flow versus an iteration.

Forward Euler is just the first-order Taylor step (Part I): it follows the tangent for $\Delta t$ and re-reads the slope, so its local error is $O(\Delta t^2)$ and its global error $O(\Delta t)$. Runge–Kutta samples the slope at several intermediate points to cancel low-order error, buying $O(\Delta t^4)$ accuracy per step. But accuracy isn't the only risk: for a stiff system (a fast mode with large $|\lambda|$) Euler is unstable unless $\Delta t<2/|\lambda|$ — too big a step doesn't just lose precision, it diverges. This is the discrete-map stability condition ($|\lambda|<1$ for the update) showing up in simulation.

Linearise: near $x^*$ the perturbation obeys $\dot\eta=f'(x^*)\eta$, whose solution is $\eta(t)=\eta_0 e^{f'(x^*)t}$ — the leaky integrator again. A negative slope makes the exponent negative, so perturbations decay back to $x^*$ (stable); a positive slope makes them grow away (unstable). At $f'(x^*)=0$ the linear term vanishes and the nonlinear terms decide.

Real part of eigenvalues = growth/decay; nonzero imaginary part = rotation. The discriminant $\tau^2-4\Delta$ decides real-vs-complex (node vs spiral), the sign of $\tau$ decides growing-vs-decaying, and the sign of $\Delta$ flags the saddle. One table classifies every linear planar flow. For a beginner: the stable spiral and the saddle are the two cases that appear most in neuroscience models — internalise those first; the others exist and matter, but these are the ones to know cold.

Linearization fails at non-hyperbolic fixed points. When an eigenvalue has zero real part — a linear center being the classic case — the linearization is silent on stability, and the nonlinear terms can turn that center into a slow inward or outward spiral. The clean "read it off the Jacobian" rule only holds away from this borderline; a center on the table is a warning, not an answer.

In FitzHugh–Nagumo $\dot v=v-\tfrac{v^3}{3}-w+I,\ \dot w=\varepsilon(v+a-bw)$, the $v$-nullcline is a cubic "N" and the $w$-nullcline is a line. How they intersect — and how injecting current $I$ slides the cubic up — is the explanation of spiking, thresholds, and refractoriness. Learn to read this picture and excitability stops being mysterious.

A nullcline is not a set of fixed points. It is where one derivative vanishes — the system still moves along it (the other variable keeps changing). Only the intersection of both nullclines, where both derivatives are zero, is a fixed point. Treating a whole nullcline as "where the system rests" is a frequent slip.

A complex-conjugate eigenvalue pair crosses the imaginary axis ($\operatorname{Re}\lambda$ changes sign): a stable spiral becomes a small oscillation. The onset of repetitive firing in many neurons is a Hopf bifurcation — the resting state loses stability and a limit cycle is born around it.

At a SNIC the limit cycle is born with infinite period (the orbit must crawl through the bottleneck where the saddle-node just vanished), so firing starts at frequency zero and ramps up continuously — a Type I curve. At a Hopf the cycle appears at finite frequency set by the imaginary part of the eigenvalues, so firing switches on abruptly at a nonzero rate — a Type II curve. The geometry of how the cycle is created is exactly what you measure as the $f$–$I$ shape.

Map stability is $|\lambda|<1$, not $\operatorname{Re}(\lambda)<0$. A discrete map's fixed point is stable when the multiplier has magnitude below 1; the real-part test belongs to continuous flows (Part VI). The two criteria are easy to swap, and the logistic map's $|f'|<1$ condition is the map version. Sensitive dependence also doesn't require randomness — chaos is fully deterministic; it only requires a positive Lyapunov exponent.

The Dirac delta $\delta(t)$ is an idealised instantaneous pulse with unit area, defined by the sifting property $\int f(t)\,\delta(t-t_0)\,dt=f(t_0)$. Any input can be written as a continuous superposition of shifted, scaled deltas — that's the sifting identity read as a sum. Linearity and time-invariance then say the output is the same superposition of shifted, scaled impulse responses, which is precisely the convolution integral $h*\text{input}$. So measuring one pulse response determines the system's reply to every possible input. (Dayan & Abbott use $\delta$ from chapter one without defining it.)

Convolution in time = multiplication in frequency: $\widehat{f*g}=\hat f\cdot\hat g$. So to filter a signal you can transform both signal and kernel, multiply pointwise, and transform back — a messy sliding integral becomes a single product. This is why frequency-domain thinking simplifies filtering, and why fast convolution is done via the FFT.

The pipeline reads stage by stage: stimulus $s$ → linear filter $h*s$ (the linear receptive field) → nonlinearity $N$ → firing rate $r$. To recover $h$, drive the cell with white noise (a stimulus whose values are uncorrelated across time) and average the stimulus snippets just before each spike. Because white noise has no internal correlations, the only structure the average can pick up is what the filter $h$ preferentially passed — so the spike-triggered average converges to a scaled copy of $h$. The monotone nonlinearity $N$ only rescales spike probability; it doesn't bend the recovered shape. This reverse-correlation logic is what makes the linear stage measurable from data.