每个3D numpy数组元素的高效1D线性回归

Question

我有一组掩膜数组的3D堆叠。我想对沿着轴0（时间）的每个行、列（空间索引）的值执行线性回归。这些堆栈的维度是不同的，但一个典型的形状可能是（50, 2000, 2000）。我的空间限制但时间密集的测试案例具有以下维度：

stack.ma_stack.shape

(1461, 390, 327)

我进行了一个快速测试，循环遍历每一行和列：

from scipy.stats.mstats import linregress
#Ordinal dates
x = stack.date_list_o
#Note: idx should be row, col
def sample_lstsq(idx):
    b = stack.ma_stack[:, idx[0], idx[1]]
    #Note, this is masked stats version
    slope, intercept, r_value, p_value, std_err = linregress(x, b)
    return slope

out = np.zeros_like(stack.ma_stack[0])
for row in np.arange(stack.ma_stack.shape[1]):
    for col in np.arange(stack.ma_stack.shape[2]):
        out[row, col] = sample_lstsq((row, col))

我知道这个方法可以运行（但速度慢）。我知道肯定有更高效的方法。我已经尝试使用索引数组和np.vectorize，但我认为这并不能提供任何真正的改进。我考虑将所有内容转储到Pandas或尝试移植到Cython，但我希望能坚持使用NumPy/SciPy。或者并行解决方案是我改善性能的最佳选择？此外，有人对NumPy/SciPy线性回归选项进行过基准测试吗？我找到了以下选项，但还没有自己进行测试：

• scipy.stats.linregress
• numpy.linalg.leastsq
• numpy.polyfit(deg=1)

我希望有一种现有的方法可以在不需要太多实现工作的情况下显著提高性能。谢谢。

---

编辑于12/3/13 @02:29

@HYRY建议的方法对于上述样本数据集非常有效（运行时间约为15秒），该数据集在所有维度（空间和时间）上均为连续（未屏蔽）。 但是，当将包含缺失数据的掩蔽数组传递给np.linalg.leastsq时，所有屏蔽值都将填充为fill_value（默认值1E20），这会导致虚假的线性拟合。

幸运的是，numpy屏蔽数组模块具有np.ma.polyfit（deg = 1），它可以处理类似于np.linalg.leastsq的2D y数组。 查看源代码（https://github.com/numpy/numpy/blob/v1.8.0/numpy/ma/extras.py#L1852），ma polyfit只是一个使用x和y掩码的组合掩码的np.polyfit包装器。 当y中的缺失数据位置恒定时，这对于2D y非常有效。

不幸的是，我的数据在空间和时间上具有可变的缺失数据位置。以下是另一个堆栈的示例：

In [146]: stack.ma_stack.shape
Out [146]: (57, 1889, 1566)

抽样单个索引会返回一个时间序列，其中有6个未屏蔽的值：

In [147]: stack.ma_stack[:,0,0]
Out [147]: 
masked_array(data = [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 519.7779541015625 -- -- -- 518.9047241210938 -- -- -- -- -- -- --
 516.6539306640625 516.0836181640625 515.9403686523438 -- -- -- --
 514.85205078125 -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --],
             mask = [ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True False  True  True
  True False  True  True  True  True  True  True  True False False False
  True  True  True  True False  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True],
       fill_value = 1e+20)

采样不同的位置会返回不同时间切片中未屏蔽值的数量：

In [148]: stack.ma_stack[:,1888,1565]
Out[148]:
masked_array(data = [-- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 729.0936889648438 -- -- -- 724.7155151367188 -- -- -- -- -- -- --
 722.076171875 720.9276733398438 721.9603881835938 -- 720.3294067382812 --
 -- 713.9591064453125 709.8037719726562 707.756103515625 -- -- --
 703.662353515625 -- -- -- -- 708.6276245117188 -- -- -- -- --],
             mask = [ True  True  True  True  True  True  True  True  True  True  True  True
  True  True  True  True  True  True  True  True  True False  True  True
  True False  True  True  True  True  True  True  True False False False
  True False  True  True False False False  True  True  True False  True
  True  True  True False  True  True  True  True  True],
       fill_value = 1e+20)

