从小学4年级的数学课开始解释线性回归

从小学4年级的数学课开始解释线性回归

flyfish

1.1 回归的由来

英国的统计学家高尔顿（Francis Galton）用统计方法研究两个变量之间关系。他研究父母身高与子女身高之间的关系，“回归”这个词由他引入的。大意是一般情况下爹妈高，孩子高；爹妈矮，孩子矮。这样拟合出了很好看的线性关系。但是，这里是有但是的，爹妈矮，孩子也不会一直比上辈矮下去。爹妈高，孩子也不会一直比上辈高，要不到处都是巨人。好像大自然有一种力量让这些人的身高朝着平均身高去长，即有“回归”到平均数的趋势。回归一词就是这么来的，今天线性回归比以前的含义要丰富些。
以预测房屋的价格（百万元）作为例子来解释线性回归。有一个开放的数据集uci-housing数据集，这是学习线性回归时用到的经典数据集，用来预测波士顿的房价。为了简单起见，这里假设价格只取决于房屋的面积（百平方米），当然我们知道价格不仅仅取决一个因素，还有它所在的城市是一线还是二线，房屋周边的地铁，医院等等等等。

1.2 模型

1.2.1 直线

先从线说起，主要是最简单的直线
小学数学四年级上册对直线非常直观的描述
线段、直线、射线
一根拉紧的线， 绷紧的弦都可以看作线段。 线段有两个端点。
把线段向两端无限延伸，就得到一条直线。
直线没有端点，是无限长的。

关于直线我们看《几何原本》是怎么说的
1、点是没有部分的
2、线只有长度而没有宽度
3、一线的两端是点
4、直线是它上面的点一样地平放着的线
这是一种非常让人费解的描述。
虽然上面一段话费解，但人们对《几何原本》赞誉有加的，其中重要的一个原因是欧几里得在书中创造了公理化的方法，而且这种方法是史无前例的
那么什么是公理化思想

对于什么是公理化思想，小平邦彦的《几何世界的邀请》书中是这样描述的
在一个数学理论体系中，我们尽可能少地选取原始概念和不加证明 的一组公理，以此为出发点，利用纯逻辑推理的规则，把该理论体 系建立成一个演绎系统，这样一种构建理论体系的思想就是公理化 思想。

按照小平邦彦的观点，现代人比前人厉害，原因之一是公理化思想。

1.2.2 平面直角坐标系

由于笛卡尔天才的想法，诞生了解析几何。
为了严谨的描述一个共同的想象。用我们常见的有x轴和y轴的平面直角坐标系，这个就是法国的笛卡尔创立的笛卡尔坐标（Cartesian Coordinates）这种坐标系统包含我们常说的平面直角坐标系。笛卡尔将平面几何图形用代数的方法来表示，直线与方程就有了关系。
我们把我们要定义的直线放在坐标系里面，代数中的二元一次方程在平面坐标上就是一根直线，那就可以用代数去解决几何问题
解析几何也就是这么来的。
他的思路是任何问题-》数学问题-》代数问题-》问题求解
就像我国著名数学家华罗庚曾说的：“数缺形时少直观，形少数时难入微；数形结合百般好，隔离分家万事休”。
解析几何的基本思想是在平面上加入“坐标”的概念，平面上的点和有序实数对(x,y)之间建立了一一对应的关系。每一对实数(x,y)都对应于平面上的一个点；反之，每一个点都对应于它的坐标(x,y) 。一句话就是代数、几何是一家，代数即几何，几何即代数，解决问题时根据需要选择。
研究解析几何还有一个人是费马，他和笛卡儿研究解析几何的方法是大相径庭的，费马主要继承了希腊人的思想，笛卡儿则是从批判希腊的传统出发。

按照笛卡尔的方式继续
我们这里重点用的是斜截式方程表示直线
首先有一个数轴

再加一个数轴

两个数轴变换一下

加上X,Y,就变成了平面直角坐标系，再把一条直线放进去，于是就有了一个直线方程$y=kx+b$，笛卡尔说要有直线方程，于是就有了直线方程。

1.2.3 直线的斜截式方程

直线方程 $y=kx+b$

1.2.4 直线的斜截式方程中的斜率 $k$

$$k=\tan \theta =\frac{y_2-y_1}{x_2-x_1}=\frac{\Delta y}{\Delta x}=\frac{2}{1}=2$$

1.2.5 二元一次方程

直线的方程都是关于$x,y$的二元一次方程
$$\mathrm{A}\mathrm{x}+\mathrm{B}\mathrm{y}+\mathrm{C}=0(\left(A^2+B^2\ne 0\right))$$

