python 求多元函数的梯度函数 多元函数梯度的计算

1.导数的实现

python中有两种常见求导的方法，一种是使用Scipy库中的derivative方法，另一种就Sympy库中的diff方法。

1.1 Scipy

scipy.misc.derivative(func, x0, dx=1.0, n=1, args=(), order=3)[source]

在一个点上找到函数的第n个导数。即给定一个函数，请使用间距为dx的中心差分公式来计算x0处的第n个导数。

参数：

• func：需要求导的函数，只写参数名即可，不要写括号，否则会报错
• x0：要求导的那个点，float类型
• dx（可选）：间距，应该是一个很小的数，float类型
• n（可选）：n阶导数。默认值为1，int类型
• args（可选）：参数元组
• order（可选）：使用的点数必须是奇数，int类型
  求一阶导数的例子：

from scipy.misc import derivative
def f(x):
    return x**3 + x**2
derivative(f, 1.0, dx=1e-6)

4.999999999921734

1.2 Sympy表达式求导

sympy是符号化运算库，能够实现表达式的求导。所谓符号化，是将数学公式以直观符号的形式输出。下面看几个例子就明白了。

from sympy import Symbol

# 符号化变量
x = Symbol('x')

func = 1/(1+x**2)

print("x:", type(x))
print(func)
print(diff(func, x))
print(diff(func, x).subs(x, 3))
print(diff(func, x).subs(x, 3).evalf())

x: <class 'sympy.core.symbol.Symbol'>
1/(x**2 + 1)
-2*x/(x**2 + 1)**2
-3/50
-0.0600000000000000

2.模拟实现梯度下降

2.1 封装函数

首先构造一个损失函数 $loss(x)=(x-2.5)^2-1$，然后创建在-1到6的范围内构建140个点，并且求出对应的损失函数值，这样就可以画出损失函数的图形。

import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import derivative

def lossFunction(x):
    return (x-2.5)**2-1

# 在-1到6的范围内构建140个点
plot_x = np.linspace(-1,6,141)
# plot_y 是对应的损失函数值
plot_y = lossFunction(plot_x)

plt.plot(plot_x,plot_y)

[<matplotlib.lines.Line2D at 0x15d3ebd4240>]

已知梯度下降的本质是多元函数的导数，这了定义一个求导的方法，使用的是scipy库中的derivative方法。

"""
算法：计算损失函数J在当前点的对应导数
输入：当前数据点theta
输出：点在损失函数上的导数
"""
def dLF(theta):
    return derivative(lossFunction, theta, dx=1e-6)

接下来我们就可以进行梯度下降的操作了。首先我们需要定义一个点 $\theta$ 作为初始值，正常应该是随机的点，但是这里先直接定为0。然后需要定义学习率 $\eta$，也就是每次下降的步长。这样的话，点 $\theta$ 每次沿着梯度的反方向移动$\eta$ 距离，即 $\theta=\theta-\eta*\nabla f$

那么还有一个问题：如何结束循环呢？梯度下降的目的是找到一个点$\theta$ ，使得损失函数值最小，因为梯度是不断下降的，所以新的点$\theta$

theta = 0.0
eta = 0.1
epsilon = 1e-6
while True:
    # 每一轮循环后，要求当前这个点的梯度是多少
    gradient = dLF(theta)
    last_theta = theta
    # 移动点，沿梯度的反方向移动步长eta
    theta = theta - eta * gradient
    # 判断theta是否达到最小值
    # 因为梯度在不断下降，因此新theta的损失函数在不断减小
    # 看差值是否达到了要求
    if(abs(lossFunction(theta) - lossFunction(last_theta)) < epsilon):
        break
print(theta)
print(lossFunction(theta))

2.498732349398569
-0.9999983930619527

下面可以创建一个用于存放所有$\theta$点位置的列表，然后将其在图上绘制出来。

为了方便测试，可以将其封装成函数进行调用。

def gradient_descent(initial_theta, eta, epsilon=1e-6):
    theta = initial_theta
    theta_history.append(theta)
    while True:
        # 每一轮循环后，要求当前这个点的梯度是多少
        gradient = dLF(theta)
        last_theta = theta
        # 移动点，沿梯度的反方向移动步长eta
        theta = theta - eta * gradient
        theta_history.append(theta)
        # 判断theta是否达到损失函数最小值的位置
        if(abs(lossFunction(theta) - lossFunction(last_theta)) < epsilon):
            break

def plot_theta_history():
    plt.plot(plot_x,plot_y)
    plt.plot(np.array(theta_history), lossFunction(np.array(theta_history)), color='red', marker='o')
    plt.show()

2.2 调整学习率

首先使用学习率 $\eta = 0.1$

eta=0.1
theta_history = []
gradient_descent(0., eta)
plot_theta_history()
print("梯度下降查找次数：",len(theta_history))

梯度下降查找次数： 35

我们发现，在刚开始时移动比较大，这是因为学习率是一定的，再乘上梯度本身数值大（比较陡），后来梯度数值小（平缓）所以移动的比较小。且经历了34次查找。

使用使用学习率 $\eta = 0.01$

eta=0.01
theta_history = []
gradient_descent(0., eta)
plot_theta_history()
print("梯度下降查找次数：",len(theta_history))

梯度下降查找次数： 310

可见学习率变低了，每一步都很小，因此需要花费更多的步数。如果我们将学习率调大，会发生什么？

eta=0.9
theta_history = []
gradient_descent(0., eta)
plot_theta_history()
print("梯度下降查找次数：",len(theta_history))

梯度下降查找次数： 35

可见在一定范围内将学习率调大，还是会逐渐收敛的。但是我们要注意，如果学习率调的过大， 一步迈到“损失函数值增加”的点上去了，在错误的道路上越走越远（如下图所示），就会导致不收敛，会报OverflowError的异常。

为了避免报错，可以对原代码进行改进：

• 在计算损失函数值时捕获一场

def lossFunction(x):
    try:
        return (x-2.5)**2-1
    except:
        return float('inf')

• 设定条件，结束死循环

def gradient_descent(initial_theta, eta, n_iters, epsilon=1e-6):
    theta = initial_theta
    theta_history.append(theta)
    i_iters = 0
    while i_iters < n_iters:
        gradient = dLF(theta)
        last_theta = theta
        theta = theta - eta * gradient
        theta_history.append(theta)
        if(abs(lossFunction(theta) - lossFunction(last_theta)) < epsilon):
            break
        i_iters += 1

总结

梯度是向量，求梯度就要求导数。在python中，除了自己手动计算以外，还有两个常用的求导方法：Scipy & Sympy。

在求出导数之后就可以模拟梯度下降的过程编写代码了，这里面还要注意退出循环的条件。当学习率过小，收敛学习速度变慢，使得算法的效率降低；学习率过大又会导致不收敛，在“错误的道路上”越走越远。我们要对异常进行进行处理。

这样，我们就手动实现了梯度下降算法。