Linear Regression from Scratch: Solvers, Geometry and Evidence

Why this session exists

Linear regression is a particularly useful place to learn optimization because a stable direct solver gives a reference for an iterative implementation. The same residual–gradient–update pattern recurs in other differentiable models. It is not true that every ML algorithm is a generalization of linear regression, or that every subsequent model lacks an analytic solution.

We will keep four questions separate: what predictions the features permit, what objective selects a fit, how the numerical algorithm finds it, and whether that fitted model is useful on unseen data. A solver comparison on one dataset answers the third question, not the fourth.

Prerequisites

Know the shapes of a matrix product, the transpose, and partial derivatives. Review the ML experiment contract for splits and leakage. The CPU checks use small local arrays. The scikit-learn continuation is also executed in the existing CPU environment, with matched intercept conventions.

(a) Intuition · 5 min

Imagine one vertical spring connecting each observed point to a candidate line. Squared stretch is an energy: larger residuals contribute disproportionately. The least-squares fit minimizes the total energy. The analogy describes vertical residuals, not shortest perpendicular distances from points to the line; those define a different fitting problem.

Let A be the design matrix, including a column of ones for the intercept, and let θ contain the intercept and feature weights. With n rows and p columns, A has shape (n, p), θ (p,), and y (n,). Predictions are $A\theta$.

The three ideas you'll own by the end

1. Model: a weighted sum of chosen features, including an optional intercept.
2. Objective: $L(\theta)=\|A\theta-y\|_2^2/(2n)$. This is half mean squared error, not MSE itself. The factor changes gradients and step-size conventions but not the minimizing θ.
3. Fit: solve the least-squares problem numerically, or iterate $\theta\leftarrow\theta-\eta\nabla L$. Compare prediction residuals and stationarity, not just rounded coefficients.

The gradient — one derivative, done once

Expand the square:

Since $A^TA$ is symmetric,

Setting the gradient to zero gives the normal equations $A^TA\theta=A^Ty$. If A has full column rank, the formal expression is $(A^TA)^{-1}A^Ty$. That expression is a derivation, not an instruction to compute an inverse. Use np.linalg.lstsq(A, y, rcond=None) for the reference solution.

A hand calculation makes the signs visible. Take rows (1, 0) and (1, 1), targets (1, 3), and initial θ (0, 0). Residuals are (−1, −3). The gradient is (−2, −1.5). With step 0.1, new θ is (0.2, 0.15), predictions (0.2, 0.35), and half-MSE falls from 2.5 to 1.915625. A gradient of the opposite sign increases the error near the starting point.

From Gauss to modern optimization

The history is useful because it separates three questions: choosing an error criterion, studying conditional association, and finding an optimizer. In Gauss's Theoria Motus (1809; the captured English translation is1857), the preface describes the1801 discovery/observations of Ceres and the orbit-estimation problem. The later least-squares discussion says he had used the principle since1795 and acknowledges Legendre's publication. This is Gauss's retrospective claim, not independent proof of priority; it cannot mean he predicted Ceres in1795 before the1801 observations. The captured OCR prints a bibliographic year as1806 in that passage, so it is not used to establish Legendre's first-publication date. The mathematical lesson is the relationship between measurement-error assumptions and a squared-residual criterion, not an invented teenage-Ceres anecdote.

Galton's1886 Regression towards Mediocrity in Hereditary Stature examines parent/offspring stature and the tendency for conditional offspring averages to be closer to the population mean. This describes conditional averages, not a law that every individual moves toward a mean or that correlation proves causation. His anthropometric framing is historical context, not a normative claim about people. Modern least squares uses the projection machinery far beyond this application.

Robbins and Monro's1951 A Stochastic Approximation Method studies finding a root of an unknown mean response from noisy experiments. To connect it to this lesson, set that mean response to an expected gradient: a noisy update can seek its zero without computing every observation's gradient each time. Convergence needs assumptions on noise, regularity and step sizes; it is not a guarantee for every neural network. BFGS-style quasi-Newton optimization instead builds an inverse-Hessian approximation from successive gradient/step information, trading memory and curvature estimates against plain gradient descent. It is not the same algorithm as the closed-form normal equation, and a single1974 invention date for the whole family was unjustified.

