Theano学习-梯度计算

1. 计算梯度

创建一个函数 \(y\) ，并且计算关于其参数 \(x\) 的微分. 为了实现这一功能，将使用函数 \(T.grad\) . 例如：计算 \(x^2\) 关于参数 \(x\) 的梯度. 注：$ d(x^2)/d(x) = 2 * x$.
以下就是计算该梯度的 Python 代码：

import numpy
import theano
import theano.tensor as T
from theano import pp
x = T.dscalar('x')
y = x ** 2
gy = T.grad(y, x)
pp(gy)
f = theano.function([x], gy)
print(f(4))  #print 8
print(numpy.allclose(f(94.2), 188.4))  #print True
print(pp(f.maker.fgraph.outputs[0])) #print (TensorConstant{2.0} * x)

我们能从 \(pp(gy)\) 看出计算出的符号梯度是正确的.
同时也能计算复杂表达式的梯度，例如 \(Logistic\) 函数的梯度.
\(Logistic\) 函数: \(s(x) = \frac{1}{1 + e^{-x}}\)
其导数： \(ds(x)/d(x) = s(x) * (1 - s(x))\)

import numpy
import theano
import theano.tensor as T
x = T.dmatrix('x')
s = T.sum(1 / (1 + T.exp(-x)))
gs = T.grad(s, x)
dlogistic = theano.function([x], gs)
print(dlogistic([[0, 1], [-1, -2]])) 
#print [[ 0.25  0.19661193]
        [ 0.19661193  0.10499359]]

总的来说，对于标量表达式 \(s\) ， Theano 提供的方法 \(T.grad(s, w)\) 即为计算 \(\frac{\partial s}{\partial w}\).
补充： \(T.grad\) 的第二个参数可以是一个列表 (list)，此时输出也为一个列表 . 不管是参数还是输出，列表的顺序是很重要的，输出列表中的第 \(i\) 个元素，就是 \(T.grad\) 中第一个参数关于第二个列表参数中第 \(i\) 个元素的导数. \(T.grad\) 的第一个参数必须是个标量.

2. 计算Jacobian

在 Theano 中，Jacobian 表示输出函数关于其输入的一阶偏导数(在数学中，就称为雅可比矩阵). Theano 中使用 \(theano.gradient.jacobian()\) 来计算需要的 Jacobian . 接下来说明怎样手动执行.
我们使用 \(scan\) 来手动计算函数 \(y\) 关于参数 \(x\) 的 Jacobian ，就是循环 \(y\) 的元素计算 \(y[i]\) 关于参数 \(x\) 的梯度.
补充： \(scan\) 在 Theano 中是一个通用操作， 允许以符号方式写入循环方程. 创建符号循环是一个艰巨的任务，目前正在努力提高 \(scan\) 的性能.

import theano
import theano.tensor as T
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)
print(f([4, 4])) 
#print [[ 8.  0.]
        [ 0.  8.]]

在此代码中，使用 \(T.arange\) 产生从 \(0\) 到 \(y.shape[0]\) 的 \(int\) 序列. 当循环序列时，计算元素 \(y[i]\) 关于参数 \(x\) 的梯度. \(scan\) 自动连接所有的行，生成一个对应于 Jacobian 的矩阵.
补充：关于 \(T.grad\) 的使用，有一些陷阱. 例如，不能将代码中的 Jacobian 表达式改为 \(theano.scan(lambda y_i, x: T.grad(y_i, x), sequences = y, non_sequences = x )\) ，即使从 \(scan\) 的文档上看似乎是可以的. 不能改写的原因在于 \(y_i\) 不再是关于 \(x\) 的函数，而 \(y[i]\) 仍然是.

3. 计算 Hessian

在 Theano 中，术语 Hessian 和数学中的概念一样：是标量输出函数关于向量输入的二阶偏微分矩阵. Theano 中使用 \(theano.gradient.hessian()\) 函数计算 Hessian . 接下来解释怎样手动执行.
手动计算 Hessian 和手动计算 Jacobian 类似，唯一的不同就是用 \(T.grad(cost, x)\) 代替 Jacobian 中的函数 \(y\), \(cost\) 通常是一个标量.

import theano
import theano.tensor as T
x = T.dvector('x')
y = x ** 2
cost = y.sum()
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)
print(f([4, 4]))
#print [[ 2.  0.]
        [ 0.  2.]]

4. Jacobian 乘以向量

在解释算法的过程中，有时我们需要表示 Jacobian 乘以向量，或者是向量乘以 Jacobian . 与先评估 Jacobian 再做乘积相比，目前有许多方式能直接计算所需结果从而避免实际的 Jacobian 评估. 这能带来显著的性能提升，具体参考下文：
[Barak A. Pearlmutter, “Fast Exact Multiplication by the Hessian”, Neural Computation, 1994]
原则上我们希望 Theano 能为我们自动识别这些模式，实际上，以通用的方式进行优化是很困难的.
因此，提供了特殊的函数专用于这些任务.

R-operator

R-operator 用于评估 Jacobian 和向量之间的乘积，写作： \(\frac{\partial f(x)}{\partial x}v\). 这个公式能够被扩展，即使 \(x\) 是一个矩阵或者张量，此时 Jacobian 成为了一个张量而其乘积成为了张量的乘积. 因此在实际中，我们需要根据权重矩阵计算这些表达式， Theano 支持这种更通用的表示形式. 为了评估表达式 \(y\) 关于参数 \(x\) 的 R-operator，将 Jacobian 与 \(v\) 右乘，你需要做下面的事：

import theano
import theano.tensor as T
W = T.dmatrix('W')
V = T.dmatrix('V')
x = T.dvector('x')
y = T.dot(x, W)
JV = T.Rop(y, W, V)
f = theano.function([W, V, x], JV)
print(f([[1, 1], [1, 1]], [[2, 2], [2, 2]], [0, 1]))
#print [ 2.  2.]

L-operator

与 R-operator 类似， L-operator是计算行向量与 Jacobian 的乘积. 数学形式为： $ v\frac{\partial f(x)}{\partial x}$ . 同样的，可以通过下面的程序执行：

import theano
import theano.tensor as T
W = T.dmatrix('W')
V = T.dmatrix('V')
x = T.dvector('x')
y = T.dot(x, W)
VJ = T.Lop(y, W, V)
f = theano.function([V, x], VJ)
print(f([2, 2], [0, 1]))
#print [[ 0.  0.]
        [ 2.  2.]]

补充：\(v\) 的评估，L-operator 与 R-operator 是不一样的. 对 L-operator 而言， \(v\) 需要和输出有同样的形状；而对 R-operator 需要和输出参数有同样的形状. 进一步说，两种运算符的结果是不一样的. L-operator 的结果与输入参数形状相同；而 R-operator 的结果具有与输出类似的形状.

5. Hessian 乘以向量

假如你需要计算 Hessian 乘以一个向量，你可以使用上面定义的算子直接计算，这比先计算精确的 Hessian ，再计算乘积更有效. 由于 Hessian 矩阵的对称性，你有两种方式得到同样的结果，尽管这两种方式可能展现出不同的性能. 下面给出这两种方式：

• 1

import theano
import theano.tensor as T
x = T.dvector('x')
v = T.dvector('v')
y = T.sum(x ** 2)
gy = T.grad(y, x)
vH = T.grad(T.sum(gy * v), x)
f = theano.function([x, v], vH)
print(f([4, 4], [2, 2]))
#print [ 4.  4.]

• 2使用 R-operator

import theano
import theano.tensor as T
x = T.dvector('x')
v = T.dvector('v')
y = T.sum(x ** 2)
gy = T.grad(y, x)
Hv = T.Rop(gy, x, v)
f = theano.function([x, v], Hv)
print(f([4, 4], [2, 2]))
#print [ 4.  4.]

6. 指示

• 1
  \(grad\) 函数以符号方式工作：接收与返回都是 Theano 变量

• 2
  \(grad\) 可以比作宏，因为它可以重复使用

• 3
  标量函数只能被 \(grad\) 直接处理，矩阵能够通过重复应用来处理

• 4
  内置函数能有效的计算向量乘以 Jacobian 和向量乘以 Hessian

• 5
  正在优化有效的计算完整的 Jacobian 和 Hessian 矩阵以及 Jacobian 乘以向量.