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

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.

• Of course: with a vector $x$ in RAM, accumulate the values of $y_i$ and $y_ix_i$ as you read in the vector $y.$ Repeat for all entries of $A^\prime A.$ Indeed, you only need $O(M^2)$ RAM if you're willing to read the columns an average of $(M+3)/2$ times each. – whuber Sep 6 '18 at 19:48
• It's not clear to me what you mean by "I'm only able to load $1 \times M$ columns of $A$"? Evidently you don't mean "I can load all $M$ columns of $A$"... And why would loading row-by-row not be possible? – jbowman Sep 7 '18 at 2:37
• Please register &/or merge your accounts (you can find information on how to do this in the My Account section of our help center), then you will be able to edit & comment on your own question. – gung - Reinstate Monica Jun 15 '19 at 11:40

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.

• (1) If you're working in batches, why do you need to read the columns any more than once? (2) Your proposals for sequential selection of columns will not work in general: they can achieve arbitrarily poor performance, as well as depending on the arbitrary ordering of the columns. (It sounds like an extremely limited form of forward stepwise regression.) – whuber Jun 15 '19 at 12:13
• This was an interesting project, where each column was compressed in a way that didn't allow partial reads. Partially reading a column for a batch required full transmission and decompression of the entire column, which was the bottleneck. The two proposal for approximation underperformed compared to the exact computation, as you mentioned, but they were the only way to process the data in a timely manner to satisfy the production requirements. With the proof-of-concept in production, we're now building a separate data store optimized for random reads to do this more effectively. – burning_question Jun 16 '19 at 13:27