L0ggm

Direct $L_0$ minimization requires searching over all $2^{p(p-1)/2}$ subsets of candidate edges. For even modest network sizes (e.g., $p = 20$ implies 190 candidate edges), this combinatorial search is computationally intractable.

{L0ggm} addresses this through continuous smooth approximations $\rho(x; \lambda, \gamma)$ that closely mimic the $L_0$ indicator and satisfy three key properties:

1. $\rho(0; \lambda, \gamma) = 0$ — zero penalty at zero
2. $\rho(x; \lambda, \gamma) \to \lambda$ as $|x| \to \infty$ — bounded, like $L_0$
3. $\rho(x; \lambda, \gamma) \to \lambda \cdot \mathbf{1}(x \neq 0)$ as $\gamma \to 0$ — converges to the $L_0$ step function

Because the approximations are differentiable, their derivatives can be used as adaptive, element-wise GLASSO penalty weights via the Local Linear Approximation (LLA; Fan & Li, 2001; Zou & Li, 2008). At iteration $t$, the penalty is linearized around the current estimate $\mathbf{K}^{(t)}$:

$$\rho(K_{ij}) \approx \rho(K_{ij}^{(t)}) + \rho'(K_{ij}^{(t)}) \cdot (K_{ij} - K_{ij}^{(t)})$$

which reduces each step to a re-weighted GLASSO problem. A single LLA pass (the default) is extremely fast and produces estimates with strong theoretical guarantees (Zou & Li, 2008). Full iterative LLA to convergence is also available via LLA = TRUE.

The five penalties available in {L0ggm} are:

Note. The Gumbel penalty subtracts $e^{-1}$ from the raw Gumbel CDF to shift the y-intercept to zero, then scales by $\lambda / (1 - e^{-1})$ so that the penalty ceiling is $\lambda$, consistent with the other approximations. The Weibull formula shown applies for $k \leq 1$; for $k > 1$ the penalty is piecewise — linear up to the mode $x^*$ and a continuity-shifted Weibull CDF above it — matching the piecewise derivative described below.

Figure 1. $L_0$ norm approximation penalties as a function of coefficient magnitude. Solid line: $L_0$ norm (step function). Dashed lines: $L_1$ norm (LASSO) and each continuous approximation penalty implemented in {L0ggm}. Gumbel is roughly equivalent to EXP (not pictured) and Weibull k = 1 is exactly EXP.

What makes it different

The Weibull penalty introduces a shape parameter $k > 0$ that continuously interpolates between qualitatively different penalty regimes:

• $k = 1$: reduces exactly to the EXP penalty
• $k < 1$: the penalty becomes less concave than exponential, with a wider rise near the origin and slower convergence to $\lambda$ — moving away from the $L_0$ step function (see Figure 1)
• $k \sim 0.40$: approximates the Atan penalty
• $k > 1$: the penalty becomes more concave than exponential, with a sharper threshold near the scale — moving closer to the $L_0$ step function

This variation in $k$ allows Weibull to control the degree of sparsity bias: smaller $k$ imposes more shrinkage on nonzero coefficients and less aggressively zeros small ones. Rather than fixing $k$ by hand, {L0ggm} estimates it from the data.

Adaptive parameter estimation

When adaptive = TRUE (the default), both the shape $k$ and scale $\gamma$ of the Weibull penalty are calibrated to the observed partial correlation distribution, making the penalty data-driven rather than fixed. The same adaptive principle applies to the EXP and Gumbel penalties.

Step 1 — Fit Weibull to empirical partial correlations. Let ${|p_{ij}|}$ denote the lower-triangular absolute partial correlations. Maximum likelihood estimates $\hat{k}$ and $\hat{\lambda}_W$ are obtained by solving:

$$\frac{\sum_i |p_i|^{\hat{k}} \log|p_i|}{\sum_i |p_i|^{\hat{k}}} - \frac{1}{\hat{k}} - \frac{1}{n}\sum_i \log|p_i| = 0$$

