如何了解正規方程

在《如何了解線性回歸》中，我們提到了使用梯度下降來求解最小二乘問題，但線上性回歸中，還有一種不需要疊代的方法來求解最小二乘問題，這就是正規方程。這是基于矩陣求導來及計算的

還是上次的資料

x = [[1, 1, 1, 1, 1], [1, 2, 3, 4, 5]]
y = [4.8, 7.2, 8.8, 11.2, 12.8]
plt.scatter(x[1], y)
plt.xlim(0,6)
plt.ylim(0,14)
plt.show()
           

和線性回歸前半部分的了解一樣，我們的目标是最小化 J J J，隻是這裡使用的 J J J是 J = 1 2 ∑ ( y − y h ) 2 J=\frac{1}{2}\sum(y-y_h)^2 J=21​∑(y−yh​)2（我們用 y y y表示真實值， y h y_h yh​表示預測值），标準的寫法是這樣的，之前我們在用梯度下降時 J J J沒有 1 / 2 1/2 1/2也沒有關系。那麼現在我們就要用矩陣的形式來表示 J J J，并取它的導數，當導數等于零時，也就可以得到最終的正規方程，下面是 J J J的矩陣表示

J = 1 2 ∑ ( y − y h ) 2 = 1 2 t r [ ( X θ − Y ) T ( X θ − Y ) ] \begin{aligned} J & = \frac{1}{2}∑(y - y_h)^2 \\ & = \frac{1}{2}tr[(Xθ - Y)^T (Xθ - Y)] \end{aligned} J​=21​∑(y−yh​)2=21​tr[(Xθ−Y)T(Xθ−Y)]​

t r tr tr是什麼， t r tr tr表示矩陣的迹，也就是矩陣主對角線元素的和，表示為 t r A = ∑ a i i trA=∑a_{ii} trA=∑aii​，其中有些重要的定理

1. t r ( A B ) = t r ( B A ) 2. t r ( A + B ) = t r ( A ) + t r ( B ) 3. d ( t r ( A B ) ) / d ( A ) = B T 4. d ( t r ( A T B ) ) / d ( A ) = B 5. t r ( A ) = t r ( A T ) 6. t r ( a ) = a （ a ∈ R ） 7. d ( t r ( A B A T C ) ) / d ( A ) = C A B + C T A B T \begin{aligned} &{1.}tr(AB) = tr(BA)\\ &{2.}tr(A+B) = tr(A) + tr(B)\\ &{3.}d(tr(AB))/d(A) = B^T\\ &{4.}d(tr(A^TB))/d(A) = B\\ &{5.}tr(A) = tr(A^T)\\ &{6.}tr(a) = a（a ∈ R）\\ &{7.}d(tr(ABA^TC))/d(A) = CAB+C^TAB^T\\ \end{aligned} ​1.tr(AB)=tr(BA)2.tr(A+B)=tr(A)+tr(B)3.d(tr(AB))/d(A)=BT4.d(tr(ATB))/d(A)=B5.tr(A)=tr(AT)6.tr(a)=a（a∈R）7.d(tr(ABATC))/d(A)=CAB+CTABT​

就這個問題，我們可以有如下的表示

(1) x i = ( 1 x ) 這 裡 的 x 就 是 上 面 的 [ 1 , 2 , 3 , 4 , 5 ] x_i= \left( \begin{matrix} 1\\ x \end{matrix} \right)\tag{1}這裡的x就是上面的[1, 2, 3, 4, 5] xi​=(1x​)這裡的x就是上面的[1,2,3,4,5](1) (2) θ = ( b a ) 這 是 y = a x + b 的 參 數 θ= \left( \begin{matrix} b\\ a \end{matrix} \right)\tag{2}這是y=ax+b的參數 θ=(ba​)這是y=ax+b的參數(2) (3) X = ( x 0 T x 1 T x 2 T x 3 T x 4 T ) X= \left( \begin{matrix} x_0^T\\ x_1^T\\ x_2^T\\ x_3^T\\ x_4^T \end{matrix} \right)\tag{3} X=⎝⎜⎜⎜⎜⎛​x0T​x1T​x2T​x3T​x4T​​⎠⎟⎟⎟⎟⎞​(3) (4) Y = ( y 0 y 1 y 2 y 3 y 4 ) 這 裡 的 y i 就 是 上 面 的 [ 4.8 , 7.2 , 8.8 , 11.2 , 12.8 ] Y= \left( \begin{matrix} y_0\\ y_1\\ y_2\\ y_3\\ y_4\\ \end{matrix} \right)\tag{4}這裡的y_i就是上面的[4.8, 7.2, 8.8, 11.2, 12.8] Y=⎝⎜⎜⎜⎜⎛​y0​y1​y2​y3​y4​​⎠⎟⎟⎟⎟⎞​這裡的yi​就是上面的[4.8,7.2,8.8,11.2,12.8](4)

