CReM

Maximum entropy model for the Conditional Reconstruction Method (CReM).

The CReM is a two-step null model for weighted, undirected networks with continuous, positive weights. The binary structure (topology) is supplied by a prior binary model — here an internally solved UBCM on the degree sequence, giving the marginal edge probability fᵢⱼ = xᵢxⱼ/(1 + xᵢxⱼ) with xᵢ = e^{-αᵢ}. Conditional on an edge existing, the weight follows an exponential distribution with rate θᵢ + θⱼ (mean 1/(θᵢ + θⱼ)); the parameters θ (one per node) constrain the strength sequence.

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

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

Constructor function for the CReM 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 CReM allows continuous, positive weights (the strength sequence need not be integer-valued).

Examples

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

julia> model = CReM(G)
CReM{SimpleWeightedGraphs.SimpleWeightedGraph{Int64, Float64}, Float64} (16 vertices, 16 unique degrees, 1.00 compression ratio)

# generating a model directly from a degree and strength sequence
julia> model = CReM(d=[1, 2, 2, 1], s=[3.0, 5.0, 4.0, 2.0])
CReM{Nothing, Float64} (4 vertices, 3 unique degrees, 0.75 compression ratio)

# generating a model with a different precision
julia> model = CReM(d=[1, 2, 2, 1], s=[3.0, 5.0, 4.0, 2.0], precision=Float32);

julia> MaxEntropyGraphs.precision(model)
Float32

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

Compute the likelihood maximising parameters of the CReM model m. This is a two-step process: first the binary (conditional) UBCM layer is solved on the degree sequence, then the weighted CReM layer is solved on the strength sequence.

By default the weighted layer is computed with the fixed-point method using the strength sequence as initial guess (the CReM fixed-point recipe is stable, unlike the UECM's).

Arguments (weighted CReM 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 UBCM 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.

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

Compute an initial guess θ₀ for the maximum likelihood parameters of the weighted (CReM) layer of the model m.

The CReM 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ᵢ > 0] / Σⱼ sⱼ.
• :strengths_minor: θᵢ = 𝟙[sᵢ > 0] / (sᵢ + 1).
• :random: random values drawn from $U(0,1)$.

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

Generate a random weighted graph from the CReM model m.

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

Arguments

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

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

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

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

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

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

Compute the Akaike Information Criterion (AIC) for the CReM model m. The parameters of the model must be computed beforehand. The (conditional) CReM has N parameters (one weighted parameter $θ$ per node; the binary layer's fᵢⱼ are fixed inputs of the conditional likelihood).

See also AICc, L_CReM.

AICc(m::CReM)

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

See also AIC, L_CReM.

BIC(m::CReM)

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

See also AIC, L_CReM.

Return the number of vertices in the CReM network

L_CReM(θ::AbstractVector, s::AbstractVector, f::AbstractVector)

Compute the log-likelihood of the CReM model given the weighted parameters θ, the strength sequence s and the per-node binary fitness f (fᵢ = xᵢ), computing the marginal edge probability fᵢⱼ = fᵢfⱼ/(1 + fᵢfⱼ) on the fly.

The log-likelihood is $\mathcal{L} = -\sum_i θ_i s_i + \sum_{i} \sum_{j<i} f_{ij} \log(θ_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 (used for the AD cross-check).

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

L_CReM(θ::AbstractVector, s::AbstractVector, f::AbstractMatrix)

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

L_CReM(m::CReM)

Return the log-likelihood of the CReM 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_CReM(::AbstractVector, ::AbstractVector, ::AbstractVector).

∇L_CReM!(∇L::AbstractVector, θ::AbstractVector, s::AbstractVector, f::AbstractVector)

Compute the gradient of the log-likelihood of the CReM model in a non-allocating manner, computing the marginal edge probability on the fly (f = per-node binary fitness).

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

See also ∇L_CReM_minus!.

∇L_CReM!(∇L::AbstractVector, θ::AbstractVector, s::AbstractVector, f::AbstractMatrix)

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

∇L_CReM_minus!(∇L::AbstractVector, θ::AbstractVector, s::AbstractVector, f::AbstractVector)

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

See also ∇L_CReM!.

∇L_CReM_minus!(∇L::AbstractVector, θ::AbstractVector, s::AbstractVector, f::AbstractMatrix)

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

CReM_iter!(θ::AbstractVector, s::AbstractVector, f::AbstractVector, G::AbstractVector)

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

The consistency equation is $θ_i = \left( \sum_{j\neq i} f_{ij}/(1 + θ_j/θ_i) \right) / s_i$.

CReM_iter!(θ::AbstractVector, s::AbstractVector, f::AbstractMatrix, G::AbstractVector)

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

set_xᵣ!(m::CReM)

Set the value of xᵣ to exp(-αᵣ) for the (binary layer of the) CReM model m.

Ĝ(m::CReM)

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

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

set_Ĝ!(m::CReM)

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

Ŵ(m::CReM)

Compute the expected (unconditional) weighted adjacency matrix for the CReM model m, i.e. ⟨wᵢⱼ⟩ = fᵢⱼ / (θᵢ + θⱼ), so that sum(Ŵ(m), dims=2) ≈ strength.

set_Ŵ!(m::CReM)

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

σˣ(m::CReM)

Compute the standard deviation for the elements of the (binary) adjacency matrix for the CReM model m, i.e. sqrt(fᵢⱼ(1 - fᵢⱼ)) (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::CReM)

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

σʷ(m::CReM)

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

$Var(w_{ij}) = \frac{f_{ij}(2 - f_{ij})}{(θ_i + θ_j)^2}$

As the network is undirected, wᵢⱼ and wⱼᵢ denote the same random variable; the corresponding covariance is accounted for by σₓ.

set_σʷ!(m::CReM)

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

precision(m::CReM)

Determine the compute precision of the CReM model m.

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

Return the expected value of the (binary) adjacency matrix for the CReM 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_CReM(xixj::T) where {T}

Helper for the CReM model computing the (binary) expected adjacency entry fᵢⱼ = xᵢxⱼ/(1 + xᵢxⱼ) from the product xᵢxⱼ of the binary maximum likelihood parameters.

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

Compute the standard deviation of metric X for the CReM 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.