Markov Random Field

Question

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.

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

0.351

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.

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