DCReM

Maximum entropy model for the directed Conditional Reconstruction Method (DCReM).

The DCReM is the directed counterpart of the CReM: a two-step null model for weighted, directed networks with continuous, positive weights (Parisi, Squartini & Garlaschelli (2020); the model is called CReM_A in the literature and DBCM+CReMa in the NuMeTriS package). The binary structure (topology) is supplied by a prior binary model — here an internally solved DBCM on the out- and in-degree sequences, giving the marginal edge probability fᵢⱼ = xᵢyⱼ/(1 + xᵢyⱼ) with xᵢ = e^{-αᵢ}, yᵢ = e^{-βᵢ}. Conditional on an edge i→j existing, its weight follows an exponential distribution with rate θᵒᵢ + θⁱⱼ (mean 1/(θᵒᵢ + θⁱⱼ)); the parameters θ = [θᵒ; θⁱ] (two per node) constrain the out- and in-strength sequences.

The object holds the maximum likelihood parameters of the weighted layer (θ), the reduced binary (conditional) DBCM parameters (αᵣ) and their exponentiated forms (xᵣ, yᵣ).

DCReM(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 DCReM type.

By default and depending on the graph type T, the definition of out/in-degree from $Graphs.jl$ and out/in-strength from $SimpleWeightedGraphs$ is applied. If you want to use different definitions, you can pass the vectors as keyword arguments. If you want to generate a model directly from degree and strength sequences without an underlying graph you can simply pass them as keyword arguments.

The DCReM allows continuous, positive weights (the strength sequences need not be integer-valued).

Examples

# generating a model from a weighted directed graph
julia> model = DCReM(MaxEntropyGraphs.rhesus_macaques())
DCReM{SimpleWeightedGraphs.SimpleWeightedDiGraph{Int64, Float64}, Float64} (16 vertices, 15 unique degree pairs, 0.94 compression ratio)

# generating a model directly from degree and strength sequences
julia> model = DCReM(d_out=[1, 1, 2, 1], d_in=[2, 1, 1, 1], s_out=[3.0, 5.0, 4.0, 2.0], s_in=[4.0, 2.0, 5.0, 3.0])
DCReM{Nothing, Float64} (4 vertices, 3 unique degree pairs, 0.75 compression ratio)

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

Compute the likelihood maximising parameters of the DCReM model m. This is a two-step process: first the binary (conditional) DBCM layer is solved on the out/in-degree sequences, then the weighted DCReM layer is solved on the out/in-strength sequences.

By default the weighted layer is computed with the fixed-point method using the strength sequences as initial guess (the CReM-family fixed-point recipe is stable).

Arguments (weighted DCReM layer)

• method::Symbol: solution method, :fixedpoint (default) or any of :BFGS, :LBFGS and :Newton.
• initial::Symbol: initial guess, :strengths (default), :strengths_minor or :random.
• 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).
• store_adjacency::Bool: cache the (binary) expected adjacency matrix m.Ĝ and use it in the weighted solve (defaults to false).

Arguments (binary conditional DBCM layer)

• method_conditional::Symbol: solution method for the binary layer (defaults to :fixedpoint).
• initial_conditional::Symbol: initial guess for the binary layer (defaults to :degrees).
• AD_method_conditional::Symbol: autodiff method for the binary layer (defaults to :AutoZygote).
• analytical_gradient_conditional::Bool: analytical gradient for the binary layer (defaults to false).

Common settings

• 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), applied to the weighted layer and to the binary layer when either is solved with :fixedpoint (passing it when neither layer uses it warns). On the weighted layer it is a relative strength tolerance (see below); on the binary layer it bounds the fixed-point increment in parameter space.
• 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 the weighted layer is its constraint residual in strength units, 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.

Examples

# default use
julia> model = DCReM(MaxEntropyGraphs.rhesus_macaques());

julia> solve_model!(model);

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

Compute an initial guess θ₀ = [θᵒ₀; θⁱ₀] for the maximum likelihood parameters of the weighted (DCReM) layer of the model m.

The DCReM parameters θ are the direct rate parameters (they appear as log(θᵒᵢ + θⁱⱼ) in the log-likelihood), so every method returns strictly positive, feasible values — a negative θᵒᵢ + θⁱⱼ would put the objective outside its domain. The forms mirror NEMtropy's crema initial guesses:

