# 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
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.

• This question appears to be off-topic because it is only about software. Morever, the syntax is incorrect--the arguments to lm have been reversed--and the problem cannot be reproduced (this machine completed the task in 18 seconds using less than 1GB RAM). – whuber Nov 13 '14 at 18:52
• Do you need all the bells & whistles of, say, lm(), or do you only need point estimates from the fit? Can you load all of the data into your software, or do you need to load & removed parts of your data at a time? – gung - Reinstate Monica Nov 13 '14 at 19:50
• @whuber ok so the computing solution is a sparse matrix package for the design matrix, the $(X^T X)^{-1} X^T Y$ should calculate ok... Still feel like there might be another statistical way to view it though, but may be not... – Corone Nov 13 '14 at 21:32
• Yes, that's precisely what I calculated. As you point out, $X$ itself contains only about ten million nonzero entries and so requires only a few hundred million bytes of storage. The calculation of $X^\prime X$ can occur in-place, creating a matrix with less than two million entries. You don't actually want to find $(X^\prime X)^{-1}$, so you just solve the system $X^\prime X \beta = X^\prime Y$, which is what took most of the two seconds. – whuber Nov 13 '14 at 21:36
• No, that's the hard way! It takes extensive computational effort to invert a $p\times p$ matrix: it amounts to solving $p$ separate systems of linear equations, whereas you need to solve only one such system. In your case that would be $1300$ times as much work as you really need to do :-). Another nice aspect of this is that you can solve these normal equations even when $X^\prime X$ is not invertible--but in such cases there are entire linear families of solutions (all of which are equally good fits). A good linear equation solver won't be bothered by such complications at all. – whuber Nov 13 '14 at 21:57

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.