Python中多变量函数的向量化偏导数

Question

在这个讨论中，alko 提供了一个非常棒的答案，讲的是如何用数字方法计算多变量函数的偏导数。

现在我有个后续问题，想要改进这个函数，让它可以接受一组输入值。我有一些代码，正在遍历一个很长的 n 维点的列表，计算每个变量的偏导数，这样做计算起来非常耗时。

用 np.vectorize 来向量化这个函数其实很简单，但这会导致 partial_derivative 包装器出现问题：

from scipy.misc import derivative
import numpy as np

def foo(x, y):
    return(x**2 + y**3)

def partial_derivative(func, var=0, point=[]):
    args = point[:]
    def wraps(x):
        args[var] = x
        return func(*args)
    return derivative(wraps, point[var], dx=1e-6)

vfoo = np.vectorize(foo)

>>>foo(3,1)
>>>10
>>>vfoo([3,3], [1,1])
>>>array([10,10])

>>>partial_derivative(foo,0,[3,1])
>>>6.0
>>>partial_derivative(vfoo,0,[[3,3], [1,1]])
>>>TypeError: can only concatenate list (not "float") to list

理想情况下，最后一行应该返回 [6.0, 6.0]。在这个例子中，传给向量化函数 vfoo 的两个数组实际上是成对组合的，所以 ([3,3], [1,1]) 被转化成了两个点，[3,1] 和 [3,1]。当这些点传给 wraps 函数时，似乎出现了问题。最终传给 derivative 函数的点是 [3,3]。此外，明显还有一个 TypeError 的错误。

有没有人有什么建议或者推荐的解决方案？有没有人遇到过类似的需求？

编辑

有时候我觉得在 Stack Overflow 上发帖就是打破思维障碍的关键。我想我已经找到了一个解决方案，供有兴趣的人参考：

vfoo = np.vectorize(foo)

foo(3,1)
X = np.array([3,3])
Y = np.array([1,1])
vfoo(X, Y)
partial_derivative(foo,0,[3,1])
partial_derivative(vfoo,0,[X, Y])

现在最后一行返回的是 array([ 6., 6.])

Answer 1

我遇到了一个小问题，就是当你用 args[var] = x 这样的方式时，这可能会永久性地改变 args[var] 的值，而且所有的值都是通过引用传递的，不管你改得多小。因此，你可能得不到你想要的准确答案。下面是一个例子：

In[67]: a = np.arange(9).reshape(3,3)
In[68]: b = a[:]
In[69]: b[0,0]=42
In[70]: a
Out[70]: 
array([[42,  1,  2],
       [ 3,  4,  5],
       [ 6,  7,  8]])

你需要通过例如以下方式来修复这个问题：

def wraps(x):
    tmp = args[var]
    args[var] = x
    ret= func(*args)
    args[var] = tmp
    return ret

另外，你可以使用 numdifftools。他们似乎知道自己在做什么。这个工具可以处理所有的偏导数：

import numpy as np
import numdifftools as nd

def partial_function(f___,input,pos,value):
    tmp  = input[pos]
    input[pos] = value
    ret = f___(*input)
    input[pos] = tmp
    return ret

def partial_derivative(f,input):
    ret = np.empty(len(input))
    for i in range(len(input)):
        fg = lambda x:partial_function(f,input,i,x)
        ret[i] = nd.Derivative(fg)(input[i])
    return ret
if __name__ == "__main__":
    f     = lambda x,y: x*x*x+y*y
    input = np.array([1.0,1.0])
    print ('partial_derivative of f() at: '+str(input))
    print (partial_derivative(f,input))

最后：如果你想让你的函数接受一个参数数组，比如：

f     = lambda x: x[0]*x[0]*x[0]+x[1]*x[1]

那么就把相应的那一行替换成（去掉了 '*'）

ret = f___(input)