DECM

Maximum entropy model for the Directed Enhanced Configuration Model (DECM).

The object holds the maximum likelihood parameters of the model (θ = [αout; αin; βout; βin]), the expected adjacency matrix (Ĝ), the expected weight matrix (Ŵ), and the variance of their elements (σ, σʷ).

The DECM constrains, per node, both the out- and in-degree sequence and the (integer) out- and in-strength sequence of a directed weighted network.

Note: this requires that the weights only assume (non-negative) integer values.

DECM(G::T; d_out::Vector, d_in::Vector, s_out::Vector, s_in::Vector, precision::Type{<:AbstractFloat}=Float64, kwargs...) where {T}

Constructor function for the DECM type.

By default and depending on the graph type T, the definition of in- and out-degree from $Graphs.jl$ and in- and out-strength from $SimpleWeightedGraphs$ is applied. If you want to use different definitions, you can pass the sequences as keyword arguments (d_out, d_in, s_out, s_in).

If you want to generate a model directly from the degree and strength sequences without an underlying graph, you can simply pass them all as keyword arguments. If you want to work from an adjacency/weight matrix, or an edge list, you can use the (weighted) graph constructors from the $JuliaGraphs$ ecosystem.

The strength sequences must be integer-valued (the DECM is defined for non-negative integer weights).

Examples

# generating a model from a directed weighted graph (the rhesus macaques grooming network)
julia> G = MaxEntropyGraphs.rhesus_macaques();

julia> model = DECM(G)
DECM{SimpleWeightedGraphs.SimpleWeightedDiGraph{Int64, Float64}, Float64} (16 vertices, 16 unique {out-degree, in-degree, out-strength, in-strength} quadruples, 1.00 compression ratio)

# generating a model directly from the degree and strength sequences
julia> model = DECM(d_out=[1, 1, 2, 1], d_in=[1, 2, 1, 1], s_out=[3, 3, 5, 2], s_in=[3, 5, 3, 2])
DECM{Nothing, Float64} (4 vertices, 4 unique {out-degree, in-degree, out-strength, in-strength} quadruples, 1.00 compression ratio)

# generating a model with a different precision
julia> model = DECM(d_out=[1, 1, 2, 1], d_in=[1, 2, 1, 1], s_out=[3, 3, 5, 2], s_in=[3, 5, 3, 2], precision=Float32);

julia> MaxEntropyGraphs.precision(model)
Float32

The degrees are taken from Graphs.outdegree/Graphs.indegree and the strengths from MaxEntropyGraphs.outstrength/MaxEntropyGraphs.instrength (the SimpleWeightedGraphs weighted degrees).

solve_model!(m::DECM; kwargs...)

Compute the likelihood maximising parameters of the DECM model m.

By default the parameters are computed using the BFGS method with the strength sequences as initial guess.

Arguments

