GCVMessagePassingRules

The Gaussian controlled variance node, GCV, and its rules for reactive message passing: a normal whose log-variance is linear in other variables of the model. Use it for a variance that is itself inferred and changes over time, as in hierarchical Gaussian filters, where the state of one layer sets the volatility of the layer below.

GCVMessagePassingRules

julia> using MessagePassingRulesBase, ExponentialFamily, BayesBase

julia> q_z, q_κ, q_ω = PointMass(0.0), PointMass(1.0), PointMass(0.0);  # a noise variance exp(0) = 1

julia> result = @call_message_update_rule(
           node = GCV, target = :y, m = (x = NormalMeanVariance(1.0, 2.0),), q = (z = q_z, κ = q_κ, ω = q_ω),
       );

julia> mean_var(getresult(result)) .≈ (1.0, 3.0)
(true, true)

Overview

A GCV node relates an output y to its mean x with a variance exp(κz + ω). Here z is a latent variable, often a random walk of its own. κ scales the effect of z on the log-variance, and ω is an offset.

The rules are variational: each rule sends the exponentiated expected log-density of the node, the expectation taken under the marginals of the variables outside the target's cluster. Every rule reads the marginals of z, κ and ω. The messages towards z, κ and ω are ExponentialLinearQuadratic densities, whose moments a Gauss–Hermite cubature computes.

The package also adds rules to the NormalMeanVariance and NormalMeanPrecision nodes of StandardMessagePassingRules. These rules accept an ExponentialLinearQuadratic message on the normal's out. With them, z or ω can itself be the output of a normal node, such as a step of a random walk.

Definition

\[p(y \mid x, z, κ, ω) = \mathcal{N}\bigl(y \mid x, \exp(κ z + ω)\bigr) = \frac{1}{\sqrt{2π \exp(κ z + ω)}} \exp\Bigl(-\frac{(y - x)^2}{2 \exp(κ z + ω)}\Bigr).\]

The rules need the expected noise precision ⟨exp(-(κz + ω))⟩ = ⟨exp(-ω)⟩ ⟨exp(-κz)⟩. For a normal q(ω) the first factor is the mean of a lognormal, exp(-⟨ω⟩ + Var(ω)/2), which is exact. The second treats κz as normal, with Var(κz) = ⟨κ⟩² Var(z) + ⟨z⟩² Var(κ) + Var(κ) Var(z), which is exact only when κ or z is a point mass.

Interfaces

using GCVMessagePassingRules, MessagePassingRulesBase, ExponentialFamily, BayesBase
MessagePassingRulesBase.nodespec(GCV)

No interface of GCV has an alias. Under the structured factorisation y and x share one cluster, and so a joint marginal q(y, x), a bivariate normal.

Algorithm

GCVApproximation(; method = GaussHermiteCubature(20)) is the node's algorithm. It carries the GaussHermiteCubature that computes the moments of the messages towards z, κ and ω. The node declares its default instance, so a model does not need to name it, and names it only to change the number of points. The rules are declared for the default instance's type, so another approximation method finds no rule.

The algorithm declares no dependencies of its own. The rules take what the default scheme delivers: messages from the target's cluster and marginals of the other clusters.

Supported rules

MessagePassingRulesBase.rule_coverage(GCV)

Each target of GCV has two rules. One is for the structured factorisation q(y, x) q(z) q(κ) q(ω), and one is for the mean field, which puts every variable in a cluster of its own. Under the structured factorisation:

• the rules towards y and x take the message on the other one and the marginals of z, κ and ω;
• the rules towards z, κ and ω take the joint q(y, x);
• the row q(y, x) is the marginal rule of the joint.

There are no belief propagation rules. The average energy is defined for both factorisations.

The rules that the package adds to the normal nodes of StandardMessagePassingRules take an ExponentialLinearQuadratic message on out, which they reduce to a normal with its moments:

Example

The message towards y from a normal message on x, with point-mass marginals that make the noise variance exp(1 ⋅ 0 + 0) = 1:

@call_message_update_rule(
    node = GCV, target = :y,
    m = (x = NormalMeanVariance(1.0, 2.0),),
    q = (z = PointMass(0.0), κ = PointMass(1.0), ω = PointMass(0.0)),
)