二元一次方程中的元，次，方程是怎么回事？
曾经有等式 $1+2=3$，后来有了未知数就有了方程，古埃及人写在纸草书上的数学问题涉及了含有未知数的等式，这时候的描述还是一般的语言文字描述的。现在看到 $1+x=3$这么简单的方程还需要在韦达，笛卡尔的工作基础之上做的改进。
方程是含有未知数的等式,方程 英语 叫 equation。
只含有一个未知数，并且未知数的指数是一次，这样的方程是一元一次方程
例如一元三次方程
$a{x}^3+b{x}^2+cx+d=0（a，b，c，d为常数，x为未知数，且a\ne 0）$
一元表示只含有一个未知数，几元就是几个未知数,我国金代数学家李冶的《测圆海镜》中说：“立天元一为某棠”天元就是现在的未知数$x$，所以用元代表未知数。立天元一为某某就是设某某为未知数$x$。
三次表示未知数的指数，如果有多个未知数就选最大的那个。
阿拉伯的花拉子密的著作《代数学》中有详细的关于方程的事，《代数学》一门也学科的名字就使用他著作的名字，可见他的伟大。

1.2.6 模型的定义

这里假设价格只取决于房屋的面积（平方米）
$y=kx+b$
$y$表示价格，是因变量
$x$表示房屋面积，是自变量。
只有一个自变量的情况称为简单线性回归（simple linear regression）或者叫简单回归。
大于一个自变量情况的叫做多元回归。

直线方程的字母要变化
房屋的面积为$x$​，价格为$y$。基于输入$x$来计算输出$y$的表达式，线性回归假设输出与各个输入之间是线性关系,这个表达式可以叫做假设函数，预测函数 或者 模型（model）
$\hat{y}=wx+b$

其中 $w$ 是权重（weight）说的系数或斜率也是它，$b$ 是偏差（bias）。它们是线性回归模型的参数（parameter）。
$y$表示真实的价格。
$y$带上帽子就是我们估计的价格即通过模型预测的价格$\hat{y} $。模型输出 $\hat{y}$是线性回归预测的价格。是我们猜测的价格。

真实价格 $y$ 和 预测的价格 $\hat{y}$ 有一定误差 $ϵ$ 或者 用 $e$ 表达。

如果我们认为房屋的价格与3种因素相关，式子这样写
$\hat{y}=w_1{x}_1+w_2{x}_2+w_3{x}_3+b$
如果我们认为房屋的价格n种因素相关，式子这样写
$\hat{y}=w_1{x}_1+w_2{x}_2+...+w_n{x}_n+b$
这里举的例子是只与一个因素有关，假如 $y=wx+b+e$有这样的数据

x = [1,2,3,4,5]
y = [6,5,7,9,12]

这5个数据点也可以这样表示$(x,y):(1,6)、(2,5)、(3,7)、(4,9)，(5，12)$
即$$(x^{(1)},y^{(1)}),(x^{(2)},y^{(2)}),...(x^{(n)},y^{(n)})$$
可视化下数据

$$\begin{gathered} 6=w\ast 1+b+e \\ 5=w\ast 2+b+e \\ 7=w\ast 3+b+e \\ 9=w\ast 4+b+e \\ 12=w\ast 5+b+e \end{gathered}$$
如果换成字母可以这样写，用上标表示第几个样本
$$\begin{gathered} y^{(1)}=w{x}^{(1)}+b+e^{(1)} \\ {y}^{(2)}=w{x}^{(2)}+b+e^{(2)} \\ {y}^{(3)}=w{x}^{(3)}+b+e^{(3)} \\ {y}^{(4)}=w{x}^{(4)}+b+e^{(4)} \\ {y}^{(5)}=w{x}^{(5)}+b+e^{(5)} \end{gathered}$$

上面关于一栋房屋的记录被称为一个样本（sample），上面是5个样本。
用来预测标签的因素叫作特征（feature）,即$x$
真实价格叫作标签（label）或目标（target）即$y$

给定一个数据集，目的是寻找模型的权重$w$和偏置$b$，
使得根据模型做出的预测的价格和真实价格差不多。

1.3 损失函数

损失函数（loss function）能够量化目标的真实值与预测值之间的差距。
通常选择非负数作为损失，数值越小表示损失越小，损失为0为理想情况。
当样本$i$的预测值为${\hat{y}}^{(i)}$，其相应的真实标签为$y^{(i)}$时，
平方误差可以定义为以下公式：
$$e=\frac{1}{2}(y-\hat{y}{)}^2$$

训练数据中所有样本的误差的和用大写的E来表示
计算如下
$$E=e^{(1)}+e^{(2)}+e^{(3)}+...+e^{(n)}$$

整个成常用的符号如下
$$\begin{array}{rl}E& =e^{(1)}+e^{(2)}+e^{(3)}+...+e^{(n)}\\ & =\sum _{i=1}^n{e}^{(i)}\\ & =\frac{1}{2}\sum _{i=1}^n(y^{(i)}-{\hat{y}}^{(i)}{)}^2\end{array}$$
累加符号$\sum$ 其他名字 Sigma Notation 、Summation Notation,它表示这样的含义

$$\begin{array}{rlr}& \sum _{i=1}^{5}i& =1+2+3+4+5\end{array}$$

常数$\frac{1}{2}$是一个技巧（trick）因为当我们对损失函数求导后常数系数为1。

为了知道模型训练的怎么样，计算在训练集$n$个样本上的损失均值，这时损失函数就变成了均方误差（MSE, mean squared error）
$$\frac{1}{n}\sum _{i=1}^{n}loss(w,b)=\frac{1}{2n}\sum _{i=1}^n{\left(\hat{y}-y\right)}^2=\frac{1}{2n}\sum _{i=1}^n{\left(wx+b-y\right)}^2$$

