<轉>如何用C++實現自動微分


作者:李瞬生
轉摘鏈接: https://www.zhihu.com/question/48356514/answer/123290631
來源:知乎
著作權歸作者所有。

實現 AD 有兩種方式,函數重載代碼生成。兩種方式的原理都一樣,鏈式法則

不難想象,任何計算都可以由第1步到第k步的序列形式,其中第 i 步計算的輸入,在之前的 i-1 步中已經計算(例如編譯器生成的匯編指令序列)。因此,任何計算都可以看作形式如下圖左側的復合函數。微積分中的鏈式法則告訴我們,符合函數的導數可寫作下圖右側的形式(假設每一步都可導)。請注意偏導數全導數的區別。

圖一

自動微分的第一個難點就在這,微積分中的鏈式法則。大家在課堂上學的鏈式法則的示例通常只有兩到三個函數,而自動微分面對的計算,有無數個函數。許多人不習慣在這樣大的規模上應用鏈式法則。不過一旦習慣,就會發現自動微分的原理十分簡單。

如果上述內容過於抽象,請參看下面這個例子以后再看一遍。

圖二

在上圖中,頂部方程是二維旋轉。輸入 x 有三個變量,旋轉角度及二維坐標。輸出 z 有兩個變量,即旋轉后的二維坐標。上圖列表第二列是該計算的序列形式,第一列是每一步對應的表達式,第三列是對應的鏈式法則(請對比圖一)。

太繁瑣了?看不出個所以然?不用擔心。
如果把計算的序列形式及其導數計算每個步驟的依賴關系表示成圖

圖三

不難發現,兩張圖是等價的。也就是說,計算序列形式的每一步都與其導數計算的步驟有一一對應的關系。源程序怎么算,其導數就可以怎么算(從順序上來說)。

以上便是自動微分的基本原理。下面我們來談實現。
如圖二,圖三,我們有兩種方式來考慮自動微分的實現。
  1. 用戶提供圖二第二列序列形式的源代碼,按順序生成第三列的微分計算。此種方法的特點是,讀一行源代碼,生成一行微分計算,因此可以動態生成。
    若源代碼這一行在做乘法,那么就依據乘法法則生成該步的微分計算。若源代碼這一行是三角函數 cos(x),那么它對應的微分計算就是 -sin(x),以此類推。每一步計算的偏導數都根據鏈式法則組合,得出該步驟的全導數。
    該種方法的常見手段是函數重載。優點是簡單直接,缺點是動態生成成本較高。
  2. 用戶提供源代碼,在編譯時生成圖三左側的程序結構圖,並生成圖三右側對應的微分程序。
    該種方法的常見手段是編譯時的代碼生成(比如用 flex-bison 做詞法、語法分析)。優點是靜態生成效率高,一次生成,多次使用。缺點是編譯原理有門檻,非計算機專業望而卻步。

力求簡單,下面給出以函數重載實現的簡單示例 SAD - Simple Automatic Differentiation。請對照 圖二中的列表閱讀。
main 函數中對應圖二表中第二列,2D 旋轉的計算。ADV 里重載的 +, -, *, sin, cos 不僅完成本來的計算,還負責圖二表中第三列導數的計算。
main 函數中 x, y, z 的序號都與圖二中對應。
 1 #include <cmath>
 2 #include <vector>
 3 #include <iostream>
 4 namespace SAD  // Simple Automatic Differentiation
 5 {
 6     class ADV
 7     {
 8     public:
 9         ADV(double v = 0, double d = 0);
10 
11         // overloaded unary and binary operators
12         ADV operator + (const ADV &x) const;
13         ADV operator - (const ADV &x) const;
14         ADV operator * (const ADV &x) const;
15         friend ADV sin(const ADV &x);
16         friend ADV cos(const ADV &x);
17 
18         double val;  // value of the variable
19         double dval;  // derivative of the variable
20     };
21 
22     ADV::ADV(double v, double d) : val(v), dval(d) {}
23 
24     ADV ADV::operator+(const ADV &x) const
25     {
26         ADV y;
27         y.val = val + x.val;
28         y.dval = dval + x.dval;
29         return y;
30     }
31 
32     ADV ADV::operator-(const ADV &x) const
33     {
34         ADV y;
35         y.val = val - x.val;
36         y.dval = dval - x.dval; // sum rule
37         return y;
38     }
39 
40     ADV ADV::operator*(const ADV &x) const
41     {
42         ADV y;
43         y.val = val*x.val;
44         y.dval = x.val*dval + val*x.dval; // product rule
45         return y;
46     }
47 
48     ADV sin(const ADV &x)
49     {
50         ADV y;
51         y.val = std::sin(x.val);
52         y.dval = std::cos(x.val)*x.dval; // chain rule
53         return y;
54     }
55 
56     ADV cos(const ADV &x)
57     {
58         ADV y;
59         y.val = std::cos(x.val);
60         y.dval = -std::sin(x.val)*x.dval; // chain rule
61         return y;
62     }
63 }
64 
65 
66 int main()
67 {
68     using namespace SAD;
69     using namespace std;
70 
71     static const double PI = 3.1415926;
72     vector<ADV> x;
73 
74     x.emplace_back(PI, 1);  // x = [PI, 2, 1]
75     x.emplace_back(2, 0);
76     x.emplace_back(1, 0);
77 
78     ADV y1 = cos(x[0]);
79     ADV y2 = sin(x[0]);
80     ADV y3 = x[1] * y1;
81     ADV y4 = x[2] * y2;
82     ADV y5 = x[1] * y2;
83     ADV y6 = x[2] * y1;
84 
85     ADV z1 = y3 + y4;
86     ADV z2 = y6 - y5;
87 
88     cout << "x = [" << x[0].val << ", " << x[1].val << ", " << x[2].val << "]" << endl;
89     cout << "z = [" << z1.val << ", " << z2.val << "]" << endl;
90     cout << "[dz1/dx0, dz2/dx0] = [" << z1.dval << "," << z2.dval << "]" << endl;
91 }

 

運行結果:

矢量 [2,1] 被旋轉 180° , 變為 [-2,-1]。關於角度的導數為 [-1,2]。


自動微分的經典教材是該題目的奠基人 Griewank 著的 Evaluating Derivatives (Society for Industrial and Applied Mathematics)
該書囊括了自動微分的所有方面,比如本文未介紹的 reverse mode, sparse Jacobian, Hessian 等。
如果不求全面,一本更通俗更面向代碼實現的書是 The Art of Differentiating Computer Programs (Society for Industrial and Applied Mathematics)

最后,自動微分是算導數的最優方法,比符號計算、有限微分更快更精確。
自動微分已經廣泛應用在優化領域,包括人工神經網絡的訓練算法 back-propagation。
要解連續優化或非線性方程,自動微分是不二的選擇。


幾何福利(自動微分的張量推導)
剛才是科普,下面給出自動微分的張量公式,及其對應的,便於編程實現的矩陣公式。

圖一的計算序列可以記作如下形式
\begin{eqnarray*}
y^{1} &=&f^{1}\left( x\right)  \\
y^{i} &=&f^{i}\left( x,y^{1},\cdots ,y^{i-1}\right) ,i=2,\cdots ,n \\
z &=&\bar{f}\left( x,y^{1},\cdots ,y^{n}\right) 
\end{eqnarray*}
其中x是自變量 y是中間變量 z是因變量 他們可以是任意維度

微分流形上的矢量有兩種:切空間 (Tangent space)的矢量,對偶空間 (Dual space)的矢量。
其中切空間的\left\{ \dfrac{\partial }{\partial x^{i}} \right\} 是關於坐標系的偏導算符
對偶空間的是關於坐標系的微分\left\{ dx^i \right\}

