如何解决加快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
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。