Linear Regression from Scratch

Make one model-system cost measurable: how much latency does added policy evidence cost when the key-rotation assistant answers stale key K10234? Linear regression is the smallest useful training problem because every moving part remains visible: feature, target, intercept, slope, residual, loss, gradient, validation, and failure.

Fit one line by hand, reproduce it with NumPy, break it on purpose, and compare it with scikit-learn only after the mechanics are clear.^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/[2]Reference 2Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/[3]Reference 3Array programming with NumPy.https://www.nature.com/articles/s41586-020-2649-2[4]Reference 4Scikit-learn: Machine Learning in Python.https://jmlr.org/papers/v12/pedregosa11a.html}

You need two prerequisites here: comfort reading a small NumPy array, and familiarity with y = mx + b as a straight-line rule. We keep one tiny request replay throughout so arithmetic, code, and debugging habits stay connected.

Predicting latency from retrieved evidence

The access assistant must answer whether stale key K10234 qualifies under policy line P-7. You replay the same question while varying retrieved policy evidence. The input feature is evidence length in hundreds of tokens. The target is end-to-end decision latency in milliseconds.

This small dataset works for beginners because each number means something you can picture:

• The feature is retrieved evidence length.
• The target is measured decision latency.
• The intercept is the latency when x = 0, which you can read as fixed server overhead.
• The slope is how many more milliseconds you expect for each extra 100 tokens.

A reasonable starting guess is:

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

Read that sentence before reading it as math: "start around 100 ms, then add about 20 ms for each extra 100 tokens of policy evidence."

That one line already lets you answer concrete questions. If a request has x = 6, the model predicts:

$$100 + 20 \times 6 = 220$$

That's the core mental model of linear regression. There isn't a hidden trick. We're finding the straight line that makes measured mistakes as small as possible.

Residuals show what the line missed

A prediction becomes useful when you compare it with reality. The difference between the actual value and the predicted value is the residual:

$$\text{residual} = y - \hat{y}$$

Using the candidate line 100 + 20x, you get:

The signs matter:

• Negative residual: the model predicted too high.
• Positive residual: the model predicted too low.
• Residual near zero: the line was close on that row.

To turn all five misses into one score, we square them and average them:

$$\text{MSE} = \frac{4 + 1 + 4 + 1 + 4}{5} = 2.8$$

This score is mean squared error. Squaring does two jobs at once: it stops positive and negative errors from canceling out, and it punishes large misses more than small ones.

For this latency model, a positive residual means the request ran slower than predicted; a negative residual means it ran faster than predicted. If release gates care about slow surprises, inspect positive residuals directly in addition to reporting MSE.

A useful beginner habit is to compare against a naive baseline. The mean-latency baseline predicts 160 for each row, which gives MSE 802.8. The fitted line is useful because it crushes the "predict the average" baseline. That ratio has a name. The standard goodness-of-fit score for regression is the coefficient of determination, R-squared, which reports the fraction of the target's variance the model explains relative to the mean baseline:^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

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

R^2 = 0 means you tied the mean baseline; R^2 = 1 means perfect fit; a negative value means you did worse than predicting the average. The formula uses sums of squared errors. Dividing numerator and denominator by the same five rows gives the identical 1 - MSE_line / MSE_baseline calculation above. Here 0.9965 says the line explains about 99.65 percent of variation in these replay measurements, not every future production request.

Why the math becomes matrix multiplication

For one feature, $\hat{y} = 100 + 20x$ is readable. For many features, you want the same idea written in a form code can scale.

Start by turning the data into a design matrix: a tidy table where each row is one training example and each column is one input the model can use.

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

Each row is one training example. The leading column is all ones so the model can learn an intercept. The second column is the feature itself.

Now bundle the weights into one vector:

$$w = \begin{bmatrix} b \\ m \end{bmatrix}$$

The prediction rule becomes:

$$\hat{y} = Xw$$

