# Obtain global linear regression estimate from subsamples

I want to estimate $$\widehat\beta$$ in a simple linear regression with scikit.

$$y = X \beta + \varepsilon$$

The problem is that the dimension of the complete $$X$$ is too large to fit into memory. Is there a way to split up the problem theoretically. Say, split up the sample into $$n$$ subamples randomly:

$$y^1, \ldots, y^n, X^1, \ldots, X^n$$

For each subsample get $$\widehat\beta^1, \ldots, \widehat\beta^n$$ and then combine these with some weighting function $$\theta$$ to get: $$\theta(\widehat\beta^1, \ldots, \widehat\beta^n) = \widehat\beta$$

Any suggestions on how to do that?

• I gave an answer without random splitting, which does not seem to add value in this context. – Christoph Hanck Dec 23 '18 at 16:11

Here is an indirect answer to the question that does not proceed by weighting coefficient estimates, but gradually computing the ingredients into the "global" estimate. I believe that biglm in R proceeds similarly.

Note that

$$\hat\beta=(X'X)^{-1}X'y$$ where $$X'X=\sum_{i=1}^Nx_ix_i'$$ with $$x_i'$$ the $$i$$th row of the regressor matrix. Likewise, $$X'y=\sum_ix_iy_i$$. Now suppose you divide $$N$$ into $$n$$ manageable blocks of size $$m$$, $$N=n\cdot m$$.

You may now load each block into memory consecutively. Then, with $$x_{j,k}'$$ the $$k$$th observation of block $$j$$, $$X'X=\sum_{j=1}^n\sum_{k=1}^mx_{j,k}x_{j,k}'$$ Each of the $$\sum_{k=1}^mx_{j,k}x_{j,k}'$$ yield an "unproblematic" matrix of size $$k\times k$$ ($$k$$ being the number of regressors, which I assume is not huge, but that the number of observations $$N$$ is).

A possible code snippet in R (sorry, I do not know python) might look like this - just for illustration, no claim to efficiency/elegance:

N <- 1000 # the "large" sample size
x <- rnorm(N)
y <- rnorm(N)
m <- 100
n <- N/m

x.sq <- x.y <- rep(NA,n)

for (j in 1:n){
x.sq[j] <- sum(x[((j-1)*m+1):(j*m)]^2)
x.y[j] <- sum(x[((j-1)*m+1):(j*m)]*y[((j-1)*m+1):(j*m)])
}
> (beta.hat <- sum(x.y)/sum(x.sq)) # recombining the samples and computing OLS (without intercept for simplicity)
[1] -0.04175772

> coef(lm(y~x-1)) # here, of course, everything fits into memory and we can check equality
x
-0.04175772