二、【机器学习作业】多变量线性回归算法（python版ex1）

学习完了机器学习的多变量线性回归课程，现在来将所学记录下来。

多变量线性回归（Linear Regression with multiple variables）

（1） 多维特征 & 多变量梯度下降

概念：
多维特征，顾名思义就是指有多个特征值或者多个变量组成的。与单变量不同的就是多了一些特征量而已。因此之前所述的数学公式只需要修改一下，即：

• 多变量的假设函数：

$$h_{\theta }(x)={\theta }_0+{\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_n{x}_n$$简化为：$$h_{\theta }(x)={\theta }^{T}X$$ 其中，$\theta$和$X$都是$n+1$维的向量，$\theta =\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right]$ ，$X=\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_n\end{array}\right]$ 。为了方便记忆，一般令$x_0=1$，通俗来说就是在已有的 n 个特征外，我们额外增加了一个 $x_0$ 特性，但这个$x_0$总是为1。

• 多变量的代价函数：
  （用于测量或者度量算法性能的好坏）
  $$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2\cdot x_j^{(i)}$$ 其中，$(h_{\theta }(x^{(i)})-y^{(i)})$是一个实数，$x_j^{(i)}$是一个向量。

• 多变量的梯度下降函数：
  （梯度下降算法不仅可以在线性回归中使用，还可以在其他最小化的计算中）
  $${\theta }_j:={\theta }_j-\alpha \frac{\partial }{\partial {\theta }_j}J(\theta )$$
  当$n=1$时，$${\theta }_0:={\theta }_0-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})$$ $${\theta }_1:={\theta }_1-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x^{(i)}$$ PS：需要同时更新${\theta }_0$和${\theta }_1。$
  当$n\ge 1$时，$${\theta }_j:={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x_j^{(i)}$$ PS：需要同时更新${\theta }_j$，j=0,…,n。

这两种梯度下降算法的思路是类似的。下面是当n=2时的例子。

变量说明：

$n$ = 特征的数量
$x^{(i)}$=第$i$ 个训练样本的输入特性值
$x_j^{(i)}$=在第$i$个特征中第$j$个特征值

（2）多维梯度下降的例子：特征缩放（Feature Scaling）

有两种常见的特征缩放：

1. 将特征除以最大值。
2. 均值归一化。

不同的特征向量值，在相近的范围内，梯度下降算法能更快的收敛。

在下面这个例子中，个人的理解是：$x_1$和$x_2$的特征值相差比较大的时候，根据代价函数绘制出的等高线图会又扁又长，在这个情况下计算梯度下降的话就会很长的时间才能达到局部的最小值，而且过程中波动的次数及角度会对结果产生影响。当采取特征缩放后，根据代价函数绘制出的等高线图看起来更加圆一些，梯度下降算法计算过程不会很久，而且来回波动的频率不会很高，那么算法的运行速率就会很快，就能快速找到局部最小值。

1. 将特征除以最大值：
   
   PS：尽量使$x_1$和$x_2$的范围定在[-1,1]中，也可以更大一些，比如[-3,3]、[-2,0.5]、[-1/3,1/3]，尽量使范围适中即可。但是不能够很大，如[-100,100]，也不能太小，如[-0.0001,0.0001]。

2. 均值归一化：
   
   PS：均值归一化后的结果，会使得 $x_i$ 的范围接近[-0.5,0.5]，均值归一化的目的和上面的方法一样，都是为了让梯度下降的运行速度更快一些。

（3）学习率（Learning rate）

$${\theta }_j:={\theta }_j-\alpha \frac{\partial }{\partial {\theta }_j}J(\theta )$$ 如何取到合适的学习率$\alpha$，在单变量线性回归中也提到过。
（参考 https://blog.csdn.net/mary_0830/article/details/96423474）


个人的理解是： 根据图中的这条下降的曲线可以判断出，梯度下降算法是否已经收敛。有一个问题是，需要多少次才知道已经收敛了呢？ 这就产生了一种算法，使用自动收敛测试算法可以告诉你梯度下降是否已经收敛了。但是选择一个适合的阈值 ε 是相当困难的，因此还是偏重选择第一种根据下降曲线判断是否收敛的方法。根据曲线可以知道，梯度下降算法是否正常的工作，上面展示了三种非正常工作的方式，要解决这类情况，可以使用较小的$\alpha$值即可。需要注意的是，当$\alpha$很小的时候，计算梯度下降的最小值就会很慢很慢。

求取合适的学习率α的步骤：
在开始时我们不确定最合适的$\alpha$值是多少，一般我们会选取等比的一系列$\alpha$值做尝试，如0.001，0.003，0.01，0.03，0.1，0.3，1，… 每次尝试间隔大约是3倍左右。
最终比较不同学习率下，算法的收敛速度，选取合适的学习率$\alpha$。

（4）正规方程（Normal equation）

最小二乘法可以将误差方程转化为有确定解的代数方程组（其方程式数目正好等于未知数的个数），从而可求解出这些未知参数。这个有确定解的代数方程组称为最小二乘法估计的正规方程。

这是一种区别于迭代方法的直接解法。

直观的理解是：

1. $J(\theta )$ 是$\theta$的一个函数，且$\theta$为实数。
   
   说明： 求这个函数的最小值则需要对这个函数进行求导，然后把导数置$0$。这是由于，导数就是函数曲线的斜率，当导数为 0 的时候，说明是斜率水平与 x 轴，也就是曲线在最低点的时候。
2. 当$\theta$不是实数时，而是一个$n+1$维的向量。
   
   说明： 逐个对对象参数${\theta }_j$ 求 $J(\theta )$ 的偏导数，然后把它们全部置$0$，求出 ${\theta }_0$，${\theta }_1$，…，${\theta }_n$ 的值，这样就能得到能够最小化代价函数$J(\theta )$的$\theta$值。

原理推导过程：
多变量线性回归代价函数为：
$$J({\theta }_0,{\theta }_1,\cdots ,{\theta }_n)=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2$$
其中：$h_{\theta }(x)={\theta }^{T}X={\theta }_0{x}_0+{\theta }_1{x}_1+\cdots +{\theta }_n{x}_n$，$\theta 、X$和$y$都是矩阵或向量。

正规方程是通过求解下面的方程来找出使得代价函数最小的参数：$$\frac{\partial }{\partial {\theta }_j}J({\theta }_j)$$

要对代价函数 $J(\theta )$中预测值与真实值的差的平方的累加进行求导。

第一步，假设有m个训练实例，每个实例有n个特征，则训练实例集为：$$X=\left[\begin{array}{ccc}{x}_0^{(1)}& \cdots & x_n^{(1)}\\ ⋮& \ddots & ⋮\\ x_0^{(m)}& \cdots & x_n^{(m)}\end{array}\right]$$其中$x_j^{(i)}$表示第i个实例第j个特征。
特征参数为：$$\theta ={\left[{\theta }_0,{\theta }_1,\cdots ,{\theta }_n\right]}^T$$
输出变量为：$$Y={\left[y^{(1)},y^{(2)},\cdots ,y^{(m)}\right]}^T$$
故代价函数为：$$\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2=(h(X)-Y{)}^T(h(X)-Y)$$对代价函数的式子进行消除累加.

第二步，对上面的式子进行化简：$$\begin{gathered} (h(X)-Y{)}^T(h(X)-Y)=（X\ast \theta -Y{)}^T（X\ast \theta -Y) \\ =（{\theta }^T\ast X^T-Y^T)（X\ast \theta -Y) \\ ={\theta }^T{X}^{T}X\theta -{\theta }^T{X}^{T}Y-Y^{T}X\theta +Y^{T}Y \end{gathered}$$

