# Markov Random Field

I look for an efficient way to implement a Pairwise Markov Random Field.

I have implemented a homogeneous (i.e., unary potential is a similar function in all nodes, and pairwise potential is a similar function too) binary Markov Random Field.

For an example probability calculation on the line graph with ten nodes.

  nodes = 10;
nodesDomain = {1, -1};
numOfObservation = 4;
g1 = GridGraph[{1, nodes}, VertexLabels -> "Name"];
possibleGraphAssignment = Tuples[nodesDomain, nodes];
weights =
ParallelMap[graphPotentialCalculation[#] &, possibleGraphAssignment];
\[ScriptCapitalD]= EmpiricalDistribution[weights -> possibleGraphAssignment]/. {b -> 0.3, a -> 0.7, q -> 0.1, p -> 0.9};
listOfRV = Parallelize[Array[RV, nodes]];


With this function I calculate the potential of each possible configuration:

graphPotentialCalculation[vertC_] :=
Module[{edgeList = EdgeList[g1], vertexList = VertexList[g1], nodeP,
edgeP},
nodeP = nodePotential[#] & /@ vertC;
edgeP = edgePotential[First[#], Last[#], vertC] & /@ edgeList;
Times @@ nodeP*Times @@ edgeP
]

edgePotential[v_, u_, graphCh_] := Module[{},
Piecewise[{{p, graphCh[[v]] == graphCh[[u]]}, {q,
graphCh[[v]] != graphCh[[u]]}}]]

nodePotential[v_] := Module[{},
Piecewise[{{a, v == 1}, {b, v == -1}}]]


I have a bottleneck in my code : calculation of the empirical distribution - to do it, I calculate the truth table, and this takes a lot of time

Any suggestion on how to speed-up code.

• Can you add a line of code that runs your code on example data? – user5601 Jul 9 '18 at 18:40
• The code in the question is run code. You can query anything about the probability in this specific graph (line graph with ten nodes and this spesific potential functions), for example the join probability of node one and nine to get value one: Probability[RV == 1 && RV == 1, listOfRV \[Distributed] \[ScriptCapitalD]] – Kiril Danilchenko Jul 10 '18 at 6:15

For maximum performance, you have to replace a, b, p, and q by numerical values at the very beginning. In the code below, I try to use integer arithmetic instead of boolean computations and to use logarithms in order to cast the many multiplications into less expensive additions. In order to keep the data transfer as low as possible, I use IntegerDigits to compute the i-th possible charge/spin/whatever distributions directly from i. I deploy the computational core into a CompiledFunction basically to exploit the low level paralellization it provides.

possGraphAssignment[i_] := Subtract[1, 2 IntegerDigits[i, 2, nodes]];

cf = Compile[{
{i, _Integer}, {nodes, _Integer}, {i1, _Integer, 1}, {i2, _Integer, 1},
{Loga, _Real}, {Logb, _Real}, {Logp, _Real}, {Logq, _Real}
},
Module[{nodeP, edgeP},
nodeP = # Logb + Subtract[1, #] Loga &[IntegerDigits[i, 2, nodes]];
edgeP = # Logq + Subtract[1, #] Logp &[Unitize[nodeP[[i1]] - nodeP[[i2]]]];
Exp[Total[nodeP] + Total[edgeP]]
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];


Let's try a simple example. First, we compute some static data.

nodes = 20;
g1 = GridGraph[{1, nodes}, VertexLabels -> "Name"];

a = 0.7;
b = 0.3;
p = 0.9;
q = 0.1;

Loga = Log[a];
Logb = Log[b];
Logp = Log[p];
Logq = Log[q];
edgeList = DeveloperToPackedArray[List @@@ EdgeList[g1]];
{i1, i2} = Transpose[edgeList];


Now the computations of the weights:

weights2 = cf[Range[0, 2^nodes - 1], nodes, i1, i2, Loga, Logb, Logp, Logq]; // RepeatedTiming // First


0.351

The empirical distribution (over the integers belonging to the possible distributions) can be obtained this way (within a second):

\[ScriptCapitalD] = EmpiricalDistribution[weights -> Range[0, 2^nodes - 1]];


In the long run (for significantly greater values of nodes), you still won't have fun with this due to inherent combinatorical explosion.

# Edit

An even faster variant is this:

cf2 = Compile[{
{i, _Integer}, {nodes, _Integer}, {i1, _Integer,
1}, {i2, _Integer, 1}, {a, _Real}, {p, _Real}
},
Module[{nodeP, m, n},
nodeP = IntegerDigits[i, 2, nodes];
m = Total[nodeP];
n = Total[Unitize[nodeP[[i1]] - nodeP[[i2]]]];
(1. - a)^m a^(nodes - m) (1. - p)^n p^(Length[i1] - n)
],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True,
RuntimeOptions -> "Speed"
];

weights3 = cf2[Range[0, 2^nodes - 1], nodes, i1, i2, a, p]; // RepeatedTiming // First
Max[Abs[weights2 - weights3]]


0.219

1.21973*10^-19

• Thanks a lot. I have some issue with Compile when I run it with the CompilationTarget -> "C" I get the error, without this option, it is run ok. CreateLibrary::nocomp: A C compiler cannot be found on your system. Please consult the documentation to learn how to set up suitable compilers. Compile::nogen: A library could not be generated from the compiled function. – Kiril Danilchenko Jul 11 '18 at 11:44
• Usually, I would advice you to install a C compiler on your computer, but this is one of the rare occasions where this would improve the speed only marginally (10 % or so). You may compile to Wolfram Virtual Machine (WVM) by using CompilationTarget -> "WVM" instead of CompilationTarget -> "C"`. – Henrik Schumacher Jul 11 '18 at 11:47
• Kiril, I wonder wether this works for you. – Henrik Schumacher Jul 14 '18 at 20:50
• Hi Henrik. It is work for me,my mistake I forget update you – Kiril Danilchenko Jul 15 '18 at 5:09
• Okay, fine!. Glad to be of help. – Henrik Schumacher Jul 15 '18 at 7:17