Kingma and Ba's Adam preprint first appeared inDecember2014, not2010. Adam tracks first and second gradient moments with bias correction; it is an adaptive stochastic optimizer, not immunity from poor learning rates or proof of universal convergence. A direct least-squares solver remains preferable when its problem structure and cost fit. Learning GD here is valuable precisely because a stable reference solution makes optimizer bugs observable.

Primary reading: Gauss1857 translation of the1809 work, Galton1886, Robbins–Monro1951, SciPy1.16.2 optimizer methods, Adam v1.

The loss surface — a convex bowl

For any vector v, $v^THv=\|Av\|^2/n\geq0$, so the objective is convex. Full column rank makes H positive definite and coefficients unique. Rank deficiency creates flat directions: there can be many coefficient vectors with identical training predictions. The fitted projection $A\theta$ is unique even when θ is not.

If $Av=0$, then $A(\theta+v)=A\theta$. For duplicate feature columns, increasing one coefficient and decreasing the other equally does not change fitted predictions. This is non-identifiability, not proof that the training data cannot be fitted. Numerical rank decisions use a tolerance, so nearly collinear columns can also cause large coefficient sensitivity.

Projection onto the column space

The columns of A span all prediction vectors the model can express on these rows. Least squares finds the point in that subspace closest to y. At the optimum, $A^T(A\theta-y)=0$: the residual is orthogonal to every column. With an intercept, residuals also sum to zero, up to numerical error.

Adding a feature can enlarge that subspace but cannot worsen the minimum training squared error on the same rows. It may leave it unchanged if the feature is redundant, and it may worsen validation error. Ordinary training $R^2$ therefore cannot by itself justify adding features. Orthogonality on the training sample also does not imply independent, homoscedastic or Gaussian errors.

Two ways to fit — same answer, different cost

There is no universal feature-count cutoff such as “10,000 is cheap.” Sparse structure, memory, hardware and accuracy requirements matter. Gradient descent is not the only answer to a singular matrix: lstsq solves rank-deficient least squares directly. Ridge also retains a linear-system solution. Lasso, by contrast, needs methods that handle its nonsmooth penalty.

NumPy's documented lstsq returns the smallest 2-norm solution when multiple minimizers exist. The scikit-learn 1.7 LinearRegression notes name SciPy lstsq for ordinary dense least squares and nnls for nonnegative fitting; its sparse path uses lsqr. It is incorrect to say every case simply inverts the Gram matrix or calls NumPy's lstsq.

Learning-rate zoo — three failure modes

Write coefficient error along a Hessian eigenvector with eigenvalue λ as e. One GD step maps it to $(1-\eta\lambda)e$. For positive curvature directions, contraction requires $0<\eta<2/\lambda_{\max}(H)$. At the upper boundary the largest-curvature component can alternate without shrinking; above it that component grows. For a positive semidefinite H, nullspace components do not move.

• Too small: the factors are close to one, so error decays slowly. Confirm a correct gradient before increasing η.
• Too large: high-curvature components oscillate or grow. Lower η, inspect scales, and check for invalid inputs; do not wait for NaNs as the only failure signal.
• Plateau: a small gradient may mean convergence, but poor conditioning can make progress slow in low-curvature directions. A model-capacity limitation can also keep predictive error high after optimization succeeds.

Feature standardization, fitted only on training rows, often improves conditioning. It does not eliminate correlations or guarantee a perfectly spherical Hessian. There is no universal safe rate interval and no rule that η=1 must diverge. The code constructs an eigenvector counterexample rather than inventing loss traces.

Linear in parameters does not mean straight in raw inputs

A model $\beta_0+\beta_1x+\beta_2x^2$ is nonlinear in x but linear in β. Once basis columns are fixed, it is still least squares. Interactions, sine/cosine terms and splines can also supply columns. Increasing basis complexity expands the function family, but extrapolation and variance can deteriorate. Fit any learned basis transforms within the training boundary and compare on held-out inputs.

Why squared error, and what it does not assume

With independent Gaussian errors of constant variance σ², the negative log-likelihood is a constant plus $\sum_i(y_i-a_i^T\theta)^2/(2\sigma^2)$. Thus least squares is maximum likelihood under that model. But least squares remains a well-defined projection objective without Gaussian errors. Non-Gaussian noise does not automatically make it useless; squared prediction loss targets a conditional mean when it exists.

