详解标准方程法（内含公式推导和代码）

标准方程法

标准方程法就是通过一些方法直接求出权值ω的值，以往的文章已将讲过了代价函数的相关内容，这里直接给出代价函数的公式：
$$J({\omega }_0,{\omega }_1,{\omega }_2,\dots ,{\omega }_n)=\sum _{i=1}^m[y_i-h_{\omega }(x_i){]}^2$$
其中$h_{\omega }(x)$的公式如下：
$$h_{\omega }(x)={\omega }_0+{\omega }_1{x}_1+{\omega }_2{x}_2+\cdots +{\omega }_n{x}_n$$
如果数据集中只有一条数据，通过矩阵可以描述成一下形式：
$$X=\left[\begin{array}{ccccc}{x}_0& x_1& x_2& \dots & x_n\end{array}\right],其中x_0=1$$

$$w=\left[\begin{array}{c}{\omega }_0\\ {\omega }_1\\ {\omega }_2\\ ⋮\\ {\omega }_n\end{array}\right]$$

$$\begin{array}{r}{h}_{\omega }(x)=X\ast w=\left[\begin{array}{ccccc}{x}_0& x_1& x_2& \dots & x_n\end{array}\right]\ast \left[\begin{array}{c}{\omega }_0\\ {\omega }_1\\ {\omega }_2\\ ⋮\\ {\omega }_n\end{array}\right]\\ ={\omega }_0+{\omega }_1{x}_1+{\omega }_2{x}_2+\cdots +{\omega }_n{x}_n\end{array}$$

一般数据集中，会有很多条数据，每一条数据通过上述公式都会得到一个$h_{\omega }(x)$，故通过矩阵可以描述成以下形式：
$$\begin{gathered} x_0{x}_1{x}_2\dots x_n \\ X={\left[\begin{array}{ccccc}1& a_1& a_2& \dots & a_n\\ 1& b_1& b_2& \dots & b_n\\ \dots & \dots & \dots & \dots & \dots \\ 1& n_1& n_2& \dots & n_n\end{array}\right]}_{n\ast (n+1)} \end{gathered}$$

$$w={\left[\begin{array}{c}{\omega }_0\\ {\omega }_1\\ {\omega }_2\\ ⋮\\ {\omega }_n\end{array}\right]}_{(n+1)\ast 1}$$

$$\begin{array}{rl}{h}_{\omega }(x)=X\ast w& =\left[\begin{array}{ccccc}1& a_1& a_2& \dots & a_n\\ 1& b_1& b_2& \dots & b_n\\ \dots & \dots & \dots & \dots & \dots \\ 1& n_1& n_2& \dots & n_n\end{array}\right]\ast \left[\begin{array}{c}{\omega }_0\\ {\omega }_1\\ {\omega }_2\\ ⋮\\ {\omega }_n\end{array}\right]\\ & =\left[\begin{array}{c}{\omega }_0+{\omega }_1{a}_1+{\omega }_2{a}_2+\cdots +{\omega }_n{a}_n\\ {\omega }_0+{\omega }_1{b}_1+{\omega }_2{b}_2+\cdots +{\omega }_n{b}_n\\ ⋮\\ {\omega }_0+{\omega }_1{n}_1+{\omega }_2{n}_2+\cdots +{\omega }_n{n}_n\end{array}\right]\\ & ={\left[\begin{array}{c}{h}_{\omega }(x_1)\\ h_{\omega }(x_2)\\ ⋮\\ h_{\omega }(x_n)\end{array}\right]}_{n\ast 1}\end{array}$$

那么代价函数如何用这种形式表示呢？别急，看下面这个公式
$$\sum _{i=1}^n[y_i-h_{\omega }(x_i){]}^2=(y_i-Xw{)}^T(y_i-Xw)$$
证明如下
$$y_i-Xw=\begin{array}{r}\left[\begin{array}{c}{y}_1-h_{\omega }(x_1)\\ y_2-h_{\omega }(x_2)\\ ⋮\\ y_n-h_{\omega }(x_n)\end{array}\right]\end{array}$$

$$\begin{array}{rl}(y_i-Xw{)}^T(y_i-Xw)& =\left[\begin{array}{cccc}{y}_1-h_{\omega }(x_1)& y_2-h_{\omega }(x_2)& \dots & y_n-h_{\omega }(x_n)\end{array}\right]\ast \left[\begin{array}{c}{y}_1-h_{\omega }(x_1)\\ y_2-h_{\omega }(x_2)\\ ⋮\\ y_n-h_{\omega }(x_n)\end{array}\right]\\ & =(y_1-h_{\omega }(x_1){)}^2+(y_2-h_{\omega }(x_2){)}^2+(y_1-h_{\omega }(x_1){)}^2+\cdots +(y_n-h_{\omega }(x_n){)}^2\\ & =\sum _{i=1}^n[y_i-h_{\omega }(x_i){]}^2\end{array}$$

证明结束。