Here, b is the intercept and m is the slope. For our latency data, the least-squares weights will turn out to be:

$$w = \begin{bmatrix} 100 \\ 20 \end{bmatrix}$$

It's still the same old line. We changed the notation so the computer can do the bookkeeping.

Look at one row of that bookkeeping by itself. For the third request, x = 3, so the design-matrix 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 selects the intercept. The 3 scales the slope. Matrix multiplication does that same dot product for each row at once.

NumPy implementation with the closed-form solve

Start with the normal equation, the closed-form condition that the least-squares optimum satisfies.^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/} It makes the derivation visible when the feature matrix is small and full rank.

Written as math, the normal equation is:

$$X^T X w = X^T y$$

Solving that system gives the same weights as the textbook formula $w = (X^T X)^{-1}X^T y$. In this small demonstration, np.linalg.solve avoids explicitly forming an inverse. Forming X^T X is useful for seeing the derivation, but it can magnify numerical conditioning problems. Later you'll use np.linalg.lstsq(X, y) directly, which also reports matrix rank and handles redundant columns without asking you to form X^T X.

Read the code in this order:

1. np.c_[np.ones_like(x), x] creates the intercept column and the feature column.
2. X.T @ X and X.T @ y are the matrix pieces of least squares.
3. np.linalg.solve(...) solves the system without explicitly computing a matrix inverse.
4. X @ w_closed turns weights back into predictions.
5. The baseline check proves the model learned a useful pattern instead of memorizing the mean target value.

For a reusable implementation, fit against X directly:

The two columns are independent, so rank equals the number of columns. Now make a bad feature table by adding double_tokens = 2 * tokens. That third column contains no new information.

The prediction remains the same because 4*x + 8*(2*x) = 20*x; the individual feature weights are no longer uniquely identifiable. Rank is the warning signal for an exact duplicate. Near-duplicate columns can be harder to catch: the solve may return weights, but small data changes can move those weights sharply. Inspect singular values or a condition number when coefficients look unstable.

When features are highly correlated but not exact copies, the practical response is usually regularization rather than hunting every near-duplicate by hand. Ridge regression (L2) shrinks correlated weights toward a shared solution; Lasso (L1) can zero some of them out. Both stabilize coefficients when X^T X is ill-conditioned. You'll meet them again in production feature pipelines; for this chapter, treat rank inspection as the diagnostic and Ridge or Lasso as the named next step when collinearity is real but incomplete.

There is also a scaling limit worth knowing. With N rows and D features, forming X^T X costs roughly O(ND^2), then solving the D x D system costs roughly O(D^3). With a handful of features that's instant. With very wide data, iterative optimization often becomes the practical path.^{[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

Gradient descent finds the same line by walking downhill

The closed-form solve is elegant, but it doesn't scale to each setting you'll care about later. Neural networks, embeddings, and large feature sets are usually trained by iterative optimization instead. Linear regression is the friendliest place to learn that loop.

Gradient descent is like walking down a mountain in thick fog. You don't see the bottom, but you can feel the slope under your feet. You take small steps in the steepest downward direction until the ground levels out.

Start with terrible weights:

$$w = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

That predicts 0 for each row. Each error is negative because each prediction is too low. The gradient points in the direction where loss grows fastest, so subtracting it moves the weights toward lower loss.

On the initial step, the average gradient works out to [-320, -1040]. The number comes from these terms when w = [0, 0]:

The slope gradient is larger because later rows have larger x values, so their errors pull harder on the slope than on the intercept.

The number 0.05 is the learning rate: it controls how large each step is. Too small and the model crawls; too large and it can overshoot the valley and bounce around. With lr = 0.05, the update becomes:

$$w = [0, 0] - 0.05 \times [-320, -1040] = [16, 52]$$

That's not the final answer. The slope overshoots. But the update is still useful because the model has learned the right direction: both the intercept and slope need to increase.

