Linear Regression from Scratch

An assistant reads retrieved context before responding. In our constructed five-row table, latency climbs from 118 to 198 milliseconds as context grows from 100 to 500 tokens. A straight line suggests 220 milliseconds for 600 tokens. How do we find that line, choose a numerical solver, and notice what a training score of 0.9965 hides?

The linear-algebra lesson introduced least-squares projections, and Training & Backpropagation demonstrated iterative updates. We'll solve the same fit both ways, then challenge it: a good optimizer can still deliver a model that misses the relationship or fails on future traffic.

Predicting latency from retrieved evidence

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 synthetic rows teach the mechanics. Benchmark performance with measurements from a deployed assistant:

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 (a possible approximation to baseline 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:

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

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.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])

print("requests", len(x))
print("evidence_range", (int(x.min()), int(x.max())), "hundreds_of_tokens")
print("latency_range", (int(y.min()), int(y.max())), "ms")

requests 5
evidence_range (1, 5) hundreds_of_tokens
latency_range (118, 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:

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

View larger

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.

The residuals run $[-2, +1, +2, +1, -2]$: negative at both endpoints and positive in the center. That arc suggests missing curvature. With only one observation at each length, it doesn't establish the true conditional mean or rule out random noise. In this constructed table, a quadratic happens to fit every row; we'll inspect that alternative below.

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

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. The evidence feature improves the fit on these rows; fresh data must establish whether that improvement transfers.

The coefficient of determination, R-squared, compares squared error with variation around the mean of the targets being scored:^{[1]Reference 1r2_scorehttps://scikit-learn.org/stable/modules/generated/sklearn.metrics.r2_score.html}

For nonconstant targets, $R^2 = 0$ matches that mean baseline, $R^2 = 1$ means zero error, and negative values mean worse predictions. Constant targets make the denominator zero; scikit-learn substitutes finite scores by default. With fewer than two observations, its r2_score returns NaN. On held-out data, the metric's denominator uses the held-out target mean; an actual deployable constant predictor must instead learn its mean from training data.

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.

Check your hand calculation against the next cell. Can a residual sum of zero coexist with nonzero MSE? Here it does: the signed misses cancel, but their squares sum to 14.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
pred = 100 + 20 * x
residual = y - pred
mse = np.mean(residual ** 2)
baseline_mse = np.mean((y - y.mean()) ** 2)
r2 = 1 - mse / baseline_mse

print("predictions", pred.astype(int))
print("residuals", residual.astype(int))
print("mse", round(float(mse), 1))
print("baseline_mse", round(float(baseline_mse), 1))
print("r2", round(float(r2), 4))

predictions [120 140 160 180 200]
residuals [-2  1  2  1 -2]
mse 2.8
baseline_mse 802.8
r2 0.9965

Why the math becomes matrix multiplication

For one feature, $\hat{y} = 100 + 20x$ already names the prediction. The same prediction for every row is one matrix-vector product once an intercept column sits beside those lengths.

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:

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:

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]$:

The leading 1 multiplies the intercept 100, while 3 multiplies the slope 20. Matrix multiplication expresses that calculation across every row; the library chooses how to schedule the arithmetic.

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}

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

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.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])

with_intercept = np.c_[np.ones_like(x), x]
without_intercept = x[:, None]

w_full, *_ = np.linalg.lstsq(with_intercept, y, rcond=None)
w_forced, *_ = np.linalg.lstsq(without_intercept, y, rcond=None)
mse_full = np.mean((with_intercept @ w_full - y) ** 2)
mse_forced = np.mean((without_intercept @ w_forced - y) ** 2)
baseline_mse = np.mean((y - y.mean()) ** 2)

print("with_intercept", w_full.round(2), "mse", round(float(mse_full), 2))
print("forced_through_zero", w_forced.round(2), "mse", round(float(mse_forced), 2))
print("mean_baseline_mse", round(float(baseline_mse), 2))
print("forced_worse_than_mean", mse_forced > baseline_mse)

with_intercept [100.  20.] mse 2.8
forced_through_zero [47.27] mse 1820.98
mean_baseline_mse 802.8
forced_worse_than_mean True

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 helps these data, even when the intercept doesn't isolate a physical source of delay. Some problems have a justified zero-intercept constraint; this table doesn't supply one.

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:

Expanding the inner product gives:

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

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:

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/}

For our latency replay, the normal system is:

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. Factoring the $(D+1) \times (D+1)$ positive-definite system with Cholesky ($X^\top X = L L^\top$) costs $O(D^3)$; the subsequent triangular solves cost $O(D^2)$. There is no need to explicitly compute an inverse.
3. Storing $X^\top X$ requires $O(D^2)$ memory.

These are dense-matrix costs, with the intercept absorbed into $D$ in big-O notation. Total time is $O(N D^2 + D^3)$, plus storage for inputs and solver workspace. At 10,000 features, a float64 Gram matrix alone has 100 million entries, roughly 800 MB; that is already substantial. Choose a solver based on rank, conditioning, sparsity, data size, and hardware. A direct solve avoids learning-rate tuning but still performs many finite-precision operations.

For a full-column-rank matrix, forming $X^\top X$ squares its 2-norm condition number: $\kappa_2(X^\top X) = \kappa_2(X)^2$. Direct least-squares methods based on QR or SVD avoid that construction. NumPy's lstsq returns a least-squares solution, numerical rank, and singular values, using a cutoff for small singular values.^{[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} We'll use the normal system below to expose the arithmetic; prefer lstsq for this dense fit in application code.

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

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])

X = np.c_[np.ones_like(x), x]

w_closed = np.linalg.solve(X.T @ X, X.T @ y)
pred_closed = X @ w_closed
residuals = y - pred_closed

baseline = np.full_like(y, y.mean())
line_mse = np.mean((pred_closed - y) ** 2)
baseline_mse = np.mean((baseline - y) ** 2)
r2 = 1 - line_mse / baseline_mse

print("closed_form_weights", w_closed.round(2))
print("residuals", residuals.round(2))
print("line_mse", line_mse.round(2))
print("baseline_mse", baseline_mse.round(2))
print("r2", round(r2, 4))
print("residual_orthogonal", bool(np.allclose(X.T @ residuals, 0)))

closed_form_weights [100.  20.]
residuals [-2.  1.  2.  1. -2.]
line_mse 2.8
baseline_mse 802.8
r2 0.9965
residual_orthogonal 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.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
X = np.c_[np.ones_like(x), x]

w, residual_sum, rank, singular_values = np.linalg.lstsq(X, y, rcond=None)
pred = X @ w

print("weights", w.round(2))
print("rank", rank, "columns", X.shape[1])
print("mse", round(float(np.mean((pred - y) ** 2)), 2))
print("singular_values", singular_values.round(2))
print("cond_X", round(float(np.linalg.cond(X)), 2))
print("cond_XTX", round(float(np.linalg.cond(X.T @ X)), 2))

weights [100.  20.]
rank 2 columns 2
mse 2.8
singular_values [7.69 0.92]
cond_X 8.37
cond_XTX 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.

Numerical rank is also a solver decision. In NumPy 2.0 and later, rcond=None uses machine precision times the larger matrix dimension as a relative singular-value cutoff. A singular value below that fraction of the largest one is treated as zero. Increasing the cutoff can discard a nearly dependent direction, changing both coefficients and fitted predictions; it doesn't prove exact redundancy in the underlying features.^{[4]Reference 4numpy.linalg.lstsqhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html}

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.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
X_bad = np.c_[np.ones_like(x), x, 2 * x]

try:
    np.linalg.solve(X_bad.T @ X_bad, X_bad.T @ y)
except np.linalg.LinAlgError:
    print("normal_equation", "singular")
else:
    print("normal_equation", "returned weights; rank still needs checking")

w, residual_sums, rank, _ = np.linalg.lstsq(X_bad, y, rcond=None)
print("lstsq_rank", rank, "columns", X_bad.shape[1])
print("lstsq_weights", w.round(2))
print("weight_norm", round(float(np.linalg.norm(w)), 2))
print("alt_norm_100_20_0", round(float(np.linalg.norm([100.0, 20.0, 0.0])), 2))
print("returned_residual_sums", residual_sums.tolist())
print("recomputed_sse", round(float(np.sum((X_bad @ w - y) ** 2)), 2))
print("mse", round(float(np.mean((X_bad @ w - y) ** 2)), 2))