在训练模型时，找到参数（$w^{\ast },b^{\ast }$），
这组参数能最小化在所有训练样本上的总损失。如下式：

$$w^{\ast },b^{\ast }={\mathrm{argmin}}_{w,b}loss(w,b).$$
$argmin$ 就是argument of the minimum 的英文缩写,使后面这个式子达到最小值时的 $w$ 和 $b$ 的取值。至于是不是最小值是不确定的，期望是全局最小值，结果可能是局部最小值。
均方误差对应了欧式距离也叫欧几里得距离，找一条直线，让所有的点到直线上的欧式距离之和最小。

欧几里得距离可以用勾股定理算
勾股定理 用一句话说是直角三角形的两条直角边的平方和等于斜边的平方。

勾股定理也叫毕达哥拉斯定理（Pythagorean theorem）
《周髀算经》中记述了(3,4,5)这一组勾股数，例如里面有周公问商高，荣方问陈子等就像子贡问曰一样。商高答周公里面就有3,4,5的事。里面还有赵爽的附录。《九章算术》有勾股定理怎么用的事，刘徽用勾股定理求圆周率。
《周髀算经》的开头里说关于周公与商高的对话，大意是
周公问商高天没有梯子可以上去,地也没法用尺子去一段一段的丈量,那么怎样知道天有多高，地有多大呢

商高回答说：“数的产生来源于对方和圆这些形体认识.其中有一条原理：当直角三角形一条直角边-勾等于3,另一条直角边-股等于4的时候,那么它的斜边-弦就必定是5.这个原理是大禹在治水的时候就总结出来的。
那时候的人都已经认为形和数是一回事，有关联的。周公与商高的对话在公元前1100年左右的西周时期，大禹都距今4000多年了。

西方的东西很多来自古希腊，古希腊古老的著作是欧几里得的《几何原本》，再往前就是毕达哥拉斯，在之前就是泰勒斯。西方人的参考来自古希腊所以叫毕达哥拉斯定理，毕达哥拉斯认为数学可以解释世界上的一切事物，同时认为一切真理都可以用比例、平方及直角三角形去反映和证实。但看历史毕达哥拉斯学派是宗教,他害死了发现根号2的学生希帕索斯,跟宗教裁判害死捍卫太阳中心说的布鲁诺是一模一样的,他们研究的很多都是巫术,数学也算其中一部分.
三国时期的赵爽在对《周髀算经》注释中证明了勾股定理,这些是记录在《九章算术》中的.在《九章算术》中看赵爽写那么多字就为了描述勾股定理，还是文言文。后人用字母、数字、符号等记录简单方便了许多。
法国的韦达在研究数学的时候，发现很多资料都是文字叙述的，跟我们看到赵爽用文言文描述勾股定理是一个感觉，太麻烦了，不管描述还是阅读效率很低的。因为韦达觉得麻烦，他就有意识的想解决这个问题，所以他用字母来表示已知数、未知数。他的后辈笛卡尔当然还是法国的笛卡尔用字母表的前几位表示已知数，后几位表示未知数。所以后来的式子abc表示的是已知数，xyz就是用来表示未知数的，如果只有一个未知数那就是x。
《九章算术》里的赵爽发明勾股定理,《几何原本》的毕达哥拉斯也发明了勾股定理, 距今4000多年前的大禹发现的是勾股数,不是证明了勾股定理, 从周公与商高的对话看,商高也没有直面回答周公的问题.猜想商高意思是天高地大没法直接测,需要借助一些人们发现规律去寻找答案,至于天多高,地多大,我现在也不知道.
在发现勾股数的情况下就能应用而不是勾股定理被发现之后才开始应用，古埃及人的金字塔中的墓室就是按照勾股数设计的尺寸。

二维的欧氏距离计算方法

