Cleve Moler 的複數步長數值微分方法

01-29

在沒有符號計算的環境下求函數的微分常用的方法是自動微分，這也是現代主流機器學習系統使用的方法。這裡介紹一種不同於此的數值微分方法，這種方法在 50 年前被 Cleve Moler，也就是 Mathworks 的創始人發明，並且現如今被使用在 MATLAB 和 Simulink 的組件 SimBiology 中。因為有一些限制所以實用性或許沒有自動微分那麼高，不過個人認為比較有趣而且原理簡單。

考慮一個解析函數 $F(z)$ ，這意味著該函數無限可微且能夠光滑的延拓的複平面上，那麼對於實數 $x_0$ 和 $h$ ，根據 Taylor 級數有：

$F(x_0+ih)=F(x_0)+ihF(x_0)-h^2F(x_0)/2!-ih^3F(x_0)/3!+...$

取虛數部分可以得到：

$F(x_0)=mathrm{Im}(F(x_0+ih))/h+O(h^2)$

由此我們可以計算 $F(x_0)$ ，其具有 $O(h^2)$ 階的精度，即

$F(x_0)approxmathrm{Im}(F(x_0+ih))/h$

舉一個例子，取

$F(x)=frac{e^x}{cos^3x+sin^3x}$

在 MATLAB 中的定義該函數並做圖：

F = @(x) exp(x) ./ (cos(x) .^ 3 + sin(x) .^ 3);nfplot(F, [-pi / 4, pi / 2]);naxis([-pi / 4, pi / 2, 0, 6])nset(gca, xtick, [-pi / 4, 0, pi / 4, pi / 2])nline(pi / 4, F(pi / 4), Marker, ., MarkerSize, 18)n

接著用我們上邊得到的方法求 $F(pi/4)$ ,也就是上圖中那點的導數值

Fpc = @(x, h) imag(F(x + 1i * h)) ./ h;n

作為對比，我們同樣用有限微分方法求導數值：

Fpd = @(x, h) (F(x + h) - F(x - h)) ./ (2 * h);n

並且用符號計算的結果用來驗證誤差：

syms Fps(x)nFps(x) = simplify(diff(F(x)));nexact = Fps(pi / 4);n

之後對 $h$ 取不同大小時兩種方法的誤差作圖：

h = 10 .^ (-(1 : 16));nerrs = [abs(Fpc(pi / 4, h) - exact) abs(Fpd(pi / 4, h) - exact)];nloglog(h, errs)nset(gca, xdir, rev)naxis([eps 1 eps 1])nlegend(complex step, finite difference, location, southwest)nxlabel(step size h)nylabel(error)n

可以看到有限微分方法在 $h$ 取較小值時，誤差反而會變大到不可接受的地步。

從例子可以看出，這種方法相對於有限差分法在步長較小時有更高的精度，相對於自動微分法的優勢則是幾乎不需要額外的編程，而是依賴於數值系統的複數值函數運算。當然該方法明顯的限制是要求函數為解析函數，且所在環境可以對函數的複數輸入求值，且只能求函數在實值位置的導數。另外，上述過程也只適用於求一階導數（有推廣求高階導數的方法，有興趣可以自己查閱）。

上述內容來自於 Cleve Moler 的博客：Complex Step Differentiation

根據相同原理寫了計算 Jacobian 矩陣的程序：Complex step Jacobian - File Exchange - MATLAB Central