# From calculation, it is expected that the local minimum occurs at x=9/4

cur_x = 6 # The algorithm starts at x=6
gamma = 0.01 # step size multiplier
precision = 0.00001
previous_step_size = cur_x

def df(x):
    return 4 * x**3 - 9 * x**2

while previous_step_size > precision:
    prev_x = cur_x
    cur_x += -gamma * df(prev_x)
    previous_step_size = abs(cur_x - prev_x)

print("The local minimum occurs at %f" % cur_x)

The local minimum occurs at 2.249965

注：关于“gradient descent”还可以参考以下资料：
https://www.analyticsvidhya.com/blog/2017/03/introduction-to-gradient-descent-algorithm-along-its-variants/
http://ruder.io/optimizing-gradient-descent/


3，梯度下降法与最小二乘法
  “机器学习”中有六个经典算法，其中就包括“最小二乘法”和“梯度下降法”，前者用于“搜索最小误差”，后者用于“用最快的速度搜索”，二者常常配合使用。代码演示如下：

# y = mx + b
# m is slope, b is y-intercept
def compute_error_for_line_given_points(b, m, coordinates):
    totalerror = 0
    for i in range(0, len(coordinates)):
        x = coordinates[i][0]
        y = coordinates[i][1]
        totalerror += (y - (m * x + b)) ** 2
    return totalerror / float(len(coordinates))

# example
compute_error_for_line_given_points(1, 2, [[3, 6], [6, 9], [12, 18]])

22.0

  以上就是用“最小二乘法”来计算误差，当输入为 $(1,2)$ 时，输出为 $22.0$

  很显然，最小二乘法需要不停地调整（试验）输入来找到一个最小误差。而应用“梯度下降法”，可以加快这个“试验”的过程。
  以上面这段程序为例，误差是斜率 m 和常数 b 的二元函数，可以表示为

$$e = g(m,b)$$
  那么，对最小二乘法的参数调优就转变为了求这个二元函数的极值问题，也就是说可以应用“梯度下降法”了。
  “梯度下降法”可以用于搜索函数的局部极值，如下，求下列函数的局部极小值
$$f(x) = x^5 - 2x^3 -2$$
  分析：这是一个一元连续函数，且可导，其导函数是：
$$f'(x) = 5x^4 - 6x^2$$
  根据“一阶导数极值判别法”：若函数f(x)可导，且f’(x)在 $x_0$ 的两侧异号，则 $x_0$ 是f(x)的极值点。那么，怎么找到这个 $x_0$ 呢？
  很简单，只需要沿斜率（导数值）的反方向逐步移动即可，如下图：导数为负时，沿x轴正向移动；导数为正时，沿x轴负方向移动。

current_x = 0.5 # the algorithm starts at x=0.5
learning_rate = 0.01 # step size multiplier
num_iterations = 60 # the number of times to train the function

# the derivative of the error function (x ** 4 = the power of 4 or x^4)
def slope_at_given_x_value(x):
    return 5 * x ** 4 - 6 * x ** 2

# Move X to the right or left depending on the slope of the error function
x = [current_x]
for i in range(num_iterations):
    previous_x = current_x
    current_x += -learning_rate * slope_at_given_x_value(previous_x)
    x.append(current_x)   #print(previous_x)

print("The local minimum occurs at %f, it is %f" % (current_x, current_x ** 5 - 2 * current_x ** 3 - 2))

The local minimum occurs at 1.092837, it is -3.051583

import numpy as np
import matplotlib.pyplot as plt
plt.plot(x, marker='*')
plt.show()

  沿梯度（斜率）的反方向移动，这就是“梯度下降法”。如上图所示，不管初始化值设为什么，在迭代过程只会越来越接近目标值，而不会偏离目标值，这就是梯度下降法的魅力。
  上面这张图是表示的是一个一元函数搜索极值的问题，未必能很好展示梯度下降法的魅力，你再返回去看上面那张“势能梯度图”，那是一个二元函数搜索极值的过程。左边的搜索路径很简洁，而右边的搜索路径，尽管因为初始值的设定，导致它的路径很曲折，但是，你有没有发现，它的每一次迭代事实上离目标都更近一步。我想，这就是梯度下降法的优点吧！
注：这段代码是一元函数求极值，如果是二元函数，则需要同时满足两个分量的偏导数的值为零，下面的线性回归程序算的就是二元偏导数。
  通过组合最小二乘法和梯度下降法，你可以得到线性回归，如下：

# Price of wheat/kg and the average price of bread
wheat_and_bread = [[0.5,5],[0.6,5.5],[0.8,6],[1.1,6.8],[1.4,7]]

def step_gradient(b_current, m_current, points, learningRate):
    b_gradient = 0
    m_gradient = 0
    N = float(len(points))
    for i in range(0, len(points)):
        x = points[i][0]
        y = points[i][1]
        b_gradient += -(2/N) * (y -((m_current * x) + b_current))
        m_gradient += -(2/N) * x * (y -((m_current * x) + b_current))
    new_b = b_current -(learningRate * b_gradient)
    new_m = m_current -(learningRate * m_gradient)
    return [new_b, new_m]

def gradient_descent_runner(points, starting_b, starting_m, learning_rate, num_iterations):
    b = starting_b
    m = starting_m
    for i in range(num_iterations):
        b, m = step_gradient(b, m, points, learning_rate)
    return [b, m]

gradient_descent_runner(wheat_and_bread, 1, 1, 0.01, 1000)

[3.853945094921183, 2.4895803107016445]

  上面这个程序的核心思想就是：在内层迭代的过程中，算出每一步误差函数相当于 m 和 b 的偏导数（梯度），然后沿梯度的反方向调整 m 和 b ；外层迭代执行梯度下降法，逐步逼近偏导数等于0的点。
  其中需要注意偏导数的近似计算公式，已知误差函数

import numpy as np
import matplotlib.pyplot as plt
a = np.array(wheat_and_bread)
plt.plot(a[:,0], a[:,1], 'ro')
b,m = gradient_descent_runner(wheat_and_bread, 1, 1, 0.01, 1000)
x = np.linspace(a[0,0], a[-1,0])
y = m * x + b
plt.plot(x, y)
plt.grid()
plt.show()

对比Numpy

import numpy as np
import matplotlib.pyplot as plt
a = np.array(wheat_and_bread)
plt.plot(a[:,0], a[:,1], 'ro')
m, b = np.polyfit(a[:,0], a[:,1], 1)
print([b,m])
x = np.linspace(a[0,0], a[-1,0])
y = m * x + b
plt.plot(x, y)
plt.grid()
plt.show()

[4.1072992700729891, 2.2189781021897814]