$$\begin{array}{l}d=\sqrt{{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2}\end{array}$$
三维的欧氏距离计算方法
$$\begin{array}{l}d=\sqrt{{\left(x_1-x_2\right)}^2+{\left(y_1-y_2\right)}^2+{\left(z_1-z_2\right)}^2}\end{array}$$
n维的欧氏距离计算方法
$$d(p,q)=\sqrt{{\left(p_1-q_1\right)}^2+{\left(p_2-q_2\right)}^2+\dots +{\left(p_i-q_i\right)}^2+\dots +{\left(p_n-q_n\right)}^2}$$
$$\text{ distance }=\sqrt{\sum _{i=1}^n{\left(p_i-q_i\right)}^2}$$

1.4 优化算法

可视化$L(w)$和$w$之间的关系,这里将$b$赋值为0
结果是个抛物线，找极值的话，一眼就看出来$L(w)$的极值大概在什么位置。
关于函数的这个“凸”，各种书籍资料都没有统一，这个凹凸正好是中国的象形字，外文书籍用的是字母。所以看书的时候要注意书中的定义和上下文，例如在陈宝林《最优化理论与算法》中用的抛物线的开口向上为凸函数和抛物线的开口向下为凹函数，张景中使用的是上凸和下凸，直观无歧义。在他的著作《直来直去的微积分》中描述了林群的不使用极限概念定义导数。

解决这个问题有两种方式一种是解析解（analytical solution） 即误差最小化问题的解可以直接用公式表达出来。还有一种是数值解（numerical solution）即通过优化算法有限次迭代模型参数来尽可能降低损失函数的值。

1.4.1 梯度下降法解决问题

1.4.1.1 找极值

假设有这样的一个函数是$y=(x-5{)}^2+2$，看看是如何找极值的

x轴代表的是我们待学习的参数 $weight$，y轴代表的是损失函数的值（即Loss值），
曲线y代表的是损失函数。我们的目标是希望损失函数的值最小。可以通过求导数的方式，达到一元二次方程的最小值点，使得导数为0即可。
什么是梯度呢？
说梯度之前先白话说极限，这里没有用教科书的严谨的定义，是从直觉上理解极限。

我们先从无穷说起。
怎么理解无穷？
我们的分析是在标准分析范围内，如果是非标准分析又是另一套方法
从$\sqrt{2}$说起
相同的两个数相乘，结果是2

一个数和自己相乘结果是2,那么这个数是多少呢？根据这个问题了解无穷小数是怎样产生的，从而了解无穷的样子。

因为
$1^2=1$，$1小于2$，
$2^2=4$，$4大于2$,
所以
这个小数肯定在$1和2$之间。

如果一一计算是$1.1^2$ ,$1.2^2$,$1.3^2$，一直到$1.9^2$，
因为
$1.4^2=1.96$,$1.96小于2$，
而$1.52=2.25$,$2.25大于2$。
所以肯定在$1.4和1.5$之间。
所以小数展开一定是以$1.4$开头

现在假设已经以此方法算出的数字是1.4142135。那么可以继续往下计算

$$\begin{array}{l}1.41421350^2=1.9999998236822500\\ 1.41421351^2=1.9999998519665201\\ 1.41421352^2=1.9999998802507904\\ 1.41421353^2=1.9999999085350609\\ 1.41421354^2=1.9999999368193316\\ 1.41421355^2=1.9999999651036025\\ 1.41421356^2=1.9999999933878736\\ 1.41421357^2=2.000000022\end{array}$$

这样我们得到下一位数字是6。

这样的计算永远无法结束

是什么让我们能够相信存在这样的一个无穷小数 $X^2=2$呢？
我们下面的规律

$$\begin{array}{c}{1}^2=1\\ 1.4^2=1.96\\ 1.41^2=1.9881\\ 1.414^2=1.999396\\ 1.4142^2=1.99996164\\ 1.41421^2=1.9999899241\\ 1.414213^2=1.999998409469\\ 1.4142135^2=1.99999982368225\\ 1.41421356^2=1.9999999933878736\end{array}$$

先确定个概念，因数和积，用于我们后面的表达。
10（因数） ×20（因数） = 200（积）
我们根据上面的列表看到因数的小数位数越多，结果积的小数点后面的数字9就越多。因此，如果把该列表永远的写下去，那么后面的9有无穷多位。这样我们感知到了无穷。
那么问题又来了。
1.99999999…（9无限循环）等于2吗？
极限来了。
$\sqrt{2}$表示为有理数组成的无穷数列，不只一种方法表示。
不足近似值构成的无穷数列
{ 1.4 , 1.41 , 1.414 , … } \{1.4,1.41,1.414,…\} {1.4,1.41,1.414,…}
过剩近似值构成的无穷数列
{ 1.5 , 1.42 , 1.415 , … } \{1.5,1.42,1.415,…\} {1.5,1.42,1.415,…}
这里的过剩和不足从词义理解就行。
我们再来看$1.99999999\dots 和2$的问题
我们可以看成数列
{ 1.9 , 1.99 , 1.999 , … } \{1.9,1.99,1.999,…\} {1.9,1.99,1.999,…}

{ 2 , 2 , 2 , 2 , … } \{2,2,2,2,…\} {2,2,2,2,…}
这两个数列都收敛于2

于是我们可以这样写 { 1.9 , 1.99 , 1.999 , … } = { 2 , 2 , 2 , 2 , … } \{1.9,1.99,1.999,…\}=\{2,2,2,2,…\} {1.9,1.99,1.999,…}={2,2,2,2,…}
简写
$$1.99999999\dots =2$$
$1.99999999\dots$看作有理数列$1.9,1.99,1.999,\dots$，其极限是2
$2,2,2,2,\dots$的极限也是2,他们两个叫等价类。
柯西将极限的定义引入了微积分，很创新。因为微积分的创造者们没有把一件事说明白，他用新的概念把这件事说明白了。问题变了，思维方式就变了，之前的人的问题是什么什么是多少，他的问题是如何定义“什么什么”的问题。他定义的框架把微积分的严谨度向前迈了一大步，所以问什么样的问题最重要。微积分的现代体系就是他建的。
牛顿和莱布尼茨没有说明白的微积分的底层也引发了数学危机。
从芝诺用不是数学的语言-大白话描述它，到现在都2500多年了。19世纪末期清代的李善兰把这个知识引入到中国。别人是先进的，那就无论好坏对错全盘吸收，然后持续积累改进迭代。如果后人把问题再变了，那又会离“微积分真理”又进了一步。人的语言可以表达不存在的东西，因相信而这不存在东西也就变得存在。

说梯度的时候，还得理解方向导数，偏导数。
最小化 $f(x)$的问题，如果是一元的，一个 导数的概念就够了；如果是多元的，就出现了偏导数，方向导数，梯度的概念。
导数精确描述了函数变化率，变化率可理解为变量的变化“快慢”问题。研究变化率的问题之前的 $y=kx+b$,就一个未知数$x$，$x$也叫自变量，可以说一元函数，研究多元函数的时候就出现了偏导。偏可以理解成部分，多元就是一个自变量固定，在编程里叫常量。三维可以可视化易于理解，超过三维在我看来只能逻辑推理，所以用一元，二元来可视化理解。假设一个二元函数的偏导数，它反应的是函数沿坐标轴方向的变化率。一个点画线的时候可以朝着坐标轴方向，也可以朝着其他方向，其他方向就成了方向导数.
用杨立昆在《科学之路》白话的说法是 “将其他变量视为常数，对某一个变量进行求导时，我们会得到偏导数。偏导数计算的是函数在这个变量方向上的斜率。由所有方向上的偏导数形成的向量就是梯度。
红叉表示梯度下降的起点
$x=0.5050505050505051$
$y=22.204570962146718$

理想情况是随机找一个点，走两步就到了全局最小值。

每个参数的迭代过程
$$\begin{array}{rl}w& \leftarrow w-\eta \sum _i\frac{\partial L^{(i)}(w,b)}{\partial w}=w-\eta \sum _i{x}^{(i)}\left(x^{(i)}w+b-y^{(i)}\right),\\ b& \leftarrow b-\eta \sum _i\frac{\partial L^{(i)}(w,b)}{\partial b}=b-\eta \sum _i\left(x^{(i)}w+b-y^{(i)}\right).\end{array}$$

如果再加一个参数
$$\begin{array}{rl}{w}_1& \leftarrow w_1-\eta \sum _i\frac{\partial L^{(i)}(w_1,w_2,b)}{\partial w_1}=w_1-\eta \sum _i{x}_1^{(i)}\left(x_1^{(i)}{w}_1+x_2^{(i)}{w}_2+b-y^{(i)}\right),\\ w_2& \leftarrow w_2-\eta \sum _i\frac{\partial L^{(i)}(w_1,w_2,b)}{\partial w_2}=w_2-\eta \sum _i{x}_2^{(i)}\left(x_1^{(i)}{w}_1+x_2^{(i)}{w}_2+b-y^{(i)}\right),\\ b& \leftarrow b-\eta \sum _i\frac{\partial L^{(i)}(w_1,w_2,b)}{\partial b}=b-\eta \sum _i\left(x_1^{(i)}{w}_1+x_2^{(i)}{w}_2+b-y^{(i)}\right).\end{array}$$
梯度下降的代码实现

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(0, 10, 100)   #从0到6选取100个点
y = ( x -5) ** 2 + 2 #我们自己定义的损失函数
plt.scatter(x[5], y[5], color='r',marker='x')#我们随机选择的起点
print(x[5], y[5])#0.5050505050505051 22.204570962146718
plt.plot(x, y,color='g')

plt.xlabel('weight', fontsize=10)
plt.ylabel('loss function', fontsize=10)
plt.show()

#定义损失函数
def L(weight):
    return (weight-5)**2 + 2

#损失函数的梯度，因为是单一变量，这里可以说是损失函数的导数
def dL(weight):
    return 2 * (weight - 5)

#初始点
weight = x[5]
weight_record = [weight]

#超参数-步长，也叫学习率
eta = 0.1

#在计算机中处理关于0的问题，由于精度或者学习率的设置无法得到导数的结果是0,只能是近似0.
epsilon = 1e-6
while True:
    gradient = dL(weight)                     #求导数
    last_weight = weight                       #上一个weight的值
    weight = weight - eta * gradient           #得到一个新的weight
    weight_record.append(weight)
    if(abs(L(weight) - L(last_weight)) < epsilon):
        break                                #处理0的问题终止循环
plt.plot(x,L(x),color="g")
plt.plot(np.array(weight_record),L(np.array(weight_record)),'-.',color='r',marker='o')
plt.show()
print(len(weight_record)) #38

那么梯度下降算法怎么实现呢？通过上面的代码总结如下：

1. 函数1：待优化的函数$f(x)$，即代码中的$L(weight)$，它可以根据给定的输入返回函数值。
2. 函数2：待优化函数的导数$g(x)$，即代码中的$dL(weight)$，它可以根据给定的输入返回函数的导数值。
3. 变量x：保存当前优化过程中的参数值，优化开始时该变量将被初始化成某个数值，优化过程中这个变量会不断变化，直到它找到最小值。
4. 变量grad：保存变量$x$点处的梯度值。即代码中的$gradient$
5. 变量step：表示沿着梯度下降方向行进的步长，也被称为学习率,即代码中的$eta$。

整个程序运行如下，当$x=5$时找到了最小值，实际计算$x=4.998833044392905$找到了最小值,共找了38步。

1.4.1.2 学习率

为什么学习率是我们要调整的参数呢？
我们把学习率设置成不同的值，看图就能找到原因
当学习率eta = 0.001 ，计算步骤是2824步，你的机器要浪费资源了。
看点多的都连成线了

当学习率eta = 0.9 ，计算步骤是38步

一句话更新过程
$\nabla$梯度符号，$\eta$ 表示学习率
$${\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}_{new}={\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t}}_{old}-\eta \nabla L(\mathrm{w}\mathrm{e}\mathrm{i}\mathrm{g}\mathrm{h}\mathrm{t})$$
梯度下降可以根据训练时每次迭代更新参数采用的样本数分为
批量梯度下降（Batch Gradient Descent）
随机梯度下降（Stochastic Gradient Descent）
小批量梯度下降（Mini-Batch Gradient Descent）

