UECM

Maximum entropy model for the Undirected Enhanced Configuration Model (UECM).

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

The UECM constrains both the degree sequence and the (integer) strength sequence of an undirected weighted network.

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

UECM(G::T; d::Vector, s::Vector, precision::Type{<:AbstractFloat}=Float64, kwargs...) where {T}

Constructor function for the UECM type.

By default and depending on the graph type T, the definition of degree from $Graphs.jl$ and strength from $SimpleWeightedGraphs$ is applied. If you want to use a different definition of degree or strength, you can pass the vectors as keyword arguments.

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

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

Examples

# generating a model from a weighted graph (here: the symmetrised rhesus macaques network)
julia> G = MaxEntropyGraphs.SimpleWeightedGraphs.SimpleWeightedGraph(MaxEntropyGraphs.rhesus_macaques());

julia> model = UECM(G)
UECM{SimpleWeightedGraphs.SimpleWeightedGraph{Int64, Float64}, Float64} (16 vertices, 16 unique {degree,strength} pairs, 1.00 compression ratio)

# generating a model directly from a degree and strength sequence
julia> model = UECM(d=[1, 2, 2, 1], s=[3, 5, 4, 2])
UECM{Nothing, Float64} (4 vertices, 4 unique {degree,strength} pairs, 1.00 compression ratio)

# generating a model with a different precision
julia> model = UECM(d=[1, 2, 2, 1], s=[3, 5, 4, 2], precision=Float32);

julia> MaxEntropyGraphs.precision(model)
Float32

The degree is taken from Graphs.degree and the strength from MaxEntropyGraphs.strength (the SimpleWeightedGraphs weighted-degree).

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

Compute the likelihood maximising parameters of the UECM model m.

By default the parameters are computed using the BFGS method with the strength sequence 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::UECM; method::Symbol=:strengths)

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

The methods available are:

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

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

Generate a random weighted graph from the UECM model m.

Arguments

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

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

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

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

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

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

Compute the Akaike Information Criterion (AIC) for the UECM model m. The parameters of the model must be computed beforehand. The UECM has 2N parameters (an $\alpha$ and a $\beta$ per node).

See also AICc, L_UECM_reduced.

AICc(m::UECM)

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

See also AIC, L_UECM_reduced.

BIC(m::UECM)

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

See also AIC, L_UECM_reduced.

Return the reduced number of {degree, strength} pairs in the UECM network

L_UECM_reduced(θ::Vector, d::Vector, s::Vector, F::Vector, n::Int=length(d))

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

Arguments

• θ: the maximum likelihood parameters of the model (θ = [α; β])
• d: the reduced degree sequence
• s: the reduced strength sequence
• F: the frequency of each (degree, strength) pair
• n: the number of unique (degree, strength) pairs (defaults to length(d))

The function is numerically stabilised (expm1/log1p) so that the 1 - exp(-βᵢ-βⱼ) denominator does not suffer from catastrophic cancellation, and stays automatic-differentiation friendly.

Examples

julia> θ = [1.0, 2.0, 3.0, 4.0, 1.5, 2.5];

julia> d = [1, 2, 1]; s = [2, 3, 1]; F = [1, 1, 1];

julia> L_UECM_reduced(θ, d, s, F);

L_UECM_reduced(m::UECM)

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

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

∇L_UECM_reduced!(∇L::AbstractVector, θ::AbstractVector, d::Vector, s::Vector, F::Vector, x::AbstractVector, y::AbstractVector, n=length(θ)÷2)

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

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

∇L_UECM_reduced_minus!(args...)

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

See also ∇L_UECM_reduced!.

UECM_reduced_iter!(θ, d, s, F, x, y, G, nz, n=length(θ)÷2)

Compute the next fixed-point iteration for the reduced UECM model. The pre-allocated buffers x, y and G are updated in place. Only the non-zero constraints nz are iterated.

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

set_xᵣ!(m::UECM)

Set the value of xᵣ to exp(-αᵣ) for the UECM model m

set_yᵣ!(m::UECM)

Set the value of yᵣ to exp(-βᵣ) for the UECM model m

Ĝ(m::UECM)

Compute the expected adjacency matrix for the UECM model m.

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

set_Ĝ!(m::UECM)

Set the expected adjacency matrix for the UECM model m

Ŵ(m::UECM)

Compute the expected (unconditional) weighted adjacency matrix for the UECM model m, i.e. ⟨wᵢⱼ⟩ = pᵢⱼ / (1 - yᵢyⱼ), so that sum(Ŵ(m), dims=2) ≈ strength.

set_Ŵ!(m::UECM)

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

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

Compute the standard deviation for the elements of the (binary) adjacency matrix for the UECM 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::UECM)

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

σʷ(m::UECM)

Compute the standard deviation for the elements of the weighted adjacency matrix for the UECM model m. The weight wᵢⱼ follows a Bernoulli–geometric mixture: with pᵢⱼ the connection probability and yᵢyⱼ 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 undirected, wᵢⱼ and wⱼᵢ denote the same random variable; the corresponding covariance is accounted for by σₓ.

set_σʷ!(m::UECM)

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

precision(m::UECM)

Determine the compute precision of the UECM model m.

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

Return the expected value of the adjacency matrix for the UECM model m at the 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_UECM(xixj::T, yiyj::T) where {T}

Helper for the UECM model computing the expected adjacency entry pᵢⱼ = (xᵢxⱼ·yᵢyⱼ)/(1 - yᵢyⱼ + xᵢxⱼ·yᵢyⱼ) from the products xᵢxⱼ and yᵢyⱼ of the maximum likelihood parameters.

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

Compute the standard deviation of metric X for the UECM 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 undirected, the entries (i,j) and (j,i) of either layer denote the same random variable, so the delta method over ordered pairs includes the within-dyad covariance cross-term (Cov(g_ij, g_ji) = σ²[g_ij]). This makes the result independent of whether X is written using one or both triangles of the matrix.

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.