从 Jacobian 矩阵、Hessian 矩阵到 Theano 实现

T.grad(cost, wrt)，一般接收两个参数，第一个参数表示需要求导的函数，放在深度学习的背景下就是代价函数，wrt（with respect to）表示代价函数所关于的参数（通俗地讲，就叫自变量，f(x)f(x)表示关于xx的函数ff）。

T.grad的第一个参数必须是标量。

>>> import thenao
>>> import theano.tensor as T
>>> x = T.dmatrix('x')
>>> y = x**2+x
>>> gy = T.grad(y, x)TypeError: cost must be a scalar.

>>> x = T.dmatrix('x')
>>> y = T.sum(x**2+x)# 这里的T.sum似乎并没有做加和的动作
>>> gy = T.grad(y, x)
>>> f = theano.function([x], gy)
>>> f([[0, 1], [2, 3]])
array([[ 1.,  3.],[ 5.,  7.]])

同理，对sigmoid型函数求导，

ds(x)dx=s(x)(1−s(x))

\frac{d\,s(x)}{dx}=s(x)(1-s(x))

>>> x = T.dmatrix('x')
>>> s = T.sum(1./(1.+T.exp(-x)))
>>> gs = T.grad(s, x)
>>> dlogistic = theano.function([x], gs)
>>> dlogistic([[0, 1], [-1, -2]])
array([[ 0.25      ,  0.19661193],[ 0.19661193,  0.10499359]])

1. Jacobian 矩阵

在向量分析中，雅克比矩阵是一阶偏导数以一定方式排列成的矩阵，其行列式称为雅克比行列式。
假设F:Rn→RmF:\, \mathbb{R}^n\rightarrow\mathbb{R}^m是一个从欧氏nn维空间转换到mm维欧氏空间的函数（如ym×1=Am×nxn×1y_{m\times 1}=A_{m\times n}x_{n\times 1}），这个函数由mm个实函数组成，y1(x1,…,xn),…,ym(x1,…,xn)y_1(x_1, \ldots,x_n),\ldots,y_m(x_1,\ldots,x_n)（刚好又可看做一个 mm 维的向量，每一个 entry 又都是一个 nn 元函数），这些函数的偏导数如果存在可以组成一个m×nm\times n的矩阵，这即是所谓的雅可比矩阵：

⎡⎣⎢⎢⎢⎢⎢⎢⎢∂y1∂x1⋮∂ym∂x1⋯⋱⋯∂y1∂xn⋮∂ym∂xn⎤⎦⎥⎥⎥⎥⎥⎥⎥

\begin{bmatrix} \frac{\partial y_1}{\partial x_1}&&\cdots&&\frac{\partial y_1}{\partial x_n}\\ \vdots&&\ddots&&\vdots\\ \frac{\partial y_m}{\partial x_1}&&\cdots&&\frac{\partial y_m}{\partial x_n} \end{bmatrix}

此矩阵表示为：JF(x1,…,xn)J_F(x_1,\ldots,x_n)，或者∂(y1,…,ym)∂(x1,…,xn)\frac{\partial(y_1,\ldots,y_m)}{\partial(x_1,\ldots,x_n)}

>>> x = T.dvector('x')
>>> y = x**2
>>> J, updates = theano.scan(lambda i, y, x: T.grad(y[i], x), sequences=T.arange(y.shape[0]), non_sequences=[y, x])
>>> f = theano.function([x], J, updates=updates)
>>> f([3, 4])
array([[ 6.,  0.],[ 0.,  8.]])

Jy(x1,x2)=∂(y1,y2)∂(x1,x2)=⎡⎣⎢⎢⎢∂y1∂x1∂y2∂x1∂y1∂x2∂y2∂x2⎤⎦⎥⎥⎥=[2x1002x2]

J_y(x_1, x_2)=\frac{\partial (y_1,y_2)}{\partial (x_1, x_2)}=\begin{bmatrix}\frac{\partial y_1}{\partial x_1}&&\frac{\partial y_1}{\partial x_2}\\ \frac{\partial y_2}{\partial x_1} &&\frac{\partial y_2}{\partial x_2}\end{bmatrix}=\begin{bmatrix}2x_1&&0 \\ 0 && 2x_2\end{bmatrix}

再考虑这样一个向量，f(x)=[3x21+x2ln(x1)sin(x2)]T\mathrm f(\mathrm x)=\left[3x_1^2+x_2 \quad \ln(x_1) \quad \sin(x_2) \right]^T，则 Jacobian 为：

∂f∂x=⎛⎝⎜⎜⎜6x11x1010cos(x2)⎞⎠⎟⎟⎟

\frac{\partial \mathbf f}{\partial \mathbf x}=\begin{pmatrix} 6x_1&1\\ \frac1{x_1}&0\\ 0&\cos(x_2) \end{pmatrix}

2. Hessian矩阵

Hessian matrix是一个自变量为向量的实值函数的二阶偏导数组成的方块矩阵，此函数如下：

f(x1,x2,…,xn)

f(x_1, x_2, \ldots, x_n)
如果 ff的所有二阶导数都存在，那么ff的 Hessian 矩阵为：

H(f)ij(x)=DiDjf(x)

H(f)_{ij}(x)=D_iD_jf(x)
其中 x=(x1,x2,…,xn)x=(x_1,x_2,\ldots,x_n)，即 H(f)H(f)为：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢∂2f∂x21∂2f∂x2x1⋮∂2f∂xnx1∂2f∂x1x2∂2f∂x22⋮∂2f∂xnx2⋯⋯⋱⋯∂2f∂x1xn∂2f∂x2xn⋮∂2f∂x2n⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

\begin{bmatrix} \frac{\partial^2f}{\partial x_1^2}&&\frac{\partial^2f}{\partial x_1x_2} &&\cdots&&\frac{\partial^2f}{\partial x_1x_n}\\ \frac{\partial^2f}{\partial x_2x_1}&&\frac{\partial^2f}{\partial x_2^2} &&\cdots&&\frac{\partial^2f}{\partial x_2x_n}\\ \vdots && \vdots && \ddots && \vdots\\ \frac{\partial^2f}{\partial x_nx_1}&&\frac{\partial^2f}{\partial x_nx_2} &&\cdots&&\frac{\partial^2f}{\partial x_n^2} \end{bmatrix}

>>> x = T.dvector('x')
>>> y = x**2
>>> cost = T.sum(y)
>>> gy = T.grad(cost, x)# 第一次求导
>>> H, updates = theano.scan(lambda i, gy, x: T.grad(gy[i], x), sequences=T.arange(gy.shape[0]), non_sequences=[gy, x])# 第二次求导
>>> f = theano.function([x], H, updates=updates)
>>> f([3, 4])
array([[ 2.,  0.],[ 0.,  2.]])

H(f)ij(x)=[2002]

H(f)_{ij}(x)= \begin{bmatrix} 2&&0\\ 0&&2 \end{bmatrix}

References

[1] Jacobian矩阵和Hessian矩阵
[2] Derivatives in Theano