[ch05-01] 正规方程法解决多变量线性回归问题

Posted 2021-03-23 woodyh5

系列博客，原文在笔者所维护的github上：https://aka.ms/beginnerAI，
点击star加星不要吝啬，星越多笔者越努力。

5.1 正规方程解法

英文名是 Normal Equations。

对于线性回归问题，除了前面提到的最小二乘法可以解决一元线性回归的问题外，也可以解决多元线性回归问题。

对于多元线性回归，可以用正规方程来解决，也就是得到一个数学上的解析解。它可以解决下面这个公式描述的问题：

[y=a_0+a_1x_1+a_2x_2+dots+a_kx_k ag{1}]

5.1.1 简单的推导方法

在做函数拟合（回归）时，我们假设函数H为：

[h(w,b) = b + x_1 w_1+x_2 w_2+...+x_n w_n ag{2}]

令(b=w_0)，则：

[h(w) = w_0 + x_1 cdot w_1 + x_2 cdot w_2+...+ x_n cdot w_n ag{3}]

公式3中的x是一个样本的n个特征值，如果我们把m个样本一起计算，将会得到下面这个矩阵：

[H(w) = X cdot W ag{4}]

公式5中的X和W的矩阵形状如下：

[ X = egin{pmatrix} 1 & x_{1,1} & x_{1,2} & dots & x_{1,n} 1 & x_{2,1} & x_{2,2} & dots & x_{2,n} \dots 1 & x_{m,1} & x_{m,2} & dots & x_{m,n} end{pmatrix} ag{5} ]

[ W= egin{pmatrix} w_0 w_1 \dots w_n end{pmatrix} ag{6} ]

然后我们期望假设函数的输出与真实值一致，则有：

[H(w) = X cdot W = Y ag{7}]

其中，Y的形状如下：

[ Y= egin{pmatrix} y_1 y_2 \dots y_m end{pmatrix} ag{8} ]

直观上看，W = Y/X，但是这里三个值都是矩阵，而矩阵没有除法，所以需要得到X的逆矩阵，用Y乘以X的逆矩阵即可。但是又会遇到一个问题，只有方阵才有逆矩阵，而X不一定是方阵，所以要先把左侧变成方阵，就可能会有逆矩阵存在了。所以，先把等式两边同时乘以X的转置矩阵，以便得到X的方阵：

[X^T X W = X^T Y ag{9}]

其中，(X^T)是X的转置矩阵，(X^T X)一定是个方阵，并且假设其存在逆矩阵，把它移到等式右侧来：

[W = (X^T X)^{-1}{X^T Y} ag{10}]

至此可以求出W的正规方程。

5.1.2 复杂的推导方法

我们仍然使用均方差损失函数：

[J(w,b) = sum (z_i - y_i)^2 ag{11}]

把b看作是一个恒等于1的feature，并把z=XW计算公式带入，并变成矩阵形式：

[J(w) = sum (x_i w_i -y_i)^2=(XW - Y)^T cdot (XW - Y) ag{12}]

对w求导，令导数为0，就是W的最小值解：

[ egin{aligned} {partial J(w) over partial w} &= {partial over partial w}[(XW - Y)^T cdot (XW - Y)] &={partial over partial w}[(X^TW^T - Y^T) cdot (XW - Y)] &={partial over partial w}[(X^TXW^TW -X^TW^TY - Y^TXW + Y^TY)] end{aligned} ag{13} ]

求导后：

第一项的结果是：(2X^TXW)

第二项和第三项的结果都是：(X^TY)

第四项的结果是：0

再令导数为0：

[ J'(w)=2X^TXW - 2X^TY=0 ag{14} ]
[ X^TXW = X^TY ag{15} ]
[ W=(X^TX)^{-1}X^TY ag{16} ]

结论和公式10一样。

以上推导的基本公式可以参考第0章的公式60-69。

逆矩阵((X^TX)^{-1})可能不存在的原因是：

1. 特征值冗余，比如(x_2=x^2_1)，即正方形的边长与面积的关系，不能做为两个特征同时存在
2. 特征数量过多，比如特征数n比样本数m还要大

以上两点在我们这个具体的例子中都不存在。

5.1.3 代码实现

我们把表5-1的样本数据带入方程内。根据公式(5)，我们应该建立如下的X,Y矩阵：

[ X = egin{pmatrix} 1 & 10.06 & 60 1 & 15.47 & 74 1 & 18.66 & 46 1 & 5.20 & 77 \dots \end{pmatrix} ag{17} ]

[ Y= egin{pmatrix} 302.86 393.04 270.67 450.59 \dots \end{pmatrix} ag{18} ]

根据公式(10)：

[W = (X^T X)^{-1}{X^T Y} ag{10}]

1. X是1000x3的矩阵，X的转置是3x1000，(X^TX)生成(3x3)的矩阵
2. ((X^TX)^{-1})也是3x3
3. 再乘以(X^T)，即(3x3)x(3x1000)的矩阵，变成3x1000
4. 再乘以Y，Y是1000x1，所以(3x1000)x(1000x1)变成3x1，就是W的解，其中包括一个偏移值b和两个权重值w，3个值在一个向量里

if __name__ == '__main__':
    reader = SimpleDataReader()
    reader.ReadData()
    X,Y = reader.GetWholeTrainSamples()
    num_example = X.shape[0]
    one = np.ones((num_example,1))
    x = np.column_stack((one, (X[0:num_example,:])))
    a = np.dot(x.T, x)
    # need to convert to matrix, because np.linalg.inv only works on matrix instead of array
    b = np.asmatrix(a)
    c = np.linalg.inv(b)
    d = np.dot(c, x.T)
    e = np.dot(d, Y)
    #print(e)
    b=e[0,0]
    w1=e[1,0]
    w2=e[2,0]
    print("w1=", w1)
    print("w2=", w2)
    print("b=", b)
    # inference
    z = w1 * 15 + w2 * 93 + b
    print("z=",z)

5.1.4 运行结果

w1= -2.0184092853092226
w2= 5.055333475112755
b= 46.235258613837644
z= 486.1051325196855

我们得到了两个权重值和一个偏移值，然后得到房价预测值z=486万元。

至此，我们得到了解析解。我们可以用这个做为标准答案，去验证我们的神经网络的训练结果。

代码位置

ch05, Level1