Laplace errors motivate absolute-error fitting; Huber loss blends quadratic and linear penalties. Heavy tails do not necessarily mean a Laplace distribution, and a central-limit argument does not establish Gaussian residuals in a particular dataset. Choose a loss with the decision costs and observed residual behaviour in mind. Even absolute loss has analytic solutions in special cases: an intercept-only absolute-error fit is a median. Smooth convex functions do not all have closed-form minimizers.

(c) Hands-on · 25 min

The following is one self-contained CPU experiment. It checks the hand step, every gradient coordinate, direct versus GD fits, a mini-batch optimization gap, rank deficiency, curvature-based divergence, polynomial features and ridge scaling. This is a solver test on synthetic data, not a reported generalization benchmark.

import numpy as np
rng = np.random.default_rng(85)
X = rng.normal(size=(160, 4))
A = np.column_stack([np.ones(len(X)), X])
y = A @ np.array([2., 1., -2., 0.5, 3.]) + rng.normal(0, 0.1, len(X))
 
def loss(B, target, theta):
    return 0.5 * np.mean((B @ theta - target) ** 2)
 
def gradient(B, target, theta):
    return B.T @ (B @ theta - target) / len(target)
 
def gd(B, target, rate, steps=800):
    theta = np.zeros(B.shape[1])
    history = [loss(B, target, theta)]
    for _ in range(steps):
        theta -= rate * gradient(B, target, theta)
        history.append(loss(B, target, theta))
    return theta, np.array(history)
 
# Hand step and finite differences.
hand_A = np.array([[1., 0.], [1., 1.]])
hand_y = np.array([1., 3.])
hand_theta = -0.1 * gradient(hand_A, hand_y, np.zeros(2))
assert np.allclose(hand_theta, [0.2, 0.15])
assert np.isclose(loss(hand_A, hand_y, hand_theta), 1.915625)
theta = rng.normal(size=A.shape[1])
numeric = np.zeros_like(theta)
for j in range(len(theta)):
    step = np.zeros_like(theta); step[j] = 1e-5
    numeric[j] = (loss(A, y, theta + step) - loss(A, y, theta - step)) / 2e-5
assert np.allclose(numeric, gradient(A, y, theta), atol=1e-8)
 
reference, _, rank, singular = np.linalg.lstsq(A, y, rcond=None)
H = A.T @ A / len(y)
eigenvalues, eigenvectors = np.linalg.eigh(H)
L = eigenvalues[-1]
fitted, history = gd(A, y, rate=1 / L)
assert np.allclose(fitted, reference, atol=1e-9)
assert np.all(np.diff(history) <= 1e-12)
assert np.linalg.norm(A.T @ (A @ reference - y)) < 1e-9
# Normal equations only as a small well-conditioned comparison, no inverse.
normal = np.linalg.solve(A.T @ A, A.T @ y)
assert np.allclose(normal, reference)
 
# Mini-batch SGD: seeded order, decaying step, check objective not rounded equality.
sgd = np.zeros(A.shape[1])
for epoch in range(160):
    order = rng.permutation(len(y))
    rate = 0.1 / (1 + epoch / 10)
    for start in range(0, len(y), 16):
        batch = order[start:start + 16]
        sgd -= rate * gradient(A[batch], y[batch], sgd)
gap = loss(A, y, sgd) - loss(A, y, reference)
assert 0 <= gap < 1e-4
 
# Duplicate one feature. Coefficients differ, fitted projection stays the same.
D = np.column_stack([A, A[:, 1]])
dup = np.linalg.lstsq(D, y, rcond=None)[0]
null = np.zeros(D.shape[1]); null[1] = 1; null[-1] = -1
assert np.allclose(D @ dup, D @ (dup + 100 * null))
assert np.allclose(D @ dup, A @ reference)
 
# Deliberately excite the largest-curvature direction above its stability bound.
error = eigenvectors[:, -1]
next_error = error - (2.1 / L) * (H @ error)
assert np.isclose(np.linalg.norm(next_error), 1.1)
# Change units of one feature; Hessian curvature grows and safe step shrinks.
wide = A.copy(); wide[:, 1] *= 1e6
assert np.linalg.eigvalsh(wide.T @ wide / len(y))[-1] > L * 1e10
 