$$\hat{\lambda}_W = \left(\frac{1}{n}\sum_i |p_i|^{\hat{k}}\right)^{1/\hat{k}}$$

Step 2 — Set adaptive scale. The scale parameter $\gamma$ is set to the 10th percentile of the fitted Weibull distribution:

$$\gamma = \hat{\lambda}_W \cdot \left(-\log(0.90)\right)^{1/\hat{k}}$$

This is simply the Weibull quantile function evaluated at $p = 0.10$: the value below which 10% of the fitted partial correlations fall, serving as a data-driven threshold between noise and signal.

The choice of the 10th percentile is grounded in a concrete calibration anchor. In practice, the Weibull MLE of absolute partial correlations tends to produce shape estimates near $\hat{k} \approx 1$ and scale estimates near $\hat{\lambda}_W \approx 0.10$. Substituting these values:

$$\gamma \approx 0.10 \cdot (-\log 0.90)^{1} \approx 0.10 \times 0.105 \approx 0.01$$

This matches $\gamma = 0.01$, the fixed default known to work well with the EXP penalty, providing a principled empirical justification for the 10th percentile rule: at typical signal levels, it recovers the well-validated fixed default, while adapting to the actual distribution of partial correlations in data where signal strength deviates from this average. The same 10th-percentile rule is applied to EXP (via its fitted exponential scale) and Gumbel (via its fitted Gumbel quantile function) when those penalties are used adaptively.

The derivative used in the LLA is defined piecewise to ensure the penalty weight is always non-increasing in $|x|$. When $k \leq 1$ the Weibull PDF is already monotonically decreasing, so the standard formula applies directly. When $k > 1$ the raw PDF has an interior mode at $x^* = \gamma \left(\tfrac{k-1}{k}\right)^{1/k}$; values below this mode would otherwise receive less weight than values above it, violating the oracle-property requirement. The derivative is therefore capped at its peak for $|x| \leq x^*$:

where $x^{\ast} = \gamma \left(\dfrac{k-1}{k}\right)^{1/k}$ is the mode of the Weibull PDF. The resulting derivative is monotonically non-increasing for all $k > 0$: as $|x|$ grows past the peak, weights decay to zero — large true edges receive vanishingly small additional penalization, directly addressing the magnitude bias of $L_1$ methods.

Distributional foundations and extreme value theory

The EXP, Weibull, and Gumbel penalties are not arbitrary constructions — they form a mathematically unified family rooted in extreme value theory. The table below shows each penalty alongside its LLA derivative, which reveals the distributional structure underlying each:

The Weibull as the general case. The EXP penalty is the exact special case $k = 1$ of the Weibull — its derivative, its adaptive scale estimation, and its convergence behavior are all inherited from the Weibull family. The Weibull strictly generalizes EXP across the full range $k \in (0, \infty)$.

The Gumbel–Weibull connection. Less obviously, the Gumbel and Weibull are also deeply related. The table shows that the Gumbel survival function $1 - G(x; \gamma)$ equals the Weibull CDF evaluated at the log-transformed argument $e^{-x}$, with shape $k = 1/\gamma$. Equivalently, if $X$ follows a Weibull distribution, then $-\log X$ follows a Gumbel distribution. This log-linear relationship is the classical link between the Type I and Type III extreme value families, and it means the Gumbel penalty is a reparametrized Weibull penalty under a change of variable — a further instance of Weibull's generality.

Generalized extreme value (GEV) distribution. The reason these distributions co-occur is the Fisher–Tippett–Gnedenko theorem (Fisher & Tippett, 1928; Gnedenko, 1943): the normalized maximum of $n$ i.i.d. random variables converges in distribution, as $n \to \infty$, to the generalized extreme value (GEV) family:

$$G(x; \mu, \sigma, \xi) = \exp\left\lbrace -\left[1 + \xi\left(\frac{x - \mu}{\sigma}\right)\right]^{-1/\xi}\right\rbrace$$

The tail index $\xi$ selects among three qualitatively distinct limiting types:

