In binary, \(0.1 = 0.0001100110011\ldots_2\) is an infinite repeating fraction — like \(1/3\) in decimal. The computer stores the nearest 64-bit value, which is not exactly \(0.1\).
So 0.1 + 0.2 returns 0.30000000000000004, and 0.1 + 0.2 == 0.3 is FALSE in R, Python, and every IEEE 754 language.
Rule: never compare floats with ==. Use a tolerance: abs(a - b) < 1e-8, all.equal() in R, np.isclose() in Python.
📝 Precision Formats at a Glance
Format
Bits
Decimal digits
Machine ε
Default in
float64 (double)
64
~15.9
2.2×10−16
R, NumPy, Stata
float32 (single)
32
~7.2
1.2×10−7
torch, GPUs, deep learning
float16 / bfloat16
16
~3.3
9.8×10−4
GPU tensor cores
Consumer GPUs run float32 30–60× faster than float64 — that is why torch defaults to float32. For econometric estimates that must match R to machine precision, request dtype=torch.float64 and accept the speed penalty (Part VIII returns to this trade-off).
IEEE 754 — How Computers Store Numbers
Every double-precision number is stored as sign, exponent, and a 52-bit fraction:
You never need the built-in constant — you can measure\(\varepsilon_{\text{mach}}\) directly. Start the increment at \(1\) and halve it repeatedly; the machine epsilon is the last increment the computer can still tell apart from \(1\):
Under round-to-nearest, \(1 + 2^{-53}\) ties back to \(1\), so the halving stops one step early — after exactly \(52\) halvings, at \(2^{-52}\). Repeat in single precision and the loop stops far sooner (\(2^{-23}\)): fewer mantissa bits, coarser spacing.
Code
# Shrink the increment until 1 + eps/2 is no longer distinguishable from 1eps <-1steps <-0while (1+ eps/2>1) { eps <- eps/2 steps <- steps +1}cat("computed eps :", eps, " after", steps, "halvings\n")
# float32 carries only ~7 digits: the same loop stops much soonere32, k = np.float32(1.0), 0while np.float32(1.0) + e32/np.float32(2.0) > np.float32(1.0): e32 /= np.float32(2.0) k +=1print(f"\nfloat32 eps : {e32:.8e} after {k} halvings (2**-23)")
float32 eps : 1.19209290e-07 after 23 halvings (2**-23)
Subtracting two nearly equal numbers destroys the leading (correct) digits and promotes rounding noise to the front:
\[x = 1.234567\underbrace{89012345}_{\text{noise}}, \quad y = 1.234567\underbrace{88999999}_{\text{noise}}, \quad x - y = 0.00000000012346\ldots\]
The relative error of the difference explodes even though \(x\) and \(y\) were each accurate to 16 digits — the leading digits that agreed simply vanish, and rounding noise is promoted to the front.
When the mean is large relative to the spread, \(\sum x_i^2\) and \(n\bar{x}^2\) are nearly equal — the subtraction cancels catastrophically and can return a negative variance.
The stable two-pass formula subtracts the mean first, while numbers are still small:
Same algebra, completely different arithmetic. Think of stock prices near 10,000 index points with daily changes of a few points — exactly this regime.
The likelihood of even a modest sample underflows: \(1000\) observations with density values around \(0.3\) give \(L = 0.3^{1000} \approx 10^{-523} \to 0\). All likelihood work happens in logs — but logs alone are not enough when you must sum exponentials (mixtures, logit denominators, particle filters):
\[\log \sum_{j=1}^{J} e^{a_j} = M + \log \sum_{j=1}^{J} e^{a_j - M}, \qquad M = \max_j a_j\]
Shifting by the max makes the largest exponent exactly \(e^0 = 1\): no overflow, no total underflow.
Code
# Log-likelihood contributions of a mixture component: large negative numbersa <-c(-1000, -1001, -1002) # log f_j(x_i): perfectly ordinary in MLE work# Naive: exponentiate first -> everything underflows to 0 -> log(0) = -Infnaive <-log(sum(exp(a)))cat("naive log(sum(exp(a))):", naive, "\n")
naive log(sum(exp(a))): -Inf
Code
# Log-sum-exp: shift by the max before exponentiatinglogsumexp <-function(a) { M <-max(a) M +log(sum(exp(a - M)))}cat("stable logsumexp(a) :", logsumexp(a), "\n")
stable logsumexp(a) : -999.5924
Code
# Same trick under the hood of a logit log-likelihood: use log1p / plogis# log(1 + exp(eta)) overflows for eta > ~710; the stable form never doeseta <-800cat("naive log(1+exp(eta)) :", log(1+exp(eta)), "\n")
naive log(1+exp(eta)) : Inf
Code
cat("stable via max trick :", max(eta, 0) +log1p(exp(-abs(eta))), "\n")
stable via max trick : 800
Code
import numpy as npfrom scipy.special import logsumexp # library version of the tricka = np.array([-1000.0, -1001.0, -1002.0]) # log-density values# Naive: exponentiate first -> underflow -> log(0) = -infprint("naive log(sum(exp(a))):", np.log(np.sum(np.exp(a))))
naive log(sum(exp(a))): -inf
Code
# Stable: shift by the max before exponentiatingM = a.max()print("manual logsumexp :", M + np.log(np.sum(np.exp(a - M))))
manual logsumexp : -999.5923940355556
Code
print("scipy.special.logsumexp:", logsumexp(a))
scipy.special.logsumexp: -999.5923940355556
Code
# Logit likelihood term log(1 + exp(eta)): naive overflows for eta > ~710eta =800.0print("naive log(1+exp(eta)) :", np.log(1+ np.exp(eta)))
\[\text{correct digits in } \hat{\mathbf{x}} \;\approx\; 16 - \log_{10} \kappa(\mathbf{A})\]
\(\kappa = 10^2\): ~14 digits — harmless
\(\kappa = 10^8\): ~8 digits — noticeable
\(\kappa = 10^{15}\): zero correct digits — the near-collinear regression case
A stable algorithm on an ill-conditioned problem still gives a bad answer — but the best possible bad answer. An unstable algorithm (normal equations!) makes it needlessly worse by squaring\(\kappa\).
Request float64 when precision matters (Part VIII)
Code
# Two quick demonstrations from the checklist# 1. Summation is not associative: same numbers, different order, different sumset.seed(14159)v <-rnorm(1e6) *1e8cat("sum ascending :", format(sum(sort(v)), digits =17), "\n")
Nearly every estimator reduces to solving a linear system:
\[\hat{\boldsymbol\beta}^{OLS} : \; (\mathbf{X}^\top\mathbf{X})\,\hat{\boldsymbol\beta} = \mathbf{X}^\top\mathbf{y}
\qquad \text{GLS, IV, GMM, Newton steps, Kalman filters — all the same shape } \mathbf{A}\mathbf{x} = \mathbf{b}\]
The textbook writes \(\hat{\boldsymbol\beta} = (\mathbf{X}^\top\mathbf{X})^{-1}\mathbf{X}^\top\mathbf{y}\) — but the textbook formula is notation, not an algorithm.
Computing \(\mathbf{A}^{-1}\) explicitly costs ~3× the flops of solving the system, and is less accurate
Forming \(\mathbf{X}^\top\mathbf{X}\)squares the condition number: \(\kappa(\mathbf{X}^\top\mathbf{X}) = \kappa(\mathbf{X})^2\)
Good software (R’s lm(), statsmodels) never inverts — it factorises: QR, Cholesky, or SVD
Rule
Instead of
Do
Never invert to solve
solve(A) %*% b
solve(A, b)
Never form X’X if avoidable
solve(t(X) %*% X, t(X) %*% y)
qr.solve(X, y) / lstsq
Exploit structure
generic solve() on SPD matrix
Cholesky: chol() / cho_solve
Check conditioning
trust any output
kappa(X), np.linalg.cond(X)
Exploit sparsity
dense dummies for 10,000 FE
Matrix::sparseMatrix, scipy.sparse
📝 Inside lm() and statsmodels.OLS
R's lm() never computes \((\mathbf{X}^\top\mathbf{X})^{-1}\). It calls LINPACK/LAPACK's QR decomposition with column pivoting: \(\mathbf{X} = \mathbf{Q}\mathbf{R}\), then solves the triangular system \(\mathbf{R}\boldsymbol\beta = \mathbf{Q}^\top\mathbf{y}\) by back-substitution.
Pivoting is why lm() returns NA coefficients (rather than garbage) for perfectly collinear columns: the rank deficiency is detected in \(\mathbf{R}\)'s diagonal.
statsmodels.OLS defaults to the Moore–Penrose pseudoinverse via SVD — even more robust, slightly slower. Both work with \(\kappa(\mathbf{X})\), not \(\kappa(\mathbf{X})^2\).
The Four Decompositions
Factor once, then solve cheaply by forward/back-substitution (\(O(n^2)\) per right-hand side):
LU with partial pivoting — general square systems:
set.seed(14159)n <-1000A <-crossprod(matrix(rnorm(n * n), n, n)) + n *diag(n) # SPD, well-conditionedb <-rnorm(n)# Three ways to get x with A x = bmb <-microbenchmark(invert =solve(A) %*% b, # explicit inverse: never do thissolve =solve(A, b), # LU solvecholesky =backsolve(chol(A), forwardsolve(t(chol(A)), b)),times =10)print(summary(mb)[, c("expr", "median")], digits =3)
A panel with \(10{,}000\) individual dummies gives a design matrix that is \(>99.9\%\) zeros. Dense storage wastes memory and flops; sparse storage keeps only the non-zeros — this is how fixest/reghdfe/pyfixest scale.
\[\text{dense: } n \times p \times 8 \text{ bytes} \qquad \text{sparse (CSC): } \approx \text{nnz} \times 12 \text{ bytes}\]
Code
set.seed(14159)n_id <-1000; t_per <-5; n <- n_id * t_per # 5,000 obs, 1,000 FE dummiesid <-factor(rep(seq_len(n_id), each = t_per))x <-rnorm(n)y <-1.5* x +rep(rnorm(n_id), each = t_per) +rnorm(n)# Dense vs sparse dummy design matrixX_dense <-model.matrix(~ x + id) # dense: hugeX_sparse <- Matrix::sparse.model.matrix(~ x + id) # sparse: tinycat(sprintf("dense size: %8.1f MB\n", as.numeric(object.size(X_dense)) /1e6))
Each iteration solves a linear system (Part II: factorise, never invert). Quasi-Newton (Broyden) replaces \(\mathbf{J}\) with a cheap rank-one update when derivatives are expensive.
Economic application — equilibrium in two interdependent markets. Demand and supply for goods 1 and 2 (substitutes, nonlinear in prices):
The optimisers and ODE solvers built here are tools without an application. The companion deck Numerical Dynamic Programming and Heterogeneous Agents is where they get used on a real economic problem — value-function iteration, policy functions, and the stationary distributions of heterogeneous-agent models.
Almost every estimator beyond OLS is defined as an optimum:
When optim() or scipy.minimize reports “convergence”, it has satisfied a numerical criterion — not found the truth. Understanding the algorithms tells you when to trust it.
Method
Uses
Cost per step
Convergence
Library
Gradient descent
\(\nabla Q\)
cheap
linear, rate depends on step size
(by hand, torch)
Newton–Raphson
\(\nabla Q\), \(\mathbf{H}\)
\(O(p^3)\) solve
quadratic
glm() (Fisher scoring)
BFGS (quasi-Newton)
\(\nabla Q\) only
\(O(p^2)\)
superlinear
optim(method="BFGS"), scipy
Nelder–Mead
\(Q\) only
cheap
slow, no rate
optim() default, scipy
L-BFGS-B
\(\nabla Q\), bounds
\(O(mp)\)
superlinear
both — large \(p\) default
Rule: derivatives available → BFGS/Newton. Noisy or non-smooth objective → Nelder–Mead. Thousands of parameters → L-BFGS or gradient methods (Part VIII).
⚠️ What "converged" Actually Certifies
An optimizer stops when its own criteria hold: small gradient norm, small step, or small objective change. This certifies a stationary point of the sample objective — nothing more. It does not certify:
a global optimum (multi-modal likelihoods: mixtures, GARCH-in-mean, random coefficients)
a well-identified parameter (flat likelihood → near-singular Hessian → huge SEs, Part II conditioning)
correctness of your gradient code — always check analytic vs numerical gradients before trusting results
Routine practice: restart from several starting points; check the Hessian is positive definite at the optimum; compare analytic and numDeriv::grad() gradients at a random point.
Optimization — The Mathematics
Gradient descent — step against the gradient with learning rate \(\gamma\):
Converges linearly for convex \(Q\) if \(\gamma < 2/L\) (\(L\) = largest Hessian eigenvalue). Too large \(\gamma\)diverges; too small crawls. The convergence rate is governed by the Hessian’s condition number \(\kappa\) (Part II again):
with \(\mathbf{W} = \text{diag}\{\mu_i(1-\mu_i)\}\). Simulated data (\(n = 5000\), true \(\boldsymbol\beta = (-0.5, 1.2, -0.8, 0.6)\)) shared by both languages via ../data/numerical-logit.csv.
Code
dat <-read.csv("../data/numerical-logit.csv")X <-cbind(1, as.matrix(dat[, c("x1", "x2", "x3")]))y <- dat$y# Stable log(1 + exp(eta)) — avoids overflow for large |eta| (Part I)log1pexp <-function(eta) pmax(eta, 0) +log1p(exp(-abs(eta)))loglik <-function(b) sum(y * (X %*% b) -log1pexp(X %*% b))# Newton-Raphson: each step solves H d = -grad (Part II: solve, never invert)b <-rep(0, 4)cat("Newton-Raphson iterations:\n")for (k in1:25) { mu <-plogis(as.numeric(X %*% b)) grad <-crossprod(X, y - mu) H <--crossprod(X * (mu * (1- mu)), X) d <-solve(-H, grad) # ascent direction b <- b +as.numeric(d)cat(sprintf(" k=%d loglik=%.6f |grad|=%.2e\n", k, loglik(b), max(abs(grad))))if (max(abs(d)) <1e-10) break}# Standard errors from the inverse negative Hessian at the optimummu <-plogis(as.numeric(X %*% b))H <--crossprod(X * (mu * (1- mu)), X)se <-sqrt(diag(solve(-H)))# Compare with glm() — same algorithm (IRLS = Fisher scoring)fit <-glm(y ~ x1 + x2 + x3, data = dat, family = binomial)tab <-data.frame(true =c(-0.5, 1.2, -0.8, 0.6),newton =round(b, 6),glm =round(as.numeric(coef(fit)), 6),se_newton =round(se, 6),se_glm =round(as.numeric(summary(fit)$coef[, 2]), 6))rownames(tab) <-c("(Intercept)", "x1", "x2", "x3")print(tab)
Code
import numpy as npimport pandas as pdimport statsmodels.api as smdat = pd.read_csv("../data/numerical-logit.csv") # same CSV as RX = np.column_stack([np.ones(len(dat)), dat[["x1", "x2", "x3"]].to_numpy()])y = dat["y"].to_numpy(dtype=float)def loglik(b): eta = X @ breturn np.sum(y * eta - np.logaddexp(0.0, eta)) # stable (Part I)# Newton-Raphson: solve H d = -grad each step (Part II)b = np.zeros(4)print("Newton-Raphson iterations:")for k inrange(1, 26): mu =1/ (1+ np.exp(-X @ b)) grad = X.T @ (y - mu) H =-(X * (mu * (1- mu))[:, None]).T @ X d = np.linalg.solve(-H, grad) b += dprint(f" k={k} loglik={loglik(b):.6f} |grad|={np.max(np.abs(grad)):.2e}")if np.max(np.abs(d)) <1e-10:break# Standard errors from the inverse negative Hessianmu =1/ (1+ np.exp(-X @ b))H =-(X * (mu * (1- mu))[:, None]).T @ Xse = np.sqrt(np.diag(np.linalg.inv(-H)))# Compare with statsmodels Logitfit = sm.Logit(y, X).fit(disp=0)true_b = [-0.5, 1.2, -0.8, 0.6]names = ["const", "x1", "x2", "x3"]print(f"\n{'param':<8}{'true':>7}{'newton':>10}{'sm.Logit':>10}{'se_newton':>10}{'se_sm':>8}")for nm, tb, bb, sb, fb, fs inzip(names, true_b, b, se, fit.params, fit.bse):print(f"{nm:<8}{tb:7.2f}{bb:10.6f}{fb:10.6f}{sb:10.6f}{fs:8.6f}")
The Same MLE with Library Optimizers — BFGS & Nelder–Mead
Real projects hand the objective to a general optimizer. Two things matter: supply the gradient when you can, and know that different algorithms take very different iteration counts to the same optimum.
Code
dat <-read.csv("../data/numerical-logit.csv")X <-cbind(1, as.matrix(dat[, c("x1", "x2", "x3")]))y <- dat$ylog1pexp <-function(eta) pmax(eta, 0) +log1p(exp(-abs(eta)))negll <-function(b) -sum(y * (X %*% b) -log1pexp(X %*% b))neggr <-function(b) -as.numeric(crossprod(X, y -plogis(as.numeric(X %*% b))))# Before optimising: verify the analytic gradient against numDerivb_test <-c(0.1, -0.2, 0.3, 0.05)cat("gradient check |analytic - numerical| :",format(max(abs(neggr(b_test) - numDeriv::grad(negll, b_test))), digits =3), "\n\n")# BFGS with analytic gradientfit_bfgs <-optim(rep(0, 4), fn = negll, gr = neggr, method ="BFGS",control =list(maxit =500), hessian =TRUE)# Nelder-Mead: derivative-free, many more function evaluationsfit_nm <-optim(rep(0, 4), fn = negll, method ="Nelder-Mead",control =list(maxit =5000))cat(sprintf("BFGS : negll = %.6f, evals(fn, gr) = %d, %d, conv = %d\n", fit_bfgs$value, fit_bfgs$counts[1], fit_bfgs$counts[2], fit_bfgs$convergence))cat(sprintf("Nelder-Mead : negll = %.6f, evals(fn) = %d, conv = %d\n\n", fit_nm$value, fit_nm$counts[1], fit_nm$convergence))# SEs from optim's numerical Hessian vs glmse_opt <-sqrt(diag(solve(fit_bfgs$hessian)))fit_glm <-glm(y ~ x1 + x2 + x3, data = dat, family = binomial)print(round(data.frame(bfgs = fit_bfgs$par,nm = fit_nm$par,glm =as.numeric(coef(fit_glm)),se_bfgs = se_opt,se_glm =as.numeric(summary(fit_glm)$coef[, 2])), 5))
Code
import numpy as npimport pandas as pdfrom scipy.optimize import minimize, check_gradimport statsmodels.api as smdat = pd.read_csv("../data/numerical-logit.csv")X = np.column_stack([np.ones(len(dat)), dat[["x1", "x2", "x3"]].to_numpy()])y = dat["y"].to_numpy(dtype=float)def negll(b): eta = X @ breturn-np.sum(y * eta - np.logaddexp(0.0, eta))def neggr(b):return-(X.T @ (y -1/ (1+ np.exp(-X @ b))))# Verify the analytic gradient before trusting any optimizer outputprint("gradient check (scipy.check_grad):",f"{check_grad(negll, neggr, np.array([0.1, -0.2, 0.3, 0.05])):.3e}\n")fit_bfgs = minimize(negll, np.zeros(4), jac=neggr, method="BFGS")fit_nm = minimize(negll, np.zeros(4), method="Nelder-Mead", options={"maxiter": 5000, "xatol": 1e-8, "fatol": 1e-8})print(f"BFGS : negll = {fit_bfgs.fun:.6f}, n_fev = {fit_bfgs.nfev}, "f"n_jev = {fit_bfgs.njev}, success = {fit_bfgs.success}")print(f"Nelder-Mead : negll = {fit_nm.fun:.6f}, n_fev = {fit_nm.nfev}, "f"success = {fit_nm.success}\n")# SEs from BFGS's inverse-Hessian approximation vs statsmodelsse_bfgs = np.sqrt(np.diag(fit_bfgs.hess_inv))fit_sm = sm.Logit(y, X).fit(disp=0)print(f"{'param':<7}{'bfgs':>10}{'nelder':>10}{'sm':>10}{'se_bfgs':>9}{'se_sm':>9}")for i, nm inenumerate(["const", "x1", "x2", "x3"]):print(f"{nm:<7}{fit_bfgs.x[i]:10.5f}{fit_nm.x[i]:10.5f} "f"{fit_sm.params[i]:10.5f}{se_bfgs[i]:9.5f}{fit_sm.bse[i]:9.5f}")
Gradient Descent — Learning Rate Is Everything
Minimise the OLS loss \(Q(\boldsymbol\beta) = \frac{1}{2n}\|\mathbf{y}-\mathbf{X}\boldsymbol\beta\|^2\) by gradient descent — the algorithm behind all of deep learning (and Part VIII’s torch code). Convergence hinges on the step size \(\gamma\) relative to the curvature \(L = \sigma_{\max}^2/n\).
Code
set.seed(14159)n <-500X <-cbind(1, rnorm(n), rnorm(n))beta_true <-c(1, 2, -1.5)y <-as.numeric(X %*% beta_true +rnorm(n))L <-max(eigen(crossprod(X) / n)$values) # Lipschitz constant of the gradientgd <-function(gamma, iters =60) { b <-rep(0, 3); path <-numeric(iters)for (k inseq_len(iters)) { b <- b - gamma *as.numeric(crossprod(X, X %*% b - y)) / n path[k] <-0.5*mean((y - X %*% b)^2) } path}rates <-c(0.05, 0.5, 1.9, 2.1) / L # fractions of the stability limit 2/Llabs <-c("0.05/L (crawls)", "0.5/L (good)", "1.9/L (fast)", "2.1/L (diverges!)")dd <-data.frame(iter =rep(1:60, 4),loss =unlist(lapply(rates, gd)),gamma =factor(rep(labs, each =60), levels = labs))cat(sprintf("stability limit: gamma < 2/L = %.4f\n", 2/ L))ggplot(dd, aes(iter, pmin(loss, 1e3), colour = gamma)) +geom_line(linewidth =1.2) +scale_y_log10() +scale_colour_manual(values =c("#BA7517", "#185FA5", "#1D9E75", "#C0132C")) +labs(x ="iteration", y ="loss (log scale, capped)",title ="Gradient descent on OLS loss: four learning rates",colour ="step size")
Code
import numpy as npimport matplotlib.pyplot as pltrng = np.random.default_rng(14159)n =500X = np.column_stack([np.ones(n), rng.standard_normal(n), rng.standard_normal(n)])beta_true = np.array([1.0, 2.0, -1.5])y = X @ beta_true + rng.standard_normal(n)L = np.linalg.eigvalsh(X.T @ X / n).max()print(f"stability limit: gamma < 2/L = {2/L:.4f}")def gd(gamma, iters=60): b = np.zeros(3); path = []for _ inrange(iters): b -= gamma * X.T @ (X @ b - y) / n path.append(0.5* np.mean((y - X @ b)**2))return np.array(path)rates = [0.05/ L, 0.5/ L, 1.9/ L, 2.1/ L]labs = ["0.05/L (crawls)", "0.5/L (good)", "1.9/L (fast)", "2.1/L (diverges!)"]cols = ["#BA7517", "#185FA5", "#1D9E75", "#C0132C"]fig, ax = plt.subplots(figsize=(9, 4.5))for g, lab, c inzip(rates, labs, cols): ax.plot(np.minimum(gd(g), 1e3), label=lab, color=c, lw=2)ax.set_yscale("log")ax.set_xlabel("iteration"); ax.set_ylabel("loss (log scale, capped)")ax.set_title("Gradient descent on OLS loss: four learning rates")ax.legend(); plt.tight_layout(); plt.show()
Part V — Numerical Differentiation & Integration
Finite Differences, Quadrature, and Monte Carlo
Numerical Differentiation — The Step-Size Dilemma
Forward and central finite differences approximate derivatives with truncation error that shrinks with \(h\):
But shrinking \(h\) revives Part I: \(f(x+h) - f(x)\) is a subtraction of nearly equal numbers — rounding error grows as \(\varepsilon_{\text{mach}}/h\). Total error is U-shaped, with the optimum at
Truncation error and rounding error pull in opposite directions — you can never have both. This is why numDeriv uses Richardson extrapolation and why marginal effects, delta-method SEs, and numerical Hessians are all sensitive to step size.
Code
# d/dx exp(x) at x = 1: true value is ef <- exp; x0 <-1; true_d <-exp(1)h <-10^seq(-1, -15, by =-0.25)err_fwd <-abs((f(x0 + h) -f(x0)) / h - true_d) / true_derr_ctr <-abs((f(x0 + h) -f(x0 - h)) / (2* h) - true_d) / true_ddd <-data.frame(h =rep(h, 2), err =c(err_fwd, err_ctr),scheme =rep(c("forward O(h)", "central O(h^2)"), each =length(h)))ggplot(dd, aes(h, err, colour = scheme)) +geom_line(linewidth =1.1) +geom_point(size =1.6) +scale_x_log10() +scale_y_log10() +scale_colour_manual(values =c("#185FA5", "#C0132C")) +geom_vline(xintercept =c(sqrt(.Machine$double.eps), .Machine$double.eps^(1/3)),linetype ="dashed", colour ="grey55") +labs(x ="step size h (log)", y ="relative error (log)",title ="The U-curve: truncation error vs rounding error",colour =NULL)
Three routes: adaptive quadrature (fast, ~machine precision for smooth 1-D integrands), fixed-node Simpson’s rule (\(O(h^4)\)), and Monte Carlo (slow \(O(1/\sqrt{S})\) — but dimension-free, which is why simulated MLE and mixed logit use it).
Code
demand <-function(p) 100* p^(-1.5)# 1. Adaptive quadrature: handles the infinite limit directlyq_adapt <-integrate(demand, lower =4, upper =Inf)cat(sprintf("integrate() : CS = %.10f (abs.error < %.1e)\n", q_adapt$value, q_adapt$abs.error))# 2. Simpson's rule on [4, 4000] with n = 1000 intervalssimpson <-function(f, a, b, n =1000) { h <- (b - a) / n x <-seq(a, b, length.out = n +1) h /3*sum(f(x) *c(1, rep(c(4, 2), (n -2) /2), 4, 1))}cat(sprintf("Simpson [4,4000] : CS = %.10f\n", simpson(demand, 4, 4000)))# 3. Monte Carlo with importance sampling: p = 4/U^2, U ~ Uniform(0,1)set.seed(14159)for (S inc(1e3, 1e5, 1e7)) { u <-runif(S) p <-4/ u^2# maps (0,1) onto (4, Inf) w <-8/ u^3# Jacobian |dp/du| cs <-mean(demand(p) * w)cat(sprintf("Monte Carlo S=%.0e: CS = %.6f (error %.1e)\n", S, cs, abs(cs -100)))}cat("\ntrue value: 100 — MC error shrinks at rate 1/sqrt(S)\n")
Code
import numpy as npfrom scipy.integrate import quad, simpsondemand =lambda p: 100* p ** (-1.5)# 1. Adaptive quadrature: handles the infinite limit directlyval, err = quad(demand, 4, np.inf)print(f"scipy quad : CS = {val:.10f} (abs.error < {err:.1e})")# 2. Simpson's rule on [4, 4000] with 1001 nodesx = np.linspace(4, 4000, 1001)print(f"Simpson [4,4000] : CS = {simpson(demand(x), x=x):.10f}")# 3. Monte Carlo with importance sampling: p = 4/U^2rng = np.random.default_rng(14159)for S in [10**3, 10**5, 10**7]: u = rng.uniform(size=S) p =4/ u**2 w =8/ u**3 cs = np.mean(demand(p) * w)print(f"Monte Carlo S=1e{int(np.log10(S))}: CS = {cs:.6f} (error {abs(cs-100):.1e})")print("\ntrue value: 100 — MC error shrinks at rate 1/sqrt(S)")
Halving \(h\) cuts Euler’s error by 2 but RK4’s by 16. Production solvers (deSolve::ode, solve_ivp) add adaptive step size and stiffness detection (LSODA switches methods automatically).
⚠️ Stiffness in One Paragraph
A system is stiff when it mixes very fast and very slow dynamics (eigenvalues of the Jacobian differ by orders of magnitude) — e.g. a fast-adjusting price coupled with slow capital accumulation.
Explicit methods (Euler, RK4) must then take tiny steps dictated by the fastest mode — even after it has died out — or the solution oscillates and explodes. Implicit methods (backward Euler, BDF) solve a small nonlinear system per step (Part III!) and remain stable at large steps.
Practice: R's deSolve::ode() default lsoda and Python's solve_ivp(method="LSODA") detect stiffness and switch automatically — one more reason to prefer library solvers over hand-rolled loops in production.
The Solow Growth Model — Setup
Capital per effective worker \(k(t)\) with Cobb–Douglas production \(f(k) = k^\alpha\):
Plan: solve on \(t \in [0, 100]\) with hand-rolled Euler and RK4 at a coarse step \(h = 5\), plus the adaptive library solver; compare against the exact path.
Solow Model — Euler vs RK4 vs Library Solver
Code
s <-0.25; alpha <-0.36; ngd <-0.1; k0 <-1f <-function(t, k) s * k^alpha - ngd * k# Exact solution (Bernoulli ODE) for the error benchmarkk_exact <-function(t) { A <- s / ngd (A + (k0^(1- alpha) - A) *exp(-(1- alpha) * ngd * t))^(1/ (1- alpha))}# Euler and RK4 with the SAME coarse step h = 5h <-5; times <-seq(0, 100, by = h); m <-length(times)k_euler <- k_rk4 <-numeric(m); k_euler[1] <- k_rk4[1] <- k0for (i in1:(m -1)) { t <- times[i]# Euler: one slope k_euler[i +1] <- k_euler[i] + h *f(t, k_euler[i])# RK4: four slopes, weighted average k1 <-f(t, k_rk4[i]) k2 <-f(t + h /2, k_rk4[i] + h /2* k1) k3 <-f(t + h /2, k_rk4[i] + h /2* k2) k4 <-f(t + h, k_rk4[i] + h * k3) k_rk4[i +1] <- k_rk4[i] + h /6* (k1 +2* k2 +2* k3 + k4)}# Library solver: adaptive LSODA via deSolvesol <- deSolve::ode(y =c(k = k0), times = times,func =function(t, y, p) list(f(t, y)), parms =NULL)cat(sprintf("steady state k* (closed form) : %.6f\n", (s / ngd)^(1/ (1- alpha))))cat(sprintf("max |error| Euler (h = 5) : %.2e\n", max(abs(k_euler -k_exact(times)))))cat(sprintf("max |error| RK4 (h = 5) : %.2e\n", max(abs(k_rk4 -k_exact(times)))))cat(sprintf("max |error| deSolve lsoda : %.2e\n", max(abs(sol[, "k"] -k_exact(times)))))dd <-data.frame(t =rep(times, 3),k =c(k_euler, k_rk4, k_exact(times)),method =rep(c("Euler h=5", "RK4 h=5", "exact"), each = m))ggplot(dd, aes(t, k, colour = method, linetype = method)) +geom_line(linewidth =1.1) +geom_hline(yintercept = (s / ngd)^(1/ (1- alpha)),linetype ="dotted", colour ="grey55") +scale_colour_manual(values =c("#1D9E75", "#C0132C", "#185FA5")) +annotate("text", x =90, y =4.35, label ="k*", colour ="grey40", size =5) +labs(x ="time", y ="capital per effective worker k(t)",title ="Solow transition path: Euler visibly off at coarse steps, RK4 on top of exact",colour =NULL, linetype =NULL)
Code
import numpy as npimport matplotlib.pyplot as pltfrom scipy.integrate import solve_ivps, alpha, ngd, k0 =0.25, 0.36, 0.1, 1.0f =lambda t, k: s * k**alpha - ngd * kdef k_exact(t): A = s / ngdreturn (A + (k0**(1- alpha) - A) * np.exp(-(1- alpha) * ngd * t))**(1/ (1- alpha))# Euler and RK4 with the SAME coarse step h = 5h =5.0times = np.arange(0, 100+ h, h)m =len(times)k_euler = np.empty(m); k_rk4 = np.empty(m)k_euler[0] = k_rk4[0] = k0for i inrange(m -1): t = times[i] k_euler[i +1] = k_euler[i] + h * f(t, k_euler[i]) k1 = f(t, k_rk4[i]) k2 = f(t + h /2, k_rk4[i] + h /2* k1) k3 = f(t + h /2, k_rk4[i] + h /2* k2) k4 = f(t + h, k_rk4[i] + h * k3) k_rk4[i +1] = k_rk4[i] + h /6* (k1 +2* k2 +2* k3 + k4)# Library solver: adaptive RK45 via solve_ivpsol = solve_ivp(f, [0, 100], [k0], t_eval=times, rtol=1e-8, atol=1e-10)print(f"steady state k* (closed form): {(s/ngd)**(1/(1-alpha)):.6f}")print(f"max |error| Euler (h = 5) : {np.max(np.abs(k_euler - k_exact(times))):.2e}")print(f"max |error| RK4 (h = 5) : {np.max(np.abs(k_rk4 - k_exact(times))):.2e}")print(f"max |error| solve_ivp RK45 : {np.max(np.abs(sol.y[0] - k_exact(times))):.2e}")fig, ax = plt.subplots(figsize=(9, 4.5))ax.plot(times, k_exact(times), color="#185FA5", lw=2.5, label="exact")ax.plot(times, k_euler, "o--", color="#C0132C", lw=1.5, ms=4, label="Euler h=5")ax.plot(times, k_rk4, "s-", color="#1D9E75", lw=1.5, ms=4, label="RK4 h=5")ax.axhline((s/ngd)**(1/(1-alpha)), ls=":", color="grey")ax.set_xlabel("time"); ax.set_ylabel("capital per effective worker k(t)")ax.set_title("Solow transition path: Euler visibly off at coarse steps, RK4 on top of exact")ax.legend(); plt.tight_layout(); plt.show()
Convergence Order in Action — Halve \(h\), Watch the Error
Theory says: Euler error \(\propto h\), RK4 error \(\propto h^4\). Verify it empirically on the Solow model — the log-log error curve’s slope is the order of the method.
Code
s <-0.25; alpha <-0.36; ngd <-0.1; k0 <-1f <-function(t, k) s * k^alpha - ngd * kk_exact_T <- (s / ngd + (k0^(1- alpha) - s / ngd) *exp(-(1- alpha) * ngd *50))^(1/ (1- alpha))solve_scheme <-function(h, scheme) { times <-seq(0, 50, by = h); k <- k0for (i in1:(length(times) -1)) { t <- times[i]if (scheme =="euler") { k <- k + h *f(t, k) } else { k1 <-f(t, k); k2 <-f(t + h/2, k + h/2* k1) k3 <-f(t + h/2, k + h/2* k2); k4 <-f(t + h, k + h * k3) k <- k + h /6* (k1 +2*k2 +2*k3 + k4) } }abs(k - k_exact_T)}hs <-c(5, 2.5, 1.25, 0.625, 0.3125)err <-data.frame(h =rep(hs, 2),error =c(sapply(hs, solve_scheme, scheme ="euler"),sapply(hs, solve_scheme, scheme ="rk4")),method =rep(c("Euler", "RK4"), each =length(hs)))# Empirical order = slope of log(error) on log(h)for (mth inc("Euler", "RK4")) { e <- err[err$method == mth, ]cat(sprintf("%-5s empirical order: %.2f\n", mth,coef(lm(log(error) ~log(h), data = e))[2]))}ggplot(err, aes(h, error, colour = method)) +geom_line(linewidth =1.1) +geom_point(size =2.5) +scale_x_log10() +scale_y_log10() +scale_colour_manual(values =c("#C0132C", "#1D9E75")) +labs(x ="step size h (log)", y ="error at t = 50 (log)",title ="Slopes on log-log axes reveal the order: ~1 for Euler, ~4 for RK4",colour =NULL)
Code
import numpy as nps, alpha, ngd, k0 =0.25, 0.36, 0.1, 1.0f =lambda t, k: s * k**alpha - ngd * kA = s / ngdk_exact_T = (A + (k0**(1- alpha) - A) * np.exp(-(1- alpha) * ngd *50))**(1/ (1- alpha))def solve_scheme(h, scheme): times = np.arange(0, 50+ h, h) k = k0for t in times[:-1]:if scheme =="euler": k = k + h * f(t, k)else: k1 = f(t, k); k2 = f(t + h/2, k + h/2* k1) k3 = f(t + h/2, k + h/2* k2); k4 = f(t + h, k + h * k3) k = k + h /6* (k1 +2*k2 +2*k3 + k4)returnabs(k - k_exact_T)hs = np.array([5, 2.5, 1.25, 0.625, 0.3125])for scheme in ["euler", "rk4"]: errs = np.array([solve_scheme(h, scheme) for h in hs]) order = np.polyfit(np.log(hs), np.log(errs), 1)[0]print(f"{scheme:5s} errors:", " ".join(f"{e:.2e}"for e in errs),f"| empirical order: {order:.2f}")
Part VII — CPU Parallelization
Processes vs Threads · Amdahl’s Law · Monte Carlo at Scale
Econometric workloads are embarrassingly parallel at the replication level: Monte Carlo draws, bootstrap resamples, cross-validation folds, grid points. Each unit is independent — perfect for process-based parallelism.
Do not stack them: 6 worker processes each spawning 12 BLAS threads = 72 threads on 16 logical CPUs — slower than serial. Pin BLAS to 1 thread inside workers.
If a fraction \(f\) of runtime is parallelizable over \(c\) cores, the maximal speedup is
\(f = 0.9\): at most 10×, ever — even with 1000 cores
\(f = 0.5\): at most 2× — parallelizing is barely worth it
The serial part (data loading, result assembly, rendering) always wins in the end
Measure before parallelizing: profile, find \(f\), then decide.
Parallel random numbers are a trap: workers must get independent, reproducible streams, not copies of the same seed.
R: mclapply(..., mc.set.seed = TRUE) with RNGkind("L'Ecuyer-CMRG") and set.seed(14159) — statistically independent streams per worker
Python: np.random.default_rng(SeedSequence(14159).spawn(n)) — one child generator per task
Never set.seed(14159)inside every worker: all workers then produce identical draws and your “10,000 replications” are 6 unique ones
Parallel Monte Carlo — Bootstrap of a Median
The median has no simple SE formula — bootstrap it. \(B = 2000\) resamples of \(n = 10{,}000\) incomes (lognormal): independent tasks, ideal for parallel workers. We time serial vs 6 processes.
Code
library(parallel)set.seed(14159)income <-rlnorm(1e4, meanlog =10, sdlog =0.75) # skewed income dataB <-2000boot_median <-function(b) {# each task: resample and add a small deterministic workload idx <-sample.int(length(income), replace =TRUE)for (j in1:30) m <-median(income[idx]) # repeat to make work visible m}# Serialt_serial <-system.time(res_s <-lapply(1:B, boot_median))["elapsed"]# Parallel: 6 forked processes with proper parallel RNG streamsRNGkind("L'Ecuyer-CMRG")set.seed(14159)t_par <-system.time( res_p <-mclapply(1:B, boot_median, mc.cores =6, mc.set.seed =TRUE))["elapsed"]ci_s <-quantile(unlist(res_s), c(0.025, 0.975))ci_p <-quantile(unlist(res_p), c(0.025, 0.975))cat(sprintf("serial : %6.2f s 95%% CI for median: [%.0f, %.0f]\n", t_serial, ci_s[1], ci_s[2]))cat(sprintf("parallel: %6.2f s 95%% CI for median: [%.0f, %.0f]\n", t_par, ci_p[1], ci_p[2]))cat(sprintf("speedup : %.1fx on 6 cores\n", t_serial / t_par))cat(sprintf("implied parallel fraction (Amdahl): f = %.2f\n", (1- t_par / t_serial) / (1-1/6)))
Code
import warningswarnings.filterwarnings("ignore")import numpy as npimport timefrom joblib import Parallel, delayedrng = np.random.default_rng(14159)income = rng.lognormal(mean=10, sigma=0.75, size=10_000)B =2000# One independent child generator per bootstrap task (reproducible streams)child_seeds = np.random.SeedSequence(14159).spawn(B)def boot_median(seed_seq): r = np.random.default_rng(seed_seq) idx = r.integers(0, income.size, income.size)for _ inrange(30): # repeat to make work visible m = np.median(income[idx])return mt0 = time.perf_counter()res_s = [boot_median(s) for s in child_seeds]t_serial = time.perf_counter() - t0t0 = time.perf_counter()res_p = Parallel(n_jobs=6)(delayed(boot_median)(s) for s in child_seeds)t_par = time.perf_counter() - t0ci_s = np.percentile(res_s, [2.5, 97.5])ci_p = np.percentile(res_p, [2.5, 97.5])print(f"serial : {t_serial:6.2f} s 95% CI for median: [{ci_s[0]:.0f}, {ci_s[1]:.0f}]")print(f"parallel: {t_par:6.2f} s 95% CI for median: [{ci_p[0]:.0f}, {ci_p[1]:.0f}]")print(f"speedup : {t_serial / t_par:.1f}x on 6 cores")print(f"implied parallel fraction (Amdahl): f = {(1- t_par/t_serial) / (1-1/6):.2f}")
Amdahl’s Law — Visualised
Code
cores <-1:64dd <-expand.grid(c = cores, f =c(0.50, 0.75, 0.90, 0.95, 0.99))dd$speedup <-1/ ((1- dd$f) + dd$f / dd$c)dd$f_lab <-factor(sprintf("f = %.2f (max %.0fx)", dd$f, 1/ (1- dd$f)))ggplot(dd, aes(c, speedup, colour = f_lab)) +geom_line(linewidth =1.2) +geom_abline(slope =1, intercept =0, linetype ="dashed", colour ="grey55") +scale_colour_manual(values =c("#C0132C", "#BA7517", "#185FA5","#1D9E75", "#7B3FA0")) +annotate("text", x =55, y =58, label ="ideal linear", colour ="grey40",angle =38, size =4) +coord_cartesian(ylim =c(0, 64)) +labs(x ="number of cores", y ="speedup S(c)",title ="Amdahl's law: the serial fraction caps the speedup",colour ="parallel fraction")
Code
import numpy as npimport matplotlib.pyplot as pltcores = np.arange(1, 65)fracs = [0.50, 0.75, 0.90, 0.95, 0.99]cols = ["#C0132C", "#BA7517", "#185FA5", "#1D9E75", "#7B3FA0"]fig, ax = plt.subplots(figsize=(9, 4.5))for f, c inzip(fracs, cols): ax.plot(cores, 1/ ((1- f) + f / cores), color=c, lw=2, label=f"f = {f:.2f} (max {1/(1-f):.0f}x)")ax.plot(cores, cores, "--", color="grey", label="ideal linear")ax.set_ylim(0, 64)ax.set_xlabel("number of cores"); ax.set_ylabel("speedup S(c)")ax.set_title("Amdahl's law: the serial fraction caps the speedup")ax.legend(); plt.tight_layout(); plt.show()
A CPU has ~8–16 powerful cores optimised for sequential logic. A GPU has thousands of simple cores (an RTX 3050: 2560 CUDA cores) optimised for doing the same operation on many numbers at once — exactly the shape of linear algebra.
\[\text{matmul: } \mathbf{C} = \mathbf{A}\mathbf{B} \quad \text{— } n^3 \text{ multiply-adds, all independent across } (i,j) \text{ cells}\]
Workload
GPU payoff
Large dense matmul, batched OLS
10–100×
Gradient-based MLE with big \(n\)
large
Neural nets, causal ML learners
the reason torch exists
Loops with branching, small data
slower than CPU
Anything dominated by data transfer
slower than CPU
Transfer tax — data must cross the PCIe bus to GPU memory and back:
\[T_{\text{total}} = \underbrace{T_{\text{transfer}}}_{\text{once per dataset}} + \underbrace{T_{\text{compute}}}_{\text{the fast part}}\]
Move data once, keep it on the device, do many operations there, bring back only results. A single small matmul is not worth the trip.
Precision tax — consumer GPUs are built for float32; float64 units are few:
float32: ~7 correct digits — fine for gradient steps, dangerous for near-singular solves
Part I’s rules bite harder: cancellation happens 10⁹ times sooner in float32
Rule: explore in float32, verify in float64, compare with a CPU reference
Same API in both languages (R torch is a native re-implementation, not a Python wrapper):
Operation
R
Python
Create tensor
torch_tensor(x, device = "cuda")
torch.tensor(x, device="cuda")
Check GPU
cuda_is_available()
torch.cuda.is_available()
Matrix product
torch_matmul(A, B) / A$matmul(B)
A @ B
Track gradients
x$requires_grad_(TRUE)
x.requires_grad_(True)
Backpropagate
loss$backward()
loss.backward()
Read gradient
x$grad
x.grad
Back to CPU/R
as.numeric(x$cpu())
x.cpu().numpy()
A systems pitfall: R torch and Python torch cannot live in one process — two libtorch builds collide on the same shared library. This deck’s Python tabs drive the real GPU, so each R torch snippet is dispatched to an isolated subprocess via callr::r() (the r_torch() helper in setup). Same code, its own process.
Matrix Multiplication — CPU vs GPU Benchmark
\(4000 \times 4000\) matrix products: ~\(1.3 \times 10^{11}\) floating-point operations each. This machine: 16-thread CPU vs an NVIDIA RTX 3050 (CUDA). Note the synchronisation: GPU calls return before the work finishes — always synchronize() before reading the clock.
Code
import torch, timen =4000A = torch.randn(n, n) # float32 on CPUB = torch.randn(n, n)def bench(fn, sync=None, reps=3): fn() # warm-up (CUDA kernels compile lazily)if sync: sync() t0 = time.perf_counter()for _ inrange(reps): fn()if sync: sync()return (time.perf_counter() - t0) / repst_cpu = bench(lambda: A @ B)print(f"CPU (float32, 16 threads): {t_cpu*1000:8.1f} ms")if torch.cuda.is_available():print("GPU:", torch.cuda.get_device_name(0))# Transfer tax: moving the inputs to the device, once t0 = time.perf_counter() Ag, Bg = A.cuda(), B.cuda() torch.cuda.synchronize()print(f"transfer to GPU (once) : {(time.perf_counter()-t0)*1000:8.1f} ms") t_gpu = bench(lambda: Ag @ Bg, sync=torch.cuda.synchronize)print(f"GPU (float32) : {t_gpu*1000:8.1f} ms")print(f"speedup (compute only) : {t_cpu / t_gpu:8.1f}x")# The precision tax: same product in float64 on the GPU Ad, Bd = Ag.double(), Bg.double() t_gpu64 = bench(lambda: Ad @ Bd, sync=torch.cuda.synchronize)print(f"GPU (float64) : {t_gpu64*1000:8.1f} ms "f"({t_gpu64/t_gpu:.0f}x slower than float32)")else:print("No CUDA device found — GPU timings skipped")
Code
# Same benchmark with R torch. The API mirrors Python torch one-to-one;# device = "cuda" works identically when R torch has a CUDA backend installed.# r_torch() runs this in an isolated subprocess (see the setup note).r_torch(function() {suppressPackageStartupMessages(library(torch)) dev <-if (cuda_is_available()) "cuda"else"cpu"cat("R torch device:", dev,if (dev =="cpu") "(this machine's R backend is CPU-only; Python tab shows the GPU)"else"", "\n\n") n <-2000 A <-torch_randn(n, n, device = dev) # float32 tensors B <-torch_randn(n, n, device = dev)# torch matmul (all 16 threads via its own BLAS)invisible(torch_matmul(A, B)) # warm-up t_torch <-system.time(for (i in1:3) invisible(torch_matmul(A, B)))["elapsed"] /3# base R %*% in float64 for reference A_r <-matrix(rnorm(n * n), n, n); B_r <-matrix(rnorm(n * n), n, n) t_base <-system.time(invisible(A_r %*% B_r))["elapsed"]cat(sprintf("torch matmul (float32, %s): %7.1f ms\n", dev, 1000* t_torch))cat(sprintf("base R %%*%% (float64, cpu): %7.1f ms\n", 1000* t_base))# The tensor API in brief: same ops, chainable, autograd-ready x <-torch_tensor(c(1, 2, 3), device = dev)cat("\nx * 2 + 1 ->", as.numeric((x *2+1)$cpu()), "\n")})
Autograd — Derivatives Without Deriving
torch records every tensor operation in a computational graph and applies the chain rule backwards — exact gradients (to float precision), no pencil work, no finite-difference U-curve (Part V):
This is reverse-mode automatic differentiation: one backward pass costs about as much as one forward pass, regardless of the number of parameters — why it powers all of deep learning and, increasingly, structural econometrics.
Code
# r_torch() runs this R torch snippet in an isolated subprocess (setup note).r_torch(function() {suppressPackageStartupMessages(library(torch))# d/dx of f(x) = x^3 + 2x at x = 2: analytic answer 3*4 + 2 = 14 x <-torch_tensor(2, requires_grad =TRUE) f <- x^3+2* x f$backward()cat("autograd df/dx at x=2 :", as.numeric(x$grad), " (analytic: 14)\n")# Gradient of a vector function: OLS loss gradient in one lineset.seed(14159) X <-torch_tensor(cbind(1, matrix(rnorm(200), 100, 2))) y <-torch_tensor(rnorm(100)) b <-torch_zeros(3, requires_grad =TRUE) loss <-torch_mean((y -torch_matmul(X, b))^2) loss$backward()cat("autograd gradient:", round(as.numeric(b$grad), 6), "\n")# Check against the analytic gradient -2/n X'(y - Xb) at b = 0 g_analytic <--2*as.numeric(torch_matmul(torch_t(X), y)) /100cat("analytic gradient:", round(g_analytic, 6), "\n")})
Code
import torch# d/dx of f(x) = x^3 + 2x at x = 2: analytic answer 14x = torch.tensor(2.0, requires_grad=True)f = x**3+2* xf.backward()print(f"autograd df/dx at x=2 : {x.grad.item():.1f} (analytic: 14)")# Gradient of a vector function: OLS loss gradient in one lineg = torch.Generator().manual_seed(14159)X = torch.cat([torch.ones(100, 1), torch.randn(100, 2, generator=g)], dim=1)y = torch.randn(100, generator=g)b = torch.zeros(3, requires_grad=True)loss = torch.mean((y - X @ b)**2)loss.backward()print("autograd gradient:", b.grad.numpy().round(6))g_analytic = (-2* X.T @ y /100).numpy()print("analytic gradient:", g_analytic.round(6))
Logit MLE on the GPU — Autograd Gradient Descent
The Part IV logit, re-estimated with torch: define the stable negative log-likelihood, let autograd differentiate it, take gradient steps on the device. Same data (../data/numerical-logit.csv), so the coefficients must match Newton’s to float precision.
Code
import torchimport pandas as pdimport numpy as npdat = pd.read_csv("../data/numerical-logit.csv")dev ="cuda"if torch.cuda.is_available() else"cpu"# Move the data to the device ONCE, in float64 for full precisionX = torch.tensor(np.column_stack([np.ones(len(dat)), dat[["x1", "x2", "x3"]].to_numpy()]), dtype=torch.float64, device=dev)y = torch.tensor(dat["y"].to_numpy(), dtype=torch.float64, device=dev)b = torch.zeros(4, dtype=torch.float64, device=dev, requires_grad=True)opt = torch.optim.Adam([b], lr=0.05) # adaptive gradient descentfor k inrange(1, 601): opt.zero_grad() eta = X @ b negll =-(y * eta - torch.logaddexp(torch.zeros_like(eta), eta)).sum() negll.backward() # autograd: exact gradient opt.step()if k %150==0:print(f" iter {k:4d} negll = {negll.item():.6f} "f"|grad| = {b.grad.abs().max().item():.2e}")b_hat = b.detach().cpu().numpy()print(f"\ndevice used: {dev}")print("torch Adam :", b_hat.round(6))print("true beta : [-0.5, 1.2, -0.8, 0.6]")print("(compare the Newton-Raphson estimates in Part IV — identical to ~6 digits)")
Code
# r_torch() runs this R torch snippet in an isolated subprocess (setup note).r_torch(function() {suppressPackageStartupMessages(library(torch)) dat <-read.csv("../data/numerical-logit.csv") dev <-if (cuda_is_available()) "cuda"else"cpu" X <-torch_tensor(cbind(1, as.matrix(dat[, c("x1", "x2", "x3")])),dtype =torch_float64(), device = dev) y <-torch_tensor(dat$y, dtype =torch_float64(), device = dev) b <-torch_zeros(4, dtype =torch_float64(), device = dev, requires_grad =TRUE) opt <-optim_adam(list(b), lr =0.05)for (k in1:600) { opt$zero_grad() eta <-torch_matmul(X, b) negll <--torch_sum(y * eta -torch_logaddexp(torch_zeros_like(eta), eta)) negll$backward() opt$step()if (k %%150==0)cat(sprintf(" iter %4d negll = %.6f\n", k, negll$item())) } b_hat <-as.numeric(b$detach()$cpu())cat("\ndevice used:", dev, "\n")cat("torch Adam :", round(b_hat, 6), "\n")cat("true beta : -0.5 1.2 -0.8 0.6\n")cat("glm() check:", round(as.numeric(coef(glm(y ~ x1 + x2 + x3,data = dat, family = binomial))), 6), "\n")})
When to Reach for Which Tool
Situation
Tool
Part
Results differ across machines / runs
check float comparisons, summation order
I
Collinear regressors, huge SEs
check \(\kappa(\mathbf{X})\), standardise, QR/SVD
II
Estimator defined by an equation
bracket + Brent; Newton with checked derivative
III
Estimator defined by an optimum
BFGS with analytic gradient; multiple starts
IV
Need derivatives / integrals of black-box \(f\)
numDeriv / quadrature; MC in high dimension
V
Continuous-time model
RK4 / lsoda — never hand-rolled Euler in production
VI
Thousands of independent replications
mclapply / joblib on 6 cores, proper RNG streams
VII
Dense linear algebra dominates; \(n\) huge
torch on GPU, float32 explore / float64 verify
VIII
One sentence per part: floats are finite (I) — so factor, never invert (II) — roots and optima are iterations built on solves (III, IV) — derivatives and integrals trade truncation against rounding (V) — dynamics compound one step’s error into a path (VI) — and when one core is not enough, replicate across cores (VII) or vectorise onto thousands of them (VIII).
Exercises — Errors, Matrices & Equations
Machine epsilon by hand — write a loop that starts at \(h = 1\) and halves until \(1 + h == 1\). Compare the result with .Machine$double.eps / np.finfo(float).eps. Repeat in float32 (np.float32 / torch_float()).
Quadratic roots stably — for \(x^2 + 10^8 x + 1 = 0\), compute both roots with the textbook formula, then with the stable variant \(x_1 = q/a\), \(x_2 = c/q\) where \(q = -\tfrac{1}{2}(b + \text{sign}(b)\sqrt{b^2 - 4ac})\). Which root loses all digits and why?
Kahan summation — implement compensated summation and compare against naive sum() on \(10^7\) values of 0.1 in float32. How many digits does each retain?
Longley by hand — standardise the Longley regressors (mean 0, sd 1) and recompute \(\kappa(\mathbf{X})\). How many digits does standardisation buy back? Does the normal-equations route now agree with lm()?
Ridge as stabiliser — add \(\lambda\mathbf{I}\) to \(\mathbf{X}^\top\mathbf{X}\) for the Longley data with \(\lambda \in \{10^{-8}, 10^{-4}, 1\}\) and track \(\kappa\) and the coefficients. Connect numerical stabilisation to shrinkage bias.
Multiple IRRs — the cash flow \((-1000, 3600, -4310, 1716)\) has three sign changes. Plot NPV\((r)\) on \([0, 1]\), find all roots with repeated bracketed uniroot / brentq calls, and explain which (if any) is economically meaningful.
Cournot equilibrium — two firms with costs \(c_i q_i^2\) face inverse demand \(P = 100 - Q\). Write the two best-response first-order conditions and solve the system with nleqslv / scipy.optimize.root for \((c_1, c_2) = (1, 2)\). Verify second-order conditions.
Exercises — Optimization, Dynamics & Acceleration
Probit from scratch — repeat the Part IV Newton–Raphson for a probit on the same CSV (gradient and Hessian involve the normal pdf/cdf ratio — mind the tails: use log(pnorm()) via pnorm(log.p=TRUE) / scipy.special.log_ndtr). Compare with glm(family = binomial(link = "probit")) / sm.Probit.
Bad starting values — start the logit Newton iteration at \(\boldsymbol\beta_0 = (20, 20, 20, 20)\). What happens and why? Add step-halving (halve \(\mathbf{d}\) until the log-likelihood improves) and show it repairs convergence.
Gradient check discipline — introduce a deliberate bug (drop the minus sign on one gradient component) and show that BFGS still reports “convergence”. What do the estimates look like? What does check_grad / numDeriv report?
Step-size U-curve for the Hessian — compute the logit Hessian by finite differences with \(h \in \{10^{-1}, \ldots, 10^{-12}\}\) and plot the error against the analytic Hessian. Where is the optimum \(h\)? Compare with Part V’s theory.
Ramsey in two equations — solve the system \(\dot{k} = k^{0.36} - c - 0.1k\), \(\dot{c} = c\,(0.36 k^{-0.64} - 0.15)\) with deSolve / solve_ivp from several initial \(c(0)\) and visualise the saddle path instability.
Stiff test — solve \(\dot{y} = -1000(y - \cos t)\) with hand-rolled Euler at \(h = 0.01\) and with lsoda / LSODA. Explain the explosion and the fix.
Parallel bootstrap of a regression — bootstrap the logit coefficient of x1 (\(B = 1000\)) serially and with 6 workers. Report the speedup and the implied Amdahl fraction. Verify both give the same CI with proper RNG streams.
GPU break-even — time torch matmul for \(n \in \{100, 500, 1000, 2000, 4000\}\) on CPU and GPU (remember to synchronise). Plot both curves: at what \(n\) does the GPU overtake, and why not sooner?
float32 vs float64 MLE — rerun the torch logit in float32. How many digits of the coefficients survive? Which part of the computation loses them?
Further Reading
The numerical-analysis canon:
Goldberg (1991). What every computer scientist should know about floating-point arithmetic. ACM Computing Surveys 23(1). doi:10.1145/103162.103163
Anderson et al. (1999). LAPACK Users’ Guide, 3rd ed. SIAM — the library underneath every solve(), qr() and chol() in Part II. doi:10.1137/1.9780898719604
IEEE (2019). Standard for Floating-Point Arithmetic (IEEE 754-2019) — the actual specification behind Part I. doi:10.1109/IEEESTD.2019.8766229
The algorithms this deck implements:
Nelder & Mead (1965). A simplex method for function minimization. Computer Journal 7(4), 308–313. doi:10.1093/comjnl/7.4.308
Broyden (1970). The convergence of a class of double-rank minimization algorithms. IMA J. Applied Mathematics 6(1), 76–90 — the B in BFGS. doi:10.1093/imamat/6.1.76
Dennis & Schnabel (1996). Numerical Methods for Unconstrained Optimization and Nonlinear Equations. SIAM — the reference for the Newton–Raphson and quasi-Newton slides. doi:10.1137/1.9781611971200
Virtanen et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods 17, 261–272. doi:10.1038/s41592-019-0686-2
Thank You
Athanassios Stavrakoudis Applied Informatics and Computational Economics Lab Department of Economics University of Ioannina, Greece