Python 实现自动微分

用了差不多一天多的时间,把自动微分基本整明白了,并且实现(抄)了一遍。
李理大神的这篇博客讲的非常清楚明白,即使我这样的小白也能看明白:http://fancyerii.github.io/books/autodiff/ ,理论讲解都在这篇里了,这里不再赘述。
代码的话是参考这篇:https://zhuanlan.zhihu.com/p/161635270 ,找了 github 上的一些自动微分的实现,都是比较成系统的集成工具,东西很多很全面,一时半会儿也看不明白,只会感觉自己智商受到了碾压,只有这篇的实现比较适合小白练手。

自动微分其实就是利用链式求导法则,构建计算图,将函数拆解为原子操作,然后利用图算法进行梯度累加,不得不说,想出这种方法进行自动微分的人简直是天才,我这种智障连看懂都需要很久的时间,别说自己提出了orz

首先需要构建一个类,这种数据类型就是计算图中的节点 Node:

import math

class Node():
    '''
    计算图的节点
    '''
    global_id = -1
    
    def __init__(self, operation, inputs):
        self.id = Node.global_id      # 该节点的id
        Node.global_id += 1
        
        self.inputs = inputs          # 该节点的输入
        self.operation = operation    # 该节点进行的运算操作
        self.grad = 0.0               # 该节点的初始化梯度
        self.calculate()              # 初始化时即计算出该节点的值
    
    def input2value(self):
        '''
        节点的输入是另一个节点,但是计算要取具体的数值
        '''
        value_input = []
        for i in self.inputs:
            if isinstance(i, Node):
                i = i.value
            value_input.append(i)
        return value_input
    
    def calculate(self):
        self.value = self.operation.compute(self.input2value())              # 该节点求输出值
        # 写到函数里是为了方便可以在构建 topo_list 的时候计算
        
    def  __str__(self):
        '''
        返回对象的描述信息
        '''
        return "Node %d : %s %s = %s, grad: %.3f" %(self.id, self.input2value(), self.operation.name, self.value, self.grad)

然后需要定义所有需要用到的原子操作类型,这里实现了:加法、减法、乘法、除法、平方、sigmoid、ln

class Operation():
    '''
    所有原子操作的基类,所有原子操作的输入都不多于两个
    '''
    def __init__(self):
        self.name = ''
        
    '''
    产生一个新的 Node,表示此次计算的结果。后面构建 topo_list 的时候会用到,每个操作必须要为输入产生孩子节点
    __call__ 表示该类可以像方法一样被调用
    '''
    def __call__(self):
        pass
    
    '''
    该操作的计算
    '''
    def compute(self, inputs):
        pass
    
    '''
    该操作的导数
    '''
    def gradient(self, output_grad):
        pass
class Identity(Operation):
    '''
    该操作将输入变量转化为接节点
    '''
    def __init__(self):
        self.name = "identity"
    
    def __call__(self, a):
        return Node(self, [a])
    
    def compute(self, inputs):
        return inputs[0]
    
    def gradient(self, inputs, output_grad):
        return [output_grad]
class Add(Operation):
    '''
    加法
    '''
    def __init__(self):
        self.name = "add"
    
    def __call__(self, a, b):
        return Node(self, [a,b])
    
    def compute(self, inputs):
        return inputs[0] + inputs[1]
    
    def gradient(self, inputs, output_grad):    # output_grad 似乎为外界对该操作整体的梯度
        return [output_grad, output_grad]

class Sub(Operation):
    '''
    减法
    '''
    def __init__(self):
        self.name = "sub"
        
    def __call__(self, a, b):
        return Node(self, [a,b])
    
    def compute(self, inputs):
        return inputs[0] - inputs[1]
    
    def gradient(self, inputs, output_grad):
        return [output_grad, -output_grad]

class Mul(Operation):
    '''
    乘法
    '''
    def __init__(self):
        self.name = "mul"
        
    def __call__(self, a, b):
        return Node(self, [a,b])
    
    def compute(self, inputs):
        return inputs[0] * inputs[1]
    
    def gradient(self, inputs, output_grad):
        return [inputs[1] * output_grad, inputs[0] * output_grad]

class Div(Operation):
    '''
    除法
    '''
    def __init__(self):
        self.name = "div"
        
    def __call__(self, a, b):
        return Node(self, [a,b])
    
    def compute(self, inputs):
        return inputs[0] / inputs[1]
    
    def gradient(self, inputs, output_grad):
        return [1.0 / inputs[1] * output_grad, -inputs[0] / inputs[1] ** 2 * output_grad]

