Python：仅使用 1 个外生变量执行数百万个简单线性回归的最快方法

Question

我正在对时间序列数据执行组件明智的回归。这基本上是在那里的，而不是倒退Ÿ对X ₁，X ₂，...，X _ñ，我们会回归Ÿ对X ₁只，对x和y ₂只，...，并采取减少的总和回归平方残差最多，并将其添加为基础学习器。这被重复 M 次，这样最终模型是许多形式 y 与 x _i（仅 1 个外生变量）的简单线性回归的总和，基本上是使用线性回归作为基础学习器的梯度提升。

问题是，由于我正在对时间序列数据执行滚动窗口回归，因此我必须进行 N × M × T 回归，这超过一百万个 OLS。虽然每个 OLS 都非常快，但在我性能较弱的笔记本电脑上运行需要几个小时。

目前，我正在使用statsmodels.OLS.fit()这样的方式来获取每个 y 的参数，而不是 x _i线性回归。的z_matrix是数据矩阵和i表示第i^个列片的回归。行数约为 100，z_matrix大小约为 100 × 500。

    ols_model = sm.OLS(endog=endog, exog=self.z_matrix[:, i][..., None]).fit()
    return ols_model.params, ols_model.ssr, ols_model.fittedvalues[..., None]

我从 2016 年的上一篇文章中阅读了在 python 中计算许多回归的最快方法？使用重复调用 statsmodels 效率不高，我尝试了建议 numpy 的答案之一，但pinv不幸的是速度较慢：

    # slower: 40sec vs 30sec for statsmodel for 100 repeated runs of 150 linear regressions
    params = np.linalg.pinv(self.z_matrix[:, [i]]).dot(endog)
    y_hat = self.z_matrix[:, [i]]@params
    ssr = sum((y_hat-endog)**2)
    return params, ssr, y_hat

有没有人有更好的建议来加速线性回归的计算？我只需要估计的参数、残差平方和和预测的 ? 价值。谢谢！

Answer 1

这是一种方法，因为您总是在没有常数的情况下运行回归。此代码在大约 0.5 秒内运行大约 900K 个模型。它保留 sse、每个 900K 回归的预测值以及估计参数。

最重要的想法是利用一个变量对另一个变量的回归背后的数学原理，即叉积与内积（模型不包含常数）的比率。可以通过使用移动窗口 demean 来修改它以包括常数来估计截距。

import numpy as np
from statsmodels.regression.linear_model import OLS
import datetime

gen = np.random.default_rng(20210514)

# Number of observations
n = 1000
# Number of predictors
m = 1000
# Window size
w = 100
# Simulate data
y = gen.standard_normal((n, 1))
x = gen.standard_normal((n, m))

now = datetime.datetime.now()
# Compute rolling covariance and variance-like terms
# These assume the model is y = x*b + e w/o a constant
c = np.r_[np.zeros((1, m)), np.cumsum(x * y, axis=0)]
v = np.r_[np.zeros((1, m)), np.cumsum(x * x, axis=0)]
c_trimmed = c[w:] - c[:-w]
v_trimmed = v[w:] - v[:-w]
# Parameters are just the ratio
params = c_trimmed / v_trimmed

# Build a selector array to quickly reshape y and the columns of x
step = np.arange(m - w + 1)
sel = np.arange(w)
locs = step[:, None] + sel
# Get the blocked reshape of y. It has n - w + 1 rows with window observations
# and looks like
# [[y[0],y[1],...,y[99]],
#  [y[1],y[2],...,y[100]],
#  ...,
#  [y[900],y[901],...,y[999]],
y_block = y[locs, 0]
# Storage for the predicted values and the sse
y_pred = np.empty((x.shape[1],) + y_block.shape)
sse = np.empty((m - w + 1, n))
# Easiest to loop over columns.
# Could do broadcasting tricks, but noth worth the trouble since number of columns is modest
for i in range(x.shape[0]):
    # Reshape a columns of x like y
    x_block = x[locs, i]
    # Get the parameters and make sure it is 2d with shape (m-w+1, 1)
    # so the broadcasting works
    p = params[:, i][:, None]
    # Get the predicted value
    y_pred[i] = x_block * p
    # And the sse
    sse[:, i] = ((y_block - y_pred[i]) ** 2).sum(1)
print(f"Time: {(datetime.datetime.now() - now).total_seconds()}s")
# Some test code
# Test any single observation
start = 124
assert start <= m - w
column = 342
assert column < x.shape[1]
res = OLS(y[start : start + 100], x[start : start + 100, [column]]).fit()
np.testing.assert_allclose(res.params[0], params[start, column])
np.testing.assert_allclose(res.fittedvalues, y_pred[column, start])
np.testing.assert_allclose(res.ssr, sse[start, column])