批量梯度下降法：每次迭代中使用所有的样本
随机梯度下降法：每次迭代中使用1个样本
小批量梯度下降法：一共有n个样本，每次迭代中使用b个样本
线性回归完整代码

import numpy as np
import matplotlib.pyplot as plt

X_train = np.array([1,2,3,4,5])
y_train = np.array([6,5,7,9,12]).reshape(-1, 1)
data_X = []
for x in X_train:
    data_X.append([1,x])
data_X = np.array((data_X))

n,feature_count = np.shape(data_X) # n, 样本数 feature_count: 特征数
max_iter = 1000 # 迭代数
weights = np.ones((feature_count,1))
alpha = 0.1 # 学习率
for i in range(0,max_iter):
    error = np.dot(data_X,weights)- y_train
    gradient = data_X.transpose().dot(error)/n
    weights = weights - alpha * gradient
print("w:",weights[1:][0]) # 输出模型参数 w
print("b",weights[0]) # 输出参数 b

X_test = np.arange(0, 7, 1).reshape(-1, 1)
plt.plot(X_test, X_test*weights[1][0]+weights[0][0], color='green')
plt.scatter(X_train,y_train,color='red')

x=X_train.tolist()
y=y_train.reshape(-1).tolist()

for i in range(0,len(X_train)):
    plt.text(int(X_train[i]),int(y_train[i]),str((x[i],y[i])), wrap=False)
