淮北做网站,黄骅贴吧超市转租信息,企业qq注册申请,网站换源码如何保留以前的文章线性回归是机器学习里面最简单也是最常用的算法#xff0c;理解了线性回归的推导之后对于后续的学习有很大帮助#xff0c;所以我决定从这里开始深入学习相关的机器学习模型。 本篇首先从矩阵求导开始切入#xff0c;然后介绍一次线性回归的推导#xff0c;再到代码实现。本…线性回归是机器学习里面最简单也是最常用的算法理解了线性回归的推导之后对于后续的学习有很大帮助所以我决定从这里开始深入学习相关的机器学习模型。 本篇首先从矩阵求导开始切入然后介绍一次线性回归的推导再到代码实现。本人也是入门哈有什么错误望轻喷。
求导
简单求导规则
首先先不管模型本身对于机器学习来说各种矩阵的运算非常重要而求导又是最基础的所以就放在最前面。主要参考
矩阵和标量之间的求导形状取决于矩阵矩阵包括向量无论是标量对矩阵求导矩阵对标量求导形状都和矩阵一致每个元素都是标量求导。
矩阵对矩阵求导矩阵求导本质上是把向量化为标量的过程假设存在变量 y → m × 1 \overrightarrow{y}_{m×1} y m×1和对应的函数 f ( y → ) n × 1 f(\overrightarrow{y})_{n×1} f(y )n×1其中每一行为 f i ( y → ) f_{i}(\overrightarrow{y}) fi(y )此时存在求导法则 ∂ f ( y → ) ∂ y → ( ∂ f 1 ( y → ) ∂ y 1 ⋯ ∂ f 1 ( y → ) ∂ y m ⋮ ⋱ ⋮ ∂ f n ( y → ) ∂ y 1 ⋯ ∂ f n ( y → ) ∂ y m ) n × m \frac{ \partial f(\overrightarrow{y}) }{ \partial \overrightarrow{y} }\begin{pmatrix} \frac{ \partial f_{1}(\overrightarrow{y}) }{ \partial {y_1} } \cdots \frac{ \partial f_{1}(\overrightarrow{y}) }{ \partial {y_m} }\\ \vdots \ddots \vdots\\ \frac{ \partial f_{n}(\overrightarrow{y}) }{ \partial {y_1} } \cdots \frac{ \partial f_{n}(\overrightarrow{y}) }{ \partial {y_m} } \end{pmatrix}_{n×m} ∂y ∂f(y ) ∂y1∂f1(y )⋮∂y1∂fn(y )⋯⋱⋯∂ym∂f1(y )⋮∂ym∂fn(y ) n×m 在矩阵运算中由于存在行列的概念因此需要区别两种布局方式分子布局和分母布局 分子布局分子布局指的是求导后的矩阵行数和分子对应行数相等。 分母布局分母布局指的是求导后的矩阵行数和分母对应行数相等。 显然上面的公式为分子布局因为分子 n n n行 ( f n × 1 ) (f_{n×1}) (fn×1)最终结果也是 n n n行下面都采用分子布局计算。
下面有两个常用的偏导式可以用于加速导数计算 假设存在方阵 A A A和向量 y → \overrightarrow{y} y ∂ A y → ∂ y → A T \frac{ \partial A\overrightarrow{y} }{ \partial \overrightarrow{y} }A^T ∂y ∂Ay AT ∂ y → T A y → ∂ y → A y → A T y → \frac{ \partial \overrightarrow{y}^TA\overrightarrow{y} }{ \partial \overrightarrow{y} }A\overrightarrow{y}A^T\overrightarrow{y} ∂y ∂y TAy Ay ATy 二次型相关,若 A A A是对称矩阵那么可以以进一步化为 2 A y → 2A\overrightarrow{y} 2Ay 。 简单证明一下1 这一段证明可能有点抽象建议可以自己写一下感受一下我简单说一下我的想法。可以看到最后对应的每个导数上面是 f f f下面是 y y y先看第一行第一行都是对 y 1 y_1 y1求导从左往右分别是 f 1 f_{1} f1到 f m f_m fm因此m对 y 1 y_1 y1的求导从左往右是 ∂ f 1 ( y → ) ∂ y 1 , ∂ f 2 ( y → ) ∂ y 2 . . . , ∂ f m ( y → ) ∂ y m \frac{ \partial f_1(\overrightarrow{y}) }{ \partial {y_1} },\frac{ \partial f_2(\overrightarrow{y}) }{ \partial {y_2} }...,\frac{ \partial f_m(\overrightarrow{y}) }{ \partial {y_m} } ∂y1∂f1(y ),∂y2∂f2(y )...,∂ym∂fm(y )这其实就对应了 A A A第一列的系数也就是 a 11 , a 21 , . . . , a m 1 a_{11},a_{21},...,a_{m1} a11,a21,...,am1相当于把列转为了行。 还有一个2我没证过也可以尝试一下熟悉熟悉规则。
链式求导
对于一般的标量函数 f ( g ( x ) ) f(g(x)) f(g(x))此时对于 f ( x ) f(x) f(x)和 x x x存在链式求导法则 ∂ f ( g ( x ) ) ∂ x ∂ f ( x ) ∂ g ( x ) ∂ g ( x ) ∂ x \frac{ \partial f(g(x)) }{ \partial x }\frac{\partial f(x)}{\partial g(x)}\frac{\partial g(x)}{\partial x} ∂x∂f(g(x))∂g(x)∂f(x)∂x∂g(x) 但是对于矩阵来说由于矩阵的乘法一般不满足交换律此时链式求导法则如下 ∂ f ( g ( x ) ) ∂ x ∂ g → ∂ x → ∂ f ∂ g → \frac{ \partial f(g(x)) }{ \partial x }\frac{\partial \overrightarrow{g}}{\partial \overrightarrow{x}}\frac{\partial f}{\partial \overrightarrow{g}} ∂x∂f(g(x))∂x ∂g ∂g ∂f 可以看到大求导在后面小求导在前面直接面向自变量的在第一个。 下面给出一个例题假设存在 x k 1 A x k B u k x_{k1}Ax_kBu_{k} xk1AxkBuk令 J x k 1 T x k 1 Jx^T_{k1}x_{k1} Jxk1Txk1J其实是一个控制系统里面的损失函数,计算 ∂ J ∂ u k \frac{ \partial J }{ \partial u_{k} } ∂uk∂J 注意其中有一步也就是 ∂ J ∂ x k 1 \frac{ \partial J }{ \partial x_{k1} } ∂xk1∂J这个利用了二次型转换相当于中间的 A A A就是单位阵 I I I I I I是一定对称的所以结果是 2 x k 1 2x_{k1} 2xk1。
线性回归
原理
首先是第一个模型线性回归。可以先考虑一元线性回归这个是最简单的回归模型表达式为 y w x b ywxb ywxb需要计算的就是两个参数 w , b w,b w,b 用矩阵计算的相对要复杂一点但是扩展性强可以先考虑用标量求导的方式来证明如下参考证明) 这个就计算到了求出导数这一步没算出解析解后面代码是用的梯度下降算的。 下面用矩阵的方式求解这种方法的扩展性比较强多元线性回归的证明也是这样本质一样假设 x , y x,y x,y均为 n × 1 n×1 n×1矩阵假设结果为一列 m × 1 m×1 m×1向量 t t t前面 m − 1 m-1 m−1个为 w w w最后一个是 b b b经过证明其结果为 t ( x T x ) − 1 y T x t(x^Tx)^{-1}y^Tx t(xTx)−1yTx 下面是简单的证明这里还是一元一次线性回归的 最后的 x T y x^Ty xTy是一个 2 × 1 2×1 2×1向量导数是 2 × 1 2×1 2×1的分别对应 w , b w,b w,b的偏导当偏导都为0也就对应了最小值所以本质还是解方程解两个偏导方程。有的证明会把损失函数前面加上 1 / 2 n 1/2n 1/2n,也就是除掉样本数但本质其实是一样的不影响求解。 我一开始其实写错了上面打叉的就是第一次写错的还是后面写代码的时候才发现的所以动手实践一下还是很重要的。 注意损失函数本身是一个标量( 1 × 1 1×1 1×1)其导数是 2 × 1 2×1 2×1的对应了求解的参数矩阵。
代码
一元一次线性回归
首先是基于解方程的思路也就是直接算出数值解对应于上面基于矩阵的证明部分。
def my_lr(x,y):n x.shape[0]x np.hstack((x,np.ones((n,1)))) # 加上一列n# print(x.shape) # (n,2)t (inv(dot(x.T,x))x.Ty) # t(xx^T)^(-1)y^Tx, 表示矩阵乘法,dot也是矩阵乘法w,bt[0][0],t[1][0]# 下面进行拟合y_hat w*x[:,0]b # 取第一列,第二列是增广的return y_hat上面的本质在解方程下面用梯度下降的思路求解 简单概括梯度下降就是 w w − α ▽ ww-α▽ ww−α▽其中 ▽ ▽ ▽是梯度多元情况的导数 ▽ ∂ J ∂ w ▽\frac{ \partial J}{ \partial w } ▽∂w∂J α α α是学习率。 首先是基于标量方式求解各个变量的也就是用 w , b w,b w,b这样的形式写出来然后求梯度迭代
def my_gradi_lr(x,y):print(--------------基于标量的梯度下降--------------)# 定义损失函数:T(w,b)1/(2m)*∑(yi_hat-yi),yi_hatwxi-b# ww-α(1/m*∑xi(yi_hat-yi)) bb-α(1/m*∑(yi_hat-yi))m x.shape[0]w,b0,0alpha0.01epochs50for i in range(epochs):y_hat w*xb # 计算拟合值w(w-alpha*(x.T(y_hat-y))/m)# ∑xi(yi_hat-yi)/mb(b-alpha*(np.sum(y_hat-y))/m) # ∑(yi_hat-yi)/my_hat w * x b # 更新新y_hat计算lossloss 1/(2*m)*(y_hat-y).T(y_hat-y) # 1/(2m)*∑(yi_hat-yi)print(epoch:%d loss:%.4f%(i1,loss))# print(w,b)return w*xb推导的过程就是链式求导因为是完全标量的还是比较简单的。 下面是基于矩阵的推导过程跟前面矩阵求解里面是一致的如果一定要比较的话上面那种标量直观一点但是参数多了就不方便了下面的适合更多参数比如之后应该会提到的多元线性回归。
def my_mat_gradi_lr(x,y):print(--------------基于矩阵的梯度下降--------------)n x.shape[0]x np.hstack((x,np.ones((n,1)))) # 加上一列nalpha0.01epochs50tnp.zeros((2,1)) # 系数矩阵for i in range(epochs):gradi 2*(x.Txt-x.Ty)/(2*n)tt - alpha * gradiloss (t.Tx.Txty.Ty-2*y.Txt)/(2*n)print(epoch:%d loss:%.4f % (i 1, loss))return xt # y_hat这两种是完全一样的具体可以看迭代的过程loss的变化是一模一样的。 下面是运算的结果sklearnsklearn model的可以视为标答数值解my model的答案是和sklearn的是一样的说明结果是对的。两个grendient都是对应的梯度下降的解法两种迭代的算法是一样的从迭代过程中的loss可以看出可以运行一下代码由于迭代次数次数比较少50所以和精确解还是有一定距离。 下面是代码
# 一元一次线性回归import numpy as np
from matplotlib import pyplot as plt
from numpy import dot
from numpy.linalg import inv
from sklearn.linear_model import LinearRegression
plt.rcParams[font.sans-serif] [SimHei]
plt.rcParams[axes.unicode_minus] Falsedef my_lr(x,y):n x.shape[0]x np.hstack((x,np.ones((n,1)))) # 加上一列n# print(x.shape)t (inv(dot(x.T,x))x.Ty) # t(xx^T)^(-1)y^Tx, 表示矩阵乘法,dot也是矩阵乘法# print(t.shape)# print(t)w,bt[0][0],t[1][0]# print(w,b)# 下面进行拟合y_hat w*x[:,0]b # 取第一列,第二列是增广的return y_hatdef skle_lr(x,y):model LinearRegression()model.fit(x, y)return model.predict(x)# 利用梯度下降求解(基于标量链式求导)
def my_gradi_lr(x,y):print(--------------基于标量的梯度下降--------------)# 定义损失函数:T(w,b)1/(2m)*∑(yi_hat-yi),yi_hatwxi-b# ww-α(1/m*∑xi(yi_hat-yi)) bb-α(1/m*∑(yi_hat-yi))m x.shape[0]w,b0,0alpha0.01epochs50for i in range(epochs):# print(epoch:,i1)y_hat w*xb # 计算拟合值w(w-alpha*(x.T(y_hat-y))/m)# ∑xi(yi_hat-yi)/mb(b-alpha*(np.sum(y_hat-y))/m) # ∑(yi_hat-yi)/my_hat w * x b # 更新新y_hat计算lossloss 1/(2*m)*(y_hat-y).T(y_hat-y) # 1/(2m)*∑(yi_hat-yi)print(epoch:%d loss:%.4f%(i1,loss))# print(w,b)return w*xb# 采用矩阵形式的梯度下降
def my_mat_gradi_lr(x,y):print(--------------基于矩阵的梯度下降--------------)n x.shape[0]x np.hstack((x,np.ones((n,1)))) # 加上一列nalpha0.01epochs50tnp.zeros((2,1)) # 系数矩阵for i in range(epochs):# print(x.T.shape)# print(t.shape)# print((x.Txt).shape)gradi 2*(x.Txt-x.Ty)/(2*n)tt - alpha * gradiloss (t.Tx.Txty.Ty-2*y.Txt)/(2*n)print(epoch:%d loss:%.4f % (i 1, loss))return xt # y_hatif __name__ __main__:xnp.array([1,2,3,4,5,6]).reshape(-1,1)ynp.array([2,4,7,9,10,11]).reshape(-1,1) # 列向量my_y_hat my_lr(x,y)skle_y_hat skle_lr(x,y)my_gradi_y_hat my_gradi_lr(x,y)my_mat_gradi_y_hat my_mat_gradi_lr(x,y)plt.scatter(x,y)plt.plot(x,my_y_hat,linestyle--,labelmy model)plt.plot(x,my_gradi_y_hat,linestyle-.,labelmy grendient)plt.plot(x,my_mat_gradi_y_hat,linestyle-.,labelmy mat grendient)plt.plot(x,skle_y_hat,linestyle-,labelsklearn model)plt.legend()plt.show()
多元一次线性回归
对于多元线性回归其实本质上还是一样的下面提供多元线性回归的代码使用的例子是一个三元一次线性回归。具体的证明和上面矩阵证明的过程是完全一致的解析求解的方式求逆的代码和一元的一样用梯度下降方式的就改一点点就可以其实一元也就是特殊的多元下面这段代码也可以用。
# 采用矩阵形式的梯度下降
def my_mat_gradi_lr(x,y):print(--------------基于矩阵的梯度下降--------------)n x.shape[0]n_samples x.shape[1]x np.hstack((x,np.ones((n,1)))) # 加上一列nalpha0.01epochs50tnp.zeros((n_samples1,1)) # 系数矩阵,唯一变化的地方for i in range(epochs):gradi 2*(x.Txt-x.Ty)/(2*n)tt - alpha * gradiloss (t.Tx.Txty.Ty-2*y.Txt)/(2*n)print(epoch:%d loss:%.4f % (i 1, loss))return xt # y_hat这里因为不方便画出图像就画出y的比较了每个点都是实际点线对应的都是拟合值x轴表示第几个点。 从结果中看解析解和sklearn得到是一样的就是精确解梯度下降还是差的比较多是因为迭代次数少50提高到5000或者改变学习率alpha就可以得到更好的结果。 贴出代码
# 多元线性回归
import numpy as np
from matplotlib import pyplot as plt
from numpy import dot
from numpy.linalg import inv
from sklearn.linear_model import LinearRegressionplt.rcParams[font.sans-serif] [SimHei]
plt.rcParams[axes.unicode_minus] Falsedef my_lr(x,y):n x.shape[0]x np.hstack((x,np.ones((n,1)))) # 加上一列n# print(x.shape)t (inv(dot(x.T,x))x.Ty) # t(xx^T)^(-1)y^Tx, 表示矩阵乘法,dot也是矩阵乘法# 下面进行拟合print(my lr:,t)return xtdef skle_lr(x,y):model LinearRegression()model.fit(x, y)print(sklearn:,model.coef_)return model.predict(x)# 采用矩阵形式的梯度下降
def my_mat_gradi_lr(x,y):print(--------------基于矩阵的梯度下降--------------)n x.shape[0]n_samples x.shape[1]x np.hstack((x,np.ones((n,1)))) # 加上一列nalpha0.01epochs50tnp.zeros((n_samples1,1)) # 系数矩阵for i in range(epochs):gradi 2*(x.Txt-x.Ty)/(2*n)tt - alpha * gradiloss (t.Tx.Txty.Ty-2*y.Txt)/(2*n)# print(epoch:%d loss:%.4f % (i 1, loss))print(matrix grandient:,t.T)return xt # y_hatif __name__ __main__:xnp.array([[1,2,3,4,5,6],[2,3,5,7,8,9],[8,7,5,4,2,1]]).Tynp.array([2,4,7,9,10,11]).reshape(-1,1) # 列向量my_y_hat my_lr(x,y)skle_y_hat skle_lr(x,y)my_mat_gradi_y_hat my_mat_gradi_lr(x,y)plot_x np.arange(y.shape[0])plt.scatter(plot_x,y)plt.plot(plot_x,my_y_hat,linestyle--,labelmy model)plt.plot(plot_x,my_mat_gradi_y_hat,linestyle-.,labelmy mat grendient)plt.plot(plot_x,skle_y_hat,linestyle-,labelsklearn model)plt.legend()plt.show()❀❀❀完结撒花❀❀❀
