Linear Regression from Scratch

An assistant reads retrieved context before responding. In a small replay table, latency climbs from 118 to 198 milliseconds as context grows from 100 to 500 tokens. A straight line suggests about 220 milliseconds for 600 tokens. How do we find that line, when does an analytical matrix inversion beat iterative steps, and what could a training score of 0.9965 be hiding?

The linear-algebra lesson introduced least-squares projections, and Training & Backpropagation demonstrated iterative parameter updates. Here we'll connect both worlds: formulate the hypothesis function and Mean Squared Error loss, solve the weights with the closed-form Normal Equation, examine its computational complexity against gradient descent, inspect the geometry of ill-conditioned loss bowls, compare Ridge and Lasso regularization, and evaluate the core Gauss-Markov assumptions that tell us when ordinary least squares can be trusted.

Predicting latency from retrieved evidence

For each request replay, the input feature is retrieved text length measured in hundreds of tokens. The target is end-to-end latency in milliseconds. Keeping units visible keeps physical intuition intact: x = 1 means 100 tokens, while y = 118 means 118 ms. These five rows represent synthetic teaching measurements:

Observed latency climbs by roughly 20 ms for each additional 100 tokens. Before setting up a formal solver, we can express our hypothesis with two parameters:

• The intercept, b, is the line's prediction at x = 0 (baseline latency overhead).
• The slope, w (or m), is the change in predicted latency for each additional 100 tokens.

For a single observation with feature vector $x \in \mathbb{R}^D$, the linear regression hypothesis function is:

$$\hat{y} = w^\top x + b$$

With one feature, this simplifies to $\hat{y} = w x + b$. A sensible first candidate line is:

$$\hat{y} = 100 + 20x$$

Read this as an operational baseline: start around 100 ms of fixed overhead, then add 20 ms for each 100 tokens of retrieved context. For x = 6 (600 tokens), the candidate line predicts $100 + 20 \times 6 = 220$ ms.

This is single-feature linear regression. We choose b and w to minimize squared prediction errors across our training dataset, a criterion called ordinary least squares (OLS). The model captures an empirical association; it doesn't prove that retrieved text alone causes this latency response.

An intercept can approximate fixed overhead, but we haven't measured x = 0. Calling the fitted 100 ms a measured server boot cost would extrapolate beyond our evidence. Likewise, predicting latency at x = 6 extends past the observed range, making it an extrapolation.

Before running the inspection script, predict five requests, an evidence range of 1 to 5, and a latency range of 118 to 198 ms.

Residuals show what the line missed

At x = 1, the candidate line predicts 120 ms while the replay measured 118. At x = 3, it predicts 160 ms while the replay measured 162. A model's predictions become informative when you name and inspect that signed gap. The residual is:

$$e_i = y_i - \hat{y}_i$$

On the raw latency scale, small 1-2 ms errors almost disappear against a 100-200 ms span. Isolating residuals on their own dedicated scale reveals both their sign and their structural pattern.

Using the candidate line $\hat{y} = 100 + 20x$, write out each prediction before combining errors into an aggregate metric:

Signs carry direct diagnostic information:

• Negative residual: the model predicted too high (the request ran faster than predicted).
• Positive residual: the model predicted too low (the request ran slower than predicted).
• Residual near zero: the observation sat directly on the fitted line.

Residual signs show how predictions miss, but the sequence of signs exposes model structure. The residuals run $[-2, +1, +2, +1, -2]$: negative at both endpoints and positive in the center. That low-high-low arc indicates that the true latency curve bends, even though every individual miss is small.

To summarize performance in a single number, square every residual and compute the mean:

$$\text{MSE} = \frac{1}{N} \sum_{i=1}^N (y_i - \hat{y}_i)^2 = \frac{4 + 1 + 4 + 1 + 4}{5} = 2.8\,\text{ms}^2$$

That quantity is mean squared error. Squaring stops positive and negative misses from canceling each other out, while penalizing large blunders much more heavily than small discrepancies.

MSE is measured in squared units: milliseconds squared. Taking its square root gives the root mean squared error: $\text{RMSE} = \sqrt{2.8} \approx 1.67$ milliseconds. That returns the error scale back to the original units of observed latency.

