Messages by hand

This page runs belief propagation on a small model by hand. You compute every message with the rules RxInfer runs, multiply the messages at a variable to get its posterior, and check the posterior against the closed form. Then infer runs the same model and returns the same posteriors.

A message is a function of one variable, sent along an edge of a factor graph. A rule computes one message of one factor node from the messages on the node's other edges. Message passing introduces the algorithm; this page computes it step by step.

using RxInfer, BayesBase

A message is an integral

A factor node $f(x, y_1, \dots, y_n)$ sends a message towards $x$ that integrates the node function against the messages on its other edges:

\[\mu_{f \to x}(x) = \int f(x, y_1, \dots, y_n) \prod_{i=1}^{n} \mu_{y_i \to f}(y_i) \, \mathrm{d}y_1 \cdots \mathrm{d}y_n .\]

This is the sum-product rule of belief propagation. The posterior of a variable is the normalized product of the messages that arrive at it. In a chain, every variable sits between two nodes, so its posterior is the product of a forward and a backward message.

The addition node

The node + relates three variables, $\mathrm{out} = \mathrm{in}_1 + \mathrm{in}_2$. Its node function is a Dirac delta, $f(\mathrm{out}, \mathrm{in}_1, \mathrm{in}_2) = \delta(\mathrm{out} - \mathrm{in}_1 - \mathrm{in}_2)$, and the function + itself is the node:

MessagePassingRulesBase.nodespec(+)

The interfaces are the node's edges: out, and in, a group with one member per term. A model creates this node when it writes s := z + u: z is the member (:in, 1) and u the member (:in, 2).

The message towards out integrates the delta out, which leaves a convolution of the two terms' messages:

\[\mu(\mathrm{out}) = \int \mu_{\mathrm{in}_1}(\mathrm{in}_1)\, \mu_{\mathrm{in}_2}(\mathrm{out} - \mathrm{in}_1) \, \mathrm{d}\mathrm{in}_1 .\]

For two normal messages $\mathcal{N}(m_1, v_1)$ and $\mathcal{N}(m_2, v_2)$, the convolution is $\mathcal{N}(m_1 + m_2, v_1 + v_2)$. @call_message_update_rule finds the rule for a node, a target and the inputs you give, and runs it. The messages on a group are a tuple, one per member:

@call_message_update_rule(
    node = +, target = :out,
    m = (in = (NormalMeanVariance(1.0, 1.0), NormalMeanVariance(2.0, 1.0)),),
)

args->sum_of_messages(args.m[:in])

The result draws the node. The input messages are solid arrows into the node, and the target is the double arrow out of it. The message is $\mathcal{N}(3, 2)$: the means add and the variances add. The value is a RuleResult, and getresult returns the message itself.

The message towards a term runs the sum backwards. Towards the second, (:in, 2), the delta sets $\mathrm{in}_2 = \mathrm{out} - \mathrm{in}_1$, and the message is $\mathcal{N}(m_{\mathrm{out}} - m_1, v_{\mathrm{out}} + v_1)$. The rule reads the other terms only, so the tuple holds nothing at the target's own position, as a graph delivers it:

@call_message_update_rule(
    node = +, target = (:in, 2),
    m = (out = NormalMeanVariance(5.0, 1.0), in = (NormalMeanVariance(1.0, 1.0), nothing)),
)

args->difference_message(args.m[:out], sum_of_messages(args.m[:in]))

The multiplication node with a known gain

The node * relates $\mathrm{out} = A \cdot \mathrm{in}$. When the gain $A = a$ is known, its message is a point mass, PointMass(a), a distribution with all its mass at $a$. The message towards out scales the input: a normal $\mathcal{N}(m, v)$ becomes $\mathcal{N}(a m, a^2 v)$.

@call_message_update_rule(
    node = *, target = :out,
    m = (A = PointMass(2.0), in = NormalMeanVariance(1.0, 1.0)),
)

args->scaled(mean(args.m[:A]), args.m[:in])

Towards in, the delta evaluates the message on out at $a \cdot \mathrm{in}$:

\[\mu(\mathrm{in}) = \int \delta(\mathrm{out} - a \cdot \mathrm{in})\, \mu_{\mathrm{out}}(\mathrm{out}) \, \mathrm{d}\mathrm{out} = \mathcal{N}(a \cdot \mathrm{in} \mid m, v) .\]

As a function of $\mathrm{in}$, this is a normal with precision $w = a^2 / v$ and weighted mean $\xi = a m / v$, the precision times the mean. It integrates to $1 / |a|$, not to one.

@call_message_update_rule(
    node = *, target = :in,
    m = (out = NormalMeanVariance(1.0, 1.0), A = PointMass(2.0)),
)

(ctx, args)->unscaled(ctx, args.m[:out], mean(args.m[:A]))

The message is $\mathcal{N}(0.5, 0.25)$ in weighted-mean form. Its log scale, $-\log 2$, is the logarithm of the constant $1/|a|$ that the message carries beyond a normalized distribution. The last section adds such constants up to the model's evidence. The card also lists matrix_correction, a service that the rule reads to keep a precision matrix invertible; its default serves here.

A model: one step of a tracking filter

A hidden position $x$ moves to $s = a x + u$, where $u$ is a noisy control input, and a sensor observes $s$ with noise:

\[\begin{aligned} x &\sim \mathcal{N}(0, 10), & u &\sim \mathcal{N}(1, 1), \\ s &= a x + u, & y &\sim \mathcal{N}(s, 0.5), \end{aligned}\]

