自动求导程序的设计与实现（Python)

最新推荐文章于 2023-01-03 15:44:22 发布

TaiJi1985

最新推荐文章于 2023-01-03 15:44:22 发布

阅读量1.1w

点赞数 14

分类专栏：机器学习 python 文章标签： python 求导自动求导表达式表达式树

本文链接：https://blog.csdn.net/taiji1985/article/details/72942947

版权

机器学习同时被 2 个专栏收录

21 篇文章 3 订阅

订阅专栏

python

10 篇文章 1 订阅

订阅专栏

动机

作者 Yangtf

最近一直在求各种导数，于是就想写一个自动求导的算法。其实python中的theano就有这个功能，但想了想，思路不难，于是就动手实现了一个。

本来想用c++实现了，但发现c++写各种问题，内存管理、操作符重载都不尽人意。花费了不少时间后，决定换语言。 Java是第一熟练语言，但不支持操作符重载，奈何？于是转战python。

源代码路径

最新的源代码在这里。
http://git.oschina.net/yangtf/python_exp

思路

函数的表示

将函数表达式表示为一个表达式树。

这里写图片描述

那个这个表达式树如何构建呢？　要自己写语法分析么？　太麻烦，有种比较简单的办法，就是使用操作符重载来实现。

定义一个类Ｅ，重载它的 + - * / **（乘方）操作，在重载中，进行二叉树的构建。

节点类型

在这个表达式树中，主要应有三种节点类型。
其一，常数节点。如 2，3
其二，变量节点，如 a,b,x,y之类。
其三，操作节点。如 + , - ,* , / ,乘方等。

求导方法

有了表达式构成的二叉树，下面就是求导了。

对常数节点求导，结果为0 。
对变量节点求导，有两种情况。如

f (a, b) = a 2 + 3 b

$f(a,b) = a^2 + 3b$

这个函数对 $a$ 求偏导，那么就将b节点看成是一个常数,求导结果为0。
对于保存了a的节点，求导结果为1。

求导的方法就是那些求导公式，举例：

(x + y)' = x' + y'

$(x+y)' = x' + y'$

求导看这篇文章 http://blog.csdn.net/taiji1985/article/details/72857554

上面的公式，对于一个根为‘+’的二叉树，分别对其左子树和右子树进行求导，然后将求导得到的和相加。

那么如何求导左子树呢？，递归的调用这个求导方法就可以了。

对乘方节点的处理时比较难的。
这里写图片描述

先对左子树f求导，对右子树g求导。
如果f求导为0，说明是指数函数，如果g求导为0，说明是幂函数，分别套用公式。
至于 $f(x)^g(x)$ 这种形式，求导公式有点复杂，还要去请教一些数学方面的高手。还没有做。

化简

求导不是最难的，最难的是化简。比如对 1 / ( 1 + e ^ ( - ( w * x + b ) ) ) 按照上述算法求导，得到的结果是：

( 0 * ( 1 + e ^ ( - ( w * x + b ) ) ) - 1 * ( 0 + e ^ ( - ( w * x + b ) ) * 1 * ( 0 * ( w * x + b ) + - ( 1 * x + w * 0 + 0 ) ) ) ) / ( 1 + e ^ ( - ( w * x + b ) ) ) * ( 1 + e ^ ( - ( w * x + b ) ) )

这就需要化简。我实现了化简的几个思路：

（1） 0+x，x+0 x-0 这种化简为 x 。0*x x*0 0/x 化简为 0
这里写图片描述

在上图中，左图c节点为0，则应让a直接指向d。删除c和b节点。右图为1*x的图，应让a直接指向d。
（2）x*1 1*x x/1 这种直接简化为x
（3）两个常量进行运算，F+F， F-F， F*F， F/F 都简化为单一节点。
（4）较为复杂的节点合并。
这里写图片描述

在上图中，右子树有个3，左子树有一个4，算法

如果右子树是一个常量节点，则在左子树中查找与p指向节点符号相同的节点。经过三个星号，找到了4，然后3*4 ->12 ,随后删除原本p指向的节点，让p直接指向原本的左子树。

(5) $x*x => x^2$

(6) $0-x => -1*x$

(7) x^1 => x

(8) log e - > 1

代码实现


# -*- coding: UTF-8 -*-

