跳转至

05.1 正规方程法

5.1 正规方程解法⚓︎

英文名是 Normal Equations。

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

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

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

5.1.1 简单的推导方法⚓︎

在做函数拟合(回归)时,我们假设函数 \(H\) 为:

\[H(w,b) = b + x_1 w_1+x_2 w_2+ \dots +x_n w_n \tag{2}\]

\(b=w_0\),则:

\[H(W) = w_0 + x_1 \cdot w_1 + x_2 \cdot w_2 + \dots + x_n \cdot w_n\tag{3}\]

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

\[H(W) = X \cdot W \tag{4}\]

公式5中的 \(X\)\(W\) 的矩阵形状如下:

\[ X = \begin{pmatrix} 1 & x_{1,1} & x_{1,2} & \dots & x_{1,n} \\\\ 1 & x_{2,1} & x_{2,2} & \dots & x_{2,n} \\\\ \vdots & \vdots & \vdots & \ddots & \vdots \\\\ 1 & x_{m,1} & x_{m,2} & \dots & x_{m,n} \end{pmatrix} \tag{5} \]
\[ W= \begin{pmatrix} w_0 \\\\ w_1 \\\\ \vdots \\\\ w_n \end{pmatrix} \tag{6} \]

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

\[H(W) = X \cdot W = Y \tag{7}\]

其中,Y的形状如下:

\[ Y= \begin{pmatrix} y_1 \\\\ y_2 \\\\ \vdots \\\\ y_m \end{pmatrix} \tag{8} \]

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

\[X^{\top} X W = X^{\top} Y \tag{9}\]

其中,\(X^{\top}\)\(X\) 的转置矩阵,\(X^{\top}X\) 一定是个方阵,并且假设其存在逆矩阵,把它移到等式右侧来:

\[W = (X^{\top} X)^{-1}{X^{\top} Y} \tag{10}\]

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

5.1.2 复杂的推导方法⚓︎

我们仍然使用均方差损失函数(略去了系数\(\frac{1}{2m}\)):

\[J(w,b) = \sum_{i=1}^m (z_i - y_i)^2 \tag{11}\]

\(b\) 看作是一个恒等于 \(1\) 的feature,并把 \(Z=XW\) 计算公式带入,并变成矩阵形式:

\[J(W) = \sum_{i=1}^m \left(\sum_{j=0}^nx_{ij}w_j -y_i\right)^2=(XW - Y)^{\top} \cdot (XW - Y) \tag{12}\]

\(W\) 求导,令导数为 \(0\),可得到 \(W\) 的最小值解:

\[ \begin{aligned} \frac{\partial J(W)}{\partial W} &= \frac{\partial}{\partial W}[(XW - Y)^{\top} \cdot (XW - Y)] \\\\ &=\frac{\partial}{\partial W}[(W^{\top}X^{\top} - Y^{\top}) \cdot (XW - Y)] \\\\ &=\frac{\partial}{\partial W}[(W^{\top}X^{\top}XW -W^{\top}X^{\top}Y - Y^{\top}XW + Y^{\top}Y)] \end{aligned} \tag{13} \]

求导后(请参考矩阵/向量求导公式):

第一项的结果是:\(2X^{\top}XW\)(分母布局,denominator layout)

第二项的结果是:\(X^{\top}Y\)(分母布局方式,denominator layout)

第三项的结果是:\(X^{\top}Y\)(分子布局方式,numerator layout,需要转置\(Y^{\top}X\)

第四项的结果是:\(0\)

再令导数为 \(0\)

\[ \frac{\partial J}{\partial W}=2X^{\top}XW - 2X^{\top}Y=0 \tag{14} $$ $$ X^{\top}XW = X^{\top}Y \tag{15} $$ $$ W=(X^{\top}X)^{-1}X^{\top}Y \tag{16} \]

结论和公式10一样。

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

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

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

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

5.1.3 代码实现⚓︎

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

\[ X = \begin{pmatrix} 1 & 10.06 & 60 \\\\ 1 & 15.47 & 74 \\\\ 1 & 18.66 & 46 \\\\ 1 & 5.20 & 77 \\\\ \vdots & \vdots & \vdots \\\\ \end{pmatrix} \tag{17} \]
\[ Y= \begin{pmatrix} 302.86 \\\\ 393.04 \\\\ 270.67 \\\\ 450.59 \\\\ \vdots \\\\ \end{pmatrix} \tag{18} \]

根据公式(10):

\[W = (X^{\top} X)^{-1}{X^{\top} Y} \tag{19}\]
  1. \(X\)\(1000\times 3\) 的矩阵,\(X\) 的转置是 \(3\times 1000\)\(X^{\top}X\) 生成 \(3\times 3\) 的矩阵;
  2. \((X^{\top}X)^{-1}\) 也是 \(3\times 3\) 的矩阵;
  3. 再乘以 \(X^{\top}\),即 \(3\times 3\) 的矩阵乘以 \(3\times 1000\) 的矩阵,变成 \(3\times 1000\) 的矩阵;
  4. 再乘以 \(Y\)\(Y\)\(1000\times 1\) 的矩阵,所以最后变成 \(3\times 1\) 的矩阵,成为 \(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