作为矩阵向量积的数值微分答案

【问题标题】：Numerical Differenciation as Matrix-Vector Product作为矩阵向量积的数值微分
【发布时间】：2015-06-11 16:29:44
【问题描述】：

我有以下代码使用公式逼近函数f() 的二阶导数：

我想比较两种不同的方法；使用循环和矩阵向量乘积，希望表明numpy 版本更快：

def get_derivative_loop(X):                             
    DDF = []
    for i in range(1,len(X)-1):
         DDF.append((f(X[i-1]) - 2*f(X[i]) + f(X[i+1]))/(h**2))
    return DDF
def get_derivative_matrix(X):
    A = (np.diag(np.ones(m)) + 
         np.diag(-2*np.ones(m-1), 1) +  
         np.diag(np.ones(m-2), 2))/(h**2)
    return np.dot(A[0:m-2], f(X))

正如预期的那样，构建矩阵A 会消耗大量时间。在numpy中构造三对角矩阵是否有更好的解决方案？

分析这两个函数会产生：

Total time: 0.003942 s
File: diff.py
Function: get_derivative_matrix at line 17

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    17                                           @profile
    18                                           def get_derivative_matrix(X):
    19         1         3584   3584.0     90.9      A = (np.diag(np.ones(m)) + np.diag(-2*np.ones(m-1), 1) + np.diag(np.ones(m-2), 2))/(h**2)
    20         1          358    358.0      9.1      return np.dot(A[0:m-2], f(X))

Total time: 0.004111 s
File: diff.py
Function: get_derivative_loop at line 22

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    22                                           @profile
    23                                           def get_derivative_loop(X):
    24         1            1      1.0      0.0      DDF = []
    25       499          188      0.4      4.6      for i in range(1, len(X)-1):
    26       498         3921      7.9     95.4          DDF.append((f(X[i-1]) - 2*f(X[i]) + f(X[i+1]))/(h**2))
    27                                           
    28         1            1      1.0      0.0      return DDF
    A = (np.diag(np.ones(m)) +
         np.diag(-2*np.ones(m-1), 1) + 
         np.diag(np.ones(m-2), 2))/(h**2)
    return np.dot(A[0:m-2], f(X))

编辑

虽然初始化只进行一次是正确的，所以不需要优化，但是我发现想出一种很好且快速的方法来设置该矩阵是很有趣的。

这是使用Divakar 方法的配置文件结果

Timer unit: 1e-06 s

Total time: 0.006923 s
File: diff.py
Function: get_derivative_matrix_divakar at line 19

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    19                                           @profile
    20                                           def get_derivative_matrix_divakar(X):
    21                                           
    22                                               # Setup output array, equivalent to A
    23         1           48     48.0      0.7      out = np.zeros((m, 3+m-2))
    24                                           
    25                                               # Setup the triplets in each row as [1,-2,1]
    26         1         1485   1485.0     21.5      out[:, 0:3] = 1
    27         1           22     22.0      0.3      out[:, 1] = -2
    28                                           
    29                                               # Slice and perform matrix-multiplication
    30         1         5368   5368.0     77.5      return np.dot(out.ravel()[:m*(m-2)].reshape(m-2, -1)/(h**2), f(X))


Total time: 0.019717 s
File: diff.py
Function: get_derivative_matrix at line 45

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    45                                           @profile
    46                                           def get_derivative_matrix(X):
    47         1        18813  18813.0     95.4      A = (np.diag(np.ones(m)) + np.diag(-2*np.ones(m-1), 1) + np.diag(np.ones(m-2), 2))/(h**2)
    48         1          904    904.0      4.6      return np.dot(A[0:m-2], f(X))



Total time: 0.000108 s
File: diff.py
Function: get_derivative_slice at line 41

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
    41                                           @profile
    42                                           def get_derivative_slice(X):
    43         1          108    108.0    100.0      return (f(X[0:-2]) - 2*f(X[1:-1]) + f(X[2:]))/(h**2)

新方法更快。但是，我不明白为什么 21.5% 花费在这个初始化上 out[:, 0:3] = 1

【问题讨论】：

标签： python numpy

【解决方案1】：

对于m = 9，没有h 缩放的三对角矩阵将如下所示 -

array([[ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  1., -2.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  1., -2.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  1., -2.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  1., -2.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  1., -2.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  1., -2.,  1.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  1., -2.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  1.]])

