Efficient linear regression given columns of A rather than rows of A

Question

I'm analyzing a large problem with a large $N \times M$ data matrix $A$, where $N$ is the number of observations, $M$ is the number of explanatory variables, and $N \gg M$. I'd like to perform single-pass linear regression on this data set against a scalar response variable. But the challenge is that I'm only able to load $N \times 1$ columns of $A$ at a time (and loading row-by-row is not possible for this problem).

So the question is, is there some way to compute $A^T A$ for linear regression using less than $O(NM)$ memory by loading columns of $A$ at a time?

Edit: Thanks for your interest. I'm looking for an efficient algorithm that makes a single pass through the data. I'd like to compute $A^T A$ exactly if possible.

Each column contains 200 million rows and is stored in separate compressed files that disallow partial read. Production system needs the data in this format. There are 60 thousand columns. The data is on a storage server, and it takes 10 seconds to load one column on a compute server in the same data center and a week to make a single pass through the data. Each compute server has 1 TB of memory. I'm currently working in parallel batches of 2 M rows. Finally, it is not possible to duplicate the data for analytics.

Answer 1

You can use some form of block matrix multiplication. Then you divide the matrix $A$ into blocks

$$A = \begin{bmatrix} B_{11} & B_{12} & \dots & B_{1n} \\ B_{21} & B_{22} & \dots & B_{2n}\\ \vdots & \vdots & \ddots & \vdots \\ B_{m1} & B_{m2} & \dots & B_{mn} \end{bmatrix}$$

and the product $A^TA$ can be expressed by a block matrix as well

$$A^TA = \begin{bmatrix} C_{11} & C_{12} & \dots & C_{1m} \\ C_{21} & C_{22} & \dots & C_{2m}\\ \vdots & \vdots & \ddots & \vdots \\ C_{m1} & C_{m2} & \dots & C_{mm} \end{bmatrix}$$

Where you compute each block as a sum of multiplications of the blocks $B$

$$C_{ij} = C_{ji}^T = \sum_{k=1}^n B_{ki}^TB_{kj}$$

And you need the memory to compute $B_{ki}^TB_{kj}$.

If you use a small number of blocks horizontally (a small $n$) then you need to reload the memory less often (because the blocks occur less often in a term $B_{ki}^TB_{kj}$). But possibly you want to use a large $n$ such that the blocks $B$ have larger columns, and reading large columns might be faster than reading large rows

Answer 2

Sharing my findings for posterity:

1. You can work in batches of $N_b$ instead of $N$, where $N_b$ is the maximum that will fit into memory. This uses $O(N_b M) + O(M^2)$ memory and reads each column $N / N_b$ times. This is the simplest method.

If you are willing to trade off accuracy for speed, you can also try:

1. Coordinate descent methods are well suited for iterating a data column at a time. There are diminishing returns on addition precision and you can stop the optimization after a few passes with an appropriate method. For this particular problem only about 10 epochs were needed to get a reasonable model using elastic net regression.

2. If you can turn the regression problem to a classification problem, you can use a slightly modified boosting algorithm: consider adding one column at a time, decide whether to accept or reject the addition, and if accepted, then modify all the weights. In our case, the response variable was standardized, so it was straightforward to convert the problem to sign classification and apply this method. This one only requires one pass through the data but gives the worst accuracy.

Lastly, keep in mind that this is a very peculiar situation. If you can build a separate data store optimized for reading rows, you can use simpler methods that don't require you to make this kind of tradeoff.