Before interpreting a model's error, compare it with the simplest plausible fallback: predicting the target mean $\bar{y} = 160$ ms for every request. That mean baseline produces an MSE of 802.8, whereas the candidate line produces 2.8. This comparison verifies that the evidence feature supplies real predictive signal beyond target averaging.

The coefficient of determination, R-squared, measures the proportion of target variance explained by the model:^{[1]Reference 1r2_scorehttps://scikit-learn.org/stable/modules/generated/sklearn.metrics.r2_score.html}

$$R^2 = 1 - \frac{\sum (y_i - \hat{y}_i)^2}{\sum (y_i - \bar{y})^2} = 1 - \frac{14}{4014} = 1 - \frac{2.8}{802.8} \approx 0.9965$$

An $R^2$ of 0 indicates performance matching the constant mean baseline. An $R^2$ of 1 indicates zero residual error. Negative values indicate predictions worse than the target mean.

Dividing the numerator and denominator by $N = 5$ reproduces the ratio 1 - MSE_line / MSE_baseline. Here, $0.9965$ describes training fit quality across these five points. It doesn't make the underlying residual curve vanish, nor does it guarantee accuracy on future requests outside this sample.

Why square residuals before averaging them instead of taking their raw sum?

Answer

Taking raw sums lets positive and negative residuals cancel out, making large over-predictions and under-predictions look like zero error. Squaring ensures all errors contribute positively, penalizing large blunders disproportionately.

Before running the evaluation cell, predict residuals [-2, 1, 2, 1, -2], MSE 2.8, baseline MSE 802.8, and $R^2 \approx 0.9965$.

Why the math becomes matrix multiplication

For a single feature, writing $\hat{y} = 100 + 20x$ is straightforward. When adding cache hit ratios, prompt complexity, and system load, code needs a vectorized representation that scales cleanly.

We construct a design matrix $X \in \mathbb{R}^{N \times (D+1)}$: each row represents one observation, and each column supplies one model term. An all-ones leading column supplies the intercept term, while subsequent columns hold feature values:

$$X = \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \\ 1 & 4 \\ 1 & 5 \end{bmatrix}$$

Packing the parameters into a single vector $w = [b, m]^\top \in \mathbb{R}^2$, the prediction rule across all $N$ requests becomes a matrix-vector product:

$$\hat{y} = X w$$

For this replay table, the optimal weights are $w = [100, 20]^\top$. Vectorization doesn't alter the math; it lets linear algebra libraries execute dot products simultaneously across rows and features.

Check the third request ($x = 3$), where the design row is $[1, 3]$:

$$\begin{bmatrix} 1 & 3 \end{bmatrix} \begin{bmatrix} 100 \\ 20 \end{bmatrix} = 1 \times 100 + 3 \times 20 = 160$$

The leading 1 multiplies the intercept 100, while 3 multiplies the slope 20. Matrix multiplication repeats that calculation across every row simultaneously.

Geometrically, the vector of predictions $\hat{y} = X w$ is a linear combination of $X$'s columns. At the ordinary least squares minimum, $X w^*$ is the orthogonal projection of target vector $y$ onto the column space $\text{col}(X)$. The resulting residual vector $e = y - X w^*$ is orthogonal to every column of $X$:^{[2]Reference 2Numerical Linear Algebra.https://people.maths.ox.ac.uk/trefethen/books.html}

$$X^\top (y - X w^*) = 0$$

Arbitrary candidate weights don't satisfy this orthogonality condition. We'll verify it in code after solving for the weights.

What happens if you remove the all-ones column from the design matrix?

Answer

The model loses its intercept parameter b and can only predict w * x, forcing the regression line through the origin (0, 0). If the true process has a non-zero baseline, this constraint severely degrades fit quality.

Before running the comparison script, predict the result of dropping the intercept: the line will be forced through zero and its MSE will climb past 802.8.

Forcing the line through zero yields an MSE of 1820.98, more than double the mean baseline's error of 802.8. Allowing a fitted offset is essential for these data, even when the offset doesn't isolate an isolated physical source of delay.

The Normal Equation: closed-form ordinary least squares

How do we solve for $w^*$ without guessing combinations of $b$ and $m$? We express the total Mean Squared Error loss in matrix notation and apply matrix calculus:

$$J(w) = \frac{1}{N} \|X w - y\|_2^2 = \frac{1}{N} (X w - y)^\top (X w - y)$$

Expanding the inner product gives:

$$J(w) = \frac{1}{N} \left( w^\top X^\top X w - 2 y^\top X w + y^\top y \right)$$

Taking the gradient with respect to the parameter vector $w$:

$$\nabla_w J(w) = \frac{2}{N} \left( X^\top X w - X^\top y \right)$$

Some textbooks define the loss with a leading $\frac{1}{2N}$ factor to cancel the 2 in the derivative. Both formulations share identical stationary points. Setting the gradient to zero produces the Normal Equation:

$$X^\top X w = X^\top y$$

When $X$ has full column rank, the Gram matrix $X^\top X \in \mathbb{R}^{(D+1) \times (D+1)}$ is symmetric positive definite and invertible. Solving yields the closed-form ordinary least squares estimator:^{[3]Reference 3The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

$$w^* = (X^\top X)^{-1} X^\top y$$

For our latency replay, the normal system is:

$$\begin{bmatrix} 5 & 15 \\ 15 & 55 \end{bmatrix} \begin{bmatrix} b \\ m \end{bmatrix} = \begin{bmatrix} 800 \\ 2600 \end{bmatrix}$$

Row-reducing gives $55m - 3(15m) = 10m = 2600 - 2400 = 200$, so $m = 20$. Back-substitution yields $5b + 15(20) = 800 \implies b = 100$. Because the loss is a convex quadratic, this stationary point is the global minimum.

Evaluating the computational complexity reveals important scaling limits:

1. Forming $X^\top X$ requires multiplying a $(D+1) \times N$ matrix by an $N \times (D+1)$ matrix, costing $O(N D^2)$ arithmetic operations.
2. Inverting or solving the $(D+1) \times (D+1)$ linear system via Cholesky decomposition ($X^\top X = L L^\top$) costs $O(D^3)$ operations.
3. Storing $X^\top X$ requires $O(D^2)$ memory.

Total time complexity for the Normal Equation is $O(N D^2 + D^3)$. When feature count $D$ is small (for example, $D \le 10{,}000$) and data fits in memory, the closed-form solve provides the exact optimal weights in one analytical step without tuning learning rates or step counts. When $D$ reaches hundreds of thousands (such as wide bag-of-words or embedding tables), $O(D^3)$ inversion and $O(D^2)$ memory become intractable.

Forming $X^\top X$ explicitly also squares the condition number of the system: $\kappa(X^\top X) = \kappa(X)^2$. For production code, np.linalg.lstsq solves $(X, y)$ directly using singular value decomposition (SVD) or QR factorization, avoiding explicit Gram matrix construction while handling rank deficiency cleanly.^{[2]Reference 2Numerical Linear Algebra.https://people.maths.ox.ac.uk/trefethen/books.html[4]Reference 4numpy.linalg.lstsqhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html}

Before running the closed-form script, predict weights [100, 20], MSE 2.8, and an orthogonality check of True.

The check X.T @ residuals confirms geometric projection: the residuals are orthogonal to the column span of $X$. No linear adjustment of the weights can reduce squared error further.

Next, inspect the conditioning and numerical stability of $X$ with np.linalg.lstsq. Expect rank 2, condition number cond_X near 8.37, and cond_XTX near 69.99.

The 2-norm condition number of $X$ is $\sigma_{\max} / \sigma_{\min} \approx 7.69 / 0.92 \approx 8.37$. Squaring that condition number when forming $X^\top X$ produces $69.99$. On five well-behaved rows that solves cleanly, but on ill-conditioned data that squaring can wipe out numerical precision.

What happens when feature columns contain redundant information? Create an artificially flawed feature matrix by appending double_tokens = 2 * tokens. The third column provides zero new information. Predict a singular matrix error from np.linalg.solve, rank 2 for 3 columns from lstsq, and an unchanged MSE of 2.8.

The prediction remains identical because $4x + 8(2x) = 20x$. However, the division of weight between the collinear columns isn't unique. Among the infinite set of solutions that achieve minimal MSE, np.linalg.lstsq returns the solution with the smallest Euclidean norm $\|w\|_2$.^{[4]Reference 4numpy.linalg.lstsqhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html} Here that produces $[100, 4, 8]$ with norm 100.4, compared to 101.98 for $[100, 20, 0]$. SVD picks a stable minimum-norm vector, but it can't separate the physical contributions of duplicate features.

Gradient descent finds the same line by walking downhill

When datasets contain millions of observations or thousands of features, constructing and inverting $X^\top X$ is impractical. Iterative optimization sidesteps matrix inversion by taking successive steps downhill along the negative gradient.

We walk the gradient descent update loop: initialize weights, evaluate predictions, compute derivatives, and adjust weights by a step size called the learning rate ($\alpha$):

$$w \leftarrow w - \alpha \nabla_w \text{MSE}$$

For Mean Squared Error, the parameter gradients derive from the chain rule:

$$\nabla_w \text{MSE} = \frac{2}{N} X^\top (X w - y)$$

Breaking this into separate scalar derivatives for the intercept and feature weights clarifies the arithmetic:

$$\frac{\partial \text{MSE}}{\partial b} = \frac{2}{N} \sum_{i=1}^N (\hat{y}_i - y_i), \qquad \frac{\partial \text{MSE}}{\partial w_j} = \frac{2}{N} \sum_{i=1}^N (\hat{y}_i - y_i) x_{i,j}$$

Starting with initial weights $w_0 = [0, 0]^\top$, every prediction is 0. Because zero sits below every observed latency value, each prediction error $\hat{y}_i - y_i$ is negative:

The slope gradient magnitude is over three times larger than the intercept gradient because higher evidence values pull more aggressively on the slope derivative.

With learning rate $\alpha = 0.05$, the first parameter update is:

$$w_1 = \begin{bmatrix} 0 \\ 0 \end{bmatrix} - 0.05 \begin{bmatrix} -320 \\ -1040 \end{bmatrix} = \begin{bmatrix} 16 \\ 52 \end{bmatrix}$$

The slope immediately overshoots its target of 20, landing at 52. Yet total MSE drops from 26402.8 to 2194.8. Optimization progress reflects joint loss reduction across the parameter space, not whether every individual coordinate moves monotonically toward its final value.

Each iteration of full-batch gradient descent requires computing $X w$ ($O(N D)$ operations) and $X^\top (\hat{y} - y)$ ($O(N D)$ operations). For $K$ iterations, total runtime is $O(K N D)$ while requiring only $O(D)$ additional storage for the gradient vector. No $D \times D$ matrix is ever constructed or inverted.

Before running the single-step verification, predict gradient [-320, -1040], next weights [16, 52], and an MSE drop from 26402.8 to 2194.8.

Now run 1000 steps of gradient descent and verify convergence toward the OLS solution $[100, 20]$.

After 1000 iterations, the iterative weights match the closed-form solution to within $10^{-4}$, confirming that both solvers converge to the identical least-squares minimum.

Convex loss bowl geometry and step size selection

Why did the first gradient step jump past the optimal slope of 20 to 52, and why did later iterations take dozens of steps to crawl from $[84.05, 24.42]$ to $[100, 20]$? The answer lies in the quadratic geometry of the loss surface.

The second-derivative matrix (Hessian) of the Mean Squared Error loss is constant:

$$H = \nabla_w^2 \text{MSE} = \frac{2}{N} X^\top X$$

Because $X^\top X$ is symmetric and positive definite when $X$ has full rank, all eigenvalues of $H$ are strictly positive: $\lambda_i > 0$. This confirms that the MSE loss surface is a strictly convex paraboloid bowl with a single global minimum. Local minima and saddle points don't exist.

The eigenvalues of $H$, $\lambda_{\max}$ and $\lambda_{\min}$, govern the surface curvature along its principal axes:

• The condition number of the Hessian is $\kappa = \frac{\lambda_{\max}}{\lambda_{\min}}$.
• For our replay data, the Hessian matrix is $\begin{bmatrix} 2 & 6 \\ 6 & 22 \end{bmatrix}$, yielding $\lambda_{\max} \approx 23.66$ and $\lambda_{\min} \approx 0.34$.
• The condition number is $\kappa \approx \frac{23.66}{0.34} \approx 70$.

A condition number of 70 means the loss bowl is 70 times steeper along its sharpest axis than along its flattest floor. The level contours are elongated ellipses rather than circles. Gradients point almost perpendicular to the long valley floor, forcing gradient descent to oscillate across the steep walls while creeping slowly along the base.

This geometry determines the mathematical bound for learning rate stability. For gradient descent on a quadratic surface to converge without diverging, the step size must satisfy:

$$0 < \alpha < \frac{2}{\lambda_{\max}}$$

For our replay, $\frac{2}{\lambda_{\max}} \approx \frac{2}{23.66} \approx 0.0845$. Our chosen learning rate of $0.05$ satisfies this stability bound ($0.05 \lt 0.0845$), ensuring eventual convergence. If we increase $\alpha$ to $0.2$, the update step exceeds the valley diameter, causing updates to bounce outward with exponentially growing loss.

The optimal theoretical fixed step size for pure quadratic descent is $\alpha^* = \frac{2}{\lambda_{\max} + \lambda_{\min}}$, which yields an asymptotic convergence rate of $\frac{\kappa - 1}{\kappa + 1}$. When features are standardized (zero mean, unit variance), the columns of $X$ become nearly orthogonal, reducing $\kappa \to 1$ and transforming narrow ravines into circular bowls where gradient descent reaches the minimum rapidly.

Test what happens when the stability bound is violated by setting $\alpha = 0.2$. Predict an exploding loss sequence where each step overshoots further away from the minimum.

Within five iterations, loss explodes from $2.6 \times 10^4$ to $9.4 \times 10^8$. Numerical divergence is an optimization step-size failure, not evidence that linear regression is unsuited for the problem.

Stabilizing ill-conditioned fits with Ridge and Lasso

When feature columns are strongly correlated, sample sizes are limited, or feature counts exceed observation counts ($D > N$), ordinary least squares suffers from severe variance. Small shifts in training observations produce wild swings in estimated coefficients. Regularization counters this by introducing a parameter penalty into the optimization objective.^{[3]Reference 3The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

Ridge Regression (L2 regularization) adds the sum of squared feature weights to the loss:

$$J_{\text{ridge}}(w) = \frac{1}{N} \|X w - y\|_2^2 + \lambda \|w_{1:D}\|_2^2$$

Setting the gradient to zero yields an analytical closed-form solution:

$$w^*_{\text{ridge}} = \left( X^\top X + \lambda N I \right)^{-1} X^\top y$$

Adding $\lambda N I$ shifts every eigenvalue of $X^\top X$ upward by $+\lambda N > 0$. Even when $X^\top X$ is strictly rank-deficient or singular, $X^\top X + \lambda N I$ is guaranteed to be strictly positive definite and invertible. Ridge regression smoothly shrinks coefficients toward zero, stabilizing correlated parameters and reducing prediction variance at the expense of introducing slight bias.

Lasso Regression (L1 regularization) penalizes the sum of absolute feature weights:

$$J_{\text{lasso}}(w) = \frac{1}{N} \|X w - y\|_2^2 + \lambda \|w_{1:D}\|_1 = \frac{1}{N} \|X w - y\|_2^2 + \lambda \sum_{j=1}^D |w_j|$$

Because the absolute value function has a non-differentiable corner at zero, Lasso can't be solved with standard matrix inversion. Instead, it's solved via subgradient optimization or coordinate descent using the soft-thresholding operator:

$$S_{\lambda}(z) = \text{sign}(z) \max(|z| - \lambda, 0)$$

The geometry of the L1 penalty explains its distinctive behavior. In parameter space, the L1 constraint boundary forms a diamond (cross-polytope) with sharp vertices positioned along the coordinate axes. The elliptical contours of the MSE loss bowl typically contact these sharp corners first. When an elliptical contour touches a corner, the corresponding parameter is driven exactly to zero. Lasso induces sparsity, performing automated feature selection by eliminating uninformative inputs.

Gauss-Markov foundations and residual diagnostic checks

Why has ordinary least squares remained the default baseline for continuous regression for two centuries? The theoretical foundation rests on the Gauss-Markov Theorem.^{[3]Reference 3The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

The theorem states that under five core assumptions, the OLS estimator $\hat{w}_{\text{OLS}}$ is BLUE: the Best Linear Unbiased Estimator. Among all linear unbiased estimators, OLS achieves the minimum sampling variance.

The five Gauss-Markov assumptions establish that guarantee:

1. Linearity in parameters: The true data-generating process takes the form $y = X w + \epsilon$. Features can undergo non-linear transformations (such as $x^2$ or $\log x$), but parameters must enter linearly.
2. Strict exogeneity (Zero conditional mean): $\mathbb{E}[\epsilon \mid X] = 0$. The noise carries zero expected value regardless of feature values. This assumption rules out omitted variable bias and measurement errors in $X$.
3. No perfect multicollinearity (Full rank): $\text{rank}(X) = D+1 \le N$. No feature is an exact linear combination of other features, ensuring $X^\top X$ is invertible.
4. Spherical error variance:
   • Homoscedasticity: $\text{Var}(\epsilon_i \mid X) = \sigma^2$ for all $i$. Noise variance remains constant across all feature levels. If latency variance spreads out for longer prompts, errors are heteroscedastic.
   • No autocorrelation: $\text{Cov}(\epsilon_i, \epsilon_j \mid X) = 0$ for $i \neq j$. Individual observations must be independent. If replays come from sequential cache warming, residuals correlate across rows.
5. Residual normality: $\epsilon \mid X \sim \mathcal{N}(0, \sigma^2 I)$.

A common misconception is that OLS requires features or errors to follow a normal distribution. Normality is not required for the Gauss-Markov BLUE guarantee, nor is it needed for OLS consistency. By the Central Limit Theorem, parameter estimates become asymptotically normal in large samples regardless of error distributions. Normality is required only for exact small-sample hypothesis testing ($t$-tests, $F$-tests) and finite-sample prediction intervals.

Our five-row latency replay highlights where linear assumptions begin to strain:

• The residuals $[-2, +1, +2, +1, -2]$ form an inverted parabola. This non-random structure violates the zero conditional mean assumption across local feature slices ($\mathbb{E}[\epsilon \mid x=3] \approx +2$, while $\mathbb{E}[\epsilon \mid x=1] \approx -2$).
• The macro fit looks stellar ($R^2 = 0.9965$), but the residual diagnostic confirms that a straight line misses systematic curvature.

Evaluation starts on unseen requests

Training metrics measure descriptive fit on data the solver has already inspected. They don't certify how well a model generalizes to fresh traffic. A production latency service requires an explicit train, validation, and test split containing requests from evidence lengths, prompt topics, and traffic volumes the model never saw during fitting.

To demonstrate data separation mechanics, refit using only the first four evidence lengths ($x \in \{1, 2, 3, 4\}$) and score the held-out fifth request ($x = 5$). Because we have already inspected all five points, this retrospective split provides an illustrative partition rather than untouched test evidence. In real projects, hold out test partitions before exploratory analysis.

Predict weights [98, 21], a prediction of 203 ms for x = 5, and an error of -5 ms.

The baseline predicts the training mean ($150.5$ ms) without inspecting the test observation. One scored point can't establish typical future variance or compute a meaningful test $R^2$, but it demonstrates data array separation: the training split determines parameters, while held-out data evaluates predictive accuracy.

Match a trusted implementation

Writing linear regression from scratch illuminates mathematical mechanics. Once from-scratch implementations work, verify numerical output against a maintained production library. scikit-learn's LinearRegression solves ordinary least squares, learning an intercept by default when fed a 2D feature matrix of shape (5, 1) and a target vector of shape (5,).^{[5]Reference 5LinearRegressionhttps://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html}

Predict exact numerical agreement before running the verification: intercept 100.0, slope 20.0, MSE 2.8, and $R^2 \approx 0.9965$.

Matching intercept, slope, MSE, and $R^2$ confirms that our custom NumPy implementation executes the exact mathematical operations of standard OLS libraries.

Operating contracts and extrapolation boundaries

A fitted line predicts 220 ms for x = 6, but a single scalar point estimate provides zero indication of variability. A prediction interval specifies the expected range for an individual future request, incorporating both model parameter uncertainty and observation noise variance $\sigma^2$. A confidence interval captures uncertainty around the estimated mean response. Neither transforms five training replays into a formal service-level agreement.

Deploying regression models requires explicit operating boundaries:

1. Record feature units, allowable training ranges ($1 \le x \le 5$), condition numbers, and held-out error distributions alongside saved weights.
2. If an incoming request arrives with x = 50 after training on 1 through 5, route it to an explicit fallback path rather than blindly trusting linear extrapolation. Beyond the training domain, physical dynamics often transition into different regimes (such as memory swapping or context truncation).

Common failure patterns

Low training error doesn't mean an operating model is sound. Inspect these characteristic failure patterns and diagnostic remedies:

Try it yourself

Run these focused experiments using the NumPy code above. Turn each question into a measurable output by printing updated weights, losses, and residual signatures: