【问题标题】:Vectorizing nested loop with conditionals and functions使用条件和函数矢量化嵌套循环
【发布时间】:2020-04-05 07:41:13
【问题描述】:

我有以下功能:

def F(x):    #F receives a numpy vector (x) with size (xsize*ysize)

    ff = np.zeros(xsize*ysize)

    count=0 

    for i in range(xsize):
       for j in range(ysize):

           a=function(i,j,xsize,ysize)

           if (a>xsize):
               ff[count] = x[count]*a
           else
               ff[count] = x[count]*i*j

           count = count +1

    return ff      

这里有一个细微差别,即(例如 xsize =4, ysize=3)

c=count
x[c=0] corresponds to x00 i=0,j=0
x[c=1]   x01   i=0, j=1
x[c=2]   x02   i=0, j=2   (i=0,  j = ysize-1)
x[c=3]   x10   i=1, j=0
...   ...  ...
x[c=n]   x32    i=3 j=2  (i=xsize-1, j=ysize-1)  

我的代码很简单

ff[c] = F[x[c]*a (condition 1)
ff[c] = F[x[c]*i*j (condition 2)

我可以避免使用广播的嵌套循环,如此链接中所述:

Python3: vectorizing nested loops

但在这种情况下,我必须调用函数(i,j,xsize,ysize),然后我有条件。 我真的需要知道 i 和 j 的值。

这个函数可以向量化吗?

编辑:function(i,j,xsize,ysize) 将使用 sympy 执行符号计算以返回浮点数。所以a 是一个浮点数,而不是一个符号表达式。

【问题讨论】:

  • function 有什么作用?看看这是否可以矢量化是相当重要的。
  • 应该将count += 1缩进到内部循环中吗?现在它只在外循环中,这意味着你的大部分任务都覆盖了自己,并不是所有的 ff 都被初始化了。
  • 功能比较复杂。它涉及点积和矩阵乘法的符号计算。
  • symbolic calculations?和sympy一样?如果是,请调整标签和示例。
  • 大家,function 只取决于输入的大小。这里只需要一个给定大小的函数结果的缓存,然后将其重用于相同大小的 x 的不同值。

标签: python python-3.x numpy vectorization sympy


【解决方案1】:

首先要注意的是,您的函数F(x) 可以描述为每个索引的x(idx) * weight(idx),其中权重仅取决于x 的维度。所以让我们根据函数get_weights_for_shape 来构造我们的代码,这样F 就相当简单了。为简单起见,weights 将是一个 (xsize by size) 矩阵,但我们也可以让 F 用于平面输入:

def F(x, xsize=None, ysize=None):
    if len(x.shape) == 2:
        # based on how you have put together your question this seems like the most reasonable representation.
        weights = get_weights_for_shape(*x.shape)
        return x * weights
    elif len(x.shape) == 1 and xsize * ysize == x.shape[0]:
        # single dimensional input with explicit size, use flattened weights.
        weights = get_weights_for_shape(xsize, ysize)
        return x * weights.flatten()
    else:
        raise TypeError("must take 2D input or 1d input with valid xsize and ysize")


# note that get_one_weight=function can be replaced with your actual function.
def get_weights_for_shape(xsize, ysize, get_one_weight=function):
    """returns weights matrix for F for given input shape"""
    # will use (xsize, ysize) shape for these calculations.
    weights = np.zeros((xsize,ysize))
    #TODO: will fill in calculations here
    return weights

所以首先我们要为每个元素运行你的function(我在参数中使用了别名get_one_weight),你说这个函数不能向量化,所以我们可以只使用列表推导。我们想要一个具有相同形状(xsize,ysize) 的矩阵a,因此对于嵌套列表的理解有点倒退:

# notice that the nested list makes the loops in opposite order:
# [ROW for i in Xs]
#  ROW = [f() for j in Ys]
a = np.array([[get_one_weight(i,j,xsize,ysize)
                    for j in range(ysize)
              ] for i in range(xsize)])

有了这个矩阵a > xsize 会给出一个布尔数组进行条件赋值:

case1 = a > xsize
weights[case1] = a[case1]

对于另一种情况,我们使用索引ij。要矢量化二维索引,我们可以使用np.meshgrid

[i,j] = np.meshgrid(range(xsize), range(ysize), indexing='ij')
case2 = ~case1 # could have other cases, in this case it's just the rest.
weights[case2] = i[case2] * j[case2]

return weights #that covers all the calculations

把它们放在一起就可以得到一个完全向量化的函数:

# note that get_one_weight=function can be replaced with your actual function.
def get_weights_for_shape(xsize, ysize, get_one_weight=function):
    """returns weights matrix for F for given input shape"""
    # will use (xsize, ysize) shape for these calculations.
    weights = np.zeros((xsize,ysize))

    # notice that the nested list makes the loop order confusing:
    # [ROW for i in Xs]
    #  ROW = [f() for j in Ys]
    a = np.array([[get_one_weight(i,j,xsize,ysize)
                        for j in range(ysize)
                  ] for i in range(xsize)])

    case1 = (a > xsize)
    weights[case1] = a[case1]

    # meshgrid lets us use indices i and j as vectorized matrices.
    [i,j] = np.meshgrid(range(xsize), range(ysize), indexing='ij')
    case2 = ~case1
    weights[case2] = i[case2] * j[case2]
    #could have more than 2 cases if applicable.

    return weights

这涵盖了大部分内容。对于您的具体情况,因为这种繁重的计算仅依赖于输入的形状,如果您希望使用类似大小的输入重复调用此函数,您可以缓存所有先前计算的权重:

def get_weights_for_shape(xsize, ysize, _cached_weights={}):
    if (xsize, ysize) not in _cached_weights:
        #assume we added an underscore to real function written above
        _cached_weights[xsize,ysize] = _get_weights_for_shape(xsize, ysize)
    return _cached_weights[xsize,ysize]

据我所知,这似乎是您将获得的最优化。唯一的改进是对function 进行矢量化(即使这意味着只是在多个线程中并行调用它),或者如果.flatten() 制作了一个可以改进但我不完全确定如何改进的昂贵副本。

【讨论】:

  • 我无法执行此操作:(xsize,ysize) = dims。它给了我一个错误:没有足够的值来解包(预期 2,得到 1)
  • 你没有展示你是如何得到xsizeysize的所以我认为它是直接从x.shape获取的,如果你的数据只有一维那么是@987654349 @ = 1 ?
  • 我的解释很糟糕。我的向量 x 总是一维的。它的大小将取决于只是两个整数的 xsize 和 ysize。如果 ysize = 2 且 xsize =3,则一维数组的大小将为 ysize*xsize。这里的问题是,对于 x[0],我们有 (i=0, j=0) x[1] (i=0, j=1)。这就是为什么我有一个嵌套循环的原因。如果我可以直接从 c 中获取 ij,我就不需要循环了。
  • 如所写,您只需传递F(x.reshape((xsize, ysize)) 即可获得您想要的。我想提请注意我在哪里使用np.meshgrid,因为这为我提供了一个矢量化索引以用于计算,尽管当您的输入与循环的形状匹配时效果会更好。
  • 它有效。它比第一个代码慢,但我认为它以更好的方式扩展。对于一定的大小,它会比我猜的第一个代码更快。但是你能解释一下重用是什么意思吗?
猜你喜欢
  • 1970-01-01
  • 2017-02-01
  • 1970-01-01
  • 1970-01-01
  • 2014-09-21
  • 2016-03-20
  • 2015-11-01
  • 1970-01-01
  • 2022-12-07
相关资源
最近更新 更多