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, 
  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.

  • 1
    $\begingroup$ Can you add a line of code that runs your code on example data? $\endgroup$
    – user5601
    Commented Jul 9, 2018 at 18:40
  • $\begingroup$ 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[9] == 1 && RV[5] == 1, listOfRV \[Distributed] \[ScriptCapitalD]] $\endgroup$ Commented Jul 10, 2018 at 6:15

1 Answer 1


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 = Developer`ToPackedArray[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


This is about 300 times faster than your original code.

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.


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]]



  • $\begingroup$ 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. $\endgroup$ Commented Jul 11, 2018 at 11:44
  • 1
    $\begingroup$ 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". $\endgroup$ Commented Jul 11, 2018 at 11:47
  • 1
    $\begingroup$ Kiril, I wonder wether this works for you. $\endgroup$ Commented Jul 14, 2018 at 20:50
  • $\begingroup$ Hi Henrik. It is work for me,my mistake I forget update you $\endgroup$ Commented Jul 15, 2018 at 5:09
  • $\begingroup$ Okay, fine!. Glad to be of help. $\endgroup$ Commented Jul 15, 2018 at 7:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.