引言

上一篇筆記中已經記錄了，如何對一個無解的線性方程組\(Ax=b\)求近似解。在這裏，我們先來回顧兩個知識點：

1. 如何判斷一個線性方程組無解：如果拿上面那個方程組\(Ax=b\)舉例，那就是向量\(b\)不在矩陣A對應的列空間中，至於列空間的概念，可以參考四個基本子空間那篇筆記
2. 如何對無解的方程組求近似解：根據上一篇筆記如何尋找一個投影矩陣可以有這麽一個思路，將向量\(b\)往矩陣\(A\)所在的列空間投影得到向量\(f\)，得到新的方程組\(A\hat{x}=f\)，這個\(\hat{x}\)便為近似解了。如果僅僅為了求近似解可以直接在\(Ax=b\)等式左右兩側同時左乘\(A^{\mathrm{T}}\)
   ，即\(A^{\mathrm{T}}Ax=A^{\mathrm{T}}b\)。這個和上面先求投影向量再求解是一樣的。

這篇筆記將會探究在機器學習的線性回歸如何求解損失函數。

\(Ax=b\)無解時求近似解

今天我們需要求一個線性方程組，長成這樣\[ \begin{equation} \left \{ \begin{array}{lr} 2 * w_1 + 2 * w_2 + b = 14 \ 4 * w_1 - 1 * w_2 + b = 5 \ 4 * w_1 + 0 * w_2 + b = 4 \ 4 * w_1 - 2 * w_2 + b = 3 \ 0 * w_1 - 3 * w_2 + b = -20 \end{array} \right. \end{equation} \]

線性回歸

拿預測房價舉例，談談什麽是最小二乘法。比如我們假設房價(price)與2個特征即面積(x1)、樓層(x2)有關。那麽我們的目標是找到一張三維空間中的平面去擬合一些數據(假設這些數據都經過歸一化處理)。先來看看平面怎麽定義的？\[ \begin{equation} price = w_1 * x_1 + w_2 * x_2 + b \end{equation} \]
我們希望所有的數據點都在這個平面上，那樣可以通過解線性方程組來算出這個平面的參數\(w_1,w_2,b\)，這正是線性代數中學到過的。

好的，我們現在有5筆數據(2, 2, 14)、(4, -1, 5)、(4, 0, 4)、(4 -2 3)、(0 -3 -20)，將它們代入(5)式得到我們的方程組吧，解出來\(w_1,w_2,b\)這樣一個線性模型就ok了。
可是這一步我們之前已經做過了，將這些數據代入方程組是無解的，即給出的這些數據根本不在一個平面上。那麽，現在我們放松條件，既然找不到一個平面能令所有的點都在它上面，我們找一個最優的平面總可以吧。

最優平面如何定義

假設我們已經有n組數據，每一組數據都是\((x_1,x_2,y)\)的集合。將一組數據\((x_1, x_2)\)代入(1)中求出price，我們認為每一組數據產生的誤差為\((price-y)^2\)，將每一組數據產生的誤差累加起來就是(6)式。即:\[ \begin{equation} J(w_1, w_2, b) = \sum_{i=1}^{n}(price_i-y_i)^2 \end{equation} \]
使得\(J(w_1, w_2, b)\)最小的那組參數\((w_1,w_2,b)\)，可以認為是最優平面的參數。

下面會給出一個動圖來展示最優平面是怎麽樣的一個情況（畫了好久才畫出滿意的效果，畫圖的代碼也會在末尾給出）：

可以看到圖中，紅色的點是我們實際的數據，這個藍色的透明平面是我畫出來的認為能擬合這些數據的最好平面。

如何能找到最優平面？

這仍然是一個數學問題，我們認為使得\(J(w_1, w_2, b)\)最小的那組參數\((w_1,w_2,b)\)，就是最終要尋找的最優平面的參數。

這樣的話，我們記\(J(w_1, w_2, b)\)為一個函數，求一個多元函數的最值我們在微積分中學到過就是求\(\frac{\partial J}{\partial w_1}, \frac{\partial J}{\partial w_2}, \frac{\partial J}{\partial b}\)，並且令它們都等於0，就能求出最終的解了。

這裏已經涉及到矩陣微積分的內容，我試著寫幾步：\[ \begin{equation} J(w_1, w_2, b) = (price-y)^{\mathrm{T}}(price-y) \end{equation} \]
\(price\)和\(y\)都是向量，再將\(price\)用參數\(w_1,w_2,b\)表示：\[ \begin{equation} J(w) = (Xw-y)^{\mathrm{T}}(Xw-y) \end{equation} \]
(8)式中，\(X\)的每一行是1組數據，它是一個nx3的矩陣；\(w\)是個向量\[ X=\left[ \begin{matrix} 第一筆數據的 \ x1 & x2 & 1 \第二筆數據的 \ x1 & x2 & 1 \. \. \. \第n筆數據的 \ x1 & x2 & 1 \\end{matrix} \right] \ \ \ \ \ \ \ \ w =\left[ \begin{matrix} w_1\w_2\b \end{matrix} \right] \]
繼續將(8)式化簡\[ \begin{equation} J(w) = (w^{\mathrm{T}}X^{\mathrm{T}}-y^{\mathrm{T}})(Xw-y) \end{equation} \]
接著去括號\[ \begin{equation} J(w) = w^{\mathrm{T}}X^{\mathrm{T}}Xw-y^{\mathrm{T}}Xw-w^{\mathrm{T}}X^{\mathrm{T}}y+y^{\mathrm{T}}y \end{equation} \]
其中，\(y^{\mathrm{T}}Xw\)與\(w^{\mathrm{T}}X^{\mathrm{T}}y\)是相等的，都是一個數，所以最終可以寫為\[ \begin{equation} J(w) = w^{\mathrm{T}}X^{\mathrm{T}}Xw-2w^{\mathrm{T}}X^{\mathrm{T}}y+y^{\mathrm{T}}y \end{equation} \]
下面就要進行矩陣微積分了，講實話我不會。但是我學會兩個trick能求出最終的\(w\)。

1. 第一個trick來自臺大的林軒田老師，我記得他很輕松地說可以把上面這個等式變換成我們會的一元二次等式，我當時帶著滿腹的懷疑按照他說的做了，不過真的得到了結果(驚嚇！可能這就是數學的魅力)。我們將(11)式變為\[ \begin{equation} J(x) = w^{\mathrm{T}}Aw - 2w^{\mathrm{T}}b + c \\ subject \ to \ A = X^{\mathrm{T}}X \\ \ \ \ \ \ \ \ \ \ \ \ b=X^{\mathrm{T}}y\\ \ \ \ \ \ \ \ \ \ \ \ c=y^{\mathrm{T}}y \end{equation} \]
   當然，這不是嚴格意義上的轉換，但是真的能讓我們像解熟悉的一元二次方程一樣求出解。對(12)求導令其為0，再將原來的值代入回去能得到\[ \begin{equation} 2X^{\mathrm{T}}Xw - 2X^{\mathrm{T}}y = 0 \end{equation} \]
   最終\[ \begin{equation} w = (X^{\mathrm{T}}X)^{-1}X^{\mathrm{T}}y \end{equation} \]

2. 第二種求解的辦法就是記住矩陣微積分的公式：

   y	\(\frac{\partial y}{\partial X}\)
   \(AX\)	\(A^{\mathrm{T}}\)
   \(X^{\mathrm{T}}A\)	\(A\)
   \(X^{\mathrm{T}}X\)	\(2X\)
   \(X^{\mathrm{T}}AX\)	\(AX+A^{\mathrm{T}}X\)

等等，(14)式好熟悉。這不就是求解線性方程組\(Ax=b\)這個方程組無解時的最優近似解麽。所以，機器學習的線性回歸其實就是最小二乘中的擬合問題。一開始就將這個問題看為求解線性方程組問題的話：\[ \left[ \begin{matrix} 第一筆數據的 \ x1 & x2 & 1 \第二筆數據的 \ x1 & x2 & 1 \. \. \. \第n筆數據的 \ x1 & x2 & 1 \\end{matrix} \right] \left[ \begin{matrix} w_1\w_2\b \end{matrix} \right]=\left[ \begin{matrix} 第一筆數據的 \ price \第二筆數據的 \ price \. \. \. \第n筆數據的 \ price \\end{matrix} \right] \]
不就是求這個方程組有沒有解麽？如果沒有解，我們就求近似解。這個近似解的求解方法就是上一篇筆記中一直強調的部分，在等式左右兩邊左乘矩陣的轉置，我們馬上能得到近似解。

畫圖代碼

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

x1 = np.linspace(-5, 5, 5)
x2 = x1
x1, x2 = np.meshgrid(x1, x2)
price = x1 * 3 + x2 * 4 - 5
np.random.seed(325)
data_x = np.random.randint(-5, 5, 5)
data_y = np.random.randint(-5, 5, 5)
data_z = data_x * 3 + data_y * 4 - 5
bias = np.array([5, 2, -3, 4, -3])
data_z = data_z + bias
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_wireframe(x1, x2, price, rstride=10, cstride=10)
for i in range(len(data_x)):
    ax.scatter(data_x[i], data_y[i], data_z[i], color='r')
ax.set_xlabel('x1')
ax.set_ylabel('x2')
ax.set_zlabel('price')
ax.set_xticks([-5, 0, 5])
ax.set_yticks([-5, 0,10])
ax.set_zticks([ -40, 0, 40])
plt.show()

線性回歸損失函數求解