每个索引的最小未遮盖值为6，最大值为45。因此，每个位置至少有一些被遮盖的值。 作为参考，我的x（时间序数）值都没有被遮盖：

In [150]: stack.date_list_o
Out[150]:
masked_array(data = [ 733197.64375     733962.64861111  733964.65694444  733996.62361111
  733999.64236111  734001.63541667  734033.64305556  734071.64722222
  734214.675       734215.65694444  734216.625       734226.64722222
  734229.63819444  734232.65694444  734233.67847222  734238.63055556
  734238.63055556  734245.65277778  734245.65277778  734255.63125
  734255.63125     734307.85        734326.65138889  734348.63888889
  734348.63958333  734351.85        734363.70763889  734364.65486111
  734390.64722222  734391.63194444  734394.65138889  734407.64652778
  734407.64722222  734494.85        734527.85        734582.85
  734602.65486111  734664.85555556  734692.64027778  734741.63541667
  734747.85        734807.85555556  734884.85555556  734911.65763889
  734913.64375     734917.64236111  734928.85555556  734944.71388889
  734961.62777778  735016.04583333  735016.62777778  735016.85555556
  735036.65347222  735054.04583333  735102.63125     735119.61180556
  735140.63263889],
             mask = False,
       fill_value = 1e+20)

所以我重塑了stack.ma_stack并运行了polyfit:

newshape = (stack.ma_stack.shape[0], stack.ma_stack.shape[1]*stack.ma_stack.shape[2])
print newshape
#(57, 2958174)
y = stack.ma_stack.reshape(newshape)
p = np.ma.polyfit(x, y, deg=1)

但是在第1500列左右，y中的每一行都被“累积”掩码，导致一些投诉和空输出：

RankWarning: Polyfit may be poorly conditioned
** On entry to DLASCL, parameter number  4 had an illegal value
...

因此，似乎在不同位置使用具有缺失数据的二维y是行不通的。我需要一个leastsq拟合，该拟合使用每个y列中的所有可用未屏蔽数据。也许可以通过仔细压缩x和y并跟踪未屏蔽索引来实现这一点。

还有其他想法吗？Pandas看起来可能是一个不错的解决方案。

---

编辑于12/3/13 @20:29

@HYRY提供了一个解决方案，适用于时间（轴=0）维度中的缺失值。 我必须稍微修改一下以处理空间（轴=1,2）维度中的缺失值。 如果特定的空间索引在时间上仅有一个未屏蔽的条目，我们肯定不想尝试进行线性回归。 这是我的实现：

def linreg(self):
    #Only compute where we have n_min unmasked values in time
    n_min = 3
    valid_idx = self.ma_stack.count(axis=0).filled(0) >= n_min
    #Returns 2D array of unmasked columns
    y = self.ma_stack[:, valid_idx]

    #Extract mask for axis 0 - invert, True where data is available
    mask = ~y.mask
    #Remove masks, fills with fill_value
    y = y.data
    #Independent variable is time ordinal
    x = self.date_list_o
    x = x.data

    #Prepare matrices and solve
    X = np.c_[x, np.ones_like(x)]
    a = np.swapaxes(np.dot(X.T, (X[None, :, :] * mask.T[:, :, None])), 0, 1)
    b = np.dot(X.T, (mask*y))
    r = np.linalg.solve(a, b.T)

    #Create output grid with original dimensions
    out = np.ma.masked_all_like(self.ma_stack[0])
    #Fill in the valid indices
    out[valid_idx] = r[:,0]

运行时非常快 - 仅 ~5-10 秒钟处理这里讨论的数组维度。

- David Shean

---

你的时间点是规律的（例如0秒，1秒，2秒，3秒），还是不规律的（例如0秒，1秒，3秒，10秒）？这将影响你选择回归函数的方式。 - Aaron Lockey

@Aaron，针对我上面描述的堆栈，是的，时间片是每天的，并且在时间或空间上没有缺失值。然而，我的大多数堆栈在时间和空间上都是不规则采样的，同一个堆栈包含2年、2天和2小时的间隔。你有什么建议？ - David Shean