• method::Symbol: solution method, :BFGS (default) or any of :BFGS, :LBFGS and :Newton, plus :fixedpoint.
• initial::Symbol: initial guess, :strengths (default), :strengths_minor, :random, or :uniform.
• maxiters::Int: maximum number of iterations (defaults to 1000).
• verbose::Bool: show log messages (defaults to false).
• ftol::Union{Real, Nothing}: tolerance for the fixedpoint method (defaults to nothing, i.e. 1e-8). It bounds the fixed-point increment $\|G(\theta) - \theta\|_\infty$ in parameter space; it is not the constraint residual. ❗ It applies to the :fixedpoint method only, and so is ignored on this model's default :BFGS path (passing it there warns). Use constraint_residual to measure how well the expected degrees and strengths actually match the observed ones.
• abstol, reltol: absolute/relative tolerances for the optimisation methods (default nothing).
• g_tol::Union{Number, Nothing}: gradient tolerance for the gradient-based methods (maps to Optim's g_abstol, default nothing). The gradient of this model is its constraint residual (up to the multiplicities), and it is the tolerance to reach for on the default path, but g_abstol is a stopping criterion rather than a guarantee: Optim can also stop on its function or parameter convergence checks and report success without the gradient ever reaching g_tol. Verify what was actually achieved with constraint_residual.
• AD_method::Symbol: autodiff method, any of :AutoZygote, :AutoReverseDiff, :AutoForwardDiff and :AutoFiniteDiff (defaults to :AutoZygote).
• analytical_gradient::Bool: use the analytical gradient instead of autodiff (defaults to false).

Notes

• the fixed-point method is very unstable for this model and should not be used. From an acceptable solution it can be used to fine-tune an existing one.
• the L-BFGS method is known to be unstable for this model and should not be used.

initial_guess(m::DECM; method::Symbol=:strengths)

Compute an initial guess θ₀ = [α_out; α_in; β_out; β_in] for the maximum likelihood parameters of the DECM model m.

The methods available are:

• :strengths (default): degrees normalised by the number of edges, strengths normalised by the total (reduced) strength.
• :strengths_minor: 1/(dᵣ+1) and 1/(sᵣ+1) per block.
• :random: random values drawn from $U(0,1)$.
• :uniform: uniformly set to -log(0.001).

rand(m::DECM; precomputed=false, rng=Random.default_rng())

Generate a random directed weighted graph from the DECM model m.

Arguments

• precomputed::Bool: not implemented yet for the DECM (the parameters are always used to generate the graph on the fly).
• rng::AbstractRNG: random number generator to use (defaults to Random.default_rng()).

Examples

julia> model = DECM(MaxEntropyGraphs.rhesus_macaques()); # generate a DECM model

julia> solve_model!(model, method=:BFGS); # compute the maximum likelihood parameters

julia> sample = rand(model); # sample a random directed weighted graph

julia> typeof(sample)
SimpleWeightedGraphs.SimpleWeightedDiGraph{Int64, Int64}

rand(m::DECM, n::Int; precomputed=false, rng=Random.default_rng())

Generate n random directed weighted graphs from the DECM model m. If multithreading is available, the graphs are generated in parallel; per-sample seeds are drawn from rng so the result is reproducible and independent of the thread schedule.

AIC(m::DECM)

Compute the Akaike Information Criterion (AIC) for the DECM model m. The parameters of the model must be computed beforehand. The DECM has 4N parameters (an $\alpha^{out}$, $\alpha^{in}$, $\beta^{out}$ and $\beta^{in}$ per node).

See also AICc, L_DECM_reduced.

AICc(m::DECM)

Compute the corrected Akaike Information Criterion (AICc) for the DECM model m.

See also AIC, L_DECM_reduced.

BIC(m::DECM)

Compute the Bayesian Information Criterion (BIC) for the DECM model m.

See also AIC, L_DECM_reduced.

Return the reduced number of {out-degree, in-degree, out-strength, in-strength} quadruples in the DECM network

L_DECM_reduced(θ::AbstractVector, d_out::Vector, d_in::Vector, s_out::Vector, s_in::Vector, F::Vector, n::Int=length(d_out))

Compute the log-likelihood of the reduced DECM model using the exponential formulation in order to maintain convexity.

Arguments

• θ: the maximum likelihood parameters of the model (θ = [α_out; α_in; β_out; β_in])
• d_out: the reduced outdegree sequence
• d_in: the reduced indegree sequence
• s_out: the reduced outstrength sequence
• s_in: the reduced instrength sequence
• F: the frequency of each (outdegree, indegree, outstrength, instrength) quadruple
• n: the number of unique quadruples (defaults to length(d_out))

The interaction term runs over the ordered pairs of reduced classes with multiplier F[i]·(F[j] - (i==j)) (a directed network has both (i,j) and (j,i) node pairs). The function is numerically stabilised (expm1/log1p) so that the 1 - exp(-β_out,ᵢ-β_in,ⱼ) denominator does not suffer from catastrophic cancellation, and stays automatic-differentiation friendly.

Examples

julia> θ = collect(range(0.1, step=0.1, length=12));

julia> d_out = [1, 2, 1]; d_in = [2, 1, 1]; s_out = [2, 3, 1]; s_in = [3, 2, 1]; F = [1, 1, 1];

julia> L_DECM_reduced(θ, d_out, d_in, s_out, s_in, F);

L_DECM_reduced(m::DECM)

Return the log-likelihood of the DECM model m based on the computed maximum likelihood parameters.

See also L_DECM_reduced(::AbstractVector, ::Vector, ::Vector, ::Vector, ::Vector, ::Vector)

∇L_DECM_reduced!(∇L::AbstractVector, θ::AbstractVector, d_out::Vector, d_in::Vector, s_out::Vector, s_in::Vector, F::Vector, x_out::AbstractVector, x_in::AbstractVector, y_out::AbstractVector, y_in::AbstractVector, n=length(θ)÷4)

Compute the gradient of the log-likelihood of the reduced DECM model in a non-allocating manner. The pre-allocated buffers x_out/x_in (xᵢ = exp(-αᵢ)) and y_out/y_in (yᵢ = exp(-βᵢ)) and the gradient ∇L are updated in place.

Each row i accumulates both orientations in a single inner loop: i→j feeds the out-blocks (α_out, β_out) and j→i feeds the in-blocks (α_in, β_in). The inner loop is branch-free (the diagonal correction is folded into the multiplier F[j] - (i==j)) so that it vectorises (@simd).

See also ∇L_DECM_reduced_minus!.

∇L_DECM_reduced_minus!(args...)

Compute minus the gradient of the log-likelihood of the reduced DECM model (used for the minimisation carried out by Optimization.jl). Non-allocating: updates the pre-allocated buffers x_out, x_in, y_out, y_in and ∇L in place.

See also ∇L_DECM_reduced!.

DECM_reduced_iter!(θ, d_out, d_in, s_out, s_in, F, nz_out, nz_in, x_out, x_in, y_out, y_in, G, n=length(θ)÷4)

Compute the next fixed-point iteration for the reduced DECM model. The pre-allocated buffers x_out, x_in, y_out, y_in and G are updated in place. Only the live channels are iterated: the out-blocks over i ∈ nz_out (with the inner sum over j ∈ nz_in) and the in-blocks over i ∈ nz_in (inner sum over j ∈ nz_out) — a dead channel contributes an exact zero to every pair probability.

Note: the fixed-point recipe is unstable for the DECM (as it is for the UECM); :BFGS/:Newton are preferred (see solve_model!).

set_xᵣ!(m::DECM)

Set the values of the degree-channel parameters xᵣout to exp(-αout) and xᵣin to exp(-αin) for the DECM model m

set_yᵣ!(m::DECM)

Set the values of the strength-channel parameters yᵣout to exp(-βout) and yᵣin to exp(-βin) for the DECM model m

Ĝ(m::DECM)

Compute the expected adjacency matrix for the DECM model m. The matrix is asymmetric: entry (i,j) is the probability of the directed edge i→j.

Note: The expected weights can be computed separately with Ŵ.

set_Ĝ!(m::DECM)

Set the expected adjacency matrix for the DECM model m

Ŵ(m::DECM)

Compute the expected (unconditional) weighted adjacency matrix for the DECM model m, i.e. ⟨wᵢⱼ⟩ = pᵢⱼ / (1 - yᵢ_out·yⱼ_in), so that sum(Ŵ(m), dims=2) ≈ outstrength and sum(Ŵ(m), dims=1)' ≈ instrength.

set_Ŵ!(m::DECM)

Set the expected weighted adjacency matrix for the DECM model m.

σˣ(m::DECM{T,N}) where {T,N}

Compute the standard deviation for the elements of the (binary) adjacency matrix for the DECM model m, i.e. sqrt(pᵢⱼ(1 - pᵢⱼ)) (the adjacency entries are Bernoulli distributed).

Note: this is the standard deviation of the binary layer; the standard deviation of the weights is available via σʷ. Read as "sigma star".

set_σ!(m::DECM)

Set the standard deviation for the elements of the (binary) adjacency matrix for the DECM model m

σʷ(m::DECM)

Compute the standard deviation for the elements of the weighted adjacency matrix for the DECM model m. The weight wᵢⱼ of the directed edge i→j follows a Bernoulli–geometric mixture: with pᵢⱼ the connection probability and yᵢ_out·yⱼ_in the geometric parameter, $⟨w_{ij}⟩ = p_{ij}/(1 - y_iy_j)$, $⟨w_{ij}^2⟩ = p_{ij}(1 + y_iy_j)/(1 - y_iy_j)^2$ and

$Var(w_{ij}) = \frac{p_{ij}(1 + y_iy_j - p_{ij})}{(1 - y_iy_j)^2}$

As the network is directed, wᵢⱼ and wⱼᵢ are distinct, independent random variables (unlike the UECM).

set_σʷ!(m::DECM)

Set the standard deviation for the elements of the weighted adjacency matrix for the DECM model m.

precision(m::DECM)

Determine the compute precision of the DECM model m.

A(m::DECM, i::Int, j::Int)

Return the expected value of the adjacency matrix for the DECM model m at the (ordered) node pair (i,j).

❗ For performance reasons, the function does not check:

• if the node pair is valid.
• if the parameters of the model have been computed.

f_DECM(xixj::T, yiyj::T) where {T}

Helper for the DECM model computing the expected adjacency entry pᵢⱼ = (xᵢxⱼ·yᵢyⱼ)/(1 - yᵢyⱼ + xᵢxⱼ·yᵢyⱼ) for the ordered pair i→j from the products xᵢxⱼ = xᵢ_out·xⱼ_in and yᵢyⱼ = yᵢ_out·yⱼ_in of the maximum likelihood parameters.

σₓ(m::DECM, X::Function; layer::Symbol=:binary, gradient_method::Symbol=:ReverseDiff)

Compute the standard deviation of metric X for the DECM model m via error propagation (the delta method of Squartini & Garlaschelli (2011), Eq. B.16).

Arguments

• layer::Symbol:
  • :binary (default): propagate over the binary adjacency matrix — X is a function of the adjacency matrix, the gradient is evaluated at m.Ĝ and weighted by m.σ (requires set_Ĝ! and set_σ!).
  • :weighted: propagate over the weighted adjacency matrix — X is a function of the weight matrix, the gradient is evaluated at m.Ŵ and weighted by m.σʷ (requires set_Ŵ! and set_σʷ!).
• gradient_method::Symbol: :ForwardDiff, :ReverseDiff (default) or :Zygote.

As the network is directed, the entries (i,j) and (j,i) of either layer are independent random variables, so the delta method carries no within-dyad covariance cross-term (unlike the undirected UECM).

Metrics mixing the two layers (functions of both the adjacency and the weight matrix) are not supported by this per-layer propagation; use ensemble sampling instead.