Python/Numpy 代码优化答案

【问题标题】：Python/Numpy Code OptimisationPython/Numpy 代码优化
【发布时间】：2014-06-24 19:09:46
【问题描述】：

我在大迭代中有以下两个矩阵代数计算。因此，我希望优化计算。

1：

    F = np.matrix(np.zeros(shape=(n+1,1)))
    F[0:n] = x - np.diag(np.array(theta)[0:n].flatten())*self.W*(theta[0:n]-self.m) + theta[0:n]*theta[n]
    F[n] = np.sum(theta[0:n]) - 1; #Lagrange multiplier term

2：

    J = np.matrix(np.zeros(shape=(n+1,n+1)))

    #Now add the dF_lamba/d(theta_i) = 1 and dF_lamba/d(lambda) = 0
    J[n,n] = 0

    #The following is correct for the off diagonal elements
    J[:n,:n] = -np.diag(np.array(theta)[0:n].flatten()) * self.W * np.diag(np.array(theta)[0:n].flatten())

    #We now update for the on diagonal elements
    J[:n,:n] = (J[:n,:n] - np.diag(np.diag(J[:n,:n])) +
           np.diag(np.array(-np.multiply(np.diag(np.diag(self.W)),np.diag(np.array(theta)[0:n].flatten())) * self.W * (theta[0:n] - self.m) + theta[n]).flatten()))

    #Finally adjust for the final columns
    J[:n,n] = theta[:n]
    J[n,:n] = 1

我不确定哪些numpy 调用的计算成本很高。是否可以在 Python 中对其进行优化以接近 C 的速度，还是我必须用 C 本身对其进行编程？

1 的测试功能

import numpy as np

def _nonLinEq(m, W, x, theta):
    #This outputs the nonlinear equations in theta
    #resulting from a the partial log derivative of a multivariate 
    #normal prior with covariance matrix E, means m and a multiinomial
    #likelihood with observations x.

    #F = [F_1, ... F_n, F_lambda]' ouput values where F_i = F(theta_i)
    n = len(m)
    F = define_F(n)
    F[0:n] = assign_values_to_F(x, theta, W, m, n)
    F[n] = assign_lagrange_multiplier_term(theta, n) #Lagrange multiplier term

    return F

def define_F(n):
    return np.matrix(np.zeros(shape=(n+1,1)))

def diag_theta(theta, n):
    return np.diag(np.array(theta)[0:n].flatten())

def multiply_terms(theta, W, m, n):
    return diag_theta(theta, n)*W*(theta[0:n]-m)

def assign_values_to_F(x,theta,W,m,n):
    return x - multiply_terms(theta, W, m, n) + theta[0:n]*theta[n]

def assign_lagrange_multiplier_term(theta, n):
    return np.sum(theta[0:n]) - 1 

def test_nonLinEq():
    n = 100
    temp = np.random.rand(n)
    m = np.transpose(np.matrix(temp/np.sum(temp)))
    W = np.matrix(np.diag(np.random.rand(n)))
    x = np.transpose(np.matrix(np.floor(np.random.rand(n)*10)))
    theta = np.transpose(np.matrix(np.append(np.ones(n)/n, -1)))

    for i in range(1000):
        _nonLinEq(m, W, x, theta)

2 的测试函数

def _jacNLE(m, W, x, theta):
    #This finds the Jacobian of our non-linear equations

    #J = (J_ij) ouput values where F_ij = dF_i/d(theta_j)
    n = len(m);

    J = define_J(n)

    #The following is correct for the off diagonal elements
    diag_theta = convert_theta_to_diagonal(theta, n)
    J[:n,:n] = input_off_diagonal_J(diag_theta, W)

    #We now update for the on diagonal elements
    J[:n,:n] = remove_J_diagonal(J, n) + new_diagonal(W, theta, m, diag_theta, n)

    #Finally adjust for the final columns
    J[:n,n] = theta[:n]
    J[n,:n] = 1

    return J

def define_J(n):
    return np.matrix(np.zeros(shape=(n+1,n+1)))

