Linear Regression from Scratch

The previous chapter traced an LLM feature from training data to a release gate. Now make one part measurable: how much latency does added policy evidence cost when the returns assistant answers order A10234? Linear regression is the smallest useful training problem because every moving part remains visible: feature, target, intercept, slope, residual, loss, gradient, validation, and failure.

You will 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][2][3][4]

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 returns assistant must answer whether damaged order A10234 is eligible 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 gives you 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]

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

In this chapter, 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.

Here is 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] 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 will 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.

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 is instant. With very wide data, iterative optimization often becomes the practical path.^[1]

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]. Here is where that number comes from 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 distance in miles alongside package weight in grams. 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][1]

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.

Here is 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. By the end, 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 robust 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]

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 logistics when you mix delivery distance in miles with package count or weight in grams without normalization.

Two interview traps 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.

What you have built and where it leads

You can now do six useful things without a library hiding the mechanics:

1. Turn a small dataset into a design matrix.
2. Fit a line and explain what each weight means.
3. Measure error with residuals and mean squared error.
4. Compare a direct least-squares solve with gradient descent.
5. Detect a redundant-feature failure through rank.
6. Hold out a request and verify the fit against scikit-learn.

That's enough to fit a model. It isn't enough to trust one.

The next step is to move from numeric latency to a binary return decision: "Should this request be routed to human review?" Logistic regression gives you probabilities, log-loss, threshold tuning, ROC, and calibration, tools later classifiers reuse.

Mastery check

Key concepts

• intercept and slope
• design matrix
• residuals
• mean squared error
• closed-form least squares
• gradient descent
• learning rate
• baseline comparison
• R-squared
• rank deficiency
• held-out evaluation

Evaluation rubric

• Foundational: Explains intercept, slope, residuals, and mean squared error with a concrete worked example before introducing matrix notation
• Intermediate: Implements linear regression in NumPy with direct least squares and a gradient descent loop, then compares against a mean baseline and a held-out request
• Advanced: Diagnoses rank deficiency, non-linearity, unstable learning rates, outliers, and training-row evaluation as distinct failure modes with concrete fixes

Follow-up questions

Common pitfalls

• Symptom: fitted line looks consistently shifted even though slope seems plausible. Cause: the design matrix forgot the all-ones intercept column, so the model is forced through zero. Fix: rebuild X with the intercept column and refit before changing anything else.
• Symptom: adding a second token-count feature makes the direct solve fail or produces hard-to-explain weights. Cause: the new column repeats information already in x. Fix: inspect rank, remove the redundant column, and use np.linalg.lstsq while diagnosing the data design.
• Symptom: training MSE looks low, but predictions still fail on new rows. Cause: training loss was treated as enough evidence, even though residual pattern or held-out performance was never checked. Fix: compare against a baseline, inspect residuals, and evaluate on rows the model did not train on.
• Symptom: gradient descent loss bounces upward or turns unstable. Cause: learning rate is too high, or feature scales are badly mismatched. Fix: lower lr, normalize features, and watch whether the loss curve starts descending smoothly.