最小二乘法求导-公式推导

多元线性回归模型

1. 建立模型：模型函数

$\hat{Y}=X{W}^T$

如果有 n+1 条数据，每条数据有 m+1 种x因素（每种x因素都对应 1 个权重w），则
👉已知数据：实际Y值= [ y 0 y 1 y 2 y 3 . . . y n ]

​y0​y1​y2​y3​...yn​​ ​ ，X= [ x 00 , x 10 . . . x m 0 x 01 , x 11 . . . x m 1 x 02 , x 12 . . . x m 2 x 03 , x 13 . . . x m 3 . . . x 0 n , x 1 n . . . x m n ]

[x00,x10...xm0x01,x11...xm1x02,x12...xm2x03,x13...xm3...x0n,x1n...xmn]" role="presentation" style="position: relative;">⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢x00,x10...xm0x01,x11...xm1x02,x12...xm2x03,x13...xm3...x0n,x1n...xmn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{x}_{00},x_{10}...x_{m0}\\ x_{01},x_{11}...x_{m1}\\ x_{02},x_{12}...x_{m2}\\ x_{03},x_{13}...x_{m3}\\ ...\\ x_{0n},x_{1n}...x_{mn}\end{array}\right]$$

​x00​,x10​...xm0​x01​,x11​...xm1​x02​,x12​...xm2​x03​,x13​...xm3​...x0n​,x1n​...xmn​​ ​

👉未知数据：模型$\hat{Y}$值= [ y 0 ^ y 1 ^ y 2 ^ . . . y n ^ ]

[y0^y1^y2^...yn^]" role="presentation" style="position: relative;">⎡⎣⎢⎢⎢⎢⎢⎢y0^y1^y2^...yn^⎤⎦⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{c}\widehat{y_0}\\ \widehat{y_1}\\ \widehat{y_2}\\ ...\\ \widehat{y_n}\end{array}\right]$$

​y0​^​y1​^​y2​^​...yn​^​​ ​ 模型参数 W= [ w 0 , w 1 , w 2 , w 3 , . . . , w m ]

[w0,w1,w2,w3,...,wm]" role="presentation" style="position: relative;">[w0,w1,w2,w3,...,wm]$$\left[\begin{array}{c}{w}_0,w_1,w_2,w_3,...,w_m\end{array}\right]$$

[w0​,w1​,w2​,w3​,...,wm​​]

2. 学习模型：损失函数

2.1 损失函数-最小二乘法

Loss = $\sum ({\hat{y}}_{i计算}-y_{i实际}{)}^2$

$\widehat{Y_{计算}}$= [ y 0 ^ y 1 ^ y 2 ^ . . . y n ^ ]

​y0​^​y1​^​y2​^​...yn​^​​ ​， 实际Y值= [ y 0 y 1 y 2 . . . y n ]

[y0y1y2...yn]" role="presentation" style="position: relative;">⎡⎣⎢⎢⎢⎢⎢⎢y0y1y2...yn⎤⎦⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{c}{y}_0\\ y_1\\ y_2\\ ...\\ y_n\end{array}\right]$$

​y0​y1​y2​...yn​​ ​ ，$\widehat{Y_{计算}}-Y$= [ y 0 ^ − y 0 y 1 ^ − y 1 y 2 ^ − y 2 . . . y n ^ − y n ]

[y0^−y0y1^−y1y2^−y2...yn^−yn]" role="presentation" style="position: relative;">⎡⎣⎢⎢⎢⎢⎢⎢y0^−y0y1^−y1y2^−y2...yn^−yn⎤⎦⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{c}\widehat{y_0}-y_0\\ \widehat{y_1}-y_1\\ \widehat{y_2}-y_2\\ ...\\ \widehat{y_n}-y_n\end{array}\right]$$

​y0​^​−y0​y1​^​−y1​y2​^​−y2​...yn​^​−yn​​ ​
则Loss = [ y 0 ^ − y 0 , y 1 ^ − y 1 , y 2 ^ − y 2 , . . . , y n ^ − y n ] [ y 0 ^ − y 0 y 1 ^ − y 1 y 2 ^ − y 2 . . . y n ^ − y n ]

[y0^−y0,y1^−y1,y2^−y2,...,yn^−yn]" role="presentation" style="position: relative;">[y0^−y0,y1^−y1,y2^−y2,...,yn^−yn]$$\left[\begin{array}{c}\widehat{y_0}-y_0,\widehat{y_1}-y_1,\widehat{y_2}-y_2,...,\widehat{y_n}-y_n\end{array}\right]$$