plt.show()

可视化结果

线性回归的结果是

系数 w：[1.60000002]
截距 b： [2.99999993]

1.4.2 代数法解决问题

下面包含了3种不同的解法

1.4.2.1 代数法实现1

分别对参数$(w,b)$分别求偏导数，并令其等于0，得到一个关于$(w,b)$的二元方程组。求解这个方程组就可以得到$(w,b)$的值.

求导后得到以下方程组:

$$\frac{\partial L_{(w,b)}}{\partial w}=2\left(w\sum _{i=1}^n{x}_i^2-\sum _{i=1}^n\left(y_i-b\right)x_i\right)=0$$
$$\frac{\partial L_{(w,b)}}{\partial b}=2\left(nb-\sum _{i=1}^n\left(y_i-w{x}_i\right)\right)$$

这个方程组是一个二元一次方程组

求解过程$$\frac{\partial L_{(w,b)}}{\partial w}=2\left(w\sum _{i=1}^n{x}_i^2-\sum _{i=1}^n\left(y_i-b\right)x_i\right)$$
[过程]：已知$L_{(w,b)}=\sum _{i=1}^n{\left(y_i-w{x}_i-b\right)}^2$，所以
$$\begin{array}{rl}\frac{\partial L_{(w,b)}}{\partial w}& =\frac{\partial }{\partial w}\left[\sum _{i=1}^n{\left(y_i-w{x}_i-b\right)}^2\right]\\ & =\sum _{i=1}^n\frac{\partial }{\partial w}\left[{\left(y_i-w{x}_i-b\right)}^2\right]\\ & =\sum _{i=1}^n\left[2\cdot \left(y_i-w{x}_i-b\right)\cdot (-x_i)\right]\\ & =\sum _{i=1}^n\left[2\cdot \left(w{x}_i^2-y_i{x}_i+b{x}_i\right)\right]\\ & =2\cdot \left(w\sum _{i=1}^n{x}_i^2-\sum _{i=1}^n{y}_i{x}_i+b\sum _{i=1}^n{x}_i\right)\\ & =2\left(w\sum _{i=1}^n{x}_i^2-\sum _{i=1}^n\left(y_i-b\right)x_i\right)\end{array}$$

求解过程
$$\frac{\partial L_{(w,b)}}{\partial b}=2\left(nb-\sum _{i=1}^n\left(y_i-w{x}_i\right)\right)$$
[过程]：已知$L_{(w,b)}=\sum _{i=1}^n{\left(y_i-w{x}_i-b\right)}^2$，所以
$$\begin{array}{rl}\frac{\partial L_{(w,b)}}{\partial b}& =\frac{\partial }{\partial b}\left[\sum _{i=1}^n{\left(y_i-w{x}_i-b\right)}^2\right]\\ & =\sum _{i=1}^n\frac{\partial }{\partial b}\left[{\left(y_i-w{x}_i-b\right)}^2\right]\\ & =\sum _{i=1}^n\left[2\cdot \left(y_i-w{x}_i-b\right)\cdot (-1)\right]\\ & =\sum _{i=1}^n\left[2\cdot \left(b-y_i+w{x}_i\right)\right]\\ & =2\cdot \left[\sum _{i=1}^{n}b-\sum _{i=1}^n{y}_i+\sum _{i=1}^{n}w{x}_i\right]\\ & =2\left(nb-\sum _{i=1}^n\left(y_i-w{x}_i\right)\right)\end{array}$$