可以看到，按行排列，7 (=m-2) 个零将[1,-2,1] 的两个三元组分开。因此，作为一种 hack，可以创建一个常规的二维数组，前三列是复制的三元组，像这样 -

array([[ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1., -2.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

上述矩阵创建的优点是易于索引，这必须非常有效。所以，为了得到我们想要的输出，剩下的工作就是切片以限制我们使用 m**2 元素并在最后处理三元组。

最后，我们会有这样的东西来得到三对角矩阵 -

def three_diag_mat(m,h):
    # Initialize output array
    out = np.zeros((m,3+m-2))

    # Setup the triplets in each row as [1,-2,1]
    out[:,:3] = 1
    out[:,1] = -2

    # Reset the ending "1" of the second last row as zero.
    out[m-2,2] = 0

    # Slice upto m**2 elements in a flattened version.
    # Then, scale down the sliced output by h**2 for the final output. 
    return (out.ravel()[:m**2].reshape(m,m))/(h**2)

运行时测试和验证结果

案例#1：

In [8]: m = 100; h = 10

In [9]: %timeit (np.diag(np.ones(m)) + 
np.diag(-2*np.ones(m-1), 1) + np.diag(np.ones(m-2), 2))/(h**2)
10000 loops, best of 3: 119 µs per loop

In [10]: %timeit three_diag_mat(m,h)
10000 loops, best of 3: 51.8 µs per loop

In [11]: np.array_equal((np.diag(np.ones(m)) + np.diag(-2*np.ones(m-1), 1) +
 np.diag(np.ones(m-2), 2))/(h**2),three_diag_mat(m,h))
Out[11]: True

案例#2：

In [12]: m = 1000; h = 10

In [13]: %timeit (np.diag(np.ones(m)) + 
np.diag(-2*np.ones(m-1), 1) + np.diag(np.ones(m-2), 2))/(h**2)
100 loops, best of 3: 16.2 ms per loop

In [14]: %timeit three_diag_mat(m,h)
100 loops, best of 3: 5.66 ms per loop

In [15]: np.array_equal((np.diag(np.ones(m)) + np.diag(-2*np.ones(m-1), 1) + 
np.diag(np.ones(m-2), 2))/(h**2),three_diag_mat(m,h))

具体用例：对于使用A[0:m-2]的用例，您可以避免少量计算并修改get_derivative_matrix，如下所示：

def get_derivative_matrix(X):

    # Setup output array, equivalent to A
    out = np.zeros((m,3+m-2))

    # Setup the triplets in each row as [1,-2,1]
    out[:,:3] = 1
    out[:,1] = -2

    # Slice and perform matrix-multiplication
    return np.dot(out.ravel()[:m*(m-2)].reshape(m-2,-1)/(h**2), f(X))

【讨论】：

看起来“丑陋” :-) 但肯定更快。但是，这会使标量积最终变得毫无用处，对吧？
@Tengis 没错！因此，您将得到类似于最后编辑代码中显示的内容。一探究竟！让我知道你得到了什么样的加速？
@Tengis 啊，这有点令人惊讶。我认为为了提高效率，您可以执行两个单独的索引来替换out[:, 0:3] = 1：out[:, 0] = 1; out[:, 2] = 1，保持其余代码相同。
我之前尝试过，因为我认为这样会减少冗余。不幸的是，并没有使整体功能更快。无论如何，我认为你的解决方案是最好的，初始化不应该是分析的一部分，因为它只完成一次。
out[:, 2] = 1 吃掉 22.4%，而其他两个任务分别吃掉 0.4% 和 0.1%。这很奇怪，不是吗？

【解决方案2】：

不需要构造矩阵。您可以直接使用向量 f。例如下个版本就好了

def get_derivative(x,f,h):
    fx=f(x)
return (fx[:-2]-2*fx[1:-1]+fx[2:])/h**2

矩阵方法在重复导数计算的情况下很有用。您存储矩阵并每次都重复使用它。它对于更高的阶精度变得更有用。

【讨论】：

我不想使用切片来解决这个问题，即使它更快（~100）。我只是认为必须有一种方法可以构造该矩阵而无需使用三次diag()
可能有一种更快的方法来设置矩阵，但是根据我的经验，如果您只设置一次矩阵而不是在每次函数调用时，初始化从来都不是我计算的瓶颈，所以这似乎是一个过早的优化。