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

Question

我正在对时间序列数据执行组件明智的回归。这基本上是我们将 y 仅对 x1、y 仅对 x2 进行回归，而不是将 y 对 x1、x2、...、xN 进行回归，并采用最能减少平方残差总和的回归，并且将其添加为基础学习器。这被重复 M 次，这样最终模型是许多形式 y 对 xi(仅 1 个外生变量)的简单线性回归的总和，基本上是使用线性回归作为基础学习器的梯度提升。
问题是，由于我正在对时间序列数据执行滚动窗口回归，因此我必须进行 N × M × T 回归，这超过一百万个 OLS。虽然每个 OLS 都非常快，但在我性能较弱的笔记本电脑上运行需要几个小时。
目前，我正在使用 statsmodels.OLS.fit()作为针对 xi 线性回归获取每个 y 参数的方法。 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 年的上一篇文章 Fastest way to calculate many regressions in 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 秒内运行了大约 90 万个模型。它保留了 sse、每个 900K 回归的预测值和估计参数。
一个重要的想法是利用一个变量对另一个变量的回归背后的数学，即叉积与内积的比率(模型不包含常数)。这可以修改为还包括一个常数，通过使用移动窗口降阶来估计截距。

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])