Computing OLS When the Design Matrix Is Too Large for Memory

In Ordinary Least Squares, you need the matrix product $X^T X$ (where $X$ is $n \times p$) to compute the coefficient estimate $\hat{\beta} = (X^T X)^{-1} X^T Y$. But suppose $n$ is so large that you cannot fit the full matrix $X$ in memory. How would you compute $X^T X$ (and $X^T Y$) without storing the entire design matrix? Describe a specific, numerically accurate method -- not an approximate algorithm like SGD.