矩阵向量积的数值微分

Question

我有以下代码来使用公式近似函数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需要花费大量的时间。是否有更好的解决方案来构造一个三对角矩阵？

分析这两种功能的结果如下：

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上

Answer 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]的两个三重态分开的。因此，作为一个黑客，我们可以创建一个普通的2D数组，将复制的三胞胎作为前三列，如下所示-

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))

Answer 2

没有必要构造矩阵。您可以直接使用向量f。例如，下面的版本就可以了。

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

在重复导数计算的情况下，矩阵方法是有用的。您存储矩阵并每次重用它。对于更高的精度来说，它变得更加有用。