allenlu2007

A great WordPress.com site

Month: August, 2014

EM, Variational Message Passing

by allenlu2007

 

之前 EM algorithm 有幾何解釋,代數解釋。

幾何基本上 based on concave (or convex) function;  代數大多 based on Jensen’s inequality.

雖然可以說明,但缺少一個大的 picture.

 

 

另外的方法是用 variational message passing 來說明。

請參考 Tom and Bishop paper

For Bayesian DAG graph

NewImage 

 

 

 

 

 

 

 

 

 

EM algorithm:  Q(H) = ??    Lower bound, L(Q), 就變成??

NewImage

 

 

Variational message passing:  Q(H)  (i) factorization form; (ii) conjugate exponential model

NewImage

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 – 電路應用

by allenlu2007

 

基本型

NewImage

進階型

所有 linear circuit (Kalman) 都會 preserve Gaussian message 

a. small signal analysis

 

b. large signal but using other approximation technique 

(i) extended Kalman (or extended factor graph model)

(ii) unscented Kalman (or unscented fator graph model)

(iii) MCMC approximation? 

 

 

Application

1. Noise analysis

2. Monte Carlo analysis of mismatch (Cadence DC mismatch?)

   bandgap, DAC, reference?

3. Monte Carlo analysis of yield? (very difficult? depends on performance criterion)

Factor Graph Misc Q&A

by allenlu2007

 

Q: What’s half edge?  Can it convert to normal edge?

A: Half edge 就是只有一邊有 factor node, 另一邊沒有任何 factor node. 

Half edge 可以變為 normal edge, 只要另一邊加上一個 factor node = 1.

 

Q: Variable 如果有兩個以上的 factor nodes, 如何表示和計算?

A: Message passing 分為兩種: variable to factor, and factor to variable

NewImage

Variable to factor message (附圖):  products from all other factor node messages (或者接收者的 message = 1)

Message from factor to variable:  sum of products except the recipient variable.

 

FFG 因為只有 factor node, variable 只有 edge.

如果 variable edge 連到一個 factor node, 稱為 half edge, 另一端的 factor 可以設為 1 (fu(x) =1 )

  –x– f(x, y, z) — y                   ===       (f=1) –x–  f(x, y, z) — y

joint pdf:  P(x, y, z, …) = f(x, y, z) * ..                    P(x, y, z, ..) = 1 * f(x, y, z) * …

v->f msg: u(x->f) = 1  (no, should be uniform dist?) 

f->v msg: u(f->x) = sum{y,z}( f(x, y, z) * u(y->f) * u(z->f)}

 

 

如果 variable edge 連到兩個 factor node, 稱為 full edge, variable to factor message 只是 passing by above definition.

如果 variable edge 連到三個 factor node, f1(x, ..), f2(x, ..), f3(x, ..). message from f1 to x

 

Q: message can be scale?  

A: Yes

 

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 演算法。