'''
Created on 2017-6-8

@author: Administrator

二元运算符    特殊方法
+    __add__,__radd__
-    __sub__,__rsub__
*    __mul__,__rmul__
/    __div__,__rdiv__,__truediv__,__rtruediv__
//    __floordiv__,__rfloordiv__
%    __mod__,__rmod__
**    __pow__,__rpow__
<<    __lshift__,__rlshift__
>>    __rshift__,__rrshift__
&    __and__,__rand__
^    __xor__,__rxor__
|    __or__,__ror__
+=    __iaddr__
-=    __isub__
*=    __imul__
/=    __idiv__,__itruediv__
//=    __ifloordiv__
%=    __imod__
**=    __ipow__
<<=    __ilshift__
>>=    __irshift__
&=    __iand__
^=    __ixor__
|=    __ior__
==    __eq__
!=,<>    __ne__
>    __get__
<    __lt__
>=    __ge__
<=    __le__

'''

class E:
    def __init__(self):
        self.left=None;
        self.right=None;
        self.parent = None;
        self.type = 'n';
        self.f = 0;
        pass
    def isOp(self,op):
        return self.type == 'op' and self.f == op;

    def isZero(self):
        return self.type == 'float' and abs(self.f) < 1e-5;
    def isOne(self):
        return self.type == 'float' and abs(self.f -1 ) < 1e-5;
    def isNum(self):
        return self.type == 'float';
    def float(self,a): #
        self.f = a;
        self.left = self.right = None;
        self.type = 'float';
        return self;
    def sym(self,name): 
        self.type = 'sym';
        self.f = name;
        return self;
    def withOp(self,op,left,right):
        self.f = op;
        self.type = 'op';

        if type(left) == int or type(left) == float:
            left = E().float(left);

        if type(right) == int or type(right) == float:
            right = E().float(right);

        if left != None:
            self.left = left.clone();
            self.left.parent = self;
        else:
            self.left =None;
        if right != None:
            self.right = right.clone();
            self.right.parent = self;
        else:
            self.right = None;
        return self;
    def clone(self): #深度复制
        x = E();
        x.type = self.type;
        x.f = self.f;
        if self.left == None:
            x.left = None;
        else:
            x.left = self.left.clone();

        if self.right == None:
            x.right = None;
        else:
            x.right = self.right.clone();
        return x;

    def __radd__(self,x):
        #print '__radd__ ',x
        r = E().withOp('+', x,self);
        return r;
    def __rsub__(self,x):
        #print '__rsub__ ',x
        r = E().withOp('-', x,self);
        return r;
    def __rmul__(self,x):
        r = E().withOp('*', x,self);
        return r;
    def __rdiv__(self,x):
        r = E().withOp('/', x,self);
        return r;
    def __neg__(self):          
        r = E().withOp('*', E().float(-1),self);
        return r;
    def __add__(self,x):
        #print 'add ',x
        r = E().withOp('+', self, x);
        return r;
    def __sub__(self,x):
        r = E().withOp('-', self, x);
        return r;
    def __mul__(self,x):
        r = E().withOp('*', self, x);
        return r;
    def __div__(self,x):
        r = E().withOp('/', self, x);
        return r;
    def __pow__(self,x):
        r = E().withOp('^', self, x);
        return r;

    def isConstOf(self,x): # 求导时，对于x是否是一个常数
        if self.type == 'float':
            return True;
        if self.type == 'sym' :
            return self.f == x.f;

        return (self.left == None or self.left.isConstOf(x)) and (self.right == None or self.right.isConstOf(x));

    def op_diff(self,x):
        # do something with None left or right
        if self.left == None:
            d_left =None;
        else:
            d_left = self.left.diff(x);
        if self.right == None:
            d_right = None;
        else:
            d_right = self.right.diff(x);

        if self.f == '+':
            return d_left+d_right;
        if self.f == '-':
            return d_left-d_right;
        if self.f == '*':
            return d_left*self.right+self.left*d_right;
        if self.f == '/':
            return (d_left*self.right-self.left*d_right)/(self.right*self.right);
        if self.f == '^':
            left_c = d_left == E().float(0);
            right_c = d_right == E().float(0);

            if left_c and right_c :
                return E().float(0);
            elif right_c: # f(x)^a  ()' = a*f(x)^(a-1)*f'(x);
                return self.right*self.left**(self.right-1)*d_left;
            elif left_c: #指数 a^g(x)  ()' = a^g(x)*loga*g'(x)
                return self.left**self.right * self.left.log() * d_right;
            else:
                print 'unsupport f(x)^g(x) style!! now ' 
                exit(1);
        pass
    def diff(self,x): # 对x求偏导数
        if self.type == 'float':
            return E().float(0);
        elif self.type == 'sym':
            if x.f == self.f: # 是同一个变量
                return E().float(1);
            else:
                return E().float(0); #不是同一个变量。
        elif self.type == 'op':
            return self.op_diff(x);


        pass
    def eq(self,x,y):
        if x == None :
            return y == None;
        else :
            return x == y;
    def __eq__(self,x):
        if x == None:
            return False;
        if x.type != self.type:
            return False;
        if x.type == 'float':
            return abs(x.f - self.f)<1e-5;
        if x.type == 'sym':
            return x.f == self.f;
        if x.type == 'op':
            if x.f != self.f :
                return False;
            return self.eq(self.left,x.left) and self.eq(self.right,x.right);
    def printme(self):
        self.setParent();
        self._printme();
        print '';
    def _op_toi(self,op):
        if op == '+' or op == '-':
            return 10;
        if op == '*' or op == '/':
            return 20;
        if op == '^':
            return 30;
        return 40;
    def _compare_op(self,a,b): #比较两个符号，谁的优先级高
        #print 'compare ',a,b,self._op_toi(a) - self._op_toi(b);
        return self._op_toi(a) - self._op_toi(b);
    def _printme(self):
        if self.type == 'float':
            print self.f ,;
        elif self.type == 'op':
            useBrack = True;
            if self.parent == None:
                useBrack = False;
            elif self._compare_op(self.f, self.parent.f)>= 0:
                useBrack = False;


            if useBrack:
                print '(',;
            #如果是 -1*x ，直接输出 -x;
            if self.left !=None and self.left == E().float(-1) and self.isOp('*'):
                print '-',;
            else:
                if self.left !=None:
                    self.left._printme();
                print self.f ,;
            if self.right != None:
                self.right._printme();
            if useBrack:
                print ')',;
        elif self.type == 'sym':
            print self.f ,;
        pass


    def child_pattern(self,x):
        if x == None:
            return 'none';
        if x.left == None:
            lc= "N";
        elif x.left.isOne():
            lc = '1';
        elif x.left.isZero():
            lc = '0';
        elif x.left.type == 'float':
            lc = 'F';
        else:
            lc ='A';

        if x.right == None:
            rc= "N";
        elif x.right.isOne():
            rc = '1';
        elif x.right.isZero():
            rc = '0';
        elif x.right.type == 'float':
            rc = 'F';
        else :
            rc ='A';



        pt= str(lc)+str(x.f) + str(rc);
        #print "PT=",pt," -------------";
        #x.printme();

        return pt;


    def evalue(self,op,a,b):
        if op == '+':
            r= a.f+b.f;
        if op == '-':
            r= a.f-b.f;
        if op == '*':
            r= a.f*b.f;
        if op == '/':
            r= a.f/b.f;
        return r;

    def _node_op(self,r,op,v):
        # 在以r为根的树中，查找一个满足从根r到该节点整条路径上节点都与op相同的float节点，并将v中的数据应用op进去。
        if r == None :
            return False;
        if r.type == 'float' : # 如果当前节点就是一个float节点，把v的值乘在这里。
            r.f = r.evalue(op,r,v);
            return True;

        if r.type != 'op' or r.f != op: #当前节点不满足op相等条件
            return False;

        if self._node_op(r.left, op, v):
            return True;

        if self._node_op(r.right, op, v):
            return True;

        return False;

        pass
    def _node_join(self,r,x,y):
        #合并两个节点 2+(2+x) => 4+x;
        #r 如果不能合并应返回的值
        #x 判断x是否是一个数字，如果是，则看能否和y中节点合并
        if x==None or y == None or  x.type != 'float' :
            return r;

        succ = self._node_op(y, r.f, x); #如果成功将x乘进了y，则删除x，把y作为父。
        if succ:
            return y;
        return r;
        #在y中查找