# A fixed polynomial basis fits an exactly quadratic target.
x = np.linspace(-2, 2, 30)
quad = np.column_stack([np.ones(len(x)), x, x*x])
qy = 1 - 3*x + 2*x*x
qfit = np.linalg.lstsq(quad, qy, rcond=None)[0]
assert np.allclose(qfit, [1, -3, 2])
 
# Objective: half mean squared error + (lam/2)*sum(nonbias weights squared).
lam = 0.2
P = np.eye(A.shape[1]); P[0, 0] = 0
ridge = np.linalg.solve(A.T @ A + len(y)*lam*P, A.T @ y)
ridge_gradient = gradient(A, y, ridge) + lam * P @ ridge
assert np.linalg.norm(ridge_gradient) < 1e-10
print('rank:', rank, 'condition:', singular[0] / singular[-1])
print('reference half-MSE:', loss(A, y, reference), 'SGD gap:', gap)
print('hand step, gradients, solvers, rank, rates, basis and ridge checks passed')

What each block does

gradient divides by the number of rows actually supplied, so the short final mini-batch is normalized correctly. gd records the loss after each update and includes the initial value. A finite-difference check compares independent changes in the objective against the analytic derivative, catching sign, transpose and factor-of-two bugs.

The reference and batch-GD coefficients agree here because this seeded design has full rank and the step size and iteration budget are sufficient. The mini-batch estimate is stochastic and receives an objective-gap tolerance rather than a fictitious universal three-decimal guarantee. In the rank-deficient case, compare predictions and residual orthogonality; the nullspace test proves why coefficient equality is unnecessary.

The curvature test uses the error recurrence itself, so it does not run an overflowing training loop. The unit-change test makes a feature a million times larger and observes a greatly increased largest eigenvalue. A fixed rate chosen for the original units is not automatically safe for that representation.

Optional library comparison

After running the preceding block, this continuation uses scikit-learn1.7.2 in the existing CPU environment. It deliberately keeps intercept conventions aligned and compares predictions with a tolerance.

from sklearn.linear_model import LinearRegression
sk = LinearRegression().fit(X, y)
assert np.allclose(sk.predict(X), A @ reference, rtol=1e-7, atol=1e-8)
print('intercept:', sk.intercept_, 'coefficients:', sk.coef_)

This comparison is diagnostic only. If it fails, check preprocessing, weights, intercept handling, rank tolerance and convergence before assuming a library bug. For a useful predictive model, fit on training rows and measure validation/test error separately.

(d) Production reality · 15 min

Separate fitting, prediction and statistical inference. You can compute a least-squares fit without normal errors. To interpret ordinary textbook standard errors or exact small-sample tests, assumptions about the conditional mean, variance, dependence and noise distribution matter. Heteroscedasticity-robust standard errors address a variance issue under their own conditions; they do not repair confounding, leakage, a wrong outcome definition or arbitrary dependence. A coefficient in observational data is not automatically a causal effect.

Use residual-versus-prediction plots to look for curvature and changing variance, and residual-versus-time plots to inspect dependence. A funnel pattern can indicate heteroscedasticity rather than a missing mean feature. Extremely large cancelling coefficients suggest collinearity. A great training fit alongside poor validation performance suggests examining feature complexity and split design, not simply increasing optimization iterations.

The earlier company's-war-story percentages and universal ship-the-baseline thresholds are removed. The choice between a baseline and a more complex model depends on absolute error cost, uncertainty, latency, maintainability and safety. A one-percent change may be operationally important or indistinguishable from noise; percentage improvement alone does not make the decision.

Ridge stretch: keep the normalization consistent

For

the gradient is $A^T(A\theta-y)/n+\lambda P\theta$ and the normal equations are

The factor n is essential with this definition. A summed squared-error objective and a mean squared-error objective use different numerical regularization strengths. The old lesson mixed their formulas. P leaves the intercept unpenalized. With a nonempty design containing an intercept and λ positive, penalizing all other coordinates makes this system positive definite: a vector in the penalty's nullspace has only its intercept coordinate, whose nonzero value cannot also be in A's nullspace. If other directions are unpenalized, inspect their rank instead of saying ridge is “always invertible.”

Where this shows up in the rest of the plan

Logistic regression changes the response model and loss. Regularization changes the preference among fits. Neural networks compose affine transformations with nonlinearities; stacking only affine layers still yields an affine function. Basis expansion changes the input columns while keeping this least-squares solver. A Gaussian prior can yield a ridge-like maximum-a-posteriori estimate, but the likelihood and prior scales determine the precise penalty.