def convert_theta_to_diagonal(theta, n):
    return np.diag(np.array(theta)[0:n].ravel())

def input_off_diagonal_J(diag_theta, W):
    return -diag_theta * W * diag_theta

def remove_J_diagonal(J, n):
    return J[:n,:n] - np.diag(np.diag(J[:n,:n]))

def matrix_prod(W, diag_theta):
    return -np.multiply(np.diag(np.diag(W)),diag_theta)

def new_diagonal(W, theta, m, diag_theta, n):
    return np.diag(np.array(matrix_prod(W, diag_theta) * W * (theta[0:n] - m) + theta[n]).ravel())

def test_jacNLE():
    n = 2
    temp = np.random.rand(n)
    m = np.transpose(np.matrix(temp/np.sum(temp)))
    W = np.matrix(np.diag(np.random.rand(n)))
    x = np.transpose(np.matrix(np.floor(np.random.rand(n)*10)))
    theta = np.transpose(np.matrix(np.append(np.ones(n)/n, -1)))

    for i in range(1000):
        _jacNLE(m, W, x, theta)

【问题讨论】：

你有没有在你的代码上使用cProfile来查看哪个最贵？
在这种情况下，我总是首先将瓶颈与其余代码隔离开来。将这些中的每一个包装成一个您可以独立调用和分析的函数。将这些添加到一个测试脚本中，该脚本还为所有输入变量生成适当大小和数据类型的虚拟值。您的测试脚本应该调用每个函数，并传入这些虚拟值。现在，您可以对代码进行计时和分析，并衡量任何改进。最后，使用测试脚本更新您的问题。
我做了cProfile 我的代码，它指出了包含这两位代码的函数。我会按照 E 先生说的去做，努力找到问题的根源。
我已将我的首字母 cprofile 包含在此处登录 dumptext.com/hNaGFjnD。正如我所解释的那样，大部分时间都花在_nonLinEq 和_jacNLE 上，它们基本上分别由1 和2 中的代码组成。
@MrE，我已经上传了我认为你在案例 1 的评论中的意思。如果这实际上是你的意思，我也会为案例 2 做。

标签： python performance optimization numpy

【解决方案1】：

我最近想自己优化 numpy 代码。

起初我会尝试line_profiler package。它允许您测量每条线路的时间消耗。起初使用看起来有点棘手，但最终它是一个不错的工具，您可以肯定地检测到瓶颈。

另一个可以加快您的代码速度的工具是numba。如果您提供一些类型信息，它会即时编译您的代码。

我注意到 numpy 中有一些特定的函数很慢。

首先：np.sum 这很奇怪，但它甚至比 python 内置的 sum-function 还要慢。我曾经试图通过将要求和的数组与matrix(ones(N))* 相乘或编写一个函数“in”numba 来避免它。恐怕我不记得哪个解决方案更快了。

我想我读过一次，np.matrix 比简单的包装器 np.asmatrix 慢。

如果你仍然填充一个空数组，你可以使用np.empty 而不是np.zeros。

另一种方法是使用一些优化的 BLAS、LAPACK 东西再次在您的机器上编译 numpy（直到我遇到从未听说过的 numpy 问题...）它应该可以工作，但我失败了。

我还听说过使用 CUDA 和您的 GPU 功能的 gnumpy-Package，但我必须为此安装一些 CUDA 软件...

[编辑：] * 我预定义了onesN = matrix(ones(N))，所以每次我只想总结一个数组时都不必构建它。这实际上是另一个提示：查找常量并预定义并仔细阅读它们...

【讨论】：

line_profiler 包真的很有用。我也将实施 numba，看看是否能获得更多的加速。

【解决方案2】：

我同意 cmets，你应该做一些分析，但显然你可以做的一件事是分配

np.diag(np.array(theta)[0:n].flatten())

到一个局部变量并在你的计算中使用它，因为你做了几次。如果你的 n 很大，这将是一个相当昂贵的计算，因为你正在展平，复制到一个新数组，然后将其创建一个以该数组作为对角线的新矩阵。

【讨论】：