[y0^−y0y1^−y1y2^−y2...yn^−yn]" role="presentation" style="position: relative;">⎡⎣⎢⎢⎢⎢⎢⎢y0^−y0y1^−y1y2^−y2...yn^−yn⎤⎦⎥⎥⎥⎥⎥⎥$$\left[\begin{array}{c}\widehat{y_0}-y_0\\ \widehat{y_1}-y_1\\ \widehat{y_2}-y_2\\ ...\\ \widehat{y_n}-y_n\end{array}\right]$$

[y0​^​−y0​,y1​^​−y1​,y2​^​−y2​,...,yn​^​−yn​​] ​y0​^​−y0​y1​^​−y1​y2​^​−y2​...yn​^​−yn​​ ​
Loss = $(\widehat{Y_{计算}}-Y{)}^T(\widehat{Y_{计算}}-Y)$

👉$\widehat{Y_{计算}}=X^{T}W$，因此 Loss = $(X^{T}W-Y{)}^T(X^{T}W-Y)$

$(X^{T}W-Y{)}^T=W^{T}X-Y^T=W^{T}X-Y^T$

则 Loss = $(X^{T}W-Y{)}^T(X^{T}W-Y)=W{X}^{T}X{W}^T-Y^{T}X{W}^T-W{X}^{T}Y+Y^{T}Y$

2.2 求导解析解

👉$\frac{\partial (Loss)}{\partial (W)}=\frac{\partial (W{X}^{T}X{W}^T)}{\partial (W)}-\frac{\partial (Y^{T}X{W}^T)}{\partial (W)}-\frac{\partial (W{X}^{T}Y)}{\partial (W)}+\frac{\partial (Y^{T}Y)}{\partial (W)}$
根据以下矩阵求导证明：

👉$\frac{\partial (Loss)}{\partial (W)}=\frac{\partial (W{X}^{T}X{W}^T)}{\partial (W)}-\frac{\partial (Y^{T}X{W}^T)}{\partial (W)}-\frac{\partial (W{X}^{T}Y)}{\partial (W)}+\frac{\partial (Y^{T}Y)}{\partial (W)}$

👉$\frac{\partial (Loss)}{\partial (W)}=2X^{T}XW-2X^{T}Y$

👉当$\frac{\partial (Loss)}{\partial (W)}=0，则W=\frac{1}{2}\ast (X^{T}X{)}^{-1}(2X^{T}Y)=(X^{T}X{)}^{-1}(X^{T}Y)$

当$(X^{T}X{)}^{-1}$计算时，只有当$X^{T}X$为满秩矩阵时，W才有解

当$W=(X^{T}X{)}^{-1}(X^{T}Y)$时，👉$\frac{\partial (Loss)}{\partial (W)}=0$，仅仅能证明Loss取到极值，并不能说明是极小值，还是极大值！(实际最小二乘法，本质就是个凸函数-平方和函数，也不用下列证明）

要如何判断Loss是极大值还是极小值？

当Loss处于极小值点时，一阶导$Los{s}^{{}^{\prime }}=\frac{d(Loss)}{W}=0$，二阶导$Los{s}^{{}^{\prime \prime }}>0$
当Loss处于极大值点时，一阶导$Los{s}^{{}^{\prime }}=\frac{d(Loss)}{W}=0$，二阶导$Los{s}^{{}^{\prime \prime }}<0$

已知最小二乘法损失函数一阶导 $$Los{s}^{{}^{\prime }}=\frac{d(Loss)}{W}=\frac{\partial (Loss)}{\partial (W)}=2X{X}^{T}W-2X{Y}^T$$
则二阶导为 L o s s ′ ′ = d ( 2 X X T W − 2 X Y T ) W = 2 X X T = 2 ∗ [ x 00 , x 10 . . . x m 0 x 01 , x 11 . . . x m 1 x 02 , x 12 . . . x m 2 x 03 , x 13 . . . x m 3 . . . x 0 n , x 1 n . . . x m n ] [ x 00 , x 01 . . . x 0 n x 10 , x 11 . . . x 1 n x 20 , x 21 . . . x 2 n x 30 , x 31 . . . x 3 n . . . x m 0 , x m 1 . . . x m n ] = [ x 00 2 , . . . , . . . , . . . , . . . . . . , x 11 2 . . . . , . . . , . . . . . . , . . . , x 33 2 , . . . , . . . . . . , . . . , . . . , . . . , x m n 2 ] Loss^{''}=\frac{d(2XX^TW-2XY^T)}{W}=2XX^T=2*

= Loss′′=Wd(2XXTW−2XYT)​=2XXT=2∗ ​x00​,x10​...xm0​x01​,x11​...xm1​x02​,x12​...xm2​x03​,x13​...xm3​...x0n​,x1n​...xmn​​ ​ ​x00​,x01​...x0n​x10​,x11​...x1n​x20​,x21​...x2n​x30​,x31​...x3n​...xm0​,xm1​...xmn​​ ​= ​x002​,...,...,...,......,x112​....,...,......,...,x332​,...,......,...,...,...,xmn2​​ ​

由于主元全为正数，且矩阵对称，因此二阶导数矩阵为正定实对称矩阵，特征值全大于0

马马虎虎…地…对于正定矩阵、实对称矩阵已经懵圈

2.3 迭代近似解-梯度下降法

参数迭代公式：$w_{k+1}=w_k-\eta \ast \frac{d(Los{s}_k)}{d{w}_k}$

沿着梯度的方向，正常情况下，Loss值会逐渐减小。
但当达到最小值后，下一次迭代，Loss值就会逐渐变大，因此当迭代后的Loss值比上次迭代的Loss值大，即$Los{s}_{k+1}>Los{s}_k$就停止迭代。并取上次迭代的Loss值（即$Los{s}_k$)作为最终损失函数值，以上次迭代的$W_k$参数，为模型参数。