class Square(Operation):
    '''
    平方
    '''
    def __init__(self):
        self.name = "square"
    
    def __call__(self, a):
        return Node(self, [a])
    
    def compute(self, inputs):
        return inputs[0] ** 2
    
    def gradient(self, inputs, output_grad):
        return [2 * inputs[0] * output_grad]

class Ln(Operation):
    '''
    取对数
    '''
    def __init__(self):
        self.name = "ln"
        
    def __call__(self, a):
        return Node(self, [a])
    
    def compute(self, inputs):
        return math.log(inputs[0])
    
    def gradient(self, inputs, output_grad):
        return [1.0 / inputs[0] * output_grad]

class Sigmoid(Operation):
    '''
    sigmoid 函数
    '''
    def __init__(self):
        self.name = "sigmoid"
        
    def __call__(self, a):
        return Node(self, [a])
    
    def compute(self, inputs):
        return 1.0 / (1 + math.exp(-inputs[0]))
    
    def gradient(self, inputs, output_grad):
        res = self.compute(inputs)
        return [res * (1 - res) * output_grad]

Executor 类是最重要最核心的一个类,它实现了对计算图进行前向计算、对所有节点进行拓扑排序、以及对输入变量进行求微分的操作。

class Executor():
    '''
    计算图的执行和自动微分
    '''
    def __init__(self, root):
        self.topo_list = self.topological_sort(root)    # 拓扑排序的顺序是孩子节点在前,父节点在后
        self.root = root
    
    '''
    给一个节点,将该节点下面的所有孩子节点从低到高加入到 topo_list 中
    '''
    def dfs(self, topo_list, node):
        if node is None or not isinstance(node, Node):
            return
        for term in node.inputs:
            self.dfs(topo_list, term)
        topo_list.append(node)    # 先添加孩子节点,再添加自己;某孩子节点可能会被添加多次
        # Q1:node 为空的判别条件?
        # Q2:能不能先添加自己,再添加孩子节点?
    
    '''
    给定根节点,将根节点下面的所有节点加入 topo_list
    '''
    def topological_sort(self, root):
        li = []
        self.dfs(li, root)
        return li
    
    '''
    计算每个节点的前向值。其实前面生成节点的时候已经计算过了
    '''
    def node_calculate(self):
        node_calcued = set()
        for term in self.topo_list:
            if term not in node_calcued:    # 保证取到的节点是还没有被计算过的
                term.calculate()
                node_calcued.add(term)
        return self.root.value              # 返回最终前向计算的结果
    
    '''
    反向对各输入变量求微分
    '''
    def derivation(self):
        reverse_topo_list = list(reversed(self.topo_list))
        reverse_topo_list[0].grad = 1.0
        
        for node in reverse_topo_list:        # Q:有重复节点怎么办?
            # 节点对其输入各变量的梯度
            print(node)
            
            grads = node.operation.gradient(node.input2value(), node.grad)
            for child, child_grad in zip(node.inputs, grads):
                if isinstance(child, Node):
                    child.grad += child_grad
        
        for node in reverse_topo_list:
            print(node)

举几个栗子试验一下:

# 原子操作实例化
identity = Identity()
add = Add()
sub = Sub()
mul = Mul()
div = Div()
square = Square()
ln = Ln()
sigmoid = Sigmoid()
# 第一个例子:y = ln(x1) + x1*x2 
x1, x2 = identity(2), identity(5)
y = add(ln(x1), mul(x1, x2))           # y 是根节点
executor = Executor(y)
print(executor.node_calculate())
executor.derivation()
print("gradient of x1 is %.3f" % x1.grad)
print("gradient of x2 is %.3f" % x2.grad)
# 第二个例子:y = (x1 + x2) * (x2 + 1)
x1, x2 = identity(2), identity(1)
y = mul(add(x1, x2), add(x2, 1))
executor = Executor(y)
print(executor.node_calculate())
executor.derivation()
print("gradient of x1 is %.3f" % x1.grad)
print("gradient of x2 is %.3f" % x2.grad)
# 第三个例子:用= (x1 + σ(x2)) / (σ(x1) + (x1 + x2)^2)
x1, x2 = identity(3), identity(-4)
numerator = add(x1, sigmoid(x2))
denominator = add(sigmoid(x1), square(add(x1, x2)))
y = div(numerator, denominator)
executor = Executor(y)
print(executor.node_calculate())
executor.derivation()
print("gradient of x1 is %.3f" % x1.grad)
print("gradient of x2 is %.3f" % x2.grad)

这里实现的是静态图,现在的深度学习框架里用的都是动态图了。

上一篇:线性回归从零实现


下一篇:RuntimeError element 0 of tensors does not require grad and does not have a grad_fn..