只想在此帖子中添加一个链接，该链接基本上询问相同的问题，但针对无缺失数据的未屏蔽数组：https://dev59.com/5XfZa4cB1Zd3GeqPWMxW?rq=1 - David Shean

Answer 1

如果我理解正确，您想要进行多元线性回归 y = k * x + b，只有一个 x，但有许多个 y，对于每个 y 您都想计算出 k 和 b。

如果 x 的形状为 50，y 的形状为 (50,1000)，您可以使用 numpy.linalg.lstsq，以下是一些演示：

import numpy as np

x = np.random.rand(50)
k = np.random.rand(1000)
b = np.random.rand(1000)

y = np.outer(x, k) + b + np.random.normal(size=(50, 1000), scale=1e-10)

r = np.linalg.lstsq(np.c_[x, np.ones_like(x)], y)[0]

print np.allclose(r[0], k)
print np.allclose(r[1], b)

对于形状为(50, 2000, 2000)的Y，您可以将其重新塑造为(50, 2000*2000)。

编辑

这是我对掩码数组的解决方案。公式如下：

将Y准备为形状为(1889*1566, 57)的2维数组，X准备为形状为(57, 2)的2维数组。掩码是一个布尔数组，与Y具有相同的形状，True表示Y中的值可用。

计算形状为(1889*1566, 2, 2)的数组a，形状为(1889*1566, 2)的数组b，然后调用numpy.linalg.solve(a, b)，以下是一些演示代码：

import numpy as np

N = 50
M = 1000

x = np.random.rand(N)
X = np.c_[x, np.ones_like(x)]
beta = np.random.rand(M, 2)
Y = np.dot(beta, X.T)
Y += np.random.normal(scale=0.1, size=Y.shape)
mask = np.random.randint(0, 2, size=Y.shape).astype(np.bool)

a = np.swapaxes(np.dot(X.T, (X[None, :, :] * mask[:, :, None])), 0, 1)
b = np.dot(X.T, (mask*Y).T)
beta2 = np.linalg.solve(a, b.T)

i = 123
print "real:", beta[i]
print "by solve:", beta2[i]

m = mask[i]
x2 = X[m]
y2 = Y[i, m]
print "by lstsq:", np.linalg.lstsq(x2, y2)[0]

输出第123个系数：

real: [ 0.35813131  0.29736779]
by solve: [ 0.38088499  0.30382547]
by lstsq: [ 0.38088499  0.30382547]

你也可以按照以下代码计算 a 数组，我认为它会比上面的方法使用更少的内存:

a2 = np.empty((M, 2, 2))
xm = mask * x
a2[:, 0, 0] = (xm*xm).sum(1)
a2[:, 1, 0] = (xm*mask).sum(1)
a2[:, 0, 1] = a2[:, 1, 0]
a2[:, 1, 1] = (mask).sum(1)

print np.allclose(a2, a)

- HYRY

---

哇，这就是我所期望的答案。不幸的是，np.linalg.lstsq无法处理缺失数据。我将编辑原帖并发布一些关于np.ma.polyfit测试结果的信息。 - David Shean

你使用的numpy版本是什么，在numpy 1.8中，numpy.linalg.solve是一个广义ufunc，你可以通过一次调用它来解决M个线性方程。而且lstsq可以通过solve()计算得出。稍后我会发布一些代码。 - HYRY

我正在使用“1.9.0.dev-Unknown”。感谢您的帮助。 - David Shean

numpy.linalg.solve看起来可以工作，如果我生成另一个57x2958174的数组，其中包含与压缩的y列对应的压缩x值。这是你想要的吗？感谢你的帮助。 - David Shean

你的解决方案在应用掩码时与np.linalg.solve方法很好地配合使用。对于我的掩码数组输入，我不得不反转现有的掩码（np.ma掩码在数据缺失时为True）。我还必须将我的x和y从np.ma转换为np.array以去除掩码，防止numpy在后续操作中自动填充填充值（1E20）。在某个时候，我会尝试回过头来稍微整理一下，但是你的解决方案很有效，而且我现在面临着截止日期的压力。非常感谢你。 - David Shean