Arithmetic

The arithmetic nodes +, -, * and dot are deterministic nodes: their output is a function of their inputs, such as out = in[1] + in[2]. Base's and LinearAlgebra's functions are the nodes themselves, so an RxInfer model writes y ~ x + z or y ~ A * x, and rule_coverage(+) takes the function. The package re-exports LinearAlgebra's dot for models to use.

A deterministic node's clusters are always out and the joint over its inputs, whatever the factorisation. Every rule here is therefore belief propagation over messages. The message towards out is the pushforward of the input messages, the distribution of f(in...). The message towards an input inverts the function against the other messages. For normal messages, with one factor of * or dot known, these messages are exact and normal. Between two uncertain factors of * the result has no closed form, and the rules return log-densities or samples instead. None of the nodes has an average energy of its own.

+, - and dot run under DefaultAlgorithm, and * runs under MultiplicationSampling, its default. No model names an algorithm for them.

Example

The sum of a normal and a known value shifts the normal. A known factor of * scales the mean, and it scales the variance by its square:

julia> using StandardMessagePassingRules, MessagePassingRulesBase, ExponentialFamily, BayesBase

julia> sum_message = getresult(@call_message_update_rule(
           node = +, target = :out,
           m = (in = (NormalMeanVariance(1.0, 1.0), PointMass(2.0)),),
       ));

julia> mean(sum_message) ≈ 3.0 && var(sum_message) ≈ 1.0
true

julia> scaled = getresult(@call_message_update_rule(
           node = *, target = :out,
           m = (A = PointMass(2.0), in = NormalMeanVariance(1.0, 1.0)),
       ));

julia> mean(scaled) ≈ 2.0 && var(scaled) ≈ 4.0
true

Addition, +

\[\mathrm{out} = \mathrm{in}_1 + \mathrm{in}_2 + \dots + \mathrm{in}_n, \quad n \ge 2\]

One node sums any number of terms, so a + b + c, which Julia parses as one call of +, is one node, with no variable for an intermediate sum. That matters in a model with loops, where an intermediate sum needs initial messages of its own and changes the schedule.

MessagePassingRulesBase.nodespec(+)

MessagePassingRulesBase.rule_coverage(+)

The message towards out adds the means and the variances of the terms:

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

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

Towards a term, the rule subtracts the others from out. With out observed at 5, the message towards the second of two terms is a normal around 5 - 1. The group's tuple holds nothing at the target's own position:

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

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

The messages are normal whenever any term is normal, and a PointMass when every term is known. A normal shifted by a known value keeps its mean–precision or weighted-mean form. Any other result is in mean–variance or mean–covariance form. Towards out, Distributions' convolve sums any other distributions, and it throws a MethodError for a pair it does not know. Every rule has log scale zero.

The joint marginal of the terms, which the free energy takes the entropy of, depends on out's message:

• with a normal out, the normal terms are one joint normal, their messages times out's message of their sum, and each known term is a block of its own;
• with out known, the terms lie on the plane where they sum to it. Every normal term but the last is free, and the last is the rest of the sum: an StandardMessagePassingRules.InputsGivenSum.

@call_marginal_update_rule(
    node = +, target = (:in,),
    m = (out = PointMass(3.0), in = (NormalMeanVariance(2.0, 2.0), PointMass(2.0), NormalMeanVariance(0.0, 1.0))),
)

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

StandardMessagePassingRules.InputsGivenSum

Subtraction, -

\[\mathrm{out} = \mathrm{in}_1 - \mathrm{in}_2\]

MessagePassingRulesBase.rule_coverage(-)

Its rules are +'s, rearranged. Towards out, the message is the difference of the inputs. Towards in1, it is the sum of out and in2, and it takes any two distributions through convolve, as + does towards out. Towards in2, it is the difference of in1 and out.

Multiplication, *

\[\mathrm{out} = A \cdot \mathrm{in}\]

MessagePassingRulesBase.nodespec(*)

MessagePassingRulesBase.rule_coverage(*)

With a known matrix A, the message towards out maps the normal message on in through A: the mean becomes A μ and the covariance A Σ Aᵀ:

@call_message_update_rule(
    node = *, target = :out,
    m = (A = PointMass([1.0 1.0; 0.0 1.0]), in = MvNormalMeanCovariance([1.0, 2.0], diageye(2))),
)

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

Towards in, the message is a normal in weighted-mean form, with precision Aᵀ Σ⁻¹ A and weighted mean Aᵀ Σ⁻¹ m from the message N(m, Σ) on out. The card lists the service the rule reads, matrix_correction:

@call_message_update_rule(
    node = *, target = :in,
    m = (out = MvNormalMeanCovariance([3.0, 2.0], diageye(2)), A = PointMass([1.0 1.0; 0.0 1.0])),
)

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

With one factor known, the rules are in closed form. Towards out, the message is a scaled normal or Gamma. Towards the unknown factor, it is a normal in weighted-mean form or a Gamma, with the log scale -d log |c| of a known scalar c. A known factor may be a scalar, a vector or a matrix. With a vector, the other factor is a scalar. With a matrix, the rules compute only A * in, never in * A, so they refuse a matrix in towards A. The precision matrices that the rules towards an input build may be singular. The rules correct them with the context's matrix_correction, ReplaceZeroDiagonalEntries(tiny) by default.

Between two univariate normals, the messages are log-densities, a ContinuousUnivariateLogPdf. Towards an input, the rule computes the exact integral. Towards out, it sums a Bessel series truncated at ten terms. Between two other univariate distributions, the rules draw samples from the context's rng, as many as the algorithm says. These messages are unnormalised and declare no log scale.

MultiplicationSampling(; samples::Int = 3000)

julia> MultiplicationSampling().samples
3000

julia> MultiplicationSampling(samples = 500).samples
500

Dot product, dot

\[\mathrm{out} = \mathrm{in}_1^\top \mathrm{in}_2\]

MessagePassingRulesBase.rule_coverage(dot)

Limitations. One input must be known. Two normal inputs have no closed form. The rules for that case exist only to throw an error that points to the SoftDot node, whose package handles it. The precision a w aᵀ towards an input has rank one, and the rule corrects it with the context's matrix_correction, ReplaceZeroDiagonalEntries(tiny) by default.