One common mistake is using wildly different feature scales, such as cache age in milliseconds alongside retrieved evidence measured in millions of tokens. The loss surface becomes stretched, so one coordinate can drive much larger updates than another. Gradient descent may zig-zag or diverge before it reaches the minimum. Normalize features before training when their ranges differ.

For mean squared error, the gradient is:

$$\nabla_w \text{MSE} = \frac{2}{n}X^T(Xw - y)$$

Where does the 2 come from? It's the derivative of the squared term. If you write out MSE as the average of $(Xw - y)^2$ and differentiate with respect to $w$, the chain rule brings down a 2 and leaves $X^T$ in front because each weight touches each prediction through the design matrix.^{[2]Reference 2Pattern Recognition and Machine Learning.https://www.microsoft.com/en-us/research/publication/pattern-recognition-machine-learning/[1]Reference 1The Elements of Statistical Learning.https://hastie.su.domains/ElemStatLearn/}

Plain English:

• Xw - y is the vector of prediction errors.
• X.T turns row-wise errors back into "how much should each weight change?"
• 2 / n averages the result across rows.

The full training loop:

The initial update looks dramatic because the model started absurdly wrong. Each progress line reports loss and weights after the same update, so you can compare snapshots without mixing states. Eventually, gradient descent reaches the same answer as the closed-form solve.

Remember this distinction:

• Least squares solve: jumps straight to the answer when the matrix is small enough and stable enough.
• Gradient descent: improves the answer step by step and scales to settings where a direct solve isn't practical.

The one-feature data also exposes a failed optimization run cleanly. Increase lr from 0.05 to 0.2; the steps overshoot and loss grows:

When to use each approach

This comparison matters because later AI systems still reuse the same question: "Should I solve this directly, or optimize it iteratively?" Linear regression is the smallest place to learn that tradeoff with numbers you can still inspect.

Evaluation starts on unseen requests

The five rows above were chosen to make mechanics visible. Their training R^2 is descriptive, not release evidence. A real latency model needs a train, validation, and test split with held-out replays from evidence lengths, policy documents, and load conditions the fit never saw.

For a first clean check, fit the first four evidence lengths and reserve the fifth request as a tiny test set:

One held-out row isn't an uncertainty estimate, but it teaches the correct boundary: fit on training rows, judge on data held back from fitting. A later evaluation chapter will make that split stronger with cross-validation and leakage controls.

Match a trusted implementation

From-scratch code teaches mechanics. Once it works, compare it with a maintained library on the same features and target. LinearRegression fits ordinary least squares, while mean_squared_error and r2_score apply the same metrics used above.^{[4]Reference 4Scikit-learn: Machine Learning in Python.https://jmlr.org/papers/v12/pedregosa11a.html}

Common failure patterns

Low error doesn't mean the model is healthy. Beginners usually get tricked by one of these patterns:

Our one-feature example hides one practical issue: feature scaling. With one feature, the optimizer has a simple job. Once one column is in milliseconds and another is in millions of tokens, gradient descent can zig-zag badly until you rescale the inputs. The same problem appears in AI systems when you mix latency in milliseconds with token counts, queue depth, or embedding dimensions without normalization.

Two misconceptions are worth separating:

• Linear regression doesn't require the raw feature x or target y to be normally distributed. Normal residuals matter for certain statistical inference claims, not for fitting the least-squares line itself.
• A straight line is most trustworthy inside the range you trained on. Predicting x = 6 after training on 1 through 5 is a small extension. Predicting x = 50 would be extrapolation, and the model may be confidently wrong.

Try it yourself

Start with the NumPy code above, then use the library comparison as a check rather than a substitute for understanding the fit.

If a practice answer feels hand-wavy, make it measurable. Print the new weights, print the new loss, or sketch the new residual pattern. Beginners improve faster when each guess turns into an observable result.

The fitted line is not release evidence. Next, logistic regression turns numeric latency into a binary review decision and adds probabilities, log loss, threshold tuning, ROC, and calibration.