梯度下降法 (Gradient descent method)
香港中文大學數學系
2026年5月16日
1 單變量函數的梯度下降法 (GDM)
我們現在介紹梯度下降法，以數值方式求得最小二乘法問題的近似解。
尋找可微分的單變量函數 \(f(x)\) 最小值之問題如下。
import numpy as np
import matplotlib.pyplot as plt
# 建立 x 數值 (對應 xmin=-3, xmax=3)
x = np.linspace(-3, 3, 100)
# 定義函數 f(x) = x² + 1
y = x ** 2 + 1
# 建立圖形
plt.figure(figsize=(6, 4))
# 繪製函數 (藍色粗線)
plt.plot(x, y, 'blue', linewidth=2)
# 在 x* = 0 繪製垂直虛線 (從 y=0 到 y=1)
plt.plot([0, 0], [0, 1], 'k--') # 黑色虛線
# 在 x 軸下方標示 x*
plt.text(0, -0.3, r'$x^*$', ha='center', fontsize=12)
# 座標軸設定 (完全對應 TikZ)
plt.axhline(0, color='black', linewidth=0.8) # x 軸
plt.axvline(0, color='black', linewidth=0.8) # y 軸
plt.xlabel(r'$x$', fontsize=12)
plt.ylabel(r'$f(x)$', fontsize=12)
plt.xlim(-3, 3)## (-3.0, 3.0)## (-1.0, 5.0)## ([], [])## ([], [])
根據費馬臨界點定理，使函數達到最小值（稱為最佳解）的點 \(x^*\) 滿足以下條件： \[ f'(x^*) = 0. \] 因此，我們可以確定滿足方程式 \(f'(x) = 0\) 的哪個解是最佳解。然而，當函數十分複雜時，要透過解方程式來找到臨界點並不容易。在這種情況下，我們會以數值方式求得臨界點。在本節中，我們將探討梯度下降法，這是尋找給定函數最小值的最常用數值優化方法。梯度下降法的基本概念是找出函數的斜率，向「下坡」方向移動，並重複此過程直到達到極值。
梯度下降法 (Gradient Descent Method) 演算法
- 步驟 1 設定初始迭代值 \(x_1\)、容差 \(0 \leq \epsilon \ll 1\)、初始學習率 \(\eta\) (eta) 以及迭代次數 \(k := 1\)。
- 步驟 2 計算 \(g_k = f'(x_k)\)。如果 \(|g_k| \leq \epsilon\)，則停止。
- 步驟 3 設定 \(x_{k+1} = x_k - \eta g_k\)，\(k := k + 1\) 並回到 步驟 2。
讓我們解釋一下這個 GDM 演算法。設定 \(x_k\) 並計算 \(g_k = f'(x_k)\)。如果 \(|g_k| < \epsilon\)（在我們允許的誤差範圍內滿足臨界點定理），演算法就會停止並給出 \(x_k\) 作為最佳解。這裡，\(\epsilon\)（epsilon，容差）滿足 \(0 \leq \epsilon \ll 1\)。不等式中的 \(``\ll"\) 符號表示 epsilon 遠小於 1。因此，我們設定一個非常接近零的微小容差 \(\epsilon = 10^{-6}\)，以及一個學習率 \(\eta\)。同時將迭代次數 \(k\) 設為 1。如果 \(|g_k| > \epsilon\)，我們則使用公式 \(x_{k+1} = x_k - \eta g_k\) 來決定下一個值。接著再次計算 \(g_k = f'(x_k)\) 以判斷是否滿足臨界點定理。如果不滿足，則以類似方式繼續計算下一個點 \(x_1\) (意即進行下一輪更新)。透過這樣的方式，我們產生 \(x_1, x_2, \dots (\to x^*)\)。在經過多次重複步驟 \(x_{k+1} = x_k - \eta f'(x_k)\) 後，如果 \(g_k = f'(x_k)\) 的值落在給定的容差範圍內，演算法就會停止。我們預期某個 \(x_k\) 或極限 \(x^*\) 將會滿足 \(f'(x^*) = 0\)。
讓我們詳細看看 GDM 是如何運作的。假設經過第 \(k\) 次迭代後，函數 \(f\) 在近似解 \(x_k\) 處的切線斜率為負，即如下圖所示的 \(g_k = f'(x_k) < 0\)。這代表當 \(x_k\) 從左向右移動時，函數值會下降，所以我們可以預期 \(x^*\) 位於 \(x_k\) 的右側。如果 \(x_k\) 從左向右移動，那麼 \(x_{k+1}\) 就會更接近最佳值 \(x^*\)。因此，我們往 \(-\eta g_k > 0\) 的方向，從 \(x_k\) 移動至 \(x_{k+1} = x_k - \eta g_k\)。
import numpy as np
import matplotlib.pyplot as plt
# ----------------------
# 函數的數據
# ----------------------
x = np.linspace(-3, 3, 100)
y = x ** 2 + 1
# ----------------------
# 建立圖表
# ----------------------
plt.figure(figsize=(7, 5))
# 繪製 f(x) = x² + 1 (藍色粗線)
plt.plot(x, y, color='blue', linewidth=3)
# 在 x* = 0 繪製垂直虛線
plt.plot([0, 0], [0, 1], 'k--', linewidth=1.5)
# 紅色箭頭: (-1.5, 3.25) → (-0.75, 1.5625)
plt.arrow(
-1.5, 3.25, # 起點
0.75, -1.6875, # dx, dy
color='red',
linewidth=2,
head_width=0.12,
head_length=0.15,
length_includes_head=True
)
# ----------------------
# 所有標籤 (精確位置)
# ----------------------
plt.text(-1.7, 3.5, r"$f'(x_k) < 0$", color='red', fontsize=14)
plt.text(-1.5, -0.4, r"$x_k$", ha='center', fontsize=12)
plt.text(-0.75, -0.4, r"$x_{k+1}$", ha='center', fontsize=12)
plt.text(0, -0.4, r"$x^*$", ha='center', fontsize=12)
# ----------------------
# 座標軸設定 (對應 TikZ)
# ----------------------
plt.axhline(0, color='black', linewidth=1)
plt.axvline(0, color='black', linewidth=1)
plt.xlabel(r'$x$', fontsize=12)
plt.ylabel(r'$f(x)$', fontsize=12)
plt.xlim(-3, 3)## (-3.0, 3.0)## (-1.0, 5.0)## ([], [])## ([], [])
同樣地，假設經過第 \(k\) 次迭代後，函數 \(f\) 在近似解 \(x_k\) 處的切線斜率為正，即如下圖所示的 \(g_k = f'(x_k) > 0\)。這代表當 \(x_k\) 從左向右移動時，函數值是上升的，所以我們可以預期 \(x^*\) 位於 \(x_k\) 的左側。如果 \(x_k\) 向左移動，那麼 \(x_{k+1}\) 就會更接近最佳值 \(x^*\)。因此，我們往 \(-\eta g_k < 0\) 的方向，從 \(x_k\) 移動至 \(x_{k+1} = x_k - \eta g_k\)。如果我們重複這個過程直到 \(x_k\) 移動到 \(x^*\)，\(f'(x_k)\) 將會收斂至 0。透過這種方式，我們便能得到近似解 \(x_k\)。
import numpy as np
import matplotlib.pyplot as plt
# 建立曲線的 x 數值
x = np.linspace(-3, 3, 100)
y = x ** 2 + 1
# 建立圖形
plt.figure(figsize=(7, 5))
# 繪製函數: f(x) = x² + 1 (藍色粗線)
plt.plot(x, y, 'blue', linewidth=3)
# 在 x* = 0 繪製垂直虛線
plt.plot([0, 0], [0, 1], 'k--', linewidth=1.5)
# 紅色箭頭: (1.5, 3.25) → (0.75, 1.5625)
plt.arrow(1.5, 3.25, -0.75, -1.6875,
color='red', linewidth=2,
head_width=0.12, head_length=0.15,
length_includes_head=True)
# 所有標籤 (精確位置)
plt.text(1.7, 3.5, r"$f'(x_k) > 0$", color='red', fontsize=14)
plt.text(1.5, -0.4, r"$x_k$", ha='center', fontsize=12)
plt.text(0.75, -0.4, r"$x_{k+1}$", ha='center', fontsize=12)
plt.text(0, -0.4, r"$x^*$", ha='center', fontsize=12)
# 座標軸樣式 (完全對應 TikZ)
plt.axhline(0, color='black', linewidth=1)
plt.axvline(0, color='black', linewidth=1)
plt.xlabel(r'$x$', fontsize=12)
plt.ylabel(r'$f(x)$', fontsize=12)
# 座標軸極限
plt.xlim(-3, 3)## (-3.0, 3.0)## (-1.0, 5.0)## ([], [])## ([], [])
這說明了 GDM 將會找出滿足 \[f(x_1) > f(x_2) > f(x_3) > \dots\] 的 \(x_1, x_2, x_3, \dots\)。在這裡，\(\eta\) 決定了移動的幅度。我們稱之為「學習率」。如果學習率太大，\(x_{k+1}\) 可能會越過 \(x^*\)，甚至可能導致函數值不減反增。因此，我們必須謹慎選擇一個合理的學習率 \(\eta\)。
import numpy as np
import matplotlib.pyplot as plt
# 建立有兩個子圖的圖形 (並排)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# --------------------------
# 共用的函數數據
# --------------------------
x = np.linspace(-4, 4, 100)
y = x**2 + 1
# ==========================
# 左圖: 不收斂 (步長太大)
# ==========================
ax1.plot(x, y, 'blue', linewidth=3)
# 紅色水平箭頭: (-2,5) → (2,5)
ax1.arrow(-2, 5, 4, 0,
color='red', linewidth=2,
head_width=0.2, head_length=0.15, length_includes_head=True)
# 標籤
ax1.text(-2, -0.8, r'$x_k$', ha='center', fontsize=12)
ax1.text(2, -0.8, r'$x_{k+1}$', ha='center', fontsize=12)
ax1.text(0, -1.5, 'Do not converge', ha='center', fontsize=12)
# 座標軸樣式 (對應 TikZ)
ax1.axhline(0, color='black', linewidth=1)
ax1.axvline(0, color='black', linewidth=1)
ax1.set_xlabel(r'$x$')
ax1.set_ylabel(r'$f(x)$')
ax1.set_xlim(-4, 4)## (-4.0, 4.0)## (-1.0, 10.0)## []## []# ==========================
# 右圖: 太慢 (步長太小)
# ==========================
ax2.plot(x, y, 'blue', linewidth=3)
# 三個紅色小箭頭
ax2.arrow(-2, 5, 0.25, -0.9375,
color='red', linewidth=2, head_width=0.15, head_length=0.1, length_includes_head=True)
ax2.arrow(-1.75, 4.0625, 0.1875, -0.6211,
color='red', linewidth=2, head_width=0.15, head_length=0.1, length_includes_head=True)
ax2.arrow(-1.5625, 3.4414, 0.15625, -0.4687,
color='red', linewidth=2, head_width=0.15, head_length=0.1, length_includes_head=True)
# 標籤
ax2.text(-2, -0.8, r'$x_k$', ha='center', fontsize=12)
ax2.text(-1.75, -0.8, r'$x_{k+1}$', ha='center', fontsize=12)
ax2.text(-1.5625, -0.8, r'$x_{k+2}$', ha='center', fontsize=12)
ax2.text(-1.40625, -0.8, r'$x_{k+3}$', ha='center', fontsize=12)
ax2.text(-1, -1.5, 'Too slow', ha='center', fontsize=12)
# 座標軸樣式
ax2.axhline(0, color='black', linewidth=1)
ax2.axvline(0, color='black', linewidth=1)
ax2.set_xlabel(r'$x$')
ax2.set_ylabel(r'$f(x)$')
ax2.set_xlim(-4, 4)## (-4.0, 4.0)## (-1.0, 10.0)## []## []
import numpy as np
import matplotlib.pyplot as plt
# ==============================================
# 1. 定義函數及其導數 (來自 LaTeX)
# f(x) = x² + 1 f'(x) = 2x
# ==============================================
def f(x):
return x**2 + 1
def df(x):
return 2 * x
# ==============================================
# 2. 梯度下降演算法 (與 LaTeX 完全相符)
# ==============================================
def gradient_descent(x_start, eta=0.1, tol=1e-6, max_iter=1000):
x = x_start
path = [x] # 儲存迭代路徑以便繪圖
for k in range(max_iter):
grad = df(x)
# 收斂條件: |f'(x_k)| ≤ ε
if abs(grad) <= tol:
print(f"演算法成功！ (於第 {k+1} 次迭代收斂)")
break
# 更新規則: x_{k+1} = x_k - η * f'(x_k)
x = x - eta * grad
path.append(x)
else:
print("未能在最大迭代次數內收斂")
return x, np.array(path)
# ==============================================
# 3. 執行 GDM (使用 LaTeX 中的參數)
# ==============================================
x_opt, path = gradient_descent(x_start=1.5, eta=0.1, tol=1e-6)## 演算法成功！ (於第 68 次迭代收斂)##
## 最佳解 x* = 0.000000## 最小值 f(x*) = 1.000000## 迭代次數: 68# ==============================================
# 4. 繪圖 1: 函數 + 收斂路徑 (對應 LaTeX 圖形)
# ==============================================
x_plot = np.linspace(-3, 3, 500)
y_plot = f(x_plot)
plt.figure(figsize=(10, 5))
plt.plot(x_plot, y_plot, 'b-', linewidth=2, label=r'$f(x) = x^2 + 1$')
plt.scatter(path, f(path), color='red', s=30, zorder=5)
plt.plot(path, f(path), 'r--', linewidth=1.5, label='GD Path')
plt.axvline(x=0, color='k', linestyle='dashed', linewidth=1, label=r'$x^* = 0$')
plt.xlabel(r'$x$')
plt.ylabel(r'$f(x)$')
plt.title('Gradient Descent on $f(x) = x^2 + 1$')
plt.grid(alpha=0.3)
plt.legend()
plt.xlim(-3, 3)## (-3.0, 3.0)# ==============================================
# 5. 示範: 學習率太大 (發散)
# ==============================================
print("\n--- 測試過大的學習率 (η=1.0) ---")##
## --- 測試過大的學習率 (η=1.0) ---## 未能在最大迭代次數內收斂# ==============================================
# 6. 示範: 學習率太小 (太慢)
# ==============================================
print("\n--- 測試過小的學習率 (η=0.01) ---")##
## --- 測試過小的學習率 (η=0.01) ---## 未能在最大迭代次數內收斂
2 範例 1：單變量函數的梯度下降法
求 \(f(x) = 2x^2 - 3x + 2\) 的最小值。在下列程式碼中，取 \(x_1 = 0\)、\(\eta = 0.1\) 及 \(\epsilon = 10^{-6}\)。
解答
由 \(f'(x) = 4x - 3 = 0\)，臨界點為 \(x = \frac{3}{4} = 0.75\)。由於 \(f(x)\) 為凸函數（開口向上），因此在 \(x = \frac{3}{4}\) 處取得最小值 \[ f\left(\frac{3}{4}\right) = \frac{7}{8} = 0.875 \]
演算法的輸出為：
- 演算法成功！ (Algorithm Succeed!)
- \(x^* = 0.749999834194560\)
- \(|g^*| = 6.63221759289456e-7\)
- \(f(x^*) = 0.875000000000055\)
- 迭代次數 (Iteration number) = \(31\)
在上述範例中，梯度下降法所產生的點 \((x_1, f(x_1)),\ (x_2, f(x_2)),\ (x_3, f(x_3)),\ \dots\) 與函數 \(y = f(x)\) 的圖形一起顯示於座標平面上。從圖中可以輕易看出，該數列收斂到了曲線上具有最小值的點 \((x^*, f(x^*))\)。
可以看出該數列收斂到了最佳解 \(x^*\)。
\[ \begin{array}{cc} \textbf{參數 (Parameter)} & \textbf{數值 (Value)} \\ \hline \textrm{函數} & f(x) = 2x^2 - 3x + 2 \\ \textrm{導數} & f'(x) = 4x - 3 \\ \textrm{初始迭代值} ~x_1 & 0.0 \\ \textrm{容差}~ \epsilon & 10^{-6} \\ \textrm{學習率} ~\eta & 0.1 \\ \textrm{收斂解} ~ x^* & 0.75 \\ \textrm{最小值}~ f(x^*) & 0.875 \\ \textrm{收斂所需迭代次數} & 31 \\ \hline \end{array} \]
範例 1 梯度下降法結果總結
import numpy as np
import matplotlib.pyplot as plt
# ==============================================
# 1. 定義函數及其導數 (來自 LaTeX)
# f(x) = 2x² - 3x + 2
# f'(x) = 4x - 3
# ==============================================
def f(x):
return 2 * x**2 - 3 * x + 2
def df(x):
return 4 * x - 3
# ==============================================
# 2. 梯度下降法 (參數完全與 LaTeX 相符)
# ==============================================
x0 = 0.0 # 初始迭代值
tol = 1e-6 # 容差
eta = 0.1 # 學習率
max_iter = 300
path = [] # 儲存迭代路徑
for k in range(max_iter):
path.append(x0)
g0 = df(x0)
# 停止條件
if abs(g0) <= tol:
print("Algorithm Succeed!")
break
# 更新步驟
x0 = x0 - eta * g0## Algorithm Succeed!# ==============================================
# 3. 列印結果 (對應 LaTeX 輸出)
# ==============================================
print(f"x* = {x0}")## x* = 0.7499998341945602## |g*| = 6.632217592894563e-07## f(x*) = 0.8750000000000553## iteration number = 31# ==============================================
# 4. 繪圖 (完全對應 LaTeX 圖形)
# ==============================================
x_plot = np.linspace(-0.2, 1.2, 200)
y_plot = f(x_plot)
plt.figure(figsize=(8, 6))
# 繪製藍色曲線
plt.plot(x_plot, y_plot, 'b-', linewidth=2, label=r'$f(x)=2x^2-3x+2$')
# 繪製紅色迭代點
x_path = np.array(path)
y_path = f(x_path)
plt.scatter(x_path, y_path, color='red', s=25, zorder=5)
# 最佳點 x* = 0.75 處的虛線
plt.plot([0.75, 0.75], [0, 0.875], 'k--', linewidth=1)
plt.text(0.75, -0.05, r'$x^*$', ha='center', fontsize=12)
# 座標軸設定 (與 LaTeX 相同)
plt.axhline(0, color='black', linewidth=0.8)
plt.axvline(0, color='black', linewidth=0.8)
plt.xlabel(r'$x$')
plt.ylabel(r'$f(x)$')
plt.xlim(-0.2, 1.2)## (-0.2, 1.2)## (0.8, 2.2)plt.grid(alpha=0.3)
plt.title("Gradient Descent Path on $f(x) = 2x^2 - 3x + 2$")
plt.tight_layout()
plt.show()
3 範例 2：單變量函數的梯度下降法
求 \(f(x) = 9x^2 - 7x + 6\) 的最小值。利用上述的 GDM 程式碼，並設定 \(x_1 = 0\)、\(\eta = 0.1\) 及 \(\epsilon = 10^{-6}\)。
解答
首先，計算導數： \[ f'(x) = 18x - 7. \] 設定 \(f'(x) = 0\) 得到臨界點： \[ 18x - 7 = 0 \implies x^* = \frac{7}{18} \approx 0.3889. \] 由於 \(f''(x) = 18 > 0\)，此函數為凸函數，所以 \(x^* = \frac{7}{18}\) 即為極小值點。最小值為： \[ f\left(\frac{7}{18}\right) = 9\left(\frac{7}{18}\right)^2 - 7\left(\frac{7}{18}\right) + 6 = \frac{167}{36} \approx 4.6389. \]
範例 2 梯度下降法結果總結
\[ \begin{array}{cc} \textbf{參數 (Parameter)} & \textbf{數值 (Value)} \\ \hline \textrm{函數} & f(x) = 9x^2 - 7x + 6 \\ \textrm{導數} & f'(x) = 18x - 7 \\ \textrm{初始迭代值}~ x_1 & 0.0 \\ \textrm{容差}~ \epsilon & 10^{-6} \\ \textrm{學習率}~ \eta & 0.1 \\ \textrm{真實最小值}~ x^* & 7/18 \approx 0.3889 \\ \textrm{最小值}~ f(x^*) & 167/36 \approx 4.6389 \\ \hline \end{array} \]
import numpy as np
import matplotlib.pyplot as plt
# --------------------------
# 定義函數及其導數
# --------------------------
def f(x):
return 9 * x**2 - 7 * x + 6
def df(x):
return 18 * x - 7
# --------------------------
# 梯度下降法參數
# --------------------------
x = 0.0 # 初始值 x1 = 0
eta = 0.1 # 學習率 η = 0.1
epsilon = 1e-6 # 容差 ε = 10⁻⁶
max_iter = 500
path = [] # 儲存 x 值的路徑
# --------------------------
# 執行梯度下降法
# --------------------------
for k in range(max_iter):
path.append(x)
gradient = df(x)
if abs(gradient) <= epsilon:
print("Algorithm Succeed!")
break
x = x - eta * gradient## Algorithm Succeed!## Optimal x* = 0.38888894## Gradient at x* = 9.21e-07## Minimum f(x*) = 4.63888889## Iterations taken: 72# --------------------------
# 繪製函數與路徑圖
# --------------------------
x_plot = np.linspace(-0.2, 1.0, 200)
y_plot = f(x_plot)
plt.figure(figsize=(8, 5))
plt.plot(x_plot, y_plot, 'b-', linewidth=2, label=r'$f(x)=9x^2-7x+6$')
plt.scatter(path, [f(p) for p in path], color='red', s=25, zorder=5)
plt.axvline(x=7/18, color='k', linestyle='--', label=r'True minimum $x^*=7/18$')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Gradient Descent on $f(x)=9x^2-7x+6$')
plt.legend()
plt.grid(alpha=0.3)
plt.show()
4 梯度下降法的應用
第一節介紹的 GDM 是用於最小化單變量函數 \(f(x)\)。讓我們將此演算法推廣至最小化多變量函數。我們所學的最小二乘法問題是一個至少包含兩個自變數的多變量函數 \(f(\mathbf{x})\)。最小二乘法問題可以透過使用多變量轉換後的誤差函數 \(E(\mathbf{u})\) 來解決。這類問題同樣也能夠使用梯度下降法 (GDM) 來求解。
用於最小化多變量函數的 GDM 演算法
步驟 1 設定初始迭代值 \(\mathbf{u}_1 \in \mathbb{R}^n\)、容差 \(0\leq\epsilon<1\)、初始學習率 \(\eta\) (eta) 以及迭代次數 \(k:=1\)。
步驟 2 計算 \(\mathbf{g}_k = \nabla E(\mathbf{u}_k)\)。如果 \(\|\mathbf{g}_k\| \leq \epsilon\)，則停止。
步驟 3 設定 \(\mathbf{u}_{k+1} = \mathbf{u}_k - \eta \mathbf{g}_k\)，\(k:=k+1\) 並回到 步驟 2。
用於最小化多變量函數的 GDM 演算法的所有步驟與單變量函數的步驟完全相同。
單變量函數與多變量函數的對應關係
\[ \begin{array} \hline \textbf{單變量函數} & \textbf{多變量函數} \\ \hline \textrm{純量}~ x & \textrm{向量}~ \mathbf{u} \\ \textrm{絕對值}~ |a| & \textrm{向量範數}~ \|\mathbf{g}\| \\ \textrm{導數}~ f'(x) & \textrm{梯度}~ \nabla E(\mathbf{u}) \\ \hline \end{array} \]
令 \(x,y\) 為自變數，\(z\) 為第三個變數。現在我們可以考慮一個依賴於 \(x\) 和 \(y\) 的函數 \(z=f(x,y)\)。多變量函數也可以用同樣的方式定義。
在三維（座標）空間中，\(z=f(x,y)\) 可以被視為一個點 \((x,y,z)\)。隨著 \(x\) 和 \(y\) 的變動，\(z=f(x,y)\) 的圖形會形成一個曲面。舉例來說，雙變量函數 \(z=f(x,y)=-xye^{-x^2-y^2}\) 的圖形可以使用以下程式碼來繪製。
# 匯入所需函式庫
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# --------------------------
# 範例 1: f(x,y) = -xye^(-x²-y²) 的 3D 曲面圖
# --------------------------
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
# 建立網格
x_surf = np.linspace(-2, 2, 100)
y_surf = np.linspace(-2, 2, 100)
X_surf, Y_surf = np.meshgrid(x_surf, y_surf)
# 定義曲面函數
f = -X_surf * Y_surf * np.exp(-X_surf**2 - Y_surf**2)
# 繪製 3D 曲面
surf = ax.plot_surface(X_surf, Y_surf, f, cmap='viridis',
alpha=0.6, edgecolor='none')
ax.set_title('3D Surface of $f(x,y) = -xye^{-x^2 - y^2}$', fontsize=14)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.view_init(30, 30) # Same viewing angle as your LaTeX figure
# Show all plots
plt.tight_layout()
plt.show()
上圖中的 \(f\) 函數圖形直觀地顯示了 \(f\) 在山峰處有一個局部極大值，在山谷處有一個局部極小值。為了找出這些臨界點，我們需要引入多變量函數梯度的概念。
雙變量函數 \(f(x,y)\) 的梯度定義如下： \[ \operatorname{grad} f(x,y) = \nabla f(x,y) = \begin{bmatrix} \frac{\partial f}{\partial x} \\[6pt] \frac{\partial f}{\partial y} \end{bmatrix}. \]
\(f(x,y)\) 的梯度是一個 \(2\times1\) 向量，其中 \(\frac{\partial f}{\partial x}\) 代表 \(f\) 對 \(x\) 的偏導數。對 \(x\) 的偏導數意味著將 \(x\) 以外的其他變數視為常數。同樣地，\(\frac{\partial f}{\partial y}\) 代表 \(f\) 對 \(y\) 的偏導數。
找出 \(z = f(x,y) = -xye^{-x^2-y^2}\) 的梯度過程如下。
輸出為： \[ (2x^2 y e^{(-x^2 - y^2)} - y e^{(-x^2 - y^2)},\ 2x y^2 e^{(-x^2 - y^2)} - x e^{(-x^2 - y^2)}) \]
答案：
\[ \operatorname{grad} f(x,y) = \nabla f(x,y) = \left( 2x^2 y e^{-(x^2+y^2)} - y e^{-(x^2+y^2)},\ 2x y^2 e^{-(x^2+y^2)} - x e^{-(x^2+y^2)} \right) \]
import sympy as sp
import numpy as np
import matplotlib.pyplot as plt
# 1. 符號梯度計算 (SymPy)
x, y = sp.symbols('x y')
f = -x * y * sp.exp(-x**2 - y**2)
# 計算偏導數
df_dx = sp.diff(f, x)
df_dy = sp.diff(f, y)
# 列印符號梯度 (與圖片上的公式相符)
print("Symbolic Gradient (SymPy):")## Symbolic Gradient (SymPy):## ∂f/∂x = 2*x**2*y*exp(-x**2 - y**2) - y*exp(-x**2 - y**2)## ∂f/∂y = 2*x*y**2*exp(-x**2 - y**2) - x*exp(-x**2 - y**2)# 2. 轉換為數值函數 (NumPy)
df_dx_num = sp.lambdify((x, y), df_dx, 'numpy')
df_dy_num = sp.lambdify((x, y), df_dy, 'numpy')
f_num = sp.lambdify((x, y), f, 'numpy')
# 3. 視覺化梯度向量場
x_vals = np.linspace(-2, 2, 20)
y_vals = np.linspace(-2, 2, 20)
X, Y = np.meshgrid(x_vals, y_vals)
# 評估網格點的梯度
U = df_dx_num(X, Y)
V = df_dy_num(X, Y)
# 繪製向量場與曲面等高線
plt.figure(figsize=(8, 6))
contour = plt.contour(X, Y, f_num(X, Y), cmap='viridis')
plt.clabel(contour, inline=True, fontsize=8)
plt.quiver(X, Y, U, V, color='red', scale=50, width=0.003)
plt.title(r'Gradient Vector Field of $f(x,y) = -x y e^{-x^2 - y^2}$')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha=0.3)
plt.axis('equal')## (np.float64(-2.0), np.float64(2.0), np.float64(-2.0), np.float64(2.0))
5 多變量函數的梯度
同理，我們也能求出具有三個或更多變數之多變量函數 \(f\) 的梯度。一般而言，\(n\) 變量函數 \(f(x_1, x_2, \dots, x_n)\) 的梯度定義為： \[ \operatorname{grad} f(x_1, x_2, \dots, x_n) = \nabla f(x_1, x_2, \dots, x_n) = \begin{bmatrix} \frac{\partial f}{\partial x_1} \\[6pt] \frac{\partial f}{\partial x_2} \\[6pt] \vdots \\[6pt] \frac{\partial f}{\partial x_n} \end{bmatrix}. \]
6 範例 3：利用梯度下降法解決最小二乘法問題
已知：
\[ \begin{array}{ccccccc} \hline x_i & u_i \\ \hline 0 & 1 \\ 1 & 3 \\ 2 & 4 \\ 3 & 4 \\ \hline \end{array} \]
第一個最小二乘法問題的範例如下。對於每一筆數據 \((x_i, u_i)\)，令 \(\hat{u}_i\) 為將 \(x_i\) 代入線性函數 \(u = a + bx\) 所得的值。因此： \[ \hat{u}_i = a + bx_i. \] 如果這個方程式的解不存在，我們可以找出使得誤差平方 \((u_i - \hat{u}_i)^2\) 達到最小的 \(a, b\)。誤差函數 \(E(\mathbf{u})\) 定義為誤差平方和 \((u_i - \hat{u}_i)^2\)（包含 \(\frac{1}{2}\) 係數是為了計算方便，並不影響最終結果）：
\[\begin{align*} E(\mathbf{u}) & = E(a,b) = \frac{1}{2} \left\{ (a-1)^2 + (a+b-3)^2 + (a+2b-4)^2 + (a+3b-4)^2 \right\} \\ & = \frac{1}{2} \| A\mathbf{u} - \mathbf{y} \|^2. \end{align*}\]
讓我們使用 GDM 來找出使 \(E(\mathbf{u})\) 達到最小的近似解 \(\hat{\mathbf{u}}\)。
設定初始迭代值 \(\mathbf{u}_1 = (-2.5, -2.5)\)、容差 \(\epsilon = 10^{-6}\) 及初始學習率 \(\eta = 0.1\)。
用於最小二乘法問題的 GDM 參數
\[ \begin{array}{ccc} \textbf{參數 (Parameter)} & \textbf{數值 (Value)} \\ \hline \textrm{初始迭代值} ~\mathbf{u}_1 & (-2.5, -2.5) \\ \textrm{容差} ~ \epsilon & 10^{-6} \\ \textrm{學習率} ~ \eta & 0.1 \\ \hline \end{array} \]
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# ==============================
# 1. 定義誤差函數 E(a,b)
# ==============================
def error_func(a, b):
return 0.5 * (
(a - 1)**2 +
(a + b - 3)**2 +
(a + 2*b - 4)**2 +
(a + 3*b - 4)**2
)
# ==============================
# 2. 定義 E(a,b) 的梯度
# ==============================
def gradient(a, b):
dE_da = (a - 1) + (a + b - 3) + (a + 2*b - 4) + (a + 3*b - 4)
dE_db = (a + b - 3) + 2*(a + 2*b - 4) + 3*(a + 3*b - 4)
return np.array([dE_da, dE_db])
# ==============================
# 3. 梯度下降法參數 (與您的 LaTeX 相符)
# ==============================
u = np.array([-2.5, -2.5]) # 初始點
tol = 1e-6 # 容差
eta = 0.1 # 學習率
max_iter = 300
path = [] # 儲存迭代路徑
# ==============================
# 4. 執行梯度下降法
# ==============================
for k in range(max_iter):
grad = gradient(u[0], u[1])
grad_norm = np.linalg.norm(grad)
path.append((u[0], u[1], error_func(u[0], u[1])))
if grad_norm <= tol:
print("Algorithm Succeed!")
break
u = u - eta * grad## Algorithm Succeed!## u* = [1.49999929 1.00000033]## E(x*) = 0.500000000000364## iteration number = 118# ==============================
# 6. 繪圖 1: 3D 誤差曲面 + 下降路徑
# ==============================
a_vals = np.linspace(-3, 5, 50)
b_vals = np.linspace(-3, 5, 50)
A, B = np.meshgrid(a_vals, b_vals)
Z = error_func(A, B)
path_np = np.array(path)
fig = plt.figure(figsize=(12, 8))
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(A, B, Z, cmap='viridis', alpha=0.6, edgecolor='none')
ax.plot(path_np[:,0], path_np[:,1], path_np[:,2], 'r-', linewidth=2)
ax.scatter(path_np[:,0], path_np[:,1], path_np[:,2], c='black', s=15)
ax.set_xlabel('a')
ax.set_ylabel('b')
ax.set_zlabel('E(a,b)')
ax.view_init(30, 30)
ax.set_box_aspect([5, 5, 1])
plt.title("3D Error Function & Gradient Descent Path")
plt.show()
# ==============================
# 7. 繪圖 2: 最小二乘法直線擬合
# ==============================
x_data = np.array([0, 1, 2, 3])
y_data = np.array([1, 3, 4, 4])
x_fit = np.linspace(-0.5, 3.5, 100)
y_fit = u[0] + u[1] * x_fit
plt.figure(figsize=(8, 6))
plt.scatter(x_data, y_data, c='red', s=60, label='Data points')
plt.plot(x_fit, y_fit, 'b-', linewidth=2, label=f'Fit: y = {u[0]:.2f} + {u[1]:.2f}x')
plt.xlabel('x')
plt.ylabel('y')
plt.grid(alpha=0.3)
plt.legend()
plt.title('Least-Squares Linear Fit')
plt.xlim(-1, 4)## (-1.0, 4.0)## (-1.0, 5.0)
如我們稍早所見，我們透過演算法從給定數據中找出了最小二乘法直線。演算法求得了最小二乘法解，並以解答作為係數得出最小二乘法直線 \(y = \frac{3}{2} + x\)。
GDM 最小二乘法線性擬合總結
\[ \begin{array}{ccc} \textbf{參數 (Parameter)} & \textbf{數值 (Value)} \\ \hline \textrm{初始迭代值}~\mathbf{u}_1 & (-2.5, -2.5) \\ \textrm{容差}~ \epsilon & 10^{-6} \\ \textrm{學習率}~ \eta & 0.1 \\ \textrm{收斂解}~ \mathbf{u}^* & (1.5, 1.0) \\ \textrm{最小誤差}~ E(\mathbf{x}^*) & 0.5 \\ \textrm{收斂所需迭代次數}~ & 118 \\ \textrm{最小二乘法直線}~ & y = x + \frac{3}{2} \\ \hline \end{array} \]
最小二乘法直線解
最小二乘法直線為 \[ y = \frac{3}{2} + x. \]
因此，我們從最小二乘法解得到的直線是 \[ y = \frac{3}{2} + x. \]
7 範例 4：利用梯度下降法解決最小二乘法問題
我們也可以找出最符合給定數據的二次函數 \(y = a + bx + cx^2\)。在這種情況下，誤差函數 \(E(\mathbf{u})\) 如下： \[ \begin{aligned} E(a,b,c) &= \frac{1}{2}\left\{(a-1)^2 + (a+b+c-3)^2 + (a+2b+4c-4)^2 + (a+3b+9c-4)^2\right\} \\ &= \frac{1}{2} \| A\mathbf{u} - \mathbf{y} \|^2. \end{aligned} \]
最小二乘法問題總結
\[ \begin{array}{ccc} \textbf{問題類型 (Problem Type)} & \textbf{形式 / 結果 (Form/Result)} \\ \hline \textrm{線性最小二乘法模型} & y = a + bx \\ \textrm{線性最小二乘法解} & y = x + \frac{3}{2} \\ \textrm{二次最小二乘法模型} & y = a + bx + cx^2 \\ \textrm{二次擬合的誤差函數} & E(a,b,c) = \frac{1}{2}\|A\mathbf{u}-\mathbf{y}\|^2 \\ \hline \end{array} \]
import numpy as np
import matplotlib.pyplot as plt
# --------------------------
# 1. 定義誤差函數 E(a,b,c)
# --------------------------
def E(a, b, c):
return 0.5 * (
(a - 1)**2 +
(a + b + c - 3)**2 +
(a + 2*b + 4*c - 4)**2 +
(a + 3*b + 9*c - 4)**2
)
# --------------------------
# 2. 計算 E(a,b,c) 的梯度 (解析形式)
# --------------------------
def grad_E(a, b, c):
# 對 a 的偏導數
dE_da = (a - 1) + (a + b + c - 3) + (a + 2*b + 4*c - 4) + (a + 3*b + 9*c - 4)
# 對 b 的偏導數
dE_db = (a + b + c - 3) + 2*(a + 2*b + 4*c - 4) + 3*(a + 3*b + 9*c - 4)
# 對 c 的偏導數
dE_dc = (a + b + c - 3) + 4*(a + 2*b + 4*c - 4) + 9*(a + 3*b + 9*c - 4)
return np.array([dE_da, dE_db, dE_dc])
# --------------------------
# 3. 梯度下降演算法
# --------------------------
# 初始參數 (與您的程式碼相符)
u = np.array([1.0, 1.0, 1.0]) # 初始值 (a,b,c)
tol = 1e-6 # 容差等級
eta = 0.01 # 學習率
max_iter = 5000 # 最大迭代次數
for k in range(max_iter):
g = grad_E(u[0], u[1], u[2])
gn = np.linalg.norm(g)
if gn <= tol:
print("Algorithm Succeed!")
break
u = u - eta * g## Algorithm Succeed!## u* = [ 1.00000138 2.49999724 -0.49999918]## E(x*) = 1.5966473717177356e-12## iteration number = 4207# --------------------------
# 5. 繪製數據點與最小二乘法二次曲線
# --------------------------
# 原始數據點
x_data = np.array([0, 1, 2, 3])
y_data = np.array([1, 3, 4, 4])
# 產生曲線點
x_curve = np.linspace(-0.5, 3.5, 100)
y_curve = u[0] + u[1] * x_curve + u[2] * (x_curve ** 2)
# 繪圖
plt.figure(figsize=(8, 6))
plt.scatter(x_data, y_data, color='red', zorder=5, label='Data Points')
plt.plot(x_curve, y_curve, color='blue', label=f'Quadratic Fit: $y = {u[0]:.2f} + {u[1]:.2f}x + {u[2]:.2f}x^2$')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Least-Squares Quadratic Curve via Gradient Descent')
plt.legend()
plt.grid(alpha=0.3)
plt.ylim(-1, 5)## (-1.0, 5.0)## (-0.5, 4.0)