Paper deep dive
Geometry-Aware Probabilistic Circuits via Voronoi Tessellations
Sahil Sidheekh, Sriraam Natarajan
Intelligence
Status: succeeded | Model: google/gemini-3.1-flash-lite-preview | Prompt: intel-v1 | Confidence: 94%
Last extracted: 3/22/2026, 6:29:34 AM
Summary
The paper introduces 'Geometry-Aware Probabilistic Circuits' (PCs) by incorporating Voronoi Tessellations (VT) to enable local geometric structure modeling. It addresses the resulting intractability of inference by proposing two solutions: a certified approximate inference framework using axis-aligned box bounds and a structural condition (Hierarchical Factorized Voronoi) for exact inference. The approach includes a differentiable relaxation for gradient-based learning.
Entities (4)
Relation Signals (3)
Voronoi Tessellations → integratedinto → Probabilistic Circuits
confidence 95% · We propose Voronoi tessellations (VT) as a natural way to incorporate geometric structure directly into the sum nodes of a PC.
Certified Approximate Inference → provides → Lower and Upper Bounds
confidence 95% · provides guaranteed lower and upper bounds for inference
Hierarchical Factorized Voronoi → enables → Exact Tractable Inference
confidence 90% · a structural condition for VT under which exact tractable inference is recovered.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Probabilistic circuits (PCs) enable exact and tractable inference but employ data independent mixture weights that limit their ability to capture local geometry of the data manifold. We propose Voronoi tessellations (VT) as a natural way to incorporate geometric structure directly into the sum nodes of a PC. However, naïvely introducing such structure breaks tractability. We formalize this incompatibility and develop two complementary solutions: (1) an approximate inference framework that provides guaranteed lower and upper bounds for inference, and (2) a structural condition for VT under which exact tractable inference is recovered. Finally, we introduce a differentiable relaxation for VT that enables gradient-based learning and empirically validate the resulting approach on standard density estimation tasks.
Tags
Links
- Source: https://arxiv.org/abs/2603.11946v1
- Canonical: https://arxiv.org/abs/2603.11946v1
Trouble viewing inline? Open PDF directly →
Full Text
102,313 characters extracted from source content.
Expand or collapse full text
Geometry-Aware Probabilistic Circuits via Voronoi Tessellations Sahil Sidheekh Sriraam Natarajan Abstract Probabilistic circuits (PCs) enable exact and tractable inference but employ data independent mixture weights that limit their ability to capture local geometry of the data manifold. We propose Voronoi tessellations (VT) as a natural way to incorporate geometric structure directly into the sum nodes of a PC. However, naïvely introducing such structure breaks tractability. We formalize this incompatibility and develop two complementary solutions: (1) an approximate inference framework that provides guaranteed lower and upper bounds for inference, and (2) a structural condition for VT under which exact tractable inference is recovered. Finally, we introduce a differentiable relaxation for VT that enables gradient-based learning and empirically validate the resulting approach on standard density estimation tasks. Machine Learning, ICML 1 Introduction Probabilistic circuits (PCs) have emerged as a powerful class of generative models that can learn and reason about complex data distributions under uncertainty. By enforcing structural properties, PCs enable exact and linear time inference of likelihoods, marginals and conditionals (Darwiche, 2003; Poon and Domingos, 2011; Choi et al., 2020), making them valuable for applications requiring reliable probabilistic reasoning, such as density estimation, out-of-distribution detection (Braun et al., 2025), causal/counterfactual reasoning (Zečević et al., 2021), multimodal fusion (Sidheekh et al., 2025, 2024) and structured prediction, among others. Despite recent advances in building expressive PCs, most existing architectures share an important limitation: the mixture weights associated with sum nodes are typically data-independent. This means that the routing decisions within the circuit are fixed globally and do not adapt to individual inputs. While this design choice is often important for preserving tractability, it restricts the ability of PCs to capture and adapt to the local geometric structure in the data manifold. In many real-world distributions, the underlying structure varies across regions of the input space, exhibiting piecewise behavior or locality that cannot be adequately modeled using globally shared mixture weights. A natural question that arises is: Can we introduce geometry-aware, input dependent routing into PCs while maintaining tractable inference? Voronoi tessellations (VT) (Aurenhammer, 1991) offer an appealing geometric mechanism for such routing and have been studied in the context of other generative models such as normalizing flows (Chen et al., 2022). By partitioning the input space into convex polyhedral regions based on proximity to a learned set of centroids, Voronoi cells provide a principled way to assign the inputs to local experts. Incorporating VT within PCs could thus help achieve geometric interpretability and adaptive routing, making it a promising approach for modeling distributions with spatially varying structure. However, naively incorporating Voronoi-based gating into probabilistic circuits creates a fundamental conflict with tractable inference: Voronoi cells are defined by oblique half-space intersections that couple multiple input dimensions. Computing integrals over such convex polyhedron is #P\#P-hard (Dyer and Frieze, 1988) and do not decompose along the variable partitions encoded by the circuit’s product nodes. Even when expert distributions are fully factorized, marginalization over Voronoi-gated regions cannot be performed recursively. In deep circuits, these geometric constraints compound across layers, making exact inference intractable even for simple queries. In this work, we formalize this incompatibility between Voronoi-based routing and tractable inference in PCs, and develop two complementary strategies to address it. First, we present a certified approximate inference framework that preserves the reliability guarantees of PCs by computing provable lower and upper bounds on partition functions, marginals, and conditionals. We achieve this by replacing Voronoi cells with tractable axis-aligned box approximations and propagating bounds through the circuit. Second, we identify a structural condition under which exact tractable inference can be recovered. By factorizing the Voronoi tessellation in a manner that aligns with the circuit decomposition, we obtain a class of geometry-aware circuits that support exact inference. This construction, which we refer to as Hierarchical Factorized Voronoi (HFV) probabilistic circuits, enforces a shared factorization between the gating mechanism and the expert distributions, thereby restoring recursive integrability. Learning Voronoi tessellated PCs introduces an additional challenge: hard region assignments are typically non-differentiable. To enable end-to-end gradient-based learning, we introduce a soft gating mechanism based on temperature-scaled distance weighting and employ annealing during training. At test time we revert to hard Voronoi assignments, recovering the exact inference guarantees. We prove that this soft-to-hard transition is well-behaved, with exponentially fast convergence as temperature decreases. Overall, we make the following contributions: (1) We present the first geometric approach based on Voronoi tessellations for training PCs; (2) We formalize the incompatibility between Voronoi-based routing and tractable inference in PCs and outline two different solutions: one based on certified approximate inference and the other based on hierarchical factorizations; (3) For both cases and general geometric modeling, we theoretically analyze the properties of the algorithms and identify the potential gaps; (4) Finally, we perform proof-of-concept experiments to validate the effectiveness and efficiency of the algorithms. Next, we present the necessary background and position our work in the context of the related work. We then outline the geometric formulation and present the two different strategies for learning. Finally, we present our experiments before concluding by outlining areas of future research. 2 Background & Related Work We begin by introducing the core concepts behind tractable probabilistic models (Poon and Domingos, 2011; Darwiche, 2003; Kisa et al., 2014; Rahman et al., 2014) collectively known as probabilistic circuits(Choi et al., 2020).Although we focus on continuous random variables =X1,…,XDX=\X_1,…,X_D\ with joint domain Ω⊆ℝD ^D, the definitions naturally extend to discrete and mixed-variable settings. Definition 2.1. A probabilistic circuit C over variables X is a rooted DAG in which each node n is one of the following. A leaf node represents a tractable univariate distribution pn(xi)p_n(x_i) over a single variable XiX_i. A product node computes fn()=∏c∈ch(n)fc()f_n(x)= _c (n)f_c(x). A sum node computes fn()=∑c∈ch(n)πn,cfc()f_n(x)= _c (n) _n,cf_c(x) with πn,c≥0 _n,c≥ 0 and ∑cπn,c=1 _c _n,c=1. The circuit output is f()=fr()f_C(x)=f_r(x) where r is the root. Each node n in a PC is associated with a scope, the set of variables it depends on. Sum nodes represent mixtures while product nodes represent factorizations. Definition 2.2 (Scope). The scope of node n, denoted scope(n)⊆scope(n) , is defined recursively. If n is a leaf over XiX_i then scope(n)=Xiscope(n)=\X_i\. If n is an internal node then scope(n)=⋃c∈ch(n)scope(c)scope(n)= _c (n)scope(c). Exact tractable inference in a PC hinges on enforcing structural constraints that make marginalization compatible with the circuit factorization. Two key structural properties are: Definition 2.3 (Smoothness). A PC is smooth if for every sum node n, all children share the same scope so that scope(c)=scope(c′)scope(c)=scope(c ) for all c,c′∈ch(n)c,c (n). Definition 2.4 (Decomposability). A PC is decomposable if for every product node n with children c1,…,ckc_1,…,c_k the scopes are disjoint so that scope(ci)∩scope(cj)=∅scope(c_i) (c_j)= for all i≠ji≠ j. Smoothness ensures that mixtures are well defined over consistent variable sets, while decomposability ensures that product nodes combine distributions over independent variable sets, which allows integrals to factor. Together, decomposability and smoothness enable efficient exact marginal and conditional inference (Poon and Domingos, 2011; Darwiche, 2003). A third property is: Definition 2.5 (Determinism). A PC is deterministic if for every sum node (s)(s) and input ()(x) at most one of its child has positive output i.e. |c∈ch(s):fc()>0|≤1. |\\,c (s):f_c(x)>0\,\ |≤ 1. Determinism enables efficient MAP inference and yields sparse and interpretable routing behavior, and is achieved, for example in cutset networks (Rahman et al., 2014), via axis aligned splits of the input space. One of the key research themes within the field of PCs is improving their expressivity to match the performance of deep generative models, while preserving tractability (Sidheekh and Natarajan, 2024). This has led to tensorized formulations of PCs (Peharz et al., 2020b, a; Liu et al., 2024; Loconte et al., 2025a; Zhang et al., 2025) that can be scaled to millions of parameters and trained efficiently via parallelized computations on GPUs using backpropagation, similar to deep neural networks. However, achieving expressivity through larger circuits or more mixture components often leads to diminishing returns (Liu et al., 2023a). As the number of mixture components grow, optimization becomes harder, parameter redundancy increases, and improvements in likelihood saturate, even though inference remains tractable. To address this, recent works have explored alternative learning paradigms such as latent variable distillation (Liu et al., 2023a, b), where structural or semantic information about the data manifold, extracted using a more expressive teacher model (often a deep generative model), is used as auxiliary supervisory signal to learn a student PC, guiding the latent variables associated with its sum node to be meaningful. Along similar lines, better regularization (Vergari et al., 2015; Shih et al., 2021; Liu and Van den Broeck, 2021; Dang et al., 2022) and optimization strategies (Suresh et al., 2026; Karanam et al., 2025) have also been proposed to improve generalization. A complementary direction to increase the expressivity of PCs involves relaxing classical assumptions and extending their representational language. This has resulted in hybrid models that integrate PCs with neural components (Correia et al., 2023; Gala et al., 2024), invertible transformations (Sidheekh et al., 2023), or non-standard mixture constructions (Loconte et al., 2024, 2025b). While these models have expanded the representational scope of PCs, they also reveal a recurring challenge: introducing additional dependencies or operations inside the circuit can silently break tractability unless they are carefully aligned with the circuit’s factorization structure. For example, in Sidheekh et al. 2023, tractability is recovered by enforcing the neural transformations to satisfy decomposability. Allowing mixture weights to depend on the input is a natural way to increase model expressivity and enable local specialization, for example, mixture-of-experts (MoE) models (Jacobs et al., 1991; Shazeer et al., 2017), where a learned gating function routes inputs to specialized subnetworks. Within the PC literature, prior approaches have explored similar directions: CSPNs (Shao et al., 2022) parameterize sum-node weights as neural functions of observed features, enabling conditional density estimation over target variables, but sacrificing tractable inference over the conditioning variables. SPQNs (Sharir and Shashua, 2018) introduce quotient nodes to encode conditionals directly, gaining expressive efficiency but restricting tractable marginalization to subsets agreeing with an induced variable ordering (Sharir and Shashua, 2018). Probabilistic neural circuits (PNCs) (Dos Martires, 2024) generalize this idea further by defining the mixing weights as neural functions of ancestor variables, again trading off general tractable marginalization for increased expressivity. While conceptually related and similar in essence, these approaches differ fundamentally from the setting we study – CSPNs condition on external observed features rather than the modeled variables themselves. Thus, routing is driven by auxiliary information that remains fixed during inference, not on the spatial structure of the data being modeled. SPQNs and PNCs impose implicit variable orderings without geometric interpretation, so routing decisions follow graph-theoretic dependencies rather than spatial proximity or geometric regions, and they sacrifice general any order marginalization capabilities of a PC. Similarly, integral circuits and continuous mixtures extend mixture representations using latent variables that are integrated out, rather than inducing explicit geometric partitions of the input space. However, many real-world distributions exhibit strong locality and piecewise structure: different regions of the input space may follow distinct statistical patterns and dependencies. From a modeling perspective, this suggests routing inputs to local experts based on geometry, rather than relying on globally shared mixture weights. We posit that such a geometry-aware routing can offer additional benefits beyond likelihood improvement. It can enable interpretability through explicit regions of responsibility, support editability and knowledge incorporation by modifying local components without retraining the entire model, and is naturally suited for online or continual learning scenarios (Veness et al., 2021) where new regions of the space may appear and require minimal adaptation. This motivates the development of geometry-aware PCs, and we aim to build principled theoretical foundation in this work. 3 Geometry-Aware Probabilistic Circuits A natural way to equip PCs with geometry awareness is to replace the constant (global) sum node weights with geometry-aware gating, so that different expert subcircuits specialize to different regions of the data manifold. Voronoi tessellations (Aurenhammer, 1991) provide a principled mechanism for such routing. Formally, given centroids 1,…,K⊂ℝd\c_1,…,c_K\ ^d, the Voronoi cell of kc_k is defined as Vk=∈ℝd:‖−k‖22≤‖−j‖22∀jV_k\;=\; \u ^d:\ \|u-c_k\|_2^2≤\|u-c_j\|_2^2\ ∀ j \ (1) Voronoi cells partition space into convex polyhedra with disjoint interiors and boundaries of measure zero, defined by the intersection of half spaces. This assigns inputs to regions based on proximity to learned prototypes, naturally capturing spatial structure, and has been successfully used in clustering (Du et al., 1999), density estimation (Polianskii et al., 2022; Marchetti et al., 2023), and likelihood based generative models (Chen et al., 2022). The resulting routing is deterministic, interpretable, and naturally connects to mixture-of-experts formulations. Thus, we define a geometry-aware sum node by gating mixture components using a Voronoi partition defined over the node’s scope. Definition 3.1. A Voronoi-gated sum node (S) over scope SX_S consists of centroids 1,…,K⊂ℝ|S|\c_1,…,c_K\ ^|X_S| inducing Voronoi cells Vk\V_k\, mixture weights π1,…,πK\ _1,…, _K\ with πk≥0 _k≥ 0 and ∑kπk=1 _k _k=1, and child subcircuits p1,…,pK\p_1,…,p_K\ each with scope SX_S. It computes f(S)=∑k=1Kgk(S)πkpk(S)f(x_S)= _k=1^Kg_k(x_S)\, _k\,p_k(x_S), where gk(S)=[S∈Vk]g_k(x_S)=I[x_S∈ V_k] Since the Voronoi cells partition ℝ|S|R^|X_S|, exactly one gate is active for almost every Sx_S, so the node is deterministic up to measure-zero boundaries. However, as we show below, even a single Voronoi-gated sum node can break the factorization of integrals in a PC, making inference intractable. Proposition 3.2 (Single-Layer Intractability). Let f()=∑k=1Kgk()πkpk()f(x)= _k=1^Kg_k(x)\, _k\,p_k(x) be a Voronoi-gated sum node over =X1,…,XDX=\X_1,…,X_D\. Suppose each child is fully factorized, pk()=∏i=1Dpk(i)(xi)p_k(x)= _i=1^Dp_k^(i)(x_i). Then the partition function Z=∫f()=∑k=1Kπk∫Vkpk()Z= f(x)\,dx= _k=1^K _k _V_kp_k(x)\,dx requires integrating pkp_k over Voronoi cells VkV_k, which are convex polytopes with oblique boundaries. In general, ∫Vk∏ipk(i)(xi)d _V_k _ip_k^(i)(x_i)\,dx does not factor into a product of one-dimensional integrals. Remark 3.3. Definition 3.1 places centroids directly in the input space, which corresponds to an identity embedding. We can also use a learned embedding ϕ:ℝ|S|→ℝdφ:R^|X_S| ^d with centroids in ℝdR^d and gates gk(S)=[ϕ(S)∈Vk]g_k(x_S)=I[φ(x_S)∈ V_k]. The negative results shows that the obstruction is the geometry of the gating regions rather than embedding complexity. Exact tractable inference in smooth decomposable PCs rests on a simple recursion: integrals factor at product nodes because scopes are disjoint, and integrals distribute over sums because mixture weights are constant. Voronoi gating disrupts this recursion by introducing cell-restricted integrals of the form ∫Vkpk(S)S _V_kp_k(x_S)\,dx_S. Even if pkp_k factorizes across variables, the region VkV_k generally does not, as its oblique facets couple variables that the circuit attempts to separate. When applied to deep PCs (see appendix), Voronoi gating at multiple scopes induces intersections and projections of such polyhedral constraints across the circuit hierarchy, further preventing the bottom-up factorization needed for exact inference. This is the core incompatibility between geometric routing and circuit factorization. However, geometrically aligning the regions w.r.t a PC’s variable decomposition can help retain tractability, as we show next. Definition 3.4 (Geometric Alignment). Consider a partition S=S1⊔S2X_S=X_S_1 _S_2. A collection of gating regions Rk\R_k\ is aligned w.r.t this partition if each region decomposes as Rk=Rk(1)×Rk(2)R_k=R_k^(1)× R_k^(2) with Rk(i)⊆ℝ|Si|R_k^(i) ^|X_S_i|. Equivalently, membership in RkR_k can be decided independently as S∈Rkx_S∈ R_k if and only if S1∈Rk(1)x_S_1∈ R_k^(1) and S2∈Rk(2)x_S_2∈ R_k^(2). Theorem 3.5. Consider a voronoi gated sum node f(S)=∑k[S∈Rk]πkpk(S)f(x_S)= _kI[x_S∈ R_k]\; _k\;p_k(x_S) on a variable partition S=S1⊔S2X_S=X_S_1 _S_2, where each expert factors as pk(S)=pk(1)(S1)pk(2)(S2)p_k(x_S)=p_k^(1)(x_S_1)\,p_k^(2)(x_S_2). If Rk\R_k\ is aligned with the partition, then the partition function decomposes as ∫f(S)S=∑kπk(∫Rk(1)pk(1)(S1)S1)(∫Rk(2)pk(2)(S2)S2). f(x_S)\,dx_S= _k _k ( _R_k^(1)p_k^(1)(x_S_1)\,dx_S_1 ) ( _R_k^(2)p_k^(2)(x_S_2)\,dx_S_2 ). A simple way to satisfy Definition 3.4 is to use axis-aligned partitions, where each RkR_k is a Cartesian product of intervals (or, more generally, a product of lower-dimensional sets). This includes rectangular boxes, and recovers the tractable region decompositions used implicitly by cutset-style splits (Rahman et al., 2014). However, axis-aligned regions are substantially less expressive than general convex polytopes and often require many rectangles to approximate a single slanted facet. This motivates two complementary directions. We can (1) design geometry-aware gating mechanisms that align with the circuit decomposition so that the induced constraints factor compatibly with the circuit structure and exact tractability is recovered. Alternatively, we can (2) accept intractability to obtain additional expressiveness and derive certified lower and upper bounds on the partition function, which enables an approximate inference with guarantees. We now present both approaches. 3.1 Certified Approximate Inference Though geometric gating breaks exact tractability, it is possible retain the reliability of PCs via inference procedures that produce certificates, i.e. provable lower and upper bounds on partition functions, marginals, and conditionals. In this section, we develop a general certified inference framework for Voronoi-Tessellated PCs. The main idea is to replace intractable polyhedral regions with tractable axis-aligned regions for which integration is compatible with decomposability, and to propagate the resulting local bounds through the circuit. Bounding Cell-Restricted Integrals with Boxes. Let Vk⊂ℝdV_k ^d be a Voronoi cell and Ω⊆ℝd ^d a bounded domain. An inner box Bk−B_k^- and outer box Bk+B_k^+ satisfy Bk−⊆Vk∩Ω⊆Bk+,B_k^- V_k∩ B_k^+, where both are axis-aligned: Bk±=∏i=1d[ai±,bi±]B_k^±= _i=1^d[a_i^±,b_i^±]. In practice, the domain Ω may be the full space or a data-dependent bounding box (e.g., per-variable min/max or high-probability truncation). The nesting property yields immediate bounds for any non-negative integrand. Lemma 3.6 (Cell Integral Bounds). Let p:ℝd→ℝ≥0p:R^d _≥ 0 be a non-negative density and let Bk−⊆Vk⊆Bk+B_k^- V_k B_k^+. Then ∫Bk−p()≤∫Vkp()≤∫Bk+p() _B_k^-p(x)\,dx≤ _V_kp(x)\,dx≤ _B_k^+p(x)\,dx Thus, the hard part becomes: (i) constructing useful boxes Bk±B_k^±, and (i) integrating circuit outputs over axis-aligned boxes efficiently. We now describe how we construct box approximations for Voronoi cells, optionally restricted to a bounded domain Ω=∏i[ℓi,ui] = _i[ _i,u_i]. Proposition 3.7 (Outer Box Computation). Let Pk=Vk∩Ω⊆ℝdP_k=V_k∩ ^d where Ω=∏i=1d[ℓi,ui] = _i=1^d[ _i,u_i] is a box domain. The tightest axis-aligned outer box containing PkP_k is given by Bk+=∏i=1d[min∈Pkxi,max∈Pkxi].B_k^+= _i=1^d [ _x∈ P_kx_i,\ _x∈ P_kx_i ]. Each bound is thus the optimum of a linear program over the polytope PkP_k. The LP constraints follow directly from the half-space representation of Voronoi cells together with the box constraints from Ω . We can construct a valid inner box Bk−B_k^- by centering it at the centroid kc_k and choosing per-coordinate radii so that the box remains inside all Voronoi half-space constraints. We provide a closed-form conservative construction below. Proposition 3.8 (Inner Box Construction). Assume Ω=ℝd =R^d for simplicity. Let δk:=minj≠k‖k−j‖2 _k:= _j≠ k\|c_k-c_j\|_2 be the nearest-centroid distance. Then the axis-aligned box Bk−=∏i=1d[k,i−r,k,i+r],B_k^-= _i=1^d[c_k,i-r,\ c_k,i+r], with r=δk2d,r= _k2 d, satisfies Bk−⊆VkB_k^- V_k. When Ω is bounded, we can simply intersect Bk−B_k^- with Ω to preserve containment within Vk∩ΩV_k∩ :Bk−←Bk−∩ΩB_k^-← B_k^-∩ . This construction trades tightness for simplicity and robustness. If tighter inner boxes are needed we can optimize radii rir_i per dimension subject to the Voronoi half-space constraints. 3.1.1 Anytime Bound Refinement Algorithm 1 Adaptive Anytime Bound Refinement 0: Voronoi-PC C, target gap ϵε, max iters T, dom. Ω 0: Bounds (Z−,Z+)(Z^-,Z^+) with Z+−Z−≤ϵZ^+-Z^-≤ε when achievable 1: Initialize partition ←ΩP←\ \ and classify all boxes for each cell 2: (Z−,Z+)←CertifiedBounds(,)(Z^-,Z^+)← CertifiedBounds(C,P) 3: for t=1t=1 to T do 4: if Z+−Z−≤ϵZ^+-Z^-≤ε then 5: return (Z−,Z+)(Z^-,Z^+) 6: end if 7: // Select boundary box with largest gap contribution 8: (n∗,k∗,B∗)←argmax(n^*,k^*,B^*)← gap contribution over Boundary boxes 9: // where wn,kw_n,k is weight in global bound computation 10: j∗←argmaxj(bj−aj)j^*← _j(b_j-a_j) for B∗=∏j[aj,bj]B^*= _j[a_j,b_j] // Longest dimension 11: Bisect B∗B^* along dimension j∗j^* into BL,BRB_L,B_R 12: Remove B∗B^* from P; add BL,BRB_L,B_R to P 13: Reclassify BL,BRB_L,B_R for all cells using memb. tests 14: (Z−,Z+)←CertifiedBounds(,)(Z^-,Z^+)← CertifiedBounds(C,P) 15: end for 16: return (Z−,Z+)(Z^-,Z^+) The basic box approximations above provide valid certified bounds but may be loose, particularly in high dimensions. We thus develop an anytime refinement algorithm that monotonically tightens the bounds through recursive box subdivision, enabling flexible trade-offs between computational cost and bound quality. The core idea is to recursively bisect boxes and test sub-boxes for containment or intersection with Voronoi cells. For outer boxes, we can discard sub-boxes that don’t intersect the cell, and for inner boxes, we retain sub-boxes fully contained within the cell. Formally, we maintain a disjoint axis-aligned partition P of domain ΩS=∏i∈S[ℓi,ui] _S= _i∈ S[ _i,u_i] where each box B∈B is labeled for each Voronoi cell VkV_k as: Inside if B⊆VkB V_k, Outside if B∩Vk=∅B∩ V_k= , or Boundary otherwise. This induces approximations Vk−():=⋃B:labk(B)=InsideBV_k^-(P):= _B:\,lab_k(B)= InsideB and Vk+():=⋃B:labk(B)≠OutsideBV_k^+(P):= _B:\,lab_k(B)≠ OutsideB satisfying Vk−()⊆Vk∩ΩS⊆Vk+()V_k^-(P) V_k∩ _S V_k^+(P). Since P is disjoint, integration decomposes additively with node bounds ∑kπkIk±() _k _kI_k^±(P) propagating via Theorem 3.10. Classification can be done using the half-space representation Vk=⋂j≠k:(j−k)⊤≤12(‖j‖2−‖k‖2)V_k= _j≠ k\x:(c_j-c_k) x≤ 12(\|c_j\|^2-\|c_k\|^2)\. For box B=∏i[li,ui]B= _i[l_i,u_i] and half-space normal a, the extrema max∈Ba⊤ _x∈ Ba x and min∈Ba⊤ _x∈ Ba x are computed by selecting uiu_i (resp. lil_i) when ai≥0a_i≥ 0 for the maximum, and vice versa for the minimum. Then B⊆VkB V_k iff all half-space maxima satisfy the constraint, and B∩Vk=∅B∩ V_k= if any half-space minimum violates its constraint. Refinement bisects Boundary boxes at midpoints, replaces them with disjoint children, reclassifies, and recomputes bounds, prioritizing boxes with largest gap contribution. Algorithm 1 outlines the refinement, and the below theorem establishes its monotone tightening and convergence properties. Theorem 3.9. Let tP_t denote the partition after t refinement steps with bounds (Zt−,Zt+)(Z_t^-,Z_t^+). Then (i) Zt−≤Z≤Zt+Z_t^-≤ Z≤ Z_t^+ ∀t∀\ t (i) Zt−≤Zt+1−Z_t^-≤ Z_t+1^- and Zt+1+≤Zt+Z_t+1^+≤ Z_t^+ for all t (i) if refinement drives the boundary volume μ(Vk+(t)∖Vk−(t))→0μ(V_k^+(P_t) V_k^-(P_t))→ 0 for each cell, then limt→∞Zt±=Z _t→∞Z_t^±=Z. Under uniform refinement, the gap scales as Zt+−Zt−=O(2−t/d)(Z0+−Z0−)Z_t^+-Z_t^-=O(2^-t/d)(Z_0^+-Z_0^-), requiring depth O(dlog(1/ϵ))O(d (1/ε)) to achieve target gap ϵε. 3.1.2 Propagating Bounds Through the Circuit. We next show how to propagate local cell bounds through the sum and product structure of the circuit. We first focus on the partition function (Z)(Z), and later discuss how the same machinery applies to marginals and conditionals. Theorem 3.10 (Bound Propagation). Let C be a Voronoi-gated PC. For each node n let In=∫fn()I_n= f_n(x)\,dx denote the integral of the subcircuit rooted at n, and let (In−,In+)(I_n^-,I_n^+) denote certified bounds on InI_n. We can compute bounds bottom up as follows. If n is a leaf, then In−=In+=∫pn(x)xI_n^-=I_n^+= p_n(x)\,dx. If n is a product node with children cj\c_j\ and disjoint scopes, then In−=∏jIcj−,In+=∏jIcj+.I_n^-= _jI_c_j^-, I_n^+= _jI_c_j^+. If n is a standard sum node with children cj\c_j\ and weights πn,cj\ _n,c_j\, then In−=∑jπn,cjIcj−,In+=∑jπn,cjIcj+.I_n^-= _j _n,c_jI_c_j^-, I_n^+= _j _n,c_jI_c_j^+. If n is a Voronoi-gated sum node with children pk\p_k\ and cell boxes (Bk−,Bk+)\(B_k^-,B_k^+)\, then In−=∑kπk∫Bk−pk();In+=∑kπk∫Bk+pk()I_n^-= _k _k _B_k^-p_k(x)\,dx;\ I_n^+= _k _k _B_k^+p_k(x)\,dx. At the root, we obtain bounds (Z−,Z+)(Z^-,Z^+) satisfying Z−≤Z≤Z+Z^-≤ Z≤ Z^+. Algorithm 2 Certified Bound Computation 0: Voronoi-gated PC C, box approximations (Bk−,Bk+)\(B_k^-,B_k^+)\ for each Voronoi cell 0: Bounds (Z−,Z+)(Z^-,Z^+) on Z=∫f()Z= f_C(x)\,dx 1: for each node n in reverse topological order do 2: if n is a leaf then 3: In−←∫pnI_n^-← p_n; In+←∫pnI_n^+← p_n 4: else if n is a product node then 5: In−←∏c∈ch(n)Ic−I_n^-← _c (n)I_c^-; In+←∏c∈ch(n)Ic+I_n^+← _c (n)I_c^+ 6: else if n is a standard sum node then 7: In−←∑cπn,cIc−I_n^-← _c _n,cI_c^-; In+←∑cπn,cIc+I_n^+← _c _n,cI_c^+ 8: else if n is a Voronoi-gated sum node then 9: for each cell k do 10: Jk−←IntegrateBox(pk,Bk−)J_k^-← IntegrateBox(p_k,B_k^-) 11: Jk+←IntegrateBox(pk,Bk+)J_k^+← IntegrateBox(p_k,B_k^+) 12: end for 13: In−←∑kπkJk−I_n^-← _k _kJ_k^-; In+←∑kπkJk+I_n^+← _k _kJ_k^+ 14: end if 15: end for 16: return (Iroot−,Iroot+)(I_root^-,I_root^+) Algorithm 2 summarizes the bound computation. The only non-standard operation is integrating a decomposable subcircuit over a box, which as we show in Theorem3.14 can be implemented recursively by applying interval integration at leaves and factorization at product nodes. We can apply the same machinery to bound marginals and conditionals. Corollary 3.11 (Marginal Bounds). For p(A)=∫p()A¯p(x_A)= p(x)\,dx_ A we obtain bounds (p−(A),p+(A))(p^-(x_A),p^+(x_A)) by propagating box restricted bounds through the circuit after restricting each box approximation to the A¯ A dimensions. Corollary 3.12 (Conditional Bounds). For disjoint A,BA,B and p(A∣B)=p(A,B)/p(B)p(x_A _B)=p(x_A,x_B)/p(x_B), if we have bounds p−(A,B)≤p(A,B)≤p+(A,B)p^-(x_A,x_B)≤ p(x_A,x_B)≤ p^+(x_A,x_B) and p−(B)≤p(B)≤p+(B)p^-(x_B)≤ p(x_B)≤ p^+(x_B) with p−(B)>0p^-(x_B)>0, then p−(A,B)p+(B)≤p(A∣B)≤p+(A,B)p−(B). p^-(x_A,x_B)p^+(x_B)≤ p(x_A _B)≤ p^+(x_A,x_B)p^-(x_B). Certified bounds provide guaranteed inference in settings where exact computation is intractable, but they have clear limitations. In high dimensions the initial inner boxes may capture only a small fraction of each Voronoi cell, leading to loose initial bounds. Tight bounds require many refinements, and the refinement cost grows quickly with dimension and the number of cells. Finally, while refinement can approximate exact inference arbitrarily well, it does not restore exactness at finite computation. This motivates designing Voronoi gating mechanisms that align with circuit decomposition and recover exact tractable inference. 3.2 Hierarchical Factorized Voronoi PCs Next, we show how to enforce geometric alignment by design through hierarchical factorization of Voronoi tessellations to retain tractable inference. The key idea is to partition the Voronoi centroids in a manner that mirrors the circuit’s vtree structure, ensuring that at each product node, the induced Voronoi cells decompose into independent factors over disjoint variable subsets. Let a scope SX_S be partitioned into disjoint blocks S=⨆i=1mSiX_S= _i=1^mX_S_i. For each block i, choose centroids 1(i),…,Ki(i)⊂ℝ|Si|\c^(i)_1,…,c^(i)_K_i\ ^|X_S_i| and let Vki(i)ki=1Ki\V^(i)_k_i\_k_i=1^K_i be the induced Voronoi cells in ℝ|Si|R^|X_S_i|. These induce a joint partition of ℝ|S|R^|X_S| into product cells indexed by =(k1,…,km)k=(k_1,…,k_m) as: V=Vk1(1)×⋯×Vkm(m).V_k\;=\;V^(1)_k_1×·s× V^(m)_k_m. The corresponding hard gate factors along blocks: g(S)=[S∈V]=∏i=1m[Si∈Vki(i)]=∏i=1mgki(i)(Si),g_k(x_S)=I[x_S∈ V_k]= _i=1^mI[x_S_i∈ V^(i)_k_i]= _i=1^mg^(i)_k_i(x_S_i), so membership can be decided independently within each block. The defining feature is that each joint cell is a Cartesian product of lower-dimensional Voronoi cells, which is exactly the geometric analogue of decomposability. We can now define a gated sum node that uses a factorized Voronoi partition and couples it to a matching factorized expert. Definition 3.13 (HFV-Gated Sum Node). Let S=⨆i=1mSiX_S= _i=1^mX_S_i. An HFV-gated sum node over scope SX_S consists of a factorized Voronoi partition with KiK_i cells on factor i, mixture weights π∈[K1]×⋯×[Km]\ _k\_k∈[K_1]×·s×[K_m] with ∑π=1 _k _k=1, and factor subcircuits pki(i)\p^(i)_k_i\ where pki(i)p^(i)_k_i has scope SiX_S_i. The node computes f(S)=∑g(S)π∏i=1mpki(i)(Si),f(x_S)= _kg_k(x_S)\, _k\, _i=1^mp^(i)_k_i(x_S_i), where g(S)=∏i=1mgki(i)(Si)g_k(x_S)= _i=1^mg^(i)_k_i(x_S_i). The critical design choice is that the gate and the expert share the same factorization pattern. This alignment restores the ability to apply Fubini’s theorem 111Proofs are given in the appendix. to reduce high-dimensional integrals into products of lower-dimensional integrals. To obtain a full circuit over =X1,…,XDX=\X_1,…,X_D\, we align HFV gating with a variable tree (vtree), so that at each internal vtree node the gate factors across its left/right child scopes. This yields a hierarchical (multi-resolution) geometric partition: coarse routing happens at higher scopes, while finer routing refines decisions within smaller scopes. We now show that HFV restores exact tractable inference. The result mirrors standard PC tractability, but incurs an additional multiplicative factor that reflects the number of joint Voronoi cell combinations at each gated sum. Theorem 3.14 (Tractability of HFV-PCs). Let C be an HFV-PC with |||C| nodes, maximum factorization degree m, and at most K Voronoi cells per factor. Then the partition function, marginals, and conditionals are computable exactly in time O(||Km)O(|C|K^m). For binary HFV-PCs the time is O(||K2)O(|C|K^2). HFV intentionally restricts geometry to preserve decomposability: it can represent rich piecewise structure within each scope block and refine it hierarchically, but it cannot realize arbitrary oblique polytopes spanning variables across different blocks at a given node. When such geometry is essential, we can revert to general Voronoi gating and still preserve reliability via certified approximate inference. Figure 1: Mean Test Log-likelihood (↑ ) on synthetic 2D and 3D density estimation tasks achieved by EinsumNet and HCLT along with their geometry-aware extensions using Voronoi tessellations (VT) and hierarchical factorized Voronoi (HFV), averaged across 33 trials. For VT, values correspond to the lower bound on the log-likelihood obtained via our certified approximate inference framework. 3.3 Learning via Soft Gating Both geometry-aware constructions we have introduced in this paper rely on hard routing, where inputs activate a single region. This yields crisp locality and interpretability and enables exact inference in HFV-PCs as well as certified bounds for general VT-PCs. However, hard Voronoi assignments are non-differentiable, preventing gradient-based learning of the centroids together with the PC parameters. We thus introduce a smooth relaxation that supports standard backpropagation during training, and revert to hard gating at test time to recover the desired inference guarantees. Definition 3.15 (Soft Voronoi Gate). Given centroids 1,…,K⊂ℝd\c_1,…,c_K\ ^d and an inverse temperature α>0α>0, we define the soft Voronoi gate wk(;α)=exp(−α‖−k‖2)∑j=1Kexp(−α‖−j‖2).w_k(u;α)= \! (-α\|u-c_k\|^2 ) _j=1^K \! (-α\|u-c_j\|^2 ). (2) For any u, the weights satisfy wk(;α)>0w_k(u;α)>0 and ∑k=1Kwk(;α)=1 _k=1^Kw_k(u;α)=1 and are smooth in both u and k\c_k\ The temperature parameter controls the sharpness of routing. When α is small the weights are diffuse, and when α is large the weights concentrate on the nearest centroid. Correspondingly, a soft Voronoi-gated sum node computes: f(S;α)=∑k=1Kwk(S;α)πkpk(S).f(x_S;α)= _k=1^Kw_k(x_S;α) _kp_k(x_S). For HFV-PCs, factorization can be maintained by applying soft gates per factor. Given partition S=⨆i=1mSiX_S= _i=1^mX_S_i, define a factorized soft gate w(S;α)=∏i=1mwki(i)(Si;α)w_k(x_S;α)= _i=1^mw_k_i^(i)(x_S_i;α) where each factor uses centroids in ℝ|Si|R^|X_S_i|. This ensures gradient signals to factor-i centroids depend only on variables in SiX_S_i, preserving the decomposition structure during optimization. For any α>0α>0, the soft gates are smooth and the circuit likelihood is differentiable with respect to all parameters. The centroid gradient has the form ∇kwk(;α)=2αwk(;α)(1−wk(;α))(−k) _c_kw_k(u;α)=2α w_k(u;α)(1-w_k(u;α))(u-c_k) The factor wk(1−wk)w_k(1-w_k) peaks when the gate is uncertain (near decision boundaries), and vanishes when routing is already confident (wk≈0w_k≈ 0 or 11). Thus centroids are primarily updated in regions where the current tessellation is contested. Soft-to-Hard Convergence. We train VT-PCs and HFV-PCs via soft gates but ultimately require the exact inference guarantees of hard HFV gating. The next result highlights that increasing α recovers hard Voronoi assignments, with an exponential rate governed by a geometric margin. Theorem 3.16 (Soft-to-Hard Convergence). Let gk()=[∈Vk]g_k(u)=I[u∈ V_k] be the hard Voronoi gate induced by k\c_k\ and let wk(;α)w_k(u;α) be the soft gate in (2). Then for any u not lying on a Voronoi boundary, limα→∞wk(;α)=gk() _α→∞w_k(u;α)=g_k(u). Moreover, if k∗()=argmink‖−k‖k^*(u)= _k\|u-c_k\| and the margin γ()=minj≠k∗()(‖−j‖2−‖−k∗‖2)γ(u)= _j≠ k^*(u) (\|u-c_j\|^2-\|u-c_k^*\|^2 ) is positive, then 1−wk∗(;α)≤(K−1)e−αγ()1-w_k^*(u;α)≤(K-1)e^-α\,γ(u). Further, if p is integrable, then wk(⋅;α)→gk(⋅)w_k(·;α)→ g_k(·) also implies limα→∞∫wk(;α)p()=∫Vkp() _α→∞ w_k(u;α)\,p(u)\,du= _V_kp(u)\,du. The gap γ()γ(u) quantifies how much closer u is to its nearest centroid than to the runner-up. Larger margins mean the softmax ratios exp(−α(dj−dk∗)) (-α(d_j-d_k^*)) decay faster, so routing becomes effectively hard at smaller α. The exponential bound formalizes this geometric picture and justifies annealing: as training progresses, increasing α sharpens routing while remaining stable on points that are already well-separated by the current centroids. Training. We thus train soft HFV-PCs and VT-PCs via maximum-likelihood, gradually increasing α so that routing sharpens over time. This annealing schedule avoids early training instabilities where centroids collapse or assignments become overly brittle before the experts have adapted. We use softmax projection to enforce π≥0 _k≥ 0 and ∑π=1 _k _k=1. After training we use the learned centroids to define a Voronoi tessellation and perform inference in one of two modes. Hard-gated inference, which replaces wki(i)w^(i)_k_i with gki(i)g^(i)_k_i and recovers the exact tractability guarantees. This is the default mode we use in experiments because it preserves exact partition functions and exact marginals. Soft-gated inference uses a finite α. This yields a fully smooth model, but exact marginalization generally requires integrating smooth gate functions, which typically forfeit the exact HFV guarantees. We therefore treat soft gating primarily as a training device and hard gating as the inference-time model. 4 Experiments & Results To validate our theoretical framework, we consider eight synthetic distributions where geometric structure is explicit and exact verification is feasible: four 2D2D (Alphabet, CheckerBoard, Pinwheel, Spiral) and four 3D3D (BentLissajous, InterlockedCircles, Knotted, TwistedEight) (Sidheekh et al., 2022, 2023), each with 1010k train, 55k validation, and 55k test samples. We compare two notable base PC architectures: EinsumNet (random binary region graph) (Peharz et al., 2020a) and HCLTs (Chow-Liu trees) (Liu and Van den Broeck, 2021) against their geometry-aware variants: VT-EinsumNet/HCLT using Voronoi tessellations at root sum nodes with certified approximate inference, and HFV-EinsumNet/HCLT using hierarchical factorized Voronoi gating aligned with circuit decomposition for exact inference. We use Gaussian leaves and Tucker sum-product layers in all models (Loconte et al., 2025a), keeping the base architecture same, with 1010 input and sum units for 3D3D (and 55 each for 2D2D) datasets, for a fair comparison. We employ maximum likelihood via stochastic gradient descent to learn the parameters of the PC and VT, by training using an Adam optimizer with a leatning rate of 0.010.01 and batch size of 500500 for 100100 epochs. For VT and HFV variants, we employ soft gating with linear temperature annealing for and k-means centroid initialization (100100 iterations), and switch to hard gating at test time. VT models report certified lower bounds using inner box approximations without refinement. Figure 2: Visualization of the distribution and voronoi tessellations learned by a VT-EinsumNet (left) and HFV-EinsumNet (right) on the 2D2D pinwheel dataset. The axis aligned boxes in the left figure represent the Inner Boxes computed for estimating the lower bound on the partition function using our conservative construction. Results. Figure 1 depicts the mean test log-likelihoods across all datasets and model variants. VT models (triangles) consistently achieve strong performance, with their certified lower bounds often exceeding baseline exact log-likelihoods, demonstrating that the increased expressivity from unconstrained geometry-aware routing captures structure missed by input-independent weights, even when accounting for conservative approximation gaps. HFV models (squares) achieve performance comparable to baselines (circles). This is expected in our low-dimensional, shallow-circuit regime, as the alignment constraints required for exact tractability in HFV induce fully deterministic, factorized partitions, which can reduce expressive power relative to smooth decomposable PCs. The advantage of HFV here however, is conceptual and algorithmic: it retains exact tractable inference while providing an explicit geometric interpretation that can be useful for downstream tasks such as continual learning. Figure 2 visualizes the learned routing structure on the 2D pinwheel dataset. In VT (left), the Voronoi cells adapt to the arms of the distribution and assign regions of responsibility to local experts; the axis-aligned inner boxes shown are the conservative subsets used by our certified lower-bound computation, and we see that in practice they capture the majority of the modeled probability mass within each cell on this dataset. In HFV (right), the partition is hierarchical and factorized, yielding axis-aligned regions that preserve exact tractability by construction. Together, these overlays offer an interpretable view of where specialization occurs and help explain why geometry-aware routing is effective on distributions with strong locality. Figure 3 shows validation learning curves of EinsumNet and HCLT on the 2D spiral dataset. For VT we plot the certified lower bound (solid red) together with the certified interval induced by the partition function Z∈[Z−,Z+]Z∈[Z^-,Z^+] (shaded region). As training proceeds and the temperature annealing sharpens the gates, the lower bound increases steadily and the envelope stabilizes reflecting improved fit of local experts and tighter normalization bounds as the learned partitions align with the data support. HFV remains tractable throughout, and its curve corresponds to exact likelihood under the aligned factorized gating. Overall, these results suggest that VT offers substantial gains while remaining certifiable, and that soft-gate training with temperature annealing yields stable learning dynamics. Figure 3: Learning curves on the 2D spiral dataset. Validation log-likelihood across epochs for EinsumNet (left) and HCLT (right), averaged over 33 trials. VT models report a certified log-likelihood lower bound while baselines and HFV use exact tractable evaluation. 5 Conclusion We developed a principled foundation for introducing geometric awareness into PCs, by replacing constant sum-node weights with geometry driven gates induced via Voronoi tessellations. Our key message is that retaining tractable inference requires geometric alignment with the circuits factorization, as general Voronoi partitioning breaks recursive marginalization. We addressed this tension in two ways: first, we developed a theoretically grounded certified approximate inference framework that preserved geometric expressivity while providing reliable inference. Second, we identified structural properties on the tessellation under which tractable inference was recovered. Finally, we presented a soft relaxation with convergence guarantees that enabled stable gradient based learning, and empirically validated the resulting model on low dimensional geometric manifolds. Future work includes extending the framework to learned embeddings, developing tighter certification schemes for scaling to higher dimensional data, and leveraging geometric awareness for continual learning, controlled generation and interpretable anomaly detection using PCs. Acknowledgements The authors gratefully acknowledge the generous support by the AFOSR award FA9550-23-1-0239, the ARO award W911NF2010224 and the DARPA Assured Neuro Symbolic Learning and Reasoning (ANSR) award HR001122S0039. References F. Aurenhammer (1991) Voronoi diagrams—a survey of a fundamental geometric data structure. ACM Computing Surveys (CSUR) 23 (3), p. 345–405. Cited by: §1, §3. S. Braun, S. Sidheekh, A. Vergari, M. Mundt, S. Natarajan, and K. Kersting (2025) Tractable representation learning with probabilistic circuits. Transactions on Machine Learning Research. External Links: ISSN 2835-8856 Cited by: §1. R. T. Chen, B. Amos, and M. Nickel (2022) Semi-discrete normalizing flows through differentiable tessellation. Advances in Neural Information Processing Systems 35. Cited by: §1, §3. Y. Choi, A. Vergari, and G. Van den Broeck (2020) Probabilistic circuits: a unifying framework for tractable probabilistic models. UCLA, Tech. Rep. Cited by: §1, §2. A. H. Correia, G. Gala, E. Quaeghebeur, C. De Campos, and R. Peharz (2023) Continuous mixtures of tractable probabilistic models. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 37, p. 7244–7252. Cited by: §2. M. Dang, A. Liu, and G. Van den Broeck (2022) Sparse probabilistic circuits via pruning and growing. Advances in Neural Information Processing Systems 35, p. 28374–28385. Cited by: §2. A. Darwiche (2003) A differential approach to inference in bayesian networks. Journal of the ACM (JACM) 50 (3), p. 280–305. Cited by: §1, §2, §2. P. Z. Dos Martires (2024) Probabilistic neural circuits. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 38, p. 17280–17289. Cited by: §2. Q. Du, V. Faber, and M. Gunzburger (1999) Centroidal voronoi tessellations: applications and algorithms. SIAM review 41 (4), p. 637–676. Cited by: §3. M. E. Dyer and A. M. Frieze (1988) On the complexity of computing the volume of a polyhedron. SIAM Journal on Computing 17 (5), p. 967–974. Cited by: §1. G. Gala, C. de Campos, R. Peharz, A. Vergari, and E. Quaeghebeur (2024) Probabilistic integral circuits. In International Conference on Artificial Intelligence and Statistics, p. 2143–2151. Cited by: §2. R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton (1991) Adaptive mixtures of local experts. Neural computation 3 (1), p. 79–87. Cited by: §2. A. Karanam, S. Mathur, S. Sidheekh, and S. Natarajan (2025) A unified framework for human-allied learning of probabilistic circuits. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39, p. 17779–17787. Cited by: §2. D. Kisa, G. Van den Broeck, A. Choi, and A. Darwiche (2014) Probabilistic sentential decision diagrams. In Fourteenth International Conference on the Principles of Knowledge Representation and Reasoning, Cited by: §2. T. A. Lab (2024) cirkit External Links: Link Cited by: Appendix C. A. Liu, K. Ahmed, and G. Van den Broeck (2024) Scaling tractable probabilistic circuits: a systems perspective. In Proceedings of the 41th International Conference on Machine Learning (ICML), Cited by: §2. A. Liu and G. Van den Broeck (2021) Tractable regularization of probabilistic circuits. Advances in Neural Information Processing Systems 34, p. 3558–3570. Cited by: §2, §4. A. Liu, H. Zhang, and G. V. den Broeck (2023a) Scaling up probabilistic circuits by latent variable distillation. In The Eleventh International Conference on Learning Representations, Cited by: §2. X. Liu, A. Liu, G. Van den Broeck, and Y. Liang (2023b) Understanding the distillation process from deep generative models to tractable probabilistic circuits. In International Conference on Machine Learning, p. 21825–21838. Cited by: §2. L. Loconte, A. Mari, G. Gala, R. Peharz, C. de Campos, E. Quaeghebeur, G. Vessio, and A. Vergari (2025a) What is the relationship between tensor factorizations and circuits (and how can we exploit it)?. Transactions on Machine Learning Research. External Links: ISSN 2835-8856 Cited by: §2, §4. L. Loconte, S. Mengel, and A. Vergari (2025b) Sum of squares circuits. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39, p. 19077–19085. Cited by: §2. L. Loconte, A. M. Sladek, S. Mengel, M. Trapp, A. Solin, N. Gillis, and A. Vergari (2024) Subtractive mixture models via squaring: representation and learning. In The Twelfth International Conference on Learning Representations, Cited by: §2. G. L. Marchetti, V. Polianskii, A. Varava, F. T. Pokorny, and D. Kragic (2023) An efficient and continuous voronoi density estimator. In International Conference on Artificial Intelligence and Statistics, p. 4732–4744. Cited by: §3. R. Peharz, S. Lang, A. Vergari, K. Stelzner, A. Molina, M. Trapp, G. Van den Broeck, K. Kersting, and Z. Ghahramani (2020a) Einsum networks: fast and scalable learning of tractable probabilistic circuits. In International Conference on Machine Learning, p. 7563–7574. Cited by: §2, §4. R. Peharz, A. Vergari, K. Stelzner, A. Molina, X. Shao, M. Trapp, K. Kersting, and Z. Ghahramani (2020b) Random sum-product networks: a simple and effective approach to probabilistic deep learning. In UAI, Cited by: §2. V. Polianskii, G. L. Marchetti, A. Kravberg, A. Varava, F. T. Pokorny, and D. Kragic (2022) Voronoi density estimator for high-dimensional data: computation, compactification and convergence. In Uncertainty in Artificial Intelligence, p. 1644–1653. Cited by: §3. H. Poon and P. Domingos (2011) Sum-product networks: a new deep architecture. In 2011 IEEE International Conference on Computer Vision Workshops (ICCV Workshops), p. 689–690. Cited by: §1, §2, §2. T. Rahman, P. Kothalkar, and V. Gogate (2014) Cutset networks: a simple, tractable, and scalable approach for improving the accuracy of chow-liu trees. In ECML PKDD, Cited by: §2, §2, §3. X. Shao, A. Molina, A. Vergari, K. Stelzner, R. Peharz, T. Liebig, and K. Kersting (2022) Conditional sum-product networks: modular probabilistic circuits via gate functions. International Journal of Approximate Reasoning. Cited by: §2. O. Sharir and A. Shashua (2018) Sum-product-quotient networks. In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, A. Storkey and F. Perez-Cruz (Eds.), Proceedings of Machine Learning Research, Vol. 84, p. 529–537. Cited by: §2. N. Shazeer, A. Mirhoseini, K. Maziarz, A. Davis, Q. Le, G. Hinton, and J. Dean (2017) Outrageously large neural networks: the sparsely-gated mixture-of-experts layer. In International Conference on Learning Representations, Cited by: §2. A. Shih, D. Sadigh, and S. Ermon (2021) Hyperspns: compact and expressive probabilistic circuits. Advances in neural information processing systems 34, p. 8571–8582. Cited by: §2. S. Sidheekh, C. B. Dock, T. Jain, R. Balan, and M. K. Singh (2022) VQ-flows: vector quantized local normalizing flows. In Uncertainty in Artificial Intelligence, p. 1835–1845. Cited by: §B.2, §4. S. Sidheekh, K. Kersting, and S. Natarajan (2023) Probabilistic flow circuits: towards unified deep models for tractable probabilistic inference. In Uncertainty in Artificial Intelligence, p. 1964–1973. Cited by: §B.2, §2, §4. S. Sidheekh and S. Natarajan (2024) Building expressive and tractable probabilistic generative models: A review. In Proceedings of the Thirty-Third International Joint Conference on Artificial Intelligence, 2024, p. 8234–8243. Cited by: §2. S. Sidheekh, P. Tenali, S. Mathur, E. Blasch, K. Kersting, and S. Natarajan (2025) Credibility-aware multimodal fusion using probabilistic circuits. In The 28th International Conference on Artificial Intelligence and Statistics, Cited by: §1. S. Sidheekh, P. Tenali, S. Mathur, E. Blasch, and S. Natarajan (2024) On the robustness and reliability of late multi-modal fusion using probabilistic circuits. In 2024 27th International Conference on Information Fusion (FUSION), p. 1–8. Cited by: §1. H. Suresh, S. Sidheekh, V. Shreeram M.P, S. Natarajan, and N. C. Krishnan (2026) Tractable sharpness-aware learning of probabilistic circuits. In Proceedings of the 40th AAAI Conference on Artificial Intelligence, Cited by: §2. J. Veness, T. Lattimore, D. Budden, A. Bhoopchand, C. Mattern, A. Grabska-Barwinska, S. Ozair, S. Gómez, A. Guez, D. Hassabis, et al. (2021) Gated linear networks. Proceedings of the AAAI Conference on Artificial Intelligence 35 (11), p. 10015–10023. Cited by: §2. A. Vergari, N. Di Mauro, and F. Esposito (2015) Simplifying, regularizing and strengthening sum-product network structure learning. In Joint European conference on machine learning and knowledge discovery in databases, p. 343–358. Cited by: §2. M. Zečević, D. Dhami, A. Karanam, S. Natarajan, and K. Kersting (2021) Interventional sum-product networks: causal inference with tractable probabilistic models. Advances in neural information processing systems 34, p. 15019–15031. Cited by: §1. H. Zhang, M. Dang, B. Wang, S. Ermon, N. Peng, and G. V. den Broeck (2025) Scaling probabilistic circuits via monarch matrices. In Forty-second International Conference on Machine Learning, Cited by: §2. Appendix A Proofs of Main Results A.1 Intractability of Voronoi-Gated PCs We begin by establishing why unconstrained Voronoi gating breaks tractable inference, even in the simplest case of a single sum node with fully factorized experts. Proposition A.1 (Single-Layer Intractability). Let f()=∑k=1Kgk()πkpk()f(x)= _k=1^Kg_k(x)\, _k\,p_k(x) be a Voronoi-gated sum node over variables =X1,…,XDX=\X_1,…,X_D\ where gk()=[∈Vk]g_k(x)=I[x∈ V_k] with Voronoi cells Vk\V_k\ induced by centroids 1,…,K⊂ℝD\c_1,…,c_K\ ^D. Suppose each expert is fully factorized: pk()=∏i=1Dpk(i)(xi)p_k(x)= _i=1^Dp_k^(i)(x_i). Then the partition function Z=∫f()=∑k=1Kπk∫Vkpk()Z= f(x)\,dx= _k=1^K _k _V_kp_k(x)\,dx requires integrating pkp_k over Voronoi cells VkV_k, which are convex polytopes with oblique boundaries. In general, the integral ∫Vk∏ipk(i)(xi)d _V_k _ip_k^(i)(x_i)\,dx does not factor into a product of one-dimensional integrals. Proof. Each Voronoi cell VkV_k is defined as the region of points closer to centroid kc_k than to any other centroid: Vk=∈ℝD:‖−k‖2≤‖−j‖2 for all j.V_k= \x ^D:\|x-c_k\|^2≤\|x-c_j\|^2 for all j \. Equivalently, VkV_k is the intersection of K−1K-1 half-spaces: Vk=⋂j≠kHkj−,V_k= _j≠ kH_kj^-, where Hkj−=:(j−k)⊤≤12(‖j‖2−‖k‖2)H_kj^-=\x:(c_j-c_k) x≤ 12(\|c_j\|^2-\|c_k\|^2)\. The boundary hyperplane between cells VkV_k and VjV_j is Hkj=:(j−k)⊤=12(‖j‖2−‖k‖2).H_kj= \x:(c_j-c_k) x= 12(\|c_j\|^2-\|c_k\|^2) \. For generic centroids (i.e., centroids not satisfying special symmetries), the normal vector j−kc_j-c_k has multiple nonzero components. This means the boundary hyperplane is not aligned with any coordinate axis, it is oblique with respect to the standard coordinate system. Consider the integral: ∫Vk∏i=1Dpk(i)(xi)d. _V_k _i=1^Dp_k^(i)(x_i)\,dx. Even though the integrand factors across variables, the region VkV_k does not. Specifically, VkV_k is not a Cartesian product of the form I1×⋯×IDI_1×·s× I_D where each Ii⊆ℝI_i is an interval. The oblique boundaries couple the coordinates: whether a point x lies in VkV_k depends on the joint configuration of all coordinates, not on independent conditions on each coordinate separately. To apply Fubini’s theorem and factor the integral as ∏i∫Iipk(i)(xi)xi _i _I_ip_k^(i)(x_i)\,dx_i, we would need the integration domain to be a product of univariate domains. Since VkV_k is not such a product, Fubini does not reduce the D-dimensional integral to a product of D one-dimensional integrals. Concrete Example (2D): Let D=2D=2 with centroids 1=(0,1)c_1=(0,1) and 2=(1,0)c_2=(1,0). The Voronoi boundary between V1V_1 and V2V_2 satisfies ‖(x1,x2)−(0,1)‖2=‖(x1,x2)−(1,0)‖2,\|(x_1,x_2)-(0,1)\|^2=\|(x_1,x_2)-(1,0)\|^2, which simplifies to x12+(x2−1)2=(x1−1)2+x22x_1^2+(x_2-1)^2=(x_1-1)^2+x_2^2, yielding x1=x2x_1=x_2. Thus V1=(x1,x2):x2>x1V_1=\(x_1,x_2):x_2>x_1\ (the region above the diagonal). Now consider ∫V1p1(1)(x1)p1(2)(x2)x1x2. _V_1p_1^(1)(x_1)p_1^(2)(x_2)\,dx_1\,dx_2. Integrating first over x2x_2 for fixed x1x_1: ∫V1p1(1)(x1)p1(2)(x2)x1x2=∫−∞p1(1)(x1)(∫x1∞p1(2)(x2)x2)x1. _V_1p_1^(1)(x_1)p_1^(2)(x_2)\,dx_1\,dx_2\\ = _-∞^∞p_1^(1)(x_1) ( _x_1^∞p_1^(2)(x_2)\,dx_2 )dx_1. The inner integral has lower limit x1x_1, which depends on the outer integration variable. This prevents factorization into (∫p1(1)(x1)x1)(∫p1(2)(x2)x2) ( p_1^(1)(x_1)dx_1 ) ( p_1^(2)(x_2)dx_2 ), and the dependence persists regardless of how simple the univariate densities p1(i)p_1^(i) are. Thus, the oblique geometry of Voronoi cells couples variables in a way that is fundamentally incompatible with the factorization exploited by decomposable probabilistic circuits. ∎ A.1.1 Deep Circuits and Constraint Accumulation Proposition A.2 (Constraint Accumulation). In a deep circuit of depth L with Voronoi gating at multiple levels, each contribution is integrated over a region determined by the active Voronoi cell at every ancestor sum node. The feasible set is an intersection of Voronoi regions restricted to different scopes. At the root, routing selects VkrootV_k_root over full scope X. At an internal node with scope S⊂S , routing selects VkSV_k_S over SX_S. After marginalizing out variables not in S, the root constraint projects onto SX_S in a way that generally doesn’t align with VkSV_k_S. The resulting region is an intersection of projected polytopes, which can be arbitrarily complex. Example A.3 (Two-Level Interaction). Consider vtree root scope X1,X2,X3,X4\X_1,X_2,X_3,X_4\ with children X1,X2\X_1,X_2\ and X3,X4\X_3,X_4\. • Root gate imposes x1+x2+x3+x4≤c1x_1+x_2+x_3+x_4≤ c_1. • Left child gate imposes x1−x2≤c2x_1-x_2≤ c_2. To integrate over X1,X2\X_1,X_2\, we need both the local constraint x1−x2≤c2x_1-x_2≤ c_2 and the projected root constraint. After marginalizing out X3,X4\X_3,X_4\, the root constraint becomes context-dependent, preventing independent marginalization. Thus, in deep circuits, unconstrained Voronoi gating introduces hierarchical geometric constraints that interact in complex ways after marginalization, compounding intractability. HFV-PCs avoid this by ensuring geometric constraints factor along circuit variable partitions at every level. A.2 Geometric Alignment and Tractability Recovery The preceding result shows that unconstrained Voronoi cells break tractability. We now prove that alignment with the circuit’s variable decomposition is sufficient to restore it. Theorem A.4 (Alignment Enables Factorization). Consider a gated sum node f(S)=∑k=1Kgk(S)πkpk(S)f(x_S)= _k=1^Kg_k(x_S)\, _k\,p_k(x_S) over scope S=S1⊔S2X_S=X_S_1 _S_2 (disjoint partition). Suppose each gating region RkR_k satisfies the alignment condition: there exist Rk(1)⊆ℝ|S1|R_k^(1) ^|X_S_1| and Rk(2)⊆ℝ|S2|R_k^(2) ^|X_S_2| such that Rk=Rk(1)×Rk(2),R_k=R_k^(1)× R_k^(2), and suppose each expert factors as pk(S)=pk(1)(S1)pk(2)(S2)p_k(x_S)=p_k^(1)(x_S_1)p_k^(2)(x_S_2). Then the partition function decomposes: ∫f(S)S=∑k=1Kπk(∫Rk(1)pk(1)(S1)S1)(∫Rk(2)pk(2)(S2)S2). f(x_S)\,dx_S= _k=1^K _k ( _R_k^(1)p_k^(1)(x_S_1)\,dx_S_1 ) ( _R_k^(2)p_k^(2)(x_S_2)\,dx_S_2 ). Proof. Start with the definition of f and substitute the gate and expert factorizations: ∫f(S)S f(x_S)\,dx_S =∫∑k=1K[S∈Rk]πkpk(S)dS = _k=1^KI[x_S∈ R_k]\, _k\,p_k(x_S)\,dx_S =∑k=1Kπk∫Rkpk(S)S. = _k=1^K _k _R_kp_k(x_S)\,dx_S. By the alignment condition, S∈Rkx_S∈ R_k if and only if S1∈Rk(1)x_S_1∈ R_k^(1) and S2∈Rk(2)x_S_2∈ R_k^(2) simultaneously. Since S1X_S_1 and S2X_S_2 are disjoint, we can write dS=dS1dS2dx_S=dx_S_1\,dx_S_2 and apply Fubini’s theorem: ∫Rkpk(S)S _R_kp_k(x_S)\,dx_S =∫Rk(1)×Rk(2)pk(1)(S1)pk(2)(S2)S1S2 = _R_k^(1)× R_k^(2)p_k^(1)(x_S_1)p_k^(2)(x_S_2)\,dx_S_1\,dx_S_2 =(∫Rk(1)pk(1)(S1)S1)(∫Rk(2)pk(2)(S2)S2). = ( _R_k^(1)p_k^(1)(x_S_1)\,dx_S_1 ) ( _R_k^(2)p_k^(2)(x_S_2)\,dx_S_2 ). Substituting back completes the proof. ∎ Remark A.5. This result generalizes immediately to m-way partitions S=⨆i=1mSiX_S= _i=1^mX_S_i with regions Rk=∏i=1mRk(i)R_k= _i=1^mR_k^(i) and experts pk(S)=∏ipk(i)(Si)p_k(x_S)= _ip_k^(i)(x_S_i). The key requirement is that both the gating regions and the expert distributions respect the same factorization structure. A.3 Certified Approximate Inference A.3.1 Inner Box Construction Proposition A.6 (Conservative Inner Box). Let Vk\V_k\ be Voronoi cells induced by centroids k⊂ℝd\c_k\ ^d. Define the nearest-centroid distance δk:=minj≠k‖k−j‖2. _k:= _j≠ k\|c_k-c_j\|_2. Then the axis-aligned box centered at kc_k with radius r=δk2dr= _k2 d satisfies Bk−:=∏i=1d[k,i−r,k,i+r]⊆VkB_k^-:= _i=1^d[c_k,i-r,c_k,i+r] V_k. Proof. Take any point ∈Bk−x∈ B_k^-. By construction, |xi−k,i|≤r|x_i-c_k,i|≤ r for all i, so ‖−k‖2=∑i=1d(xi−k,i)2≤d⋅r2=dr=δk2.\|x-c_k\|_2= _i=1^d(x_i-c_k,i)^2≤ d· r^2= d\,r= _k2. Now consider any other centroid jc_j with j≠kj≠ k. By the triangle inequality: ‖−j‖2≥‖k−j‖2−‖−k‖2≥δk−δk2=δk2.\|x-c_j\|_2≥\|c_k-c_j\|_2-\|x-c_k\|_2≥ _k- _k2= _k2. Thus ‖−k‖2≤δk2≤‖−j‖2\|x-c_k\|_2≤ _k2≤\|x-c_j\|_2 for all j≠kj≠ k, which means x is closer to kc_k than to any other centroid. By the definition of Voronoi cells, ∈Vkx∈ V_k. Since this holds for every ∈Bk−x∈ B_k^-, we have Bk−⊆VkB_k^- V_k. ∎ Remark A.7. The factor 1/(2d)1/(2 d) ensures containment but is conservative. A tighter construction can be obtained by optimizing the per-dimension radii rir_i subject to the half-space constraints defining VkV_k, but this requires solving a constrained optimization problem for each cell. The simple formula above trades tightness for computational convenience and robustness. Remark A.8 (Extension to bounded domains). If the domain is a box Ω=∏i[ℓi,ui] = _i[ _i,u_i], intersect: Bk−←Bk−∩ΩB_k^-← B_k^-∩ . This preserves Bk−⊆Vk∩ΩB_k^- V_k∩ because intersection with Ω only removes points outside the domain. A.3.2 Anytime Refinement: Convergence Analysis Theorem A.9 (Monotone Tightening and Convergence). Let tP_t denote the partition after t refinement steps with bounds (Zt−,Zt+)(Z_t^-,Z_t^+). Then (i) Zt−≤Z≤Zt+Z_t^-≤ Z≤ Z_t^+ ∀t∀\ t (i) Zt−≤Zt+1−Z_t^-≤ Z_t+1^- and Zt+1+≤Zt+Z_t+1^+≤ Z_t^+ for all t (i) if refinement drives the boundary volume μ(Vk+(t)∖Vk−(t))→0μ(V_k^+(P_t) V_k^-(P_t))→ 0 for each cell, then limt→∞Zt±=Z _t→∞Z_t^±=Z. Under uniform refinement, the gap scales as Zt+−Zt−=O(2−t/d)(Z0+−Z0−)Z_t^+-Z_t^-=O(2^-t/d)(Z_0^+-Z_0^-), requiring depth O(dlog(1/ϵ))O(d (1/ε)) to achieve target gap ϵε. Proof. We formalize the refinement scheme and then prove each claim. Fix a Voronoi-gated sum node n with scope S and bounded domain ΩS⊂ℝ|S| _S ^|S|. Let Vkk=1K\V_k\_k=1^K be the Voronoi cells intersected with ΩS _S (so each VkV_k is measurable and ⋃kVk=ΩS _kV_k= _S up to boundaries of measure zero). At refinement step t, we maintain a disjoint axis-aligned partition tP_t of ΩS _S into boxes, i.e., ΩS=⨄B∈tB _S= _B _tB (disjoint union). For each cell k, define the inner/outer approximations induced by tP_t: Vk−(t) V_k^-(P_t) :=⋃B∈t:B⊆VkB, := _B _t:\,B V_kB, Vk+(t) V_k^+(P_t) :=⋃B∈t:B∩Vk≠∅B. := _B _t:\,B∩ V_k≠ B. (3) By construction, Vk−(t)⊆Vk⊆Vk+(t)V_k^-(P_t) V_k V_k^+(P_t) for every t. Let pk(S)≥0p_k(x_S)≥ 0 denote the expert density (subcircuit) attached to cell k. Define the node-level cell integrals and their bounds Ik:=∫Vkpk(S)S,Ik,t−:=∫Vk−(t)pk(S)S,Ik,t+:=∫Vk+(t)pk(S)S. I_k:= _V_kp_k(x_S)\,dx_S, I_k,t^-:= _V_k^-(P_t)p_k(x_S)\,dx_S, I_k,t^+:= _V_k^+(P_t)p_k(x_S)\,dx_S. (4) Because tP_t is disjoint and Vk−(t)V_k^-(P_t) (resp. Vk+(t)V_k^+(P_t)) is a union of boxes from tP_t, we can write these integrals as sums of box integrals (no overlap): Ik,t−=∑B∈t:B⊆VkIntegrateBox(pk,B), I_k,t^-= _B _t:\,B V_k IntegrateBox(p_k,B), Ik,t+=∑B∈t:B∩Vk≠∅IntegrateBox(pk,B). I_k,t^+= _B _t:\,B∩ V_k≠ IntegrateBox(p_k,B). The node-level mixture integral at n is I:=∑k=1KπkIkI:= _k=1^K _kI_k and its bounds are It−:=∑kπkIk,t−I_t^-:= _k _kI_k,t^- and It+:=∑kπkIk,t+I_t^+:= _k _kI_k,t^+, where πk≥0 _k≥ 0 and ∑kπk=1 _k _k=1. Circuit-level bounds (Zt−,Zt+)(Z_t^-,Z_t^+) are obtained by propagating node-level bounds upward via the same sum/product rules as in Theorem 3.10. Claim (i): Validity. Zt−≤Z≤Zt+Z_t^-≤ Z≤ Z_t^+. We first establish node-level validity at every Voronoi-gated node n. For each k, since Vk−(t)⊆Vk⊆Vk+(t)V_k^-(P_t) V_k V_k^+(P_t) and pk≥0p_k≥ 0, monotonicity of integration gives Ik,t−≤Ik≤Ik,t+.I_k,t^-≤ I_k≤ I_k,t^+. Multiplying by πk≥0 _k≥ 0 and summing over k yields It−≤I≤It+I_t^-≤ I≤ I_t^+ at node n. For standard sum nodes, linearity and nonnegative weights preserve inequalities; for product nodes, decomposability implies integrals factor and multiplying nonnegative bounds preserves inequalities. Thus, by induction in reverse topological order (as in Theorem 3.10), every node integral is sandwiched by its bounds, in particular at the root Zt−≤Z≤Zt+Z_t^-≤ Z≤ Z_t^+. Claim (i): Monotonic tightening. By definition, t+1P_t+1 is obtained from tP_t by refining (bisecting) a subset of boxes and reclassifying. Refinement preserves disjointness and refines the partition: every box B′∈t+1B _t+1 is contained in some box B∈tB _t, and every B∈tB _t is the disjoint union of its descendants in t+1P_t+1. We show set inclusions: Vk−(t)⊆Vk−(t+1),Vk+(t+1)⊆Vk+(t).V_k^-(P_t) V_k^-(P_t+1), V_k^+(P_t+1) V_k^+(P_t). (5) For the inner sets: take any B∈tB _t with B⊆VkB V_k. Under refinement, B is either unchanged or replaced by disjoint children Br′⊆B\B _r\ B. In either case, every descendant Br′B _r satisfies Br′⊆B⊆VkB _r B V_k, hence all of B’s mass remains Inside and is included in Vk−(t+1)V_k^-(P_t+1). Therefore the union of Inside boxes can only expand, proving the first inclusion in (5). For the outer sets: take any B′∈t+1B _t+1 such that B′∩Vk≠∅B ∩ V_k≠ . Let B∈tB _t be its (unique) ancestor box with B′⊆B B. Then necessarily B∩Vk≠∅B∩ V_k≠ as well, hence B⊆Vk+(t)B V_k^+(P_t). Since B′B is contained in B, we have B′⊆Vk+(t)B V_k^+(P_t), and taking unions over all such B′B gives the second inclusion in (5). Now apply monotonicity of integration (using pk≥0p_k≥ 0) to (5): Ik,t−=∫Vk−(t)pk≤∫Vk−(t+1)pk=Ik,t+1−I_k,t^-= _V_k^-(P_t)p_k≤ _V_k^-(P_t+1)p_k=I_k,t+1^- Ik,t+1+=∫Vk+(t+1)pk≤∫Vk+(t)pk=Ik,t+.I_k,t+1^+= _V_k^+(P_t+1)p_k≤ _V_k^+(P_t)p_k=I_k,t^+. Weighting by πk≥0 _k≥ 0 and summing yields node-level monotonicity It−≤It+1−I_t^-≤ I_t+1^- and It+1+≤It+I_t+1^+≤ I_t^+. Finally, the circuit-level bounds follow by propagating inequalities through the circuit: standard sums preserve monotonicity by linearity with nonnegative weights, and products preserve monotonicity because all integrals are nonnegative and multiplication is monotone in each argument. Hence Zt−≤Zt+1−Z_t^-≤ Z_t+1^- and Zt+1+≤Zt+Z_t+1^+≤ Z_t^+. Claim (i): Convergence under vanishing boundary volume. Fix a Voronoi-gated node n and a cell k. Define the undecided region Uk,t:=Vk+(t)∖Vk−(t),U_k,t:=V_k^+(P_t) V_k^-(P_t), so that Vk−(t)⊆Vk⊆Vk+(t)V_k^-(P_t) V_k V_k^+(P_t) implies 0≤Ik,t+−Ik,t−=∫Uk,tpk(S)S.0≤ I_k,t^+-I_k,t^-= _U_k,tp_k(x_S)\,dx_S. (6) Assume that the boundary volume μ(Uk,t)→0μ(U_k,t)→ 0 as t→∞t→∞. Since pkp_k is integrable over ΩS _S (it is a density component on a bounded domain), the Lebesgue integral is absolutely continuous with respect to the measure: for any ε>0 >0 there exists δ>0δ>0 such that μ(A)<δμ(A)<δ implies ∫Apk<ε _Ap_k< . Applying this to A=Uk,tA=U_k,t yields ∫Uk,tpk→0 _U_k,tp_k→ 0, hence Ik,t+−Ik,t−→0I_k,t^+-I_k,t^-→ 0 by (6). Therefore Ik,t−→IkI_k,t^-→ I_k and Ik,t+→IkI_k,t^+→ I_k (squeezing with validity from (i)). Multiplying by πk _k and summing over k shows that the node-level bounds converge to the true node integral. Applying the same inductive argument up the circuit—using continuity of addition and multiplication on ℝ≥0R_≥ 0 and the fact that all intermediate quantities are finite—yields convergence at the root: Zt−→Z_t^-→ Z and Zt+→Z_t^+→ Z. The convergence-rate statement Zt+−Zt−=O(2−t/d)(Z0+−Z0−)Z_t^+-Z_t^-=O(2^-t/d)(Z_0^+-Z_0^-) under uniform refinement is an informal geometric bound that depends on regularity of the density and the surface-area-to-volume behavior of Voronoi facets. A sufficient condition is, for example, that each expert density is bounded on ΩS _S by ‖pk‖∞<∞\|p_k\|_∞<∞ and that the undecided region Uk,tU_k,t lies in a tubular neighborhood of the cell boundary whose thickness scales with the maximum box diameter. Under uniform refinement, after t steps the maximum box side length is O(L⋅2−t/d)O(L· 2^-t/d) where L is the initial domain diameter. The undecided region Uk,tU_k,t lies in a neighborhood of thickness O(2−t/d)O(2^-t/d) around the (d−1)(d-1)-dimensional boundary ∂Vk∂ V_k. Thus μ(Uk,t)=O(surface area×thickness)=O(2−t/d).μ(U_k,t)=O(surface area×thickness)=O(2^-t/d). Assuming bounded density ‖pk‖∞<∞\|p_k\|_∞<∞: Ik,t+−Ik,t−≤‖pk‖∞⋅μ(Uk,t)=O(2−t/d).I_k,t^+-I_k,t^-≤\|p_k\|_∞·μ(U_k,t)=O(2^-t/d). Summing over k and propagating through the circuit gives Zt+−Zt−=O(2−t/d)(Z0+−Z0−).Z_t^+-Z_t^-=O(2^-t/d)(Z_0^+-Z_0^-). To achieve gap ϵε, we need t≥dlog2((Z0+−Z0−)/ϵ)=O(dlog(1/ϵ))t≥ d _2((Z_0^+-Z_0^-)/ε)=O(d (1/ε)). ∎ Remark A.10 (Curse of dimensionality). The factor d in the refinement depth makes certified bounds impractical in high dimensions. For example, achieving ϵ=0.01ε=0.01 in d=10d=10 requires roughly 10log2(100)≈6610 _2(100)≈ 66 refinement levels, leading to exponential blowup. This motivates the HFV-PC construction in Section 3.2. A.3.3 Computational Complexity of Refinement Proposition A.11 (Per-Iteration Complexity). Each iteration of Algorithm 1 requires: • O(||⋅Kmax)O(|C|· K_ ) to compute certified bounds via bottom-up propagation. • O(K⋅d⋅2d)O(K· d· 2^d) to test box-polytope containment/intersection for K cells in dimension d. Total per-iteration cost: O(||⋅Kmax+K⋅d⋅2d)O(|C|· K_ +K· d· 2^d). Proof. Bound computation: The circuit has |||C| nodes. At each Voronoi-gated node, we integrate the expert subcircuit over at most K boxes. Assuming box integration over a subcircuit is linear in subcircuit size, the total cost is O(||⋅Kmax)O(|C|· K_ ). Box classification: To determine whether box B=∏i[li,ui]⊆ℝdB= _i[l_i,u_i] ^d satisfies B⊆VkB V_k or B∩Vk=∅B∩ V_k= , we use the half-space representation: Vk=⋂j≠k:(j−k)⊤≤cjk.V_k= _j≠ k\x:(c_j-c_k) x≤ c_jk\. For each half-space, compute extremal values of a⊤a x over box B. This takes O(d)O(d) per half-space. There are K−1K-1 half-spaces, so testing one box against one cell takes O(K⋅d)O(K· d). Testing all K cells takes O(K2⋅d)O(K^2· d). Alternatively, testing all 2d2^d corners of B against all K centroids takes O(2d⋅K⋅d)O(2^d· K· d). The dominant term is O(2d)O(2^d) for large d. ∎ Remark A.12. The 2d2^d factor in corner evaluations reinforces the curse of dimensionality. In practice, for d≤5d≤ 5, adaptive refinement is effective; for higher dimensions, HFV-PCs are preferable. A.3.4 Marginals and Conditionals The certified bound framework extends naturally to marginals and conditionals. Marginal bounds. To compute p(A)=∫p()A¯p(x_A)= p(x)\,dx_ A where A¯=∖A A=X _A, propagate bounds through the circuit while marginalizing out A¯ A. At each Voronoi-gated node: • If the scope S includes variables in A¯ A, restrict box approximations to dimensions being integrated (project boxes onto subspace A¯∩SX_ A∩ S) and integrate over those dimensions. • If the scope is entirely within A, no change needed. The result is bounds (p−(A),p+(A))(p^-(x_A),p^+(x_A)) satisfying p−(A)≤p(A)≤p+(A)p^-(x_A)≤ p(x_A)≤ p^+(x_A) for all Ax_A. Conditional bounds. For disjoint A,BA,B and p(A∣B)=p(A,B)/p(B)p(x_A _B)=p(x_A,x_B)/p(x_B), we have bounds on numerator and denominator. Assuming p−(B)>0p^-(x_B)>0, apply interval arithmetic: p−(A,B)p+(B)≤p(A∣B)≤p+(A,B)p−(B). p^-(x_A,x_B)p^+(x_B)≤ p(x_A _B)≤ p^+(x_A,x_B)p^-(x_B). These bounds are valid because division is monotone in the numerator and antitone in the denominator when all quantities are positive. Practical considerations. If p−(B)p^-(x_B) is very small or zero, the upper bound becomes loose or undefined. Refining the partition to tighten p−(B)p^-(x_B) away from zero is necessary for meaningful conditional bounds. A.4 Hierarchical Factorized Voronoi PCs A.4.1 Factorized Voronoi Cells The key idea of HFV-PCs is to replace a single high-dimensional Voronoi tessellation with a product of lower-dimensional tessellations that align with the circuit’s variable decomposition. Definition A.13 (Factorized Voronoi Partition). Let S=⨆i=1mSiX_S= _i=1^mX_S_i be a partition of scope S into m disjoint blocks. For each block i, choose centroids 1(i),…,Ki(i)⊂ℝ|Si|\c^(i)_1,…,c^(i)_K_i\ ^|X_S_i| and let Vki(i)ki=1Ki\V^(i)_k_i\_k_i=1^K_i be the induced Voronoi cells in ℝ|Si|R^|X_S_i|. Define the joint product cell indexed by =(k1,…,km)∈[K1]×⋯×[Km]k=(k_1,…,k_m)∈[K_1]×·s×[K_m]: V:=Vk1(1)×⋯×Vkm(m)⊆ℝ|S|.V_k:=V^(1)_k_1×·s× V^(m)_k_m ^|X_S|. The corresponding hard gate factors as: g(S)=[S∈V]=∏i=1m[Si∈Vki(i)]=∏i=1mgki(i)(Si).g_k(x_S)=I[x_S∈ V_k]= _i=1^mI[x_S_i∈ V^(i)_k_i]= _i=1^mg^(i)_k_i(x_S_i). Membership in V_k can be decided independently for each block. This is exactly the geometric analogue of decomposability. Proposition A.14 (Factorized Gate Decomposition). For a factorized cell V_k, the hard gate decomposes as g(S)=∏i=1mgki(i)(Si).g_k(x_S)= _i=1^mg^(i)_k_i(x_S_i). Proof. By definition, S∈Vx_S∈ V_k if and only if Si∈Vki(i)x_S_i∈ V^(i)_k_i for all i. Since the conditions involve disjoint variable sets, the indicator of their conjunction factors: [S∈V]=[⋂iSi∈Vki(i)]=∏i[Si∈Vki(i)].I[x_S∈ V_k]=I [ _i\x_S_i∈ V^(i)_k_i\ ]= _iI[x_S_i∈ V^(i)_k_i]. ∎ A.4.2 HFV-Gated Sum Nodes and Tractability Theorem A.15 (Tractability of HFV-PCs). Let C be an HFV-PC with |||C| nodes, maximum factorization degree m, and at most K Voronoi cells per factor. Then the partition function, all marginals, and all conditionals are computable exactly in time O(||Km)O(|C|K^m). Proof. We prove the partition function claim. The marginal and conditional claims follow by the same bottom-up evaluation used for standard smooth decomposable PCs. For each node n we consider In=∫fn(scope(n))scope(n)I_n= f_n(x_scope(n))\,dx_scope(n). If n is a leaf, InI_n is a univariate integral and is tractable by assumption. If n is a product node, decomposability implies that children have disjoint scopes, so the integral factors as In=∏c∈ch(n)IcI_n= _c (n)I_c. If n is a standard sum node, linearity gives In=∑cπcIcI_n= _c _cI_c. It remains to handle an HFV-gated sum node n with scope S=⨆i=1mSiX_S= _i=1^mX_S_i. Using (LABEL:eq:hfv_sum) and exchanging summation and integration we obtain In I_n =∫∑g(S)π∏i=1mpki(i)(Si)dS = _kg_k(x_S)\, _k _i=1^mp^(i)_k_i(x_S_i)\,dx_S =∑π∫∏i=1mgki(i)(Si)pki(i)(Si)dS. = _k _k _i=1^mg^(i)_k_i(x_S_i)\,p^(i)_k_i(x_S_i)\,dx_S. (7) Since the scopes SiX_S_i are disjoint and each factor depends only on Six_S_i, we apply Fubini to factor the integral, ∫∏i=1mgki(i)(Si)pki(i)(Si)dS _i=1^mg^(i)_k_i(x_S_i)\,p^(i)_k_i(x_S_i)\,dx_S (8) =∏i=1m∫gki(i)(Si)pki(i)(Si)Si. = _i=1^m g^(i)_k_i(x_S_i)\,p^(i)_k_i(x_S_i)\,dx_S_i. (9) Substituting gki(i)(Si)=[Si∈Vki(i)]g^(i)_k_i(x_S_i)=I[x_S_i∈ V^(i)_k_i] yields In=∑π∏i=1mIki(i),Iki(i)=∫Vki(i)pki(i)(Si)Si.I_n= _k _k _i=1^mI^(i)_k_i, I^(i)_k_i= _V^(i)_k_ip^(i)_k_i(x_S_i)\,dx_S_i. Each Iki(i)I^(i)_k_i is the integral of an HFV subcircuit over a lower-dimensional Voronoi cell, so we compute it recursively by the same argument. This establishes exact tractability. For complexity, at an HFV-gated sum node we compute ∑iKi _iK_i factor integrals and then evaluate the sum over at most ∏iKi≤Km _iK_i≤ K^m joint indices, each term requiring O(m)O(m) arithmetic operations. Summing over all circuit nodes yields O(||Km)O(|C|K^m) time. ∎ Cell-Restricted Integrals and the Recursion Base Case The recursion requires that we can evaluate integrals of the form ∫Vki(i)pki(i) _V^(i)_k_ip^(i)_k_i at every level. In an HFV-PC the same factorization property ensures that cell restrictions are handled locally at the appropriate scope. The recursion bottoms out at univariate leaves, where Voronoi cells are intervals and integration is straightforward. Proposition A.16 (Univariate Voronoi Cells). In ℝ1R^1, let c1<⋯<cKc_1<·s<c_K be centroids. The induced Voronoi tessellation consists of intervals Vk=(−∞,c1+c22]k=1(ck−1+ck2,ck+ck+12]1<k<K(cK−1+cK2,∞)k=K.V_k= cases(-∞, c_1+c_22]&k=1\\ ( c_k-1+c_k2, c_k+c_k+12]&1<k<K\\ ( c_K-1+c_K2,∞)&k=K. cases For standard univariate leaf families such as Gaussians, mixtures of Gaussians, exponentials, and bounded distributions, integrating over an interval is tractable by closed form CDF evaluation or simple numerical routines. Definition A.17 (Binary HFV-PC). A binary HFV-PC over variables X with vtree T is a probabilistic circuit in which each internal vtree node v with children vLv_L and vRv_R induces HFV-gated sum nodes over scope v=vL⊔vRX_v=X_v_L _v_R. Each HFV-gated sum uses the binary partition (vL,vR)(X_v_L,X_v_R) and computes: f(v)=∑kL,kRgkL(L)(vL)gkR(R)(vR)πkL,kRpkL(L)(vL)pkR(R)(vR).f(x_v)= _k_L,k_Rg^(L)_k_L(x_v_L)\,g^(R)_k_R(x_v_R)\, _k_L,k_R\,p^(L)_k_L(x_v_L)\,p^(R)_k_R(x_v_R). A.5 Learning via Soft Gating A.5.1 Differentiability and Gradients Proposition A.18 (Differentiability). Let C be a soft Voronoi gated PC with parameters Θ=ki(i),π,θleaf =\\c^(i)_k_i\,\ _k\, _leaf\. For any finite temperature α>0α>0, the likelihood p(;Θ,α)p(x; ,α) is differentiable with respect to all parameters. Proof. Each soft gate wk(;α)=exp(−α‖−k‖2)∑jexp(−α‖−j‖2)w_k(u;α)= (-α\|u-c_k\|^2) _j (-α\|u-c_j\|^2) is a composition of smooth operations: squared Euclidean distance (polynomial), exponential, and softmax normalization (ratio of positive smooth functions). The circuit output is built from sums and products of smooth gates and leaf densities, hence differentiable in all parameters. ∎ Proposition A.19 (Centroid Gradient). For soft gate wk(;α)w_k(u;α), the gradient with respect to centroid kc_k is: ∇kwk(;α)=2αwk(;α)(1−wk(;α))(−k). _c_kw_k(u;α)=2α\,w_k(u;α) (1-w_k(u;α) )\,(u-c_k). Proof. Let dj=‖−j‖2d_j=\|u-c_j\|^2 and Z=∑jexp(−αdj)Z= _j (-α d_j). Then wk=exp(−αdk)/Zw_k= (-α d_k)/Z. Taking the derivative with respect to kc_k: ∇kwk=∇k[exp(−αdk)]Z−exp(−αdk)∇k[Z]Z2. _c_kw_k= _c_k[ (-α d_k)]Z- (-α d_k)\, _c_k[Z]Z^2. We have ∇kdk=−2(−k) _c_kd_k=-2(u-c_k), so: ∇k[exp(−αdk)]=2αexp(−αdk)(−k). _c_k[ (-α d_k)]=2α (-α d_k)(u-c_k). Since only the k-th term in Z depends on kc_k: ∇k[Z]=2αexp(−αdk)(−k). _c_k[Z]=2α (-α d_k)(u-c_k). Substituting: ∇kwk _c_kw_k =2αexp(−αdk)(−k)Z−exp(−αdk)⋅2αexp(−αdk)(−k)Z2 = 2α (-α d_k)(u-c_k)Z- (-α d_k)· 2α (-α d_k)(u-c_k)Z^2 =2αwk(1−wk)(−k). =2α\,w_k\,(1-w_k)\,(u-c_k). ∎ The factor wk(1−wk)w_k(1-w_k) is largest when wk≈1/2w_k≈ 1/2 (near decision boundaries where the model is uncertain). The gradient magnitude peaks in regions of ambiguity and vanishes where routing is confident, naturally focusing centroid updates on improving contested region geometry. A.5.2 Soft-to-Hard Convergence Theorem A.20 (Soft-to-Hard Convergence). Let Vk\V_k\ be Voronoi cells induced by centroids k\c_k\ and let gk()=[∈Vk]g_k(u)=I[u∈ V_k] be the hard gate. Let wk(;α)w_k(u;α) be the soft gate. Then: (i) For any u not on a Voronoi boundary, limα→∞wk(;α)=gk() _α→∞w_k(u;α)=g_k(u). (i) If k∗()=argminj‖−j‖k^*(u)= _j\|u-c_j\| and margin γ():=minj≠k∗(‖−j‖2−‖−k∗‖2)>0γ(u):= _j≠ k^*(\|u-c_j\|^2-\|u-c_k^*\|^2)>0, then 1−wk∗(;α)≤(K−1)exp(−αγ()).1-w_k^*(u;α)≤(K-1) (-αγ(u)). (i) If p is integrable, then limα→∞∫wk(;α)p()d=∫Vkp()d _α→∞ w_k(u;α)p(u)\,du= _V_kp(u)\,du. Proof. (i) Pointwise convergence: Let u be a point not on any Voronoi boundary. Then there exists unique nearest centroid k∗c_k^* with dk∗:=‖−k∗‖2<dj:=‖−j‖2d_k^*:=\|u-c_k^*\|^2<d_j:=\|u-c_j\|^2 for all j≠k∗j≠ k^*. Rewrite: wk∗(;α)=11+∑j≠k∗exp(−α(dj−dk∗)).w_k^*(u;α)= 11+ _j≠ k^* (-α(d_j-d_k^*)). Since dj−dk∗>0d_j-d_k^*>0 for all j≠k∗j≠ k^*, as α→∞α→∞, each term exp(−α(dj−dk∗))→0 (-α(d_j-d_k^*))→ 0 exponentially. Therefore: limα→∞wk∗(;α)=1=gk∗(). _α→∞w_k^*(u;α)=1=g_k^*(u). For j≠k∗j≠ k^*, we have limα→∞wj(;α)=0=gj() _α→∞w_j(u;α)=0=g_j(u). (i) Exponential rate: Define margin γ()=minj≠k∗(dj−dk∗)>0γ(u)= _j≠ k^*(d_j-d_k^*)>0. Then for all j≠k∗j≠ k^*: exp(−α(dj−dk∗))≤exp(−αγ()). (-α(d_j-d_k^*))≤ (-αγ(u)). Summing: ∑j≠k∗exp(−α(dj−dk∗))≤(K−1)exp(−αγ()). _j≠ k^* (-α(d_j-d_k^*))≤(K-1) (-αγ(u)). Therefore: wk∗≥11+(K−1)exp(−αγ),w_k^*≥ 11+(K-1) (-αγ), implying 1−wk∗≤(K−1)exp(−αγ()).1-w_k^*≤(K-1) (-αγ(u)). (i) Integral convergence: We have 0≤wk≤10≤ w_k≤ 1 for all ,αu,α. Pointwise limit gives wk→gkw_k→ g_k except on measure-zero boundaries. By dominated convergence: limα→∞∫wk(;α)p()=∫gk()p()=∫Vkp(). _α→∞ w_k(u;α)p(u)\,du= g_k(u)p(u)\,du= _V_kp(u)\,du. ∎ The margin γ()γ(u) quantifies how much closer u is to its nearest centroid compared to second-nearest. Large margins yield faster convergence. Small margins (near decision boundaries) require larger α for hard-like behavior. Algorithm 3 Soft Gating Training with Annealing 0: Dataset =(n)n=1ND=\x^(n)\_n=1^N, schedule αtt=1T\ _t\_t=1^T, learning rate η 0: Trained parameters Θ 1: Initialize: 2: Centroids ki(i)\c^(i)_k_i\ via k-means (per factor for HFV) 3: Mixture weights π\ _k\ uniformly on simplex 4: Leaf parameters θleaf _leaf randomly or from prior 5: for t=1t=1 to T do 6: Set inverse temperature α←αtα← _t 7: for each minibatch ℬ⊂B do 8: ℒ←−1|ℬ|∑∈ℬlogp(;Θ,α)L←- 1|B| _x p(x; ,α) 9: Θ←Θ−η∇Θℒ ← -η\, _ L 10: Enforce π on simplex (projection or softmax) 11: end for 12: end for 13: return Θ (a) Alphabets (2D) (b) Checkerboard (2D) (c) Pinwheel (2D) (d) Spiral (2D) (e) Bent Lissajous (3D) (f) Interlocked Circles (3D) (g) Knotted (3D) (h) Twisted Eight (3D) Figure 4: Synthetic 2D/3D density estimation benchmarks. Top row: 2D datasets with diverse local geometry and disconnected support. Bottom row: 3D manifold-like datasets exhibiting crossings, interlocks, and knotting. These benchmarks stress-test whether geometry-aware routing can specialize locally while maintaining reliable inference (VT via certification; HFV via alignment). Appendix B Synthetic Dataset Construction Figure 4 visualizes the synthetic 2D/3D benchmarks used in our experiments. These datasets were designed to emphasize geometric structure (e.g., curved manifolds, crossings, knots, and disconnected supports) that can be difficult to capture with input-independent mixture weights, making them well-suited for evaluating geometry-aware routing. All datasets were generated following a unified protocol: (1) generate raw samples from a geometric structure as defined below, (2) add Gaussian noise (0,σ2)N(0,σ^2I) with σ=0.01σ=0.01, (3) standardize to zero mean and unit variance per dimension, (4) split into 10k/5k/5k10k/5k/5k train/val/test sets. B.1 Two-Dimensional Datasets CheckerBoard. Nine Gaussian clusters at grid positions −1.5,0,1.52\-1.5,0,1.5\^2. Each cluster samples uniformly from a square [k−0.6,k+0.6][c_k-0.6,c_k+0.6] then adds Gaussian noise σ=0.15σ=0.15 and clips to bounds. Pinwheel. Five radial arms at angles 2πk/52π k/5 for k=0,…,4k=0,…,4. Each sample: radius r∼(1.0,0.32)r (1.0,0.3^2), angular offset δθ∼(0,0.22)δθ (0,0.2^2), then convert to Cartesian (rcos(θk+δθ),rsin(θk+δθ))(r ( _k+δθ),r ( _k+δθ)). Spiral. Two interleaved Archimedean spirals. For each spiral: sample θ∼Uniform(0,2π)θ (0,2π), set θ←θ⋅2πθ← θ· 2π (arc-length correction), radius r=2θr=2θ, convert to (±rcosθ,±rsinθ)(± r θ,± r θ) (opposite signs for two spirals), add noise σ=0.1σ=0.1. Alphabet. Samples uniformly from pixel-based 7×57× 5 binary grids representing uppercase letters arranged spatially. The results reported are for the letter W. Active pixels are converted to continuous coordinates (cell size 0.20.2), positioned in a grid layout with letter gap 0.40.4. B.2 Three-Dimensional Datasets All 3D datasets were generated by sampling a parameter t∼Uniform([−π,π])t ([-π,π]), apply a parametric curve as defined below, scaling by 44, adding noise σ=0.01σ=0.01, foloowing (Sidheekh et al., 2022, 2023) BentLissajous. Lissajous curve: (x,y,z)=(sin(2t),cos(t),cos(2t))(x,y,z)=( (2t), (t), (2t)). InterlockedCircles. Two circles in orthogonal planes: Circle 1 in xyxy-plane (sint,cost,0)( t, t,0), Circle 2 in xzxz-plane (1+sint,0,cost)(1+ t,0, t). Knotted. Trefoil knot: (x,y,z)=(sint+2sin2t,cost−2cos2t,sin3t)(x,y,z)=( t+2 2t, t-2 2t, 3t). TwistedEight. Two circles in orthogonal planes: (sint,cost,0)( t, t,0) and (2+sint,0,cost)(2+ t,0, t). Appendix C Implementation Details All models were implemented in PyTorch using the CirKit package (Lab, 2024) and trained with Adam (lr 0.010.01, batch size 500500, 100100 epochs) on a single 2424GB NVIDIA L4 GPU. VT and HFV gating were implemented as drop-in replacements for sum layers, with learnable centroids initialized via k-means and optimized jointly with circuit parameters. We annealed the soft-gating inverse temperature linearly, α:1→50α:1→ 50, to transition from smooth routing to near-hard assignments. For VT models, evaluation and model selection used the certified normalization bounds produced by our box-based inference routine (reporting the log-likelihood lower bound), whereas HFV models preserved exact tractable inference throughout. For all models, the best performing epoch in terms of validation performance was saved and loaded for testing.