Edited at

Belief Propergation(確率伝播法) をpythonで。ついでにQuantum LDPC-code


. Belief Propergation(確率伝播法) をpythonで

ここでは、BPをやってみようと思います。


0. Belief Propergation(確率伝播法)

概念的には機械学習やら画像処理、符号理論などでよく使われる手法です。

多分色々役に立つよ!


1. お願い

プログラムがすごい汚いのでアドバイスあったらください...

命名規則とかもぐちゃぐちゃなので...

データ構造とか、Pythonのスマートな書き方とかそういう基本的なことを教えてください...


2. 何をする

BPを使って量子誤り訂正符号の通信路情報を更新して復号する。

基礎的な量子誤り訂正符号の復号に関しては前の記事を。

量子誤り訂正符号の復号に関するプログラムは今回本質ではないのでいつかgithubにあげときます。


3. プログラム解説


3.1 補題的な関数の準備やらインポートやら

まずnumpyをインポートする。

import numpy as np

次に同時分布の計算と正規化の計算を行う関数。

正規化、時間かかるからあんまやりたくないんだけど、しないとメッセージが途中桁落ちしちゃう...

#任意の同時分布を返す。

def jointDistribution(prob):
result=np.ones(1)
reshape=()
for p in prob:
result=np.kron(result,p)
reshape=reshape+(len(p),)
return result.reshape(reshape)

#正規化関数
def normalize(Prob):
nProb=np.ones(len(Prob),dtype='float128')
sum=np.array(Prob).sum()
if sum==0:
return np.array([1/len(Prob)]*len(Prob))
return Prob/sum


3.2 Nodeクラス

Nodeクラスは、変数ノードのデータ構造とメッセージの計算を行う。関数ノードと変数ノードで分けようと思ったけど、なんか関数いじれば同じ感じにできそうだったので一緒にしてしまった。

メッセージ計算の処理が違うし、なんかもっとスマートにできる気がするので、

class Node:

def __init__(self,id,bit,func):
self.id = id
self.neighbor = {}
self.message = {}
self.bit = bit
self.func = func
def addNeighbor(self, node):
self.neighbor[node.id]=node
def getMessage(self,key):
return self.message[key]
def setMessage(self,message,key):
if message.shape!=self.message[key].shape:
print("メッセージの型が違います")
print("セットしたいメッセージの型",message.shape)
print("メッセージの型",self.message[key].shape)
exit()
self.message[key]=message
def initializeMessage(self):
for key in self.neighbor:
bit= self.bit if self.neighbor[key].bit==None else self.neighbor[key].bit
self.message[key] = np.ones(bit)

def calc(self,target):
prob=[]
number=()
k=0
for key in self.neighbor:
if key!=target:
prob=prob+[self.neighbor[key].message[self.id]]
number=number+(k,)
else:
bit= self.bit if self.neighbor[key].bit==None else self.neighbor[key].bit
prob=prob+[np.ones(bit)]
k+=1

#関数ノードだったら
if None==self.bit:
prob=jointDistribution(prob)
self.message[target]=normalize((self.func*prob).sum(number))

#変数ノードだったら
else:
p=np.ones(self.bit)
for i in range(len(prob)):
if i in number:
p*=prob[i]
self.message[target]=normalize(p)

def calcAll(self):
for key in self.neighbor:
self.calc(key)


3.2.1 変数や関数の説明

・self.id

ノードのidを表す。文字列にしてる。絶対被らないようにしてください。チェックする機能ない

・self.neighbor

隣接するノードを表す。

・self.message

隣接するノードに対するメッセージを表す。

・self.bit

変数ノードのbitを表す。関数ノードの場合はNoneにする。

・self.func

関数を行列表現で渡す。全パターンの出力を行列にして渡す。ポテンシャルともいうっぽい。

・addNeighbor(Node node)

自分のノードに隣接するノードを加える。例えばグラフ構造が木になっており、単方向にしかメッセージを交換しないのであれば、自分のノードに加えるだけで良い。ただloopy BPなど、グラフ構造が木でないとき、両方のグラフにaddNeighborで加える必要がある。

なお、ちゃんと接続できているかチェックする機能はない...

・getMessage(String key)

ノードのIDに応じて、メッセージを取得する。

・setMessage(float *message,String key)

メッセージをセットする。ちょっと例外処理などを行う時用。

・initializeMessage()

メッセージを一様分布で初期化

・calc(String key)

id方向へのメッセージを計算

・calcAll()

全方向のメッセージを計算


3.3 factorGraph クラス

隣接行列を与えてグラフを一気に作ったり、メッセージを計算する順番を与えたりする管理クラス。

クラスの使い方あってるのかな...?

Javaを授業でやった時に色々「クラス」の使い方を勉強したけどもう何も覚えてない...

class factorGraph:

def __init__(self):
self.node={}
def adjacencyMatrixToFactorGraph(self,v,f,matrix,vbit,ffunc):
for i in range(len(v)):
for j in range(len(f)):
if(matrix[i][j]==1):
if not v[i] in self.node:
self.node[v[i]]=Node(v[i],vbit[i],None)
if not f[j] in self.node:
self.node[f[j]]=Node(f[j],None,ffunc[j])
self.node[v[i]].addNeighbor(self.node[f[j]])
self.node[f[j]].addNeighbor(self.node[v[i]])
for key in self.node:
self.node[key].initializeMessage()
def calcSequential(self,message,uniqueProcess):
for m in message:
if m[1]==None:
self.node[m[0]].calcAll()
else:
self.node[m[0]].calc(m[1])
if not uniqueProcess is None:
#例外処理をする場合
for u in uniqueProcess:
if uniqueProcess[0]==m[0] and uniqueProcess[1]==m[1]:
self.node[m[0]].setMessage(uniqueProcess[2](self.node[m[0]].getMessage(m[1])),m[1])