第三步，将化简后的式子，对$\theta$求导：（简化一下过程，直接把对上式的每一项的导数写出来）
第一项：$${\theta }^T{X}^{T}X\theta$$是一个标量，对标量求导得到$X^{T}X\theta$。
该矩阵求导为分母布局下的标量/向量形式。

第二项：同样是对标量求导，得到$X^{T}Y$

该矩阵求导为分母布局下的标量/向量形式。
第三项：得到$X^{T}Y$
该矩阵求导为分母布局下的标量/向量形式。

第四项：导数为0
说明：矩阵推导公式在下面会介绍。

综上，对上式求导的结果是：$-2X^{T}Y+X^{T}X\theta$

为了取到代价函数的最小值，所以让导函数等于零。
就得到了：$$\theta =(X^{T}X{)}^{-1}{X}^{T}Y$$

---

看到上面的求导公式，我刚开始是一脸懵。为什么是这样的呢？原来当我查找相关资料了之后才发现，原来之前学过的求导公式是不一样的。下面进行知识的补充部分。

一般来讲，我们约定$x=(x_1,x_2,...x_n{)}^T$，这是分母布局。这就是说，$X$和$\theta$也属于分母布局。常见的矩阵求导方式有：向量对向量求导，标量对向量求导，向量对标量求导。

