When should you use gradient descent over the normal equation to fit a linear regression?
Use gradient descent when the feature matrix is too wide, sparse, or continuously arriving for a direct least-squares solve to fit comfortably in memory, and use mini-batch or stochastic updates when you need online learning. For a modest, fixed, well-conditioned dense dataset, a direct solver is usually simpler and faster; there is no universal feature-count cutoff.
How to think about it
Use gradient descent when the feature matrix is too wide, sparse, or continuously arriving for a direct least-squares solve to fit comfortably in memory. For a modest, fixed, dense dataset, I would usually choose a direct solver instead, because it is simpler and often faster; there is no universal cutoff such as p = 10,000.
Why the choice exists
Let X be the design matrix, with n training examples and p features. Let y be the target vector and β the coefficients. Ordinary least squares minimises the squared-error objective:
J(β) = (1/n) ||Xβ - y||²
Taking the derivative and setting it to zero gives the normal equations:
XᵀXβ = Xᵀy
If XᵀX is invertible, the familiar formula is:
β = (XᵀX)⁻¹Xᵀy
That formula hides the bill.
For dense data, constructing XᵀX costs O(np²), because every example contributes to roughly p² feature pairs. Solving the resulting p by p system costs about O(p³), and storing the matrix costs O(p²) memory.
The cubic term is not the only concern. A matrix with p = 100,000 has 10¹⁰ entries. At eight bytes per floating-point entry, XᵀX alone needs about 80 GB before temporary workspaces or other model data. A p³ calculation is on the order of 10¹⁵ arithmetic operations. A closed-form formula can still have a very non-closed-form memory bill.
There is also a numerical issue. In exact mathematics, the normal equation is fine when the columns of X are linearly independent. In floating-point arithmetic, forming XᵀX squares the condition number: roughly, κ(XᵀX) = κ(X)². The condition number describes how much small numerical errors can be amplified. Highly correlated features can therefore make the normal-equation route less accurate.
For that reason, production libraries often solve least squares with QR factorisation or singular value decomposition rather than explicitly computing the inverse or even forming XᵀX. When people say “use the normal equation,” they often mean “use a direct least-squares solve.” That distinction is a useful senior-level clarification.
What gradient descent changes
Gradient descent never constructs the p by p Gram matrix. It repeatedly moves the coefficients in the direction that reduces the loss:
βₜ₊₁ = βₜ - α (2/n) Xᵀ(Xβₜ - y)
Here, α is the learning rate, which controls the step size.
A full-batch update touches every training example, so one update costs O(np) for dense data. If the model needs 50 passes over the data, the rough cost is O(50np). That can be much cheaper than O(np² + p³) when p is large, although the number of passes is not known in advance.
There are three common update patterns:
| Method | Data per update | Main trade-off |
|---|---|---|
| Batch gradient descent | All n examples | Stable updates, but every step scans the dataset |
| Stochastic gradient descent | One example | Cheap updates and online learning, but noisy |
| Mini-batch gradient descent | A batch of b examples | Usually the practical compromise |
A mini-batch update costs about O(bp). If one epoch processes every example once, the total work per epoch is still roughly O(np). The advantage is that the model only needs a manageable batch in memory, and sparse matrix operations can avoid touching zero features.
Linear regression with squared loss is convex. Its Hessian is 2XᵀX/n, which is positive semidefinite. Therefore, full-batch gradient descent has no bad local minima to get trapped in. With a suitable fixed learning rate, commonly expressed as 0 < α < 2/L, where L is the largest eigenvalue of the Hessian, it converges to a global minimiser.
But “convex” does not mean “any learning rate works.”
Common misconception: gradient descent is guaranteed to find the global optimum only under appropriate step sizes and enough iterations. If the learning rate is too large, the loss can oscillate or explode. If it is too small, training may appear frozen. If X is rank-deficient, the optimum may not be unique, even though every optimum has the same minimum loss.
The condition number matters here too. Strongly correlated features create directions in which the loss is very steep and directions in which it is almost flat. Gradient descent then zigzags through the steep direction while barely moving through the flat one. Scaling features to comparable ranges usually helps, but scaling cannot remove genuine correlation.
A concrete example
Suppose a delivery company predicts delivery time from 50,000 historical orders. The first version has 20 numeric features: distance, package weight, traffic estimates, weather measurements, and so on.
For this model, forming XᵀX involves roughly:
50,000 × 20² = 20,000,000
feature-pair contributions. The resulting matrix has only 400 entries. A direct least-squares solve is easy, and it will probably beat gradient descent that needs dozens of passes to reach the same answer.
The implementation might look like this for a small, dense problem:
import numpy as np
# Normal-equation calculation for a small, dense, full-rank problem.
beta_direct = np.linalg.solve(X.T @ X, X.T @ y)
# Full-batch gradient descent on the same squared-error objective.
beta = np.zeros(X.shape[1])
n = X.shape[0]
for epoch in range(epochs):
residual = X @ beta - y
beta -= alpha * (2.0 / n) * (X.T @ residual)The first line avoids explicitly calling a matrix-inverse function, but it still forms X.T @ X. For production code, I would normally use a library least-squares routine based on QR or SVD when numerical stability matters.
Now expand the delivery model. It includes one-hot encoded store IDs, product categories, postcode segments, and thousands of other categorical values. The model has p = 200,000 sparse features and n = 5,000,000 orders. Each order contains, on average, 30 nonzero feature values.
A dense XᵀX would contain 40 billion entries and require about 320 GB in float64. Its cubic solve would involve roughly 8 × 10¹⁵ operations. Even if the original X is sparse, feature co-occurrences can make XᵀX much less sparse, and sparse factorisations can suffer substantial fill-in.
A mini-batch method avoids that matrix. A batch of 4,096 orders contains about 4,096 × 30 = 122,880 nonzero feature visits. One pass through all five million orders involves about 150 million such visits. Five passes involve about 750 million. The actual runtime depends on hardware and implementation, but the shape of the decision is clear: sparse gradient updates are plausible; a dense normal-equation matrix is not.
The nuance that earns the senior signal
Feature count alone is not the decision rule.
If p = 20 but n = 100,000,000, even the linear scan through the data may dominate. A direct method can still be attractive because XᵀX is tiny and can be accumulated in chunks without keeping every row in memory. If exact least-squares accuracy matters, a chunked direct solve may beat SGD, which may need many passes and still has stochastic noise.
Conversely, p = 500 is not automatically small if the matrix is badly conditioned, highly sparse, or part of a latency-sensitive pipeline. The representation and the required accuracy matter as much as the dimensions.
A useful decision pattern is:
| Situation | Usual first choice | Reason |
|---|---|---|
| Small, dense, fixed dataset | QR or SVD direct solve | Accurate and no learning-rate tuning |
| Wide sparse feature matrix | Mini-batch gradient descent or SGD | Works with X without forming a dense Gram matrix |
| Continuously arriving data | SGD | Updates the model one example or batch at a time |
| Very high numerical accuracy | QR, SVD, or a suitable iterative least-squares solver | More reliable on ill-conditioned data |
p is greater than n | Ridge or another regularised iterative method | The unregularised normal matrix is singular |
Regularisation changes the objective, not the basic scaling problem. Ridge regression is commonly written as:
(XᵀX + λI)β = Xᵀy
where λ is greater than zero and I is the identity matrix. The exact factor multiplying λ depends on how the loss is normalised. Adding the ridge term makes the penalised coefficient block positive definite even when XᵀX is singular, but it does not remove the O(p²) storage or O(p³) direct-solve cost. Gradient descent can apply the same penalty through an additional term in the gradient.
One more practical detail: standard implementations usually do not penalise the intercept. That means the regularisation matrix is not literally λI across every parameter, even though the simplified formula is written that way.
A failure mode to recognise
The first symptom of a bad gradient-descent setup is often a training loss that jumps upward, alternates wildly, or becomes inf or nan after a few updates. The usual causes are an excessive learning rate, features with dramatically different scales, or overflow in the feature values.
I would inspect the loss after every epoch, check that inputs and targets are finite, standardise numeric features using training-set statistics, and reduce the learning rate. Sparse matrices should not be centred blindly, because subtracting a mean can turn a sparse matrix into a dense one.
The opposite symptom is a loss that decreases but remains far above a direct-solver baseline. That usually means too few epochs, a learning rate that is too small, or poor conditioning from correlated features. With SGD, also check whether the model is using L2 regularisation or a learning-rate schedule, since those settings mean it may no longer be optimising exactly the same objective as ordinary least squares.
What they’ll ask next
Do you explicitly compute the inverse in production?
Usually no. I would solve XᵀXβ = Xᵀy with a linear-system solver, or preferably use QR or SVD for better numerical stability. Explicitly computing an inverse does extra work and can amplify floating-point error.
Does gradient descent always reach the same answer as the normal equation?
Full-batch gradient descent reaches a global least-squares minimiser under suitable learning-rate conditions and enough iterations. It may not reach it exactly in finite time. SGD with a fixed learning rate generally keeps fluctuating around the optimum, and regularisation or early stopping intentionally produces a different solution. If the design matrix is rank-deficient, several coefficient vectors may be equally optimal.
What would you do when the number of features is greater than the number of examples?
The unregularised XᵀX matrix is singular, so the ordinary normal-equation inverse does not exist. I would use ridge regression, a sparse iterative solver, or gradient-based optimisation. The choice would depend on sparsity, the required accuracy, and whether the data arrives in batches or as a stream.
Say this in the interview
Use gradient descent when p² memory, cubic direct-solve cost, sparsity, or streaming data makes the normal equation impractical; for small, fixed, well-conditioned dense data, use a stable direct least-squares solver instead.