(e) Recall + stretch · 10 min

1. Derive the normal equations for half-MSE and state when coefficients are unique.
2. For the two-row hand example, compute the gradient and one update without running Python.
3. Why should a library use a least-squares solver rather than explicitly invert the Gram matrix?
4. Can two different coefficient vectors both be correct? Produce a nullspace example.
5. What is the fixed-step stability bound, and why is η=1 not universally wrong?
6. Distinguish a batch, a mini-batch update and an epoch.
7. Fit a curved function using a linear-in-parameters model. What changes in the design?
8. Does Gaussian maximum likelihood mean OLS is invalid for every non-Gaussian target?
9. Derive the ridge system above, including n and the unpenalized intercept.
10. Why does decreasing training residual error not establish production value?

Explain-out-loud test

Explain projection geometry, draw a flat loss direction from duplicate columns, and diagnose a loss curve before touching the learning rate. Then explain why correct gradients, good conditioning and honest evaluation are three different checks.

Additional recall answers and production experiment

Why the column of ones? It makes the intercept one coordinate of theta. Without it, the fitted function is constrained through the origin in feature space. Do not both append ones and ask a solver to add an independent intercept without checking rank.

What assumptions matter for inference? A correctly specified conditional mean and exogeneity concern coefficient interpretation; dependence and variance assumptions concern standard errors; Gaussian errors give particular exact finite-sample distributions, not the existence of a least-squares solution. Robust covariance estimators address specified variance/dependence conditions, not omitted confounders. Normality is not a blanket prerequisite for prediction, but violating mean/split assumptions can certainly harm prediction.

What caused a singular matrix? Rank deficiency may arise from duplicate columns, an all-zero column, too few rows or other linear dependence—not only two perfectly correlated columns. SVD/lstsq can still produce a minimum-norm solution. Ridge with an unpenalized intercept is invertible under the conditions derived above, not for every arbitrary partially penalized matrix.

Feature-scaling practice solution. Multiply one feature by10^6 and recompute the largest eigenvalue of A-transpose-A/n before selecting a GD step; the safe bound shrinks dramatically. Fit scaling on training rows only. To report coefficients in original units for z=(x-mu)/scale, use beta_original=beta_scaled/scale and intercept_original=intercept_scaled-sum(beta_scaled*mu/scale). Standardizing does not remove collinearity. Validate predictions after inverse parameter transformation.

A/B-test and baseline case, explicitly hypothetical. Randomize treatment at the appropriate user/cluster level; freeze pre-treatment covariates and outcome window; compare a treatment coefficient with a design-consistent unadjusted estimate; use variance estimation consistent with randomization/dependence; and audit balance, attrition and interference. A later business reversal is not uniquely caused by heteroscedasticity. CUPED-style adjustment can reduce variance when covariates predict outcomes, but post-treatment leakage or a wrong estimand can invalidate it. Netflix-style rating baselines and linear feature baselines are useful comparators, not proof that linear models contribute80% of every production score or that a1% gain should always be rejected.

Gaussian-prior stretch. With Gaussian observation variance sigma² and independent zero-mean Gaussian nonbias coefficients of variance tau², the negative log posterior is residual-sum-squares/(2sigma²) plus nonbias squared norm/(2tau²). Dividing by n/sigma² gives this lesson's half-MSE plus lambda/2 penalty with lambda=sigma²/(n*tau²). The corresponding linear system adds sigma²/tau² * P; this recovers the earlier n-lambda relationship rather than mixing objective conventions.

Sources and execution limits

• NumPy 2.4 lstsq: rank handling, minimum-norm solution and singular-value cutoff.
• scikit-learn 1.7.2 linear models: OLS objective, collinearity and solver complexity.
• scikit-learn 1.7.2 LinearRegression: SciPy solver branches and intercept API.

The exact NumPy program and scikit-learn continuation are checked locally on CPU. The historical reconstruction above distinguishes original authors’ claims from priority proof. No GPU, external training data or predictive-performance benchmark is claimed.

Optional reading

Use these links for additional background. Version-specific behavior is documented in the primary references above.

• Optional resource 1
• Optional resource 2
• Optional resource 3
• Optional resource 4
• Optional resource 5
• Optional resource 6

Course navigation

Next: Logistic regression · Learning-plan hub.