Projection

CVIProjection handles the functions and message families that the Gaussian methods cannot. It draws samples of the inputs, pushes the samples through f, and projects the result onto an exponential family with ExponentialFamilyProjection. Projection here means finding the member of a chosen family, such as the normals or the gamma distributions, that is closest to the result.

Loading the rules

The rules of CVIProjection are in a package extension. Until ExponentialFamilyProjection is loaded, DeltaApproximation refuses the method, and its error says which package to load:

using DeltaMessagePassingRules, MessagePassingRulesBase, ExponentialFamily, BayesBase

try
    DeltaApproximation(method = CVIProjection())
catch err
    showerror(stdout, err)
end

ArgumentError: `CVIProjection(…)` is not an approximation method of the Delta node. It takes `Unscented()` and `Linearization()` from MessagePassingRulesApproximations, and `CVIProjection()` once ExponentialFamilyProjection is loaded. `CVIProjection`'s rules are in an extension: load it with `using ExponentialFamilyProjection`.

Loading the package loads the rules:

using ExponentialFamilyProjection

algorithm = DeltaApproximation(method = CVIProjection())

The method reads different inputs from the Gaussian methods. Towards out it reads the message on out itself, the marginal q(out) and the joint over the inputs:

MessagePassingRulesBase.dependencies_spec(DeltaFn, algorithm)

Example

Take out = exp(in), a standard normal message on in, and a gamma message on out, a positive quantity. The joint over the single input is the input's message reweighted by the likelihood of out, projected back onto the normals:

using Random: Xoshiro

struct WithFunction{F}
    f::F
end

MessagePassingRulesBase.getnodefn(node::WithFunction, ::MessagePassingRulesBase.Target{:out}) = node.f

ctx = MessagePassingRulesBase.RuleContext(node = WithFunction(exp), rng = Xoshiro(42))

joint = @call_marginal_update_rule(
    node = DeltaFn, target = (:in,), algorithm = algorithm, ctx = ctx,
    m = (out = Gamma(2.0, 1.0), in = (NormalMeanVariance(0.0, 1.0),)),
)

(algo, ctx, args)->begin
        g = getnodefn(ctx.node, Target(:out))
        (method, m_out, m_ins) = (algo.method, args.m[:out], args.m[:in])
        if length(m_ins) == 1
            m_in = only(m_ins)
            logp = generic_logpdf(m_in, (z->logpdf(m_out, g(z))))
            FactorizedJoint((project_to(create_project_to_ins(method, m_in, 1), logp, m_in),))
        else
            (proposal, strategy) = (method.proposal_distribution, method.sampling_strategy)
            columns = Any[input_samples(ctx.rng, source, strategy) for source = sampling_sources(proposal.distribution, m_ins)]
            loglikelihood(z, i) = mean((sample->logpdf(m_out, g(replace_at(sample, i, z)...))), zip(columns...))
            projections = ntuple(length(m_ins)) do i
                    m_in = m_ins[i]
                    prj = create_project_to_ins(method, m_in, i)
                    matches = ExponentialFamilyProjection.get_projected_to_type(prj) === ExponentialFamily.exponential_family_typetag(m_in) && ExponentialFamilyProjection.get_projected_to_dims(prj) == size(m_in)
                    projection = if matches
                            project_to(prj, generic_logpdf(m_in, (z->loglikelihood(z, i))), m_in)
                        else
                            project_to(prj, generic_logpdf(m_in, (z->loglikelihood(z, i) + logpdf(m_in, z))))
                        end
                    columns[i] = input_samples(ctx.rng, projection, strategy)
                    projection
                end
            result = FactorizedJoint(projections)
            proposal.distribution = result
            result
        end
    end

The rules draw their samples from the rule context's generator, ctx.rng, which an engine owns and a caller sets here. The node object is the same stand-in as on the overview. The message towards out samples that joint, pushes the samples through exp, and projects them onto the family of q(out), here the gamma distributions:

q_in = getresult(joint)

@call_message_update_rule(
    node = DeltaFn, target = :out, algorithm = algorithm, ctx = ctx,
    m = (out = Gamma(2.0, 1.0),), q = (out = Gamma(3.0, 1.0),), clusters = ((:in,) => q_in,),
)

(algo, ctx, args)->begin
        f = getnodefn(ctx.node, Target(:out))
        method = algo.method
        samples = map((q_in->cvilinearize(rand(ctx.rng, sampling_optimized(q_in), method.outsamples))), components(args.q[(:in,)]))
        out_samples = map((x->f(x...)), zip(samples...))
        estimate = project_to(create_project_to(method, args.q[:out], out_samples), out_samples)
        DivisionOf(estimate, args.m[:out])
    end

The result is the projection divided by the message that arrived on out, left unevaluated: its product with that message, at the variable, is the projection.

For a node with several inputs, the rule for the joint projects each input in turn, against samples of the others. A CVISamplingStrategy decides how those samples are drawn:

DeltaApproximation(method = CVIProjection(outsamples = 200, sampling_strategy = MeanBased()))

The rule for the joint keeps its result as the next call's proposal, in a ProposalDistributionContainer, so the method carries state and that rule is impure.

Limitations

• The rules need ExponentialFamilyProjection loaded.
• Sampling makes the messages random: two runs agree only with the same generator.
• The message towards out reads q(out) and the joint over the inputs, which depend on the node's own messages; a model whose out is not observed needs initial marginals for them.
• One CVIProjection value given to several nodes shares one proposal among them.
• A known inverse is not used, and is ignored with a warning.

API

The method:

CVIProjection(; outsamples = 100, out_prjparams = nothing, in_prjparams = nothing,
                proposal_distribution = ProposalDistributionContainer{Any}(nothing),
                sampling_strategy = FullSampling(10))

using DeltaMessagePassingRules, ExponentialFamilyProjection

algorithm = DeltaApproximation(method = CVIProjection(sampling_strategy = MeanBased()))

How the other inputs are sampled while one input of a node with several is projected:

CVISamplingStrategy

FullSampling(samples::Int = 10)

MeanBased()

The proposal the inputs are sampled from:

ProposalDistributionContainer(distribution)
ProposalDistributionContainer{Any}(nothing)