上面的 x x x是一個 2 × 1 2\times1 2×1的向量，這樣 x T θ x^Tθ xTθ就是一個 1 × 2 1\times2 1×2的向量乘一個 2 × 1 2\times1 2×1，就得到一個答案，就是 a x + b ax+b ax+b，這個很眼熟吧。是以 X T θ − Y X^Tθ-Y XTθ−Y就包含了所有的 y − y h = a x + b − y h y-y_h=ax+b-y_h y−yh​=ax+b−yh​，轉置相乘再去迹就可以得到所有式子的平方和。這個式子比較好了解

麻煩的就是下面 d ( J ) / d ( θ ) d(J)/d(θ) d(J)/d(θ)。首先暴力拆開 ( X θ − Y ) T ( X θ − Y ) (Xθ-Y)^T(Xθ-Y) (Xθ−Y)T(Xθ−Y)

( X θ − Y ) T ( X θ − Y ) = ( θ T X T − Y T ) ( X θ − Y ) = θ T X T X θ − θ T X T Y − Y T X θ + Y T Y \begin{aligned} (Xθ - Y)^T (Xθ - Y) & = (θ^TX^T - Y^T)(Xθ - Y)\\ &= θ^TX^TXθ - θ^TX^TY - Y^TXθ + Y^TY \end{aligned} (Xθ−Y)T(Xθ−Y)​=(θTXT−YT)(Xθ−Y)=θTXTXθ−θTXTY−YTXθ+YTY​

是以這時要解 d ( J ) / d ( θ ) d(J)/d(θ) d(J)/d(θ)，我們算一下次元看， θ T X T X θ θ^TX^TXθ θTXTXθ是 ( 1 × 2 × 2 × 5 × 5 × 2 × 2 × 1 ) = ( 1 × 1 ) (1\times2 \times 2\times5 \times 5\times2 \times 2\times1)=(1\times1) (1×2×2×5×5×2×2×1)=(1×1)， θ T X T Y θ^TX^TY θTXTY是 ( 1 × 2 × 2 × 5 × 5 × 1 ) = ( 1 × 1 ) (1\times2 \times 2\times5 \times 5\times1)=(1\times1) (1×2×2×5×5×1)=(1×1)，還有 Y T X θ Y^TXθ YTXθ是 ( 1 × 5 × 5 × 2 × 2 × 1 ) = ( 1 × 1 ) (1\times5 \times 5\times2 \times 2\times1)=(1\times1) (1×5×5×2×2×1)=(1×1)， Y T Y Y^TY YTY是 ( 1 × 5 × 5 × 1 ) = ( 1 × 1 ) (1\times5 \times 5\times1)=(1\times1) (1×5×5×1)=(1×1)，都是些 1 × 1 1\times1 1×1的陣。苦逼沒學過矩陣的求導，隻有馬上去學了，至少可以看出最後的 Y T Y Y^TY YTY與 θ θ θ半毛錢關系都沒有，是以就可以忽略了，得到下面的式子

d ( J ) / d ( θ ) = 1 2 d ( t r ( θ T X T X θ − θ T X T Y − Y T X θ + Y T Y ) ) / d ( θ ) = 1 2 d ( t r ( θ T X T X θ − θ T X T Y − Y T X θ ) ) / d ( θ ) \begin{aligned} d(J)/d(θ) &= \frac{1}{2}d(tr(θ^TX^TXθ - θ^TX^TY - Y^TXθ + Y^TY)) / d(θ)\\ &=\frac{1}{2}d(tr(θ^TX^TXθ - θ^TX^TY - Y^TXθ)) / d(θ) \end{aligned} d(J)/d(θ)​=21​d(tr(θTXTXθ−θTXTY−YTXθ+YTY))/d(θ)=21​d(tr(θTXTXθ−θTXTY−YTXθ))/d(θ)​