结果得到w和b的解析解
一个字母头上加了横线（bar）就是平均值
计算结果1
$$\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i$$
$$\bar{y}=\frac{1}{n}\sum _{i=1}^n{y}_i$$
$$w=\frac{{\sum }_{i=1}^n{y}_i(x_i-\bar{x})}{{\sum }_{i=1}^n{x}_i^2-\frac{1}{n}({\sum }_{i=1}^n{x}_i{)}^2}$$
$$b=\frac{1}{n}\sum _{i=1}^n(y_i-w{x}_i)$$
或者
$$b=\bar{y}-w\bar{x}$$

代数法完整的代码实现
完整的单纯的代码实现，计算过程不用库，才看的通透。
按照上面的计算结果代码如下

import numpy as np
x= np.array([1,2,3,4,5])
y= np.array([6,5,7,9,12])
n=5
m_x, m_y = np.mean(x), np.mean(y)
s1 = np.sum((x - m_x) * (y))
s2 = np.sum(x ** 2) - np.sum(x) ** 2 / n

w = s1 / s2
b = m_y - w * m_x
print("w: {} b: {}".format(w, b))
#w: 1.6 b: 2.999999999999999

1.4.2.2 代数法实现2

不同的结果有不同的编写代码方式
计算结果2
$$\begin{array}{r}w=\frac{{\sum }_{i=1}^n\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)}{{\sum }_{i=1}^n{\left(x_i-\bar{x}\right)}^2}\end{array}$$
$$\begin{array}{l}{S}_1={\sum }_{i=1}^n\left(x_i-\bar{x}\right)\left(y_i-\bar{y}\right)\\ S_2={\sum }_{i=1}^n{\left(x_i-\bar{x}\right)}^2\\ w=\frac{S_1}{S_2}\\ b=\bar{y}-w\bar{x}\end{array}$$
实现代码如下

import numpy as np

x= np.array([1,2,3,4,5])
y= np.array([6,5,7,9,12])

m_x, m_y = np.mean(x), np.mean(y)
s1 = np.sum((x - m_x) * (y - m_y))
s2 = np.sum((x - m_x) ** 2)

w = s1 / s2
b = m_y - w * m_x
print("w: {} b: {}".format(w, b))
#w: 1.6 b: 2.999999999999999

1.4.2.3 代数法实现3

计算结果3
$$\begin{array}{rl}w& =\frac{n\Sigma x_i{y}_i-\Sigma x_i\Sigma y_i}{n\sum x_i^2-{\left(\Sigma x_i\right)}^2}\\ & \\ b& =\frac{\sum x_i^2\Sigma y_i-\Sigma x_i(\Sigma x_i{y}_i)}{n\sum x_i^2-{\left(\Sigma x_i\right)}^2}\end{array}$$

import numpy as np
x= np.array([1,2,3,4,5])
y= np.array([6,5,7,9,12])
n=5
t1 = 0
t2 = 0
t3 = 0
t4 = 0
for  i in range(5):
    t1 = t1 + x[i] * x[i]
    t2 = t2 + x[i]
    t3 = t3 + x[i] * y[i]
    t4 = t4 + y[i]

w = (t3*n - t2*t4) / (t1*n - t2*t2);
b = (t1*t4 - t2*t3) / (t1*n - t2*t2);
print("w: {} b: {}".format(w, b))

# w: 1.6 b: 3.0

并不是所有的问题都存在解析解，所以在深度学习里并没有广泛应用。

1.4.3 矩阵法解决问题

简单的通过直线方程，建立方程组,然后求方程组的解,就不需要按照有预测值或者估计值方法去做了。

1.4.3.1 方程组变矩阵

先看一个方程组是如何变成矩阵的。
$$\begin{array}{rl}-x+2y-z& =-1\\ -2x+2y-z& =-3\\ 3x-y+z& =8\end{array}$$
三维向量的单个方程式
$$\left(\begin{array}{r}-x+2y-z\\ -2x+2y-z\\ 3x-y+z\end{array}\right)=\left(\begin{array}{r}-1\\ -3\\ 8\end{array}\right)$$
上面的向量方程（vector equation）不省略等于1的系数的时候是这样子的
$$\left(\begin{array}{r}-1x+2y-1z\\ -2x+2y-1z\\ 3x-1y+1z\end{array}\right)=\left(\begin{array}{r}-1\\ -3\\ 8\end{array}\right)$$
方程组可以写成矩阵形式的矩阵方程（matrix equation）
$$\left(\begin{array}{rrr}-1& 2& -1\\ -2& 2& -1\\ 3& -1& 1\end{array}\right)\left(\begin{array}{l}x\\ y\\ z\end{array}\right)=\left(\begin{array}{r}-1\\ -3\\ 8\end{array}\right)$$