但是呢，到这里肯定是还没有结束的，我们知道了代价函数，接下来一步就是去求最小值。提到函数和最小值，用我为数不多的数学知识那就是求导了，我现在脑海中只剩了一顿求导，令其为零了。好了，接下来就让我们一起来求导吧！！！
$$\begin{array}{rl}& \frac{\partial (y-Xw{)}^T(y-Xw)}{\partial w}\\ & =\frac{\partial (y^T-w^T{X}^T)(y-Xw)}{\partial w}\\ & =\frac{\partial (y^{T}y-y^{T}Xw-w^T{X}^{T}y+w^T{X}^{T}Xw)}{\partial w}\\ & =\frac{\partial y^{T}y}{\partial w}-\frac{\partial y^{T}Xw}{\partial w}-\frac{\partial w^T{X}^{T}y}{\partial w}+\frac{\partial w^T{X}^{T}Xw}{\partial w}\end{array}$$
对其各自求偏导得：
$$\begin{array}{rl}\frac{\partial y^{T}y}{\partial w}& =0\\ \frac{\partial y^{T}Xw}{\partial w}& =X^{T}y\\ \frac{\partial w^T{X}^{T}y}{\partial w}& =\frac{(\partial w^T{X}^{T}y{)}^T}{\partial w}=\frac{\partial y^{T}Xw}{\partial w}=X^{T}y\\ \frac{\partial w^T{X}^{T}Xw}{\partial w}& =2X^{T}Xw\end{array}$$
整合一下：
$$\frac{\partial y^{T}y}{\partial w}-\frac{\partial y^{T}Xw}{\partial w}-\frac{\partial w^T{X}^{T}y}{\partial w}+\frac{\partial w^T{X}^{T}Xw}{\partial w}=0-X^{T}y-X^{T}y+2X^{T}Xw$$
令其为零：
$$\begin{array}{rl}-2X^{T}y+2X^{T}Xw& =0\\ X^{T}Xw& =X^{T}y\\ (X^{T}X{)}^{-1}{X}^{T}Xw& =(X^{T}X{)}^{-1}{X}^{T}y\\ w& =(X^{T}X{)}^{-1}{X}^{T}y\end{array}$$
到这里$w$就是我们苦苦寻求的值了，而其中$(X^{T}X{)}^{-1}$是$X^{T}X)$的逆矩阵。

不可逆的情况

这里再给出两种不可逆的情况

1. 线性相关特征（多重共线性）
   • x1为房屋面积，单位是平方英尺
   • x2为房屋面积，单位是平方米
2. 特征数据太多（样本数m$\le$特征数量n）

梯度下降和标准方程的比较

	梯度下降法	标准方程法
优点	当特征值非常多的时候也可以很好的工作	不需要学习率 不需要迭代 可以得到全局最优解
缺点	需要选择合适的学习率 需要迭代很多个周期 最能得到最优解的近似值	需要计算$(X^{T}X{)}^{-1}$ 时间复杂度大约是$O(n^3)$ n是特征数量

代码实现

import numpy as np
from numpy import genfromtxt
import matplotlib.pyplot as plt

# 载入数据
data = np.genfromtxt("data.csv", delimiter=",")
x_data = data[:,0,np.newaxis]
y_data = data[:,1,np.newaxis]
plt.scatter(x_data,y_data)
plt.show()

print(np.mat(x_data).shape)
print(np.mat(y_data).shape)
# 给样本添加偏置项
X_data = np.concatenate((np.ones((100,1)),x_data),axis=1)
print(X_data.shape)

# 标准方程法求解回归参数
def weights(xArr, yArr):
    xMat = np.mat(xArr)
    yMat = np.mat(yArr)
    xTx = xMat.T*xMat # 矩阵乘法
    # 计算矩阵的值,如果值为0，说明该矩阵没有逆矩阵
    if np.linalg.det(xTx) == 0.0:
        print("This matrix cannot do inverse")
        return
    # xTx.I为xTx的逆矩阵
    ws = xTx.I*xMat.T*yMat
    return ws
    
#使用numpy本身函数求解
#相关函数可以点击链接查看
#https://www.numpy.org.cn/article/basics/numpy_matrices_vectors.html#%E7%94%A8numpy%E6%B1%82%E8%A7%A3%E6%96%B9%E7%A8%8B%E7%BB%84
def numpy_weights(xArr, yArr):
    xMat = np.mat(xArr)
    yMat = np.mat(yArr)
    xT = np.transpose(xMat)#转置矩阵
    xTx = np.dot(xT, xMat)#矩阵乘
    if np.linalg.det(xTx) == 0.0:
        print("This matrix cannot do inverse")
        return
    xTy = np.dot(xT, yMat)
    ws = np.linalg.solve(xTx,xTy)
    return ws

#调用函数
ws = weights(X_data,y_data)
ws1 = numpy_weights(X_data,y_data)
print(ws)
print(ws1)
    
 # 画图
x_test = np.array([[20],[80]])
y_test = ws[0] + x_test*ws[1]
plt.plot(x_data, y_data, 'b.')
plt.plot(x_test, y_test, 'r')
plt.show()