自動微分

在數學和電腦代數中，自動微分有時稱作演算式微分，是一種可以藉由電腦程式計算一個函數導數的方法。兩種傳統做微分的方法為：

- 對一個函數的表示式做符號上的微分，並且計算其在某一點上的值。
- 使用差分。

使用符號微分最主要的缺點是速度慢及將電腦程式轉換成表示式的困難。此外，很多函數在要計算更高階微分時會變得複雜。

使用差分的兩個重要的缺點是捨棄誤差及數值化過程和相消誤差。此兩者傳統方法在計算更高階微分時，都有複雜度及誤差增加的問題。自動微分則解決上述的問題。

自動微分使用這個事實：任何實作一個向量函數 y=F(x)的電腦程式，一般而言，可以被分解成由基本指定運算所成的序列，而其中每一個都可以藉由查表而輕易地微分。這些計算某一特定項的 "基本偏微分" 是依照微積分中的複合函數求導法則來合併成某個 F 的微分資訊（如梯度、切線、雅可比矩陣等）。這個過程會產生確實(數值上準確)的導數。因為只在最基礎的層面做符號轉換，自動微分避免了複雜的符號運算的問題。

## 複合函數求導法則，前向和反向積累

自動微分的基礎是，根據複合函數求導法則來合併微分值。以$f(x)=g(h(x))$為例，根據複合函數求導法則，我們有：

$$\frac{df}{dx} = \frac{dg}{dh} \frac{dh}{dx}$$

通常有兩個不同的模式：「前向積累」（或「前向模式」）和「反向積累」（或反向模式）。

前向積累由右到左地使用複合函數求導法則，即先計算 $dh/dx$，然後才$dg/dh$。

反向積累則是由左到右。

### 前向積累

前向積累式的自動微分是最容易理解和實作的。

$f(x_1,x_2) = x_1 x_2 + \sin(x_1)$這個函數是可被電腦（或程式設計師）

解釋成一連串對變數$w_i$的運算。

前向積累式自動微分的工具則會增加相對應的作用於第二項上的運算。

原本程式碼敘述｜為導數而新增的敘述

$w_1 = x_1$｜$w'_1 = 1$ (種子)

$w_2 = x_2$｜$w'_2 = 0$ (種子)

$w_3 = w_1 w_2$｜$w'_3 = w'_1 w_2 + w_1 w'_2 = 1 x_2 + x_1 0 = x_2$

$w_4 = \sin(w_1)$｜$w'_4 = \cos(w_1)w'_1 = \cos(x_1) 1$

$w_5 = w_3 + w_4$｜$w'_5 = w'_3 + w'_4 = x_2 + \cos(x_1)$

計算$f(x_1,x_2) = x_1 x_2 + \sin(x_1)$的導數需要初始化。

以區別是要對$x_1$或$x_2$來求導數。

上述表格則以$w'_1=1$和$w'_2=0$來初始化。

並且我們可以看到其結果$x_2 + \cos(x_1)$正是對$x_1$的導數。

注意，雖然表格列出了符號微分。

但以電腦的角度而言，電腦總是儲存數值。圖2則以圖形表明上述的敘述。

為了計算這個例子的導數，其分別為$\partial f/\partial x_1$和$\partial f / \partial x_2$。

需要計算兩次，一次是以$w'_1 = 1$和$w'_2 = 0$做初始值。

另一次則以$w'_1 = 0$和$w'_2 = 1$做為初始值。

前向積累的計算複雜度則正比於原來程式的計算複雜度。

對於函數$f:\mathbb{R} \rightarrow \mathbb{R}^m$且 $m \gg 1$來說。

前向積累只要計算一次，優於需要計算 m 次的反向積累。

## 前向式積累自動微分的由來與二元數

前向式積累自動微分可藉由擴充實數中的代數並得到一個新的算術系統來達成。

每一個數都會新增另一數，用來表示一函數在這數上導數的數。

而每一個算術運算都被擴充於此新的代數。

這個擴充後的代數就是二元數的代數。

將每一個數$\,x$替換成數$x+x'\varepsilon$，其中$x'$是一個實數，但$\varepsilon$則只是一個據有$\varepsilon^2=0$這個特性的符號。

使用這特性，我們可以有運算

$$(x + x'\varepsilon) + (y + y'\varepsilon) = x + y + (x' + y')\varepsilon$$

$$(x + x'\varepsilon) \cdot (y + y'\varepsilon) = xy + xy'\varepsilon + yx'\varepsilon + x'y'\varepsilon^2 = xy + (x y' + yx')\varepsilon$$

減法和除法則類似。

現在，我們可以計算多項式。

如果$P(x) = p_0 + p_1 x + p_2x^2 + \cdots + p_n x^n$，則

$P(x + x'\varepsilon)$｜$=\,$｜$p_0 + p_1(x + x'\varepsilon) + \cdots + p_n (x + x'\varepsilon)^n$

$=\,$｜$p_0 + p_1 x + \cdots + p_n x^n$

$$\, {} + p_1x'\varepsilon + 2p_2xx'\varepsilon + \cdots + np_n x^{n-1} x'\varepsilon$$

$=\,$｜$P(x) + P^{(1)}(x)x'\varepsilon$

其中 $P^{(1)}$ 表示 $P$ 對第一個參數的導數。

而 $x'$ 則稱作「種子」，可以任意選擇。

新的算術是由有序對、寫成$\langle x, x' \rangle$ 及對第一項的運算和對第二項的第一階微分運算所組成。

將上述結果應用於多項式的解析函數上。

我們可以得到一系列關於基本算術和一些標準函數的新算術：

$$\langle u,u'\rangle +\langle v,v'\rangle = \langle u+v, u'+v' \rangle$$

$$\langle u,u'\rangle -\langle v,v'\rangle = \langle u-v, u'-v' \rangle$$

$$\langle u,u'\rangle *\langle v,v'\rangle = \langle u v, u'v+uv' \rangle$$

$$\langle u,u'\rangle /\langle v,v'\rangle = \left\langle \frac{u}{v}, \frac{u'v-uv'}{v^2} \right\rangle \quad ( v\ne 0)$$

$$\sin\langle u,u'\rangle = \langle \sin(u) , u' \cos(u) \rangle$$

$$\cos\langle u,u'\rangle = \langle \cos(u) , -u' \sin(u) \rangle$$

$$\exp\langle u,u'\rangle = \langle \exp u , u' \exp u \rangle$$

$$\log\langle u,u'\rangle = \langle \log(u) , u'/u \rangle \quad (u>0)$$

$$\langle u,u'\rangle^k = \langle u^k , k u^{k-1} u' \rangle \quad (u \ne 0)$$

$$\left| \langle u,u'\rangle \right| = \langle \left| u \right| , u' \mbox{sign} u \rangle \quad (u \ne 0)$$

並且，一般而言，對於一個函數 $g$，我們會有

$$g(\langle u,u' \rangle , \langle v,v' \rangle ) = \langle g(u,v) , g^{(1)}(u,v) u' + g^{(2)}(u,v) v' \rangle$$

其中$g^{(1)}$和$g^{(2)}$分別是$g$對其第一項和第二項的導數。

對一個二元算術運算作用於混合的參數時（數對$\langle u, u' \rangle$和實數$c$）。

實數會先被轉成$\langle c, 0 \rangle$。

函數$f : \mathbb{R}\rightarrow\mathbb{R}$在$x_0$上的導數則為以上述算術計算$f(\langle x_0, 1 \rangle)$，其結果為$\langle f ( x_0 ) , f' ( x_0 ) \rangle$。

### 向量參數和函數

藉由採取方向導數的運算。

多變數函數也可以同單變數函數的效率和機制來處理。

亦即，函數$f:\mathbb{R}^n\rightarrow\mathbb{R}^m$在$x \in \mathbb{R}^n$這點。

和$x' \in \mathbb{R}^n$這個方向上的方向導數$y' \in \mathbb{R}^m$。

可以使用上述相同的算術來計算$$(\langle y_1,y'_1\rangle, \ldots, \langle y_m,y'_m\rangle) = f(\langle x_1,x'_1\rangle, \ldots, \langle x_n,x'_n\rangle)$$ 而求得。

### 更高階微分

以上的算術可以被一般化，以用於二階及三階導數。

然而，此算術的規則將會迅速變得複雜。

其複雜度將與最高階導數階數成平化。

取而代之的是使用限縮泰勒級數。

這是可行的，因為函數的泰勒級數中的通項為己知係數和函數導數的乘積。

使用自動微分來計算黑塞矩陣在某些最佳化已被證明是可行的。

## 實作

前向式積累是由對程式的非標準化轉譯程序來實作。

即將實數替換成二元數，常數則換成有第二項為零係數的二元數。

而數值上基本運算則被換成二元數的運算。

非標準化轉譯程序一般使用兩者策略之一：程式碼轉換和運算符重載。

### 程式碼轉換

<image>

一個函數的程式碼會被自動產生的程式碼所替換。

新生成用來計算導數的程式碼則會插入原程式碼中。

程式碼轉換可實作在所有的程式語言上，且它對編譯器而言，是容易最佳化的。

然而，實作自動微分的工具則是比較困難的。

例子：

- ADIC (C/C++, 前向積累)
- ADIFOR (Fortran77)
- OpenAD (Fortran77, Fortran95, C/C++)
- TAPENADE (Fortran77, Fortran95)

### 運算符重載

<image>

如果所使用的程式語言支持，運算符重載是個可行的方法。

實數的物件跟基本數學運算必須重載以滿足上述 augmented 算術。

這不須要改變要被微分的函數的程式碼。

運算符重載對前向積累是容易實作的，並且可能對反向積累亦如此。

然而，與前向積累相比，現有的編譯器在最佳化程式碼方面則是較為落後。

例子：

- ADC Version 4.0 (C/C++)
- ADF Version 4.0 (Fortran 77, Fortran 95)
- ADOL-C (C/C++)
- FADBAD++ (C/C++)
- CppAD (C/C++)
- MAD (Matlab)
- Sacado (Trilinos的一部分) (C++, forward/reverse)
- JAX（Python，自動微分、JIT編譯）

## 外部連結

- www.autodiff.org , An "entry site to everything you want to know about automatic differentiation"
- Automatic Differentiation of Parallel OpenMP Programs
- Automatic Differentiation, C++ Templates and Photogrammetry
- Automatic Differentiation, Operator Overloading Approach
- Compute analytic derivatives through a web-based interface Automatic differentiation of nonlinear models
- Compute analytic derivatives of any Fortran77 or Fortran95 through a web-based interface Automatic Differentiation of Fortran programs