但有时迭代次数过大，会导致训练很久很久很久，因此，可以设置一个最大迭代次数。

当迭代次数超过最大迭代次数后，仍未找到极值点$Los{s}_k$，则停止迭代。（可修改参数初始值或学习率后，重新训练）

因此，迭代停止的条件为两个：
$Los{s}_{k+1}>Los{s}_k$ 或 超过最大迭代次数

不过！！！$Los{s}_{k+1}>Los{s}_k$ 并不意味着取到极小值，因为又可能迭代步长过长，使Loss超过极值点，因此，一般可以限制一个Loss差值范围，但这个差值范围要设置为多少呢？？？
这又是一门学问了，我还不懂、、、、

另外，有可能出现反复震荡，无法收敛到极小值的现象，那怎么办。。。凉拌

简单点：迭代停止的条件为-超过最大迭代次数，或迭代后的$Los{s}_{k+1}>Los{s}_k$

3. 手动代码实现

# 1. 建立模型:多元线性回归模型
rows,columns = X.shape
X = np.array(X)
W = np.zeros(columns+1)
X_0 = np.ones((rows,1))
X = np.concatenate((X_0,X),axis=1)
Y_hat = np.matmul(X,W.T) # 创建多元线性回归模型
# print(Y.shape,Y_hat.shape)
# 2. 学习模型：损失函数模型+优化算法
# 求导法
def learn_model(X,Y,Y_hat):
    W = np.matmul((np.matmul(X.T,X))**-1,np.matmul(X.T,Y))
    Y_hat = np.matmul(X,W.T)
    Loss = np.sum((Y - Y_hat) *(Y - Y_hat))
    print("——————————求导法——————————")
    print(f'模型参数W为：{W}')
    print(f"最终损失值为{round(Loss, 2)}")

    return W,Loss,Y_hat
# 梯度下降法
def gradient_down(Y_hat,W):
    Loss1 = np.sum((Y - Y_hat) *(Y - Y_hat))
    Loss = np.sum((Y - Y_hat) *(Y - Y_hat))
    η = 0.01 # np.arange([0.1 for i in range(columns)]) # 设置学习率η
    times = 30
    # Y_hat = X*(W.T)
    while Loss1 <= Loss and times!=0:
        Loss = Loss1
        gradient = 2*np.matmul(np.matmul(X.T,X),W) - 2*np.matmul(X.T,Y)
        W = W - η * gradient
        Y_hat = np.matmul(X,W.T)
        Loss1 = np.sum((Y - Y_hat) *(Y - Y_hat))
        # print(f"第{31-times}次迭代，损失值为：{round(Loss,2)}")
        times -= 1
    print("——————————迭代法——————————")
    print(f"最终迭代次数为{30-times}，损失值为{round(Loss,2)}")
    print(f'模型参数W为：{W}')
learn_model(X,Y,Y_hat)
gradient_down(Y_hat,W)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

3.1 求导法的实践问题

我不理解！！！！！！！为什么求导法得出的LOSS值，居然大于梯度下降法迭代1次的Loss值
why???
what’s wrong with my 程序。。。。

虽然没有进行归一化处理。。。但也不至于这么个结果啊！
于是，进行一下归一化处理

# 归一化处理
datas_1 = (datas_1-datas_1.min())/(datas_1.max()-datas_1.min())
1
2

I don’t understand