分子布局： 分子为 y 或者分母为 x^T (即，分子为列向量或者分母为行向量)
分母布局： 分子为 y^T 或者分母为 x(即，分子为行向量或者分母为列向量)

---

例如，上面第三步中的所有式子都是对$\theta$求导，那么由于$\theta$是列向量，符合分母布局的要求，那么上面的所有式子都是分母布局。

按照不同的布局方式，有几种情形，比如在分子布局方式下计算：标量/向量，向量/标量，向量/向量，标量/矩阵，矩阵/标量。

布局方式	分子布局	分母布局
标量/向量	分母是向量，且是分子布局，则把分母的向量按照行向量铺开	分母是向量，且是分母布局，则把分母的向量按照列向量铺开
向量/标量	分子是向量，且是分子布局，则把分子按照列向量铺开	分子是向量，且是分母布局，则把分子按照行向量铺开
向量/向量	分子分母都是向量，且是分子布局，则分子向量按照列向量铺开，分母向量按照行向量铺开	分子分母都是向量，且是分母布局，则分子向量按照行向量铺开，分母向量按照列向量铺开
标量/矩阵	分子布局下，X矩阵是转置后铺开的	分母布局下，X矩阵无需转置，就是原矩阵

第三步中的分母布局解释清楚了，那么为什么是标量\向量的布局形式呢？

在草稿纸上可以写出，${\theta }^T{X}^{T}X\theta$的矩阵表达式，这里就具体罗列出来了。简单来说，就是${\theta }^T$是1×（n+1）维的向量，$X^T$和$X$都是（n+1）×（n+1）维，$\theta$是（n+1）×1维的向量。根据矩阵的乘法法则，可以计算出，${\theta }^T{X}^{T}X\theta$最后的结果是1×1的形式，即结果是一个标量。第二项和第三项以此类推即可。

下面举出一个简单的例子：（向量对向量求导）

在向量对向量求导中，复合矩阵的求导，对分子布局来说，复合矩阵求导就和复合函数求导类似，从外面开始“扒皮”，一层接一层往里面“扒”；但是对于分母布局来说，复合矩阵是从里面开始“扒皮”的，从里向外开始计算。

在标量对向量的求导中，与复合函数求导法则类似，从外向里面“扒皮”。

下面这个公式可以用来计算第一项。将$X^{T}X$看成是A，下面是分母布局下的标量\向量的形式，那么对应得到的是$2X^{T}X\theta$。

下面这个公式可以用于解释第二项和第三项的结果。将第二项中的$X^{T}Y$看成是$a$，在分母布局下的标量\向量中，第二项的求导后的结果就是$X^{T}Y$；第三项可以将$X^{T}Y$看成是$a^T$（$Y^{T}X=a$，根据转置的公式，就可以得到$a^T=(Y^{T}X{)}^T=X^{T}Y$）就可以得到结果为$X^{T}Y$。

