How to run linear regression in a parallel/distributed way for big data setting?

Question

I am working on a very large linear regression problem, with data size so large that they have to be stored on a cluster of machines. It will be way too big to aggregate all the samples into one single machine's memory (even disk)

To do regression these data, I am thinking about a parallel approach, i.e., run regression on each individual box, and then calculate the beta based on the statistics of each individual beta (probably a mean or a median)

does this make any sense ? if so, how should I get the total expected $R^2$ from each individual $R^2$?

Answer 1

Short Answer:

Yes, running linear regression in parallel has been done. For example, Xiangrui Meng et al. (2016) for Machine Learning in Apache Spark. The way it works is using stochastic gradient descent (SGD). In section 3, core features, the author mentioned:

Generalized linear models are learned via optimization algorithms which parallelize gradient computation, using fast C++-based linear algebra libraries for worker computations.

An example on how SGD works can be found in my answer here: How could stochastic gradient descent save time comparing to standard gradient descent?

---

Long Answer:

Note, the notation is not consistent with the link I provided, I feel matrix notation is better in this question.

To do a linear regression we are trying to do

$$\text{minimize}~\|X\beta-y\|^2$$

The derivative is

$$2X^T(X\beta-y)$$

In small data settings, we can set the derivative to $0$ and solve it directly. (e.g., QR decomposition in R.) In big data settings, the data matrix $X$ is too big to be stored in memory, and may be hard to solve directly. (I am not familiar with how to do QR decomposition or Cholesky decomposition for huge matrices).

One way to parallelize this is by trying to use an iterative method: stochastic gradient descent, where we can approximate the gradient using a subset of the data. (If we use $X_s$, $y_s$ to represent a subset of the data, the gradient can be approximated by $2X_s^T(X_s\beta-y_s)$, and we can update $\beta$ with the approximated gradient).

In addition, for the $R^2$ statistic, we can compute $R^2$ for all data in parallel or approximate it by using a subset of the data.

Intuition on how it works (mapreduce paradigm):

I keep saying approximation using a subset; the intuition for why this works can be described in the following example: suppose I have 100 billion data points and we want to calculate the average of all data points. Suppose conducting such an operation takes a very long time, and further that the whole data cannot be stored in memory.

What we can do is to just take a subset, say 1 billion items, and calculate the average of these. The approximation thus produced should not be far away from the truth (i.e., using the whole data).

To parallelize, we can use 100 computers, with each of them taking a different subset of the 1 billion data points and calculating the average of these. (Commonly referred to as the MAP step). Finally, run another average on these 100 numbers (a.k.a. the REDUCE step).

Note the "mapreduce paradigm" would work well in some cases, but not well in others. For the example, the "average" operation mentioned earlier is very easy, because we know $\text{mean}(<x,y>)=\text{mean}(x)+\text{mean(y)}$, (assuming the length of $x$ and $y$ are the same). For some iterative methods, i.e., the current iteration is dependent on previous iteration results, it is hard to parallelize. Stochastic gradient descent solves this problem by approximating the gradient using a subset of data. And details can be found in @user20160 's answer.

References:

Xiangrui Meng et al. (2016). MLlib: Machine Learning in Apache Spark

Answer 2

As @hxd1011 mentioned, one approach is to formulate linear regression as an optimization problem, then solve it using an iterative algorithm (e.g. stochastic gradient descent). This approach can be parallelized but there are a couple important questions: 1) How should be problem be broken into subproblems? 2) Given that optimization algorithms like SGD are inherently sequential, how should solutions to the subproblems be combined to obtain a global solution?

Zinkevich et al. (2010) describe some previous approaches to parallelizing across multiple machines:

• 1) Parallelize SGD as follows: Split the data across multiple machines. At each step, each local machine estimates the gradient using a subset of the data. All gradient estimates are passed to a central machine, which aggregates them to perform a global parameter update. The downside of this approach is that it requires heavy network communication, which reduces efficiency.

• 2) Partition the data evenly across local machines. Each machine solves the problem exactly for its own subset of the data, using a batch solver. Final parameter estimates from the local machines are averaged to produce a global solution. The benefit of this approach is that it requires very little network communication, but the downside is that the parameter estimates can be suboptimal.

They propose a new approach:

• 3) Allow each local machine to randomly draw data points. Run SGD on each machine. Finally, average the parameters across machines to obtain a global solution. Like (2), this method requires little network communication. But, the parameter estimates are better because each machine is allowed to access a larger fraction of the data.

The parallelized optimization approach is very general, and applies to many machine learning algorithms (not just linear regression).

Another alternative would be to use parallel/distributed matrix decomposition algorithms or linear solvers. Least squares linear regression has special structure that allows it to be solved using matrix decomposition methods. This is how you'd typically solve it in the case of a smaller data set that fits in memory. This can be parallelized by distributing blocks of the matrix across multiple machines, then solving the problem using parallel/distributed matrix computations. Given that this approach is more specialized to solving linear systems, it would be interesting to see how its performance compares to the more general distributed optimization approach. If anyone can provide more information about this, I'd be glad to hear.

References:

Zinkevich et al. (2010). Parallelized Stochastic Gradient Descent.

Answer 3

Long, long, before map reduce I solved this. Below is reference to an old paper of mine in Journal of Econometrics 1980. It was for parallel nonlinear maximum likelihood and would work for M-estimation.

The method is exact for regressions. Split data into k subsets on k processors/units (could be done sequentially as well.) Do k regressions keep the regression coefficients an the X'X matrix for each. Call these b1,...,bk and W1,...,Wk respectively then the overall regression coefficients is given by b=inverse(W1+..+Wk)*(W1*b1+...+Wk*bk) one needs another pass through the data to calculate the residuals using b for the parameters to get sigma^2 the estimated error variance, R^2 overall F and the like. Then the covariance matrix of b is given exactly by sigma^2 (inverse(W1+..+Wk)). Above * indicates matrix multiplication.

https://www.sciencedirect.com/science/article/pii/0304407680900950

Answer 4

As far as I understand, the formulas for the intercept and the slope in a simple linear regression model can be rewritten as expressions on various sums that don't contain the mean value, so that these sums can be calculated in parallel during a map phase or similar:

• sum(x)
• sum(y)
• sum(x*x)
• sum(x*y)

In a final phase that is done in a single instance, the actual regression is calculated:

slope = beta = (sum(x*y) - sum(x)*sum(y)/n)/(sum(x*x) - sum(x)*sum(x)/n)

intercept = alpha = (sum(y) - beta*sum(x))/n

A similar approach could be used to calculate linear regression parameters using SQL aggregates.

Sorry, I haven't tested this in great detail, but simulating a few cases in an Excel spreadsheet resulted in the same regression line that Excel produced.