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?

  • $\begingroup$ I gave an answer without random splitting, which does not seem to add value in this context. $\endgroup$ – 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.