常用的矩阵求导公式：

矩阵求导公式参考 https://blog.csdn.net/daaikuaichuan/article/details/80620518

---

下面使用正规函数的一个例子：

其中，$X$是$m\ast (n+1)$维矩阵，$y$ 是一个 $m$ 维向量。
$m$是样本个数，$n$是特征量个数。

根据上面图中的公式求得的$\theta$值就是使得代价函数$J(\theta )$取到最小的值。

PS：如果你使用“正规方程算法”，那么就不需要特征缩放。（但是使用梯度下降算法，则需要特征缩放，需要让特征值靠近一个区间范围）。

若$X^{T}X$是不可逆的情况下，那产生的原因和解决办法分别是什么呢？

产生的原因：

1. 存在多余的特征，可能存在线性相关的量。
2. 过多的特征（$m⩽n$）。

分别解决的办法：

1. 删除其中一个线性相关的特征量。
2. 减少特征量或者使用一个叫“正则化”的方法。（“正则化”可以让你使用很多的特征来配置很多参数，即使你有一个相对较小的训练集。）

上面两种算法的比较：

说明： 梯度下降算法在特征量$n$很大的情况下，通常还是会使用的。当$n$是百位数或者千位数时，是可以使用“正规方程算法”的；但是如果$n$是万位数或以上时，通常不会使用“正规方程算法”，因为该算法的计算效率就会很慢，为了计算还是会使用梯度下降算法或者其他算法来进行操作。具体地说，只要特征变量数量小于一万，通常使用标准方程法，而不使用梯度下降法。

（5）python代码实现

单变量线性回归算法的代码在上一条博客已经说明了。
（网址： https://blog.csdn.net/mary_0830/article/details/96423474）

下面编写 多变量线性回归 的代码：
（ex1data2.txt 包含的数据有 房子大小、卧室数量和房子的价格）

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
iterations = 1500

def computeCost(X, y, theta):
    sum_all = np.power(((X.dot(theta.T))-y),2) #数组元素求n次方
    J = np.sum(sum_all)/(2*X.shape[0])
    return J

def GradientDescent(X, y, theta): # 初始化一个theta临时矩阵
    #temp = np.array(np.zeros(theta.shape))
     # 参数theta的数量
    #parameters = int(theta.flatten().shape[0])
     # 初始化一个ndarray，包含每次 iterations 的 cost
    #cost = np.zeros((iterations,1))
    alpha = 0.001
    #iterations = 1000
    # cost存放每次修改theta后代价函数的值
    cost=[]
    #样本数量 m
    m = X.shape[0]
    for i in range(iterations): #进一步向量化求解
        temp = theta - (alpha/m) * (X.dot(theta.T) - y).T.dot(X) # theta 是一个行向量，计算得出的结果也是行向量
        theta = temp
        cost.append(computeCost(X, y, theta)) # theta 是未知数，X和y是已知数
    return theta,cost

def normalEqution(X,y):
    theta = np.linalg.inv(X.T@X)@X.T@y   #X.T@X等价于X.T.dot(X)
    #np.linalg.inv()：矩阵求逆矩阵求逆;
    return theta

path = 'ex1data2.txt'
# names添加列名，header用指定的行来作为标题，若原无标题且指定标题则设为None
#读取CSV（逗号分割）文件到 DataFrame
data1 = pd.read_csv(path, header=None, names=['Size', 'Bedrooms', 'Price'])  
#查看数据集的前 5 行
data1.head()
print(data1.head())
#查看数据的分布情况
data1.describe()
print(data1.describe())
#绘制数据的图形
#data.plot(kind='scatter', x='Population', y='Profit', figsize=(8,5))
data1 = (data1 - data1.mean()) / data1.std()
print(data1.head())
data1.head()
# add ones column
data1.insert(0, 'Ones', 1)
iterations = 1000
# set X (training data) and y (target variable)
cols = data1.shape[1]#X.shape[1]图片水平尺寸,X.shape[0]图片垂直尺寸,X.shape[2]图片通道数
X1 = data1.iloc[:,0:cols-1] #通过行号来取行数据
y1 = data1.iloc[:,cols-1:cols]

