How to fit OLS with many categorical levels, on more than one category This question is not meant to be a software question, but I will illustrate the issue using R a bit.
My Understanding of the Simple Case
If I have a simple linear model with a categorical variable on $k$ levels, my understanding of the approach is to replace that categorical variable with $(k-1)$ dummy variables.
So if $A_i$ is the categorical variable for subject $i$ that can take three levels then one way of doing the dummy variables is:
$Y=\beta_1 X_1 + \beta_2 X_2 + c + \epsilon$
Where $X_1$ and $X_2$ are $1$ or $0$ to flag if $A_i$ is in the second or third level (and both are $0$ to indicate the first level).
This means the design matrix is $3 \times N$.
Of course, with just one category, I can just as well split the data up by category and average each set independently, thus avoiding the linear algebra completely.  This clearly scales as big as I want without problems, since I can even scan through my data if $N$ is truly gigantic.
The Problem
The problem comes if I have two (or more) categorical variables, each with a large number of possible levels.
Consider the following R code:
library(data.table)
N = 4000000
C1 = 2000
C2 = 500
DT = data.table(A = factor(sample(C1, N, rep = T)), B = factor(sample(C2, N, rep=T)), Y = rnorm(N))
modelFit <- lm(Y~ A+ B, data = DT)

This toy example is a 4m row dataset with two categorical variables, one of 2000 levels and the other of 500.  That should result in a total of 2500 parameters, but fitted with 4m datapoints, so this should be sufficient data.
The problem is by the method I described above this needs a $2500 \times 4000000$ design matrix.  Even at 1 byte per cell, that matrix would take 9.3gb!
It seems that lm in R hits a problem like this, since when I run the above code I get:
Error: cannot allocate vector of size 74.5 Gb
In addition: Warning messages:
1: In model.matrix.default(mt, mf, contrasts) :
  Reached total allocation of 8087Mb: see help(memory.size)

The Question
So, is there an alternative approach, like in the single category case, that can get the result without this huge design matrix?  Some sort of clever averaging procedure?
It seems like such a simple model having only categorical variables.
 A: I would consider the following two approaches: dimensionality reduction and sparse regression. The traditional methods for the dimensionality reduction are projective and manifold, as mentioned in my answer here and the reference within, as well as latent variables modeling (LVM) methods, such as exploratory factor analysis (EFA) and confirmatory factor analysis (CFA).
Recently I've run across two R packages - ClustOfVar and clere -  that essentially implement dimensionality reduction, but differently than traditional methods like PCA, and, in my view, are closer in their nature to the LVM approach. This alternative approach is referred to as variable clustering and you can find more details in the packages' JSS vignettes: this paper and this paper (unpublished yet?), correspondingly.
Additionally, since high-dimensional data is frequently sparse, the second main approach is to use sparse regression. R ecosystem offers many packages, useful for sparse regression analysis, such as Matrix, SprseM, MatrixModels, glmnet and flare. For links and more relevant resources, please see my related answer on DS SE site: https://datascience.stackexchange.com/a/918/2452.
For some overview and more specific examples of support for categorical explanatory variables by PCA, MDS and MCA methods as well as latent variable modeling approaches, please see this paper, these presentation slides (starting from slide 15), this paper, this paper and this paper.