with the gain $a = 2$ and the observation $y = 3$. Name the product $z = a x$. The factor graph is a chain of five nodes:

Belief propagation sends messages forwards from the priors to $s$, and backwards from the observation to $x$.

Forward messages

A prior node with known parameters sends its own distribution towards out:

a, y = 2.0, 3.0

prior_x = @call_message_update_rule(
    node = NormalMeanVariance, target = :out,
    m = (μ = PointMass(0.0), v = PointMass(10.0)),
)

args->NormalMeanVariance(mean(args.m[:μ]), mean(args.m[:v]))

The message on $x$ then passes through *, and meets the message on $u$ at +:

m_x = getresult(prior_x)
m_z = getresult(@call_message_update_rule(node = *, target = :out, m = (A = PointMass(a), in = m_x)))
m_u = getresult(@call_message_update_rule(node = NormalMeanVariance, target = :out, m = (μ = PointMass(1.0), v = PointMass(1.0))))
m_s = getresult(@call_message_update_rule(node = +, target = :out, m = (in = (m_z, m_u),)))

NormalMeanVariance{Float64}(μ=1.0, v=41.0)

The forward message on $s$ is the prediction of $s$ before the observation: $\mathcal{N}(2 \cdot 0 + 1,\; 4 \cdot 10 + 1)$.

Backward messages

The observation enters as a point mass on the likelihood's out. The message towards μ, the mean, is the likelihood of the observation as a function of $s$:

likelihood_s = @call_message_update_rule(
    node = NormalMeanVariance, target = :μ,
    m = (out = PointMass(y), v = PointMass(0.5)),
)

args->NormalMeanVariance(mean(args.m[:out]), mean(args.m[:v]))

The backward message then crosses + towards its first term, (:in, 1), which is $z$, and * towards in, which is $x$:

b_s = getresult(likelihood_s)
b_z = getresult(@call_message_update_rule(node = +, target = (:in, 1), m = (out = b_s, in = (nothing, m_u))))
b_x = getresult(@call_message_update_rule(node = *, target = :in, m = (out = b_z, A = PointMass(a))))

NormalWeightedMeanPrecision{Float64}(xi=2.6666666666666665, w=2.6666666666666665)

Posteriors

The posterior of $x$ is the product of its forward and backward messages. BayesBase's prod(ClosedProd(), left, right) multiplies two densities in closed form, and returns a normalized distribution:

q_x = prod(ClosedProd(), m_x, b_x)
q_s = prod(ClosedProd(), m_s, b_s)
mean_var(q_x)

(0.963855421686747, 0.3614457831325301)

The closed form agrees. Given $x$, the observation is $y \sim \mathcal{N}(a x + 1, 1.5)$, since the variances of $u$ and of the sensor noise add. The posterior precision of $x$ is the prior's plus $a^2 / 1.5$, and its mean is $a (y - 1) / 1.5$ divided by that precision:

w = 1 / 10 + a^2 / 1.5
(mean = a * (y - 1) / 1.5 / w, var = 1 / w)

(mean = 0.963855421686747, var = 0.3614457831325301)

The same model in RxInfer

The model below is the same chain. := defines a deterministic relationship, so z and s are the outputs of * and +:

@model function tracking_step(y, a)
    x ~ NormalMeanVariance(0.0, 10.0)
    u ~ NormalMeanVariance(1.0, 1.0)
    z := a * x
    s := z + u
    y ~ NormalMeanVariance(s, 0.5)
end

result = infer(model = tracking_step(a = a), data = (y = y,), free_energy = true)

Inference results:
  Posteriors       | available for (s, z, u, x)
  Free Energy:     | [2.82998]

infer runs the rules you ran, in the order their inputs arrive, and multiplies the messages at each variable:

(
    x = (by_hand = mean_var(q_x), rxinfer = mean_var(result.posteriors[:x])),
    s = (by_hand = mean_var(q_s), rxinfer = mean_var(result.posteriors[:s])),
)

(x = (by_hand = (0.963855421686747, 0.3614457831325301), rxinfer = (0.963855421686747, 0.3614457831325301)), s = (by_hand = (2.975903614457831, 0.4939759036144578), rxinfer = (2.975903614457831, 0.4939759036144578)))

RxInfer gives the rules of a stochastic node its observations and constants as point-mass marginals, q = (…), rather than as messages. For a point mass the two are the same distribution, and the rules return the same message.

The evidence

The model's evidence is the probability of the observation, $p(y) = \mathcal{N}(y \mid 1, a^2 \cdot 10 + 1 + 0.5)$. On a graph without loops, belief propagation is exact, and the Bethe free energy that free_energy = true computes equals $-\log p(y)$:

(free_energy = result.free_energy[end], minus_log_evidence = -logpdf(NormalMeanVariance(1.0, a^2 * 10 + 1.5), y))

(free_energy = 2.8299780179073366, minus_log_evidence = 2.8299780179073366)

The log scales of the messages add up to the same number. With logscales = true, every message and posterior carries its log scale, and the log scale of a posterior is $\log p(y)$:

with_logscales = infer(model = tracking_step(a = a), data = (y = y,), logscales = true)
getlogscale(with_logscales.posteriors[:x])

-2.8299780179073366

Next steps

• Variational message passing by hand runs a model where belief propagation has no closed form, and minimizes the free energy with a factorized posterior.
• Understanding rules explains how RxInfer finds the rule for each message.
• Creating your own custom nodes writes a node and its rules.