• $\xi = 0$ (Gumbel, Type I): the limit of exponential-class tails; $G(x) = \exp\left(-e^{-(x-\mu)/\sigma}\right)$
• $\xi > 0$ (Fréchet, Type II): heavy-tailed, unbounded support
• $\xi < 0$ (Weibull, Type III): bounded support, arising from distributions with a finite right endpoint

Only Type I (Gumbel) and Type III (Weibull) are used as penalties in {L0ggm}. This is not coincidental: these are precisely the two GEV types whose CDFs are bounded above and approach 1 as $x \to \infty$. The Fréchet CDF grows without bound, violating the fundamental $L_0$ approximation requirement that $\rho(x) \to \lambda$ as $|x| \to \infty$, and is therefore unsuitable as a penalty.

Why extreme value CDFs are ideal $L_0$ approximations. The $L_0$ penalty is a Heaviside step function: zero below a threshold, constant above. Extreme value CDFs share exactly this structural character — they are bounded, concave, and rise from 0 to a finite ceiling — because they represent the limiting distribution of extreme order statistics under the most general conditions. This is no coincidence from the perspective of sparsity: recovering sparse structure amounts to discriminating between the very smallest coefficients (noise) and all others (signal), which is precisely the problem that extreme value distributions are designed to describe.

The Weibull CDF in particular describes the distribution of the minimum of a large sample from distributions with polynomial density near zero. Because partial correlations are bounded in $[0, 1]$, the relevant domain of the penalty lies in the range $[0, \hat{\lambda}_W]$ where $\hat{\lambda}_W < 1$. In this regime, the shape parameter $k$ controls the sharpness of the threshold at the scale $\gamma$: for $x < \gamma$, raising to a larger power $k$ shrinks $(x/\gamma)^k$ toward zero (since $x/\gamma < 1$), while for $x > \gamma$, the same large $k$ amplifies $(x/\gamma)^k$ toward infinity. As $k \to \infty$, the Weibull CDF converges pointwise to $\lambda \cdot \mathbf{1}(x > \gamma)$ — a step function at the scale — the $L_0$ indicator itself.

The same structure governs the LLA derivative. For $k > 1$, the piecewise derivative holds a constant plateau near the origin (up to the mode $x^*$) and then decays — producing a consistent penalty on all near-zero coefficients and vanishing penalization on large ones. This is the continuous analog of $L_0$ hard thresholding: aggressively zero small edges, leave large edges essentially unpenalized. For $k < 1$, the derivative decays from a very large value at the origin, imposing the largest penalty weights on the smallest coefficients — an increasingly proportional, $L_1$-like profile. The adaptive fitting procedure estimates $k$ from the data without restriction; in dense or moderate-signal settings it will select $k > 1$, driving the penalty toward the $L_0$ limit, while sparse settings with many near-zero partial correlations tend to produce $k < 1$.

Generates GGM data from a small-world topology via a three-stage process closely aligned with empirical psychometric network properties.

Stage 1 — Ring lattice. A ring lattice is constructed with neighbors + 1 nearest-neighbor connections per node, where neighbors is derived from nodes and density. The extra neighbor ensures the lattice always has more edges than the target, guaranteeing pruning can proceed.

Stage 2 — Pruning. The lattice is randomly pruned to the target density by removing edges uniformly at random while maintaining graph connectivity. This creates degree heterogeneity before any rewiring, providing a non-uniform prior for the next stage.

Stage 3 — Degree-weighted rewiring. Each edge is independently rewired with probability rewire. Unlike standard Watts-Strogatz, the new endpoint is selected proportionally to $\sqrt{d_i + d_k}$ (combined degree of the candidate pair), moderating the extreme hubs of linear preferential attachment while still producing realistic degree heterogeneity.

Edge weights. Retained lattice edges receive weights by neighbor-distance priority (nearer = stronger); rewired edges receive lowest priority. Absolute weights are drawn from a Weibull distribution via weibull_parameters().

Returns: a named list with data, parameters (including omega smallworldness), population (R, Omega), and convergence.

result <- simulate_smallworld(
  nodes = 20, density = 0.30, rewire = 0.20, sample_size = 500
)

