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). On this model it bounds the constraint residual — the largest absolute mismatch between an expected and an observed reduced degree or strength — 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 check any solve.
• 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 (UECM_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 networks it converged on 0. The replacement solves each block exactly and is feasible by construction; on the 100 of those networks that are well posed it converges on 98, against 86 for :BFGS, and it is both the most accurate and much the cheapest path — on a 250-vertex weighted network it reaches a degree/strength residual of 9·10⁻⁹ in 0.28 s, against 1.4·10⁻⁵ in 57 s for :BFGS. On the 50 networks carrying a runaway constraint (see Conditioning below) it reaches ftol on only 8; there the optimum is at an infinite parameter and :BFGS (38/50) is the method to use — solve_model! says so, naming the offending constraint, rather than failing silently. Such a solve runs the full maxiters before giving up, so lower it if you want the answer sooner. The legacy map is still exported as UECM_reduced_iter!.
• the L-BFGS method is known to be unstable for this model and should not be used.
• :Newton requires the default initial = :strengths. The UECM likelihood lives on the open box βᵢ > 0; the first-order methods carry that box explicitly (via Fminbox), but Fminbox does not accept Newton, so :Newton solves the box-constrained problem without box protection. From :strengths that is fine (residuals 1e-9–1e-14 on every network tested), but from a far start — :uniform sets β = 6.9, i.e. y ≈ 1e-3, against an optimum near β ≈ 8e-3 — the Newton step overshoots past β = 0, where the objective is NaN, and the solve aborts with a ConvergenceError. That honest failure is deliberate: IPNewton with the box removes the failures but reports success far from the optimum instead (measured off by up to 5·10³ in -L), which is worse. See validation/bicm_uecm_solver_geometry.md.

The feasible set

The UECM is defined wherever every pair that exists has yᵢyⱼ < 1, i.e. βᵢ + βⱼ > 0 — a pairwise condition. The only per-coordinate case is the same-class pair, needing 2βᵢ > 0, and a class of multiplicity Fᵢ = 1 has no same-class pair. The Fminbox used by the first-order methods therefore floors only the βᵢ of classes with Fᵢ ≥ 2 (see MaxEntropyGraphs._UECM_β_FLOOR); flooring every class, as v0.7.0 did, pinned singleton classes against a wall that is not part of the model and reported Success from it (measured: 5 of 128 networks, degree/strength residual 0.17–2.31). :fixedpoint carries the exact feasible set and is unaffected either way.

Conditioning

Unlike the DECM, the UECM has no gauge freedom: its pair terms couple αᵢ + αⱼ within one block, so a shift adds 2c rather than cancelling, and the Hessian has no null directions. On clean data it is well conditioned (κ ≈ 250–7000). It does share the DECM's runaway constraints, though: a saturated degree (k = n-1) drives α → -∞, and a node all of whose links carry weight 1 (s = k) drives β → +∞. In the extreme where every weight is 1, s = k holds for every node, the whole β block is unidentifiable and the model degenerates to a UBCM — the strength sequence then carries no information beyond the degree sequence.

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

UECM_reduced_coordinate_iter!(θ, d, s, F, x, y, live, n=length(θ)÷2)

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

Why not the Picard recipe

The classical recipe (still available as UECM_reduced_iter!) writes ⟨kᵢ⟩ = xᵢ·Aᵢ and ⟨sᵢ⟩ = yᵢ·Bᵢ and updates xᵢ ← dᵢ/Aᵢ, yᵢ ← sᵢ/Bᵢ with Aᵢ, Bᵢ frozen at the previous iterate. With g = xᵢxⱼ, t = yᵢyⱼ and D(t) = (1-t)(1-t+gt):

• Aᵢ is strictly decreasing in xᵢ (∂/∂xᵢ [c/(d₀+cxᵢ)] = -c²/(d₀+cxᵢ)² < 0), so the degree step undershoots: it always lands in (xᵢ, xᵢ*). It is safe.
• Bᵢ = Σⱼ wⱼ gⱼyⱼ/D(yᵢyⱼ) is strictly increasing in yᵢ whenever D'(0) = xᵢxⱼ - 2 < 0, which holds throughout the sparse regime, so the strength step overshoots — and without bound: as yᵢ → 0, Bᵢ → xᵢΣⱼ wⱼxⱼyⱼ, giving yᵢ ← sᵢ/(xᵢΣⱼ wⱼxⱼyⱼ) = O(1/(x²y)).

The model is only defined for yᵢyⱼ < 1, and nothing in that recipe enforces it. From the :strengths cold start on the symmetrised rhesus network the strength step overshoots its own root by a factor 208–1383 on every class, leaving the domain on the first iteration; the optimum there sits at max yᵢyⱼ = 0.974, so there is no slack to absorb it. Measured over 150 random weighted networks the Picard/Anderson path converged on 0.

What this map does instead

Each class is solved exactly, one block at a time, which is possible because both blocks are monotone:

• ⟨kᵢ⟩(xᵢ) increases from 0 to Σⱼ wⱼ, so a root exists iff dᵢ is below that ceiling (it is not for a saturated degree, which is a genuine runaway) and is unique;
• ⟨sᵢ⟩(yᵢ) increases from 0 to ∞ on (0, ȳᵢ) — d/dt [gt/D] = g(1-(1-g)t²)/D² > 0 for every t ∈ (0,1) and g > 0 — so a root always exists, is unique, and is feasible by construction.

A two-dimensional Newton polish (_uecm_polish!) then follows runaway ridges that the alternating solves only crawl along. Both steps are exact block maximisations of a concave objective, so the sweep increases the likelihood monotonically and can never leave the domain.

See also _monotone_root, UECM_reduced_iter! and validation/uecm_decm_fixedpoint.md.

_uecm_polish!(i, x, y, d, s, F, live, ȳ)

Safeguarded two-dimensional Newton step on class i's own pair of constraints, in (log xᵢ, log yᵢ).

The alternating one-dimensional solves converge for a well-posed problem but crawl along a runaway ridge — a class with sᵢ = dᵢ (all its links carry weight 1) has its optimum at xᵢ → ∞, yᵢ → 0 jointly, a direction no single coordinate move can follow; measured rate 0.9998 per sweep, i.e. some 7·10⁴ sweeps to reach 10⁻⁹. Solving the class's two equations together follows that ridge directly.

Every trial point must stay inside the domain (log yᵢ < log ȳᵢ) and reduce this node's residual. That per-node test is deliberately local — it is what lets the step travel far along a ridge — so it is not a convergence guarantee on its own: near a runaway the node Jacobian is nearly singular and a step that helps this node can wreck the others. The caller therefore accepts the whole polish pass only when it lowers the global residual (see UECM_reduced_coordinate_iter!).

_UECM_β_FLOOR

Lower bound imposed on the β entries of degree/strength classes of multiplicity Fᵢ ≥ 2 when the UECM is solved with a first-order method.

The UECM is defined wherever every pair that actually exists satisfies yᵢyⱼ < 1, i.e. βᵢ + βⱼ > 0 — a pairwise condition, not a per-coordinate one. The only per-coordinate case is the same-class pair, which needs 2βᵢ > 0; a class of multiplicity Fᵢ = 1 contains a single vertex and therefore has no same-class pair, so nothing in the model bounds the sign of its own βᵢ.

Applying this floor to every class — as versions up to v0.7.0 did — is therefore too tight, and it is not a harmless over-restriction: on 5 of 128 random weighted networks the ML optimum genuinely has a βᵢ < 0 on a singleton class, and a solver pinned at the floor reported Success with a degree/strength residual between 0.17 and 2.31. Singleton classes now get -Inf and are held inside the domain by the objective itself (L_UECM_reduced returns NaN off it), exactly as the DECM does.

The value is deliberately tiny — a feasibility guard, not a regularisation, so it does not move a well-posed solution. A class whose ML solution genuinely pushes βᵢ to the boundary will report that βᵢ sitting at this floor.

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.