1. 计算梯度
创建一个函数
y
,并且计算关于其参数
以下就是计算该梯度的 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)=11+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
的第二个参数可以是一个列表 (list),此时输出也为一个列表 . 不管是参数还是输出,列表的顺序是很重要的,输出列表中的第
i
个元素,就是
2. 计算Jacobian
在 Theano 中,Jacobian 表示输出函数关于其输入的一阶偏导数(在数学中,就称为雅可比矩阵). Theano 中使用
theano.gradient.jacobian()
来计算需要的 Jacobian . 接下来说明怎样手动执行.
我们使用
scan
来手动计算函数
y
关于参数
补充:
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
到
补充:关于
T.grad
的使用,有一些陷阱. 例如,不能将代码中的 Jacobian 表达式改为
theano.scan(lambdayi,x:T.grad(yi,x),sequences=y,nonsequences=x)
,即使从
scan
的文档上看似乎是可以的. 不能改写的原因在于
yi
不再是关于
x
的函数,而
3. 计算 Hessian
在 Theano 中,术语 Hessian 和数学中的概念一样:是标量输出函数关于向量输入的二阶偏微分矩阵. Theano 中使用
theano.gradient.hessian()
函数计算 Hessian . 接下来解释怎样手动执行.
手动计算 Hessian 和手动计算 Jacobian 类似,唯一的不同就是用
T.grad(cost,x)
代替 Jacobian 中的函数
y
,
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 和向量之间的乘积,写作:
∂f(x)∂xv
. 这个公式能够被扩展,即使
x
是一个矩阵或者张量,此时 Jacobian 成为了一个张量而其乘积成为了张量的乘积. 因此在实际中,我们需要根据权重矩阵计算这些表达式, Theano 支持这种更通用的表示形式. 为了评估表达式
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∂f(x)∂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 而言,
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 和向量乘以 Hessian5
正在优化有效的计算完整的 Jacobian 和 Hessian 矩阵以及 Jacobian 乘以向量.