Generates GGM data with explicit community (block) structure, where edge density is controlled separately within and between communities via a blocks × blocks density matrix.

Structure. Nodes are partitioned into blocks communities of size nodes (scalar or per-block vector). For each pair of communities $(i, j)$, edges are included independently with probability density_matrix[i, j]. The resulting graph is required to be connected; disconnected draws are rejected and resampled.

Edge weights and diffusion. Absolute weights are drawn from a Weibull distribution via weibull_parameters(). The diffusion parameter controls the minimum fraction of high-weight edges reserved for within-community positions: 1 - diffusion of within-community slots are filled from the top-ranked weight draws, with the remaining slots drawn randomly. Because random draws can also land within communities by chance, the actual within-community weight advantage typically exceeds 1 - diffusion. Newman-Girvan modularity $Q$ (returned in parameters) summarizes the community contrast actually achieved.

Between-community sign. Each between-community edge is independently signed negative with probability negative_proportion, modeling inhibitory connections. Within-community edges are always positive.

Returns: a named list with data, parameters (including Q modularity and weibull edge weight parameters), population (R, Omega, membership), and convergence.

dm <- matrix(0.20, nrow = 3, ncol = 3)
diag(dm) <- 0.90

result <- simulate_sbm(
  nodes = 6, blocks = 3, density_matrix = dm, sample_size = 500
)

Predicts the Weibull shape $k$ and scale $\lambda_W$ that characterize absolute partial correlation edge weights for a network of given size and sample size. Predictions come from a Seemingly Unrelated Regression (SUR) model fitted to 194 empirical psychometric networks (Huth et al., 2025), where absolute partial correlations were found to follow a Weibull distribution more consistently than Beta, Gamma, or log-normal alternatives.

Predictors. The two SUR equations have an asymmetric structure:

• Shape is predicted from snr and rlp $= 1/\log(p)$, capturing the diminishing marginal effect of network size on weight concentration. Shape governs the spread of the edge weight distribution and is not affected by sampling precision.
• Scale is predicted from snr, rlp, and scaling $= \sqrt{1/(n - p - 2)}$ (the standard error of partial correlations). Scale governs the typical magnitude of edge weights and is directly affected by estimation precision.

Model fit. Shape: $R^2 = 0.887$, RMSE $= 0.048$. Scale: $R^2 = 0.885$, RMSE $= 0.011$. Residual correlation between equations: $0.261$ (motivating SUR). Empirical ranges: shape $\in [0.72, 1.63]$ (M $= 1.07$, SD $= 0.14$); scale $\in [0.03, 0.19]$ (M $= 0.10$, SD $= 0.03$).

# Predicted parameters for a 15-node network, n = 300
weibull_parameters(nodes = 15, sample_size = 300)

# With bootstrapped residuals for Monte Carlo simulation
weibull_parameters(nodes = 15, sample_size = 300, bootstrap = TRUE)

The shared simulation pipeline is implemented in helpers_simulation.R and handles tasks common to both simulate_smallworld() and simulate_sbm().

generate_edges() draws absolute edge weights from a Weibull distribution parameterized by weibull_parameters(). Draws are repeated until no weight exceeds the empirical maximum ($\approx 0.845$), and are then sampled proportionally to $w^2$ (preferring larger weights for occupied edges). The Weibull parameters are restricted to empirical bounds (shape $\in [0.70, 1.70]$, scale $\in [0.03, 0.19]$); draws outside these bounds are discarded.

condition_network() recovers positive definiteness by adding a minimal ridge penalty $\lambda$ to the diagonal of the precision matrix via uniroot(), targeting a user-specified condition number. Maximum shrinkage is capped at approximately 23% following Peeters et al. (2020).

simulate_data() generates multivariate normal observations via a Cholesky decomposition of the population correlation matrix $\mathbf{R}$, then applies variable-level skew transformations. Skew values are rounded to the nearest 0.05 increment within $[-2, 2]$; if skew_range is provided, per-variable skew is drawn uniformly from the interval.