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 10_000 — larger than the other models, see the notes on conditioning below).
• verbose::Bool: show log messages (defaults to false).
• ftol::Union{Real, Nothing}: tolerance for the fixedpoint method (defaults to nothing, i.e. 1e-8). On this model it bounds the constraint residual — the largest absolute mismatch over the four constrained sequences — rather than the parameter-space increment used by the binary models. (The CReM/DCReM/CRWCM layers already use ftol this way.) The increment is not a usable test here: a runaway constraint has no finite fixed point, so the iterate keeps moving while the residual it is chasing falls, and an orbit that cycles can dip under an increment threshold and report success at an arbitrary 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

• :fixedpoint is block coordinate ascent (DECM_reduced_coordinate_iter!), not the Picard recipe that this model used up to v0.7.0. That recipe's strength step provably overshoots its own root without bound and left the model's domain on the first iteration from any cold start: measured over 150 random weighted digraphs it converged on 0. The replacement solves each block exactly and is feasible by construction; on the 71 of those digraphs that are well posed it converges on 71/71, with a median residual of 1.5·10⁻⁹ against 3.1·10⁻⁸ for :BFGS, and at a fraction of the cost — a 150-vertex weighted digraph takes 0.19 s against 593 s for :BFGS, at a residual 170× smaller. On the 79 carrying a runaway constraint (see Conditioning below) it reaches ftol on 22; there the optimum is at an infinite parameter and :BFGS (78/79) is the method to use — solve_model! says so, naming the offending constraint, rather than failing silently. Such a solve runs the full maxiters (10_000 by default on this model) before giving up, so lower it if you want the answer sooner. The legacy map is still exported as DECM_reduced_iter!.
• :BFGS (the default) and :Newton are the recommended gradient methods; :LBFGS converges to the same optimum but needs an order of magnitude more iterations (see below).

Conditioning of the DECM likelihood

The DECM objective is ill-conditioned in two distinct ways, and both are properties of the model and the data rather than defects — every solve below still reproduces the constraints to ~1e-9.

1. A two-fold gauge freedom, always present. The likelihood is exactly invariant under (α_out, α_in) → (α_out + c, α_in - c) and (β_out, β_in) → (β_out + c, β_in - c), because the pair terms depend only on the sums α_out,ᵢ + α_in,ⱼ and β_out,ᵢ + β_in,ⱼ while the linear part shifts by -c·(Σ F·k_out - Σ F·k_in) resp. -c·(Σ F·s_out - Σ F·s_in), both identically zero. θᵣ is therefore only determined up to that orbit, and the Hessian is singular along it. :Newton factorises the Hessian, so it needs the flat directions removed; solve_model! does that internally for :Newton only, which changes no gauge-invariant quantity (Ĝ, Ŵ and every metric are unaffected).
2. Runaway parameters from constraints at the edge of their feasible range. A saturated degree (k = N-1, the node links to everyone) drives α → -∞, and a minimum strength (s = k, every incident link carrying weight exactly 1) drives β → +∞. Each such constraint adds a near-null Hessian direction; condition numbers reach 1e15-1e17 on affected networks versus ~1e3 on clean ones. The corresponding fitness is simply not identifiable — the constraint sits at the boundary of what any ensemble can realise. (The rhesus_macaques network shipped with the package has one such s = k node.)

Measured iterations to convergence across a spread of networks, from the default :strengths guess and from the deliberately far :uniform one:

The limited-memory approximation simply cannot represent a Hessian of this condition number, which is why :LBFGS is not recommended here — it reaches the same optimum, but may exceed maxiters from a far start. :BFGS and :Newton are insensitive to the initial guess. Hence the larger maxiters default, which is only a cap and so costs the faster methods nothing.

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!).

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

One sweep of block coordinate ascent on the reduced DECM log-likelihood. θ is updated in place and returned; the four parameter buffers are overwritten. This is the map behind solve_model!(m, method=:fixedpoint).

The directed twin of UECM_reduced_coordinate_iter!, and it replaces the Picard recipe (DECM_reduced_iter!) for the same reason: with g = x_out,ᵢx_in,ⱼ, t = y_out,ᵢy_in,ⱼ, the degree step divides by a factor strictly decreasing in its own x and therefore undershoots safely, while the strength step divides by Bᵢ = Σⱼ wⱼ gⱼy_in,ⱼ/D(t) — strictly increasing in its own y whenever g < 2 — and therefore overshoots without bound, leaving the domain t < 1 on the first iteration. Measured over 150 random weighted digraphs, the Picard/Anderson path converged on 0.

Each block is instead solved exactly: ⟨kᵢ⟩(xᵢ) rises monotonically from 0 to Σⱼ wⱼ and ⟨sᵢ⟩(yᵢ) from 0 to ∞ on (0, ȳᵢ), so the strength root always exists, is unique, and is feasible by construction. A two-dimensional Newton polish (_decm_polish) follows runaway ridges.

See also validation/uecm_decm_fixedpoint.md.

_decm_polish(a₀, b₀, bmax, i, x_opp, y_opp, F, live, d, s)

Safeguarded two-dimensional Newton step on one orientation of class i, in (log xᵢ, log yᵢ); returns the polished (a, b).

The directed twin of _uecm_polish!, and it matters more here: the rhesus macaques network has an out-class with k = s = 1, whose optimum is at β_out → +∞. With the one-dimensional solves alone the residual stalls at 2·10⁻⁴ with a per-sweep rate of 0.999823 — some 7·10⁴ sweeps to reach 10⁻⁹. Every trial point must stay in the domain and reduce this channel's residual; that test is local, so the caller accepts the pass as a whole only when it lowers the global residual (see DECM_reduced_coordinate_iter!).

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.