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 乘以向量.