normal_equation singular
lstsq_rank 2 columns 3
lstsq_weights [100.   4.   8.]
weight_norm 100.4
alt_norm_100_20_0 101.98
returned_residual_sums []
recomputed_sse 14.0
mse 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.

The returned residual_sums is empty even though squared errors sum to 14. NumPy leaves this array empty when the design is rank deficient or has no more rows than columns. Empty means that this diagnostic wasn't returned, not that the fit is exact. Compute y - X @ w directly when you need errors for every solver case.^{[4]Reference 4numpy.linalg.lstsqhttps://numpy.org/doc/stable/reference/generated/numpy.linalg.lstsq.html}

Gradient descent finds the same line by walking downhill

When the Gram matrix or direct factorization is too expensive, iterative methods avoid constructing it. Millions of rows with a handful of features can still be manageable for a direct solve; many sparse features can favor sparse iterative solvers. Full-batch gradient descent is one option, not an automatic winner at a fixed data size.

The direct solve formed a Gram matrix and factored it. Gradient descent never builds that matrix. It starts from a guess and steps against the loss slope, scaled by a learning rate $\alpha$:

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

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

The NumPy shapes must preserve those row-by-row errors. Our single-target loop uses X.shape == (N, D+1), w.shape == (D+1,), and y.shape == (N,). Subtracting an (N, 1) target from an (N,) prediction is valid broadcasting, but creates an (N, N) matrix comparing every target with every prediction.^{[5]Reference 5Broadcasting.https://numpy.org/doc/stable/user/basics.broadcasting.html} Predict whether the following bug raises an exception or produces a plausible-looking number:

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y_column = np.array([118.0, 141.0, 162.0, 181.0, 198.0])[:, None]
pred = 100 + 20 * x
wrong_error = pred - y_column
y = y_column[:, 0]  # One target per request, matching this loop's contract.
error = pred - y
assert error.shape == y.shape == (len(x),)

print("broadcast_error_shape", wrong_error.shape)
print("broadcast_mse", round(float(np.mean(wrong_error ** 2)), 1))
print("paired_error_shape", error.shape)
print("paired_mse", round(float(np.mean(error ** 2)), 1))

broadcast_error_shape (5, 5)
broadcast_mse 1602.8
paired_error_shape (5,)
paired_mse 2.8

There is no shape error: the wrong loss is 1602.8 rather than 2.8. Consistent column vectors are also valid if weights, predictions, targets, and gradients all use that convention. For multiple targets, preserve the output axis instead of flattening it. Check shapes before allowing an optimizer to minimize the wrong objective.

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:

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 full-batch iteration computes $X w$ and $X^\top (\hat{y} - y)$, each costing $O(N D)$ for dense data. For $K$ iterations, runtime is $O(K N D)$. The code below stores $O(D)$ weights and gradients plus $O(N)$ predictions/errors, in addition to $X$ and $y$. Processing rows in batches can reduce that temporary $O(N)$ storage, but changes the implementation; stochastic mini-batch updates also change the optimization dynamics.

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

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
X = np.c_[np.ones_like(x), x]
w = np.zeros(2)

error = X @ w - y
grad = (2 / len(X)) * X.T @ error
w_next = w - 0.05 * grad

print("gradient", grad)
print("next_weights", w_next)
print("old_mse", round(float(np.mean(error ** 2)), 1))
print("new_mse", round(float(np.mean((X @ w_next - y) ** 2)), 1))

gradient [ -320. -1040.]
next_weights [16. 52.]
old_mse 26402.8
new_mse 2194.8

After one step, the slope is farther from 20 even though loss fell. Watch the two coordinates recover together over 1,000 steps, then compare with the direct fit.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
X = np.c_[np.ones_like(x), x]

w = np.zeros(2)
lr = 0.05

for step in range(1000):
    pred = X @ w
    error = pred - y
    grad = (2 / len(X)) * (X.T @ error)
    w -= lr * grad

    if step in {0, 99, 999}:
        loss = np.mean((X @ w - y) ** 2)
        print(f"after_step={step + 1} loss={loss:.2f} w={w.round(2)}")

w_direct, *_ = np.linalg.lstsq(X, y, rcond=None)
print("matches_direct", bool(np.allclose(w, w_direct, atol=1e-4, rtol=0)))

after_step=1 loss=2194.80 w=[16. 52.]
after_step=100 loss=49.09 w=[84.05 24.42]
after_step=1000 loss=2.80 w=[100.  20.]
matches_direct True

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:

For full column rank, $X^\top X$ is positive definite and all eigenvalues of $H$ are positive. The loss is strictly convex with one minimum and no nonglobal local minima or saddle points. A rank-deficient design still gives a convex loss, but flat directions allow multiple minimizing weight vectors.

The duplicate-feature design above has one such direction: $v=[0,2,-1]^\top$ satisfies $Xv=0$. Adding $v$ changes coefficients without changing any prediction. Every gradient is a combination of the rows of $X$, so it is orthogonal to $v$: gradient descent can't change the initial component along this flat direction. With a convergent step size and zero initialization, it reaches the minimum-norm solution. A nonzero initial component in the null space survives.^{[6]Reference 6Linear Regression (Cornell CS 3/5780, Spring 2026)https://www.cs.cornell.edu/courses/cs3780/2026sp/slides/lec12.pdf}

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
X = np.c_[np.ones_like(x), x, 2 * x]
null_direction = np.array([0.0, 2.0, -1.0])
assert np.allclose(X @ null_direction, 0)

def fit_from(start):
    w = start.copy()
    for _ in range(5_000):
        w -= 0.01 * (2 / len(x)) * X.T @ (X @ w - y)
    return w

from_zero = fit_from(np.zeros(3))
from_null = fit_from(null_direction)
print("zero_start_weights", from_zero.round(2).tolist())
print("null_start_weights", from_null.round(2).tolist())
print("same_predictions", bool(np.allclose(X @ from_zero, X @ from_null)))
print("surviving_difference", (from_null - from_zero).round(2).tolist())

zero_start_weights [100.0, 4.0, 8.0]
null_start_weights [100.0, 6.0, 7.0]
same_predictions True
surviving_difference [0.0, 2.0, -1.0]

Both runs minimize the same training MSE, but their coefficients differ. Solver agreement should compare predictions and loss as well as coefficients; rank deficiency changes what coefficient agreement can mean.

View larger

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 near 70 means a curvature ratio near 70; contour semiaxes have ratio $\sqrt{70} \approx 8.37$ in these parameter coordinates. With our step size, the steep direction oscillates while the flat direction decays slowly. A sufficiently smaller step would avoid that oscillation, but still make slow progress. Conditioning and the learning rate jointly determine the path.

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:

For our table, $2/\lambda_{\max} \approx 0.0845$, so $\alpha = 0.05$ converges in exact arithmetic. In each eigenvector direction, parameter error is multiplied by $1-\alpha\lambda_i$ per step. At $\alpha=0.2$, the steep direction's multiplier has magnitude greater than 1, producing exponential growth when that direction has nonzero initial error.

The fixed step $\alpha^* = 2/(\lambda_{\max}+\lambda_{\min})$ minimizes the worst-case parameter-error contraction factor, $(\kappa-1)/(\kappa+1)$, for this positive-definite quadratic. The corresponding worst-case loss-gap factor is its square. Neural networks need learning-rate choices tested against their own optimization dynamics.

Centering and scaling our single nonconstant feature to unit population variance makes it orthogonal to the intercept column and gives $H=2I$, so $\kappa=1$. With several features, standardization preserves their correlations. Two duplicate standardized columns remain duplicates.^{[7]Reference 7StandardScalerhttps://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.StandardScaler.html}

Try to predict both results before running: the one-feature Hessian becomes isotropic, but the duplicate-feature design stays rank deficient.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
z = (x - x.mean()) / x.std(ddof=0)
X_one = np.c_[np.ones_like(z), z]
H_one = (2 / len(x)) * X_one.T @ X_one

features = np.c_[x, 2 * x]
scaled = (features - features.mean(axis=0)) / features.std(axis=0, ddof=0)
X_duplicate = np.c_[np.ones_like(x), scaled]

# Display tiny roundoff as zero; retain the original Hessian for computations.
H_display = np.where(np.abs(H_one) < 1e-12, 0.0, H_one)
print("single_feature_H", H_display.round(2).tolist())
print("single_feature_condition", round(float(np.linalg.cond(H_one)), 2))
print("scaled_duplicates_equal", bool(np.allclose(scaled[:, 0], scaled[:, 1])))
print("duplicate_design_rank", np.linalg.matrix_rank(X_duplicate), "columns", X_duplicate.shape[1])

single_feature_H [[2.0, 0.0], [0.0, 2.0]]
single_feature_condition 1.0
scaled_duplicates_equal True
duplicate_design_rank 2 columns 3

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.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
X = np.c_[np.ones_like(x), x]
w = np.zeros(2)
losses = []

for _ in range(5):
    error = X @ w - y
    losses.append(float(np.mean(error ** 2)))
    grad = (2 / len(X)) * X.T @ error
    w -= 0.2 * grad

print("losses", [round(loss, 1) for loss in losses])
print("loss_increased", losses[-1] > losses[0])

losses [26402.8, 349474.8, 4852473.8, 67584328.5, 941482657.6]
loss_increased True

The five recorded losses, from initialization through the fourth update, grow from $2.6 \times 10^4$ to $9.4 \times 10^8$. The loop performs a fifth update after its last recording. This example isolates a step-size failure; it doesn't tell us whether a straight line is the right statistical model.

Stabilizing ill-conditioned fits with Ridge and Lasso

Strongly correlated columns can make coefficients sensitive to small data changes. If there are more parameters than rows, the OLS weights aren't unique without an additional choice such as minimum norm. Neither situation alone proves that prediction variance is severe. Regularization adds a parameter penalty and can improve stability; its predictive value needs held-out evaluation.^{[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:

Our objective excludes the intercept $w_0=b$. Define $R=\operatorname{diag}(0,1,\ldots,1)$, then set the gradient to zero:

Use a linear solve rather than explicitly computing that inverse. With at least one row, an all-ones intercept column, all other coefficients penalized, and $\lambda>0$, the system is positive definite even when feature columns are redundant. For any vector $v$, its quadratic form is $\|Xv\|^2+\lambda N\sum_{j=1}^D v_j^2$. It can be zero only when all feature coordinates vanish and the intercept does too.

The familiar $X^\top X+\lambda N I$ formula penalizes every coordinate, including the intercept. That different objective does shift every eigenvalue by the same amount; the feature-only penalty generally doesn't. Ridge shrinks feature weights overall and introduces bias in their estimates. Every correlated coefficient needn't shrink monotonically, and better future prediction isn't guaranteed.

What changes if we accidentally replace $R$ with $I$? In this table, using $\lambda=1$ changes both the baseline and slope. The feature-only fit preserves the training target mean because its intercept remains free.

import numpy as np
from sklearn.linear_model import Ridge

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
X = np.c_[np.ones_like(x), x]
lam = 1.0
R = np.diag([0.0, 1.0])
w_feature = np.linalg.solve(X.T @ X + lam * len(x) * R, X.T @ y)
w_all = np.linalg.solve(X.T @ X + lam * len(x) * np.eye(2), X.T @ y)

# sklearn Ridge uses summed squared error + alpha * squared feature norm.
model = Ridge(alpha=lam * len(x)).fit(x[:, None], y)
w_library = np.r_[model.intercept_, model.coef_]
assert np.allclose(w_feature, w_library)

print("feature_only_weights", w_feature.round(2).tolist())
print("all_coordinates_weights", w_all.round(2).tolist())
print("feature_only_matches_sklearn", bool(np.allclose(w_feature, w_library)))
print("feature_only_mean_matches_target", bool(np.isclose((X @ w_feature).mean(), y.mean())))
print("all_coordinates_mean_matches_target", bool(np.isclose((X @ w_all).mean(), y.mean())))

feature_only_weights [120.0, 13.33]
all_coordinates_weights [24.0, 37.33]
feature_only_matches_sklearn True
feature_only_mean_matches_target True
all_coordinates_mean_matches_target False

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

Lasso also leaves the intercept unpenalized here. Its absolute-value corner requires different optimization machinery, such as coordinate descent or a proximal method. The soft-thresholding operator is:

The threshold depends on the loss convention and feature norm. For our MSE-plus-penalty objective, a coordinate update for a nonzero feature column is

The factor $\lambda/2$ comes from differentiating MSE; using half-MSE would change it. The intercept uses an unpenalized update. In two penalized dimensions, the L1 constraint is a diamond whose corners have zero coordinates. This geometry favors sparse solutions, but contact can occur along an edge too. Lasso doesn't identify which features are truly uninformative: correlated inputs can be interchangeable, selected features may be spurious, and weak useful features may be omitted.

Fit scaling statistics on training data and reuse them on validation/test data. Both penalties depend on feature units; choose strength with a validation procedure that matches future traffic.

Library parameters need the same care. scikit-learn's Ridge uses summed squared error, so our objective maps to alpha = lambda * N. Its Lasso uses half-MSE, so our Lasso objective maps to alpha = lambda / 2. Identically named strength parameters aren't necessarily on the same scale.^{[8]Reference 8Linear Modelshttps://scikit-learn.org/stable/modules/linear_model.html}

Gauss-Markov foundations and residual diagnostic checks

The Gauss-Markov theorem gives a precise, limited reason to value OLS. Conditional on a full-column-rank design and the assumptions below, it is BLUE, the Best Linear Unbiased Estimator. "Linear" describes an estimator that is a linear function of the observed targets $y$; "best" means its covariance is no larger than that of another linear unbiased estimator.^{[9]Reference 9Notes on the Gauss-Markov theoremhttps://www.stat.berkeley.edu/~census/GaussMar.pdf}

The conditions are:

1. Linearity in parameters: $y=Xw+\epsilon$. Columns may include $x^2$ or $\log x$; the coefficients still enter linearly.
2. Zero conditional error mean: $\mathbb{E}[\epsilon\mid X]=0$. Omitted effects correlated with included features or measurement error in features can violate this condition; an omitted variable alone needn't do so.
3. Full column rank: $\operatorname{rank}(X)=D+1\le N$, so the weights are identifiable.
4. Constant conditional error variance: $\operatorname{Var}(\epsilon_i\mid X)=\sigma^2<\infty$.
5. Uncorrelated conditional errors: $\operatorname{Cov}(\epsilon_i,\epsilon_j\mid X)=0$ for distinct rows. Uncorrelated doesn't require independent. Sequential traffic may violate this condition through shared load or cache state.

Together, the last two conditions say $\operatorname{Cov}(\epsilon\mid X)=\sigma^2I$. No normality assumption is part of the BLUE guarantee. Gaussian errors give the familiar exact small-sample $t$/$F$ results under the rest of the classical model; other inference methods use other assumptions. OLS consistency and asymptotic normality need sampling and moment conditions as well as identification. A larger sample alone doesn't overcome dependence, endogeneity, or heavy tails that break those conditions.

Don't substitute fitted residuals $e=y-X\hat w$ for unobservable population errors $\epsilon$. OLS residuals share constraints such as $X^\top e=0$ and can be correlated even when the true errors aren't. Five residuals aren't five estimates of a conditional expectation.

For our table, adding $x^2$ gives coefficients $[93,26,-1]$ and zero training error. The model is curved in $x$ but still linear in its parameters. Compare its extrapolation with the straight line before deciding which would generalize.

import numpy as np

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])
X_line = np.c_[np.ones_like(x), x]
X_curve = np.c_[np.ones_like(x), x, x ** 2]
w_line, *_ = np.linalg.lstsq(X_line, y, rcond=None)
w_curve, *_ = np.linalg.lstsq(X_curve, y, rcond=None)

print("quadratic_weights", w_curve.round(2).tolist())
print("quadratic_training_fit_exact", bool(np.allclose(X_curve @ w_curve, y)))
print("line_x6_ms", round(float(np.array([1.0, 6.0]) @ w_line), 1))
print("quadratic_x6_ms", round(float(np.array([1.0, 6.0, 36.0]) @ w_curve), 1))

quadratic_weights [93.0, 26.0, -1.0]
quadratic_training_fit_exact True
line_x6_ms 220.0
quadratic_x6_ms 213.0

Both predictions are outside the observed range. Exact training fit doesn't establish a physical law. Collect fresh data, compare held-out errors, and inspect whether feature transforms remain sensible over the intended operating range.

Evaluation starts on unseen requests

Training metrics describe data used by the solver. To estimate performance on fresh traffic, plan a train, validation, and test split that matches intended use. A random split may suit independent requests; a time split tests future traffic, and a grouped split can keep repeated users or replay families together. Test novel length/topic/load ranges separately when those shifts matter. Neither random splitting nor deliberately unfamiliar inputs automatically represent deployment.

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.

Will removing the final row make the slope steeper or flatter? The first four increments are larger; the fit becomes $98+21x$ and predicts 203 ms for the held-out x=5.

import numpy as np

x_train = np.array([1.0, 2.0, 3.0, 4.0])
y_train = np.array([118.0, 141.0, 162.0, 181.0])
x_test = np.array([5.0])
y_test = np.array([198.0])

X_train = np.c_[np.ones_like(x_train), x_train]
X_test = np.c_[np.ones_like(x_test), x_test]
w, *_ = np.linalg.lstsq(X_train, y_train, rcond=None)
prediction = X_test @ w
baseline = np.full_like(y_test, y_train.mean())

print("train_weights", w.round(2))
print("test_prediction_ms", prediction.round(2))
print("test_error_ms", (y_test - prediction).round(2))
print("line_squared_error", round(float(np.mean((prediction - y_test) ** 2)), 2))
print("baseline_squared_error", round(float(np.mean((baseline - y_test) ** 2)), 2))

train_weights [98. 21.]
test_prediction_ms [203.]
test_error_ms [-5.]
line_squared_error 25.0
baseline_squared_error 2256.25

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

The NumPy fit is in hand. The check that remains is whether a maintained solver prints the same intercept, slope, and errors on these five rows. 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,).^{[10]Reference 10LinearRegressionhttps://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$.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

x = np.array([1.0, 2.0, 3.0, 4.0, 5.0]).reshape(-1, 1)
y = np.array([118.0, 141.0, 162.0, 181.0, 198.0])

model = LinearRegression().fit(x, y)
pred = model.predict(x)

print("intercept", round(float(model.intercept_), 2))
print("slope", np.round(model.coef_, 2))
print("mse", round(float(mean_squared_error(y, pred)), 2))
print("r2", round(float(r2_score(y, pred)), 4))

intercept 100.0
slope [20.]
mse 2.8
r2 0.9965

Agreement supports the correctness of our calculation on this fixture. It doesn't prove every edge case or imply identical numerical algorithms: a library can center data and use a least-squares factorization instead of solving a Gram system. Compare rank/tolerances and ill-conditioned cases when those matter.

As checked on October 10, 2026, scikit-learn 1.9's LinearRegression uses tol as the singular-value cutoff cond for dense least squares; earlier releases only applied it to the sparse iterative solver. This differs from NumPy's lstsq(rcond=None) default cutoff. The well-conditioned fixture agrees, but a near-collinear design may give different numerical ranks unless you align these choices. Record library versions and tolerances when comparing solvers.

Operating contracts and extrapolation boundaries

A fitted line predicts 220 ms for x = 6, but the point estimate doesn't express uncertainty. A classical regression prediction interval includes both uncertainty in the fitted mean and variation of an individual future observation. A confidence interval for the mean excludes that extra individual-observation variation. State a coverage level and the method's assumptions; nominal coverage can fail under misspecification or distribution shift. Neither interval turns five constructed rows into a service-level agreement.

Even a careful interval doesn't tell the service what to do with a length it never saw. Write that boundary down with the weights:

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

These symptoms suggest checks, not unique causes. Use evidence to choose the next experiment:

Try it yourself

The five latencies show a stable step, a diverging one, and a bend that $R^2$ hides. Print the weights, losses, and residuals for each change below: