allenlu2007

A great WordPress.com site

Category: math

平行公理和微分幾何平行移動(Parallel Transport)

by allenlu2007

垂直和平行都是幾何最重要的基石。

  • 垂直的局部性質很直接,兩直線相交,夾角為 90 為垂直。一般人隨手可以畫出垂直相交線,即使在曲面也不例外。垂直的微分定義是兩直線切向量內積為 0.
  • 平行的概念看似直接,兩直線永不相交。在歐幾里德平面幾何很容易想像。兩個問題:(A) 這樣平行直線唯一嗎?或一定存在嗎? (B)平行是否有局部的微分定義
  • (A) 的回應就是平行公理;(B)的回應就是平行移動 (parallel transport).

歐幾里得平面幾何五大公理。最後一個公理就是平行公理:一條直線線外一點只能做出唯一平行直線

俄國數學家羅巴切夫斯基 (Lobachevsky) 一輩子想要其他公理證明平行公理或者用更直觀的公理取代平行公理都未能成功。後來用反證法:一條直線外一點一定會有兩條或兩條以上的平行直線,試圖找出矛盾。但卻發現邏輯自洽。因而開創非歐幾何學。我們現在知道這其實對應雙曲面幾何如下。這是一個馬鞍面,立起來就是一個雙曲面。雙曲面上的"直線“線外一點可以找到無數 diverged “直線”,和原來的“直線”永不相交,因此可以有無數條平行線。

黎曼則另闢蹊徑開拓黎曼幾何學。對應(橢圓)球面幾何。球面上所有的“直線”都是 中心在球心的封閉大圓。因此所有的“直線”都相交。平行公理變成一條直線外一點沒有任何平行直線。常見的誤區是球面的緯線是平行線。緯線除了赤道是“直線”,其他的緯線都是球面上的“曲線”,不能視為平行直線!

下圖很好的總結三種幾何學。

曲面幾何和曲率

曲面幾何可以定義曲率。曲率乍看之下是把一個曲面放在更高維空間(外)顯出的性質。

例如上圖是把二維的曲面放在(嵌入)三維的歐氏空間,可以明顯看出不同的曲面。同時可以定量的描述這個曲面的曲率。例如球面的曲率非常直觀和球的半徑相反。越小的球,曲率越大。但球面曲率是一個整體的性質,還是每一點(和其鄰近點)的性質?曲率是和半徑成反比,還是和半徑平方成反比?

先回顧二維歐氏平面上的一維曲線。很明顯曲線上每一點(和其鄰近點)都可以定義“線曲率” \kappa 是相切圓(不是拋物線!)半徑的倒數,i.e. \kappa = 1/r . 當然可以根據相切圓是在曲線的上方或是下方定義 \kappa 是正或是負。\kappa 明顯和泰勒展開的二階導數成正比1。直線的二階導數為 0, 曲率為 0.

再回到二維曲面,曲面上每一點都有無限條曲面上的曲線相交於這一點。每一條曲線都有本身的“線曲率”,可以是正或負。高斯曲率定義二維曲面每一點的“曲面曲率” R 是相交該點所有曲線(更正確應該是 geodesics)最大的線曲率乘最小的線曲率(含正負號), i.e. R = \kappa_1 \kappa_2 ; \kappa_1 是最大線曲率;\kappa_2 是最小線曲率。

以完美球面為例,每一點所有相交曲線的線曲率都一樣。都是大圓(圓心即是球心)半徑的倒數。所以球面每一點的高斯曲率 R = \kappa_1 \kappa_2 = 1/r^2 . 這符合直覺。

如果是非完美球面,例如橢圓球面半徑為 (a, b, c) , 曲面每一點的高斯曲率都不同,還是可以找到對應的 \kappa_1 \,, \kappa_2 . (完美或非完美)球面上同一點 \kappa_1 \,, \kappa_2 一定同正或同負,所以球面幾何每一點的高斯曲率永遠為正

再考慮雙曲面幾何。每一點相交曲線的”線曲率“必定有正有負,\kappa_1 > 0 \,, \kappa_2 < 0 . 所以雙曲面幾何每一點的高斯曲率永遠為負

外顯和內稟曲率

So far so good!
檢查高斯曲率在平面和圓柱面的曲率。平面是 trivial case, 所有 \kappa=0 . 圓柱面就很有意思,x,y 軸切面是圓形,半徑為 d; z 軸切面是(兩條)直線。帶入高斯曲率,\kappa_1 = 1/d \,,\, \kappa_2=0 , R = \kappa_1 \kappa_2 = 0 . 也就是圓柱面每一點的高斯曲率為 0, 和平面一模一樣!

顯然這不直觀。直觀上會認為圓柱面和平面不同,具有曲率,至少在 x,y 方向有曲率。為什麼高斯會定義高斯曲率為 \kappa_1 \kappa_2 而非 (\kappa_1 + \kappa_2)/2 ? (事實上,每一點的平均曲率定義為 mean curvature = (\kappa_1 + \kappa_2)/2 )。

圓柱面是一個很好的例子,說明外顯曲率和內稟曲率的不同

  • 外顯曲率(extrinsic curvature)就是把曲面放在更高維的歐氏空間,從更高維歐氏空間測量和定義的曲率。
  • 內稟曲率(intrinsic curvature)就是在(比嵌入歐氏空間低維)曲面本身測量和定義的曲率。更接地氣的說法,就是在曲面上的生物,可以測量和定義的曲率。
  • 舉例而言,上圖二維歐氏空間上的一維曲線,對於生活在曲線上一隻只會前進和後退(無限小)的蚯蚓,它能測量每一段曲線的長度,但是否能測量曲線上每一點的(內稟)曲率?直觀上顯然不能。唯一的機會是封閉曲線。如果蚯蚓往一個方向前進,最後回到同一點,代表這是一條曲線。假設曲線是圓形,甚至可以算出整條曲線的平均直徑,以及整條線平均曲率。實務上沒有用途:(1)大多數的曲線都是開放而非封閉曲線。(2)即使是封閉曲線,很少是圓形,算出平均半徑,或是平均曲率沒有意義。注意曲率是每一點局部的特性。除非是完美圓形或球形,一條曲線的平均曲率沒有太大意義。上文的 mean curvature 是每一點相交曲線的平均曲率,是有意義。
  • 簡單而言,一維曲線的“線曲率” \kappa 是嵌入在二維或是更高維歐氏空間,所以“線曲率”是外顯曲率。
  • 再看二維曲面,上文由“線曲率”所定義的高斯曲率 R = \kappa_1\kappa_2 或是 mean curvature = (\kappa_1+\kappa_2)/2 應該也是外顯曲率。NO!
  • 結果很意外。高斯曲率是內稟曲率,mean curvature 是外顯曲率。
  • 對於二維曲面上的一隻螞蟻,不用線曲率,它能測量每一點的曲率嗎?能!見下文。具體的方法是在每一點畫一個圓(或三角形),測量圓的圓周率(或是三角形內角和)。就可以得到該點的高斯曲率 R , 結果剛好等於最大和最小線曲率的乘積 R = \kappa_1\kappa_2 . 高斯自己似乎也覺得不可思議,稱為"絕妙定理“[@wikiTheoremaEgregium2019] 如下。
  • Gaussian curvature can be determined entirely by measuring angles, distances and their rates on a surface, without reference to the particular manner in which the surface is embedded in the ambient 3-dimensional Euclidean space. In other words, the Gaussian curvature of a surface does not change if one bends the surface without stretching it. Thus the Gaussian curvature is an intrinsic invariant of a surface.
  • 所以平面和圓柱面的高斯曲率是相同的。一隻曲面上的螞蟻無法分辨(locally)是一個平面或是圓柱面。因為圓柱面展開是一個平面。但是 mean curvature 兩者不同。因為 mean curvature 是外顯曲率。
  • 兩個曲面如果每一點高斯曲率相同,則為 isometric. 反之,無法 global isometric, 最多是 local isometric. 球面和平面無法等距 mapping. 只要看地圖就知道。赤道部分比例可以很準,但是兩極比例差很大。
  • 三維空間如何定義“內稟“曲率?因為很難直觀想像三維彎曲空間嵌入一個四維空間。三維彎曲空間包含三個 basis (x, y, z), 每一點(和鄰近點形成的空間)可以分解成 xy, yz, xz 三組曲面,各自有高斯曲率。所以曲率是一個三維 vector. Feymann lecture said there are 6 values to fully describe curvature.
  • 圓周率不是很好的曲面特徵,因為不是定值。Use N-dimension unit volume instead. S(r) = dV(r)/dr or dV(r) = S(r) dr (Gray and Vanhecke 1979).

曲面幾何的其他特徵

除了平行公理可以用來分類不同的曲面,是否有其他更直觀的特徵分類?
Yes! 可以用三角形內角和是否大於,等於,或小於 180^o \,(or \,\pi) 如下圖。

另外一個特徵是圓周率。在平面幾何任意圓的周長除以直徑永遠等於\pi .
但在曲面則會小於或大於\pi .

以上特徵(三角形內角和,圓周率)不只是定性的描述,都和高斯曲率有定量的關係。更重要的是這是曲面內稟的特性!就是在曲面的生物,就可以測量到高斯曲率。

下表總結三種幾何學:

Geometry Sphere 球面 Plane 平面/柱面 Hyperbolic 雙曲面
Creator Riemann Euclid Lobachevsky
平行公理平行線 0 1 \ge 2
三角內角和 >\pi \pi <\pi
圓周率 <\pi \pi >\pi
高斯曲率 > 0 0 < 0

上表的平行線,三角內角和,圓周率都是 global (全域) 特性。但是高斯曲率卻是 local (局部) 特性。如何連結 global and local?

  • 圓周率 Bertrand–Diguet–Puiseux theorem (Wifi 2019a). 是否有積分形式?
    K = \lim_{r\to0^+} 3  \frac{2\pi r - C(r)}{\pi r^3}  = \lim_{r\to0^+} 12 \frac{\pi r^2 - A(r)}{\pi r^4}

  • 三角形內角和 (Gauss-Bonnet theorem (Wiki 2019a)
    \sum_{i=1}^3 \theta_i = \pi + \iint K dA \quad\text{or}\quad K = \lim_{A\to0^+}\frac{\sum^3 \theta_i - \pi}{A}

  • A good article to integrate everything including Laplacian using unit volume (Gray and Vanhecke 1979).

  • 一個明顯的問題是平行公理是 global 的特性(直線永不相交),是否有直接對應局部曲率的公式? Yes,平行移動 parallel transport 就是平行的微分定義!! 和曲面微分幾何有非常深層的連結。

  • 反向問題:既然可以用三角內角和以及圓周率極限定義曲率,還需要用平行定義取率嗎? Absolutely Yes, 平行觀念和座標 grid line and local basis vector 定義吻合(local basis vector = grid line 的切線),可以直接融入純量場,向量場,張量場的微分運算。 相反三角內角和或圓周率基於極限的定義在各種微分幾何的運算並不實際。

(Keng 2018) 剛好和我想的一致,平行公理和 parallel transport.

Parallel Transport (平行移動) 的定義

(Wiki 2019a) 微分幾何對平行移動的定義:
平行移動是將流形上的幾何特性沿著光滑曲線移動的一種方法。如果流形上有一個 affine connnection(covariant derivative),那麼 affine connection 保證我們可以將流形上的向量沿著曲線移動使得它們關於這個 connection 保持「平行」。

在某種意義上說,關於 connection 的平行移動提供了將流形的局部幾何沿著曲線移動的方法:即「連接」了鄰近點的幾何。有許多種平行移動的概念,但其中一種方式等同於提供了一個 connection。事實上,通常的 connection 是平行移動的無窮小類比。反之,平行移動是 connection 的局部實現。

因為平行移動給出了 connection 的一種局部實現,它也提供了曲率的一種局部實現(holonomy)。Ambrose-Singer 定理明確了曲率與 holonomy 的關係。
\frac{D}{dx}\frac{D}{dy}V - \frac{D}{dy}\frac{D}{dx}V = R(\frac{\partial \sigma}{\partial x},\frac{\partial \sigma}{\partial y})V \quad\text{where R is curvature tensor}

簡言之,平行移動基本是一種 local connection, 以及實現 curvature 的方式。
問題是:平行移動和平行公理有什麼關係?

平行公理和 Parallel Transport 和 geodesic 三角形內角和

平行公理:一條直線和線外一點平行直線的關係。
什麼是直線?最直觀的定義就是最短路徑,在非平面的“直線”看起來並不直。
數學的表示是: \min \int ds . 使用變分法的 Euler-Lagrange equation 可以得到 Geodesics (測地線或地直線)。平面幾何的 geodesic 就是直線,推廣到微分幾何如下:

  • 一條 Geodesics (地直線)上所有點的切向量(以及和切線固定夾角的向量)都互相平行(如下圖三條地直線,都是球面上的大圓)。
  • Geodesic (地直線)線外一點和其最短距離的路徑也是 geodesic, 並且垂直夾角。
  • 其實是否直角不重要。如果不是直角但仍是 geodesic,parallel transport 就會轉該夾角。沿著同一條 geodesic 的 parallel transport 和 geodesic 的夾角不變。
  • How to prove? 如果 geodesic 構成的三角形內角和為 180, parallel transport 的夾角為 0. 大於 180, parallel transport 會產生正夾角。小於 180, parallel transport 產生負夾角。

平面和曲面幾何的 Parallel Transport (Holonomy)

平面幾何:平行移動先定義 Geodesic A, 接著線外一點 P, 可以找到一條 Geodesic B 垂直於 A at T. 再來根據平行公理,只有一條平行線 Geodesic C 平行於 A, 且垂直於 B. 因為 A 和 C 是平行線,永不相交。為了讓 parallel transport 形成封閉迴路 (holonomy), 可以在 A 上找一點 Q, 並且做垂直線到 C, 稱為 geodesic D.
此時從 A_P 的切線,一路到 A_Q 都是 parallel transport. 和 A_P 經由 B \to C \to D \to A_P 的 parallel transport vector 會是平行。也就是夾角為 0!這是一個 trivial case. 沿著長方形兩邊是切向量,兩邊是垂直向量的 holonomy, 在 close loop 的夾角為 0.

柱面幾何:直接剖面切開就是平面。parallel transport 和平面一樣。如果沿著三角形 holonomy, 因為內角和為 180, 回到起點的 parallel transport 夾角為 0. Gauss curvature 為 0.

錐面幾何:圓錐面如下圖(a)的曲率?圓錐面可以展開成平面缺角圓形如下圖(b). 除了 A 點之外任何點的曲率為 0. 因為錐面上一個 close loop, 如果不包含 A, parallel transport 夾角為 0. 因此曲率為 0.

如果包含 A 的圓,parallel transport 夾角 = \alpha , 面積 = (\pi-\alpha/2) a^2 , A 的曲率趨近無窮大。
R_A = \lim_{a\to 0}\frac{2\alpha}{a^2(2\pi-\alpha)}\to\infty

另一個問題,上圖球面非大圓 R 的 parallel transport 夾角? 常見的回答是 0. 錯誤的答案!正確的答案是用圓錐面和球面的切線為 R,再展開如上圖 b.

球面幾何:平行移動先定義 geodesic A (上圖左邊弧線), 赤道的 geodesic B 垂直於 A at T. 由平行定理我們知道球面幾何不存在和 A 平行(不相交)的 geodesic. 我們只能選擇和 B 垂直的 geodesic C, given we know A 和 C 會相交於 Q 點(其實在北極)。
此時從 A_T 的切線,一路到 A_Q 都是 parallel transport. 另一路是A_T 經由 B \to C \to A_Q 的 C 的切線(也是 parallel transport) 一定有一個夾角 \alpha . 球面幾何的三個內角和 =90+90+\alpha>180 . 這個夾角和 parallel transport 所包含的面積有關。
如果只是在局部做一個封閉迴路(holonomy),局部看起來就像歐氏空間。Parallel transport 的 holonomy 夾角為 0. Use Gauss curvature
R = angle / area \quad \text{where}\quad area = 4\pi r^2 /2 * \alpha /{2 \pi} = {\alpha r^2}
R= \frac{1}{r^2} \quad \text{independent of }\alpha

雙曲線幾何:直接考慮鞍面的三角形,內角和小於 180. parallel transport 為負夾角。

平行移動的數學表示

Parallel Transport \vec{V} 沿著一條線(不需要直線 or geodesic) \vec{r}(\lambda) 的局部變化為 \vec{0} , \lambda 是參數 normalized to 弧長,所以和座標系無關。

Parallel transport 更數學而且和座標系無關的定義:\vec{V} \vec{r}(\lambda) 切線 \vec{u} 的方向導數(微分)為 \vec{0} . (注意向量場的所有 basis 方向的導數是張量,某個特定方向, \vec{u} , 的導數是向量,也就是“affine connection". 在曲面空間稱為 Levi-Civita connection.)

\nabla_{\vec{u}} \vec{V} = \vec{u}\cdot\nabla\vec{V} = \vec{0} \quad \text{where}\quad \vec{u} = \frac{d\vec{r}}{d\lambda}

一個類比:如果不是一個向量 \vec{V} , 而是一個純量 \phi 滿足下式。(注意純量場的所有 basis 方向的導數是向量,稱為梯度 gradient. 某個特定方向, \vec{u} , 的導數是純量。)\vec{r} 就是等高線或等位線。

\nabla_{\vec{u}} \phi = \vec{u}\cdot\nabla\phi = 0 \quad \text{where}\quad \vec{u} = \frac{d\vec{r}}{d\lambda}

平行移動和等高線的數學公式雖然非常類似。問題的形式 (formulation) 卻不同。
平行移動是給定一條任意線 \vec{r}(\lambda) (open or close) 以及一點 P 的 vector \vec{V}(P) , 用 (1) 找出沿線所有點平行移動的 \vec{V} .
等高線是給定一個純量場 \phi 以及任一點 P , 用 (2) 找出通過 P 的等位線 \vec{r}(\lambda) (open or close).

張量分析用於幾何

張量(tensor) 最重要的特性就是座標系無關,可用於平面或曲面。張量分為三類:

  1. Basis vector invariant: 純量 也就是 0 階張量。
  2. Basis vector contra-variant. 為什麼 contra-variant? 因為 basis vector scale up, 張量分量會 scale down. 一階張量是 vector (n維), 二階張量 (nxn維), 三階張量(nxnxn維),依此類推。e.g. Einstein curvature tensor 是四階張量(4x4x4x4)維。 Ricci curvature tensor 是二階張量 (4×4)維。
  3. Basis vector co-variant. 也稱為 one-form. 像是 gradient, \nabla , basis vector scale up, dual basis vector scale down, 張量分量也會 scale up. 其他特性和一般張量一樣。

其他的張量運算下面展開。

先用歐氏空間的方向導數為例:

在卡氏座標系(笛卡爾直角座標系):\{\mathbf{e_1,\,e_2}\} , global constant orthonormal basis vectors.
曲線和切向量:
\vec{u} = \frac{d\vec{r}}{d\lambda} = u_1\mathbf{e_1} + u_2\mathbf{e_2}

注意:千萬不要用 [u_1, u_2] 表示法,只有在卡氏座標系可用。對於非直角或非歐氏空間一定要如上公式帶上 local basis vectors!!

先考慮 gradient operator, \nabla . 在卡氏座標系:

\nabla = \frac{\partial}{\partial e_1} \mathbf{e_1} +  \frac{\partial}{\partial e_2} \mathbf{e_2}

一個純量場(0階張量)的梯度, \nabla \phi , 是一個向量。更準確的說是一個一階張量
一個一維量場的梯度, \nabla \vec{V} , 是一個二階張量。以此類推。

純量場的方向導數,基本上是純量場的梯度 (gradient) 和方向向量的內積。

\nabla_{\vec{u}}\phi = \vec{u}\cdot\nabla\phi = (u_1\mathbf{e_1} + u_2\mathbf{e_2}) \cdot (\frac{\partial \phi}{\partial e_1}\mathbf{e_1}+\frac{\partial \phi}{\partial e_2}\mathbf{e_2})\\ = u_1 \frac{\partial \phi}{\partial e_1} + u_2 \frac{\partial \phi}{\partial e_2} \quad\text{where}\quad\mathbf{e_i}\cdot \mathbf{e_j}=\delta_{ij} \quad(1)\,\text{!wrong for non-Euclidean!}

\nabla_{\vec{u}}\phi =  \vec{u}\cdot\nabla\phi = 0 代表 \vec{u} 和梯度 \nabla \phi 垂直。因此 \vec{u} 形成 \phi 的等位線。

(1) 只能用於卡氏座標系。接下來幾個重點:1.座標系無關;2.推廣到曲面;3.簡化公式
現在引入 covariant/contravariant basis and tensor, 和愛因斯坦 summation notation.

Basis vectors \{\mathbf{e_1,\,e_2}\} , 配合 contra-variant 分量。
Dual basis vectors \{\mathbf{e^1,\,e^2}\} , 配合 co-variant 分量。

\mathbf{e_i}\cdot\mathbf{e^j}=\mathbf{e^i}\cdot\mathbf{e_j}=\delta_i^j=\delta_j^i \qquad(2)

注意 \mathbf{e_i}\cdot \mathbf{e_j}\ne\delta_{ij} , 有兩個原因。如果 \mathbf{e_i} 放大兩倍,\mathbf{e_i}\cdot \mathbf{e_i}=4 , 就是座標系無法 scale. 另外非直角座標系,\mathbf{e_1}\cdot \mathbf{e_2}\ne 0 , 就是座標系無法非直角。正確的結論如下,可以用於非歐非直角座標系。g_{ij} 是二階張量 (metric tensor) 的分量。
\mathbf{e_i}\cdot\mathbf{e_j}=g_{ij} \qquad \mathbf{e^i}\cdot\mathbf{e^j}=g^{ij} \qquad(3)

對於歐氏空間非直角座標系,(1) 可以修改為:
\nabla_{\vec{u}}\phi = \vec{u}\cdot\nabla\phi = (u^1\mathbf{e_1} + u^2\mathbf{e_2}) \cdot (\frac{\partial \phi}{\partial e_1}\mathbf{e^1}+\frac{\partial \phi}{\partial e_2}\mathbf{e^2})\\ = u^1 \frac{\partial \phi}{\partial e_1} + u^2 \frac{\partial \phi}{\partial e_2} =  \quad\text{where}\quad\mathbf{e^i}\cdot \mathbf{e_j}=\delta^i_j \qquad(4)

使用愛因斯坦 notation 重寫 (4)
\nabla_{\vec{u}}\phi = \vec{u}\cdot\nabla\phi = u^i\mathbf{e_i} \cdot\partial_j \phi \mathbf{e^j} = u^i \partial_i \phi\qquad(5)

曲面空間的方向導數:

曲面空間的純量場方向導數和 (5) 完全相同。因為純量場是 basis invariant.

但是曲面空間的向量場是 basis covariant or contravariant, 所以曲面的向量場的方向導數會比歐氏空間複雜的多。

曲面空間向量場 \vec{V} 的方向 \vec{u} 導數

  • 似乎可視每個分量為純量場,求每個純量場的方向導數。Wrong! 原因是 basis vector 在非歐空間不是 constant vector, 甚至歐氏空間的極座標也不是 constant vector. 因此 basis vector 的微分會產生新的 component! (Christopher symbol, or connection).
  • 向量場的梯度是二階張量,直接求梯度太複雜。比較好的作法是分解方向導數的“方向向量”為 basis vector linear combination, 最後再結合為真正方向導數。

$latex \vec{V} = V^1 \mathbf{e_1} + V^2 \mathbf{e_2}
= V^{\alpha} \mathbf{e_\alpha} \\
\vec{u} = u^1 \mathbf{e_1} + u^2 \mathbf{e_2} $

$latex \nabla_{\vec{u}}\vec{V} = \vec{u}\cdot\nabla\vec{V}
= (u^1 \mathbf{e_1} + u^2 \mathbf{e_2})\cdot\nabla\vec{V} \\
\nabla_{e_\beta}\vec{V}=\mathbf{e_\beta} \cdot\nabla\vec{V}
= \mathbf{e_\beta} \cdot\nabla (V^\alpha \mathbf{e_\alpha})
= \mathbf{e_\beta} \cdot \frac{\partial(V^\alpha \mathbf{e_\alpha})}{\partial e_i} \mathbf{e^i} \\
= \frac{\partial(V^\alpha \mathbf{e_\alpha})}{\partial e_\beta}
= \frac{\partial V^\alpha}{\partial e_\beta}\mathbf{e_\alpha} +
V^\alpha \frac{\partial \mathbf{e_\alpha}}{\partial e_\beta}
= \frac{\partial V^\alpha}{\partial e_\beta}\mathbf{e_\alpha} +
V^\alpha \Gamma^k_{\alpha\beta}\mathbf{e_k} \\
= \frac{\partial V^\alpha}{\partial e_\beta}\mathbf{e_\alpha} +
V^i \Gamma^{\alpha}_{i\beta}\mathbf{e_\alpha}
= (\frac{\partial V^\alpha}{\partial e_\beta} +
V^i \Gamma^{\alpha}_{i\beta})\mathbf{e_\alpha} $

此處利用 Christopher symbol 以及 \alpha\to i \,,\, k\to\alpha
$latex \frac{\partial \mathbf{e_\alpha}}{\partial e_\beta} =
\Gamma^k_{\alpha\beta}\mathbf{e_k} $
如果是歐氏空間且直角座標系,Christopher symbol \Gamma^a_{bc}=0
如果是歐氏空間但極座標系,因為每一點的 local basis vector 方向都不同,Christopher symbol 不為 0, 一共有 2x2x2=8 個 components [@cyrilChristoffelSymbol2016].

\Gamma_{i j}^{r}=\left(\begin{array}{cc}{\Gamma^r_{\theta\theta}=-r} & {0} \\ {0} & {0}\end{array}\right) \qquad \Gamma_{i j}^{\theta}=\left(\begin{array}{cc}{0} & {\frac{1}{r}} \\ {\frac{1}{r}} & {0}\end{array}\right)

愛因斯坦 notation and convention: 向量場的梯度微分是二階張量(樓上加樓下)。
\nabla_{e_\beta}\vec{V}= \nabla_\beta V^\alpha  = \partial_\beta V^\alpha+\Gamma^{\alpha}_{i\beta} V^i = \partial_\mu V^\nu+\Gamma^{\nu}_{\lambda\mu} V^\lambda

最後 $latex \nabla_{\vec{u}}\vec{V} = \vec{u}\cdot\nabla\vec{V}
= (u^1 \mathbf{e_1} + u^2 \mathbf{e_2})\cdot\nabla\vec{V}\\
= u^\beta (\frac{\partial V^\alpha}{\partial e_\beta} +
V^i \Gamma^{\alpha}_{i\beta})\mathbf{e_\alpha}
= (u^i \frac{\partial V^k}{\partial e_i} +
u^i V^j \Gamma^{k}_{ij})\mathbf{e_k}
= (u^i {\partial_i V^k} +
u^i V^j \Gamma^{k}_{ij})\mathbf{e_k} $

愛因斯坦 notation and convention: 向量場的方向導數是一階張量(樓上)。
\nabla_{\vec{u}}\vec{V} = u^i {\partial_i V^k} + u^i V^j \Gamma^{k}_{ij}

整理愛因斯坦 notation 的原則:

  • \vec{u}=u^i\mathbf{e_i}=u^i , 只有一個樓上 index,代表一階張量 contra-variant.
  • \nabla\phi = \partial_i \phi \mathbf{e^i} = \partial_i \phi , 只有一個樓下 index,代表一階張量 co-variant (one-form).
  • \vec{u}\cdot\nabla\phi = u^i \partial_i \phi , 張量和 1-form 張量內積結合,同一個 index i 樓上樓下抵銷,變成 0 階張量(純量)。
  • 張量微分(gradient) 階數+1\nabla (0 階張量)得到一階 one form 張量;\nabla (一階張量)得到二階張量;以此類推。\nabla_\beta V^\alpha  = \partial_\mu V^\nu+\Gamma^{\nu}_{\lambda\mu} V^\lambda 是二階(樓上加樓下)tensor.
  • 張量和 one-form 張量內積 階數-1,index 樓上樓下抵銷。
  • 張量場的方向導數是先微分張量(+1)場再和方向向量內積(-1),階數不變,e.g. 純量場方向導數是純量。向量場的方向導數是一階張量:\nabla_{\vec{u}}\vec{V} = u^i {\partial_i V^k} + u^i V^j \Gamma^{k}_{ij} 是一階(樓上) 張量。
  • 因為張量場的方向導數階數不變。因此可以重覆這個運算。最常見是二次方向導數(沿不同的方向,如 basis vectors)。除了在平面歐氏座標系以外,一般是不能交換。\nabla_{\vec{w}}\nabla_{\vec{u}}\vec{V} \ne \nabla_{\vec{u}}\nabla_{\vec{w}}\vec{V} . 事實上,\nabla_{\vec{w}}\nabla_{\vec{u}}\vec{V} - \nabla_{\vec{u}}\nabla_{\vec{w}}\vec{V} = R\vec{V} ? (TBC)

回到曲線的方向導數=絕對導數(Absolute Derivative:等高線,地直線,平行移動)

曲線和切向量回顧:
\vec{u} = \frac{d\vec{r}}{d\lambda} = u_1\mathbf{e_1} + u_2\mathbf{e_2}
純量場的線方向導數
\nabla_{\vec{u}}\phi(\vec{r}) = \frac{d\vec{r}}{d\lambda}\cdot\nabla\phi(\vec{r}) = \frac{d\phi(\vec{r})}{d\lambda}

其實這就是著名的 gradient theorem, 或是線積分基本定理推廣到曲面。重點是曲面線積分是路徑無關!或是封閉迴路線積分為 0. 但在向量場不是如此。
$latex \int_p^q \nabla_{\vec{u}}\phi(\vec{r}) d\lambda
= \int_p^q \nabla\phi(\vec{r})\cdot d\vec{r} = \phi(q) – \phi(p) $

以純量場為例

Let \phi(\vec{r}) = |\vec{r}|^2 座標系無關:
射線運動卡氏座標系:
\phi = x^2 + y^2 and \nabla \phi = 2x \mathbf{e_x} + 2y \mathbf{e_y}
\vec{r} = c_0\lambda\mathbf{e_x} + c_1\lambda\mathbf{e_y} , where c_0^2+c_1^2=1
\frac{d\vec{r}}{d\lambda} = c_0\mathbf{e_x} + c_1\mathbf{e_y}

\nabla_{\vec{u}}\phi(\vec{r}) = \frac{d\vec{r}}{d\lambda}\cdot\nabla\phi(\vec{r}) = 2x c_0 + 2y c_1= 2\lambda c_0^2 + 2\lambda c_1^2 = 2\lambda

\phi(\vec{r}) = \phi(\lambda) = c_0^2 \lambda^2 + c_1^2 \lambda^2 = \lambda^2
\frac{d\phi(\vec{r})}{d\lambda} = 2\lambda

射線運動極座標系:
\phi = r^2 and \nabla \phi = 2r \mathbf{e_r}
\vec{r}= \lambda \mathbf{e^r} + \theta_o \mathbf{e^\theta}
\frac{d\vec{r}}{d\lambda} = \mathbf{e_r}

\nabla_{\vec{u}}\phi(\vec{r}) = \frac{d\vec{r}}{d\lambda}\cdot\nabla\phi(\vec{r}) = 2r = 2\lambda

\phi(\vec{r}) = \phi(\lambda) = \lambda^2
\frac{d\phi(\vec{r})}{d\lambda} = 2\lambda

以向量場為例

向量場的線方向導數,稱為絕對導數(absolute derivative)
$latex \nabla_{\vec{u}}\vec{V}(\vec{r}) = \frac{d\vec{r}}{d\lambda}\cdot\nabla\vec{V}(\vec{r}) \\
= u^i {\partial_i V^k} + u^i V^j \Gamma^{k}_{ij}
= \frac{d V^k}{d\lambda} + \frac{d x^i}{d\lambda} V^j \Gamma^{k}_{ij} $

Parallel transport

\vec{u} = \frac{d\vec{r}}{d\lambda} or u^k = \frac{dx^k}{d\lambda} , parallel transport 就是向量場 \vec{V} 沿 \vec{r} 方向導數為 \vec{0} ,上式變成:
\frac{d V^k}{d\lambda} + \frac{dx^i}{d\lambda} V^j \Gamma^{k}_{ij}=0

積分形式: V^k(q) - V^k(p) = -\int_p^q V^j\Gamma^{k}_{ij}{dx^i} . 注意 p 點和 q 點的 basis vectors 不一定相同。即使 V^k(p) \ne V^k(q) 也不代表 \vec{V}(p) \ne \vec{V}(q) . 下面極坐標例子可以看出。
另一個重點是 V^k(q) 有很多 dx^i 可以到達,結果不一定相等!(path dependent)

卡氏座標系:\Gamma=0 , V^k = c^k , constant as expected.
極座標系:k =\{r \,, \theta \} ; ds^2 = dr^2 + r^2 d\theta^2 .

|\mathbf{e^r}| = \frac{\partial s}{\partial r} =1
|\mathbf{e^\theta}| = \frac{\partial s}{\partial \theta} =r

\vec{x} = x^r \mathbf{e^r} + x^\theta \mathbf{e^{\theta}}
\frac{d V^r}{d\lambda} + \frac{dx^i}{d\lambda} V^j \Gamma^{r}_{ij} = \frac{d V^r}{d\lambda}-\frac{dx^\theta}{d\lambda}V^\theta x^r= 0

\frac{d V^\theta}{d\lambda} + \frac{dx^i}{d\lambda} V^j \Gamma^{\theta}_{ij} = \frac{d V^\theta}{d\lambda} + \frac{dx^r}{d\lambda} \frac{V^\theta}{x^r} + \frac{dx^\theta}{d\lambda} \frac{V^r}{x^r} = 0

Example 1: 射線運動 (起點在r=1上)
Let x^r = \lambda+1\,,\, x^\theta = c . \vec{V}(\lambda=0) = r_o \mathbf{e^r} + \theta_o \mathbf{e^\theta}
\frac{d V^r}{d\lambda} = 0

\frac{d V^\theta}{d\lambda} + \frac{V^\theta}{\lambda+1} = 0

V^r = r_o \,,\, V^\theta = \frac{\theta_0}{\lambda+1} = \frac{\theta_0}{x^r} \text{ Yes, scale down the basis vector } e^\theta

Example 2: 圓周運動
Let x^r = 1\,,\, x^\theta = \lambda . 弧長等於夾角! \vec{V}(\lambda=0) = r_o \mathbf{e^r}
\frac{d V^r}{d\lambda}-V^\theta = 0

\frac{d V^\theta}{d\lambda} + V^r = 0

V^r = r_o\cos\lambda \,,\, V^\theta = -r_o\sin\lambda as expected below.

平行移動和曲率的關係

Geodesics 是 parallel transport 的特例

\vec{V} = \vec{u} = \frac{d\vec{r}}{d\lambda} or V^k = u^k = \frac{dx^k}{d\lambda} , 上式變成:
\nabla_{\vec{u}}\vec{u}(\vec{r}) = 0

\frac{d^2 x^k}{d\lambda^2} + \Gamma^{k}_{ij}\frac{dx^i}{d\lambda} \frac{dx^j}{d\lambda}=0

Reference

Gray, A., and L. Vanhecke. 1979. “Riemannian Geometry as Determined by
the Volumes of Small Geodesic Balls.” Acta Mathematica 142: 157–98.
https://doi.org/10.1007/BF02395060.

Keng, Brian. 2018. “Hyperbolic Geometry and Poincaré Embeddings.”
Bounded Rationality. June 17, 2018.
http://bjlkeng.github.io/posts/hyperbolic-geometry-and-poincare-embeddings/.

Wiki. 2019a. “Gaussian Curvature.” Wikipedia.
https://en.wikipedia.org/w/index.php?title=Gaussian_curvature&oldid=910677726.

———. 2019b. “Parallel Transport.” Wikipedia.
https://en.wikipedia.org/w/index.php?title=Parallel_transport&oldid=911069387.

———. 2019c. “Theorema Egregium.” Wikipedia.
https://en.wikipedia.org/w/index.php?title=Theorema_Egregium&oldid=895212511.


  1. \kappa = \frac{y''}{(1+y'^2)^{3/2}}  

相對論和張量分析 – Coordinate Covariant, Contravariant, Invariant (座標系協變,逆變,不變)

by allenlu2007

黎曼微分幾何是廣義相對論的數學基礎。張量分析是黎曼微分幾何的必備數學工具。

  • 黎曼幾何是彎曲幾何(curved geometry),不同於歐式平直幾何(straight geometry)。
  • 黎曼幾何研究內稟(intrinsic)的特性以及嵌入(embedded)歐式幾何的特性。例如 2D Swiss Roll 的內稟幾何特性和一張平紙無異 (i.e. intrinsic curvature = 0), 但是 embedded 在高維(3D)歐式空間卻有曲率。

由於相對論的影響,不同(運動)觀察者所量測的物理定律 (e.g. least action principle) 都要一樣。這剛好對應幾何不同座標系“不變性”。不變的 scalar 可以使用 invariant. 不變的 vector/tensor, 因為分量仍然和座標系相關,我們用 covariant 表述。

觀察者和座標系容易混淆的重點:

  1. 不同(運動)觀察者對應幾何不同座標系. 不要執著一個是動態,一個是靜態。這裏所說的幾何座標系都包含時間軸 (t -axis), 從下到上代表過去到未來。任何運動都可以用一條靜態的世界線(worldline) 表示。

  2. 一條 wordline 是相對於一個觀察者(座標系)。例如 A 觀察到一個靜止物體 S,在 A 的(直角)座標系 S 的位置就是一條垂直的 world line. B 是一個相對 A 定速 v 移動的觀察者。B 要如何得到 S 的 world line?

    有兩個方法:

    (i) B 創造自己的(直角)座標系以及 S 的 world line. 顯然在 B 的直角座標系 S world line 是一條斜線。斜率是 -v (assuming v \ll c ). 兩個不同的直角座標系 for S 運動,對應 A 和 B 觀察者 (observer oriented)。這不是我們所要的。如果有無窮多觀察者(等速或等加速觀察者)就需要無窮多直角座標系以及無窮多的 world lines for the same S 的運動。
    我們關心 S 的運動 (object oriented) 一次看到 A 和 B 座標系的詮釋,見 (ii). 這對廣義相對論研究時空內稟幾何特性很重要。

    (ii) 在 A 的直角座標系 (grid line) 重劃 B 的(非直角)座標系 (grid line)。同一條垂直的 world line, 在 A 座標系 grid line 是靜止;但在 B 座標系 grid line 則是運動。

  3. 下圖是一個例子。A 觀察者是直角座標系 (x, t) , 網格線 (grid line) 是正方形。但對 B 觀察者(相對速度為 v )對應座標系 (x', t') , 網格線是平行四邊形。此處是用 Lorentz 變換,所以 t'=0,1,2 格線是斜線。如果是 Galilean 變換t'=0,1,2 會是水平線。對於 A 的靜止物體 S, 其 world line 就是 x=0 垂直線。對 B 座標系 (x',t') 而言,x=0 垂直線 world line 對應一個後退的運動。

  4. 注意這種重劃座標系網格並不限於相對速度固定的觀察者。原則上適用任何運動觀察者。假設觀察者 C 是相對于 A 有固定等加速度。也可以定義 (x'', t'') 網格線。下圖是 Lorentz 變換 C 的座標系。 x''=const 變成雙曲線。t''=const 變成聚集在 x''=-1 的投射線 (假設光速=1). 網格線變成更奇怪的四邊形。注意在 Lorentz transform 所有速度無法超越光速。觀察者 C 無法持續加速超越光速。 Galilean 變換 C 的座標系則不同。x''=const 變成拋物線。t''=const 還是水平線。

回到物理和幾何的對照:

  • 運動的物體 \iff 幾何的一條世界線 C (不是面或其他形狀)
  • 不同慣性運動或是非慣性運動觀察者 \iff 幾何的不同座標系以及網格
  • 狹義相對論:物理定律在不同慣性座標系不變 (scalar:invariant; vector/tensor:covariant)
  • 廣義相對論:物理定律考慮等效原理在不同座標系不變
  • 物理定律 \iff 幾何不變量 (如何對應 scalar: invariant; tensor: covariant)

    • 微分幾何的 1st fundmental form ds^2 = c^2 dt^2-dx^2-dy^2-dz^2 = c^2 dt'^2-dx'^2-dy'^2-dz'^2 幾何量是微分弧長 coordinate invariant 對應物理的 Lorentz transform. 如果少了 c^2 dt^2 就是 (-1)*Galilean transform.

    • 注意微分幾何的 1st fundamental form 可以推廣到 curved space or coordinate. 這裏使用愛因斯坦 summation notation. ds^2 = g_{\alpha\beta}{dx^\alpha}{dx^\beta}

    • 牛頓第二定律: \vec{F} = m\vec{a} \, 在 Galilean transform 不變。因為是 vector, vector components 在不同座標系不同,稱為 covariant. 對應的幾何不變量是“curvature vector". 但在 Lorentz transform, \vec{F} \,and\, \vec{a} 甚至不是 4-vector/tensor, 連 covariant 也不是。如何得到符合 Lorentz tranform 4-vector/tensor 的力學定律呢?下面描述。

    • 最小作用力原理 (Least action principle) 是比牛頓運動定律更基本的定理。分析力學首先找到 Lagrangian L (單位是能量),differential action dS = L dt (單位是能量*時間,大自然兩者都要省)。Total action 就是空間兩點 (p_1, p_2) 和時間兩點 (t_1, t_2) 路徑積分。最小作用原理用變分法找到最小 S 的路徑。 S = \int_{t_1}^{t_2} L dt

    • 舉一個例子:古典力學 Lagrangian of free particle (no potential field) in Cartesian and polar coordinate. [@ActionPhysics2019]

      $latex L =\frac{1}{2} m v^2=\frac{1}{2} m (\dot{x}^2 + \dot{y}^2)
      \quad L =\frac{1}{2} m (\dot{r}^2 + r^2 \dot{\varphi}^2) $

    • 最小作用原理推廣到 4D 時空 for two events (x_1, x_2, x_3, x_4) and (x'_1, x'_2, x'_3, x'_4) . 因為時間 t 已經變成 x_1 , 需要用 dummy variable \lambda . 上述的 S 公式完全適用。

    • 但在相對論力學 4D 時空的 Lagrangian of free particle (no potential field). [@RelativisticLagrangian2019] L = -\frac{m_o c^2}{c} \sqrt{g_{\alpha\beta}\frac{dx^\alpha} {d\lambda}\frac{dx^\beta}{d\lambda}} = -\frac{m_o c^2}{c} \frac{ds}{d\lambda} \\ \approx m_o c^2 + \frac{1}{2} m_o v^2! \quad \text{assuming } v \ll c g_{\alpha\beta} 沒有單位。根號的單位是速度,和光速 c 單位抵銷。m_o c^2 的單位是能量。所以 L 的單位在 4D 時空仍然是能量。\lambda 只是 dummy variable (時間單位)。ds 是 1st fundmental form 的根號,或是微分弧長。滿足 Lorentz transform invariant. In a word, L 對應幾何量是弧速。 S = \int L d\lambda = \int -\frac{m_o c^2}{c} \sqrt{g_{\alpha\beta}\frac{dx^\alpha} {d\lambda}\frac{dx^\beta}{d\lambda}} d\lambda  重點是 S (action) 對應的是 4D 時空兩點之間的總弧長 (with a proportional constant). 最小作用原理對應幾何量是兩點之間最短距離的弧線,稱為 geodesic (測地線或地直線) which is coordinate invariant. [@cyrilGeodesicEquation2018] 也就是說,不同慣性運動的觀察者的最小作用原理都是一樣。但如何推廣到非慣性觀察者?

UTM. 2018. “Proper Time, Coordinate Systems, Lorentz Transformations | Internet Encyclopedia of Philosophy.” 2018. https://www.iep.utm.edu/proper-t/.

Coordinate Invariant, Covariant, Contravariant

Coordinate invariant 基本是 scalar quantity or field, 例如 ds , geodesic, scalar curvature, etc.

Contravariance 基本是 vector/tensor, 例如 curvature tensor, …. 什麼是 contravariant? (A^x e_x) base vector (e_x) 變大,component (A^x) 變小。所以是 contra-variance.

Covariance 也是 vector/tensor, 為什麼是 covariance? 一個例子是 gradient of scalar field. (\frac{\partial f}{\partial x} e^x, \frac{\partial f}{\partial y} e^y, \frac{\partial f}{\partial z} e^z) , base vector e^x 越大, component 也越大。

Covariance and Contra-variance of Vector/Tensor

歐式空間(或是 tangent space on manifold)

如何在座標系表達一個向量 \vec{v} ?
選擇一組 basis vectors (假設不正交). 然後分解成 basis vectors 的平行分量 (下圖右 typo v_x \to v^x \,\,and\,\, v_y \to v^y )。
\vec{v} = v^x \mathbf{e_x}+v^y \mathbf{e_y}
我們稱 v^x and v^y contra-variance, 因為 v^x \,,v^y \mathbf{e_x} \,,\mathbf{e_y} 反比的關係。
問題是如何用數學式計算 v^x and v^y ? (見下文)
如果是在直角座標系如下,可以看出 contra-variant!

v^x = \frac{\vec{v} \cdot e_x}{\lVert\mathbf{e_x}\rVert^2} \quad  v^y = \frac{\vec{v} \cdot e_y}{\lVert\mathbf{e_y}\rVert^2}

顯然在非直角座標系是錯的!正確的數學表述在後面。

此時稱 \mathbf{e_x} \,,\mathbf{e_y} (tangent) basis vectors 可以使用 Einstein summation notation:
\vec{v} = v^x \mathbf{e_x}+v^y \mathbf{e_y}= v^i \mathbf{e_i}\quad \text{or } v^i \text{ in short}

下一個問題,是否還有其他方法表示一個向量?Yes!
我們可以用垂直投影分量表示一個向量 (下圖左, v_x = \vec{v} \cdot e_x \,\, v_y = \vec{v} \cdot e_y ).
我們稱 v_x and v_y covariance, 因為 v_x \,,v_y \mathbf{e_x} \,,\mathbf{e_y} 正比的關係。
\vec{v} = v_x \mathbf{e^x}+v_y \mathbf{e^y}= v_i \mathbf{e^i} \quad \text{or } v_i \text{ in short}

問題是如何找到對應的 (cotangent) basis vector e^x and e^y ?

為了滿足 v_x = \vec{v} \cdot e_x v_y = \vec{v} \cdot e_y ,可以得出 e^i e_j = \delta^i_j !
也就是說 e^x \perp e_y \,, e^y \perp e_x . 同時 normalize e^x e_x = e^y e_y = 1 . 也就是反矩陣關係!

注意 e^x \nparallel e_x \,, e^y \nparallel e_y 除非 e_x \,, e_y 是直角座標系。

此時可以用新的 (cotangent) basis vector 表示:
\vec{v} = v_x \mathbf{e^x}+v_y \mathbf{e^y}= v_i \mathbf{e^i} = v^i \mathbf{e_i}

v_x = \vec{v} \cdot e_x \quad v_y = \vec{v} \cdot e_y \quad or \quad v_i = \vec{v} \cdot e_i

v^x = \vec{v} \cdot e^x \quad v^y = \vec{v} \cdot e^y \quad or \quad v^i = \vec{v} \cdot e^i

因為 e^i e_i 是反比。e_i 越大,e^i 越小,v^i 就愈小,稱為 contra-variance. 但是 v_i 卻是愈大,和 e_i 正比,稱為 covariance.

Inner Product Using Einstein Notation
\vec{u} \cdot \vec{v} = u^i \mathbf{e_i} v_j \mathbf{e^j} = u^i v_j \delta^i_j = u^i v_i = u_i v^i \\ = u^i \mathbf{e_i} v^j \mathbf{e_j} = u^i v^j g_{ij} = u_i v_j g^{ij}

Scalar field is invariant.
Vector/tensor field is contravariant.
Gradient of scalar field is a vector field, which is co-variant.

https://www.quora.com/What-are-%E2%80%9Ccovariant%E2%80%9D-and-%E2%80%9Ccontravariant%E2%80%9D-vectors-as-intuitive-representations

為了直觀了解 curved coordinate and space, 用以下的順序理解:直角座標系\to 斜角座標系 \to 歐式空間極座標系 (curved coordinate) \to 球面座標系 (curved space)

  1. orthonormal 座標系
    e^1 = (1, 0) \,\,  e^2= (0, 1)
    e_1 = (1, 0) \,\,  e_2= (0, 1)

  2. scaled 座標系
    e^1 = (a, 0) \,\,  e^2= (0, b)
    e_1 = (\frac{1}{a}, 0) \,\,  e_2= (0, \frac{1}{b})

  3. Non-orthogonal 座標系
    e^1 = (a, 0) \,\,  e^2= (d, b)
    e_1 = (\frac{1}{a}, \frac{-d}{ab}) \,\,  e_2= (0, \frac{1}{b})

  4. Polar 座標系 (r,\theta) \to (dr, r d\theta)
    e^1 = (a, 0) \,\,  e^2= (0, b)
    e_1 = (\frac{1}{a}, 0) \,\,  e_2= (0, \frac{1}{b})

三維空間
推廣上述二維的垂直觀念到三維空間。
e^1 的方向是由 e_2 \,, e_3 所形成平面的法線 (normal) 決定,i.e. e^1 e_2 = e^1 e_3 = 0 .

e^1 的大小是由 e_1 決定,i.e. e^1 e_1 = 1

同理 e^i 的方向是是由 e_j \,, e_k 所形成平面的法線 (normal) 決定。e^i e_j = \delta^i_j

推廣到 n 維空間:
Tangent basis: e1, e2, …, en, C(n,1)=n
Reciprocal basis: e1.,… C(n, n-1) = n

Curved Space or Coordinate System
Basis Vector (e_1, e_2, ...)
Reciprocal Basis Vector (e^1, e^2, ...)
e_i e_j = g_{ij} metric tensor g matrix
e^i e^j = g^{ij} inverse metric tensor, 真的是 inverse g matrix.
e^i = g^{ij} e_j \quad\quad e^i e_k = g^{ij}e_j e_k = g^{ij}g_{jk} = \delta^i_j

Covariant Derivative of Scalar Field and Vector Field

Covariant derivative 是 directional derivative 的推廣。同樣用以下的順序理解:直角座標系\to 斜角座標系 \to 歐式空間極座標系 (curved coordinate) \to 球面座標系 (curved space)

  1. Scalar field: f(x, y) on orthonormal 座標系
    e^1 = (1, 0) \,\,  e^2= (0, 1)
    e_1 = (1, 0) \,\,  e_2= (0, 1)
    \nabla f = \frac{\partial f}{\partial x} e^x + \frac{\partial f}{\partial y} e^y

  2. Scalar field: f(x, y) on scaled 座標系
    e^1 = (a, 0) \,\,  e^2= (0, b)
    e_1 = (\frac{1}{a}, 0) \,\,  e_2= (0, \frac{1}{b})
    \nabla f = \frac{\partial f}{\partial x} e^x + \frac{\partial f}{\partial y} e^y

  3. Scalar field: f(r, \theta) on 極座標系
    e_i\cdot e_j = g_{ij} , e^i\cdot e^j = g^{ij} , e_i\cdot e^j = \delta_{ij} .
    e_i = g_{ij} e^j , e^i = g^{ij} e_j
    ds^2 = dr^2 + r^2 d\theta^2 = g_{ij}dx^idx^j
    g_{ij}=\left(\begin{array}{cc}{g_{rr}=1} & {g_{r\theta}=0} \\ {g_{\theta r}=0} & {g_{\theta\theta}=r^2}\end{array}\right) \,\, g^{ij}=\left(\begin{array}{cc}{g^{rr}=1} & {g^{r\theta}=0} \\ {g^{\theta r}=0} & {g^{\theta\theta}=\frac{1}{r^2}}\end{array}\right)
    \det(g_{ij})=|g|= r^2 , \det(g^{ij})=|g|^{-1}=1/r^2
    $latex \nabla f(r,\theta) = \frac{\partial f}{\partial r} e^r + \frac{\partial f}{\partial \theta} e^\theta \quad\text{(1- form, covariant)}\\
    = \frac{\partial f}{\partial r} e_r + \frac{1}{r^2}\frac{\partial f}{\partial \theta} e_\theta \qquad\text{(contra-variant)} $

    另外常見的表示法是 Euclidean metric, normalized form, |g|=1 :
    \nabla f(r,\theta) = \frac{\partial f}{\partial r} \hat{e}_r + \frac{1}{r}\frac{\partial f}{\partial \theta} \hat{e}_\theta

一個純量場(0階張量)的梯度, \nabla \phi , 是一個向量。更準確的說是一個一階張量(1-form)。在求方向導數或是散度都很方便。
一個一維量場的梯度, \nabla \phi , 是一個二階張量(index 在樓上和樓下)。以此類推。

重新回到物理和幾何的對照

  1. (三維空間 + 一維時間): 滿足 Galilean tranform. F = ma , 最小作用原理~?。
  2. 四維時空:滿足 Lorentz tranform. 最小作用原理 \iff Geodesic.
  3. (n-維空間 + 一維時間):多粒子運動。constrain 變為 manifold. Tangent space (contra-variant basis): Lagrangian mechanics. Cotangent space (co-variant basis): Hamiltonian mechanics. How to expand this?

Reference

Wiki. 2019. “Covariance and Contravariance of Vectors.” In Wikipedia.
https://en.wikipedia.org/w/index.php?title=Covariance\_and\_contravariance\_of\_vectors&oldid=904269873

“Action (Physics).” 2019. Wikipedia.
https://en.wikipedia.org/w/index.php?title=Action_(physics)&oldid=904671703

Cyril. 2018. “Geodesic Equation from the Principle of Least Action.” September 18, 2018.
http://einsteinrelativelyeasy.com/index.php/general-relativity/97-geodesic-equation-from-the-principle-of-least-action.

“Relativistic Lagrangian Mechanics.” 2019. Wikipedia.
https://en.wikipedia.org/w/index.php?title=Relativistic_Lagrangian_mechanics&oldid=884947845

Appendix

Example 1: 牛頓空間以及不同慣性運動觀察者

  • 3D 歐式空間 + 1D 時間 (空間時間各自獨立,互不影響)
  • 不同慣性觀察者:
    • x' = x - v t
    • t' = t \,\,;\,\, y' = y \,\,;\,\, z' = z
    • (x,y,z,t) 座標系。為了簡化在 2D plan,只考慮(x, t) . 垂直軸是 t 軸。
    • 靜止在原點的物體。它的 world line 是一條垂直線
  • 對應座標轉換:x' = constant grid line 變成斜線,斜率由速度決定
  • 物理定律不變 (invariant):F = m a = m \ddot{x} = m \ddot{x'}

Example 2: 狹義相對論慣性運動觀察者

  • 空間:4D 閔式幾何 (Minkovsky space-time 4D space)
  • 運動方程式:least action principle <-> geodesic
  • 不同慣性觀察者: y = x – v*t

Example 2: 廣義相對論

愛因斯坦“時空場”方程式

by allenlu2007

愛因斯坦受了 Michelson 實驗光速在“以太”不同相對速度都不變,以及 Lorentz 對此提出的 Lorentz 變換 (相對於 Gelileo 變換) 的影響。提出狹義相對論的兩大定理(假設),而推導出狹義相對論。

  1. 所有物理定律和慣性座標系無關 (covariant principle),不同慣性座標系對應不同慣性相對運動的觀察者。
  2. 光速在所有慣性座標系是不變且是極速 (invariant)

定理一稱為(狹義)相對性原理。意即所有慣性座標系(系統)都是相對而且等價。沒有哪一個座標系比較獨特。一般靜止的座標系只是相對地球是靜止。對於月球,太陽,或是宇宙中的觀察者,地球上靜止的座標系事實上是近似等速直線運動。(狹義)相對性原理在 Galileo 和牛頓就提出。牛頓力學本身就符合(狹義)相對性原理。

愛因斯坦真正的突破是提出定理二。時間和空間因而連結成為四維時空。推導出一系列令人驚訝的結果:時間膨脹,長度收縮,質能互換。愛因斯坦的老師 Minkovsky 引入閔式空間: (x_0, x_2, x_3, x_4) 替代 (ict, x, y, z) . 一開始愛因斯坦並沒有 appreciate 閔式空間。閔式時空 (space-time) 引入微分幾何的觀念,為後來的廣義相對論開路。閔式空間和歐式空間一樣是平直空間。其微分第一基本式 (1st fundmental form)
ds^2 = -dx_0^2+dx_1^2+dx_2^2+dx_3^2 (平直空間 cross terms = 0)
也就是時空(ds )座標不變 (invariant)! 所以時間膨脹和長度收縮就是保持 ds 不變。相對於歐式空間時間和空間是互相獨立。

後來廣義相對論再把閔式空間推向黎曼空間,從平直空間推廣到彎曲空間。

廣義相對論的定理

愛因斯坦深受馬赫所有運動都是相對的影響,對於狹義相對論的定理只適用於慣性座標系,也就是慣性座標系的特殊地位並不滿意。由此提出廣義相對論的三個定理。

  1. 所有物理定律和(慣性和非慣性)四維時空座標系無關 (invariant->covariant?)
  2. 光速在所有(慣性和非慣性)四維時空座標系是不變且是極速 (invariant)
  3. 非慣性或加速座標系觀察的慣性力和慣性座標系的引力等效(等效原理)

廣義相對論的重點:

  1. 移除慣性座標系的特殊地位。
  2. 牛頓的引力定律 (F = \frac{GMm}{r^2}) 不滿足四維時空座標系不變,需要修正。
  3. 牛頓的引力是超距力違背光速極速。(=>需要場論)

一個想法是導入引力場滿足 Lorentz 變換,類似電磁場方程式。引力波的速度等於或小於光速。如果只是這樣,愛因斯坦就無法名垂千古。愛因斯坦思考的是更深層的問題。
電磁場波動方程式是一個嵌入在時空框架的波動。愛因斯坦時空場方程式是時空框架內禀的彎曲,把引力場的運動方程式變成一個靜止幾何彎曲空間的地直線理論。

這個思考的邏輯如下:
始於等效原理。比薩斜塔丟下的鉛球,鐵球,鋁球,棒球,同時落地,為什麼?牛頓的解釋是萬有引力定律的重力質量 m_g 和運動定律的慣性質量 m_i 相等,F = \frac{GMm_g}{r^2} = m_i a or a = \frac{GM}{r^2} . 也就是運動本身和球的特性(質量)無關。愛因斯坦認為這不是巧合,其中有非常深刻的意義。

如果你不覺得這有什麼特別(其實我中學完全不懂有什麼特別),可以比較電磁理論的庫倫定律,F = \frac{KQq}{r^2} = m_i a . 形式看起來一樣,但電磁力和慣性質量毫無關係。其他的各種力(磁力,弱力,核力)也都和慣性質量無關。這也是愛因斯坦時空場理論並不是另一個類似電磁場 Maxwell equations 的波動方程式。

愛因斯坦對於比薩斜塔的實驗的解釋非常符合 Occam razor 原則。所有的球都同時落地,最簡單的解釋就是所有的球都沒動,而是觀察者做了一個反向的加速運動。愛因斯坦是用封閉電梯或是火箭的思想實驗,在沒有引力的外太空火箭頂端放了鉛球,鐵球,鋁球,棒球。火箭一加速對於關在火箭內的觀察者會認為所有的球同時落地。

因此我們可以把引力等效轉換成加速座標系的觀察者 A。相對於火箭外的靜止觀察者 B,自始自終這些球都是靜止不動。

我們先以 B 劃一個四維平直時空座標。靜止的球形成一條幾何上直的世界線(worldline), 對應這個平直空間是一條直線。在這個四維空間再劃上 A 的時空座標。因為 A 是等加速運動的觀察者。A 的座標系格線是雙曲線(t)和射線(x). A 對應的座標系不是平直空間,而是彎曲空間!更重要對 B 靜止的球的世界線在 A 的(彎曲時空)座標系的世界線也是幾何的 geodesic (地直線)。運動方程式對應幾何彎曲時空的 geodesic!

簡單結論:地球質量 \to 引力場的運動方程式 \to 加速運動觀察者相對靜止(等效原理+馬赫相對性原理)\to 加速運動觀察者的座標系對應彎曲時空(符合四維時空座標無關:Lorentz transform?),運動方程式變成幾何的 geodesic (地直線)。

可以跳過引力:地球質量 \to 彎曲時空。運動方程式變成幾何的 geodesic (地直線)。我懷疑愛因斯坦根本不認為引力存在!只有時空彎曲。

簡單而言:質能造成時空彎曲。時空彎曲又會影響質能分佈。質能的波動(如黑洞合併)造成時空漣漪,也稱為引力波。時空場方程式一邊是質能,另一邊是時空曲率。

時空曲率是四維張量 (rank-4 to rank-2)

黎曼曲率是 R_{\rho\sigma\mu\nu} 四維而且是四階張量。亦即一共有 4x4x4x4 = 256 個分量。相反質能只是一個純量。愛因斯坦面臨的問題是如何讓方程式兩邊能對接。也就是說:
curvature = \kappa (mass)  \kappa 是比例常數。包含萬有引力常數 G .

Ricci 之前把四維四階張量轉換成四維二階張量 (4×4). 因此方程式左邊是 4×4 張量。
質能被 energy-momentum (density) 四維向量取代 (4×1). 愛因斯坦進一步把 4×1 vector 改造成 4×4 energy-momentum tensor 如下。(Wiki 2019)

愛因斯坦一開始的引力方程式如下:
R_{\mu\nu} = \kappa T_{\mu\nu}  \kappa = \frac{8\pi G}{c^4} 是比例常數。 G 是牛頓萬有引力常數. 在弱場的近似解就是牛頓萬有引力方程式。

上式一個問題是無法滿足四維時空座標系 invariant. With some help from Hilbert, 愛因斯坦減去一個二階張量,得到座標系 invariant 的結果。
G_{\mu\nu} = R_{\mu\nu} - \frac{1}{2} g_{\mu\nu}R = \kappa T_{\mu\nu}

是在非慣性座標,非慣性座標系可以轉換為帶引力的慣性座標系,反之亦真。下一個問題是如何描述帶引力的四維時空座標系(NO)。
invariant!

  1. 在狹義相對論中質量和距離也和運動有關
  2. 牛頓的引力是超距力違背光速極速

Reference

Wiki. 2019. “Introduction to the Mathematics of General Relativity.” Wikipedia.
https://en.wikipedia.org/w/index.php?title=Introduction_to_the_mathematics_of_general_relativity&oldid=914452016.

語音信號處理–Mel Frequency Cepstrum Coefficients (MFCCs)

by allenlu2007

承上文 PhysioNet 2016 challenge, 最重要的 features 是 MFC (Mel Frequency Cepstrum) 係數的前 13 項。再加上 sample kurtosis and sample entropy.   到底 Mel Frequency Cepstrum Coefficients 是什麽?爲什麽有效?

本文主要參考: http://research.cs.tamu.edu/prism/lectures/sp/l9.pdf

以及 http://mirlab.org/jang/books/audioSignalProcessing/speechFeatureMfcc_chinese.asp?title=12-2%20MFCC

This lecture is based on [Taylor, 2009, ch. 12; Rabiner and Schafer, 2007, ch. 5; Rabiner and Schafer, 1978, ch. 7 ] 

Cepstrum

Spectrum 的定義是 20log10|F{x[n]}| or 10log10|F{x[n]}|.  Spectrum 把字首反排成 Cepstrum, 定義也很特別。簡單來説就是 Cepstrum = IDFT{Spectrum}.   

image

1. abs and log operations 自動把 phase information 消除 (but frequency information keeps). 對於 communication 是罪惡,但對 voice 似乎無妨。也許是耳朵對 phase 不敏感。但對 log (frequency response) 敏感。

2. abs and log operations 都是 nonlinear operations.  因此 c[n] 不會是 x[n] 的 replica 或是 linear filtering.  c[n] 的 harmonics 部分會加强。 c[n] is complex sequence?  如果改爲 DCT (discrete cosine transform) 就會變成 real sequence.

Examples

1. x[n] 是 pure sine wave.  After DFT, ABS, and LOG 仍然得到 c[n] 是 sine wave 如下圖左。rceps(sin[x])

很明顯一般的聲音不是如此。如果是周期 waveform with harmonics, 在 abs and log operations 之後,加强 harmonics 部分。在 DCT 之后加强 waveform 的 distortion.  例如 sawtooth waveform 最後得到的 cepstrum 如下圖右。基本上還是保持類周期性質,衹是 distortion 更嚴重。rceps(sawtooth[x]) 

image  image

Cepstrum 的新詮釋

前段提到 Cepstrum = IDFT{Spectrum}, 是從比較 x[n] and c[n] 角度出發。 

但在 reference  給了一個完全不同的詮釋如下。重點是: log spectra (power spectrum) 視爲 waveform!  IDFT 是用來分開 power spectrum 的 frequency and amplitude modulation.

image

下圖就是一個例子。log spectrum 在 log 后變得像 periodic waveform.  在 IDFT 后對應高頻部分是 power spectrum 上的 ripple.  如果在 cepstrum 加上 filter, 可以移除 ripple, 稱爲 liftering.  Liftering 有實際的用途嗎?

image

image

System Convolution

OFDM 一個最大的優點就是把 system convolution 轉換爲 multiplication (FEQ). 

同樣 Cepstrum 的一個很大的優點也是如此。把 system convolution 轉換爲 addition.

這對聲音很重要。可以分開 source and filter.  Source: glottal excitation, 對應 high frequency coefficients.  Filter: vocal tract, 對應 low frequency coefficients.

image

Comparison of LPC, STFT, Cepstrum smoothing, and Homomorphic smoothing

回到 frequency domain (spectrum), 我們可以比較幾種常見的語音算法,包含 Short-time Fourier transform (STFT), LPC (Linear Predictive Coefficients), Homomorphic smoothing, Mel-cepstral smoothing.  How about wavelet?  結果如下:

STFT: 似乎最大的問題是 ripple, 也就是無法分開 source and filter?

LPC, Homomorphic, and Mel-Cepstrum 看來 spectrum 都很像。重點還是在 parameter domain (time domain for STFT, qerfreny for cepstraum, etc.)  是否有明顯可辨識的特徵。

image

Voice detection 是常見的應用。以下是常見參數 for speech signal.

50-ms window, 12.5-ms shift,

Fs=8KHz, Nmfcc = 14 (number of cestral coefficients) , and R = 22 (cepstral sine lifter parameter).

In summary, Mel Cepstrum works 的主要因素是 harmonic rich signal!  經由 log operation 放大 harmonics, 再經由 DCT 得到 cepstrum spike.

image

SGD and SGD with momentum

by allenlu2007

 

Gradient Descent (GD, aka steepest descent), Stochastic Gradient Descent (SGD), Newton method, etc.  都是 optimization 的算法。Optimization 問題就是 (1) define a cost function C(x); (2) find C(x) 的 minimum (or maximum).

我們可以有不同的分類: 

1st order optimization vs. 2nd order optimization

用一個 example, C(x) = a x^2  (a > 0)

GD, SGD, etc. 都是 1st order optimization algorithm.  主要是利用 1st order derivative (gradient) information 做 optimization.

Gradient descent 如下:

C'(x) = 2ax  =>  x(t+1) = x(t) – r * C'(x) = x(t) – r * (2ax(t)) = x(t) (1-2ar)    r 稱為 learning rate

只要 -1< 1-2ar < 1   or  -2 < -2ar < 0  =>  1/a > r > 0   x(t+1) 就可以 converge 到 minimum.

一般而言,r 愈小 (r~0) 愈容易 converge.  但 converging speed 很慢。r 變大 converge 會變快,特別當 r ~ 1/2a 時收歙非常快。但 r 再大收歙又變慢 (overshoot oscillation), 如果 r > 1/a,又變成不收歙。

因此如何選 r 是 GD 最大的問題。r(opt) = 1/2a

由上文 r 的選取 highly depending on a.   也就是 2nd order derivative (曲率)。

如果是一個二維的 cost function, C(x, y) = a (x+y)^2 + b (x-y)^2  而且 a >> b or a << b 代表很扁的斜橢圓。這時 GD 就會有問題。因為只是 x and y 方向的 optimization 會被最 max(a, b) 所限制住。必須讓 r 很小,因此收歙速度會很慢。

如果更高 dimension 會更容易出現 bottleneck.  

Newton method 則大幅改善這個問題。基本的觀念就是使用 2nd order derivative 再加上 1st order derivative (gradient + curvature).  就可以大幅增加收歙速度。據 Stephen Boyd 說,任何 convex optimization 問題,使用 Newton method 只需要 12 次 iteration 就收歙到非常非常精準。實務上大約用 5-7次就夠解決 90% 的 convex optimization 問題。如果 cost function 是 parabolic, 只需要 1 次 iteration 就 ok.

Newton method or 2nd order optimization 最大的問題就是無法處理 non-convex optimization.  要不然就是 stuck 在 local minimum/maximum.  要不然就是無法收歙 (oscillation).

對於 machine learning 的問題,大多數是 non-convex optimization (e.g. neural network).  因此 Newton method or other 2nd order optimization 大多無用武之地。

 

Deterministic vs. Stochastic

前面提到的 GD or Newton method, 都是 deterministic optimization.  一般單純 deterministic optimization 問題的 dimension 是固定且不大,可以用 GD or Newton method.  

但對一些很多 training samples (high dimension input) 但 low parameter dimension 的 optimization 問題。如下圖的 curve fitting (or linear regression) 問題。cost function 就是 mean square error function.   可以看出來 input sample (dimension) 可能很大 (e.g. > 100), 但 parameter space 只是 2 dimensions (a and b).  

可以有兩種方式來解 optimization 問題。一是 GD, 不過要等到所有 samples (>100, e.g.) 收到,再做 iteration.

另一種方式就是 sample-by-sample iteration.  就是從 one sample (or two samples) 先算 (a, b), 每多一個 sample 就 iterate (a, b) 一次。 這就是 stochastic gradient descent (SGD) 的觀念。

NewImage

當然這也有假設 input samples 是 stationary 且 random.   可以想像如果一開始的 samples 都是集中在某一區,算出來的 (a, b) 可能就會有很大的誤差。

下圖是 (a, b) 經過 iteration 的 trajectory.  Red trajectory 是 GD,  Pink trajectory 是 SGD.  

SGD 的好處是不用等所有的 samples 就開始 learning, 因此收歛速度可能一開始比較快。但最後似乎 jittering a lot.  因此另一種 algorithm 是 mini-batch GD.  故名思義就是介於 all samples iteration 和 one sample iteration 之間。

NewImage

下圖是 mini-batch GD 的收歙 trajectory.  

NewImage

SGD 另外一個很大的優點是可以 on-line learning.  如果有新的 training set 進來,可以 incremental 改善原來的 weight.  而不需要全部重新 train/learn.  例如 ImageNet 是從幾張 sample photos 成長到現在 1.5M 張,之後還會再增加。如果用 GD 很難每次更新都重算,實務上不可能。

 

綜上所述,machine learning 所需要的 optimizer 是 SGD

 1. First order optimization,  因為 machine learning 的 optimization 絕大多數是 non-convex, 存在很多的 local miin/max.  

2. Stochastic optimization 而非 deterministic optimization.  因為 training set 太多而且會變動。

SGD 最主要的問題就是收歙太慢。特別是扁斜橢圓的 error contour.  因此有再針對如何加速 1st order optimization 的改進算法。 另一個問題是如何選 r (learning rate)?

 

Momentum Base Gradient Descent

以下參考此文 : “On the momentum term in gradient descent learning algorithms” by Ning Qian Columbia.

Momentum 

 

1. GD 基本上如之前所述

NewImage

NewImage

 

2.  GD 加上 momentum term (Rumelhart et al. 1986)

以下這段話的意思是: 加上 p*dW(t-1) 的修正項是避免 error vector (dW) 在橢圓短軸方向 oscillation 而影響收歙的速度。但如何選取 p and epsilon?

NewImage 

 

上文是 hand-waving 的說法,Qian 給了一個力學的想法。

GD 在 dt 趨近 0 時,可以把 iteration 的過程變成力學問題如下:

NewImage

剛才說過 epsilon 如果太大會造成收歙的問題,所以 epsilon 要小。當然以 Eq(3) 來看,epsilon 不管大小都會收歙,收歙速度是 exponential decay.  這是因為假設 dt 趨近 0, 所以 (epsilon * dt) 還是無窮小。

我們想做的是如何在同樣的 epsilon 可以改善收歙速度,特別是在 epsilon 很小的 case. 

在力學的經驗中,eq(3) 就是粒子在 highly damped 的力學公式。如果還原為完整的牛頓力學公式如下:

NewImage

m 是質量,u 是磨擦系數 (u ~ 1/epsilon or epsilon ~ 1/u).   

所以 (3) 是 (4) 的 special case when u >> m.  

此時因為磨擦力太大,gradient 所提供的吸力無法讓粒子加速,只能慢慢的滑入最低點。

不論是圓坑,橢圓坑,扁橢圓坑都可以收歙。唯一的問題就是慢。如何加速收歙?

 

有兩個方法: (1) 之前談過讓 epsilon 變小 (或增加 learning rate), 不過會有收歙的問題。

(2) 增加 m (多一個參數) 讓 gradient 提供的吸力可以加速粒子。

不過直覺想如果 m 太大,也會造成粒子 circle around 而不收歙。因此 m 也需要在一個範圍內。

m=0 就回到 GD,  m > 0 但必須小於某個 value 才有更快的收歙。

 

其實 (4) 的 discrete time 公式就是 Eq (2).

NewImage

NewImage

NewImage

如果 p = 0 ==>  m = 0.   epsilon ~ 1/(u*dt)

因此 GD with a momentum term 相當於一個牛頓力學粒子在 viscous medium under a conservative force field.

 

Energy Function

再從能量的觀點來看。之前的 Energy function (or cost function) = E(w).

加入 momentum term 之後, total energy function (or cost function) 變為: 動能+位能

NewImage

加上動能項可以改變原來的 energy function (cost function) 的 profile.  例如之前很難收歙的扁橢圓在加入動能後可能就不是扁橢圓而加快收歙速度。另外有 m (m > 0) 之後也會改變每次 weight change (=velocity=dw/dt) 的方向。

當然最後 dw/dt = 0,   Et –>  E(w),  所以改變 cost function (energy function) 對最後的 optimization 並無影響 。但是可以改變收歙的速度。

如果 m 讓動能 >> 位能,那就像沒磨擦力的 case, 也要很久才收歙。反過來如果動能 << 位能,又回到磨擦力 dominant case.  所以最佳的情況似乎是讓(initial)動能 ~ 位能。最終當然動能=0, 位能-> min 位能。

 

Qian 的 paper 給了一個定量的解釋如下:

NewImage

假設 H (Hessian, symmetry and positive definite) 可以被分解成 H=QKQ’

NewImage

NewImage

Equation (18) 對應是 GD 的 exponential decay time constant (m=0).

GD with momentum (m <> 0) 則變成 2nd order differential equations.

NewImage

NewImage

以下是結果。如果選 m 滿足 eq(23) 可以加速 exponential decay rate. 

NewImage

 

可以快多少。下面 Result2 告訴我們 eigenvalue 可以最大變 X2.  當然 exponential decay rate X2 代表平方 (e.g.  原來的 GD error = 0.1,  GD with momentum error 最好可以 (0.1)^2 = 0.01. 

NewImage

 

Discrete Case

上節是用 continuous and 力學 analogy 來解釋 momentum GD. 

回到 real discrete case, 也可以得到收歙加速的效果。

NewImage

NewImage

和 continuous time 不同。discrete time 的 eigenvalue 愈小愈好,同時需要小於1.

NewImage

NewImage

 

NewImage

結論就是 0 < p < 1  以及  p < 1- sqrt(epsilon*ki).   同樣最小 ki 也是 bottleneck. 

NewImage

 

可以參考 Hinton 的 optimzation note

NewImage 

p = beta = 0.95

(1-beta) epsilon = 0.05*epsilon

看來 epsilon 一般會設的很小。 

 

————————————————————————– 

以上部份在另文討論。

SGD 有 noise average 好處? 

Conjugate Gradient Descent or Hamiltonian Dynamic GD

Conjugate GD –> Solve Energy Optimization problem –> Leverage Hamiltonian dynamic

 

SGD  Check Zaid’s talk in CVPR 2013 (but no momentum)

SGD with momentum

Hamiltonian dynamics 

MCMC 

 

The other thing is Energy function

Stanford Machine Learning and Deep Learning Lecture and Others

by allenlu2007

 

Stanford 可說是 machine learning 的大本營。

原因是 (1) 幾大公司急需相關 knowlege or expertise 都在附近, 如 Google, Facebook, Apple, Microsoft (research center). 學生老師開公司, 找工作, 或 consulting 都容易。(2) 對於新領域,特別是有具大商業價值的 area.  Stanford 一向都敢為天下先。

 

CS229  Machine Learning by Andrew Ng 不錯的基本課程,適合初接觸 ML.

Andrew Ng (CMU->MIT->Berkeley).  之前和 Google 合作 Google Brain.  目前在百度,是 Baidu Brain 掌門人。

Cover learning 的三個領域: supervised learning, unsupervised learning, reinforcement learning 以及相關的方法。

 

CS228/228T  Probabilistic Graph Models (PGM) by Daphne Koller.  

PGM 的進階課程。我並不是那麼喜歡 Koller 的上課方式。Koller no doubt 是 graph 的 master, 但常常沒有給一個明確目標 why to do this 或解釋可能的應用 (如 variable elimination lecture), 很容易讓聽者陷在 tedious 的公式推導中,有點像是 textbook 的 replay.  後來發現她的確寫了一本 PGM textbook and on-line readable.太重 details.  缺乏比較直覺的說明。

相較之下,我比較推薦先看 MLSS 的 Prof. Sam Roweis 在台北的 tutorial lecture about factor graph model.

Roweis 四個 lectures (MLSS in Taipei) 是我聽過最好的 PGM lectures.

 

Lecture1: introduction of PGM; conditional independent == factorization

Lecture2: 

Lecture3: latent variables, clustering

Lecture4: inference, belief propagation, junction tree: 非常棒,除非你已經精通 PGM, 應該會有豁然開朗 的感覺。  p(x|e) 就是用 query node as root (x, single variable), evident node (e, multiple variables) as leaves, trim out other leave nodes (sum to 1 for non-evident leaves).  Then use message passing to from leaf to root and root to leaf.   For multiple queries (again with single variable), do the messaging passing for every node back and forth (x2) obtains the entire tree CDP. 

http://videolectures.net/mlss06tw_roweis_mlpgm/

Roweis 是 1999 Caltech PhD, 待過 MIT and University of Toronto, then NYU professor.  I am so surprised to know he passed away in 2010.  What a pity to such a brilliant person!

Lecun 和 Roweis 都是 NYU 的 machine learning professors and collaborators. 

 

Neural Network and Deep Learning

Ng and Koller 的課程對 Neural network and deep learning 都沒有著墨。Ng 之後的 CS229 似乎加入了 neural network and deep learning 的介紹。不過新的課程還沒有放在 youtube.  Koller 主要 focus on graph model, 對 neural network 或 deep learning 似乎沒有興趣。

Deep learning (DL) 其實是 neural network 的更新版, 是 machine learning 的一個 branch, 也是近年的當紅炸子雞。在 (automatic) speech recognition (ASR/SR), computer vision (CV), nature language processing (NLP) 都有不錯 or 突破性 (?) 的進展。配合 Internet 公司在這方面的需求。Neural network 可說是鹹魚翻身,又變成顯學。不過之前 neural network 的 reputation 太差 (眼高手低)。現在只稱為 deep learning 表示和之前 neural network 不同。

Neural network 可分為三階段 (每階段約 20 年)

Gen1 1960 Cornell:  Perceptron, no hidden layer, hard-decision.  Feature by hand.  Learning for weight.

Perceptron 之後改良為 Logistic regression. 差別是 soft-decision instead of hard-decision.

In general, kernel function y = x -> linear regression.   y = sgn(x) -> perceptron.   y = sigmoid(x) -> logistic regression.   其實不同的 kernel function 對應不同的 cost function.   Ng CS229 note 的 machine learning 有提到可以 generalize to exponential family distribution.   What’s the relationship between deterministic classifier function with pdf distribution?   linear regression, perceptron, logistic regression,  softmax/multinomial regression (generalized logistic regression, multi-label classification).

 

Gen2 1980:  Shallow learning neural network.  Back-propagation learning on weight and feature?

 

Gen3 2005:  Deep learning neural network.  Both feature!! and weight by learning.  因為之前 neural network 的 reputation 太差 (眼高手低)。大家多半只稱為 deep learning 和之前 neural network 區別。

不知道 Gen3 是否能突破性的解決之前無法解決的問題。不然可能又要等 Gen4 的突破。   

Deep learning 主要 contributors 包括 Geoffrey Hinton, Yann Lecun, and a Japanese (Fukushima?).

 

Coursera:  Neural Network for Machine Learning by Geoffrey Hinton (Google Brain 掌門人)

Hinton (U. of Toronto) 本身就是 Gen3 deep learning 的主要 contributor.  他的 coursera lecture 說實在我沒有很深的印像。可能是太循序漸近,還沒到 DL 我就失了耐心。等我全部看完再給 comment.  另外 MLSS 也有很多 DL tutorial.  之後再 comment.

 

NYU Deep Learning Lecture by Yann Lecun link (Facebook AI 掌門人)

Yann Lecun 也是 Gen3 deep learning 的主要 contributor.  他的 NYU lecture 以及 techtalks 聽起來比較有趣也比較快。Lecun 和 Roweis 都是 NYU 的 machine learning professors and collaborators.  不過 Lecun 的 lecture assume 你熟悉 machine learning 的基本技 (regression, classification, kernel method, K-mean).  最好已經聽過 CS229 再聽 Lecun lecture 會收穫更多。

 

 

CVPR 2013 “Large Scale Visual Recognition With Deep Learning” by Marc Ranzato

Lecture link.   另外 Youtube 也有其他的 video lecture.

Google Deep Learning, Marc Ranzato, NYU (PhD with LeCun) 和 Toronto (postdoc with Hinton) 都待過,目前在 Google。他是義大利人,不過英文通暢。他的 slide 最合我胃口,由淺入深。E.g. CVPR 2013 slide,不像有一些 deep learning or machine learning 學者,一下就跳的很深。同時他對 probabilistic (graph model) 也非常熟,和 deep learning 可以做深入比較。

 

Google (NYU -> U. Toronto)  Marc Ranzato

Youtube lecture 2012

and other video?

Ranzato 的 slide 非常清楚,明白易懂。之後看完再 post.

 CVPR 2013 slide,不像有一些 deep learning or machine learning 學者,一下就跳的很深。同時他對 probabilistic (graph model) 也非常熟,和 deep learning 可以做深入比較。另外他談到很多 practical 問題,非常有用!   非常推薦 watch his video.

Video: http://techtalks.tv/talks/large-scale-visual-recognition-part-vi/58589/

關於 why deep, difference between PGF and neural

 

 

CS231n: Convolutional Neural Network for Visual Recognition by ??

CNN for visual recognition or computer vision 是近來很熱的 topic, 可以用於無人車和無人機的駕駛。可惜 Stanford 這門課沒有 video broadcasting.  我又太懶沒去看 lecture note.  雖然是我非常有興趣的題目,只好之後再說。

 

CS224d: Deep Learning for Natural Language Processing by Richard Socher

我聽了 lecture 1, 應該很有趣。Richard 的講課也很生動。不過因為 NLP 和我目前想看的應用無關。先 skip.

 

Nvidia Deep Learning Tutorial

 

Factor Graph – Kalman Filter

by allenlu2007

之前花了很多時間了解 Kalman filter (見前文)。主要在 recursive equation 的推導,physical insight (least square minimisation, state space model of Gaussian distribution).  

進階包含: Kalman filter 和 Kalman smoother 的差異.  以及當 variance 未知時用 EM algorithm 來求解。

用 HMM 來類比 Kalman filter:  alpha-beta (forward-backwar) algorithm 相當於 Kalman filter/smoother.  Baum-Welch algorithm 相當於 EM algorithm.  另外還有 Viterbi algorithm for HMM sequential optimisation.

 

之後 study probabilistic graph model (PGM) 以及 factor graph.  如 Bishop 所言其實 Kalman filter 只是一個 PGM/FG 的簡單的應用。Kalman filter 放在 PGM framework 下是一個單純的 chain (special case of tree).

HMM 的 forward-backward algorithm 和 Viterbi algorithm 正好對應 PGM 中 message passing 的 sum-product 和 max-product algorithm.  同時在 Gaussian message 和 Gaussian factor (linear combination) 的情況下,sum-product algorithm 等價於 min-sum algorithm, 也就相等於 max-product algorithm.

PGM 可以把所有的 algorithm 都整合在同一 framework, 例如 EM algorithm, 例如 parameter estimation. 更重要的是 PGM 更 general 的 message passing algorithm 可以處理各種特別 model.  例如 switched state space model, 例如 time varying model, 例如 input control model, etc.  Kalman filter 原來主要亮點是 recursive equation; 在 PGM message passing 只是最基本的 sum-product message passing 結果。

 

參考 Bishop video

先用一個非常簡單的例子來說明不同 graph 的優缺點。

Example 1: 在一個 noisy (zero mean with known variance w) 環境,estimate dc 訊號值。

Solution 1A: 直覺來說,最佳值就是所有 observables 的平均值,empirical mean:  u = (x1+x2+ … + xn) / n.

如果 noise variance 是未知,最佳的 variance estimation 是 empirical variance:  sum(xi – u)^2 / (n-1)

乍看之下不需要更複雜的 probabilistic graph model 或演算法。後面的例子會 show 這些 model 和算法的價值。

 

Solution 1B: 假設 noise 是 Gaussian distribution with zero mean.  Maximum likelihood parameter estimation 對 likelihood function 微分,可以得到同樣結果。另一個假設是 xi 是 independence (i.i.d.).  

 NewImage

NewImage 

 

Solution 1C (Plate model): 以下我們用不同的 approach, 就是 Bayesian approach, 基本上就是有 prior distribution.   我們用以下的 model:    假設 mean 是 Gaussian distribution p(u) with mean u and variance v.  我們最後的目標就是 estimate u (maybe w in some case?) 根據觀察到的 observables.   Observables 包含 sensor noise (noise with zero mean and variance w)  p(x|u) 同樣是 Gaussian distribution.

注意 noise variance w 和 mean 的 variance v 是兩個不同的概念。Frequentist 或 maximum likelihood lover 可能會困惑為什麼需要假設 mean 是一個 Gaussian distribution 而非一個 fixed parameter.  Maximum likelihood 對於 parameter estimation 沒有 prior information 假設。

我們用以下的 factor graph 來 model: mean and observable 的關係。

NewImage

 

Include observables: it becomes joint Gaussian distribution.  Given u, all xi are conditional independent.   

NewImage 

 

u 的 total message = product of incoming messages = exp((u-uo)^2/v]  exp[(x1-u)^2/w] * exp[(x2-u)^2/w] * … * exp[(xn-u)^2/w)]

其實可以看到就是 u 的 marginal pdf:  N(uo, v) * N(x1, w) * N(x2, w) * .. * N(xn, w)

mean = (uo * w + x1 * v + x2 * v + .. + xn * v) / (w + n * v) = (uo * w/v + x1 + x2 + … + xn) / (w/v + n)

variance = (v // (w/n)) = ( n/w + 1/v )^(-1) 

uo and v 是 prior information (mean and variance) of u.  如果事先知道 u 的 distribution, v 很小: mean ~ uo

相反如果不知道 u 的 distribution, let uo=0 and v 很大:  mean ~ (x1+x2+ ..+xn)/n 和 empirical mean 相同。

如果隨便選 uo 和 v, 當 n 夠大時,仍然會收歛到 empirical mean.

 

另外也可以從 conjugate prior 出發,得到相同的結果。

Conjugate prior 已經開始引入 sequence 的概念。新的 observables 會變為 prior information 的一部份。

NewImage

 

意即 Bayesian interpretation 多引入一個 prior distribution, N(uo, v), 的彈性。傳統的統計學,可能很難了解為什麼需要 prior information, 應該 “讓証據說話”  

實務上 prior information 是一個有用的觀念,很多情況我們對要 estimate 的 parameter/distribution 有 prior information.

如果真的沒有 prior information, 就讓 prior distribution 的 variance 變大,代表 prior information 不清楚。

另外一個好處是在傳統的 machine learning, inferring 和 learning 是兩件不同的事。Learning 多半是 parameter estimation.

在 Bayesian interpretation 中,inferring 和 learning 都是 inferring.

Bayesian 還有其他好處嗎?

 

 

Solution 1D:  導入 Hidden Markov Model (HMM)

Prior distribution —> initial distribution of the estimated parameter

引入 recursive 觀念。新的 observables 持續 update inferring.  並導入 conjugate prior. 

意即 prior information 可視為之前 observables 所得到的 information.  同一種 pdf family.

  

Or hidden markov model (HMM).  HMM 是更 general 並且有用的 model.

(i) sequential observables, 可以建立 recursive equation.

(ii) 在 DC model 中, u(k+1) = u(k).  更 general case, u(k+1) = A u(k) 可以是 linear transformation (或加上 input signal).

 

NewImage

 

Forney factor graph 的表示方式: 灰字代表 unknown. 黑字代表 known or observables.

NewImage

u 是 Gaussian distribution with mean uo and variance w.  這是 initial distribution 或是 prior information.

Based on 之前的公式:  

u1 的 mean : (uo w + x1 v) / (v + w);   variance = v w / (v + w) = v // w (v and w 並聯,小於 v and w)

u2 的 mean : mean(u1) w + x2 * v // w ) / (w + v // w ) = (uo w + x1 v + x2 v ) / (w + 2v);  variance = v // w // w

…  

結果和 Solution 1C 一致。比起 1C 的好處是 (i) recursive; (ii) 可以引入更 general 的 linear transformation chain model. 

補充說明: Solution 1D 相當 forward message passing, 也就是 Kalman filter.  如果做 backward message passing.  

再結合 forward and backward messages, 就得到 Kalman smoother.  以本例來說:

Backward: x2 的 mean and variance 回傳給 u1

u1_back mean = x2;  variance = w  

u1_forward mean = (uo w + x1 v) / (v+w);   variance = v // w

Combine forward and backward message ==> equivalent to Kalman smoother

u1_total mean = [ x2 * v // w + (uo w + x1 v) / (v+w) * w ] / ( w + v // w) = (uo w + x1 v + x2 v) / (w + 2v)

u1_total variance = v // w // w

以上的結果和 u2 的 mean 和 variance 一樣。這也是我們預期。因為 P(u1 | x1, x2, u) = P(u2 | x1, x2, u).

再次強調 Kalman smoother 並不是 backward message.  而是 combine forward and backward message!!

結論是我們不用再記複雜的 Kalman filter or Kalman smoother 的 recursive equation!  只要 apply forward (and backward if we want to retrospect) message 的 sum-product (equivalent to max-product) rules!!

 

 

Example 2: 在一個 noisy (zero mean with unknown variance) 環境estimate dc 訊號值以及 noise variance。

為什麼需要知道 noise variance?  在一些應用如通訊,估計 SNR (signal to noise ratio) 是非常重要且有用的參數。

 

Solution 2A and 2B 和 1A and 1B 一樣。用 empirical variance 或是 maximum likelihood 可以估計 noise variance.

從 PGM 或是 Bayesian approach, 也有兩種 model.  一種是 lumped 的 plate model, 另外是 HMM (不過我們稱 EM  algorithm!)

 

Solution 2C (Plate model)

下圖,意即除了 mean 的 prior information P(μ) 之外,還包括 variance (σ^2) 的 prior information.  比較方便是定義 precision λ = 1/σ^2.  λ 的 pdf 是 P(λ).

P(μ) 假設為 Gaussian distribution with mean uo and variance v.

P(λ) 假設為 Gamma distribution, 因為 variance (1/precision) 只能為 positive number.  Square of Gaussian distribution is Chi-squared distribution.  Gamma distribution 是更 general distribution.   Chi-squared distribution 是 Gamma distribution 的特例。Needless to say, Gamma distribution 仍是 exponential family.  

Gamma(α, β)  α 稱為 shape parameter, β 稱為 rate parameter

NewImage

E(X or λ or Φ) = α/β;   Var(X or λ or Φ) = α/β^2  這是 precision 的 pdf

Variance (σ^2) 的 distribution 稱為 Inverse Gamma distribution.

E(σ^2) = β / (α-1) ~ β/α     Var(σ^2) = β^2 / (α-1)^2(α-2) ~ β^2/α^3

在 α >> 1 時,  precision and variance 的 mean 基本上為倒數。但 variance 顯然不是。

NewImage

因此 factor graph 如上圖。Joint pdf P(μ, λ, x) = P(μ) * P(λ) * P(x|μ, λ)

乍看之下,μ and λ 是完全 independent random variables.  不過這只有在 x 和 x 的 descendents 都是 unknowns 才如此。如果 x 或 x descendents 是 observables, μ and λ 就變為 dependent.  這可由 Bayesian directed graph explaining away 解釋。

Joint pdf P(μ, λ, x) = P(μ) * P(λ) * P(x|μ, λ)

P(μ, λ) = P(μ) * P(λ)  (是 independent)

P(μ, λ|x) ≠ P(μ|x) * P(λ|x)  (Given x, μ and λ 變為 dependent) 

以上的 model 需要四個 parameters fully characterized, 分別是 μ 的 mean and variance; λ 的 mean and variance.  不過會有一個問題是無法得出一個 general conjugate prior, why? 因為 real data 的 mean and variance 是 dependent? 如果用 independent mean and variance 來 model prior information 會有問題?

應該是 conjugate prior 有問題,但仍可以用以上的 model for message passing?

另一種 model 的出發點是既然 P(μ) 是一個 Gaussian distribution.  同時 P(x|μ, λ) 也是 Gaussian distribution with mean μ.  就讓 P(μ) 儘量接近 P(x|μ,λ).   σ2 和 λ 剛好是倒數。所以 σ2 是 inverse-Gamma distribution.

另一種 model 是讓 μ and λ 有如下的 dependency.  同樣有四個 parameters,  μo and no (for μ), νo, and σo (for σ).  這是 conjugate prior 的做法。

NewImage

NewImage

no 是 prior sample size 還蠻有 physical insight!  代表 prior information 相當於 no 個 observables.

這 no observables 的 mean 是 μo, variance 是 σ2.  不過 σ2 pdf 又引入另一個 νo (degrees of freedom) 以乎和 no 有重覆。

假設 no 和 νo 一致,那就簡化到三個參數: μo (prior mean), σo^2 (prior variance), (no,νo) (prior sample size).  換言之讓 u, σ dependent, 可以簡化到三個參數,但 dependency 會讓 message passing 更複雜。

 

另一篇 paper 更清楚, 同時加上了 n observations.  最後的 prior (p(μ,Φ))+ posterior (p(Y)) joint pdf 如下。就 graph model 而言,就是 plate model.

NewImage

NewImage 

NewImage

NewImage

當 samples 愈大 (n):

Normal pdf 逼近 delta function 收歛到 empirical mean (y_bar).

Gamma pdf 的兩個參數: shape (α) and rate (β).  當 samples 愈大 (n), shape 愈大,rate 愈大(scale 愈小)。Gamma pdf 會逼近 delta function (可用 python or matlab function plot).  Φ (precision) 會收歛到 Gamma function 的 mean = α/β = νn/SSn —> 1 / (empirical variance)

所以 Inverse Gamma pdf 的 mean 收歛到 (SSn/νn) ~ SSo/n + SS/n + po/pn ( y_bar – mo)^2  ~ SS/n  意即 empirical variance!

Gamma pdf 的 variance 收歛到 SS/n^2 ~ 0 意即最後變為 delta function.

 

上述公式的詮釋

po 相當於 prior information 等價於 po 個 samples

mn 是 prior information 和 observables 平均的結果

νn degrees of freedom, 應該和 po 有相同的功效

SSn 包含 3 個部份:

SSo: prior information of variation

SS: observables 的 variation (sum of squares)

最後一項是: prior mean 和 sample mean 的 variation

 

Solution 2D 導入 Hidden Markov Model (HMM)

上節 (Solution 1D) 從 plate 到 chain 提供更 general 的 model.  (a) recursive equation; (b) linear transformation.
 
 

下圖的 μ 為 x(k),  Φ 為 Σ(k),  x 為 y(k)

NewImage

 

其實上圖就是 Kalman filter 的推廣。以下為更 general Kalman filter 的 problem statement.
Estimate DC signal with unknown variance,  只是 A=1, Q=0, H=1 的 special case!
NewImage

 

我們主要先考慮 VB inference 因為可以把 noise variance 和 observations 連結一起。

目標是計算 posterior distribution:

 

NewImage

 

Initialisation and prediction 基本上都和 Kalman filter 一致。困難的是 eq. (3) 的 integration.

密訣是假設 mean 是 normal pdf;  variance 是 Gamma distribution 同時 mean and variance 是 independent.

NewImage 

By Bayesian directed graph rule, 上述的 conditional mean and variance 非 independent.  

我們仍用這樣的近似。

NewImage

NewImage 

NewImage 

NewImage 

注意以上的 equations 是聯立方程式 (mk, Pk, Σk are unknowns 但有三個 equations to solve).

為了避免解聯立方程式,我們再把每個 update phase 做 N 次 iteration 來解如下。

Predict phase 和之前 equation 相同,不過引入了 noise 的 heuristic dynamic parmeter ρ. 在 time varying noise 時 ρ 的角色很重要,更詳細的說明請見後文。 

Update phase 就是用 fixed N iteration 來解上述 equation (7). 

NewImage 

ρ = 1 : stationary noise variance.  ρ < 1 表示 precision Φ (1/σ^2) 的 variance 增加。

NewImage 

  

可以看到引入 HMM 比 plate model 更 general.   (1) 是 recursive method 而非 batch method (需要等到所有的 measurements).  Recursive method 非常適合 real time processing.   (2) 更 general 可以加上 process noise (Q > 0), dynamic model (A != 1).  以下是 simulation results.

Q=0 or R >> Q   (mean=-3.8; Q=1e-5; R=0.81)  

下圖是 mean 的 true value (state), measurements, estimation. 和普通 Kalman filter 不同的是 noise variance R is unknown. 是由 equation (7) estimate.  

 NewImage

下圖是 noise variance R 的 variance 以及 break down to empirical variance and P (residue variance).

前文解釋過: 開始時 K 接近 1, estimate mean 接近 measurement, empirical variance 非常小。

residue variance (P) ~ R/n.  這由 log plot 可以看出。

隨著 estimate mean 接近 true mean, K 接近 0, empirical variance dominate.  Residue variance -> 0. 

NewImage

Log scale. 

NewImage

 更 general case 請參考 Adaptive Kalman filter.  包含 time varying noise and time varying signal. 

 

本文主要強調使用 PGM 的威力。可以處理更複雜的問題。只要能用 PGM model, 就可 inference model 中的 parameters.   How to explain EM algorithm?

Forney Factor Graph – Simple FEC (part 2)

by allenlu2007

前文用 FFG 來 decode FEC,主要重點

1. forward and backward message passing to compute 所有 marginal probability

2. Iterative the message passing: probability 可以收歛 (desire case: probability 收歛到 0 or 1. 但也有可能收歛到 0.5!

為什麼需要 iteration? 如果是 Y=(0, 0, 1, 0), 一次的 message passing 就相當正確判斷 codeword 是 (0, 0, 0, 0). 其實可以不用 iteration.  但在某些 initial condition (e.g. Y=(0, 0, 1, 1) 或是 graph structure (loop?), 一次的 message passing 無法判斷valid codeword, 更不用說收歛。藉著 iteration, 可以讓 probability 收歛。

3. 如果有必要,initial probability 加上 perturbation 可以破壞對稱性而幫助 probability 收歛。

最大的問題是不管 max-product 或 sum-product with iteration 即使收歛也並不保証收歛到 valid codeword!

因此是否有其他的 algorithm 可以保確收歛到 valid codeword.  答案是肯定的,就是用 state transition graph or trellis graph/diagram.  

 

State Transition Graph

以 Hamming code (7, 4) 為例。對應的 parity check matrix 如下。

NewImage

H 對應的 valid trellis diagram and sequence state transition graph 如下圖。暫時先不管下圖如何產生。

NewImage

State transition graph 也是一個 FFG.  Variables 包含 X1, X2, .. X7, and S1, .. S6.

S1, S2, … S6 代表 states:  S1 and S6 只有兩個 states {0, 1}, S2/S3/S4 有四個 states {0, 1, 2, 3} (由下而上編號)

Factors 包含 f(X1, S1), f(S1, X4, S2), f(S2, X3, S3), f(S3, X2, X5, S4), f(S4, X6, S6), f(S6, X7).

f(X1, S1) = 1 if (X1, S1) = (0, 0) or (1, 1);  else 0

f(S1, X4, S2) = 1 if (S1, X4, S2) = (0,0,0), (0,1,2), (1,0,3), (1,1,1); else 0

f(S2, X3, S3) = 1 if (S2, X3, S3) = (0,0,0), (0,1,1), (1,0,1), (1,1,0), (2,0,2), (2,1,3), (3,0,3),(3,1,2); else 0

f(S3, X2, X5, S4) = 1 if (S3, X2, X5, S4) = (0,0,0,0),(0,1,1,2),(1,0,1,1),(1,1,0,3),(2,0,0,2),(2,1,1,0),(3,1,0,1),(3,0,1,3); else 0

f(S4, X6, S6) = 1 if (S4, X6, S6) = (0,0,0), (1,1,0), (2,1,1), (3,0,1); else 0

f(S6, X7) =1 if (0,0) or (1,1); else 0

Total branches = 8+8 = 16 = 2^4 matches with 2^4 codeword space! 

In summary,  state transition graph 的 FFG 是 chain structure; 和前文由 H generation matrix 直接產生的 web structure FFG (Tanner graph?) 不同。目的是 always decode valid codeword.   How to do it?

 

Channel Model

Chanel model 的方式和前文相同,分為 memoryless 或同樣有 state-space channel models 如下圖。

 NewImage

NewImage

很明顯如果是 memoryless channel, 結合 FEC state transition graph 和 memoryless channel graph 的 combined FFG 變成一個 Hidden Markov Chain (HMM).

如果是 state-space channel model (e.g. partial response channel or ISI), 則變成 dual chain 的 FFG.  BTW, dual chain 存在 loop, 只有近似解。

本文只考慮 memoryless channel. 理論上 HMM 的 algorithm 都可以用在解 state transition graph/trellis graph.

最終是要算出 P(branch | Y1, Y2, ..Yn) = P(b|y) 而得到 valid codeword.  

It is proved by Wiberg that:

Using sum-product to compute trellis graph P(b|y):  BCJR algorithm

Using max-product to compute trellis graph P(b|y): soft output Viterbi algorithm (SOVA)

 

Recap Baum Welch Algorithm (1963)

NewImage

Baum Welch algorithm 是一種 EM (expectation maximization) algorithm.  在 (A, B, Pi) unknown 情況下,given Y1, Y2, .. Yt, 找出最有可能的 (A, B, Pi).   EM algorithm 的 iteration, 就是在增加 likelihood.

表面上 Baum Welch algorithm 和 trellis decode 不同。在 trellis graph, A 是 given by trellis structure (如上圖), either 0 or 1.  BCJR 卻發現 Baum Welch algorithm 其實可以用在 trellis decoding.  

Iterative decoding 就如同 EM 的 iteration (?)  似乎 No!

EM algorithm is equivalent to Iterative decoding?  似乎No!

BCJR and Baum-Welch algorithm 相同嗎? 似乎 No!

 

Algorithm summary (from Wiki)

NewImage

 

 

BCJR Algorithm Review (1974)

BCJR 是用於 trellis code decoding, 主要是 convolution code, trellis code modulation, 但如前文討論,也可用在 linear block code decoding.  For practical reason, BCJR never have been adopted as widely as Viterbi algorithm (VA).

Con:

* Lots of computation (product and sum) compared with VA.  VA only needs to add-compare-select (ACS).

* VA only need “forward only” and store a path history of each survivor.  This is particularly useful for real-time communication purpose.

* VA does not need to know the SNR.  BCJR needs to know SNR a prior.  

* Linear block code decoding 一般是用 linear algebra based computation (syndrome calculation, etc.) 而不用 trellis based algorithm.  

 

以下節自 Abrantes’s “From BCJR to Turbo Decoding

首先 encoder/channel/decoder model

NewImage

再來是 prior and post probability

NewImage

NewImage

NewImage 

Key 是 log-likelihood ratio L.  可以証明 (by using Markovian property) 

NewImage

NewImage

NewImage

NewImage

NewImage

比較 BCJR and Baum-Welch

Alpha

BCJR: alpha(k-1)(s’) = P(s’, y<k)

BW: alpha(t)(i) = P(Xt=i, Y1, … Yt)

基本上相同,但 BCJR 往前一個 time step

 

Beta

BCJR:  beta(k)(s) = P(y>k | s)  

BW:  beta(t)(i) = P(Yt+1, … YT | Xt = i)

相同的定義

 

Gamma

BCJR: r(k)(s’, s) = P(yk, s | s’)

BW: P(Xt=i | Y)  

gamma 在 BCJR 和 BW 的定義完全不同。

BCJR 基本上把 Y1, Y2, … YT 拆成三份: 1.. (k-1) (Alpha), k (Gamma), (k+1), … T (Beta)  是 branch metrics.

BW 基本上把 Y1, Y2, …, YT 拆成兩份: 1..t (Alpha), and t+1, …, T (Beta); Gamma 比較像 temporary variable.

 

 

NewImage

注意只有 recursive (traverse the Markov chain), 並沒有像 EM 的 iteratively update!

 

Numerical Examples

Convolution code with code rate 1/2, 4 states, and the trellis graph shown below.

NewImage

NewImage

Ec/No = 1dB (SNR).  

順序是先 forward pass, 用 Gamma 來算 Alpha, 等到 forward pass is done, start backward pass for Beta.

NewImage

 

NewImage 

NewImage 

NewImage 

計算完所有的 P, 下一步就是算 L(uk | y).

NewImage

 

 

BCJR (or MAP) Summary

1. BCJR is strictly based on the trellis graph.  因此不會 decode 出 invalid codeword!

2. BCJR needs to know the SNR to compute Lc (?) 如果 SNR  is unknown 會有問題嗎?

3. BCJR needs lots of product and summation computation.

4. BCJR needs to do forward and backward, cannot be used for real-time computing.

NewImage

因些很少用在實際應用上。Solutons:

Use log to replace product -> log-MAP algorithm

NewImage

 

Use max function to approximate ->  max-log-MAP algorithm (Viterbi? Some yes, some No)

NewImage

Use max-product (and log) instead of sum-product (MAP) -> Viterbi algorithm (VA)

Store forward pass and not use backward -> VA

Without the need of SNR, maximizes r*v -> VA

VA 完勝 BCJR on the trellis decoding!

 

Q: Difference between VA and SOVA

A: VA is only forward and store.  SOVA does both forward and backward.  Reference

 

Q: Why BCJR?  

A: The MAP (per bit) for trellis decode.  Performance 會比 VA or SOVA 好? Reference 

傳統的 turbo decoding 的演算法則是由 Bahl,Cocke,Jelinek,及 Raviv 所提出之
所謂 BCJR 演算法。BCJR 演算法主要觀念在解一個最大事後機率(MAP)檢測問
題,因此基本上 BCJR 演算法就是一個 MAP decoder。BCJR 演算法運用到 turbo
decoding 的運作方法是利用事後機率來判斷個別 bit 應該隸屬於 0 或是 1,然後
再重建整個原來的資料序列,這一點與 Viterbi decoding 找一條最相似的的資料
序列作為原來的資料序列的解碼方式不同,由於 BCJR 演算法是找一個最佳 bit
的解碼方式,因此就平均的位元錯誤率(BER)來講 BCJR 演算法會比 Viterbi
decoding 來的好。然而,由於 BCJR 演算法的複雜度太高(詳細 BCJR 演算法請看
參考文獻一)不利於硬體實現,因此一些化簡的動作是必須的。
最常的化簡的動作可以歸納成以下 5 個規則如下:

1. 使用 likelihood ratio 取代在 trellis 中的機率值。
2. 分別在 likelihood ratio 或是機率值上取”log”值,如此一來

NewImage

3. 只考慮最鄰近的 trellis path。也就是 max* –> max
4. 在計算時只考慮部分的 trellis path。
5. 因為 BCJR 演算法必須全部的 codeword 接收後才能進行解碼,可以設定 window長度來解碼。

 

規則 1 與規則 2 並不會傷害 soft-information,然而規則 3、規則 4
與規則 5 可能就會降低整個系統效能。

如果我們採用規則 2 與規則 5 來簡化 BCJR 演算法則稱為 MAP 演算法,若我們
更進一步加入規則 3 來簡化 MAP 演算法則我們可以得到 Max-Log-MAP 演算法。

 

 

 

 

 

 

 

 

 

Forney Factor Graph – Simple Error Correction Code

by allenlu2007

Loeliger 在 IEEE Signal Processing Magazine 發表的 tutorial paper, 介紹 factor graph 和應用,非常好的文章。

本文節錄其中 error correction code example.

 

先以 Hamming code (n, k) 為例

NewImage 

H 是 parity-check matrix (n-k)xn.  G 是 generator matrix kxn.

NewImage

For Hamming code (7, 4)

H is (7-4)x7 = 3×7

G is 4×7

NewImage

 

Example 1(如何 model): Hamming code (7, 4)

NewImage

Ic 就是 contraints, 必須全部滿足。意即可分解為 products 如 (17).

NewImage

 

接下來考慮 channel/noise 的影響

NewImage

最簡單的 channel 是 memoryless channel

NewImage

或者複雜一點的 state-space channel model (with memory and crosstalk?)

NewImage

NewImage

 

Example 2 (如何 decode):  (4, 2) toy code

我們用一個更簡單例子來 demo 如何用 message passing (belief propagation) decode

假設一個 (4, 2) (非 Hamming) toy code.  

delta(x1+x2+x3).delta(x3-x4)   .. (1)

或 delta(x1+x2+x3).delta(x3+x4)  .. (2)

注意 x3 = x4 相當於 x3+x4 = 0

下圖(a)是 equation (1) 的 FFG, 尚未包含 probability error.

 

下圖 (b) 則是加上 memoryless channel noise model (example 1).

NewImage

當然 error probability 0.1 (10%) depends on SNR.  在此我們假定為 0.1.  

 NewImage

 

所有 y =  (y1, y2, y3, y4) 為已知。

我們最後的目的是計算 marginal pdf : 

NewImage

下圖 Step (c) to (f) 解釋了 message passing 的計算方式。

Step (c) 是 from factor (f1, f2, f3, f4) to variables (x1, x2, x3, x4) pass messages based on equation (30).

FFG forward message passing always starts from leaf, 這一步很容易。

NewImage

Step (d) 再來就是 factor (+, =) 之間的 message passing.  需要用到前文 sum-product message passing 公式 

NewImage

Luckily Loeliger 已經整理兩個 table 如下。Table 1 是 sum-product; Table 2 是 max-product message passing for “+” and “=” factors.   從 Table 1 很容易可以得到 step (d) 的結果 (0.5, 0.5) and (0.82, 0.18)

NewImage

NewImage

 

Step (e) 開始算 backward message passing, 同樣用 Table 1 公式。

For x1 and x2, input is (0.5, 0.5) and (0.9, 0.1) and the result (0.5*0.1+0.5*0.9, 0.5*0.9+0.5*0.1)=(0.5,0.5).

For x3, input is (0.82, 0.18) and (0.9, 0.1) and the result (0.82*0.9, 0.18*0.1) normalize to (0.976,0.024)

same for x4, (0.82*0.1, 0.18*0.9) normalize to (0.336, 0.664)

 

Step (f) 算 x1, x2, x3, x4 的 marginal pdf.  只要把所有的 message multiply and normalize 如下圖。

NewImage

從 sum-product 結果來看 (x1, x2, x3, x4) decode 為 (0, 0, 0, 0) given Y=(0, 0, 1, 0).

再觀察 message passing 的 physical meaning:

x1(0), x2(0) belief x3 應該為 0 (0.82, 0.18); 但 x3(1) 受到 x4(0) 的影響,只能送給 x1, x2 中立 (0.5, 0.5) 的 message.  最後的結果就是 x3 只能接受 x1, x2 的 belief.  

Q: 如果再 iterative 一次,應該會更收歛到 (0, 0, 0, 0)?

A: Based on the python code in Appendix, 大約三次就 marginal probability 完全收歛 1 (0, 0, 0, 0)

Error ~ p^2, 和牛頓法相當。應該是 message ~ u->(y)*u<-(y) 的關係? (上圖)

1’s iteration

Prob(x1(0)) = 0.900000
Prob(x1(1)) = 0.100000
Prob(x2(0)) = 0.900000
Prob(x2(1)) = 0.100000
Prob(x3(0)) = 0.820000
Prob(x3(1)) = 0.180000
Prob(x4(0)) = 0.820000
Prob(x4(1)) = 0.180000

2’s iteration

Prob(x1(0)) = 0.982699
Prob(x1(1)) = 0.017301
Prob(x2(0)) = 0.982699
Prob(x2(1)) = 0.017301
Prob(x3(0)) = 0.989533
Prob(x3(1)) = 0.010467
Prob(x4(0)) = 0.989533
Prob(x4(1)) = 0.010467

3’s iteration
Prob(x1(0)) = 0.999688
Prob(x1(1)) = 0.000312
Prob(x2(0)) = 0.999688
Prob(x2(1)) = 0.000312
Prob(x3(0)) = 0.999996
Prob(x3(1)) = 0.000004
Prob(x4(0)) = 0.999996
Prob(x4(1)) = 0.000004

 

Q: 如果 Y = (0, 0, 1, 1) 如下圖, 結果會變如何?  

NewImage 

Step (d):  x3, x4 “=” to “+” 的 message 從 (0.5, 0.5) 變為 (0.01, 0.99).  Insight 是 (x3, x4) 非常可能是 (1,1)

反之,從 x1, x2 來的 message 仍為 (0.82, 0.18).  Insight 是 (x3, x4) 非常可能是 (0, 0)

Step (e): 此時 backward message from “+” to x1 and x2 為 (0.9*0.01+0.1*0.99, 0.9*0.99+0.1*0.01) = (0.11, 0.89)

backward message from “=” to x3 and x4 為 (0.82*0.1, 0.18*0.9) normalize to (0.336, 0.664)

Step (f): 最後 x1 and x2 are (0.11*0.9, 0.89*0.1) normalize to (0.53, 0.47)

x3 and x4 (0.1*0.336, 0.9*0.664) normalize to (0.05, 0.95)

有趣的是 (x1, x2, x3, x4) decode 為 (0, 0, 1, 1), 並非一個正確的 code word!

Q: 如果再把 (f) 的結果 iterative 一次,可以看到 (x3, x4) 會愈來愈收歛到  (1, 1), 但 (x1, x2) 因為 symmetry 的原故會 deocd 為 (0, 0) or (1, 1) !   Again, (x1, x2, x3, x4) = (0, 0, 1, 1) or (1, 1, 1, 1) 都不是正確的 codeword!!

A: Based on the python code in the Appendix

1’s iteration (0, 0, 1, 1)…
Prob(x1(0)) = 0.525974
Prob(x1(1)) = 0.474026
Prob(x2(0)) = 0.525974
Prob(x2(1)) = 0.474026
Prob(x3(0)) = 0.053247
Prob(x3(1)) = 0.946753
Prob(x4(0)) = 0.053247
Prob(x4(1)) = 0.946753

2’s iteration (0, 0, 1, 1)…
Prob(x1(0)) = 0.500164
Prob(x1(1)) = 0.499836
Prob(x2(0)) = 0.500164
Prob(x2(1)) = 0.499836
Prob(x3(0)) = 0.003170
Prob(x3(1)) = 0.996830
Prob(x4(0)) = 0.003170
Prob(x4(1)) = 0.996830

3’s iteration (0, 0, 1, 1) or (1, 1, 0, 0) (undetermined)
Prob(x1(0)) = 0.500000
Prob(x1(1)) = 0.500000
Prob(x2(0)) = 0.500000
Prob(x2(1)) = 0.500000
Prob(x3(0)) = 0.000010
Prob(x3(1)) = 0.999990
Prob(x4(0)) = 0.000010
Prob(x4(1)) = 0.999990

下圖 summarize 結果。中間的曲線是 x1 and x2; 上下曲線是 x3 and x4.  Iteration 再多都 (x1, x2) 都收歛到 (0.5, 0.5) 而不會收歛到 (0, 1) or (1, 0).

NewImage

 

Q: 如果 (x1, x2) slightly asymmetry (perturbation), 會 decode 為 (1, 0, 1, 1) or (0, 1, 1, 1) 嗎?

A:  YES!  Make (x1, x2) initial error probability from (0.1, 0.1) to (0.1, 0.1001); 只加上小小的 perturbation (0.0001), 如下兩圖 (linear and log scale) 經過 6-7 次的 iteration, (x1, x2) 開始從 (0.5, 0.5) 分開到 (1, 0) or (0, 1).

NewImage

NewImage

 

Summary:

1. Sum-product 並不保証一次收歛到正確的 codeword. 特別 initial probability 如果沒設好情況下。 

2. Iterative decode 也不保証最後收歛到正確的 codeword.

3. Symmetry 似乎會造成收歛到 saddle point.  然而如果加入一些 random perturbation on initial error probability, 可以克服 symmetry 造成的問題。

Q: Max-prodcut 會不同嗎?

A: 下兩圖是加入小 perturbation 的收歛圖,看來和 sum product 幾乎一樣。其他所有的測試也差不多。Max product and sum product 的 tradeoff?  在 Gaussian message passing, max product and sum product 是相同的。

NewImage

NewImage

 

需要用 Viterbi decoding (sequence decode?) to make sure 收歛到 valid codeword?

Proof table 1 and 2?

Factor from edge, why forward message is like the above derivation? 

 

Appendix:  Python code for (4, 2) Toy Code (tab disappears in editor)

#!/usr/bin/env python

## comment

# — +(fp) — =(fe) —
# | | | |
# x1 x2 x3 x4
# f1 f2 f3 f4
# y1 y2 y3 y4

# (y1, y2, y3, y4) are binary observables
# based on Loeliger’s paper
# http://www.ee.cityu.edu.hk/~liping/Research/Journal/Factor%20Graph.pdf

## the following two are for plotting, comment out if not plotting
import numpy as np
import matplotlib.pyplot as plt

plotlog = True
sum_product = False

## input
y1 = 0
y2 = 0
y3 = 1
y4 = 0

## initial probability of binary channel error probability
p = 0.1
#p1 = p+0.0001
p1 = p
p2 = p
p3 = p
p4 = p

## iterative number
count = 3

## array variable for plotting
x1_0_arr = []
x1_1_arr = []
x2_0_arr = []
x2_1_arr = []
x3_0_arr = []
x3_1_arr = []
x4_0_arr = []
x4_1_arr = []

## initial forward message passing based on (y1,y2,y3,y4) and p
if y1 == 0:
x1_fp_0 = 1-p1
x1_fp_1 = p1
else:
x1_fp_0 = p1
x1_fp_1 = 1-p1

if y2 == 0:
x2_fp_0 = 1-p2
x2_fp_1 = p2
else:
x2_fp_0 = p2
x2_fp_1 = 1-p2

if y3 == 0:
x3_fe_0 = 1-p3
x3_fe_1 = p3
else:
x3_fe_0 = p3
x3_fe_1 = 1-p3

if y4 == 0:
x4_fe_0 = 1-p4
x4_fe_1 = p4
else:
x4_fe_0 = p4
x4_fe_1 = 1-p4

for i in range(count):

if i > 0:
print “\n\nIterative %d …\n” %(i)

print “\nForward message passing …”
print “Prob(f<-x1(0)) = %f” %(x1_fp_0)
print “Prob(f<-x1(1)) = %f” %(x1_fp_1)
print “Prob(f<-x2(0)) = %f” %(x2_fp_0)
print “Prob(f<-x2(1)) = %f” %(x2_fp_1)
print “Prob(f<-x3(0)) = %f” %(x3_fe_0)
print “Prob(f<-x3(1)) = %f” %(x3_fe_1)
print “Prob(f<-x4(0)) = %f” %(x4_fe_0)
print “Prob(f<-x4(1)) = %f” %(x4_fe_1)

###### Step (d)
## sum product message fp to fe
if sum_product:
fp_fe_0_p = x1_fp_0 * x2_fp_0 + x1_fp_1 * x2_fp_1
fp_fe_1_p = x1_fp_0 * x2_fp_1 + x1_fp_1 * x2_fp_0
fp_fe_0 = fp_fe_0_p / (fp_fe_0_p + fp_fe_1_p)
fp_fe_1 = fp_fe_1_p / (fp_fe_0_p + fp_fe_1_p)
else: # max_product
fp_fe_0_p = max(x1_fp_0 * x2_fp_0, x1_fp_1 * x2_fp_1)
fp_fe_1_p = max(x1_fp_0 * x2_fp_1, x1_fp_1 * x2_fp_0)
fp_fe_0 = fp_fe_0_p / (fp_fe_0_p + fp_fe_1_p)
fp_fe_1 = fp_fe_1_p / (fp_fe_0_p + fp_fe_1_p)

## sum product message fe to fp
fe_fp_0_p = x3_fe_0 * x4_fe_0
fe_fp_1_p = x3_fe_1 * x4_fe_1
fe_fp_0 = fe_fp_0_p / (fe_fp_0_p + fe_fp_1_p)
fe_fp_1 = fe_fp_1_p / (fe_fp_0_p + fe_fp_1_p)

print “\fp and fn message passing …”
print “Prob(fp->fe(0)) = %f” %(fp_fe_0)
print “Prob(fp->fe(1)) = %f” %(fp_fe_1)
print “Prob(fe->fp(0)) = %f” %(fe_fp_0)
print “Prob(fe->fp(1)) = %f” %(fe_fp_1)

###### Step (e)
# x1 backward message
if sum_product:
fp_x1_0_p = fe_fp_0 * x2_fp_0 + fe_fp_1 * x2_fp_1
fp_x1_1_p = fe_fp_1 * x2_fp_0 + fe_fp_0 * x2_fp_1
fp_x1_0 = fp_x1_0_p / (fp_x1_0_p + fp_x1_1_p)
fp_x1_1 = fp_x1_1_p / (fp_x1_0_p + fp_x1_1_p)
else: # max product
fp_x1_0_p = max(fe_fp_0 * x2_fp_0, fe_fp_1 * x2_fp_1)
fp_x1_1_p = max(fe_fp_1 * x2_fp_0, fe_fp_0 * x2_fp_1)
fp_x1_0 = fp_x1_0_p / (fp_x1_0_p + fp_x1_1_p)
fp_x1_1 = fp_x1_1_p / (fp_x1_0_p + fp_x1_1_p)

# x2 backward message
if sum_product:
fp_x2_0_p = fe_fp_0 * x1_fp_0 + fe_fp_1 * x1_fp_1
fp_x2_1_p = fe_fp_1 * x1_fp_0 + fe_fp_0 * x1_fp_1
fp_x2_0 = fp_x2_0_p / (fp_x2_0_p + fp_x2_1_p)
fp_x2_1 = fp_x2_1_p / (fp_x2_0_p + fp_x2_1_p)
else: # max product
fp_x2_0_p = max(fe_fp_0 * x1_fp_0, fe_fp_1 * x1_fp_1)
fp_x2_1_p = max(fe_fp_1 * x1_fp_0, fe_fp_0 * x1_fp_1)
fp_x2_0 = fp_x2_0_p / (fp_x2_0_p + fp_x2_1_p)
fp_x2_1 = fp_x2_1_p / (fp_x2_0_p + fp_x2_1_p)

# x3 backward message
fe_x3_0_p = fp_fe_0 * x4_fe_0
fe_x3_1_p = fp_fe_1 * x4_fe_1
fe_x3_0 = fe_x3_0_p / (fe_x3_0_p + fe_x3_1_p)
fe_x3_1 = fe_x3_1_p / (fe_x3_0_p + fe_x3_1_p)

# x4 backward message
fe_x4_0_p = fp_fe_0 * x3_fe_0
fe_x4_1_p = fp_fe_1 * x3_fe_1
fe_x4_0 = fe_x4_0_p / (fe_x4_0_p + fe_x4_1_p)
fe_x4_1 = fe_x4_1_p / (fe_x4_0_p + fe_x4_1_p)

print “\nBackward message passing …”
print “Prob(f->x1(0)) = %f” %(fp_x1_0)
print “Prob(f->x1(1)) = %f” %(fp_x1_1)
print “Prob(f->x2(0)) = %f” %(fp_x2_0)
print “Prob(f->x2(1)) = %f” %(fp_x2_1)
print “Prob(f->x3(0)) = %f” %(fe_x3_0)
print “Prob(f->x3(1)) = %f” %(fe_x3_1)
print “Prob(f->x4(0)) = %f” %(fe_x4_0)
print “Prob(f->x4(1)) = %f” %(fe_x4_1)

#### Step (e) compute marginal probability
x1_0_p = x1_fp_0 * fp_x1_0
x1_1_p = x1_fp_1 * fp_x1_1
x1_0 = x1_0_p / (x1_0_p + x1_1_p)
x1_1 = x1_1_p / (x1_0_p + x1_1_p)

x2_0_p = x2_fp_0 * fp_x2_0
x2_1_p = x2_fp_1 * fp_x2_1
x2_0 = x2_0_p / (x2_0_p + x2_1_p)
x2_1 = x2_1_p / (x2_0_p + x2_1_p)

x3_0_p = x3_fe_0 * fe_x3_0
x3_1_p = x3_fe_1 * fe_x3_1
x3_0 = x3_0_p / (x3_0_p + x3_1_p)
x3_1 = x3_1_p / (x3_0_p + x3_1_p)

x4_0_p = x4_fe_0 * fe_x4_0
x4_1_p = x4_fe_1 * fe_x4_1
x4_0 = x4_0_p / (x4_0_p + x4_1_p)
x4_1 = x4_1_p / (x4_0_p + x4_1_p)

print “\nFinal marginal probability …”
print “Prob(x1(0)) = %f” %(x1_0)
print “Prob(x1(1)) = %f” %(x1_1)
print “Prob(x2(0)) = %f” %(x2_0)
print “Prob(x2(1)) = %f” %(x2_1)
print “Prob(x3(0)) = %f” %(x3_0)
print “Prob(x3(1)) = %f” %(x3_1)
print “Prob(x4(0)) = %f” %(x4_0)
print “Prob(x4(1)) = %f” %(x4_1)

x1_fp_0 = x1_0
x1_fp_1 = x1_1
x2_fp_0 = x2_0
x2_fp_1 = x2_1
x3_fe_0 = x3_0
x3_fe_1 = x3_1
x4_fe_0 = x4_0
x4_fe_1 = x4_1

x1_0_arr.append(x1_0)
x1_1_arr.append(x1_1)
x2_0_arr.append(x2_0)
x2_1_arr.append(x2_1)
x3_0_arr.append(x3_0)
x3_1_arr.append(x3_1)
x4_0_arr.append(x4_0)
x4_1_arr.append(x4_1)

print x1_0_arr
print x1_1_arr
print x2_0_arr
print x2_1_arr
print x3_0_arr
print x3_1_arr
print x4_0_arr
print x4_1_arr

if plotlog:
plt.semilogy(range(count), x1_0_arr, range(count), x1_1_arr, marker=”+”)
plt.semilogy(range(count), x2_0_arr, range(count), x2_1_arr, marker=”*”)
plt.semilogy(range(count), x3_0_arr, range(count), x3_1_arr, marker=”^”)
plt.semilogy(range(count), x4_0_arr, range(count), x4_1_arr, marker=”o”)
plt.axis([0, count, 1e-3, 1])
else:
plt.plot(range(count), x1_0_arr, range(count), x1_1_arr, marker=”+”)
plt.plot(range(count), x2_0_arr, range(count), x2_1_arr, marker=”*”)
plt.plot(range(count), x3_0_arr, range(count), x3_1_arr, marker=”^”)
plt.plot(range(count), x4_0_arr, range(count), x4_1_arr, marker=”o”)
plt.xlabel(‘Iteration number’)
plt.ylabel(‘Marginal probablity’)
plt.grid()
plt.show()

 

Graph Model Message Passing vs. Nerual Net BackProp vs. Variational Inference

by allenlu2007

前文介紹 Gaussian message.   主要的缺點是如果通過某一些 block (非 +, =, gain) Gaussian 就會被破壞。

解決的方法是 variational inference, 其實就是 approxmiation, 假設還是 Gaussian.

但在 neural network 也有 forward path 和 error backward propagation.  似乎更 general.

也是只 pass variational message (dE/dW, or dE/dX) back prop.  但只用 1st order statistics.  因為經過 nonlinear block (e.g. sigmoid), cost function 已經非 convex.  即使用 2nd order statistics 也沒用?

 

Gaussian message –> 2nd order statistics,  fast but limited use (not Gaussian)

Variation inference –> 2nd order or 1st order statistics? 2nd order

Neural back prop –> 1st order statistics

 

–> 這些觀念應該可以 combine 在一起。雖然一個是 Bayesian, 另一個是 neural network.

 

 

前文介紹 Forney factor graph (FFG) 的基本型。FFG 可以用在更複雜的 case.

介紹完 FFG 進階型之後,討論最重要的 continuous message: Gaussian message.

主要是 Gaussian message 可以得到 close form, 同時可用為其他 continuous message 的近似。

再來應用 Gaussian message 在 Kalman filter.

 

FFG 進階型  

1. variable X 出現在兩個以上的 factor, 可以用 “=” factor, 代表 delta function.

NewImage

 

2. 可以 generalize “=” factor to block diagram factor 如下

NewImage

 Let  U = x3,  X = x3′, Z = x3”;  g(u,w) = U;  h(X,Y) = X 就可化簡為 1.

 

Discrete and Continuous

Discrete message passing (from factor to variable)

NewImage

Continuous message passing (from factor to variable)

NewImage

Special case, if Y = h(X1) is a delta function;  u(y) = u(X1) as expected

 

對於 continuous message passing 的幾種近似法

NewImage

 

Gaussian Message (reference)

故名思義,Gaussian message 的 message 的型式就是 Gaussian distribution 如下。

NewImage

 

接下來就是 Gaussian message 在 sum-product 以積分代替 summation:

NewImage

 

  • Gaussian message 的 products 還是 Gaussian message
  • (7) 的 integration (marginalization) 也是 Gaussian form
  • 更進一步,上述的 u(x1), u(x2), .. u(xm) 以及 h(y, x1…, xm) 如果是 Gaussian form, u(y) 也是 Gaussian message!!!

先考慮簡單的 case if q(x, y) is the quadratic form: 

NewImage

  • 注意原來的 integration 變成 minimization !!
更重要的是 delta function 可以變為 Gaussian function 如下:
NewImage

因此原來的

NewImage

or 

NewImage 

就變成:

 NewImage

同時比較 max-product rule:

NewImage 

如果f, u1, … un are Gaussian,  (16) 和上述公式會得到同樣的結果t!!

  • 如果 factor 和所有 input message 都是 Gaussian form, output message 也是 Gaussian form!
  • 以上,sum-product 和 min-sum 在 Gaussian message 和 Gaussian factor 是等價!!
  • 所謂的 Gaussian factor 最簡單就是 linear combination, (“=”, “+”, scaling, etc.) 可以用 delta function 表示。delta function 可以用 Gaussian function 表示。
  • 如果 factor 不是 Gaussian form, output message 就不是 Gaussian message (e.g. Extended Kalman filter, Unscented Kalman filter).  當然 sum-product 和 max-product 的結果也不同。
  • Exponential family message passing (variational message passing): if input message and factor are exponential family, output message is also exponential family?  (Yes?)  sum-product and max-product are equivalent? (No?)

 

Example

Equality constraint 如下圖 1:

NewImage 

非常漂亮簡潔的推導!  Brute force 推導如下下圖,非常 dirty 也沒有太多 insight. 

 

NewImage

NewImage 

 

NewImage

NewImage

其他的推導請參考 reference.   Some composite blocks 結果如下

 NewImage

Rule 6: Vx = A-1Vy A-H 是錯的。正確如下:

NewImage

NewImage

NewImage

Rule 5 是 Rule 6 duality. 

NewImage

 

In more details (Gaussian assumption may be relaxed for some rules)

NewImage 

NewImage

NewImage

NewImage

 

Kalman Filter (Gaussian Message + HMM)

Kalman filter 的 model 如下 (similar to HMM)

NewImage

重新看 Kalman filter 和以上公式的關係

NewImage

 

正好用到 Rule 3/4, Rule 6, and Rule5.   “+” 相當 variance/covariance 串聯; “=” 相當 variance/covariance 並聯。對於 variance/covariance^(-1) 剛好相反。

更細來說 Rule 3/4 和 Rule 6 對所有 message (distribution) 都適用 (A is invertible); 只有 Rule 5 需要 Gaussian message.  Practically speaking, 當有夠多的 Uk and Wk, Xk 最後趨近 Gaussian message (central limit theorem).

推導 Kalman Filter

Step 1:  Entry6 的 output, denoted Xk_hat (Kalman filter 的 estimation stage) and assuming mean Uk and covariance = Q, b = 1

          Xk_hat 的 forward Gaussian message

          mean = A m(k-1) + Uk

          covariance = A V(k-1) A^H + Q

Step 2: Entry 5 的 output Xk (Kalman filter 的 update stage) and assuming Wk zero mean and covariance R.

Xk 的 forward Gaussian message  如下

NewImage

 

Eq. (35) 即是 Kalman filter 中的 mean.  (Xk 的 mean) = (xk_hat 的 mean) + (Kalman gain) * innovation

 innovation = (observation – prediction) 

 Kalman gain = Vinx A^H ( VinY + AVinxA^H)-1 = Vinx A^H (R + A Vinx A^H)-1

 

Eq. (36) 則是 Kalman filter 中的 covariance

Voutz = (I – Vinx A^H G A) Vinx = ( I – Kg * A) (xk_hat 的 covariance)

 

In summary:  Eq (35) 和 Eq (36) 就是熟知的 Kalman filter = Forward Gaussian message

 

推導 Kalman Smoother

Kalman filter:  就是找 Gaussian conditional pdf: P(xk | y1, y2, …yk) given current observables (real-time)

Kalman smoother: 就是找 Gaussian conditional pdf: P(xk | y1, y2, … yT) given all observables (non real-time)

 

Q: Kalman filter 是 forward message.  那 Kalman smoother 是 backward message? 

A: No.  smoother 是 (forward message) * (backward message), i.e. product of forward and backward message.

也就是 P(x(k)) | all observed y)

 

如何推導? TBD.  以下是 high level 的說明,不是証明。

Kalman filter 是 alpha (forward message); Kalman smoother 是 alpha * beta, i.e. (forward message)*(backward message).  如果真的用 (forward message)*(backward message) 來推導 Kalman smoother 的 close form.  應該會非常 messy.  Reference 只用 forward message 來推導 Kalman filter.  同樣方法應可用在 Kalman smoother.

實際上 Kalman smoother 是用 conditional PDF 來推導 recursive close form,參考 reference.

 

下圖是用 HMM 來看 Kalman smoother (Z 是 state variable, X 是 observable)

f(zn, zn-1) 是濃縮的 factor = (transition probability) * (emission probability)

NewImage

NewImage

NewImage 

NewImage

 

Q&A

Q: 在 Gaussian message 時,Marginalization (sum product) and Maximization (max product, MAP) 得到的 message 剛好相同 

A: Yes. 但必須包含 factor.  Factor 也必須是 Gaussian form, 如 delta function 或 linear combination.

 

Q: 在 Kalman filter 中,sum-product 代表 P(xk | y1, … yk) 的 mean and variance (minimum mean square error); max-product 則是 ML estimator (most likely).  這兩個結果一樣。HMM 有同樣結果嗎?

A:  No,  HMM 非 Gaussian message.  HMM 是 discrete PMF.  Sum-product algorithm 對應 forward-backward message (alpha-beta).   類似的 Baum-Welch algorithm belongs to sum product algorithm (EM algorithm)!

Max-product algorithm 在 HMM 對應的是 Viterbi family algorithm (VA, SOVA, hard-decision VA).  

 

Viterbi in Gaussian Messge

Q: 在 HMM model, most probable sequence 和 most probable values 是不同的。

most probable sequence: Viterbi algorithm

most probable value: forward-backward algorithm (alpha * beta)

在 Gaussian message (continuous HMM) 也有類似的 algorithm 嗎?


A: No.  在 Continuous HMM Gaussian message both 

most probable value: forward-backward algorithm (i.e. Kalman smoother)

most probable sequence: forward-backward algorithm (i.e. Kalman smoother)

原因 sum-product (forward-backward) and max-product (VA) 在 Gaussian message 等價!!

 

NewImage 

Or 

 NewImage

 

 

NewImage