简短的回答：

是的，已经完成了并行运行线性回归。例如，祥瑞孟等人。(2016) 用于 Apache Spark 中的机器学习。它的工作方式是使用随机梯度下降 (SGD)。在第 3 节，核心特征中，作者提到：

广义线性模型是通过并行梯度计算的优化算法学习的，使用基于 C++ 的快速线性代数库进行工作计算。

关于 SGD 如何工作的示例可以在我的回答中找到：与标准梯度下降相比，随机梯度下降如何节省时间？

长答案：

请注意，符号与我提供的链接不一致，我觉得在这个问题中矩阵符号更好。

要进行线性回归，我们正在尝试做

导数是

在小数据设置中，我们可以将导数设为，直接求解。（例如，R 中的 QR 分解。）在大数据设置中，数据矩阵太大而无法存储在内存中，并且可能难以直接求解。（我不熟悉如何对大型矩阵进行 QR 分解或 Cholesky 分解）。$0$$X$

一种并行化的方法是尝试使用迭代方法：随机梯度下降，我们可以使用数据的子集来近似梯度。（如果我们用 ,来表示数据的一个子集，梯度可以近似为，我们可以用近似的梯度更新β）$X_s$$y_s$$2X_s^T(X_s\beta-y_s)$$\beta$

此外，对于统计量，我们可以并行计算 ^2，也可以使用数据的子集对其进行近似计算。$R^2$$R^2$

关于它如何工作的直觉（mapreduce 范式）：

我一直说使用子集进行近似；可以在以下示例中描述为什么这样做的直觉：假设我有 1000 亿个数据点，我们想要计算所有数据点的平均值。假设进行这样的操作需要很长时间，而且整个数据无法存储在内存中。

我们可以做的是只取一个子集，比如 10 亿个项目，然后计算它们的平均值。这样产生的近似值不应该远离真相（即，使用整个数据）。

为了实现并行化，我们可以使用 100 台计算机，每台计算机获取 10 亿个数据点的不同子集并计算它们的平均值。（通常称为 MAP 步骤）。最后，对这 100 个数字运行另一个平均值（也称为 REDUCE 步骤）。

请注意，“mapreduce 范例”在某些情况下效果很好，但在其他情况下效果不佳。例如，前面提到的“平均”操作非常简单，因为我们知道， (假设和的长度相同）。对于一些迭代方法，即当前迭代依赖于先前的迭代结果，很难并行化。随机梯度下降通过使用数据子集逼近梯度来解决这个问题。详细信息可以在@user20160 的回答中找到。$\text{mean}(<x,y>)=\text{mean}(x)+\text{mean(y)}$$x$$y$

参考：

祥瑞孟等。（2016 年）。MLlib：Apache Spark 中的机器学习

正如@hxd1011 提到的，一种方法是将线性回归公式化为优化问题，然后使用迭代算法（例如随机梯度下降）解决它。这种方法可以并行化，但有几个重要问题：1）应该如何将问题分解为子问题？2）鉴于像SGD这样的优化算法本质上是顺序的，应该如何组合子问题的解决方案以获得全局解决方案？

津克维奇等人。(2010) 描述了一些以前的跨多台机器并行化的方法：

• 1) 将 SGD 并行化如下：将数据拆分到多台机器上。在每一步，每台本地机器使用数据子集估计梯度。所有梯度估计都传递到中央机器，中央机器将它们聚合以执行全局参数更新。这种方法的缺点是它需要繁重的网络通信，从而降低了效率。

• 2) 在本地机器上均匀地划分数据。每台机器都使用批处理求解器精确地为自己的数据子集解决问题。来自本地机器的最终参数估计被平均以产生全局解决方案。这种方法的好处是它需要很少的网络通信，但缺点是参数估计可能不是最优的。

他们提出了一种新方法：

• 3）让每台本地机器随机抽取数据点。在每台机器上运行 SGD。最后，对机器之间的参数进行平均以获得全局解决方案。与（2）一样，这种方法需要很少的网络通信。但是，参数估计更好，因为允许每台机器访问更大比例的数据。

并行优化方法非常通用，适用于许多机器学习算法（不仅仅是线性回归）。

另一种选择是使用并行/分布式矩阵分解算法或线性求解器。最小二乘线性回归具有允许使用矩阵分解方法求解的特殊结构。在适合内存的较小数据集的情况下，这通常是您解决问题的方法。这可以通过在多台机器上分布矩阵块来并行化，然后使用并行/分布式矩阵计算来解决问题。鉴于这种方法更专门用于解决线性系统，看看它的性能如何与更一般的分布式优化方法进行比较会很有趣。如果有人可以提供有关此的更多信息，我将很高兴听到。

参考：

津克维奇等人。（2010 年）。并行随机梯度下降。

很久很久了，在 map reduce 我解决了这个问题之前。下面是我在 1980 年计量经济学杂志上的一篇旧论文的参考。它适用于并行非线性最大似然，适用于 M 估计。

该方法适用于回归。将数据拆分为 k 个处理器/单元上的 k 个子集（也可以按顺序完成）。k 个回归是否将回归系数保持为每个 X'X 矩阵。分别调用这些 b1,...,bk 和 W1,...,Wk 然后整体回归系数由 b=inverse(W1+..+Wk)*(W1*b1+...+Wk*bk) 给出需要再次通过数据来计算残差，使用 b 作为参数来获得 sigma^2 估计误差方差、R^2 总体 F 等。然后 b 的协方差矩阵完全由 sigma^2 (inverse(W1+..+Wk)) 给出。* 以上表示矩阵乘法。

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

据我了解，简单线性回归模型中的截距和斜率公式可以重写为不包含平均值的各种总和的表达式，以便可以在映射阶段并行计算这些总和或类似的：

• 总和（x）
• 总和（y）
• 总和（x*x）
• 总和（x*y）

在单个实例中完成的最后阶段，计算实际回归：

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

类似的方法可用于使用 SQL 聚合计算线性回归参数。

抱歉，我没有对此进行详细测试，但是在 Excel 电子表格中模拟一些案例会产生与 Excel 生成的相同回归线。