args->through_noise(args.m[:x], args.q[:z], args.q[:κ], args.q[:ω])

The card draws the message on x as a solid input and the marginals of z, κ and ω as dashed ones. The result adds the noise variance to the variance of x, 2 + 1 = 3:

julia> using MessagePassingRulesBase, GCVMessagePassingRules, ExponentialFamily, BayesBase

julia> result = @call_message_update_rule(
           node = GCV, target = :y,
           m = (x = NormalMeanVariance(1.0, 2.0),),
           q = (z = PointMass(0.0), κ = PointMass(1.0), ω = PointMass(0.0)),
       );

julia> mean_var(getresult(result)) .≈ (1.0, 3.0)
(true, true)

The message towards ω under the mean field is an ExponentialLinearQuadratic, whose moments the cubature of the default algorithm computes. Its coefficient b is ⟨(y - x)²⟩ = (1 - 0)² + 0.5 + 0.5 = 2:

julia> using MessagePassingRulesBase, GCVMessagePassingRules, ExponentialFamily, BayesBase

julia> result = @call_message_update_rule(
           node = GCV, target = :ω,
           q = (y = NormalMeanVariance(1.0, 0.5), x = NormalMeanVariance(0.0, 0.5), z = PointMass(0.0), κ = PointMass(1.0)),
       );

julia> message = getresult(result);

julia> message isa ExponentialLinearQuadratic
true

julia> params(message) == (1.0, 2.0, -1.0, 0.0)
true

A step of a hierarchical Gaussian filter

In a two-layer hierarchical Gaussian filter, the upper layer z is a random walk, zₜ ~ N(zₜ₋₁, 1), and it sets the volatility of the lower random walk, xₜ ~ GCV(xₜ₋₁, zₜ, κ, ω). The update of the upper layer passes through two nodes. First the GCV node sends its message towards zₜ, from the marginals of the lower layer:

towards_z = @call_message_update_rule(
    node = GCV, target = :z,
    q = (y = NormalMeanVariance(1.0, 0.5), x = NormalMeanVariance(0.0, 0.5), κ = PointMass(1.0), ω = PointMass(0.0)),
)
m_z = getresult(towards_z)
mean_var(m_z)

(1.814142762171816, 3.673218039648464)

The message is an ExponentialLinearQuadratic, and its moments come from the cubature. Then the random walk's NormalMeanVariance node carries it back towards zₜ₋₁. This is one of the rules the package adds to that node:

@call_message_update_rule(node = NormalMeanVariance, target = :μ, m = (out = m_z, v = PointMass(1.0)))

args->begin
        (m, v) = mean_var(args.m[:out])
        NormalMeanVariance(m, v + mean(args.m[:v]))
    end

The rule reduces the message on out to a normal with its moments and adds the random walk's variance, 1.

Limitations

• There are no belief-propagation rules: every rule needs the marginals of z, κ and ω, so a model must initialise them.
• The node is univariate; the messages on y and x must be univariate normals or ExponentialLinearQuadratics.
• ⟨exp(-κz)⟩ treats κz as normal, which is an approximation unless κ or z is a point mass.
• The approximation method must be a Gauss–Hermite cubature.
• The average energy needs a normal q(z), and a multivariate normal q(y, x) or normal q(y) and q(x).
• The rules declare no log scale.

API

The node and its algorithm:

GCV

GCVApproximation(; method = GaussHermiteCubature(20))
GCVApproximation(method)

julia> using MessagePassingRulesApproximations: GaussHermiteCubature

julia> GCVApproximation(method = GaussHermiteCubature(32)) isa GCVApproximation
true

The density of the messages towards z, κ and ω:

ExponentialLinearQuadratic(approximation, a, b, c, d)

julia> using MessagePassingRulesApproximations: GaussHermiteCubature

julia> using ExponentialFamily, BayesBase

julia> d = ExponentialLinearQuadratic(GaussHermiteCubature(20), 1.0, 2.0, -1.0, 0.0);

julia> logpdf(d, 0.0) ≈ -1.0
true

julia> prod(GenericProd(), NormalMeanVariance(0.0, 1.0), d) isa NormalMeanVariance
true