#         if y.type == 'op' and y.type == r.type  and y.f == r.f:
#             if y.left != None and y.left.type=='float':
#                 y.left.f = self.evalue(y.f, x, y.left);
#                 
#                 return y;
#             if y.right != None and y.right.type=='float':
#                 y.right.f = self.evalue(y.f, x, y.right);
#                 return y;
#             
#        return r;

    def _opt_node(self,x):
        #左子树 0,1检测
        r = x;
        if x == None :
            return x;

        pt = self.child_pattern(x);
        if pt == 'F-1':
            pt = pt; # for debug

        if pt == '0*A' or pt == '0/A' or pt== 'A*0':
            r = E().float(0);
        if pt == '0+A' or pt == '0+1':
            r = x.right;
        if pt == 'A+0' or pt == '1+0':
            r = x.left;
        if pt == 'A*1':
            r = x.left;

        #左子树常数化简
        pt = self.child_pattern(x);
        pt = pt.replace('0', 'F').replace('1','F');
        #print '#####', pt;

        if pt == 'F+F':
            r = E().float(x.left.f+x.right.f);
        if pt == 'F-F':
            r = E().float(x.left.f-x.right.f);
        if pt == 'F*F':
            r = E().float(x.left.f*x.right.f);
        if pt == 'F/F':
            r = E().float(x.left.f/x.right.f);
        return r;
    def optm(self): # 优化式子
        # 后续遍历，从下网上优化
        if self.left!= None:
            self.left = self.left.optm();
        if self.right!=None:
            self.right = self.right.optm();

        self.left = self._opt_node(self.left);
        self.right = self._opt_node(self.right);

        r = self._opt_node(self);

        # 0-x -> -1*x
        if self.isOp('-'):
            if self.left!=None and self.left == E().float(0):
                self.f = '*';
                self.left = E().float(-1);


        #优化常数项（多个常数项相乘，如2*3*x ->6*x)
        r = self._node_join(r,r.left,r.right);
        r = self._node_join(r,r.right,r.left);

        if r.left != None and r.left == r.right:
            if r.isOp('*'):
                r.f = '^';
                r.right = E().float(2);


        #优化乘方
        if r.isOp('^') and r.right != None and r.right.isOne():
            return r.left;


        return r;
        pass

    #求以e为底的对数
    def log(self):
        if self.type == 'sym' and self.f == 'e':
            return E().float(1);
        r = E().withOp('log', None, self);
        return r;

    #设置所有parent指针

    def setParent(self):
        if self.left !=None :
            self.left.parent = self;
            self.left.setParent();
        if self.right != None:
            self.right.parent = self;
            self.right.setParent();
pass


# class Optmer:
#     def __init__(self):
#         pass
#     def addParentPointer(self,tree):
#         if tree.left != None:
#             tree.left.parent = tree;
#             self.addParentPointer(tree.left);
#         if tree.right != None:
#             tree.right.parent = tree;
#             self.addParentPointer(tree.right);
#         
#     def optNode(self,node):
#         self.addParentPointer(node);
# 
#     def _zeroOptNode(self,node):
#         if node == None:
#             return;
#         if node.isZero():
#             node.parent.
#         pass

x = E().sym('x');
#c = 2*x**2+3*x**4+E().float(4)**x;
e = E().sym('e');
w = E().sym('w');
b = E().sym('b');
c = 1/(1+e**(-(w*x+b)));
c.printme();

d = c.diff(w);
d.printme();
d.optm().optm().printme();

运行测试

以 sigmoid函数为例，进行求导。

待求导的函数
1 / ( 1 + e ^ ( - ( w * x + b ) ) )

求导后，化简前
( 0 * ( 1 + e ^ ( - ( w * x + b ) ) ) - 1 * ( 0 + e ^ ( - ( w * x + b ) ) * 1 * ( 0 * ( w * x + b ) + - ( 1 * x + w * 0 + 0 ) ) ) ) / ( 1 + e ^ ( - ( w * x + b ) ) ) * ( 1 + e ^ ( - ( w * x + b ) ) )
化简后，中间还是有一个1在哪里，问题在哪里太晚了，不查了。结果是对的。
e ^ ( - ( w * x + b ) ) * 1 * x / ( 1 + e ^ ( - ( w * x + b ) ) ) ^ 2