• :strengths (default): θᵒᵢ = 𝟙[s^{out}_i > 0] / (Σⱼ s^{out}_j + Σⱼ s^{in}_j) and the analog for θⁱᵢ.
• :strengths_minor: θᵒᵢ = 𝟙[s^{out}_i > 0] / (s^{out}_i + 1) and the analog for θⁱᵢ.
• :random: random values drawn from $U(0,1)$.

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

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

The binary structure is drawn from the (binary) DBCM layer (fᵢⱼ = xᵢyⱼ/(1 + xᵢyⱼ)); for each realised edge i→j a continuous weight is drawn from an exponential distribution with rate θᵒᵢ + θⁱⱼ.

Arguments

• precomputed::Bool: not implemented yet for the DCReM (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 = DCReM(MaxEntropyGraphs.rhesus_macaques()); # generate a DCReM model

julia> solve_model!(model); # compute the maximum likelihood parameters

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

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

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

Generate n random weighted directed graphs from the DCReM 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::DCReM)

Compute the Akaike Information Criterion (AIC) for the DCReM model m. The parameters of the model must be computed beforehand. The (conditional) DCReM has 2N parameters (one out- and one in-rate $θ$ per node; the binary layer's fᵢⱼ are fixed inputs of the conditional likelihood).

See also AICc, L_DCReM.

AICc(m::DCReM)

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

See also AIC, L_DCReM.

BIC(m::DCReM)

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

See also AIC, L_DCReM.

Return the number of vertices in the DCReM network

L_DCReM(θ::AbstractVector, s_out::AbstractVector, s_in::AbstractVector, x::AbstractVector, y::AbstractVector)

Compute the log-likelihood of the DCReM model given the weighted parameters θ = [θᵒ; θⁱ], the strength sequences s_out/s_in and the per-node binary fitnesses x, y, computing the marginal edge probability fᵢⱼ = xᵢyⱼ/(1 + xᵢyⱼ) on the fly.

The (conditional) log-likelihood is $\mathcal{L} = -\sum_i \left[ θ^{o}_i s^{out}_i + θ^{i}_i s^{in}_i \right] + \sum_{i≠j} f_{ij} \log(θ^{o}_i + θ^{i}_j)$. The log is domain-guarded (returns NaN for θᵒᵢ + θⁱⱼ ≤ 0) so that the line search rejects steps that leave the feasible region; this keeps it automatic-differentiation friendly.

This method is used when the marginal probability matrix is not precomputed (i.e. model.Ĝ is nothing).

L_DCReM(θ::AbstractVector, s_out::AbstractVector, s_in::AbstractVector, f::AbstractMatrix)

Compute the log-likelihood of the DCReM model, using the precomputed marginal probability matrix f (i.e. model.Ĝ). See also L_DCReM(::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector).

L_DCReM(m::DCReM)

Return the log-likelihood of the DCReM model m based on the computed maximum likelihood parameters. Depending on the status of the model, the precomputed marginal probability matrix (m.Ĝ) is used, or the marginal probability is computed on the fly.

See also L_DCReM(::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector, ::AbstractVector).

∇L_DCReM!(∇L::AbstractVector, θ::AbstractVector, s_out::AbstractVector, s_in::AbstractVector, x::AbstractVector, y::AbstractVector)

Compute the gradient of the log-likelihood of the DCReM model in a non-allocating manner, computing the marginal edge probability on the fly (x, y = per-node binary fitnesses).

The inner loop is branch-free (the spurious j == i self-term is folded out afterwards) so that it vectorises (@simd). ∂L/∂θᵒᵢ = -s^{out}_i + Σⱼ≠ᵢ fᵢⱼ/(θᵒᵢ + θⁱⱼ) and ∂L/∂θⁱᵢ = -s^{in}_i + Σⱼ≠ᵢ fⱼᵢ/(θⁱᵢ + θᵒⱼ).

See also ∇L_DCReM_minus!.

∇L_DCReM!(∇L::AbstractVector, θ::AbstractVector, s_out::AbstractVector, s_in::AbstractVector, f::AbstractMatrix)

Compute the gradient of the log-likelihood of the DCReM model, using the precomputed marginal probability matrix f (i.e. model.Ĝ). The f[i,i] = 0 diagonal makes the self-terms vanish automatically.

∇L_DCReM_minus!(∇L::AbstractVector, θ::AbstractVector, s_out::AbstractVector, s_in::AbstractVector, x::AbstractVector, y::AbstractVector)

Compute minus the gradient of the log-likelihood of the DCReM model (used for the minimisation carried out by Optimization.jl), computing the marginal edge probability on the fly. Non-allocating.

See also ∇L_DCReM!.

∇L_DCReM_minus!(∇L::AbstractVector, θ::AbstractVector, s_out::AbstractVector, s_in::AbstractVector, f::AbstractMatrix)

Compute minus the gradient of the log-likelihood of the DCReM model, using the precomputed marginal probability matrix f (i.e. model.Ĝ). Non-allocating.

DCReM_iter!(θ::AbstractVector, s_out::AbstractVector, s_in::AbstractVector, x::AbstractVector, y::AbstractVector, G::AbstractVector)

Compute the next fixed-point iteration for the DCReM model, computing the marginal edge probability on the fly (x, y = per-node binary fitnesses). The pre-allocated buffer G is updated in place.

The consistency equations are $θ^{o}_i = \left( \sum_{j\neq i} f_{ij}/(1 + θ^{i}_j/θ^{o}_i) \right) / s^{out}_i$ and $θ^{i}_i = \left( \sum_{j\neq i} f_{ji}/(1 + θ^{o}_j/θ^{i}_i) \right) / s^{in}_i$.

DCReM_iter!(θ::AbstractVector, s_out::AbstractVector, s_in::AbstractVector, f::AbstractMatrix, G::AbstractVector)

Compute the next fixed-point iteration for the DCReM model, using the precomputed marginal probability matrix f (i.e. model.Ĝ). The pre-allocated buffer G is updated in place.

set_xᵣ!(m::DCReM)

Set the value of xᵣ to exp(-αᵣ[1:d_unique]) for the (binary layer of the) DCReM model m.

set_yᵣ!(m::DCReM)

Set the value of yᵣ to exp(-αᵣ[d_unique+1:end]) for the (binary layer of the) DCReM model m.

Ĝ(m::DCReM)

Compute the expected (binary) adjacency matrix for the DCReM model m, i.e. fᵢⱼ = xᵢyⱼ/(1 + xᵢyⱼ), so that sum(Ĝ(m), dims=2) ≈ outdegree and sum(Ĝ(m), dims=1) ≈ indegree.

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

set_Ĝ!(m::DCReM)

Set the expected (binary) adjacency matrix for the DCReM model m.

Ŵ(m::DCReM)

Compute the expected (unconditional) weighted adjacency matrix for the DCReM model m, i.e. ⟨wᵢⱼ⟩ = fᵢⱼ / (θᵒᵢ + θⁱⱼ), so that sum(Ŵ(m), dims=2) ≈ outstrength and sum(Ŵ(m), dims=1) ≈ instrength.

set_Ŵ!(m::DCReM)

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

σˣ(m::DCReM)

Compute the standard deviation for the elements of the (binary) adjacency matrix for the DCReM model m, i.e. sqrt(fᵢⱼ(1 - fᵢⱼ)) (the adjacency entries are Bernoulli distributed and, under the conditional DBCM layer, independent).

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::DCReM)

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

