前言

前幾天在翻舊書，看到國中以前做過的求根\極值方法: 牛頓法，以前都是做單變量函數的，現在發現其實有多變量的。研究後發現蠻有趣的，就寫了這篇。

主體

多變量牛頓法主要求函數極值，原理跟單變量差不多，只是換成更高維度的東西。

單變量的迭代公式: $x_{k+1}=x_k - \frac{f’(x)}{f’’(x)}$

多變量的迭代公式: $x_{k+1}=x_k - \nabla f(x_k)\mathbf H^{-1} f(x_k)$
這邊$x$是一堆變量$a, b, c …$。
$\nabla f(x_k)$是$f(x)$在$x_k$的梯度，$\mathbf H^{-1} f(x_k)$是$f(x)$在$x_k$的Hessian逆矩陣。
這兩個對應到單變量一階導跟二階導倒數。

梯度可以用不同的方向向量去定義，但是這裡將梯度定義為分別將$x$, $y$固定後的導數組成的向量，也就是偏導數組成的向量，例如:
$\nabla f(x, y) = \frac{\partial f}{\partial x}(x, y)i + \frac{\partial f}{\partial y}(x, y)j$

Hessian逆矩陣是從Hessian矩陣求逆得到的，Hessian矩陣是多變數函數的所有二階偏導數組成的方塊矩陣，例如:
$\mathbf H f(x, y) = \begin{bmatrix} \frac{\partial^2 f}{\partial x^2} & \frac{\partial^2 f}{\partial x \partial y} \\\ \frac{\partial^2 f}{\partial y \partial x} & \frac{\partial^2 f}{\partial y^2} \end{bmatrix}$

實作想法

考慮到視覺化的需求，用python來寫最合適。

因為只討論雙變量函數，所以Hessian矩陣只需要$2\times2$的就好，可以手動設。

而numpy提供的np.linalg.solve()也恰好可以求Hessian逆矩陣與梯度矩陣的乘積。

接著是計算偏導數，
如果你有看過我寫的曲線插值-樣條函數這篇，或許會發現我的微分是用原始定義去做的 $lim_{h \to 0} \frac{1}{h}[f(x + h) - f(x)]$ ，將$h$設成一個很小的數字，就可以逼近我們要的微分數值。
同樣的道理，偏微分也可以照著做:

$f_x(x, y) = \frac{\partial f}{\partial x}(x, y) = lim_{h \to 0} \frac{1}{h}[f(x + h, y) - f(x, y)]$

換成程式:

1
2
3

def dfdx(x0, y0):
    h = 1e-5
    return (f(x0 + h, y0) - f(x0, y0)) / h

二階偏導同樣:

$\frac{\partial}{\partial x}(\frac{\partial f}{\partial x}(x, y)) = \frac{\partial}{\partial x}f_x(x, y)$

$=lim_{g \to 0} \frac{1}{g}[f_x(x + h, y) - f_x(x, y)]$

1
2
3

def dfdxdx(x0, y0):
    g = 1e-5
    return (dfdx(x0 + g, y0) - dfdx(x0, y0)) / g

其餘偏導類似。

code

import numpy as np
import math
import matplotlib.pyplot as plt

f = lambda x, y: x * np.exp(-x * x - y * y)

eps = 1e-6

# 偏導數
def dfdx(x0, y0):
    h = 1e-5
    return (f(x0 + h, y0) - f(x0, y0)) / h

def dfdy(x0, y0):
    h = 1e-5
    return (f(x0, y0 + h) - f(x0, y0)) / h

def dfdxdx(x0, y0):
    g = 1e-5
    return (dfdx(x0 + g, y0) - dfdx(x0, y0)) / g

def dfdxdy(x0, y0):
    g = 1e-5
    return (dfdx(x0, y0 + g) - dfdx(x0, y0)) / g

def dfdydy(x0, y0):
    g = 1e-5
    return (dfdy(x0, y0 + g) - dfdy(x0, y0)) / g