实际进行多元线性回归时，是不需要进行归一化的。因为预测Y的实值 $\hat{Y}=X{W}^T$，对X归一化后，数值大小发生了改变，预测出的实值也不再是原本的数值大小。

就像是原本收益Y一般是1000的，归一化后预测就会变得非常小，例如0.9

猜想：迭代法的$W=(X^{T}X{)}^{-1}(X^{T}Y)$，要求$(X^{T}X{)}^{-1}$是满秩矩阵，但实际数据的$(X^{T}X{)}^{-1}$并不是满秩矩阵，但numpy也能求出非满秩矩阵的逆矩阵，猜测，可以是伪逆矩阵求解。。。。【不懂，我胡说八道…】

总之，求导法实际不太可行。

3.1 梯度下降法的实践问题

梯度下降法有个很大的问题，很大很大的问题，学习率η的选取问题！！！！big problem

唉，懒得复盘了

η小小的改变，迭代效果差距非常大！！！！！！！！

看网上说，η也应该随着切线斜率同步改变的，但是。。。。。懒得想了

还不如直接就多运行几次，慢慢调整一下η的值就好

对比一下手动训练模型和sklearn训练模型的速度

查了一下sklearn的官方文档，说是用随机梯度下降法，但具体的实现没有曝光
但是，sklearn的速度好快啊！！！！！我…是它的670倍。。。。。。。。

呵呵。。。调包侠岂不是更快活！！！

import pandas as pd
import numpy as np
from sklearn import linear_model
import time
import matplotlib.pyplot as plt

# 获取所需数据：'推荐分值', '专业度','回复速度','用户群活跃天数'
datas = pd.read_excel('./datas1.xlsx')
important_features = ['推荐分值', '辅导老师专业度','回复速度','群活跃天数']
datas_1 = datas[important_features]

# 明确实值Y为'推荐分值'，X分别为'专业度','回复速度','用户群活跃天数'
Y = datas_1['推荐分值']
X_original = datas_1.drop('推荐分值',axis=1)

# 1. 建立模型:多元线性回归模型
rows,columns = X_original.shape
X_original = np.array(X_original)
X_0 = np.ones((rows,1))
X = np.concatenate((X_0,X_original),axis=1)
# 计时器：这是计时装饰器函数，参数 func 是被装饰的函数
def cal_time(func):
    def wrapper(*args, **kw):
        start_time = time.time()
        func(*args, **kw)
        end_time = time.time()
        print(f'用时：{end_time-start_time}秒\n')
    return wrapper
# 2. 学习模型：损失函数模型+优化算法
# 求导法
@cal_time
def learn_model():
    W = np.matmul((np.matmul(X.T,X))**-1,np.matmul(X.T,Y))
    Y_hat = np.matmul(X,W.T)
    Loss = np.sum((Y - Y_hat) *(Y - Y_hat))
    print("——————————手动：求导法——————————")
    print(f'模型参数W为：{W}')
    print(f"最终损失值为{round(Loss, 2)}")
    return W,Loss,Y_hat


# 梯度下降法
@cal_time
def gradient_down():
    W0 = np.zeros(columns + 1)
    Y_hat = np.matmul(X, W0.T)  # 创建多元线性回归模型
    Loss0 = np.sum((Y - Y_hat) *(Y - Y_hat))
    a = 0.00002 # 设置学习率a
    times = 10000000
    num = 0
    while num <= times:
        gradient = 2*np.matmul(np.matmul(X.T,X),W0) - 2*np.matmul(X.T,Y)
        W1 = W0 - a * gradient
        Y_hat = np.matmul(X,W1.T)
        Loss1 = np.sum((Y - Y_hat) *(Y - Y_hat))
        if Loss1 > Loss0:
            break
        # print(f"第{num+1}次迭代，损失值为：{round(Loss0,2)}")
        num += 1
        Loss0 = Loss1
        W0 = W1
    print("——————————手动：梯度下降法——————————")
    print(f'模型参数W为：{W0}')
    print(f"最终迭代次数为{num+1}，损失值为{round(Loss0,2)}")

# sklearn的坐标下降法
@cal_time
def sklearn_mulreg():
    # 1. 建立模型
    reg = linear_model.LinearRegression()

    # 2. 学习模型
    reg.fit(X_original, Y)
    w = reg.coef_
    b = reg.intercept_


    print("——————————sklearn:随机梯度下降法——————————")
    print(f'模型参数W为：{w,b}')
    Y_hat = np.matmul(X_original, w)+b
    Loss = np.sum((Y - Y_hat) * (Y - Y_hat))
    print(f"加了正则化，最终损失值为{round(Loss,2)}")

learn_model()
gradient_down()
sklearn_mulreg()


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88