Why does a closed-form least-squares solve scale cubically in features but linearly in rows?

## The two costs in an exact least-squares fit The closed-form fit for a linear model with `n` rows and `p` features runs in two stages, and they scale completely differently. **Stage one: summarise the data.** The fit only needs the cross-products of the columns: the `p x p` matrix `X'X` (every feature dotted with every other feature) and the `p`-vector `X'y` (every feature dotted with the target). Building `X'X` means, for each row, adding that row's `p x p` outer product into an accumulator — about `n * p^2` multiply-adds. Rows enter this term **linearly**: doubling the rows doubles this stage and nothing else. It is also a single streaming pass, so the rows never have to be in memory at once. **Stage two: solve the system.** With `X'X` and `X'y` in hand, the weights come from solving the `p x p` linear system `(X'X) w = X'y`. Solving a dense `p x p` system by a factorisation costs on the order of `p^3` operations, and the matrix itself occupies `p^2` numbers. The number of rows does not appear anywhere in this stage. Total: roughly `n*p^2 + p^3`, plus `p^2` memory. That single expression explains every practical decision on this topic. ## Why this makes "big data" the wrong question "Ten million rows" says almost nothing about whether the exact solve is feasible. What matters is the width. **Wide-in-rows, narrow-in-features: solve exactly.** A marketing-mix model with 3 million weekly rows and 8 spend channels fits a 9 x 9 system once the intercept is included. Stage one is 3,000,000 * 81 accumulations — a fraction of a second of arithmetic and one pass over the file. Stage two is a 9 x 9 solve, which is instantaneous. Iterating here would be pointless: you would take many passes over the same 3 million rows to approach, approximately, the answer one pass already gave you exactly. **Wide-in-features: the exact solve falls apart.** A store-by-SKU one-hot design with 60,000 columns has `p^2 = 3.6 * 10^9` entries in `X'X`, which at 8 bytes each is roughly 29 GB before you have solved anything, and `p^3` is on the order of `10^14` operations. The memory wall arrives first: you cannot even hold the system, let alone factor it. This is the situation the exact solve does not survive, and it is the reason to fit iteratively instead — an iterative method never forms `X'X` at all, it only ever needs to multiply by `X` and by `X'`, and it exploits the fact that a one-hot design is mostly zeros. So the switch to an iterative fit is triggered by **`p`**, and secondarily by the sparsity of the design, far more often than by `n`. ## A detail worth getting right: nobody inverts the matrix The textbook formula is written `w = (X'X)^-1 X'y`, and it is common to hear candidates say the method "inverts a matrix". A competent implementation does not. It solves the system by **factorisation** — a Cholesky decomposition of `X'X`, or a QR decomposition of `X` directly — because: - computing an explicit inverse and then multiplying costs more arithmetic than solving once; - forming `X'X` at all **squares the condition number** of the problem, so a design that was merely awkward becomes numerically hopeless. Factorising `X` directly via QR (about `2*n*p^2`) avoids that squaring; - an explicit inverse of a nearly-singular matrix is full of enormous entries that amplify floating-point noise into the answer. The inverse in the formula is notation for "solve this system", not an instruction. ## The intermediate regime Between "9 x 9" and "60,000 columns" there is a wide band where the exact solve is perfectly fine. A few thousand features means a `p^2` matrix of a few tens of millions of entries and a `p^3` solve of a few times `10^10` operations: seconds, on a laptop, however many rows you have — because the rows only cost you the streaming pass. Practitioners routinely underestimate this band and reach for an iterative fit far earlier than they need to, giving up an exact answer for an approximate one. ## What to say when asked to size a fit Give the two terms and then the numbers. State that rows are a one-pass, linear cost that streams; that features drive a `p^2` memory footprint and a `p^3` solve; that the practical failure is usually memory rather than time; and that the answer therefore depends on `p` and on whether the design is dense or sparse. Then name the crossover you would actually check: can I hold a `p x p` matrix of doubles in RAM? If yes, solve exactly. If no, iterate.