用逆矩阵（inverse matrix）形式把$x,y,z$ 全部放在左边

$$\left(\begin{array}{l}x\\ y\\ z\end{array}\right)={\left(\begin{array}{rrr}-1& 2& -1\\ -2& 2& -1\\ 3& -1& 1\end{array}\right)}^{-1}\left(\begin{array}{r}-1\\ -3\\ 8\end{array}\right)$$
$${\left(\begin{array}{rrr}-1& 2& -1\\ -2& 2& -1\\ 3& -1& 1\end{array}\right)}^{-1}=\left(\begin{array}{rrr}1& -1& 0\\ -1& 2& 1\\ -4& 5& 2\end{array}\right)$$
$$\left(\begin{array}{l}x\\ y\\ z\end{array}\right)=\left(\begin{array}{rrr}1& -1& 0\\ -1& 2& 1\\ -4& 5& 2\end{array}\right)\left(\begin{array}{r}-1\\ -3\\ 8\end{array}\right)=\left(\begin{array}{l}2\\ 3\\ 5\end{array}\right)$$
$x=2,y=3,z=5$ 就是方程组的解

1.4.3.2 实践

采用矩阵的方法来解，已知 $x$ 和 $y$,求 $w$ 和 $b$。换成我们自己的字母
$$\{\begin{array}{c}{y}^1=w{x}^1+b\\ y^2=w{x}^2+b\\ y^3=w{x}^3+b\\ \dots \\ y^n=w{x}^n+b\end{array}$$
方程组改为矩阵表示
$$\left[\begin{array}{c}{y}^1\\ y^2\\ y^3\\ ⋮\\ y^n\end{array}\right]=w\left[\begin{array}{c}{x}^1\\ x^2\\ x^3\\ ⋮\\ x^n\end{array}\right]+b\Rightarrow \left[\begin{array}{c}{y}^1\\ y^2\\ y^3\\ ⋮\\ y^n\end{array}\right]=\left[\begin{array}{cc}{x}^1& 1\\ x^2& 1\\ x^3& 1\\ ⋮& \\ x^n& 1\end{array}\right]\left[\begin{array}{c}w\\ b\end{array}\right]$$
符号是$Y=XW$

$$Y=\left[\begin{array}{c}{y}^1\\ y^2\\ y^3\\ ⋮\\ y^n\end{array}\right]$$
$$X=\left[\begin{array}{cc}{x}^1& 1\\ x^2& 1\\ x^3& 1\\ ⋮& \\ x^n& 1\end{array}\right]$$
$$W=\left[\begin{array}{c}w\\ b\end{array}\right]$$
损失函数
$$L(W)=\frac{1}{2}(XW-Y{)}^T(XW-Y)$$
一个字母右上角加了T就是矩阵转置（matrix transpose）
$$(1,2,3,4{)}^T=\left(\begin{array}{l}1\\ 2\\ 3\\ 4\end{array}\right)$$

损失函数对 $W$ 求导
$$\frac{\partial L}{\partial W}=X^T(XW-Y)=0$$

跟上面的一样把要求的未知数放在左边
$$XW=Y\Rightarrow X^{T}XW=X^{T}Y$$
$$W=(X^{T}X{)}^{-1}{X}^{T}Y$$
代码实现，与公式一模一样

import numpy as np
x= np.array([1,2,3,4,5])
y= np.array([6,5,7,9,12])

# 对x（小写）进行扩展 X（大写）
X = np.c_[x,np.ones((len(x), 1))]
print(X)
'''
np.linalg.inv用来求矩阵的逆矩阵
T表示矩阵的转置
'''
W = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
print("w: {} b: {}".format(W[0], W[1]))
# w: 1.5999999999999996 b: 3.0000000000000053

把一个二元一次方程组中的两个方程画出直线，两条直线可能或平行了没有交点（没有解），有一个交点（一个解）或有无数个交点（无数解）

1.5 预测

通过“已学习”的线性回归模型$\hat{y}=wx+b$来预测。
这样的过程叫预测（prediction）或推断（inference）。
系数 w：[1.60000002] 约等于1.6
截距 b： [2.99999993] 约等于3
$y=1.6x+3$ 当 $x=7$ 时 ，模型的输出是 $14.2$

1.6 参考文献

《教育部审定2013·义务教育教科书 数学（四年级，上册）》人民教育出版社
《普通高中课程标准实验教科书·数学2（必修）》《直线的方程》 B版 人民教育出版社
《数学史概论（第三版）》5.3 解析几何的诞生 李文林
《几何原本》
《几何世界的邀请》小平邦彦
《机器学习 [Machine Learning]》西瓜书 周志华
《深度学习》花书 [deep learning] Ian Goodfellow、Yoshua Bengio和Aaron Courville
《最优化理论与算法》陈宝林
《直来直去的微积分》张景中
《机器学习公式详解》 南瓜书 谢文睿,秦州 编
《科学之路：人、机器与未来》杨立昆
《高等数学 第七版 下册》《方向导数与梯度》
《数学悖论与三次数学危机》 韩雪涛