使用joschu/cgt项目进行符号计算与自动微分教程

什么是joschu/cgt项目

joschu/cgt是一个用于符号计算和自动微分的Python库，它允许用户构建数学表达式图，并自动计算这些表达式的梯度。该库特别适合机器学习中的优化问题，能够高效地处理向量、矩阵运算，并支持GPU加速。

基本工作流程

1. 定义符号变量

首先需要创建符号变量，这些变量代表计算图中的输入节点：

import cgt
a = cgt.scalar(name='a') # 浮点型标量，可选名称参数
b = cgt.scalar(name='b')
n = cgt.scalar(name='n', dtype='int64') # 整型标量

这里创建了三个符号变量：两个浮点型标量a和b，一个整型标量n。这些变量将成为我们计算图的输入节点。

2. 构建表达式

利用这些变量构建更复杂的数学表达式：

c = (a**n + b**n)**(1.0/n)

CGT会自动重载基本运算符(如+、<、等)，使得我们可以像操作普通数值一样操作这些符号变量。在后台，CGT会构建一个表示这个计算过程的有向无环图(DAG)**。

3. 编译为可执行函数

将符号表达式编译为可执行函数：

f = cgt.function([a,b,n], c)
print f(8,15,2)  # 输出17.0

这个函数可以接受数值输入并返回计算结果。例如，当输入(8,15,2)时，计算的是(8² + 15²)的平方根，即17.0。

可视化计算图

CGT的一个强大功能是可以可视化计算图：

cgt.as_dot(c)

这将生成计算图的可视化表示，其中节点表示操作，边表示数据依赖关系。标签0和1表示函数参数的索引位置。

线性回归示例

让我们通过一个线性回归的例子展示更复杂的用法：

定义符号变量

X_nk = cgt.matrix("X")  # 输入特征矩阵
y_n = cgt.vector("y")    # 目标值向量
w_k = cgt.vector("w")    # 权重向量
b = cgt.scalar("b")      # 偏置项

构建预测和损失函数

ypred_n = X_nk.dot(w_k) + b  # 线性预测
L = cgt.sum(cgt.square(ypred_n - y_n))  # 平方损失
print "L = ", cgt.print_expr(L)

这里我们构建了一个标准的线性回归模型，使用平方误差作为损失函数。

计算梯度

grads = dLdw, dLdb = cgt.grad(L, [w_k, b])

CGT使用反向模式自动微分(即反向传播算法)自动计算损失函数对参数的梯度。这些梯度本身也是符号表达式。

编译为可执行函数

f_loss = cgt.function([X_nk, y_n, w_k, b], L)
f_grads = cgt.function([X_nk, y_n, w_k, b], grads)

使用数值数据进行测试

import numpy as np, numpy.random as nr
nr.seed(0)

# 生成随机数据
Xval = nr.randn(100,3)
yval = nr.randn(100)

# 初始化参数
wval = nr.randn(3)
bval = nr.randn()

# 计算损失和梯度
print "loss:", f_loss(Xval, yval, wval, bval)
print "grad:", f_grads(Xval, yval, wval, bval)

优化参数

方法1：使用SciPy优化器

def f(theta):
    return f_loss(Xval, yval, theta[0:3], theta[3])

def gradf(theta):
    gw,gb = f_grads(Xval, yval, theta[0:3], theta[3])
    return np.concatenate([gw,gb.reshape(1)])

import scipy.optimize as opt
theta_opt = opt.fmin_bfgs(f, np.zeros(4), gradf, disp=False)
print "Optimal theta:", theta_opt

方法2：使用CGT内置的优化

theta = cgt.vector('theta')
w_k_1 = theta[0:3]
b_1 = theta[3]
ypred_n_1 = X_nk.dot(w_k_1) + b_1
L_1 = cgt.sum(cgt.square(ypred_n_1 - y_n))
dLdtheta, = cgt.grad(L_1, [theta])

f = cgt.function([theta], L_1, givens=[(X_nk,Xval), (y_n,yval)])
gradf = cgt.function([theta], dLdtheta, givens=[(X_nk,Xval), (y_n,yval)])

theta_opt = opt.fmin_bfgs(f, np.zeros(4), gradf, disp=False)
print "Optimal theta:", theta_opt

方法3：使用共享变量和原地更新

这种方法特别适合大规模问题和GPU计算：

stepsize=0.001
w_k_2 = cgt.shared(wval.copy(), name='w')
b_2 = cgt.shared(bval, name='b')
ypred_n_2 = X_nk.dot(w_k_2) + b_2
L_2 = cgt.sum(cgt.square(ypred_n_2 - y_n))
dLdw_2,dLdb_2 = cgt.grad(L_2, [w_k_2, b_2])

updates = [(w_k_2, w_k_2 - stepsize*dLdw_2), (b_2, b_2 - stepsize*dLdb_2)]
givens = [(X_nk,Xval), (y_n,yval)]
f_update = cgt.function([], L_2, givens=givens, updates=updates)

for i in xrange(100):
    f_update()
print w_k_2.op.get_value(), b_2.op.get_value()

总结

joschu/cgt项目提供了一个强大的符号计算和自动微分框架，特别适合机器学习中的优化问题。通过本教程，我们学习了：

1. 如何定义符号变量和构建表达式
2. 如何可视化计算图
3. 如何自动计算梯度
4. 如何将符号表达式编译为可执行函数
5. 三种不同的参数优化方法

CGT的设计使得从简单的数学表达式到复杂的机器学习模型都能高效地实现和优化。其支持GPU加速的特性使其特别适合大规模机器学习任务。