# Hessian矩陣
def H(x0, y0):
    h = [[dfdxdx(x0, y0), dfdxdy(x0, y0)], [dfdxdy(x0, y0), dfdydy(x0, y0)]]
    return h

# 梯度矩陣
def grd(x0, y0):
    return [dfdx(x0, y0), dfdy(x0, y0)]

X = np.linspace(-5, 5, 1024)
Y = np.linspace(-5, 5, 1024)

X, Y = np.meshgrid(X, Y)
Z = f(X, Y)
plt.xlabel("x")
plt.ylabel("y")

# 等高線圖
C = plt.contour(X, Y, Z)
plt.clabel(C, inline = True)

# 初始值
x0 = 1
y0 = 0.2

# 求Hessian逆矩陣與梯度矩陣的乘積
k = -np.linalg.solve(H(x0, y0), grd(x0, y0))

# 先迭代一次
x1 = x0 + k[0]
y1 = y0 + k[1]

# 畫線
plt.plot((x0, x1), (y0, y1), 'bo-')

for i in range(10):
    x0 = x1
    y0 = y1
    k = -np.linalg.solve(H(x0, y0), grd(x0, y0))
    x1 = x0 + L * k[0]
    y1 = y0 + L * k[1]
    plt.plot((x0, x1), (y0, y1), 'ro-')

# 最終值
print(x0, y0)
print("extreme:", f(x0, y0))

測試

$(ln(|x - 1|))(5y - 1) - 1$, 起始值=(1.1, 9), 迭代次數 6
收斂值 $(x,y) = (1.9999999999995695, 0.20000000000075954), f(x, y) = -1.0$

$cos(\sqrt{x^2 + y^2})$, 起始值=(0.5, 1), 迭代次數 6
收斂值 $(x,y) = (-4.999999323762127e-06, -5.000008901435384e-06), f(x, y) = 0.999999999975$

$xe^{-x^2-y^2}$, 起始值=(1, 0.2), 迭代次數 7
收斂值 $(x,y) = (0.7071017811951483, -4.999999497366208e-06), f(x, y) = 0.4288819424481872$

$xe^{-x^2-y^2}$, 起始值=(1, 0.3), 迭代次數 21
收斂值 無法收斂

$sin(x) \times y^{\frac{2}{3}}$, 起始值=(-2, 0.9), 迭代次數 15
收斂值 $(x,y) = (-7.683814908002442e-05, -0.0020110391746036667), f(x, y) = -1.2242137360945188e-06$

$sin(x) \times y^{\frac{2}{3}}$, 起始值=(-1, 2), 迭代次數 21
收斂值 $(x,y) = (-20.42035724833452 -4848738.969554567), f(x, y) = -28647.443380356$
註解: 收斂到奇怪的地方了…

結論

看的出來，多變量牛頓法其實就跟單變量牛頓法有著一樣的缺點: 不穩定、起始點敏感
更糟糕的情況，Hessian矩陣可能沒有逆，例如:$f(x, y) = ln(x^2 + y^2)$。

但同時優點也很明顯: 收斂速度快 (上面的迭代次數有調高，以求比較精準的收斂值)，比較簡單的函數甚至可以兩步到位。

個人覺得實際應用不可能直接用純的牛頓法，弊大於利，畢竟起始點必須選的很好才有可能收斂，應該搭配其他演算法抓起始點的大略位置，再給牛頓法去迭代。

結語

這次的東西用到好多沒看過的數學方法，不過還好不太深，還可以理解，但是僅僅能寫出陽春版的牛頓法，貌似還可以通過Hessian矩陣是否正定來判斷狀態?
待研究…

Reference

https://zh.wikipedia.org/zh-tw/%E7%89%9B%E9%A1%BF%E6%B3%95
https://en.wikipedia.org/wiki/Newton's_method_in_optimization
https://en.wikipedia.org/wiki/Gradient
https://en.wikipedia.org/wiki/Hessian_matrix
https://zhuanlan.zhihu.com/p/46536960
https://zhuanlan.zhihu.com/p/37524275
https://www.zhihu.com/question/19723347