從切矢量出發 我們可以得到自動微分的正序模式 (forward mode)
從對偶矢量出發 我們可以得到自動微分的逆序模式 (reverse mode)

任意切矢量 \dfrac{d}{ds}=\dfrac{dx^{i}}{ds}\dfrac{%
\partial }{\partial x^{i}} 的定義是其對應的方向導數算符
將它依次應用在計算序列的左邊 便可獲得下圖左側的張量公式

\[
\begin{tabular}{c|c}
Tensor Form & Matrix Form \\ \hline\hline
$\dfrac{dy^{1k}}{ds}=\dfrac{\partial f^{1k}}{\partial x^{j}}\dfrac{dx^{j}}{ds%
}$ & $D_{s}^{1}=P_{x}^{1}d_{s}$ \\ 
$\dfrac{dy^{ik}}{ds}=\dfrac{\partial f^{ik}}{\partial x^{j}}\dfrac{dx^{j}}{ds%
}+\dfrac{\partial f^{ik}}{\partial y^{jh}}\dfrac{dy^{jh}}{ds}$ & $%
D_{s}^{i}=P_{x}^{i}d_{s}+\sum_{j<i}P_{^{j}}^{i}D_{s}^{j}$ \\ 
$\dfrac{dz^{i}}{ds}=\dfrac{\partial \bar{f}^{i}}{\partial x^{j}}\dfrac{dx^{j}%
}{ds}+\dfrac{\partial \bar{f}^{i}}{\partial y^{jh}}\dfrac{dy^{jh}}{ds}$ & $%
J_{s}=P_{x}^{z}d_{s}+\sum_{j<i}P_{j}^{z}D_{s}^{j}$%
\end{tabular}%
\]

對張量縮並可得上圖右側的矩陣公式 其中 J 便是一階導 Jacobian 矩陣
這種計算導數的方式與計算序列同序 故名正序模式
每一個方向導數的計算復雜度與計算序列相同 空間復雜度也相同
注意 若計算序列的自變量有 n 維 則獲得 Jacobian 矩陣需要計算 n 個方向上的方向導數

任意對偶矢量 d_{t}=t_{i}dz^{i} 的定義是其對應方向的微分
對該矢量做關於坐標基 \left\{ dx^i \right\} 的坐標變換 便可獲得下圖左側的張量公式

\[
\begin{tabular}{c|c}
Tensor Form & Matrix Form \\ \hline\hline
\multicolumn{1}{c|}{$d_{t}=t_{i}\dfrac{\partial z^{i}}{\partial x^{j}}%
dx^{j}+t_{i}\dfrac{\partial z^{i}}{\partial y^{jl}}dy^{jl}$} & $%
D_{x}+=TP_{x}^{z},D_{j}+=TP_{j}^{z}$ \\ 
\multicolumn{1}{c|}{$dy^{ik}=\dfrac{\partial y^{ik}}{\partial x^{j}}dx^{j}+%
\dfrac{\partial y^{ik}}{\partial y^{jl}}dy^{jl}$} & $%
D_{x}+=D_{i}P_{x}^{i},D_{j}+=D_{i}P_{j}^{i}$ \\ 
\multicolumn{1}{c|}{$dy^{1k}=\dfrac{\partial y^{1k}}{\partial x^{j}}dx^{j}$}
& $D_{x}+=D_{1}P_{x}^{1}$%
\end{tabular}%
\]

對漲量縮並可得上圖右側的矩陣公式 其中 D += TP 表示疊加 D = D + TP
D_x 便是關於x的一階導
這種計算導數的方式與計算序列逆序 故名逆序模式
展開每一個基底 dz^i 的計算復雜度與計算序列相同 但由於內存訪問的順序與計算序列相逆 所有中間結果都需要保存下來 因而空間復雜度與計算復雜度相當
若計算序列的因變量有 m 維 則獲得 Jacobian 矩陣需要展開 m 個基底 dz^i


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM