Theano 学习

1. 计算梯度

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

import numpyimport theanoimport theano.tensor as Tfrom theano import ppx = T.dscalar('x')y = x ** 2gy = T.grad(y, x)pp(gy)f = theano.function([x], gy)print(f(4))  #print 8print(numpy.allclose(f(94.2), 188.4))  #print Trueprint(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 numpyimport theanoimport theano.tensor as Tx = 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) 即为计算 ∂s∂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 theanoimport theano.tensor as Tx = T.dvector('x')y = x ** 2J, 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(lambdayi,x:T.grad(yi,x),sequences=y,nonsequences=x) ，即使从 scan 的文档上看似乎是可以的. 不能改写的原因在于 yi 不再是关于 x 的函数，而 y[i] 仍然是.

3. 计算 Hessian

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

import theanoimport theano.tensor as Tx = T.dvector('x')y = x ** 2cost = 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 支持这种更通用的表示形式. 为了评估表达式 y 关于参数 x 的 R-operator，将 Jacobian 与 v 右乘，你需要做下面的事：

import theanoimport theano.tensor as TW = 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 theanoimport theano.tensor as TW = 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 theanoimport theano.tensor as Tx = 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 theanoimport theano.tensor as Tx = 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 乘以向量.