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()

 

Advertisements