梯度下降
梯度下降法的原理
梯度下降法(gradient descent)是一种常用的一阶(first-order)优化方法，是求解无约束优化问题最简单、最经典的方法之一。
梯度下降最典型的例子就是从山上往下走，每次都寻找当前位置最陡峭的方向小碎步往下走，最终就会到达山下（暂不考虑有山谷的情况）。
首先来解释什么是梯度？这就要先讲微分。对于微分，相信大家都不陌生，看几个例子就更加熟悉了。
先来看单变量的微分：
再看多变量的微分：
补充:导数和微分的区别
导数是函数在某一点处的斜率，是Δy和Δx的比值；而微分是指函数在某一点处的切线在横坐标取得增量Δx以后，纵坐标取得的增量，一般表示为dy。
梯度就是由微分结果组成的向量，令
有
那么，函数f(x,y,z)在(1,2,3)处的微分为
因此，函数f(x,y,z)在(1,2,3)处的梯度为(6,11,6)。
梯度是一个向量，对于一元函数，梯度就是该点处的导数，表示切线的斜率。对于多元函数，梯度的方向就是函数在该点上升最快的方向。
梯度下降法就是每次都寻找梯度的反方向，这样就能到达局部的最低点。
那为什么按照梯度的反方向能到达局部的最低点呢？这个问题直观上很容易看出来，但严禁起见，我们还是给出数学证明。
对于连续可微函数f(x)，从某个随机点出发，想找到局部最低点，可以通过构造一个序列
那么我们就能够不断执行该过程即可收敛到局部极小点，可参考下图。
那么问题就是如何找到下一个点
其中
若想
步长
因此，有
每一步我们都按照
这里再对
在梯度前加负号就是朝梯度的反方向前进，因为梯度是上升最快的方向，所以方向就是下降最快的方向。
梯度下降的实例
一元函数的梯度下降
设一元函数为
函数的微分为
设起点为
经过4次迭代：
多元函数的梯度下降
设二元函数为
函数的梯度为
设起点为(2,3)，步长
loss function(损失函数)
损失函数也叫代价函数(cost function)，是用来衡量模型预测出来的值h(θ)与真实值y之间的差异的函数，如果有多个样本，则可以将所有代价函数的取值求均值，记做J(θ)。代价函数有下面几个性质：
- 对于每种算法来说，代价函数不是唯一的；
- 代价函数是参数θ的函数；
- 总的代价函数J(θ)可以用来评价模型的好坏，代价函数越小说明模型和参数越符合训练样本(x, y)；
- J(θ)是一个标量。
最常见的代价函数是均方误差函数，即
其中，
- m为训练样本的个数
表示估计值，表达式如下
- y是原训练样本中的值
我们需要做的就是找到θ的值，使得J(θ)最小。代价函数的图形跟我们上面画过的图很像，如下图所示。
看到这个图，相信大家也就知道了我们可以用梯度下降算法来求可以使代价函数最小的θ值。
先求代价函数的梯度
这里有两个变量
这么看公式可能很多同学会不太明白，我们把每个矩阵的具体内容表示出来，大家就很容易理解了。
矩阵
矩阵X为：
矩阵y为:
这样写出来后再去对应上面的公式，就很容易理解了。
下面我们来举一个用梯度下降算法来实现线性回归的例子。有一组数据如下图所示，我们尝试用求出这些点的线性回归模型。
首先产生矩阵X和矩阵y
# generate matrix X
按照上面的公式定义梯度函数
def gradient_function(theta, X, y):diff = np.dot(X, theta) - yreturn (1./m) * np.dot(np.transpose(X), diff)
接下来就是最重要的梯度下降算法，我们取
def gradient_descent(X, y, alpha):theta = np.array([1, 1]).reshape(2, 1)gradient = gradient_function(theta, X, y)while not np.all(np.absolute(gradient) <= 1e-5):theta = theta - alpha * gradientgradient = gradient_function(theta, X, y)return theta通过该过程，最终求出的
附录 source code
matlab一元函数的梯度下降程序
clc;
close all;
clear all;
%%
delta = 1/100000;
x = -1.1:delta:1.1;
y = x.^2;
dot = [1, 0.2, 0.04, 0.008];
figure;plot(x,y);
axis([-1.2, 1.2, -0.2, 1.3]);
grid on
hold on
plot(dot, dot.^2,'r');
for i=1:length(dot)text(dot(i),dot(i)^2,['theta_{',num2str(i),'}']);
end
title('一元函数的梯度下降过程');
python一元函数的梯度下降程序
import numpy as np
import matplotlib.pyplot as pltdelta = 1/100000
x = np.arange(-1.1, 1.1, delta)
y = x ** 2
dot = np.array([1, 0.2, 0.04, 0.008])
plt.figure(figsize=(7,5))
plt.plot(x,y)
plt.grid(True)
plt.xlim(-1.2, 1.2)
plt.ylim(-0.2, 1.3)
plt.plot(dot, dot**2, 'r')
for i in range(len(dot)):plt.text(dot[i],dot[i]**2,r'$theta_%d$' % i)
plt.title('一元函数的梯度下降过程')
plt.show()
julia一元函数的梯度下降程序
using PyPlot
delta = 1/100000
x = -1.1:delta:1.1
y = x.^2
dot = [1, 0.2, 0.04, 0.008]
plot(x, y)
grid(true)
axis("tight")
plot(dot, dot.^2, color="r")
for i=1:length(dot)text(dot[i], dot[i]^2, "$theta_$i$")
end
title("Single variable function gradient descent")
matlab二元函数的梯度下降程序
pecision = 1/100;
[x,y] = meshgrid(-3.1:pecision:3.1);
z = x.^2 + y.^2;
figure;
mesh(x,y,z);
dot = [[2,3];[1.6,2.4];[1.28,1.92];[5.09e-10, 7.64e-10]];
hold on
scatter3(dot(:,1),dot(:,2),dot(:,1).^2+dot(:,2).^2,'r*');
for i=1:4text(dot(i,1)+0.4,dot(i,2),dot(i,1).^2+0.2+dot(i,2).^2+0.2,['theta_{',num2str(i),'}']);
end
title('二元函数的梯度下降过程')
python二元函数的梯度下降程序
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
x = np.linspace(-3.1,3.1,300)
y = np.linspace(-3.1,3.1,300)
x,y = np.meshgrid(x, y)
z = x**2 + y**2
dot = np.array([[2,3],[1.6,2.4],[1.28,1.92],[5.09e-10, 7.64e-10]])
fig = plt.figure(figsize = (10,6))
ax = fig.gca(projection = '3d')
cm = plt.cm.get_cmap('YlGnBu')
surf = ax.plot_surface(x, y, z, cmap=cm)
fig.colorbar(surf,shrink=0.5, aspect=5)
ax.scatter3D(dot[:,0], dot[:,1], dot[:,0]**2 + dot[:,1]**2, marker='H',c='r')
for i in range(len(dot)-1):ax.text(dot[i,0]+0.4, dot[i,1], dot[i,0]**2 + dot[i,1]**2, r'$Theta_%d$' % i)
ax.text(dot[3,0]+0.4, dot[3,1]+0.4, dot[3,0]**2 + dot[3,1]**2-0.4, r'min')
plt.show()
julia二元函数的梯度下降程序
这个图的text死活标不上，希望知道的朋友可以告知一下。再多说一句，虽然我之前出了个Julia的教程，里面也包含4种绘图工具(Plots,GR,Gadfly & PyPlot)，但没有画过3维的图形，今天为了画这个图可真是费尽周折，Julia官网上的3D绘图的程序基本没有一个可以直接使用的，具体的绘图过程和调试中碰到的问题我还会整理篇文章到知乎和公众号，大家可以看一下。
using Plots
Plots.plotlyjs()
n = 50
x = range(-3, stop=3, length=n)
y= x
z = zeros(n,n)
for i in 1:n, k in 1:nz[i,k] = x[i]^2 + y[k]^2
endsurface(x, y, z)
dot = [[2 3]; [1.6 2.4]; [1.28 1.92]; [5.09e-10 7.64e-10]]
scatter!(dot[:,1], dot[:,2], dot[:,1].^2 .+ dot[:,2].^2)
matlab梯度下降的线性回归
m = 18;
X0 = ones(m,1);
X1 = (1:m)';
X = [X0, X1];
y = [2,3,3,5,8,10,10,13,15,15,16,19,19,20,22,22,25,28]';
alpha = 0.01;
theta = gradient_descent(X, y, alpha, m);function [grad_res] = gradient_function(theta, X, y, m)diff = X * theta - y;grad_res = X' * diff / m;
endfunction [theta_res] = gradient_descent(X, y, alpha, m)theta = [1;1];gradient = gradient_function(theta, X, y, m);while sum(abs(gradient)>1e-5)>=1theta = theta - alpha * gradient;gradient = gradient_function(theta, X, y, m);endtheta_res = theta;
end
python梯度下降的线性回归
import numpy as np
import matplotlib.pyplot as plt # y = np.array([2,3,3,5,8,10,10,13,15,15,16,19,19,20,22,22,25,28])
# x = np.arange(1,len(y)+1)
# plt.figure()
# plt.scatter(x,y)
# plt.grid(True)
# plt.show()# sample length
m = 18# generate matrix X
X0 = np.ones((m, 1))
X1 = np.arange(1, m+1).reshape(m, 1)
X = np.hstack((X0, X1))# matrix y
y = np.array([2,3,3,5,8,10,10,13,15,15,16,19,19,20,22,22,25,28]).reshape(m,1)# alpha
alpha = 0.01def cost_function(theta, X, y):diff = np.dot(X, theta) - yreturn (1./2*m) * np.dot(np.transpose(diff), diff)def gradient_function(theta, X, y):diff = np.dot(X, theta) - yreturn (1./m) * np.dot(np.transpose(X), diff)def gradient_descent(X, y, alpha):theta = np.array([1, 1]).reshape(2, 1)gradient = gradient_function(theta, X, y)while not np.all(np.absolute(gradient) <= 1e-5):theta = theta - alpha * gradientgradient = gradient_function(theta, X, y)return theta[theta0, theta1] = gradient_descent(X, y, alpha)
plt.figure()
plt.scatter(X1,y)
plt.plot(X1, theta0 + theta1*X1, color='r')
plt.title('基于梯度下降算法的线性回归拟合')
plt.grid(True)
plt.show()
julia梯度下降的线性回归
m = 18
X0 = ones(m,1)
X1 = Array(1:m)
X = [X0 X1];
y = [2,3,3,5,8,10,10,13,15,15,16,19,19,20,22,22,25,28];
alpha = 0.01;
theta = gradient_descent(X, y, alpha, m)function gradient_function(theta, X, y, m)diff = X * theta .- y;grad_res = X' * diff / m;
endfunction gradient_descent(X, y, alpha, m)theta = [1,1]gradient = gradient_function(theta, X, y, m)while all(abs.(gradient) .>1e-5)==truetheta = theta - alpha * gradientgradient = gradient_function(theta, X, y, m)endtheta_res = theta
end最后附上我出过的免费的Julia教学视频链接：
Julia 教程 从入门到进阶 - 网易云课堂study.163.com微信公众号：Quant_Times
