Message graph demo

g={1->4,2->5,3->6,4->5,5->6,6->4}//graph2mat//symmetrize

GraphPlot[{{1->4,1},{2->5,1},{3->6,1},{4->5,0},{5->6,0},{6->4,0}},VertexLabeling->True,EdgeLabeling->All,EdgeRenderingFunction->Function[{r,v,lbl},If[lbl==0,{Thickness[.01],Red,Line[r]},{Thickness[.01],Green,Line[r]}]],VertexRenderingFunction->Function[{r,v},{White,Switch[v,2,EdgeForm[{Thickness[.01],Red}],1,EdgeForm[{Thickness[.01],Dashed,Black}],_,EdgeForm[{Thickness[.005],Black}]],Disk[r,.1],Black,Text[v,r]}],VertexCoordinateRules->{1->{(Sqrt[1/2]),(Sqrt[1/2])},3->{-(Sqrt[1/2]),(Sqrt[1/2])},2->{0,-1}}]

Listing all the labelings

coloredPlot[g_,chosen_]:=Module[{},
GraphPlot[g,VertexRenderingFunction->Function[{r,v},{If[MemberQ[chosen,v],Black,Red],Disk[r,If[v==1,.3,.2]]}],EdgeLabeling->None]
]

g1=coloredPlot[g,{1}~Join~#]&/@Subsets[Range[2,6]]//Partition[#,4]&//GraphicsGrid;

g2=coloredPlot[g,{}~Join~#]&/@Subsets[Range[2,6]]//Partition[#,4]&//GraphicsGrid;

GraphicsRow[Rasterize/@{g1,g2},Spacings->{50,0}]

Message passing graph (each node is a message)

messageGraph[g_]:=Module[{},
edges=Position[g,1];
n=Length[edges];
(*Array[Boole[edges[[#1,2]]==edges[[#2,1]]]&,{n,n}]*)
Cases[Table[If[edges[[i,2]]==edges[[j,1]]&&(edges[[i,1]]≠edges[[j,2]]),edges[[i]]->edges[[j]],0],{i,n},{j,n}]//Flatten,Except[0]]
]

(* makes 0-1 matrix symmetric *)
symmetrize[mat_]:=Unitize[mat+Transpose[mat]];

g=graph2mat[{1->4,2->5,3->6,4->5,5->6,6->4}]//symmetrize;
gg=messageGraph[g];
gg//GraphPlot[#,VertexLabeling->True,DirectedEdges->True]&

J=g Array[If[#1>3&&#2>3,-1,1]&,{6,6}];
J//MatrixForm

Loopy BP implementation

loopyBP[J_,θ_]:=Module[{n=Length[J]},Function[{m},Table[ArcTanh[Tanh[J[[i,j]]] Tanh[θ[[i]]+Total[m[[All,i]]]-m[[j,i]]]],{i,1,n},{j,1,n}]]]

θ={0,-1,0,0,0,0};
f=loopyBP[J,θ];
(Total[#,{1}]+θ)&/@NestList[f,Array[0.&,{6,6}],12]//ArrayPlot

((Total[#,{1}]+θ)&/@NestList[f,Array[0.&,{6,6}],30])//Column

v=.;v[i_]:=((2Total[#,{1}]+θ)&/@NestList[f,Array[0.&,{6,6}],20])[[All,i]];
ListLinePlot[v[#]&/@Range[4],PlotRange->All,Mesh->All]

2(Total[Nest[loopyBP[J,θ],ConstantArray[0,{6,6}],30]][[1]]//N)

marginalize[J,θ,1]//N

Checking against exact marginalization on tree

J2={{0,0,0,1,0,0},{0,0,0,0,1,0},{0,0,0,0,0,1},{1,0,0,0,-1,-1},{0,1,0,-1,0,0},{0,0,1,-1,0,0}};
GraphPlot[Unitize[J2],VertexLabeling->True]

Marginalize

marginalize[J_,θ_,var_]:=Module[{n=Length[θ]},
pl=Select[Tuples[{-1,1},n],#[[var]]==1&];
nl=Select[Tuples[{-1,1},n],#[[var]]==-1&];
pot[l_]:=l.J.l/2+l.θ;
Log[Total[Exp[pot[#]]&/@pl]/Total[Exp[pot[#]&/@nl]]]
]

answer1=marginalize[J2,θ,1]//FullSimplify

Loopy BP

answer2=(2Total[Nest[loopyBP[J2,θ],ConstantArray[0,{6,6}],10]])[[1]]//TrigToExp//FullSimplify

Reduce[answer1==answer2]

{answer1//N,answer2//N}

Message passing charts

GraphPlot[{{1->"Xi","branch 1"},{2->"Xi","branch2"},{3->"Xi","branch3"}},DirectedEdges->True,VertexLabeling->True]

LayeredGraphPlot[Table["X1="<>ToString[i]->"Xi="<>ToString[j],{i,5},{j,5}]//Flatten,DirectedEdges->True,VertexLabeling->True,PlotLabel->"Marginalization"]

Final configuration

bpgraph[g_,seen_,actual_]:=Module[{rs=0.20},
GraphPlot[g,VertexLabeling->True,EdgeLabeling->All,VertexRenderingFunction->Function[{r,v},{White,EdgeForm[Black],Rectangle[r-{rs,rs},r+{rs,0}],Rectangle[r+{-rs,0},r+{rs,rs}],Black,Rectangle[r-{rs,rs},r+{-rs+(2rs)seen[[v]],0}],Rectangle[r-{rs,0},r+{-rs+(2rs)actual[[v]],rs}]}]]
]

acts=(((Exp[#]))/(Exp[#]+1))&@marginalize[J,θ,#]&/@Range[6]//N;

θ={0,-1,0,0,0,0};
f=loopyBP[J,θ];
seen=((Exp[2#])/((Exp[2#]+1)))&/@((Total[#,{1}]+θ)&@Nest[f,Array[0.&,{6,6}],200]);

bpgraph[g,acts,seen]

(* draws graph with nodes annotated (0..1) *)
statgraph[g_,seen_]:=Module[{rs=0.20},GraphPlot[g,VertexLabeling->True,EdgeLabeling->All,VertexRenderingFunction->Function[{r,v},{White,EdgeForm[Black],Rectangle[r-{rs,rs},r+{rs,rs}],Black,Rectangle[r-{rs,rs},r+{-rs+(2 rs) seen[[v]],rs}]}]]]

seenArr=((Exp[2#])/((Exp[2#]+1)))&/@((Total[#,{1}]+θ)&/@NestList[f,Array[0.&,{6,6}],20]);
images=statgraph[g,#]&/@seenArr;

BP on grid with random potentials

g=makeChain[5,5]//graph2mat//symmetrize;
n=Length[g];
θ=Array[Boole[#==1]&,25];
J=g RandomReal[{-4,4},{25,25}];
f=loopyBP[J,θ];
Total[NestList[f,ConstantArray[0.,{n,n}],100],{2}]//Transpose//ArrayPlot

g=makeChain[5,5]//graph2mat//symmetrize;
n=Length[g];
θ=Array[Boole[#==1]&,25];
J=g RandomReal[{-4,4},{25,25}];
f=loopyBP[J,θ];
seenArr=((Exp[2#])/((Exp[2#]+1)))&/@((Total[#,{1}]+θ)&/@NestList[f,ConstantArray[0.,{n,n}],100]);
seenArr//Transpose//ArrayPlot

seenArr=((Exp[2#])/((Exp[2#]+1)))&/@((Total[#,{1}]+θ)&/@NestList[f,ConstantArray[0.,{n,n}],300]);
seenArr//Transpose//ArrayPlot

images=statgraph[g,#]&/@seenArr[[;;5]];
Export["C:\\research\\ideas\\animate2.swf",images]

ListLinePlot[seenArr[[50;;200,1]]]

start=100;
end=100+14*5;
images=statgraph[g,#]&/@seenArr[[start;;end]];
Export["C:\\research\\ideas\\animate2.swf",images]

Export["E:\\research\\ideas\\sum-product\\loopy-blogpost\\loopy-blogpost.html",SelectedNotebook[]]