# convert to matrices and initialize theta
X1 = np.matrix(X1.values) #X1矩阵
y1 = np.matrix(y1.values) #y1矩阵
theta1 = np.matrix(np.array([0,0,0]))
X1.shape, theta1.shape, y1.shape

computeCost(X1,y1,theta1)
final_theta1, cost1 = GradientDescent(X1, y1, theta1)
print('GradientDescent:',final_theta1)
print('GradientDescent_Costfunction:',computeCost(X1,y1,final_theta1))
final_theta2 = normalEqution(X1,y1)
print('normalEqution :',final_theta2)
print('normalEqution_Costfunction:',computeCost(X1,y1,final_theta2.T))

fig, ax = plt.subplots(figsize=(7,5))
ax.plot(np.arange(iterations), cost1, 'y')
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Iteration-Cost')
plt.show()

实验结果：
显示数据中前5行以及数据的分布情况：

特征归一化后：

（PS：若在没有特征归一化的情况下，那数据的值以及图形是怎么样的呢？）
下面只需要将上面的代码注释掉就可以了。
#data1 = (data1 - data1.mean()) / data1.std() #print(data1.head()) #data1.head()

从上面的结果看出，梯度下降算法在没有进行特征归一化时， 计算出的$\theta$值是"nan nan nan"，经过查找相关信息，原来nan代表Not A Number（不是一个数），它并不等于0，因为nan不是一个数，所以相关计算都无法得到数字。当然，所有涉及nan的操作，返回的都是nan。在使用梯度下降算法时若不使用归一化，那么可能就计算不出最小的代价函数值。但是，在正规方程中，不经过归一化操作也可以进行损失函数$\theta$值的计算以及计算出最小的代价函数值。同时，控制台也会显示时间超时的警告，应该是在计算梯度下降算法时，因为数据过大，不能够计算出$\theta$值和代价函数最小值（即很难收敛或者说不能收敛）。通过结果还能看出，到了一定的迭代次数才能够计算出代价函数，因此该结果与下面归一化后的结果看起来是“相反的”。

进行归一化操作可以在最优解的寻优过程明显会变得平缓，加快了梯度下降求最优解的速度，更容易正确的收敛到最优解。


实验将迭代次数设置为iterations = 150000，alpha = 0.01时，算法的运行速度明显慢了很多。如图所示。

实验将迭代次数设置为iterations = 150000，alpha = 1000000时，学习率增大，迭代次数也相对的会减少，但是收敛速度相比归一化后的收敛速率慢了很多。

下面是上图放大后的图形，在第9次迭代后，能够计算出代价函数的数值了。

代价函数与迭代次数的关系图：
（$\alpha$ = 0.001）

（这是在学习率$\alpha$为0.01时的关系图）

当改变学习率时，结果如下：（$\alpha$ = 0.1）

（$\alpha$ = 0.5）

（$\alpha$ = 0.5）
（$\alpha$ = 1）

从上面的对比可以看出，在学习率不同的情况下，代价函数和迭代次数的关系呈现出不一样的关系。当学习率越小，所呈现的关系是稍微平滑一些的，随着学习率的增加，代价函数与迭代次数的关系曲线显示出“垂直”的关系，说明当学习率越大的时候，所需要的迭代次数会随之减少。

在梯度下降和正规方程算法中得出的$\theta$值以及根据$\theta$值计算出的代价函数取得最小的值：

如果您在阅读之中发现文章错误之处或者出现疑问，欢迎在评论指出。