然後可以發現 θ T X T Y θ^TX^TY θTXTY是 Y T X θ Y^TXθ YTXθ的轉置，也就是 ( Y T X θ ) T = θ T X T Y (Y^TXθ)^T=θ^TX^TY (YTXθ)T=θTXTY。由迹的一些定理（上面的2和5），我們可以這樣處理一下

d ( J ) / d ( θ ) = 1 2 d ( t r ( θ T X T X θ − θ T X T Y − Y T X θ + Y T Y ) ) / d ( θ ) = 1 2 d ( t r ( θ T X T X θ − θ T X T Y − Y T X θ ) ) / d ( θ ) = 1 2 ( d ( t r ( θ T X T X θ ) / d ( θ ) − d ( t r ( θ T X T Y ) / d ( θ ) − d ( t r ( Y T X θ ) / d ( θ ) ) = 1 2 ( d ( t r ( θ T X T X θ ) / d ( θ ) − 2 d ( t r ( θ T X T Y ) / d ( θ ) ) \begin{aligned} d(J)/d(θ) &= \frac{1}{2}d(tr(θ^TX^TXθ - θ^TX^TY - Y^TXθ + Y^TY)) / d(θ)\\ &=\frac{1}{2}d(tr(θ^TX^TXθ - θ^TX^TY - Y^TXθ)) / d(θ)\\ &=\frac{1}{2}(d(tr(θ^TX^TXθ) / d(θ) - d(tr(θ^TX^TY) / d(θ) - d(tr(Y^TXθ) / d(θ))\\ &=\frac{1}{2}(d(tr(θ^TX^TXθ) / d(θ) - 2d(tr(θ^TX^TY) / d(θ)) \end{aligned} d(J)/d(θ)​=21​d(tr(θTXTXθ−θTXTY−YTXθ+YTY))/d(θ)=21​d(tr(θTXTXθ−θTXTY−YTXθ))/d(θ)=21​(d(tr(θTXTXθ)/d(θ)−d(tr(θTXTY)/d(θ)−d(tr(YTXθ)/d(θ))=21​(d(tr(θTXTXθ)/d(θ)−2d(tr(θTXTY)/d(θ))​

然後現在需要求導了，這裡有一些有用的求導的定理

1. d ( a ) / d ( x ) = 0 2. d ( x ) / d ( x ) = I （ I 是 單 位 矩 陣 ） 3. d ( A x ) / d ( x ) = A T 4. d ( x T A ) / d ( x ) = A 5. d ( A T B A ) / d ( A ) = ( B + B T ) A \begin{aligned} &{1.}d(a)/d(x) = 0\\ &{2.}d(x)/d(x) = I （I是機關矩陣）\\ &{3.}d(Ax)/d(x) = A^T\\ &{4.} d(x^TA)/d(x) = A\\ &{5.}d(A^TBA)/d(A) = (B+B^T)A\\ \end{aligned} ​1.d(a)/d(x)=02.d(x)/d(x)=I（I是機關矩陣）3.d(Ax)/d(x)=AT4.d(xTA)/d(x)=A5.d(ATBA)/d(A)=(B+BT)A​

現在就可以繼續向下做了，由求導的5和迹的4可以得到如下的結果

d ( J ) / d ( θ ) = 1 2 d ( t r ( θ T X T X θ − θ T X T Y − Y T X θ + Y T Y ) ) / d ( θ ) = 1 2 d ( t r ( θ T X T X θ − θ T X T Y − Y T X θ ) ) / d ( θ ) = 1 2 ( d ( t r ( θ T X T X θ ) / d ( θ ) − d ( t r ( θ T X T Y ) / d ( θ ) − d ( t r ( Y T X θ ) / d ( θ ) ) = 1 2 ( d ( t r ( θ T X T X θ ) / d ( θ ) − 2 d ( t r ( θ T X T Y ) / d ( θ ) ) = 1 2 ( 2 X T X θ − 2 X T Y ) = X T X θ − X T Y \begin{aligned} d(J)/d(θ) &= \frac{1}{2}d(tr(θ^TX^TXθ - θ^TX^TY - Y^TXθ + Y^TY)) / d(θ)\\ &=\frac{1}{2}d(tr(θ^TX^TXθ - θ^TX^TY - Y^TXθ)) / d(θ)\\ &=\frac{1}{2}(d(tr(θ^TX^TXθ) / d(θ) - d(tr(θ^TX^TY) / d(θ) - d(tr(Y^TXθ) / d(θ))\\ &=\frac{1}{2}(d(tr(θ^TX^TXθ) / d(θ) - 2d(tr(θ^TX^TY) / d(θ))\\ &=\frac{1}{2}(2X^TXθ - 2X^TY)\\ &=X^TXθ - X^TY \end{aligned} d(J)/d(θ)​=21​d(tr(θTXTXθ−θTXTY−YTXθ+YTY))/d(θ)=21​d(tr(θTXTXθ−θTXTY−YTXθ))/d(θ)=21​(d(tr(θTXTXθ)/d(θ)−d(tr(θTXTY)/d(θ)−d(tr(YTXθ)/d(θ))=21​(d(tr(θTXTXθ)/d(θ)−2d(tr(θTXTY)/d(θ))=21​(2XTXθ−2XTY)=XTXθ−XTY​

這是讓 d ( J ) / d ( θ ) = 0 d(J)/d(θ)=0 d(J)/d(θ)=0，就是導數等于零，我們知道在 J J J最小時拟合最佳，而 J J J的導數為零時對應 J J J的一個極小值，這裡就是它的最小值，得到

d ( J ) / d ( θ ) = 0 X T X θ − X T Y = 0 X T X θ = X T Y θ = ( X T X ) − 1 X T Y \begin{aligned} d(J)/d(θ)& = 0\\ X^TXθ - X^TY& = 0\\ X^TXθ& = X^TY\\ θ& = (X^TX)^{-1}X^TY\\ \end{aligned} d(J)/d(θ)XTXθ−XTYXTXθθ​=0=0=XTY=(XTX)−1XTY​

這樣就得到了正規方程 θ = ( X T X ) − 1 X T Y θ=(X^TX)^{-1}X^TY θ=(XTX)−1XTY，這裡我們算算 θ θ θ的次元， ( X T X ) − 1 X T Y (X^TX)^{-1}X^TY (XTX)−1XTY得到 ( 2 × 5 × 5 × 2 × 2 × 5 × 5 × 1 ) = ( 2 × 1 ) (2\times5 \times 5\times2 \times 2\times5 \times 5\times1)=(2\times1) (2×5×5×2×2×5×5×1)=(2×1)，可以看到這裡的 θ θ θ就是一個 2 × 1 2\times1 2×1的向量，裡面包含了所有參數的值，這裡隻有 a a a和 b b b，是以這裡的 θ [ 0 ] = b θ[0]=b θ[0]=b，而 θ [ 1 ] = a θ[1]=a θ[1]=a

我們回到開頭的那個問題，看看資料和資料的次元

th = np.zeros((2, 1))  # 這是θ
X = np.array(x).T
Y = np.array(y).reshape(5, 1)
print(th.shape)
print(X.shape)
print(Y.shape)
           

這裡的輸出是

(2, 1)
(5, 2)
(5, 1)
           

現在我們構造這個正規方程并且計算，在python-numpy中，

np.dot()

是矩陣相乘，

np.linalg.inv()

是求逆矩陣，現在輸出看看

th = np.dot(np.dot(np.linalg.inv(np.dot(X.T, X)), X.T), Y)
print('a='+str(th[1])+' b='+str(th[0]))
           

結果為，也就是 y = 2 x + 2.96 y=2x+2.96 y=2x+2.96

a=[2.] b=[2.96]
           

我們畫個圖看看

plt.scatter(x[1], y)
x2 = np.linspace(0, 6, 100)
fy = th[1]*x2+th[0]
plt.plot(x2, fy)
plt.xlim(0,6)
plt.ylim(0,14)
plt.show()
           

這個資料是按照 y = 2 x + 3 y=2x+3 y=2x+3構造的，這個結果和梯度下降的結果是差不多的，這就是矩陣求導最小二乘的解法

梯度下降和正規方程都适用于求最小二乘的一種方法，相同點在于它們都求了導數，目标都是最小化J；不同在于梯度下降需要設定學習率和需要多次疊代，而正規方程不需要。一般使用時，n<10000時且為線性模型我們可以用正規方程，當n>10000時由于逆矩陣的運算複雜度高，效率很低；而且正規方程隻适用于線性模型，梯度下降法則可以适用各種模型，是更常用的做法