σʷ(m::DCReM)

Compute the standard deviation for the elements of the weighted adjacency matrix for the DCReM model m. The (unconditional) weight wᵢⱼ is a mixture of an exponential (with probability fᵢⱼ) and zero, so $⟨w_{ij}^2⟩ = 2f_{ij}/(θ^{o}_i + θ^{i}_j)^2$ and $Var(w_{ij}) = f_{ij}(2 - f_{ij})/(θ^{o}_i + θ^{i}_j)^2$.

Under the conditional DBCM layer the weights of distinct ordered pairs are independent (Cov(wᵢⱼ, wⱼᵢ) = 0), unlike for the CRWCM.

set_σʷ!(m::DCReM)

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

precision(m::DCReM)

Determine the compute precision of the DCReM model m.

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

Return the expected value of the (binary) adjacency matrix for the DCReM model m at the node pair (i,j) (fᵢⱼ = xᵢyⱼ/(1 + xᵢyⱼ)).

❗ For performance reasons, the function does not check:

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

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

Compute the standard deviation of metric X for the DCReM model m via error propagation (the delta method of Squartini & Garlaschelli (2011)).

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_σʷ!). This is the layer to use for weighted metrics such as the triadic fluxes.
• gradient_method::Symbol: :ForwardDiff, :ReverseDiff (default) or :Zygote.

Under the conditional DBCM binary layer all matrix entries are independent (no covariance terms).