加快for循环-Python

问题描述

我有一个运行良好的代码，但我希望加快收敛的时间。代码片段如下所示：

def myfunction(x,i):
    y = x + (min(0,target[i] - data[i,:]x))*data[i]/(norm(data[i])**2))
    return y
rows,columns = data.shape

start = time.time()
iterate = 0
iterate_count = []
norm_count = []
res = 5

x_not = np.ones(columns)
norm_count.append(norm(x_not))
iterate_count.append(0)

while res > 1e-8:
    for row in range(rows):
        y = myfunction(x_not,row)
        x_not = y
    iterate += 1
    iterate_count.append(iterate)
    norm_count.append(norm(x_not))
    res = abs(norm_count[-1] - norm_count[-2])

print('Converge at {} iterations'.format(iterate))
print('Duration: {:.4f} seconds'.format(time.time() - start))

我在Python中相对较新。我将不胜感激。

Ax = b是我们要解决的问题。在这里，“ A”是“数据”，“ b”是“目标”

解决方法

U！花了一段时间后，我认为无法通过设置问题的方式来解决。在该行的每次迭代中，您都修改x_not，然后传递更新后的结果以获得下一行的解决方案。这种设置无法轻易向量化。您可以从失败的尝试中学到矢量化的思考过程，因此我将其包含在答案中。我还包括另一种迭代方法来求解线性方程组。我提供了一个矢量化版本-使用矩阵乘法和向量加法来更新解决方案，以及一个循环版本-使用for循环来更新解决方案以演示您可以期望获得的收益。 / p>

1。尝试失败

让我们看看您在这里做什么。

def myfunction(x,i): 
    y = x + (min(0,target[i] - data[i,:] @ x)) * (data[i] / (norm(data[i])**2)) 
    return y

• 您减去
  • （i和data的第x_not行的点积）
  • 来自i的第target行，
  • 限制为零。
• 您将此结果乘以i的第data行除以该行的范数的平方。我们称之为part2
• 然后将其添加到i的第x_not个元素

现在让我们看一下矩阵的形状。

• data是(M,N)。
• target是(M,)。
• x_not是(N,)

您可以对整个矩阵进行操作，而不必按行进行这些操作！

1.1。简化点积。

您可以执行data[i,:] @ x，而不是执行data @ x_not，这将为第i个元素提供一个数组，该数组将第i行与{{ 1}}。现在我们有了x_not，形状为data @ x_not

然后，您可以从整个(M,)数组中减去它，因此target的形状为target - (data @ x_not)。

到目前为止，我们有

(M,)

接下来，如果大于零，则将其设置为零。

part1 = target - (data @ x_not)

1.2。查找按行规范。

最后，您想将其乘以part1[part1 > 0] = 0 的行，然后除以该行L2范数的平方。要获取矩阵每一行的范数，您可以

data

这是一个rownorms = np.linalg.norm(data,axis=1) 数组，因此我们需要将其转换为(M,)数组，以便可以划分每一行。 (M,1)执行此操作。然后用rownorms[:,None]除以这个。

data

1.3。添加到part2 = data / (rownorms[:,None]**2)

最后，我们将 {em> x_not的每一行添加到原始part1 * part2并返回结果

x_not

在这里我们被困住了。在您的方法中，对result = x_not + (part1 * part2).sum(axis=0) 的每次调用都会给出一个取决于myfunction()的{​​{1}}值，该值在上一次对part1的调用中已更改。

2。为什么要向量化？

使用target[i]的内置方法而不是循环，可以将计算卸载到C后端，因此运行速度更快。如果您的numpy链接到BLAS backend，则可以使用处理器的SIMD寄存器you can extract even more speed

conjugate gradient method是解决某些方程组的简单迭代方法。还有其他更复杂的算法可以很好地解决通用系统问题，但这对于我们的演示来说应该是正确的。再次，目的不是要有一个能够完美解决任何线性方程组的迭代算法，而是要显示出对向量进行矢量化处理后可以期望的加速程度。

提供系统

myfunction()

让我们定义一些变量：

numpy

我们将解决系统data @ x_not = target

A = data.T @ data
b = data.T @ target

---

为了与完全矢量化的方法与使用迭代来更新A @ x = b和x = np.zeros((columns,)) # Initial guess. Can be anything resid = b - A @ x p = resid while (np.abs(resid) > tolerance).any(): Ap = A @ p alpha = (resid.T @ resid) / (p.T @ Ap) x = x + alpha * p resid_new = resid - alpha * Ap beta = (resid_new.T @ resid_new) / (resid.T @ resid) p = resid_new + beta * p resid = resid_new + 0 的行的方法进行对比，让我们定义执行此操作的CG解算器的另一种实现。

x

还有我们原始的矢量方法：

resid_new

让我们解决一个简单的系统，看看它是否首先起作用：

def solve_loopy(data,target,itermax = 100,tolerance = 1e-8):
    A = data.T @ data
    b = data.T @ target
    
    rows,columns = data.shape
        
    x = np.zeros((columns,)) # Initial guess. Can be anything
    resid = b - A @ x
    resid_new = b - A @ x
    p = resid
    
    niter = 0
    while (np.abs(resid) > tolerance).any() and niter < itermax:
        Ap = A @ p
        alpha = (resid.T @ resid) / (p.T @ Ap)
        for i in range(len(x)):
            x[i] = x[i] + alpha * p[i]
            resid_new[i] = resid[i] - alpha * Ap[i] 
        # resid_new = resid - alpha * A @ p
        beta = (resid_new.T @ resid_new) / (resid.T @ resid)
        p = resid_new + beta * p
        resid = resid_new + 0
        niter += 1
    
    return x

应提供def solve_vect(data,tolerance = 1e-8): A = data.T @ data b = data.T @ target rows,columns = data.shape x = np.zeros((columns,)) # Initial guess. Can be anything resid = b - A @ x resid_new = b - A @ x p = resid niter = 0 while (np.abs(resid) > tolerance).any() and niter < itermax: Ap = A @ p alpha = (resid.T @ resid) / (p.T @ Ap) x = x + alpha * p resid_new = resid - alpha * Ap beta = (resid_new.T @ resid_new) / (resid.T @ resid) p = resid_new + beta * p resid = resid_new + 0 niter += 1 return x

的解决方案

2x1 +   x2  =   -5
−x1 +   x2  =   -2

两者都给出正确的解决方案[-1,-3]，是的！现在开始做更大的事情：

data = np.array([[ 2,1],[-1,1]])
target = np.array([-5,-2])
print(solve_loopy(data,target))
print(solve_vect(data,target))

确保解决方案仍然正确：

[-1,-3]

嗯，看来CG方法对我们创建的条件差的随机矩阵不起作用。好吧，至少两者都给出相同结果。

data = np.random.random((100,100))
target = np.random.random((100,))

但不要气disc！我们并不在乎它是否完美运行，这是要证明向量化的惊人程度。因此，让我们计时一下：

sol1 = solve_loopy(data,target)
np.allclose(data @ sol1,target)
# Output: False

sol2 = solve_vect(data,target)
np.allclose(data @ sol2,target)
# Output: False

好！只需通过在更新我们的解决方案时避免循环，即可实现〜2倍的加速！

对于更大的系统，这会更好。

np.allclose(sol1,sol2)
# Output: True

这给我们：

import timeit

timeit.timeit('solve_loopy(data,target)',number=10,setup='from __main__ import solve_loopy,data,target')
# Output: 0.25586539999994784

timeit.timeit('solve_vect(data,setup='from __main__ import solve_vect,target')
# Output: 0.12008900000000722