def getMessage(self,id1,id2):
return self.node[id1].getMessage(id2)


3.3.1 変数や関数の説明

・self.node

nodeを格納する。ノードのIDがkey

・adjacencyMatrixToFactorGraph(v,f,matrix,vbit,ffunc)

・v,f

それぞれ、nodeのIDをもつ。


vとfの例

v=['x1','x2','x3']

f=['f12','f13']

さらにmatrixは、横が関数ノードの番号、縦が変数ノードの番号で

$$\left(\begin{array}{ll}

1&1\\

1&0\\

0&1\\

\end{array}\right)$$

と書く。したがって


隣接行列の例

matrix=np.array([[1,1],[1,0],[0,1]])


となる。

・vbit

vbitは、各変数ノードの確率の要素数がいくらあるかを示している。

というかbitって書いてしまった!やってしまった。

・ffunc

ffuncは、関数ノードの関数を表している。

$$f_{12}:(x_1,x_2)$$

は、出力は、全てのペアに対して行列で与えられる。

もし$x_1$と$x_2$が2つの要素をもつ確率(1 bit)で、等しい時1を表す関数

$$\delta_{x_1,x_2}=\left\{\begin{array}{ll}

1 & ({\rm if\ } x_1=x_2)\\

0 & ({\rm if\ } x_1\neq x_2)

\end{array}\right.$$

を行列で書くと


関数の例

f12matrix=np.zeros((2,2))

f12matrix[0][0]=1
f12matrix[1][1]=1

となる。


4. 量子LDPC符号を復号してみる。

シチュエーション

復号対象:STEANE符号

最尤復号+確率伝播法で確率更新


4.1結果

Figure_1.png

悪くなったけどSTEANE符号はループがたくさんあるからかな(ループが多くあるとうまくいかない)


5. プログラム

BPを行うためのプログラム


bp.py

import numpy as np

#任意の同時分布を返す。
def jointDistribution(prob):
result=np.ones(1)
reshape=()
for p in prob:
result=np.kron(result,p)
reshape=reshape+(len(p),)
return result.reshape(reshape)

#正規化関数
def normalize(Prob):
nProb=np.ones(len(Prob),dtype='float128')
sum=np.array(Prob).sum()
if sum==0:
return np.array([1/len(Prob)]*len(Prob))
return Prob/sum

class Node:
def __init__(self,id,bit,func):
self.id = id
self.neighbor = {}
self.message = {}
self.bit = bit
self.func = func
def addNeighbor(self, node):
self.neighbor[node.id]=node
def getMessage(self,key):
return self.message[key]
def setMessage(self,message,key):
if message.shape!=self.message[key].shape:
print("メッセージの型が違います")
print("セットしたいメッセージの型",message.shape)
print("メッセージの型",self.message[key].shape)
exit()
self.message[key]=message
def initializeMessage(self):
for key in self.neighbor:
bit= self.bit if self.neighbor[key].bit==None else self.neighbor[key].bit
self.message[key] = np.ones(bit)

def calc(self,target):
prob=[]
number=()
k=0
for key in self.neighbor:
if key!=target:
prob=prob+[self.neighbor[key].message[self.id]]
number=number+(k,)
else:
bit= self.bit if self.neighbor[key].bit==None else self.neighbor[key].bit
prob=prob+[np.ones(bit)]
k+=1

#関数ノードだったら
if None==self.bit:
prob=jointDistribution(prob)
self.message[target]=normalize((self.func*prob).sum(number))

#変数ノードだったら
else:
p=np.ones(self.bit)
for i in range(len(prob)):
if i in number:
p*=prob[i]
self.message[target]=normalize(p)

def calcAll(self):
for key in self.neighbor:
self.calc(key)

class factorGraph:
def __init__(self):
self.node={}
def adjacencyMatrixToFactorGraph(self,v,f,matrix,vbit,ffunc):
for i in range(len(v)):
for j in range(len(f)):
if(matrix[i][j]==1):
if not v[i] in self.node:
self.node[v[i]]=Node(v[i],vbit[i],None)
if not f[j] in self.node:
self.node[f[j]]=Node(f[j],None,ffunc[j])
self.node[v[i]].addNeighbor(self.node[f[j]])
self.node[f[j]].addNeighbor(self.node[v[i]])
for key in self.node:
self.node[key].initializeMessage()
def calcSequential(self,message,uniqueProcess):
for m in message:
if m[1]==None:
self.node[m[0]].calcAll()
else:
self.node[m[0]].calc(m[1])
if not uniqueProcess is None:
#例外処理をする場合
for u in uniqueProcess:
if uniqueProcess[0]==m[0] and uniqueProcess[1]==m[1]:
self.node[m[0]].setMessage(uniqueProcess[2](self.node[m[0]].getMessage(m[1])),m[1])

def getMessage(self,id1,id2):
return self.node[id1].getMessage(id2)



6.改善点

グラフによっては処理速度がきつくなってくる。

そんな時、様々な近似処理が行えるように改修していきたい。


7.参考文献

確率伝播法とその応用,林和則

確率伝播法を組んでみる(Python)