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

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

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型函数求导，

$$\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:\, \mathbb{R}^n\rightarrow\mathbb{R}^m$是一个从欧氏$n$维空间转换到$m$维欧氏空间的函数（如$y_{m\times 1}=A_{m\times n}x_{n\times 1}$），这个函数由$m$个实函数组成，$y_1(x_1, \ldots,x_n),\ldots,y_m(x_1,\ldots,x_n)$（刚好又可看做一个 $m$ 维的向量，每一个 entry 又都是一个 $n$ 元函数），这些函数的偏导数如果存在可以组成一个$m\times n$的矩阵，这即是所谓的雅可比矩阵：

$$\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}$$

此矩阵表示为：$J_F(x_1,\ldots,x_n)$，或者$\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.]])

$$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}$$

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

$$\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(x_1, x_2, \ldots, x_n)$$
如果 $f$的所有二阶导数都存在，那么$f$的 Hessian 矩阵为：
$$H(f)_{ij}(x)=D_iD_jf(x)$$
其中 $x=(x_1,x_2,\ldots,x_n)$，即 $H(f)$为：

$$\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)= \begin{bmatrix} 2&&0\\ 0&&2 \end{bmatrix}$$

References

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