Paper deep dive
Behavior Learning (BL): Learning Hierarchical Optimization Structures from Data
Zhenyao Ma, Yue Liang, Dongxu Li
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 90%
Last extracted: 7/20/2026, 2:59:08 PM
Summary
The paper introduces Behavior Learning (BL), a machine learning framework that learns interpretable and identifiable hierarchical optimization structures from data. Inspired by behavioral science, BL models responses as solutions to Utility Maximization Problems (UMPs) composed into a compositional utility function. It offers three architectural variants: BL(Single), BL(Shallow), and BL(Deep), with the latter supporting deep hierarchical compositions. The Identifiable BL (IBL) variant guarantees identifiability. The framework unifies predictive performance, intrinsic interpretability, and identifiability, applicable to domains like economics, physics, and biology.
Entities (9)
Relation Signals (9)
Behavior Learning → isbasedon → Utility Maximization Problem
confidence 95% · BL parameterizes a compositional utility function built from intrinsically interpretable modular blocks... Each block represents and can be written in symbolic form as a utility maximization problem (UMP)
Behavior Learning → hasvariant → Identifiable Behavior Learning
confidence 92% · Its smooth and monotone variant (IBL) guarantees identifiability.
Behavior Learning → hasarchitecture → BL(Shallow)
confidence 90% · BL admits hierarchical structure, mainly in three architectural variants: ... BL(Shallow)...
Behavior Learning → hasarchitecture → BL(Deep)
confidence 90% · BL admits hierarchical structure, mainly in three architectural variants: ... and BL(Deep)...
Behavior Learning → hasarchitecture → BL(Single)
confidence 90% · BL admits hierarchical structure, mainly in three architectural variants: BL(Single)...
BL(Deep) → models → Hierarchical Optimization Structures
confidence 90% · BL(Deep), a deep hierarchical composition of multiple blocks. The latter two model, and can be symbolically interpreted as, hierarchical optimization structures.
Behavior Learning → hasproperty → Universal Approximation
confidence 88% · Theoretically, we establish the universal approximation property of BL
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Inspired by behavioral science, we propose Behavior Learning (BL), a novel general-purpose machine learning framework that learns interpretable and identifiable optimization structures from data, ranging from single optimization problems to hierarchical compositions. It unifies predictive performance, intrinsic interpretability, and identifiability, with broad applicability to scientific domains involving optimization. BL parameterizes a compositional utility function built from intrinsically interpretable modular blocks, which induces a data distribution for prediction and generation. Each block represents and can be written in symbolic form as a utility maximization problem (UMP), a foundational paradigm in behavioral science and a universal framework of optimization. BL supports architectures ranging from a single UMP to hierarchical compositions, the latter modeling hierarchical optimization structures. Its smooth and monotone variant (IBL) guarantees identifiability. Theoretically, we establish the universal approximation property of BL, and analyze the M-estimation properties of IBL. Empirically, BL demonstrates strong predictive performance, intrinsic interpretability and scalability to high-dimensional data. Code: this https URL ; install via pip install blnetwork.
Tags
Links
- Source: https://arxiv.org/abs/2602.20152v1
- Canonical: https://arxiv.org/abs/2602.20152v1
Trouble viewing inline? Open PDF directly →
Full Text
187,690 characters extracted from source content.
Expand or collapse full text
Behavior Learning (BL): Learning Hierarchical Optimization Structures from Data Zhenyao Ma Xiamen University &Yue Liang University of Tübingen &Dongxu Li Xi’an Jiaotong University Correspondence to: zhenyaoma@stu.xmu.edu.cn; Code contact: yue.liang@student.uni-tuebingen.de Abstract Inspired by behavioral science, we propose Behavior Learning (BL), a novel general-purpose machine learning framework that learns interpretable and identifiable optimization structures from data, ranging from single optimization problems to hierarchical compositions. It unifies predictive performance, intrinsic interpretability, and identifiability, with broad applicability to scientific domains involving optimization. BL parameterizes a compositional utility function built from intrinsically interpretable modular blocks, which induces a data distribution for prediction and generation. Each block represents and can be written in symbolic form as a utility maximization problem (UMP), a foundational paradigm in behavioral science and a universal framework of optimization. BL supports architectures ranging from a single UMP to hierarchical compositions, the latter modeling hierarchical optimization structures. Its smooth and monotone variant (IBL) guarantees identifiability. Theoretically, we establish the universal approximation property of BL, and analyze the M-estimation properties of IBL. Empirically, BL demonstrates strong predictive performance, intrinsic interpretability and scalability to high-dimensional data. Code: MoonYLiang/Behavior-Learning on GitHub; install via pip install blnetwork. Figure 1: Behavior Learning (BL). (a) Human behavior modeled as a UMP. (b) Learning scheme of BL, where CompU denotes the compositional utility function. (c) BL offers intrinsic interpretability (via symbolic form as a hierarchical optimization structure), identifiability, and inference capability. (d) Three architectural variants of BL, from single UMP to deep compositions. 1 Introduction Scientific research often grapples with phenomena that resist precise formalization (Anderson, 1972; Mitchell, 2009), including human and social domains (Simon, 1955; Arthur, 2009). Such phenomena are difficult to predict and even harder to falsify through theory alone. Interpretable machine learning (Interpretable ML) (Molnar, 2020), with its powerful approximation capabilities and built-in transparency, offers a promising alternative for modeling such phenomena. Yet a long-standing tension remains unresolved: model predictive performance and intrinsic interpretability often trade off—a challenge commonly known as the performance–interpretability trade-off (Arrieta et al., 2020). High-performing models such as deep neural networks (LeCun et al., 2015) typically lack transparency, while intrinsically interpretable models struggle to capture complex nonlinear patterns. Some efforts have been made to mitigate the performance–interpretability trade-off. For example, Hastie (2017); Alvarez Melis & Jaakkola (2018); Angelino et al. (2018); Nori et al. (2019); Koh et al. (2020); Agarwal et al. (2021); Kraus et al. (2024); Liu et al. (2024b); Plonsky et al. (2025) demonstrate varied strengths. However, two fundamental limitations remain, restricting their scientific applicability. (i) Insufficient alignment with scientific theories. Most approaches focus on extending existing machine learning methods to achieve interpretability, rather than developing a scientifically grounded framework (e.g., based on optimization problems or differential equations). This often hinders alignment with scientific theories and limits the ability to extract scientific knowledge from learned models (Roscher et al., 2020; Bereska & Gavves, 2024; Longo et al., 2024). (i) Non-uniqueness of interpretations. Most models are non-identifiable—their interpretations are not uniquely determined by observable predictions in a mathematical sense (Ran & Hu, 2017; Méloux et al., 2025). As a result, such models cannot support reliable estimation of ground-truth parameters (Newey & McFadden, 1994; Van der Vaart, 2000), and may even lack Popperian falsifiability (Popper, 2005), ultimately limiting their scientific credibility. These limitations naturally raise a key question: can we design an interpretable ML framework that mitigates the performance–interpretability trade-off while being scientifically grounded and identifiable? Inspired by behavioral science, we propose Behavior Learning (BL): a general-purpose machine learning framework that learns interpretable and identifiable (hierarchical) optimization structures from data. It unifies high predictive performance, intrinsic interpretability, and identifiability, with broad applicability to scientific domains involving optimization. As illustrated in Figure 1, BL builds on one of the most fundamental paradigms in behavioral science—utility maximization—which posits that human behavior arises from solving a utility maximization problem (UMP) (Samuelson, 1948; Debreu, 1959; Mas-Colell et al., 1995). Motivated by this paradigm, BL learns interpretable optimization structures from data. It models responses (y) as drawn from a probability distribution induced by a UMP or a composition of multiple interacting UMPs. This distribution is parameterized by a compositional utility function BL(,)BL(x,y), constructed from intrinsically interpretable modular blocks ℬ(,)B(x,y). Each block is a learnable penalty-based formulation that represents an optimization problem (UMP), which can be written in symbolic form and offers transparency comparable to linear regression. BL admits hierarchical structure, mainly in three architectural variants: BL(Single), defined by a single block; BL(Shallow), a moderately layered composition of blocks; and BL(Deep), a deep hierarchical composition of multiple blocks. The latter two model, and can be symbolically interpreted as, hierarchical optimization structures. All variants are trained end-to-end to induce a conditional Gibbs distribution for prediction and generation. By refining the penalty functions in each block into smooth and monotone forms, we develop Identifiable BL (IBL), the identifiable variant of BL. Under mild conditions, IBL guarantees unique intrinsic interpretability. This property ensures the scientific credibility of its explanations and further supports recovery of the ground-truth model under appropriate conditions. While motivated by behavioral science, BL is not domain-specific. It applies broadly to any scientific domain where observed outcomes arise as solutions to optimization problems—such as macroeconomics (Ramsey, 1928; Ljungqvist & Sargent, 2018), statistical physics (Gibbs, 1902; Landau & Lifshitz, 2013), or evolutionary biology (Wright et al., 1932; Fisher, 1999). This generality is supported by a key theoretical insight (Theorem 2.2): any optimization problem can be equivalently written as a UMP. This makes BL a general-purpose modeling framework for data-driven inverse optimization (Ahuja & Orlin, 2001) across diverse scientific disciplines. We study BL both theoretically and empirically. Theoretically, we show that both BL and IBL admit universal approximation under mild assumptions (Section 2.2). For IBL, we further establish its M-estimation properties (Section 2.3), including identifiability, consistency, universal consistency, asymptotic normality, and asymptotic efficiency. Empirically, we evaluate BL across four tasks. Standard prediction tasks (Section 3.1) demonstrate its strong predictive performance. A qualitative case study (Section 3.2) illustrates its intrinsic interpretability. Prediction on high-dimensional inputs (Section 3.3) further demonstrates its scalability to high-dimensional data. Further discussion and related work are provided in Section 5 and Section 6, respectively. We also provide guidance on how to scientifically explain BL(Deep) and architectural details in Section 4 and Section A, respectively. Overall, our key contributions are threefold. (i) We propose Behavior Learning (BL), a novel general-purpose machine learning framework inspired by behavioral science, which unifies high predictive performance, intrinsic interpretability, identifiability, and scalability. (i) For scientific research, BL offers a scientifically grounded, interpretable, and identifiable machine learning approach for modeling complex phenomena that defy precise formalization. BL applies broadly to scientific disciplines involving optimization. (i) At the paradigm level, BL learns from data the optimization structure of either a single optimization problem or a hierarchical composition of problems through distributional modeling, contributing a new methodology to data-driven inverse optimization. 2 Behavior Learning (BL) 2.1 Utility Maximization Problem (UMP) The modeling of human behavior, particularly in behavioral science and decision theory, often begins with the assumption that observed outcomes arise from a latent optimization process. A canonical formulation of this idea is the Utility Maximization Problem (UMP) (Mas-Colell et al., 1995), in which an agent selects actions ∈y in response to contextual features ∈x by solving: max∈U(,)s.t.(,)≤0,(,)=0 _y \,U(x,y) .t. (x,y)≤ 0,\;\;T(x,y)=0 (1) Here, U(⋅)U(·) denotes a subjective utility function encoding the agent’s internal preferences or goals. The inequality constraint (⋅)C(·) captures resource constraints, while the equality constraint (⋅)T(·) encodes either endogenous belief consistency or exogenous conservation laws. The UMP can be recast as a cost–benefit framework, where the agent trades off utility gains against constraint violations. Formally, under mild regularity conditions, it admits an unconstrained penalty reformulation at the level of local optimality (Han & Mangasarian, 1979), as formalized below. Theorem 2.1 (Local Exact Penalty Reformulation for UMP). Let ⊂ℝdxX ^d_x and ⊂ℝdyY ^d_y be nonempty compact sets, and let U:×→ℝU:X×Y , :×→ℝmC:X×Y ^m, and :×→ℝpT:X×Y ^p be C1C^1. Assume that for a given ∈x , the Han–Mangasarian constraint qualification holds at any strict local maximizer ⋆y of the UMP. Then there exist λ0>0 _0>0, λ1∈ℝ++m _1 ^m_++, and λ2∈ℝ++p _2 ^p_++ such that ⋆y is a local maximizer of max∈λ0ϕ(U(,))−λ1⊤ρ((,))−λ2⊤ψ((,)). _y \; _0\,φ\! (U(x,y) )- _1 ρ\! (C(x,y) )- _2 ψ\! (T(x,y) ). (2) Here ϕ:ℝ→ℝφ:R is strictly increasing and C1C^1, ρ(z):=maxz,0ρ(z):= \z,0\, and ψ(z):=|z|ψ(z):=|z|. The proof is provided in Appendix B.1. This unconstrained reformulation offers greater tractability for both theoretical analysis and model training. While motivated by behavioral modeling, the UMP formulation is not domain-specific. It applies to any setting where observed outcomes are solutions to (explicit or latent) optimization problems. This is because any optimization problem can be equivalently formulated as a UMP. We state this in the following result, while the formal statement and proof are provided in Appendix B.1. Theorem 2.2 (Universality of UMP). Any optimization problem of the form max∈f(,) _y f(x,y) or min∈f(,) _y f(x,y), subject to equality and inequality constraints, is equivalent to a UMP. 2.2 BL Architecture Figure 1(b–d) illustrates the architecture of BL. We consider samples (,)∼(x,y) , where ∈ℝdx ^d denotes contextual features and y is the response, represented as (disc,cont)∈disc×ℝmc(y^disc,y^cont) _disc×R^m_c, capturing its hybrid structure. Responses are assumed to be stochastically generated by solving multiple interacting UMPs, each with a penalty-based formulation, which together compose a compositional utility function BL(,)BL(x,y). On this basis, we model the data using a conditional Gibbs distribution (Gibbs, 1902) parameterized by BLΘ(,)BL_ (x,y): pτ(∣;Θ)=exp(BLΘ(,)/τ)Zτ(;Θ),Zτ(;Θ)=∫exp(BLΘ(,′)/τ)′p_τ(y ; )= \! (BL_ (x,y)/τ )Z_τ(x; ), Z_τ(x; )= _Y\! \! (BL_ (x,y )/τ )\,dy (3) Here the temperature parameter τ>0τ>0 controls the randomness of the response. As τ→0τ→ 0, the distribution in equation 3 converges to a Dirac measure supported on argmaxBL(,) _yBL(x,y), thereby recovering the deterministic best response obtained by solving the composed UMPs. Model Structure of BL(,)BL(x,y). To represent the composition of multiple UMPs, we build BL(,)BL(x,y) by composing fundamental modular blocks ℬ(,)B(x,y). Each block provides a penalty-based formulation of a single UMP, and together they yield the overall compositional utility function. Motivated by Theorem 2.1, we parameterize ℬ(,)B(x,y) as ℬ(,;θ):=λ0⊤ϕ(UθU(,))−λ1⊤ρ(θC(,))−λ2⊤ψ(θT(,))B(x,y;θ):= _0^\! φ (U_ _U(x,y) )- _1 ρ (C_ _C(x,y) )- _2 ψ (T_ _T(x,y) ) (4) where θ:=(λ0,λ1,λ2,θU,θC,θT)θ:=( _0, _1, _2, _U, _C, _T) denotes the complete set of learnable parameters. Following Theorem 2.1, ϕφ is an increasing function; ρ penalizes inequality violations; and ψ captures symmetric deviations. Each block can be written as a well-defined UMP. We then compose BL(,)BL(x,y) from multiple ℬB-blocks through hierarchical composition to improve its representational power for optimization structures, yielding three main architectural variants, as illustrated in Figure 1(d). 1. BL(Single) applies a single instance of ℬ(,)B(x,y) as defined in equation 4, without any additional layers. It can be viewed as learning a single UMP, and offers maximal interpretability. 2. BL(Shallow) uses ℬ(,)B(x,y) as the fundamental modular block to construct a shallow network. It introduces one or two intermediate layers of computation. Each layer ℓB_ stacks multiple parallel ℬℓ,iB_ ,i blocks to produce a vector in ℝdℓR^d_ , i.e., ℓ(,;θℓ):=[ℬℓ,1(,;θℓ,1),…,ℬℓ,dℓ(,;θℓ,dℓ)]⊤B_ (x,y; _ ):=[B_ ,1(x,y; _ ,1),…,B_ ,d_ (x,y; _ ,d_ )] . The output of ℓB_ is directly fed into the next ℓ+1B_ +1, and only the final output is passed through a learnable affine transformation. 3. BL(Deep) extends the BL(Shallow) architecture to more than two layers, enabling richer hierarchical compositions of UMPs while maintaining the same recursive structure. As before, only the final output is affine transformed. The overall structure of BL(Shallow) and BL(Deep) can be expressed in a unified form, where the shallow case corresponds to L≤2L≤ 2 and the deep case to L>2L>2: BL(,):=L⋅L(⋯2(1(,))⋯)BL(x,y):=W_L·B_L (·sB_2(B_1(x,y))·s ) (5) Learning Objective. The response y may contain both discrete and continuous components. For discrete responses, we directly apply cross-entropy (Kullback & Leibler, 1951) on discy^disc. For continuous responses, since the compositional utility function is analogous to an energy function (LeCun et al., 2006), we employ denoising score matching (Vincent, 2011) on conty^cont. The final objective combines the two with nonnegative weights γd,γc _d, _c: ℒ(θ)=γd[−logpτ(disc∣)]+γc∥∇~contlogpτ(~cont∣)+σ−2(~cont−cont)∥2L(θ)= _d\,E [- p_τ(y^disc ) ]+ _c\,E \| _ y^cont p_τ( y^cont )+σ^-2( y^cont-y^cont) \|^2 (6) Implementation Details. Here, we describe the key implementation choices for the general form of BL, taken as defaults unless otherwise noted. Further details are provided in Appendix A.3. • Function Instantiation. Following equation 4, we instantiate the function ℬ(,)B(x,y) as ℬ(,)=λ0⊤tanh(u(,))−λ1⊤ReLU(c(,))−λ2⊤|t(,)|B(x,y)= _0^\! (p_u(x,y) )- _1^\! ReLU (p_c(x,y) )- _2^\! |p_t(x,y) | (7) where u,c,tp_u,p_c,p_t are polynomial feature maps of bounded degree, providing interpretable representations of utility, inequality, and equality terms, respectively. The bounded tanh reflects the principle of diminishing marginal utility (Jevons, 2013), a commonly assumed principle in behavioral science, while ReLUReLU and |⋅||·| introduce soft penalties for constraint violations. • Polynomial Maps. In BL(Single), the structure of polynomial maps is optional. In BL(Shallow) and BL(Deep), each ℬB-block employs affine transformations as its polynomial maps, with higher-degree and interaction terms omitted by default for computational efficiency. • Skip Connections. For deep variants, skip connections can be optionally introduced to improve representational efficiency. More detailed architectural descriptions for this section are provided in Appendix A. Theoretical Guarantees. Under the given architecture, the BL framework has universal approximation power: it can approximate any continuous conditional distribution arbitrarily well, provided that BL has sufficient capacity, as stated below. The proof is given in Appendix B.2. Theorem 2.3 (Universal Approximation of BL). Let ⊂ℝdX ^d and ⊂ℝmY ^m be compact sets, and let p⋆(∣)p (y ) be any continuous conditional density such that p⋆(∣)>0p (y )>0 for all (,)∈×(x,y) ×Y. Then for any τ>0τ>0 and ε>0 >0, there exists a finite BL architecture (with depth and width depending on ε ) and a parameter θ⋆θ such that the Gibbs distribution in equation 3 satisfies sup∈KL(p⋆(⋅∣)∥pτ(⋅∣;θ⋆))<ε. _x KL (p (· )\,\|\,p_τ(· ;θ ) )< . (8) Figure 2: (a) Visualization and symbolic form of BL(Single) trained on the Boston Housing dataset, modeling the UMP (maxUs.t.≤0,=0 U\;\;s.t.\;C≤ 0,\;T=0) of a representative buyer in Boston housing (details in Section 3.2). Top: computational graphs of the polynomials inside the three penalty functions—tanh (preference), ReLUReLU (budget), and |⋅||·| (belief). Each graph is respectively centered on tanh−1(U) ^-1(U), C, and T from left to right, with surrounding nodes representing input features. Directed edges (shown only if coefficient ≥0.3≥ 0.3) indicate how each feature contributes to the corresponding term. Bottom: approximate symbolic formulation of the trained BL model as a UMP. (b) The BL[2,1] architecture. Layer 1 identifies two key micro-level preference types: the Economic-sensitive Buyer and the Location-sensitive Buyer. Layer 2 aggregates these two components into an effective representative buyer. (c) The BL(Deep) [5,3,1] architecture. Layer 1 recovers five distinct micro-level housing preference types. Layer 2 identifies three macro-level trade-off types capturing different ways these primitive preferences interact. Layer 3 aggregates them into the overall representative buyer. Table 10 provides detailed descriptions of each type. BL(Deep) provides a hierarchical explanation consistent with the coarse-graining principle (Kadanoff, 1966) in statistical physics, reconstructing the full micro-to-macro optimization hierarchy. In addition, the preference and trade-off patterns uncovered by BL(Deep) are well documented in the classical economics literature (see Table 11). (d) BL can be applied to a broad class of hierarchical optimization structures in science, including hierarchical need structures, hierarchical social–organizational structures, and renormalization-style coarse-grained structures in physics. Interpretability. Alongside its expressive power, BL also exhibits strong intrinsic interpretability. • Each ℬB-block can be expressed in symbolic form as an optimization problem (UMP): the tanh term defines the objective, the ReLUReLU term corresponds to an inequality constraint, and the absolute-value term corresponds to an equality constraint. Thus, BL(Single) can be directly expressed as a symbolic UMP, whereas deeper architectures can be interpreted as compositions of UMPs, with each block retaining interpretability. • The polynomial basis ensures a level of transparency comparable to linear regression, as both objectives and constraints can be represented as linear combinations of polynomial features. It can further be visualized as a computational graph (Figure 7), in which each input’s influence on every ℬB-block is traceable through compositional pathways. • BL(Deep) composes ℬB-blocks in a layered manner, forming a hierarchical optimization structure. Interpretation proceeds in a bottom-up fashion, where the relation between any two consecutive layers can be viewed as aggregation or coarse-grained observation. Overall, the interpretive pathway is: raw input features → micro-level optimization blocks → macro-level aggregation or coarse-grained behavioral constructs → macro-level optimization systems. Section 4 provides a detailed description of this interpretation procedure. • BL also offers multiple architectural degrees of freedom that provide flexibility but simultaneously affect the resulting interpretability. In deep variants, skip connections introduce cross-layer dependency structures that are modeled in statistical physics (Yang & Schoenholz, 2017). Replacing polynomial maps with affine transformations preserves the underlying optimization semantics but reduces symbolic granularity, yielding a more qualitative rather than symbolic interpretation of each block. • BL can be interpreted as a single UMP when the final layer contains only one ℬB-block, since all lower-layer structures aggregate into a unified optimization problem. When the final layer contains multiple ℬB-blocks, BL corresponds to a linear trade-off among multiple optimization problems. 2.3 Identifiable Behavior Learning (IBL) Beyond prediction and interpretability, the BL framework supports a third fundamental goal: the identification of ground-truth parameters, which in turn endows BL with the capacity for scientifically credible modeling. We refer to this setting as Identifiable Behavior Learning (IBL). In the IBL setting, we define the modular block as ℬid(,;θ):=λ0⊤ϕid(UθU(,))−λ1⊤ρid(θC(,))−λ2⊤ψid(θT(,))B^id(x,y;θ):= _0^\! φ^id (U_ _U(x,y) )- _1 ρ^id (C_ _C(x,y) )- _2 ψ^id (T_ _T(x,y) ) (9) Unlike BL, which uses general nonlinearities, the IBL architecture imposes stricter structural constraints: ϕidφ^id and ρidρ^id are strictly increasing, while ψidψ^id is symmetric and strictly increasing in |⋅||·|. In addition, all three functions are C1C^1. These properties ensure that each UMP block stays responsive and adjusts smoothly to objectives and constraints. In practice, we instantiate equation 9 as ℬid(,)=λ0⊤tanh(u(,))−λ1⊤softplus(c(,))−λ2⊤(t(,))⊙2B^id(x,y)= _0^\! (p_u(x,y) )- _1^\! softplus (p_c(x,y) )- _2^\! (p_t(x,y) ) 2 (10) where (⋅)⊙2(·) 2 denotes elementwise square. We design IBL in three architectural forms. Similar to BL, the IBL(Single) directly uses ℬid(,)B^id(x,y) as the compositional utility function. The IBL(Shallow) and IBL(Deep) variants are defined recursively as IBL(,):=L∘⋅Lid(⋯2id(1id(,))⋯),L≥1IBL(x,y):=W _L·B^id_L (·sB^id_2 (B^id_1(x,y) )·s ), L≥ 1 (11) where ℓidB^id_ stacks multiple parallel blocks ℬℓ,iid(,)B^id_ ,i(x,y), and L∘W _L is a learnable affine transformation without bias. All other design choices follow the BL setting. Theoretical Foundation. IBL admits favorable properties for ground-truth identification. We begin by establishing identifiability, which is fundamental for statistical inference. We first state our key assumption (see Assumption B.1 for details). Assumption 2.1. Let Ψ¯ denote the quotient space of atomic parameters. We assume that the map Ψ¯→ℝ× ^X×Y, ψ¯↦gψ¯ ψ g_ ψ, is injective, and that any finite set of distinct atoms is linearly independent. We further restrict attention to minimal representations with no duplicate atoms and a fixed canonical ordering. Theorem 2.4 (Identifiability of IBL). Under Assumption B.1, the architectures IBL(Single), IBL(Shallow), and IBL(Deep) are identifiable in the parameter quotient space Θ¯ . Theorem 2.5 (Loss Identifiability of IBL). The IBL model is parameterized by θ∈Θθ∈ . Suppose Θ is compact. Then under Assumption B.1, the population loss ℒL defined in equation 6 satisfies: • If γc>0 _c>0, it admits a unique minimizer in the quotient space Θ¯ ; • If γc=0 _c=0, it admits a unique minimizer in the scale-invariant quotient space Θ~ . Theorems 2.4 and 2.5 together establish the identifiability of IBL. Theorem 2.4 shows that if two IBL models of the same structure induce the same compositional utility, then their parameters coincide up to an equivalence class. Theorem 2.5 further extends this result to loss-based identifiability. These results jointly imply that IBL admits a unique parameter estimate up to an equivalence class, and thus yields intrinsic interpretability that is unique up to the same class. Building on identifiability, Theorem 2.6 establishes the statistical consistency of IBL: under compactness of the parameter space, the learned parameters converge in probability to a minimizer of the population loss as the sample size n→∞n→∞. If the model is correctly specified, the estimator further converges to the ground-truth parameter, recovering the true underlying model, thereby endowing IBL with the potential to recover the ground-truth model. Theorem 2.6 (Consistency of IBL). Let Ξ denote the relevant parameter quotient space: Ξ=Θ¯ = if γc>0 _c>0, and Ξ=Θ~ = if γc=0 _c=0. Let θ^n∈argminθ∈Θℳn(θ) θ_n∈ _θ∈ M_n(θ) denote the empirical minimizer, and let θ∙∈argminθ∈Θℳ(θ)θ ∈ _θ∈ M(θ) denote the population minimizer. Then under the conditions of Theorem B.5, θ^n→θ∙in Ξ,ℳ(θ^n)→ℳ(θ∙). θ_n\; p\;θ , ( θ_n)\; p\;M(θ ). Moreover, if the model is correctly specified (i.e., the data distribution is realized by some θ⋆∈Θθ ∈ ), then θ∙=θ⋆θ =θ in Ξ , and thus θ^n→θ⋆ θ_n pθ . Correct specification is a strong and often unrealistic assumption. Fortunately, the IBL framework—like BL—also enjoys a universal approximation guarantee (Theorem B.6). Building on this result, we further establish the universal consistency of IBL: even under misspecification, IBL is capable of recovering the ground-truth model with sufficiently large sample sizes. Theorem 2.7 (Universal Consistency of IBL). Under the conditions of Theorem B.7, for any admissible data-generating distribution p†p satisfying the regularity assumptions of Theorem B.6, the IBL posterior sequence pθ^n\p_ θ_n\ satisfies supx∈KL(p†(⋅∣)∥pθ^n(⋅∣))→0, _x KL\! (p (· )\,\|\,p_ θ_n(· ) ) p0, i.e., the learned conditional distributions pθ^n\p_ θ_n\ converge in KL to p†p uniformly over x. Specifically, this result implies that, even under model misspecification, the learned predictive distribution pθ^np_ θ_n, parameterized by the IBL model, converges uniformly in KL to the true conditional distribution p†p , provided that the capacity of the IBL architecture grows with the sample size n. We also establish the asymptotic normality of IBL estimators (Theorem B.9), showing that the parameter estimates converge in distribution to a normal law as the sample size increases. Furthermore, under additional regularity conditions, the asymptotic variance attains the efficient information bound (Theorem B.10), demonstrating the statistical optimality of IBL. Formal statements and proofs of all theorems in this part are deferred to Appendix B.3. 3 Experiments In this section, we conduct four groups of experiments to systematically evaluate the capabilities of BL. Due to space constraints, details are provided in Appendix C. 3.1 Standard Prediction Tasks Figure 3: Predictive performance of BL and baselines. Left/Middle: relative AUC and F1-Macro gains over DT, sorted by mean (excluding BL). Right: mean F1-Macro ranks (↓ better). BL achieves first-tier performance in both metrics. Its variants rank second and third in mean F1-Macro rank, with BL(Shallow) showing no statistically significant difference from state-of-the-art models. Is BL accurate enough for standard prediction tasks? In this part, we evaluate the predictive performance of BL on 10 datasets (Table 4), covering diverse sample sizes, feature dimensions, and scientific domains. For fair comparison, we consider two BL variants—BL(Single) and BL(Shallow)—and compare them against 10 baseline models (Table 5) drawn from five methodological families: neural networks, tree-based models, gradient boosting methods, Bayesian methods, and linear regressors. All methods share a unified preprocessing and tuning pipeline. Predictive Performance. Figure 3 shows that BL attains first-tier predictive performance overall, achieving the best results among intrinsically interpretable models. Notably, BL(Shallow) surpasses MLP, highlighting that BL delivers interpretability without sacrificing performance. 3.2 Interpreting BL: A Case Study How can BL be interpreted in practice? This part presents a case study using the Boston Housing dataset, where we train a supervised BL(Single) model with a degree-2 polynomial basis, a BL[2,1] model (i.e., a two-layer BL with two B-blocks in the first layer and one in the second layer), and a BL(Deep) model with a [5,3,1] architecture to predict median home values. We illustrate how the internal structure of BL can be interpreted as explicit optimization problems and their hierarchical versions, accompanied by complementary visualizations. Further details are provided in Appendix C.3 and C.5. Symbolic Form of BL(Single) as a UMP. As shown in Figure 2, the trained BL(Single) model can be interpreted as the UMP of a representative buyer in the Boston Housing market, comprising a single objective, inequality, and equality term. Each term is represented by an estimated quadratic polynomial. For parsimony, we extract approximate symbolic expressions by retaining only the monomials with the largest (2–5) absolute coefficients, while collecting the remaining terms (including constants) into a residual term R~ R. For example, the utility term can be written as: u=−0.56⋅P2−0.6⋅RM+0.57⋅RM⋅P+R~u≈(1−P)(1+P−RM)+R~u _u=-0.56· P^2-0.6·RM+0.57·RM· P+ R_u≈(1-P)(1+P-RM)+ R_u We similarly simplify the budget and belief terms to recover an approximate UMP for the buyer. The full symbolic form is illustrated at the bottom of Figure 2. Interpreting BL(Single) via Model Visualization. Visualizations of each term’s polynomial reveal how features constitute the UMP. Three insights emerge from the visualizations in Figure 2. (i) Median housing price (MEDV) and average number of rooms (RM) are dominant across all terms—MEDV negatively affects utility in a near-quadratic form, while RM modulates its marginal effect. (i) Proportion of lower-income residents (LSTAT) features prominently in the budget constraint, reflecting implicit resource limitations. (i) Crime rate (CRIM) appears only in the belief term, suggesting that buyers treat it as influencing others’ behavior rather than their own preferences. Figure 4: Interpreting deeper BL architectures as hierarchical structures of interacting agents. Each block ℬB represents an interpretable agent solving its own UMP, while a layer corresponds to a set of heterogeneous agents operating in parallel. The next layer then aggregates and reallocates the negative energies from the previous layer, thereby performing higher-level coordination across agents. This layered organization provides a natural compositional interpretation of deep BL: bottom-layer modules encode local objectives, while upper layers synthesize these into collective outcomes. Analogous structures arise in biological and social systems—for example, in ant colonies, individual ants (first-layer agents) follow simple local rules, yet their collective behavior is coordinated through higher-level interactions (second-layer aggregation), yielding globally efficient resource allocation and task division. Interpreting BL(Deep). (1) Figure 2 (b) illustrates the optimization problems learned by the BL[2,1] model. Layer 1 identifies two micro-level preference types: an Economic-sensitive Buyer, whose utility and constraint terms load primarily on ZN (Large-lot residential share) and LSTAT (Proportion of lower-income residents); and a Location-sensitive Buyer, driven mainly by CHAS (Charles River indicator) and RAD (Highway accessibility). Layer 2 aggregates these basic preferences, yielding an effective “representative buyer” that integrates the two preference types. (2) Figure 2 (c) presents the internal structure of the BL[5,3,1] model. In Layer 1, BL recovers five distinct micro-level preference types characterizing heterogeneous patterns in the housing market. Layer 2 identifies three macro-level representative agents, each capturing a different macro-level trade-off among the basic preferences. Layer 3 then aggregates these components into a single high-level mechanism, yielding the overall representative buyer. Table 10 provides detailed descriptions of each type. (3) Beyond interpretability, we find that each preference pattern and trade-off recovered by BL(Deep) aligns with established findings in the economics literature (see Table 11). This indicates that BL successfully reconstructs underlying scientific knowledge. 3.3 Prediction on High-Dimensional Inputs Is BL scalable to high-dimensional inputs? We evaluate BL against the energy-based MLP (E-MLP) baseline across network depths d∈1,2,3d∈\1,2,3\, with all models implemented without skip connections. Experiments are conducted on four datasets spanning both image and text domains, and are evaluated using six metrics: in-distribution accuracy, calibration metrics (ECE and NLL), and OOD robustness metrics (AUROC, AUPR, and FPR@95). For OOD evaluation, we adopt symmetric ID↔ splits, using MNIST (LeCun et al., 2002) and Fashion-MNIST (Xiao et al., 2017) as one pair, and AG News and Yelp Polarity (Zhang et al., 2015) as another. E-MLP and BL are controlled to have comparable parameters. Scalability on High-Dimensional Inputs. Figure 5 and Table 1 present results for BL and E-MLP across network depths. On image datasets, the two models exhibit comparable in-distribution accuracy, while BL generally achieves stronger out-of-distribution detection performance on Fashion-MNIST at similar accuracy levels. On text datasets, BL consistently improves ID accuracy over E-MLP across depths. However, OOD detection behavior varies by dataset: BL outperforms E-MLP on Yelp, whereas E-MLP shows better OOD discrimination on AG News. BL also achieves better calibration metrics (ECE and NLL; Table 2). Downward Shift of the Pareto Frontier. Table 13 reports the parameter counts of BL and E-MLP across four tasks, and Tables 3 summarize their runtimes. The two models have highly comparable parameter sizes. Across datasets, BL exhibits slightly higher training time than E-MLP. Combining these results with their comparable predictive performance and the intrinsic interpretability of BL, in contrast with the black-box E-MLP, indicates that BL achieves a downward shift of the Pareto frontier. Figure 5: Comparison of BL and E-MLP on image and text datasets; d denotes model depth. Table 1: ID accuracy and OOD AUROC (%) on image and text datasets. BL and E-MLP are evaluated at depths 1–3 with matched parameter counts, both without skip connections. Top-two per column are blue and red. Image Datasets Model MNIST Fashion-MNIST Accuracy OOD AUROC Accuracy OOD AUROC E-MLP (depth=1) 98.15±\;±\;0.07 88.72±\;±\;1.36 88.79±\;±\;0.29 90.57±\;±\;1.39 BL (depth=1) 97.97±\;±\;0.18 91.17±\;±\;2.68 89.26±\;±\;0.22 91.89±\;±\;0.71 E-MLP (depth=2) 98.11±\;±\;0.08 90.32±\;±\;1.74 88.88±\;±\;0.26 84.61±\;±\;2.56 BL (depth=2) 98.05±\;±\;0.12 90.57±\;±\;2.49 88.96±\;±\;0.39 89.87±\;±\;2.48 E-MLP (depth=3) 98.14±\;±\;0.11 87.76±\;±\;2.55 89.33±\;±\;0.25 83.13±\;±\;1.90 BL (depth=3) 97.93±\;±\;0.27 92.92±\;±\;1.69 88.79±\;±\;0.25 89.24±\;±\;4.18 Text Datasets Model AG News Yelp Accuracy OOD AUROC Accuracy OOD AUROC E-MLP (depth=1) 88.74±\;±\;0.26 59.24±\;±\;0.21 91.16±\;±\;0.02 57.60±\;±\;0.31 BL (depth=1) 89.52±\;±\;0.16 66.18±\;±\;0.20 91.56±\;±\;0.04 57.06±\;±\;0.10 E-MLP (depth=2) 89.29±\;±\;0.20 62.48±\;±\;0.76 91.32±\;±\;0.09 57.47±\;±\;0.21 BL (depth=2) 89.22±\;±\;0.20 63.68±\;±\;0.46 91.39±\;±\;0.06 57.31±\;±\;0.27 E-MLP (depth=3) 89.37±\;±\;0.21 66.82±\;±\;1.01 91.23±\;±\;0.07 57.36±\;±\;0.27 BL (depth=3) 88.80±\;±\;0.18 64.44±\;±\;0.52 91.13±\;±\;0.09 57.16±\;±\;0.48 Table 2: ECE and NLL on image and text datasets. BL and E-MLP are evaluated at depths 1–3 with matched parameter counts. Top-two per column are blue and red. Model MNIST Fashion-MNIST ECE NLL ECE NLL E-MLP (depth=1) 0.02±\;±\;0.00 0.20±\;±\;0.02 0.08±\;±\;0.00 0.74±\;±\;0.01 BL (depth=1) 0.02±\;±\;0.00 0.26±\;±\;0.01 0.05±\;±\;0.00 0.36±\;±\;0.01 E-MLP (depth=2) 0.02±\;±\;0.00 0.23±\;±\;0.02 0.09±\;±\;0.00 0.89±\;±\;0.03 BL (depth=2) 0.02±\;±\;0.00 0.16±\;±\;0.01 0.07±\;±\;0.00 0.44±\;±\;0.01 E-MLP (depth=3) 0.02±\;±\;0.00 0.16±\;±\;0.02 0.09±\;±\;0.00 0.85±\;±\;0.04 BL (depth=3) 0.02±\;±\;0.00 0.13±\;±\;0.02 0.07±\;±\;0.00 0.49±\;±\;0.02 Model AG News Yelp ECE NLL ECE NLL E-MLP (depth=1) 0.02±\;±\;0.00 0.40±\;±\;0.01 0.01±\;±\;0.00 0.24±\;±\;0.00 BL (depth=1) 0.02±\;±\;0.00 0.31±\;±\;0.01 0.00±\;±\;0.00 0.20±\;±\;0.00 E-MLP (depth=2) 0.02±\;±\;0.00 0.42±\;±\;0.01 0.00±\;±\;0.00 0.25±\;±\;0.00 BL (depth=2) 0.06±\;±\;0.01 0.43±\;±\;0.03 0.02±\;±\;0.00 0.23±\;±\;0.01 E-MLP (depth=3) 0.01±\;±\;0.00 0.41±\;±\;0.02 0.00±\;±\;0.00 0.25±\;±\;0.01 BL (depth=3) 0.05±\;±\;0.01 0.39±\;±\;0.02 0.02±\;±\;0.00 0.22±\;±\;0.00 Table 3: Training time (seconds) of BL vs. E-MLP on high-dimensional datasets (mean ± std). Model MNIST FashionMNIST AG News Yelp E-MLP (depth=1) 100.59±\;±\;0.29 73.57±\;±\;1.20 14.69±\;±\;0.40 179.37±\;±\;0.73 BL (depth=1) 110.63±\;±\;3.34 96.52±\;±\;2.90 17.20±\;±\;0.06 181.07±\;±\;1.80 E-MLP (depth=2) 102.64±\;±\;0.26 78.25±\;±\;0.28 15.76±\;±\;0.06 179.22±\;±\;0.66 BL (depth=2) 122.85±\;±\;3.95 114.43±\;±\;3.72 21.78±\;±\;0.08 180.38±\;±\;1.44 E-MLP (depth=3) 104.52±\;±\;0.30 85.57±\;±\;1.19 16.95±\;±\;0.05 178.99±\;±\;1.42 BL (depth=3) 140.17±\;±\;4.42 130.03±\;±\;4.96 26.29±\;±\;0.24 180.36±\;±\;0.91 3.4 Constraint Enforcement Test: High-Dimensional Energy Conservation To evaluate whether the learnable penalty terms in BL are capable of enforcing near-hard constraints under finite temperature, we isolate the penalty mechanism and test it on a high-dimensional energy-conservation constraint. This diagnostic experiment removes the utility term and focuses solely on the penalty term, providing a characterization of how the penalty term controls constraint violations as a function of temperature τ and penalty scale λ. Experiment setup. We sample x∈R64x∈ R^64 i.i.d. from a standard Gaussian x∼(0,I64)x (0,I_64) and define a pure penalty compositional utility T(x,y)=‖y‖2−‖x‖2,BL(x,y)=−λT(x,y)2,T(x,y)\;=\;\|y\|^2-\|x\|^2, (x,y)\;=\;-λ\,T(x,y)^2, which plays the role of an energy-conservation residual and its quadratic penalty. We target the Gibbs distribution p(y∣x)∝exp(BL(x,y)/τ)p(y x)\; \; \! (BL(x,y)/τ ) using overdamped Langevin dynamics with step size η=10−4η=10^-4: yk+1=yk+η∇yBL(x,yk)/τ+2ητξk,ξk∼(0,I64).y_k+1=y_k+η\, _yBL(x,y_k)/τ+ 2ητ\, _k, _k (0,I_64). For each pair (λ,τ)(λ,τ) we run 512512 parallel chains, each for 15001500 Langevin steps (500 burn-in). We sweep over temperatures τ∈2.0,1.0,0.5,0.25,0.1,0.05,0.02,0.01,0.005τ∈\2.0,1.0,0.5,0.25,0.1,0.05,0.02,0.01,0.005\ at a fixed penalty λ=25λ=25, and over penalty weights λ∈0,1,3,10,30,100,200,500λ∈\0,1,3,10,30,100,200,500\ at a fixed temperature τ=0.05τ=0.05. For each configuration we record the residual magnitude |T(x,y)||T(x,y)| from the final state of every chain. We then report three summary statistics: (i) the mean violation E[|T(x,y)|]E[|T(x,y)|], (i) the 9595th percentile of |T(x,y)||T(x,y)|, and (i) the empirical probability of near-feasible samples. We declare a sample to satisfy the constraint approximately if |T(x,y)|≤εtolwithεtol=10−1,|T(x,y)|\;≤\ _tol _tol=10^-1, and estimate P(|T(x,y)|≤εtol)P(|T(x,y)|≤ _tol) across chains. This tolerance scale is chosen to be small relative to the typical unconstrained residuals, so that the near-feasible regime corresponds to a practically tight energy-conservation constraint. Figure 6: Constraint enforcement test of the BL penalty block on an energy-conservation constraint. The figure reports violation statistics |T(x,y)||T(x,y)| when varying the temperature τ (left side of panel) and the penalty weight λ (right side of panel). Constraint enforcement. Figure 6 shows that BL achieves near-hard constraint enforcement under finite temperature and penalty scaling. Violations decrease substantially as τ decreases or λ increases. At around λ=25λ=25 and τ=0.01τ=0.01, the 64-dimensional energy-conservation constraint is enforced within 10−210^-2 error. Curves remain mostly smooth and monotone in 64 dimensions, indicating stable Langevin sampling and effective penalty enforcement. 4 Scientific Explanation of BL(Deep) BL(Deep) provides a form of interpretability that is consistent with hierarchical optimization structures. In BL, each layer performs a coarse-graining of the optimization structure implemented by the layer below. An intuitive analogy is a corporate organizational hierarchy: lower-layer managers solve their own local optimization problems, while higher-layer managers aggregate and coordinate the outcomes of many such lower-layer problems to achieve broader organizational objectives. BL(Deep) follows the same principle—higher layers summarize, reorganize, and coordinate the solutions formed at lower layers. This perspective aligns with many scientific domains characterized by multi-level complexity, including (i) the formation of representative behavioral agents in behavioral sciences, and (i) renormalization in statistical physics, where fine-scale interactions are compressed into effective coarse-scale potentials. We describe the explanation procedure below. To build intuition, let us first consider a generic hierarchical optimization structure—this may refer to a multi-layer organizational structure composed of individual agents, or a multi-scale physical system composed of interacting particles. Step 1: Bottom-layer interpretation. Each bottom-layer block is an optimization problem that directly receives inputs from the environment. These blocks correspond to micro-level behavioral mechanisms, such as the decision rules of individual agents performing environment-facing tasks in an organization, or the motion laws governing a single particle in statistical physics. Examining these bottom-layer blocks reveals the fundamental optimization principles followed by all units that directly interact with the environment. Step 2: Layer-wise coarse-graining and micro-to-macro aggregation. Blocks in the next layer aggregate the outputs of lower-layer optimization problems through a new optimization step, producing a coarse-grained behavioral summary. Each higher-level block represents the effective optimization system that emerges from the interactions among many lower-level units, thereby capturing macro-level regularities distilled from micro-level mechanisms. This micro-to-macro transition is consistent with many well-established scientific principles, including: • (i) Aggregation and coordination: in hierarchical organizations, the outputs of lower-level agents are aggregated, reallocated, and coordinated by higher-level agents to achieve improved organizational objectives. • (i) Coarse-grained observation: in hierarchical behavioral systems, individual agents are grouped into categories that share characteristic optimization patterns; in statistical physics, many particles collectively form systems whose coarse-grained behavior is governed by effective potentials induced by microscopic interactions. Step 3: Bottom-up reconstruction. A global explanation is obtained by tracing the hierarchy upward, following the model’s micro-to-macro abstraction path: raw input features → micro-level optimization blocks → macro-level aggregation and coordination or coarse-grained behavioral constructs → macro-level optimization system. At each layer, we inspect the characteristics of each block and its associated optimization objective, as well as how these optimization problems evolve across layers. This reveals how each higher layer aggregates, coordinates, or coarse-grains the outputs of the layer below. Together, these observations yield a compact multi-scale interpretation in which BL is understood as a hierarchical optimization structure. 5 Discussion In what follows, we discuss the limitations and future directions of Behavior Learning from the perspectives of theoretical foundations, architecture, and applications. Scalability of theoretical assumptions. The identifiability-related statistical theorems constitute the core theoretical pillars of IBL, ensuring uniqueness of the interpretability and supporting its scientific credibility. Although these results hold under mild conditions, their behavior in large-scale, highly over-parameterized architectures remains less well understood. This highlights the need for systematic investigations into the robustness, potential failure modes, and empirical boundaries of these guarantees when applied to modern large-scale learning systems. Choice of basis functions. Polynomial basis functions enhance expressivity while preserving symbolic interpretability in BL (Single). However, high-order polynomials may introduce optimization instability, exacerbate sensitivity to initialization and normalization, and complicate training dynamics. Future work may explore alternative basis families—such as trigonometric, spline-based, or neural basis functions—and develop conditioning or normalization strategies that improve numerical stability without sacrificing interpretability. Interpretable generative modeling. BL integrates several training techniques from energy-based models while retaining intrinsic interpretability, enabling interpretable generative modeling for vision (e.g., image or video generation) and language (e.g., large language models). Extending BL to explicitly generative architectures in which outputs correspond directly to human-understandable and scientifically meaningful blocks represents a compelling direction. Such extensions could yield generative systems with greater transparency, controllability, and scientific credibility compared to traditional black-box models. Hybrid architectures for partial interpretability. A promising direction for future work is to develop hybrid architectures that integrate BL with black-box models in a principled way to achieve partial interpretability. Three avenues are particularly worth exploring: (i) Feature-level integration. Black-box neural networks can serve as high-capacity feature extractors, while BL operates on the resulting learned representations to impose structured, optimization-based semantics. (i) Decision-critical integration. BL blocks may be inserted specifically at high-risk or decision-critical components of the model, substantially reducing the interpretability and reliability risks associated with purely black-box architectures. (i) Mechanism-level integration. Because BL provides an optimization-driven inductive bias aligned with many real-world mechanisms, selectively applying BL to the parts of the system where such inductive bias is essential may yield models that better capture the underlying ground-truth processes while retaining the flexibility of deep networks, thereby improving generalization performance. BL for scientific and social-scientific modeling. BL represents data as a composition of optimization problems, closely resonating with modeling paradigms in the natural and social sciences. Its competitive performance, intrinsic interpretability, and statistical rigor position BL as a promising framework for scientific machine learning. Future research may apply BL to domains such as statistical physics, evolutionary biology, computational neuroscience, and climate dynamics, as well as behavioral science, economics, sociology, and political science—particularly in settings involving complex, partially formalized, or cognitively meaningful structures. 6 Related Work 6.1 Interpretability Interpretability has become increasingly vital in machine learning (Lipton, 2018; Molnar, 2020), especially for scientific domains (Doshi-Velez & Kim, 2017; Roscher et al., 2020). Ensuring interpretability fosters transparency and reproducibility, and may further provide insights into underlying scientific principles. The ideal form of interpretability is intrinsic interpretability, in which a model’s structure or parameters are directly understandable to humans. However, intrinsic interpretability is challenging to achieve in some widely used high-capacity models such as deep neural networks (LeCun et al., 2015). This has motivated post-hoc interpretability methods (Ribeiro et al., 2016; Lundberg & Lee, 2017), which seek to explain a pre-trained black-box model. While more broadly applicable, such explanations are often considered less suitable for scientific research (Rudin, 2019), as they may compromise stability and faithfulness to the model’s decision process. Performance–Interpretability Trade-off. The limited intrinsic interpretability observed in high-capacity models has long been recognized as a central challenge. This is commonly framed as the performance–interpretability trade-off (Rudin, 2019; Arrieta et al., 2020), which posits a tension between predictive performance and intrinsic interpretability. High-performing models such as deep neural networks often lack transparency, whereas intrinsically interpretable models struggle to capture complex nonlinear patterns. Several efforts have sought to mitigate the performance–interpretability trade-off, which can be broadly categorized into four groups. (i) Additive models. Classical GAMs (Hastie, 2017), modern GA2Ms/EBMs (Caruana et al., 2015; Nori et al., 2019), and neural variants such as NAM (Agarwal et al., 2021) and NODE-GAM (Chang et al., 2021) preserve interpretability by decomposing predictions into main effects and low-order interactions. (i) Concept-based models. Concept Bottleneck Models (Koh et al., 2020), TCAV (Kim et al., 2018), and SENN (Alvarez Melis & Jaakkola, 2018) map inputs into human-interpretable latent concepts and use them as intermediate predictors. (i) Rule- and score-based systems. SLIM (Ustun & Rudin, 2016) and CORELS (Angelino et al., 2018) generate transparent scoring functions or rule lists with provable optimality guarantees. (iv) Shape-constrained networks. Deep Lattice Networks (You et al., 2017) and related monotonic architectures impose monotonicity and calibration constraints to encode domain priors while retaining flexibility. Limitations in Scientifically Credible Modeling. The above approaches demonstrate strengths, yet two fundamental limitations restrict their applicability in scientific research. First, most methods are tool-centric modifications of machine learning architectures rather than frameworks grounded in scientific theory (e.g., optimization, dynamical systems, conservation laws). As recent surveys emphasize (Roscher et al., 2020; Karniadakis et al., 2021; Allen et al., 2023; Bereska & Gavves, 2024; Longo et al., 2024; Mersha et al., 2024), genuine scientific insight requires models linked to mechanistic principles, yet many interpretability techniques remain detached from such principles. Second, these approaches are typically non-identifiable (Ran & Hu, 2017; Méloux et al., 2025), meaning that multiple distinct parameterizations can explain the same data. This lack of uniqueness undermines their reliability for recovering ground-truth mechanisms and, in statistical terms, complicates consistency guarantees. As a result, the trained model may fail to converge to the true data-generating process as sample size increases (Newey & McFadden, 1994; Van der Vaart, 2000). Relation to BL. BL also mitigates the performance–interpretability trade-off. Unlike prior methods, it is principle-driven and scientifically grounded, learning interpretable latent optimization structures directly from data. The framework applies broadly to domains where outcomes arise as solutions to (explicit or latent) optimization problems. It is also identifiable: its smooth and monotone variant, Identifiable Behavior Learning (IBL), guarantees identifiability under mild conditions, ensuring the scientific credibility of its explanations and supporting recovery of the ground-truth model under appropriate conditions. 6.2 Data-Driven Inverse Optimization Inverse optimization (IO) (Ahuja & Orlin, 2001; Chan et al., 2025) is a core paradigm for learning latent optimization problems from observed data. Traditional IO aims to construct objectives or constraints that exactly rationalize a small set of deterministic decisions. In contrast, data-driven IO (Keshavarz et al., 2011; Aswani et al., 2018) focuses on statistically recovering the underlying problem from large-scale, noisy observational data. Inverse optimal control (IOC) (Kalman, 1964; Freeman & Kokotovic, 1996) extends this paradigm to dynamic settings, seeking to infer sequential decision processes from expert trajectories. Within machine learning, inverse reinforcement learning (IRL) (Ng et al., 2000; Wulfmeier et al., 2015) and inverse constrained reinforcement learning (ICRL) (Malik et al., 2021; Liu et al., 2024a) are prominent instances of data-driven IOC: Typically, IRL assumes fixed constraints and learns a reward function, whereas ICRL reverses this role. Both require repeatedly solving for (near-)optimal policies and matching with expert demonstrations—incurring high computational cost. In the behavioral sciences, particularly economics, numerous studies can be viewed as instances of the data-driven IO paradigm. Foundational work (McFadden, 1972; Dubin & McFadden, 1984; Hanemann, 1984; Berry et al., 1993) and related studies typically posits theoretically grounded, parametric utility maximization problems (UMPs) and estimates their structural parameters from observed behavior. Relation to BL. The BL framework also falls under the paradigm of data-driven inverse optimization but differs notably from prior related work in both machine learning and behavioral science. Compared with IRL and ICRL, BL does not rely on matching expert-demonstrated policies with the aim of improving task-specific performance. Instead, it is proposed as a general-purpose, scientifically grounded, and intrinsically interpretable framework that operates via low-cost end-to-end training with a hybrid CE–DSM objective. It jointly learns a utility functions and constraints—a direction that has received little attention in IRL and ICRL (Park et al., 2020; Jang et al., 2023; Liu & Zhu, 2024). Meanwhile, in behavioral science, related work typically formulates distinct utility maximization models under varying assumptions for specific decision contexts, and estimate their parameters accordingly. However, to the best of our knowledge, no existing work proposes a structure-free framework for learning UMPs that generalizes across contexts. BL fills this gap with a structure-free, data-driven approach that does not rely on fixed UMP structures. 6.3 Energy-based Models (EBMs) Energy-based models (EBMs) (LeCun et al., 2006) are a prominent data-driven IO scheme, rooted in the principle of energy minimization from statistical physics. They learn an energy function Eθ(x,y)E_θ(x,y) that parameterizes the compatibility between inputs and outputs, inducing a Gibbs distribution pθ(y∣x)∝exp−Eθ(x,y)p_θ(y x) \-E_θ(x,y)\ that favors outcomes corresponding to low-energy solutions. In practice, this energy function is almost always instantiated by high-capacity neural networks, endowing the learned landscape with strong expressive power but also a black-box nature. Training EBMs typically relies on objectives that circumvent the intractable partition function, with classical approaches including contrastive divergence (Hinton, 2002), persistent contrastive divergence (Tieleman, 2008), and noise-contrastive estimation (Gutmann & Hyvärinen, 2010). A particularly influential line of work is score matching (Hyvärinen & Dayan, 2005) and its denoising variant (DSM) (Vincent, 2011), which have underpinned breakthroughs in score-based generative modeling (Song & Ermon, 2019; 2020) and laid the foundation for modern diffusion methods (Song et al., 2020). Relation to BL. BL and EBMs exhibit a principled correspondence: BL is grounded in behavioral science and rooted in utility maximization, while EBMs are grounded in statistical physics and based on energy minimization. BL adopts several training techniques common to EBMs, such as Gibbs distribution modeling and denoising score matching (DSM). However, the two frameworks differ substantially in model structure. EBMs primarily focus on generative quality and typically employ black-box neural networks to learn an opaque energy function with little regard for interpretability. In contrast, BL is built on the utility maximization problem (UMP) and its equivalence to penalty formulations, yielding a principled and scientifically grounded framework. Its architecture is composed of intrinsically interpretable blocks, each of which can be explicitly expressed in symbolic form as a UMP—a foundational paradigm in behavioral science and a universal optimization framework. These properties enable BL to jointly achieve high predictive performance, intrinsic interpretability, and identifiability, thereby supporting scientifically credible modeling that extends beyond mere generative capability. 7 Acknowledgements We would like to thank Prof. Dr. Philipp Hennig, Shu Liu, and Prof. Dr. Sen Geng for their helpful discussions and valuable suggestions. We are also grateful to participants of the Xi’an Jiaotong University seminar for their constructive feedback. We acknowledge the public computational resources provided by the University of Tübingen. Finally, we sincerely thank all anonymous reviewers for their insightful comments. In particular, we appreciate Reviewer sGAR for the highly constructive advice. The energy conservation constraint experiment was added following a suggestion from Reviewer sGAR. If any errors remain, they are solely our responsibility. References Agarwal et al. (2020) Rishabh Agarwal, Nicholas Frosst, Xuezhou Zhang, Rich Caruana, and Geoffrey E Hinton. Neural additive models: Interpretable machine learning with neural nets. arXiv preprint arXiv:2004.13912, 2020. Agarwal et al. (2021) Rishabh Agarwal, Levi Melnick, Nicholas Frosst, Xuezhou Zhang, Ben Lengerich, Rich Caruana, and Geoffrey E Hinton. Neural additive models: Interpretable machine learning with neural nets. Advances in neural information processing systems, 34:4699–4711, 2021. Ahuja & Orlin (2001) Ravindra K Ahuja and James B Orlin. Inverse optimization. Operations research, 49(5):771–783, 2001. Akiba et al. (2019) Takuya Akiba, Shotaro Sano, Toshihiko Yanase, Takeru Ohta, and Masanori Koyama. Optuna: A next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining, p. 2623–2631, 2019. Allen et al. (2023) Genevera I Allen, Luqin Gan, and Lili Zheng. Interpretable machine learning for discovery: Statistical challenges and opportunities. Annual Review of Statistics and Its Application, 11, 2023. Alvarez Melis & Jaakkola (2018) David Alvarez Melis and Tommi Jaakkola. Towards robust interpretability with self-explaining neural networks. Advances in neural information processing systems, 31, 2018. Anderson (1972) Philip W Anderson. More is different: broken symmetry and the nature of the hierarchical structure of science. Science, 177(4047):393–396, 1972. Angelino et al. (2018) Elaine Angelino, Nicholas Larus-Stone, Daniel Alabi, Margo Seltzer, and Cynthia Rudin. Learning certifiably optimal rule lists for categorical data. Journal of Machine Learning Research, 18(234):1–78, 2018. Arik & Pfister (2021) Sercan Ö Arik and Tomas Pfister. Tabnet: Attentive interpretable tabular learning. In Proceedings of the AAAI conference on artificial intelligence, volume 35, p. 6679–6687, 2021. Arrieta et al. (2020) Alejandro Barredo Arrieta, Natalia Díaz-Rodríguez, Javier Del Ser, Adrien Bennetot, Siham Tabik, Alberto Barbado, Salvador García, Sergio Gil-López, Daniel Molina, Richard Benjamins, et al. Explainable artificial intelligence (xai): Concepts, taxonomies, opportunities and challenges toward responsible ai. Information fusion, 58:82–115, 2020. Arthur (2009) W Brian Arthur. Complexity and the economy. In Handbook of Research on Complexity. Edward Elgar Publishing, 2009. Aswani et al. (2018) Anil Aswani, Zuo-Jun Shen, and Auyon Siddiq. Inverse optimization with noisy data. Operations Research, 66(3):870–892, 2018. Balseiro et al. (2019) Santiago R Balseiro, Omar Besbes, and Gabriel Y Weintraub. Dynamic mechanism design with budget-constrained buyers under limited commitment. Operations Research, 67(3):711–730, 2019. Bayer et al. (2007) Patrick Bayer, Fernando Ferreira, and Robert McMillan. A unified framework for measuring preferences for schools and neighborhoods. Journal of political economy, 115(4):588–638, 2007. Bereska & Gavves (2024) Leonard Bereska and Efstratios Gavves. Mechanistic interpretability for ai safety–a review. arXiv preprint arXiv:2404.14082, 2024. Berry et al. (1993) Steven T Berry, James A Levinsohn, and Ariel Pakes. Automobile prices in market equilibrium: Part i and i, 1993. Black (1999) Sandra E Black. Do better schools matter? parental valuation of elementary education. The quarterly journal of economics, 114(2):577–599, 1999. Caruana et al. (2015) Rich Caruana, Yin Lou, Johannes Gehrke, Paul Koch, Marc Sturm, and Noemie Elhadad. Intelligible models for healthcare: Predicting pneumonia risk and hospital 30-day readmission. In Proceedings of the 21th ACM SIGKDD international conference on knowledge discovery and data mining, p. 1721–1730, 2015. Chan et al. (2025) Timothy CY Chan, Rafid Mahmood, and Ian Yihang Zhu. Inverse optimization: Theory and applications. Operations Research, 73(2):1046–1074, 2025. Chang et al. (2021) Chun-Hao Chang, Rich Caruana, and Anna Goldenberg. Node-gam: Neural generalized additive model for interpretable deep learning. arXiv preprint arXiv:2106.01613, 2021. Chay & Greenstone (2005) Kenneth Y Chay and Michael Greenstone. Does air quality matter? evidence from the housing market. Journal of political Economy, 113(2):376–424, 2005. Debreu (1959) Gerard Debreu. Theory of value: An axiomatic analysis of economic equilibrium, volume 17. Yale University Press, 1959. Doshi-Velez & Kim (2017) Finale Doshi-Velez and Been Kim. Towards a rigorous science of interpretable machine learning. arXiv preprint arXiv:1702.08608, 2017. Dubin & McFadden (1984) Jeffrey A Dubin and Daniel L McFadden. An econometric analysis of residential electric appliance holdings and consumption. Econometrica: Journal of the Econometric Society, p. 345–362, 1984. Fisher (1999) Ronald Aylmer Fisher. The genetical theory of natural selection: a complete variorum edition. Oxford University Press, 1999. Freeman & Kokotovic (1996) Randy A Freeman and Petar V Kokotovic. Inverse optimality in robust stabilization. SIAM journal on control and optimization, 34(4):1365–1391, 1996. Gardner et al. (2018) Jacob R Gardner, Geoff Pleiss, David Bindel, Kilian Q Weinberger, and Andrew Gordon Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. In Advances in Neural Information Processing Systems, 2018. Gibbons & Machin (2005) Stephen Gibbons and Stephen Machin. Valuing rail access using transport innovations. Journal of urban Economics, 57(1):148–169, 2005. Gibbs (1902) Josiah Willard Gibbs. Elementary principles in statistical mechanics: developed with especial reference to the rational foundations of thermodynamics. C. Scribner’s sons, 1902. Glaeser & Gyourko (2002) Edward L Glaeser and Joseph Gyourko. The impact of building restrictions on housing affordability. Federal Reserve Bank of New York, Economic Policy Review, 2002:1–19, 2002. Gutmann & Hyvärinen (2010) Michael Gutmann and Aapo Hyvärinen. Noise-contrastive estimation: A new estimation principle for unnormalized statistical models. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, p. 297–304. JMLR Workshop and Conference Proceedings, 2010. Han & Mangasarian (1979) S-P Han and Olvi L Mangasarian. Exact penalty functions in nonlinear programming. Mathematical programming, 17(1):251–269, 1979. Hanemann (1984) W Michael Hanemann. Discrete/continuous models of consumer demand. Econometrica: Journal of the Econometric Society, p. 541–561, 1984. Hastie (2017) Trevor J Hastie. Generalized additive models. Statistical models in S, p. 249–307, 2017. He et al. (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun. Deep residual learning for image recognition. In Proceedings of the IEEE conference on computer vision and pattern recognition, p. 770–778, 2016. Hinton (2002) Geoffrey E Hinton. Training products of experts by minimizing contrastive divergence. Neural computation, 14(8):1771–1800, 2002. Huang et al. (2017) Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger. Densely connected convolutional networks. In Proceedings of the IEEE conference on computer vision and pattern recognition, p. 4700–4708, 2017. Hyvärinen & Dayan (2005) Aapo Hyvärinen and Peter Dayan. Estimation of non-normalized statistical models by score matching. Journal of Machine Learning Research, 6(4), 2005. Jang et al. (2023) Jaehwi Jang, Minjae Song, and Daehyung Park. Inverse constraint learning and generalization by transferable reward decomposition. IEEE Robotics and Automation Letters, 9(1):279–286, 2023. Jevons (2013) William Jevons. The theory of political economy. Springer, 2013. Kadanoff (1966) Leo P Kadanoff. Scaling laws for ising models near t c. Physics Physique Fizika, 2(6):263, 1966. Kalman (1964) Rudolf Emil Kalman. When is a linear control system optimal? 1964. Karniadakis et al. (2021) George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. Physics-informed machine learning. Nature Reviews Physics, 3(6):422–440, 2021. Kayid et al. (2020) Amr Kayid, Nicholas Frosst, and Geoffrey E Hinton. Neural additive models library, 2020. Ke et al. (2017) Guolin Ke, Qi Meng, Thomas Finley, Taifeng Wang, Wei Chen, Weidong Ma, Qiwei Ye, and Tie-Yan Liu. Lightgbm: A highly efficient gradient boosting decision tree. Advances in neural information processing systems, 30, 2017. Keshavarz et al. (2011) Arezou Keshavarz, Yang Wang, and Stephen Boyd. Imputing a convex objective function. In 2011 IEEE international symposium on intelligent control, p. 613–619. IEEE, 2011. Kim et al. (2018) Been Kim, Martin Wattenberg, Justin Gilmer, Carrie Cai, James Wexler, Fernanda Viegas, et al. Interpretability beyond feature attribution: Quantitative testing with concept activation vectors (tcav). In International conference on machine learning, p. 2668–2677. PMLR, 2018. Kingma (2014) Diederik P Kingma. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980, 2014. Koh et al. (2020) Pang Wei Koh, Thao Nguyen, Yew Siang Tang, Stephen Mussmann, Emma Pierson, Been Kim, and Percy Liang. Concept bottleneck models. In International conference on machine learning, p. 5338–5348. PMLR, 2020. Kraus et al. (2024) Mathias Kraus, Daniel Tschernutter, Sven Weinzierl, and Patrick Zschech. Interpretable generalized additive neural networks. European Journal of Operational Research, 317(2):303–316, 2024. Kullback & Leibler (1951) Solomon Kullback and Richard A Leibler. On information and sufficiency. The annals of mathematical statistics, 22(1):79–86, 1951. Landau & Lifshitz (2013) Lev Davidovich Landau and Evgenii Mikhailovich Lifshitz. Statistical Physics: Volume 5, volume 5. Elsevier, 2013. LeCun et al. (2002) Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278–2324, 2002. LeCun et al. (2006) Yann LeCun, Sumit Chopra, Raia Hadsell, M Ranzato, Fujie Huang, et al. A tutorial on energy-based learning. Predicting structured data, 1(0), 2006. LeCun et al. (2015) Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. nature, 521(7553):436–444, 2015. Lipton (2018) Zachary C Lipton. The mythos of model interpretability: In machine learning, the concept of interpretability is both important and slippery. Queue, 16(3):31–57, 2018. Liu et al. (2024a) Guiliang Liu, Sheng Xu, Shicheng Liu, Ashish Gaurav, Sriram Ganapathi Subramanian, and Pascal Poupart. A comprehensive survey on inverse constrained reinforcement learning: Definitions, progress and challenges. arXiv preprint arXiv:2409.07569, 2024a. Liu & Zhu (2024) Shicheng Liu and Minghui Zhu. Meta inverse constrained reinforcement learning: Convergence guarantee and generalization analysis. International Conference on Learning Representations, 2024. Liu et al. (2024b) Ziming Liu, Yixuan Wang, Sachin Vaidya, Fabian Ruehle, James Halverson, Marin Soljačić, Thomas Y Hou, and Max Tegmark. Kan: Kolmogorov-arnold networks. arXiv preprint arXiv:2404.19756, 2024b. Ljungqvist & Sargent (2018) Lars Ljungqvist and Thomas J Sargent. Recursive macroeconomic theory. MIT press, 2018. Longo et al. (2024) Luca Longo, Mario Brcic, Federico Cabitza, Jaesik Choi, Roberto Confalonieri, Javier Del Ser, Riccardo Guidotti, Yoichi Hayashi, Francisco Herrera, Andreas Holzinger, et al. Explainable artificial intelligence (xai) 2.0: A manifesto of open challenges and interdisciplinary research directions. Information Fusion, 106:102301, 2024. Loshchilov & Hutter (2017) Ilya Loshchilov and Frank Hutter. Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101, 2017. Lundberg & Lee (2017) Scott M Lundberg and Su-In Lee. A unified approach to interpreting model predictions. Advances in neural information processing systems, 30, 2017. Malik et al. (2021) Shehryar Malik, Usman Anwar, Alireza Aghasi, and Ali Ahmed. Inverse constrained reinforcement learning. In International conference on machine learning, p. 7390–7399. PMLR, 2021. Mas-Colell et al. (1995) Andreu Mas-Colell, Michael Dennis Whinston, Jerry R Green, et al. Microeconomic theory, volume 1. Oxford university press New York, 1995. McFadden (1972) Daniel McFadden. Conditional logit analysis of qualitative choice behavior. 1972. McFadden (1977) Daniel McFadden. Modelling the choice of residential location. 1977. Méloux et al. (2025) Maxime Méloux, Silviu Maniu, François Portet, and Maxime Peyrard. Everything, everywhere, all at once: Is mechanistic interpretability identifiable? arXiv preprint arXiv:2502.20914, 2025. Mersha et al. (2024) Melkamu Mersha, Khang Lam, Joseph Wood, Ali K Alshami, and Jugal Kalita. Explainable artificial intelligence: A survey of needs, techniques, applications, and future direction. Neurocomputing, 599:128111, 2024. Mitchell (2009) Melanie Mitchell. Complexity: A guided tour. Oxford university press, 2009. Molnar (2020) Christoph Molnar. Interpretable machine learning. Lulu. com, 2020. Newey & McFadden (1994) Whitney K Newey and Daniel McFadden. Large sample estimation and hypothesis testing. Handbook of econometrics, 4:2111–2245, 1994. Ng et al. (2000) Andrew Y Ng, Stuart Russell, et al. Algorithms for inverse reinforcement learning. In Icml, volume 1, p. 2, 2000. Nori et al. (2019) Harsha Nori, Samuel Jenkins, Paul Koch, and Rich Caruana. Interpretml: A unified framework for machine learning interpretability. arXiv preprint arXiv:1909.09223, 2019. Park et al. (2020) Daehyung Park, Michael Noseworthy, Rohan Paul, Subhro Roy, and Nicholas Roy. Inferring task goals and constraints using bayesian nonparametric inverse reinforcement learning. In Conference on robot learning, p. 1005–1014. PMLR, 2020. Plonsky et al. (2025) Ori Plonsky, Reut Apel, Eyal Ert, Moshe Tennenholtz, David Bourgin, Joshua C Peterson, Daniel Reichman, Thomas L Griffiths, Stuart J Russell, Even C Carter, et al. Predicting human decisions with behavioural theories and machine learning. Nature Human Behaviour, p. 1–14, 2025. Popper (2005) Karl Popper. The logic of scientific discovery. Routledge, 2005. Ramsey (1928) Frank Plumpton Ramsey. A mathematical theory of saving. The economic journal, 38(152):543–559, 1928. Ran & Hu (2017) Zhi-Yong Ran and Bao-Gang Hu. Parameter identifiability in statistical machine learning: a review. Neural Computation, 29(5):1151–1203, 2017. Ribeiro et al. (2016) Marco Tulio Ribeiro, Sameer Singh, and Carlos Guestrin. ” why should i trust you?” explaining the predictions of any classifier. In Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, p. 1135–1144, 2016. Roscher et al. (2020) Ribana Roscher, Bastian Bohn, Marco F Duarte, and Jochen Garcke. Explainable machine learning for scientific insights and discoveries. Ieee Access, 8:42200–42216, 2020. Rosen (1974) Sherwin Rosen. Hedonic prices and implicit markets: product differentiation in pure competition. Journal of political economy, 82(1):34–55, 1974. Rudin (2019) Cynthia Rudin. Stop explaining black box machine learning models for high stakes decisions and use interpretable models instead. Nature machine intelligence, 1(5):206–215, 2019. Samuelson (1948) Paul Anthony Samuelson. Foundations of economic analysis. Science and Society, 13(1), 1948. Simon (1955) Herbert A Simon. A behavioral model of rational choice. The quarterly journal of economics, p. 99–118, 1955. Song & Ermon (2019) Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in neural information processing systems, 32, 2019. Song & Ermon (2020) Yang Song and Stefano Ermon. Improved techniques for training score-based generative models. Advances in neural information processing systems, 33:12438–12448, 2020. Song et al. (2020) Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456, 2020. Tieleman (2008) Tijmen Tieleman. Training restricted boltzmann machines using approximations to the likelihood gradient. In Proceedings of the 25th international conference on Machine learning, p. 1064–1071, 2008. Ustun & Rudin (2016) Berk Ustun and Cynthia Rudin. Supersparse linear integer models for optimized medical scoring systems. Machine Learning, 102(3):349–391, 2016. Van der Vaart (2000) Aad W Van der Vaart. Asymptotic statistics, volume 3. Cambridge university press, 2000. Vincent (2011) Pascal Vincent. A connection between score matching and denoising autoencoders. Neural computation, 23(7):1661–1674, 2011. Wright et al. (1932) Sewall Wright et al. The roles of mutation, inbreeding, crossbreeding, and selection in evolution. 1932. Wulfmeier et al. (2015) Markus Wulfmeier, Peter Ondruska, and Ingmar Posner. Maximum entropy deep inverse reinforcement learning. arXiv preprint arXiv:1507.04888, 2015. Xiao et al. (2017) Han Xiao, Kashif Rasul, and Roland Vollgraf. Fashion-mnist: a novel image dataset for benchmarking machine learning algorithms. arXiv preprint arXiv:1708.07747, 2017. Yang & Schoenholz (2017) Ge Yang and Samuel Schoenholz. Mean field residual networks: On the edge of chaos. Advances in neural information processing systems, 30, 2017. You et al. (2017) Seungil You, David Ding, Kevin Canini, Jan Pfeifer, and Maya Gupta. Deep lattice networks and partial monotonic functions. Advances in neural information processing systems, 30, 2017. Zhang et al. (2015) Xiang Zhang, Junbo Zhao, and Yann LeCun. Character-level convolutional networks for text classification. Advances in neural information processing systems, 28, 2015. Appendix A Architecture Details A.1 Learning Scheme Details Input and output of the BL function. We formulate BL as a direct mapping from input–output pairs to compositional utility representations: BL:×→ℝdout,(x,y)↦BL(x,y)∈ℝdout,BL:\ X×Y ^d_out, (x,y) (x,y) ^d_out, where the output dimension doutd_out is chosen according to the modeling choice. This formulation intentionally allows BL to return either a scalar or a vector for each (x,y)(x,y); the following cases are most common: • Scalar per candidate (pointwise evaluation). Set dout=1d_out=1. Here BL(x,y)∈ℝBL(x,y) is a scalar compositional utility evaluated for the single candidate y. This view is natural for continuous y (regression or density estimation) or when one prefers to evaluate candidates individually. • Vectorized over a finite candidate set. If =y1,…,ymY=\y_1,…,y_m\ is finite, one can choose dout=md_out=m and define the vector-valued output by stacking evaluations over the candidate set: BL(x):=[BL(x,y1)⋮BL(x,ym)]∈ℝm.BL(x):= bmatrixBL(x,y_1)\\[4.0pt] \\[4.0pt] BL(x,y_m) bmatrix ^m. This vectorized form is convenient for classification: it evaluates all class candidates at once and yields a single compositional utility vector per x. • Flexibility and equivalence. The scalar and vector modes are compatible: the vectorized form is simply a batch of pointwise evaluations. Conversely, a scalar pointwise evaluator can be used to assemble a vector by repeated calls over a candidate set. The choice between pointwise (scalar) and vectorized outputs is therefore an engineering choice that trades off computational efficiency and convenience. Given a dataset =(xi,yi)i=1nD=\(x_i,y_i)\_i=1^n, training and inference may use either mode: vectorized computation where feasible (e.g., small finite Y), or pointwise evaluation when Y is large or continuous. Conditional Gibbs model. Let (x,y)∼(x,y) with x∈ℝdx ^d and y=(ydisc,ycont)∈disc×ℝmcy=(y^disc,y^cont) _disc×R^m_c (discrete, continuous, or hybrid). BL induces a conditional Gibbs distribution with temperature τ>0τ>0: pτ(y∣x)=expBL(x,y)/τZτ(x),Zτ(x)=∫expBL(x,y′)/τy′.p_τ(y x)= \BL(x,y)/τ\Z_τ(x), Z_τ(x)=\! _Y\! \BL(x,y )/τ\\,dy . For discrete =y1,…,ymY=\y_1,…,y_m\, if we choose the vector-output formulation, we define BL(x):=[BL(x,y1),…,BL(x,ym)]∈ℝm,BL(x):= [BL(x,y_1),…,BL(x,y_m) ] ^m, so that the conditional distribution reduces to a softmax over this compositional utility vector: pτ(y=k∣x)=softmaxk(1τBL(x)).p_τ(y=k x)=softmax_k\! ( 1τ\,BL(x) ). Behaviorally, τ encodes noisy rationality; as τ→0τ\!→ 0, pτ(⋅∣x)p_τ(· x) concentrates on argmaxyBL(x,y) _yBL(x,y), corresponding to the deterministic optimal choice implied by the learned model. Supervised, unsupervised, and generative uses. BL accommodates multiple regimes. (i) Supervised: take x as input and y as label. For discrete y, one may either (a) adopt the vector-output formulation, where BL(x)∈ℝmBL(x) ^m yields a compositional utility vector over all classes and the likelihood is given by a softmax, or (b) adopt the scalar-output formulation, where BL(x,y)BL(x,y) is evaluated separately for each candidate and then normalized across classes. For continuous y, BL naturally operates in the scalar-output mode, treating BL(x,y)∈ℝBL(x,y) as a compositional utility field. (i) Unsupervised / generative: model a marginal p(y)∝expBL(y)/τp(y) \BL(y)/τ\ (empty x) or a joint p(x,y)∝expBL(x,y)/τp(x,y) \BL(x,y)/τ\; sampling the Gibbs distribution yields a generator. Learning objective. Since the response y may contain both discrete and continuous components, we estimate θ by minimizing a type-specific risk: ℒ(θ)=γd[−logpτ(ydisc∣x)]+γc∥∇y~contlogpτ(y~cont∣x)+σ−2(y~cont−ycont)∥2,L(θ)= _d\,E [- p_τ(y^disc x) ]\;+\; _c\,E \| _ y^cont p_τ( y^cont x)+σ^-2( y^cont-y^cont) \|^2, where the first term is cross-entropy on the discrete component and the second is denoising score matching (DSM) on the continuous component with y~cont=ycont+ε y^cont=y^cont+ , ε∼(0,σ2I) (0,σ^2I). Set (γd,γc)=(1,0)( _d, _c)=(1,0) for purely discrete outputs, (0,1)(0,1) for purely continuous outputs, and (>0,>0)(>0,>0) for hybrids. A.2 Model Structure Details In the main text we adopted a compact notation for BL; here we present an equivalent, more explicit matrix/vector formulation that makes dimensions, linear maps, and the per-head parameterizations explicit, which is useful for formal proofs and for implementation details. Fixed bases and head pre-activations. For a block input z (specified below), let mu(z)∈ℝdu,mc(z)∈ℝdc,mt(z)∈ℝdtm_u(z) ^d_u, m_c(z) ^d_c, m_t(z) ^d_t denote fixed basis (e.g., monomial) vectors. Learnable linear maps produce head pre-activations: u(z):=Mumu(z)+bu∈ℝru,c(z):=Mcmc(z)+bc∈ℝrc,t(z):=Mtmt(z)+bt∈ℝrt,u(z):=M_u\,m_u(z)+b_u ^r_u, c(z):=M_c\,m_c(z)+b_c ^r_c, t(z):=M_t\,m_t(z)+b_t ^r_t, with Mu∈ℝru×duM_u ^r_u× d_u, Mc∈ℝrc×dcM_c ^r_c× d_c, Mt∈ℝrt×dtM_t ^r_t× d_t and optional biases b∙b_ . Single BL block. A single modular block is ℬ(z)=λ0⊤ϕ(u(z))−λ1⊤ρ(c(z))−λ2⊤ψ(t(z)),B(z)\;=\; _0 \,φ\! (u(z) )\;-\; _1 \,ρ\! (c(z) )\;-\; _2 \,ψ\! (t(z) ), (12) where λ0∈ℝru _0 ^r_u, λ1∈ℝrc _1 ^r_c, λ2∈ℝrt _2 ^r_t are learnable weights, and ϕ,ρ,ψφ,ρ,ψ act coordinatewise with the roles specified in Theorem 2.1 (increasing ϕφ for utility, penalty ρ for inequality violations, symmetric ψ for equalities). Identifying UθU(x,y)=u(z=(x,y)),θC(x,y)=c(z=(x,y)),θT(x,y)=t(z=(x,y)),U_ _U(x,y)=u (z=(x,y) ), _ _C(x,y)=c (z=(x,y) ), _ _T(x,y)=t (z=(x,y) ), substituting into equation 12 recovers the main-text parameterization in equation 4. Layer of parallel blocks. A layer ℓB_ stacks dℓd_ parallel copies of equation 12 with (possibly) distinct parameters θℓ,i _ ,i: ℓ(zℓ):=[ℬθℓ,1(zℓ)⋮ℬθℓ,dℓ(zℓ)]∈ℝdℓ.B_ (z_ )\;:=\; bmatrixB_ _ ,1(z_ )\\[-1.0pt] \\[-1.0pt] B_ _ ,d_ (z_ ) bmatrix ^d_ . We adopt the standard layered (feedforward) form: z1:=(x,y),zℓ+1:=ℓ(zℓ)(ℓ=1,…,L−1),z_1:=(x,y), z_ +1:=B_ (z_ ) ( =1,…,L-1), so that each layer’s input is simply the previous layer’s output. This is the canonical feedforward architecture. Optionally, one may allow each layer to explicitly access the original inputs: z1:=(x,y),zℓ+1:=ℓ((x,y),zℓ).z_1:=(x,y), z_ +1:=B_ ((x,y),\,z_ ). To improve trainability one may also use residual connections: zℓ+1:=zℓ+ℓ(zℓ).z_ +1:=z_ +B_ (z_ ). Shallow/Deep composition and final affine readout. For depth L≥1L≥ 1, the BL compositional utility is produced by a final learnable affine transformation of the top layer: BL(x,y)=WLL(zL)+bL,BL(x,y)\;=\;W_L\,B_L (z_L )+b_L, (13) with WL∈ℝ1×dLW_L ^1× d_L for scalar output or WL∈ℝm×dLW_L ^m× d_L for vector output, and bias bLb_L of matching dimension. The cases L=1L=1 (with d1=1d_1=1), L≤2L≤ 2, and L>2L>2 correspond to BL(Single), BL(Shallow), and BL(Deep), respectively, exactly as described in the main text. A.3 Implementation Details A.3.1 Function Instantiation Default instantiation. In practice, we instantiate equation 4 with the specific choice (ϕ,ρ,ψ)=(tanh,ReLU,|⋅|)(φ,ρ,ψ)=( ,ReLU,|·|): ℬ(x,y;θ)=λ0⊤tanh(UθU(x,y))−λ1⊤ReLU(θC(x,y))−λ2⊤|θT(x,y)|.B(x,y;θ)= _0^\! \! (U_ _U(x,y) )- _1^\! ReLU\! (C_ _C(x,y) )- _2^\! |T_ _T(x,y) |. (14) Here λ0,λ1,λ2 _0, _1, _2 are learnable nonnegative weights. The bounded tanh captures saturation effects and diminishing returns in the utility head (Jevons, 2013), while ReLUReLU and |⋅||·| impose asymmetric (one-sided) and symmetric (two-sided) penalties for inequality and equality violations. Variants and simplifications. Several variants of equation 14 are often useful: • Identity utility head. Set ϕ=idφ=id so the utility head uses raw polynomials: ℬ=λ0⊤UθU−λ1⊤ReLU(θC)−λ2⊤|θT|.B= _0^\! U_ _U- _1^\! ReLU(C_ _C)- _2^\! |T_ _T|. • Smooth penalty alternatives. Replace ReLUReLU with softplussoftplus to yield smooth inequality penalties, or replace |⋅||·| with Huber or squared penalties to modulate sensitivity near zero for equality terms. • Dropping heads. The framework is modular, so one may omit heads depending on the task: – No T head: ignores symmetric deviations, yielding a constrained maximization with only inequality penalties. – No C head: if the T head is retained, the model reduces to a maximization problem with only equality constraints; if T is also removed, it becomes a fully unconstrained maximization. – No U head: produces a pure (soft-)constraint model focusing on feasibility. Strikingly, removing both U and T leaves only piecewise-linear ReLUReLU penalties; when followed by a final affine readout, the resulting architecture becomes highly similar to a standard MLP—suggesting that MLPs may be viewed as a closely related special instance within the broader BL framework. A.3.2 Polynomial Feature Maps and Linear Reductions We adopt a pragmatic default: use low-degree polynomial maps for single-block models to maximize interpretability, and use affine (degree-1) maps inside blocks for shallow/deep stacks to control parameter growth and compute. Below we state the instantiations and give the final block formulas used in experiments. BL(Single) — polynomial instantiation. Let mD(x,y)m_D(x,y) denote a fixed basis of monomials up to total degree D (e.g. D≤2D≤ 2): mD(x,y)=[x,y,vec(xx⊤),vec(xy⊤),vec(yy⊤),…]⊤.m_D(x,y)= [x,\;y,\;vec(x ),\;vec(xy ),\;vec(y ),… ] . Parameterize each map as a linear map on this basis: UθU(x,y)=MUmD(x,y)+bU,θC(x,y)=MCmD(x,y)+bC,θT(x,y)=MTmD(x,y)+bT,U_ _U(x,y)=M_U\,m_D(x,y)+b_U, _ _C(x,y)=M_C\,m_D(x,y)+b_C, _ _T(x,y)=M_T\,m_D(x,y)+b_T, with learnable matrices M∙M_ and biases b∙b_ . The block becomes ℬ(x,y;θ)=λ0⊤ϕ(MUmD+bU)−λ1⊤ρ(MCmD+bC)−λ2⊤ψ(MTmD+bT).B(x,y;θ)= _0 φ(M_Um_D+b_U)- _1 ρ(M_Cm_D+b_C)- _2 ψ(M_Tm_D+b_T). BL (Shallow/Deep) — linear-by-layer instantiation. For stacked architectures (Shallow/Deep) we use affine maps inside each block to keep per-layer complexity low: UθU(x,y)=AU[x;y]+bU,θC(x,y)=AC[x;y]+bC,θT(x,y)=AT[x;y]+bT,U_ _U(x,y)=A_U\,[x;y]+b_U, _ _C(x,y)=A_C\,[x;y]+b_C, _ _T(x,y)=A_T\,[x;y]+b_T, with learnable A∙A_ and b∙b_ . The corresponding block is ℬ(x,y;θ)=λ0⊤ϕ(AU[x;y]+bU)−λ1⊤ρ(AC[x;y]+bC)−λ2⊤ψ(AT[x;y]+bT).B(x,y;θ)= _0 φ(A_U[x;y]+b_U)- _1 ρ(A_C[x;y]+b_C)- _2 ψ(A_T[x;y]+b_T). On-demand higher-order terms. If diagnostics or domain knowledge indicate underfitting, we optionally augment the affine maps with selected higher-order terms or interactions. Concretely, this is done by appending a small set of monomials (e.g. xiyjx_iy_j, xi2x_i^2, yk2y_k^2) to the input vector [x;y][x;y] and re-estimating the same affine maps A∙A_ . This targeted augmentation preserves the base affine parameterization, increases expressivity only where required, and keeps both computational and statistical costs modest while retaining interpretability. Figure 7: Visualization of polynomial feature maps as computation graphs, where nodes represent variables or outputs and edges represent their effects. The left panel illustrates the linear form ℱ=ax+bF=ax+b, in which the single edge x→ℱx\!→\!F directly encodes the marginal effect of x on ℱF. The middle panel shows the quadratic form ℱ=ax2+bx+cF=ax^2+bx+c, where x not only has a direct edge x→ℱx\!→\!F but also acts on its own edge (“x→ℱx\!→\!F”), thereby modifying the strength of its self-effect through a higher-order contribution. The right panel depicts the interaction form ℱ=ax+by+cxy+dF=ax+by+cxy+d, where y has an edge y→ℱy\!→\!F and, in addition, x acts on this edge (“y→ℱy\!→\!F”), thereby modulating the strength of y’s contribution to ℱF. Symmetrically, y may act on the edge (“x→ℱx\!→\!F”), so that each variable can reshape the other’s effect through the interaction term. A.3.3 Skip Connections Skip connections are optional in our implementation. When beneficial, we often consider two patterns tailored to BL: a DenseNet-style (concatenative) variant and a ResNet-style (additive) variant. Dense skip connections (DenseNet-style, concatenation). This variant feeds each layer with the concatenation of all preceding representations, mirroring DenseNet (Huang et al., 2017). Let z1:=[x;y],s1:=1(z1)∈ℝd1.z_1:=[\,x;\,y\,], s_1:=B_1(z_1) ^d_1. For ℓ≥2 ≥ 2, zℓ:=[x;y;s1;…;sℓ−1],sℓ:=ℓ(zℓ)∈ℝdℓ.z_ :=[\,x;\,y;\,s_1;…;s_ -1\,], s_ :=B_ (z_ ) ^d_ . The final compositional utility is read out as BL(x,y)=WLsL+bL.BL(x,y)\;=\;W_L\,s_L+b_L. Pros. By exposing all earlier block outputs explicitly as inputs to later blocks, dense skips preserve a transparent feature trail: one can trace which intermediate ℬB-block outputs enter downstream computations and the final affine readout. This often improves feature reuse and yields favorable interpretability at the block level. Residual skip connections (ResNet-style, addition). This variant adds an identity (or projected) shortcut to each layer, as in ResNet (He et al., 2016). Define z1:=[x;y],s1:=1(z1)∈ℝd1,z_1:=[\,x;\,y\,], s_1:=B_1(z_1) ^d_1, and for ℓ≥2 ≥ 2, sℓ:=ℓ(sℓ−1)+Πℓsℓ−1,Πℓ∈ℝdℓ×dℓ−1,s_ :=B_ (s_ -1)\;+\; _ \,s_ -1, _ ^d_ × d_ -1, where Πℓ _ is the identity if dℓ=dℓ−1d_ =d_ -1, or a bias-free learnable projection otherwise. The readout is again BL(x,y)=WLsL+bL.BL(x,y)\;=\;W_L\,s_L+b_L. Skip Connections and Interpretability. Skip connections introduce explicit cross-layer dependency structures, a form widely studied in statistical physics and other scientific domains. Such structures enhance scientific interpretability by making long-range influences transparent. In behavioral and organizational sciences, they capture situations in which lower-level agents directly affect higher-level decision makers without routing through intermediate layers. In physics, microscopic parameters can exert direct effects on macroscopic behaviors across multiple scales. Architecturally, ResNet-style skip connections model linear cross-layer dependencies, whereas DenseNet-style connections realize concatenative (information-replicating) dependencies. These mechanisms provide flexible yet interpretable pathways for representing hierarchical interactions. Appendix B Proofs of Theorems B.1 Utility Maximization Problem (UMP) Theorem 2.1 (Local Exact Penalty Reformulation for UMP). Let ⊂ℝdxX ^d_x and ⊂ℝdyY ^d_y be nonempty compact sets, and let U:×→ℝU:X×Y , :×→ℝmC:X×Y ^m, and :×→ℝpT:X×Y ^p be C1C^1. Consider the Utility Maximization Problem (UMP) max∈U(,)s.t.(,)≤0,(,)=0. _y \,U(x,y) .t. (x,y)≤ 0,\;\;T(x,y)=0. (15) Assume there exists a feasible point ⋆∈int()y (Y) which is a strict local maximizer of equation 15 and the Han–Mangasarian constraint qualification (2.1) holds at ⋆y (in the notation of Han & Mangasarian (1979)). Let ϕ:ℝ→ℝφ:R be strictly increasing and C1C^1, and define ρ(z):=maxz,0ρ(z):= \z,0\ and ψ(z):=|z|ψ(z):=|z| (componentwise on ℝmR^m and ℝpR^p). Then there exist λ0>0 _0>0, λ1∈ℝ++m _1 ^m_++, and λ2∈ℝ++p _2 ^p_++ such that ⋆y is a local maximizer of max∈λ0ϕ(U(,))−λ1⊤ρ((,))−λ2⊤ψ((,)). _y \; _0\,φ\! (U(x,y) )- _1 ρ\! (C(x,y) )- _2 ψ\! (T(x,y) ). (16) Proof. Fix ∈x and abbreviate g():=(,)∈ℝm,h():=(,)∈ℝp.g(y):=C(x,y) ^m, h(y):=T(x,y) ^p. Feasibility of ⋆y means g(⋆)≤0g(y )≤ 0 componentwise and h(⋆)=0h(y )=0. Step 1: Convert to a constrained local minimization problem in the ambient space. Pick any λ0>0 _0>0 and define f():=−λ0ϕ(U(,)).f(y):=- _0\,φ\! (U(x,y) ). (17) Since ϕφ is strictly increasing, for any 1,2y_1,y_2 we have U(,1)>U(,2)U(x,y_1)>U(x,y_2) if and only if f(1)<f(2)f(y_1)<f(y_2). Hence ⋆y is a strict local maximizer of equation 15 if and only if ⋆y is a strict local minimizer of min∈f()s.t.g()≤0,h()=0. _y \,f(y) .t. g(y)≤ 0,\ \ h(y)=0. (18) Now use the interior-point assumption ⋆∈int()y (Y): there exists ε0>0 _0>0 such that Bε0(⋆)⊂B_ _0(y ) . Therefore, the notion of strict local minimizer over Y at ⋆y coincides with the ambient-space notion: for any function F, there exists ε∈(0,ε0] ∈(0, _0] such that F(⋆)<F()∀∈(∩Bε(⋆))∖⋆F(y )<F(y) ∀\,y∈ (Y∩ B_ (y ) ) \y \ if and only if F(⋆)<F()∀∈Bε(⋆)∖⋆.F(y )<F(y) ∀\,y∈ B_ (y ) \y \. Hence ⋆y is a strict local minimizer of the constrained problem equation 18 in the ambient-space sense. Moreover, by assumption, the triple (f,g,h)(f,g,h) is continuously differentiable on a neighborhood of ⋆y . Step 2: Embed vector weights into a norm and build a Han–Mangasarian penalty. Define the positive part g+()∈ℝmg_+(y) ^m componentwise by (g+())i:=maxgi(),0(g_+(y))_i:= \g_i(y),0\. Let λ1∈ℝ++m _1 ^m_++ and λ2∈ℝ++p _2 ^p_++ be arbitrary for the moment, and define a norm on ℝm+pR^m+p by the weighted ℓ1 _1-norm ‖(u,v)‖λ:=λ1⊤|u|+λ2⊤|v|,(u,v)∈ℝm×ℝp,\|(u,v)\|_λ:= _1 |u|+ _2 |v|, (u,v) ^m×R^p, (19) where |⋅||·| is componentwise absolute value. Since λ1,λ2 _1, _2 have strictly positive entries, ∥⋅∥λ\|·\|_λ is indeed a norm. Choose the scalar penalty function Q:[0,∞)→[0,∞)Q:[0,∞)→[0,∞) as Q(t)=tQ(t)=t. Then Q satisfies the penalty regularity condition (1.3) in Han & Mangasarian (1979), in particular Q′(0+)=1>0Q (0+)=1>0. Define for α≥0α≥ 0 the penalty function P(,α):=f()+αQ(‖(g+(),h())‖λ).P(y,α):=f(y)+α\,Q\! ( \| (g_+(y),h(y) ) \|_λ ). (20) Expanding equation 20 using equation 19 and Q(t)=tQ(t)=t yields P(,α)=−λ0ϕ(U(,))+α[λ1⊤g+()+λ2⊤|h()|].P(y,α)=- _0\,φ\! (U(x,y) )+α [ _1 g_+(y)+ _2 |h(y)| ]. (21) Since ρ(g())=g+()ρ(g(y))=g_+(y) and ψ(h())=|h()|ψ(h(y))=|h(y)| componentwise, equation 21 can be rewritten as P(,α)=−λ0ϕ(U(,))+α[λ1⊤ρ(g())+λ2⊤ψ(h())].P(y,α)=- _0\,φ\! (U(x,y) )+α [ _1 ρ\! (g(y) )+ _2 ψ\! (h(y) ) ]. (22) Step 3: Apply Han–Mangasarian Theorem 4.4. By Steps 1–2, the functions f,g,hf,g,h are C1C^1 on a neighborhood of ⋆y , ⋆y is a strict local minimizer (ambient-space sense) of the constrained problem equation 18, and the Han–Mangasarian constraint qualification (2.1) holds at ⋆y . Therefore, by (Han & Mangasarian, 1979, Thm. 4.4), there exists α¯≥0 α≥ 0 such that for every α≥α¯α≥ α, ⋆y is a local minimizer of P(⋅,α)P(·,α). Step 4: Return to local maximization and ensure strictly positive vector weights. Choose any α>maxα¯,0,α> \ α,0\, (23) so in particular α>0α>0. Define the penalized maximization objective F~():=−P(,α)=λ0ϕ(U(,))−αλ1⊤ρ(g())−αλ2⊤ψ(h()). F(y):=-P(y,α)= _0\,φ\! (U(x,y) )-α\, _1 ρ\! (g(y) )-α\, _2 ψ\! (h(y) ). Since ⋆y is a local minimizer of P(⋅,α)P(·,α), it is a local maximizer of F~ F. Finally set λ1′:=αλ1∈ℝ++m,λ2′:=αλ2∈ℝ++p. _1 :=α\, _1 ^m_++, _2 :=α\, _2 ^p_++. Then F~() F(y) equals λ0ϕ(U(,))−(λ1′)⊤ρ((,))−(λ2′)⊤ψ((,)), _0\,φ\! (U(x,y) )-( _1 ) ρ\! (C(x,y) )-( _2 ) ψ\! (T(x,y) ), which is precisely the objective in equation 16. Hence ⋆y is a local maximizer of equation 16 over Y. Since ⋆∈int()y (Y), this is equivalent to local maximality in the ambient space sense. This completes the proof. ∎ Theorem 2.2 (Universality of UMP). Let X and Y be arbitrary nonempty sets. Let f:×→ℝf:X×Y be an objective and let gii∈I≤,g~kk∈I≥,hjj∈J\g_i\_i∈ I_≤, \ g_k\_k∈ I_≥, \h_j\_j∈ J be (possibly empty, countable, or uncountable) families of real–valued constraint functions on ×X×Y. For each fixed ∈x , consider the optimization problem sup∈f(,)s.t.gi(,)≤0(i∈I≤),g~k(,)≥0(k∈I≥),hj(,)=0(j∈J). _y f(x,y) .t. g_i(x,y)≤ 0~(i∈ I_≤),\ \ g_k(x,y)≥ 0~(k∈ I_≥),\ \ h_j(x,y)=0~(j∈ J). (24) Define (with the convention sup∅:=−∞ :=-∞ and maxima taken in the extended reals) U(,):=f(,),(,):=max 0,supi∈I≤gi(,),supk∈I≥(−g~k(,)),U(x,y):=f(x,y), (x,y):= \\,0,\ _i∈ I_≤g_i(x,y),\ _k∈ I_≥ (- g_k(x,y) ) \, (,):=max 0,supj∈J|hj(,)|.T(x,y):= \\,0,\ _j∈ J\,|h_j(x,y)| \. Then for every ∈x , problem equation 24 is equivalent to the utility–maximization problem sup∈U(,)s.t.(,)≤0,(,)=0, _y U(x,y) .t. (x,y)≤ 0, (x,y)=0, (25) in the sense that the feasible sets of equation 24 and equation 25 coincide; hence the optimal values coincide, and whenever maximizers exist, the argmax sets coincide. For minimization problems, replace U by −f-f. Proof. Fix ∈x . Let F():=∈:gi(,)≤0∀i∈I≤,g~k(,)≥0∀k∈I≥,hj(,)=0∀j∈JF(x):= \y :\ g_i(x,y)≤ 0~∀ i∈ I_≤,\ g_k(x,y)≥ 0~∀ k∈ I_≥,\ h_j(x,y)=0~∀ j∈ J \ denote the feasible set of equation 24, and let F^():=∈:(,)≤0,(,)=0 F(x):= \y :\ C(x,y)≤ 0,\ T(x,y)=0 \ denote the feasible set of equation 25. We prove that F()=F^()F(x)= F(x). (i) F()⊆F^()F(x) F(x). Let ∈F()y∈ F(x). Then gi(,)≤0g_i(x,y)≤ 0 for all i∈I≤i∈ I_≤, hence supi∈I≤gi(,)≤0. _i∈ I_≤g_i(x,y)≤ 0. Similarly, g~k(,)≥0 g_k(x,y)≥ 0 for all k∈I≥k∈ I_≥ implies −g~k(,)≤0- g_k(x,y)≤ 0 for all k, hence supk∈I≥(−g~k(,))≤0. _k∈ I_≥ (- g_k(x,y) )≤ 0. Moreover, hj(,)=0h_j(x,y)=0 for all j∈Jj∈ J implies |hj(,)|=0|h_j(x,y)|=0 for all j∈Jj∈ J, hence supj∈J|hj(,)|≤0(with the convention sup∅=−∞). _j∈ J|h_j(x,y)|≤ 0 (with the convention $ =-∞$). By definition, (,)=max 0,supi∈I≤gi(,),supk∈I≥(−g~k(,))=0,C(x,y)= \\,0,\ _i∈ I_≤g_i(x,y),\ _k∈ I_≥ (- g_k(x,y) ) \=0, and (,)=max 0,supj∈J|hj(,)|=0.T(x,y)= \\,0,\ _j∈ J|h_j(x,y)| \=0. Thus ∈F^()y∈ F(x). (i) F^()⊆F() F(x) F(x). Let ∈F^()y∈ F(x). Set A:=supi∈I≤gi(,),B:=supk∈I≥(−g~k(,)),S:=supj∈J|hj(,)|.A:= _i∈ I_≤g_i(x,y), B:= _k∈ I_≥ (- g_k(x,y) ), S:= _j∈ J|h_j(x,y)|. Then (,)=max0,A,B≤0.C(x,y)= \0,A,B\≤ 0. Since 0≤max0,A,B0≤ \0,A,B\ always holds, we have max0,A,B=0 \0,A,B\=0, and in particular A≤0A≤ 0 and B≤0B≤ 0. Using the basic property of the supremum, for every i∈I≤i∈ I_≤ we have gi(,)≤supi∈I≤gi(,)=A≤0,g_i(x,y)≤ _i∈ I_≤g_i(x,y)=A≤ 0, and for every k∈I≥k∈ I_≥ we have −g~k(,)≤supk∈I≥(−g~k(,))=B≤0,- g_k(x,y)≤ _k∈ I_≥ (- g_k(x,y) )=B≤ 0, i.e., g~k(,)≥0 g_k(x,y)≥ 0. Next, (,)=0T(x,y)=0 means 0=(,)=max0,S,0=T(x,y)= \0,S\, hence S≤0S≤ 0. Since |hj(,)|≥0|h_j(x,y)|≥ 0 for every j∈Jj∈ J and |hj(,)|≤S≤0|h_j(x,y)|≤ S≤ 0, it follows that |hj(,)|=0for all j∈J,|h_j(x,y)|=0 all j∈ J, equivalently hj(,)=0h_j(x,y)=0 for all j∈Jj∈ J. Therefore ∈F()y∈ F(x). Combining (i) and (i) yields F()=F^()F(x)= F(x). Since U(,)=f(,)U(x,y)=f(x,y) (and for minimization problems one may equivalently optimize −f-f), the two problems optimize the same objective over the same feasible set. Consequently, their optimal values coincide, and whenever maximizers exist, their argmax sets coincide. ∎ B.2 BL Architecture Theorem 2.3 (Universal Approximation of BL). Let ⊂ℝdX ^d and ⊂ℝmY ^m be compact sets, and let p⋆(∣)p (y ) be any continuous conditional density such that p⋆(∣)>0p (y )>0 for all (,)∈×(x,y) ×Y. Then for any τ>0τ>0 and ε>0 >0, there exists a finite BL architecture (with some depth and width depending on ε ) and a parameter θ⋆θ such that the Gibbs distribution pτ(∣;θ⋆)=exp(BLθ⋆(,)/τ)∫exp(BLθ⋆(,′)/τ)d′p_τ(y ;θ )= (BL_θ (x,y)/τ ) _Y (BL_θ (x,y )/τ )dy (26) satisfies sup∈KL(p⋆(⋅∣)∥pτ(⋅∣;θ⋆))<ε. _x KL (p (· )\,\|\,p_τ(· ;θ ) )< . (27) Proof. Step 0 (bounded log-density). Define f(,):=logp⋆(∣)f(x,y):= p (y ). Since p⋆p is continuous and strictly positive on the compact set ×X×Y, it attains a positive minimum and finite maximum. Hence f∈C(×)f∈ C(X×Y) and is bounded. Step 1 (the BL block contains a one-hidden-layer tanh network). Recall the elementary block ℬ(,;θ):=λ0⊤tanh(u(,))−λ1⊤ReLU(c(,))−λ2⊤|t(,)|.B(x,y;θ):= _0^\! (p_u(x,y) )- _1^\! ReLU (p_c(x,y) )- _2^\! |p_t(x,y) |. (28) Set λ1= _1=0 and λ2= _2=0. Choose u(,)p_u(x,y) to be affine in [;][x;y], i.e. u(,)=W[;]+b∈ℝkp_u(x,y)=W[x;y]+b ^k for some k∈ℕk . Then ℬ(,;θ)=λ0⊤tanh(W[;]+b),B(x,y;θ)= _0 \! (W[x;y]+b ), (29) which is a standard one-hidden-layer tanh network on the compact domain ×X×Y. If λ0∈ℝk _0 ^k is unconstrained, equation 29 is the classical universal approximation class. If instead one imposes λ0≥ _0 0 componentwise, the same expressivity is retained because tanh is odd: for any scalar a∈ℝa , write a=a+−a−a=a^+-a^- with a±≥0a^±≥ 0, and note atanh(h)=a+tanh(h)+a−tanh(−h).a (h)=a^+ (h)+a^- (-h). Since −h-h is affine whenever h is affine, negative coefficients can be realized by duplicating hidden units and keeping the corresponding output weights nonnegative. Thus, up to a constant-factor increase in width, the block class contains signed linear combinations of tanh units. Step 2 (uniform approximation of the target energy). By the universal approximation theorem for single-hidden-layer networks with nonpolynomial activation (e.g., tanh ), for any δ>0δ>0 there exist a width k and parameters θ such that sup(,)∈×|ℬ(,;θ)−τf(,)|<δ. _(x,y) ×Y |B(x,y;θ)-τ f(x,y) |<δ. (30) Define g(,):=ℬ(,;θ)/τg(x,y):=B(x,y;θ)/τ and η:=δ/τη:=δ/τ. Then equation 30 is equivalent to sup(,)∈×|g(,)−f(,)|<η. _(x,y) ×Y |g(x,y)-f(x,y) |<η. (31) Step 3 (uniform KL control). For each ∈x , define q(∣):=exp(g(,))∫exp(g(,′))′.q(y ):= (g(x,y) ) _Y (g(x,y ) )\,dy . (32) The normalizer in equation 32 is finite because g is continuous and Y is compact. Let Zg():=∫exp(g(,′))′.Z_g(x):= _Y (g(x,y ) )\,dy . Since p⋆(⋅∣)p (· ) is a density, without loss of generality we may normalize the energy so that ∫ef(,′)′=1 _Ye^f(x,y )dy =1. From equation 31, for all (,)(x,y), e−η≤eg(,)ef(,)≤eη.e^-η≤ e^g(x,y)e^f(x,y)≤ e^η. Integrating over ∈y yields e−η≤Zg()≤eη,hence|logZg()|≤η,∀∈.e^-η≤ Z_g(x)≤ e^η, | Z_g(x)|≤η, ∀\,x . (33) Moreover, logp⋆(∣)q(∣)=logef(,)eg(,)/Zg()=(f(,)−g(,))+logZg(). p (y )q(y )= e^f(x,y)e^g(x,y)/Z_g(x)= (f(x,y)-g(x,y) )+ Z_g(x). Taking expectation under p⋆(⋅∣)p (· ) and using equation 31 and equation 33 gives KL(p⋆(⋅∣)∥q(⋅∣)) \! (p (· )\, \|\,q(· ) ) =p⋆(⋅∣)[f(,)−g(,)]+logZg() =E_p (· )\! [f(x,Y)-g(x,Y) ]+ Z_g(x) ≤η+η=2η,∀∈. ≤η+η=2η, ∀\,x . (34) Step 4 (choose δ and embed into BL). Choose δ:=ετ/4δ:= τ/4, so that η=δ/τ=ε/4η=δ/τ= /4. Then equation 34 implies sup∈KL(p⋆(⋅∣)∥q(⋅∣))≤2η=ε/2<ε. _x KL\! (p (· )\, \|\,q(· ) )≤ 2η= /2< . Finally, the density q(⋅∣)q(· ) equals the Gibbs distribution equation 26 with energy BLθ⋆(,):=ℬ(,;θ)BL_θ (x,y):=B(x,y;θ) (a finite BL architecture containing a single block), and temperature τ. This proves the claim. ∎ B.3 Identifiable Behavior Learning (IBL) B.3.1 Setup and Assumption Input–output space and data. Let ⊂ℝdxX ^d_x and ⊂ℝdyY ^d_y be compact sets. Assume the data distribution PX,YP_X,Y is supported on ×X×Y, and that there exists a point z0=(x0,y0)z_0=(x_0,y_0) in the interior of its support; that is, some open neighborhood of z0z_0 has positive PX,YP_X,Y-measure. All expectations are taken with respect to PX,YP_X,Y unless otherwise specified. Parameter space and polynomial feature maps. The parameter space factorizes as Θ:=ΘU×ΘC×ΘT×∘. \;:=\; _U× _C× _T×W_ . For θU∈ΘU _U∈ _U, θC∈ΘC _C∈ _C, and θT∈ΘT _T∈ _T, we define polynomial feature maps pu:×→ℝdu,pc:×→ℝdc,pt:×→ℝdt,p_u:X×Y ^d_u, p_c:X×Y ^d_c, p_t:X×Y ^d_t, each of fixed degree and injective in their coefficients (i.e., distinct coefficients yield distinct functions). For a single block, θU,θC,θT _U, _C, _T correspond to the parameters of the U, C, and T terms together with their respective external multipliers (e.g., penalty weights λ). For a deep network composed of multiple blocks, θ=(θU,θC,θT)θ=( _U, _C, _T) denotes the collection of all block-level parameters across the hierarchy, where θU _U aggregates the parameters of all U-terms, θC _C those of all C-terms, and θT _T those of all T-terms (each including their associated multipliers). The output component ∘W_ corresponds to the affine transformation in the final layer: ∘=ℝd′W_ =R^d for single-output prediction, and ∘=ℝd′×mW_ =R^d × m for m-way classification, where d′d is the output dimension induced by the preceding network, whether shallow or deep. Identifiable base block. Let λ0∈ℝdu _0 ^d_u, λ1∈ℝdc _1 ^d_c, and λ2∈ℝdt _2 ^d_t denote nonnegative weight vectors, treated as learnable parameters. We instantiate the identifiable modular block ℬid(x,y;θ)=λ0⊤tanh(pu(x,y))−λ1⊤softplus(pc(x,y))−λ2⊤(pt(x,y))⊙2,B^id(x,y;θ)\;=\; _0^\! \! (p_u(x,y) )\;-\; _1^\! softplus\! (p_c(x,y) )\;-\; _2^\! (p_t(x,y) ) 2, (35) where (⋅)⊙2(·) 2 denotes elementwise squaring. By construction, the tanh and softplussoftplus heads are strictly monotone in their arguments, while the quadratic head is even. We assume that each polynomial feature map ∙(x,y)p_ (x,y) contains no nonzero monomial independent of y; that is, no feature is a pure function of x or a constant. This ensures that ℬid(x,y)B^id(x,y) is nonconstant in y unless all weights vanish. Architectures. We implement IBL in three architectural forms, each producing a compositional utility function over (x,y)(x,y). • IBL(Single): A single block is used as the compositional utility, IBL(x,y):=ℬid(x,y).IBL(x,y):=B^id(x,y). • IBL(Shallow): Shallow IBL uses one or two stacked layers of parallel blocks. For instance, a first layer 1id(x,y):=[ℬ1,1id(x,y),…,ℬ1,d1id(x,y)]⊤∈ℝd1B^id_1(x,y):=[\,B^id_1,1(x,y),…,B^id_1,d_1(x,y)\,] ^d_1 feeds into a bias-free affine map IBLShallow(x,y):=1∘1id(x,y),IBL_Shallow(x,y):=W _1\,B^id_1(x,y), where 1∘∈ℝm×d1W _1 ^m× d_1 for classification and 1∘∈ℝ1×d1W _1 ^1× d_1 for scalar output. • IBL(Deep): Deep IBL extends the construction to depth L>2L>2, recursively defined as IBL(x,y):=L∘⋅Lid(⋯2id(1id(x,y))⋯),IBL(x,y):=W _L·B^id_L (·s\,B^id_2(B^id_1(x,y))·s ), where each ℓidB^id_ stacks parallel blocks ℬℓ,iid(x,y)B^id_ ,i(x,y), and L∘W _L is a bias-free affine transformation. The cases L=1L=1 and L=2L=2 recover the Single and Shallow architectures, respectively. Induced conditional model. Let IBL(x,y)IBL(x,y) denote the compositional utility function produced by the chosen architecture (Single, Shallow, or Deep). It induces the conditional Gibbs distribution (Discrete y∈[m])p(y∣x) (Discrete $y∈[m]$) p(y x) =softmaxyIBL(x,y), =softmax_y\IBL(x,y)\, (36) (Continuous y)p(y∣x) (Continuous $y$) p(y x) =expIBL(x,y)/τ∫expIBL(x,y~)/τy~,τ>0 fixed. = \IBL(x,y)/τ\ _Y \IBL(x, y)/τ\\,d y, τ>0 fixed. (37) Here τ is a fixed temperature parameter. Thus, IBL predicts by defining a compositional utility landscape whose Gibbs distribution governs y given x. Quotient parameter space. Definition B.1 (Symmetry Quotient Space). Define the equivalence relation ∼ on Θ as the smallest relation satisfying θt∼θt′⟺pt(i)(x,y;θt(i))⊙2=pt(i)(x,y;θt′(i))⊙2for all i and (x,y). _t _t p_t^(i)(x,y; _t^(i)) 2=p_t^(i)(x,y; _t (i)) 2 all i and (x,y). The corresponding quotient space is Θ¯:=Θ/∼. := /\! . Explanation. The T-component is designed to encode equality constraints, which are symbolically equations. Flipping the overall sign of such a constraint leaves the equation unchanged, so different parameterizations that differ only by sign should be regarded as equivalent. Definition B.2 (Scale-Invariant Quotient Space). Define the equivalence relation ≈ on Θ¯ by θ¯≈θ¯′⟺∃c>0 such that (x,y;θ¯)=c(x,y;θ¯′). θ≈ θ ∃\,c>0\; such that \; s(x,y; θ)=c\, s(x,y; θ ). The scale-invariant quotient space is then given by Θ~:=Θ¯/≈. := /\!≈. Explanation. In classification, predictions depend only on relative compositional utility differences between candidate labels. From a technical perspective, quotienting out global shifts or uniform scalings is necessary: without this identification, the cross-entropy loss admits redundant parameterizations that differ only by such transformations. At the same time, this quotient is natural and harmless, since it does not eliminate informative ratios between classes but merely discards absolute levels or scales that play no role in the softmax decision rule. Loss Functions. We adopt a hybrid loss to simultaneously accommodate discrete and continuous outputs. Specifically, cross-entropy (CE) is applied to discrete targets, while denoising score matching (DSM) is applied to continuous targets. Let γc,γd≥0 _c, _d≥ 0 with γc+γd>0 _c+ _d>0. The population risk, defined on the quotient parameter space, is given by ℳ(θ¯)=γd[−logpθ(Y∣X)]+γc[DSM(θ;X)],θ∈π−1(θ¯),M( θ)\;=\; _d\,E\! [- p_θ(Y X) ]\;+\; _c\,E\! [S_DSM(θ;X) ], θ∈π^-1( θ), (38) where π denotes the canonical projection from the original parameter space onto its quotient. For continuous outputs Y∈⊆ℝdyY ^d_y, DSM is implemented by perturbing the target with additive Gaussian noise Y~=Y+ε Y=Y+ , ε∼(0,σ2I) (0,σ^2I), and penalizing the squared discrepancy between the model score and the corresponding denoising score: DSM(θ;X)=12σ2ε[∥∇y~logpθ(y~∣X)+1σ2(Y−Y~)∥2|X,Y].S_DSM(θ;X)= 12σ^2\,E_ \! [ \| _ y p_θ( y X)+ 1σ^2(Y- Y) \|^2\, |\,X,Y ]. (39) In classification-only settings we set γc=0 _c=0 (pure CE), while in regression-only settings we set γd=0 _d=0 (pure DSM). For a single observation Z=(X,Y)Z=(X,Y), we define the per-sample loss as ℓ(θ;Z):=γd[−logpθ(Y∣X)]+γcDSM(θ;X). (θ;Z):= _d\, [- p_θ(Y X) ]+ _c\,S_DSM(θ;X). (40) The empirical criterion then takes the standard M-estimation form Q^n(θ)=1n∑i=1nℓ(θ;Zi),Zi=(Xi,Yi). Q_n(θ)\;=\; 1n _i=1^n (θ;Z_i), Z_i=(X_i,Y_i). (41) Key Assumptions. Assumption B.1 (Global Atomic Independence and Injectivity). Let Ψ¯ be the atomic parameter quotient. 1. Injectivity on the quotient. The map Ψ¯→ℝ× ^X×Y, ψ¯↦gψ¯ ψ g_ ψ, is injective. 2. Linear Independence. Atomic linear independence. Any finite collection of pairwise distinct atoms gψ¯ii=1r\g_ ψ_i\_i=1^r with ψ¯i∈Ψ¯ ψ_i∈ is linearly independent in ℝ×R^X×Y. 3. Minimality. In all model instances we only consider minimal representations: no duplicate atoms and its corresponding linear coefficient in the mixture is nonzero. 4. Canonical ordering. For each model instance, a fixed canonical ordering is imposed on the atom list. Explanation. Assumption B.1 treats each identifiable block ℬidB^id as an atomic building unit and imposes four structural requirements on representations built from these atoms. Together, these four conditions define a non-ambiguous, non-redundant, and canonical algebra of atoms: after quotienting by the natural symmetries, every model constructed from ℬB-blocks admits a unique minimal representation (up to the prescribed equivalences). This structural regularity is the foundation on which identifiability statements are built: it guarantees that observing the model output (or the objective it optimizes) allows one, in principle, to recover the underlying atomic components and their coefficients in the appropriate quotient sense. Practical remark. In practice, these conditions can be encouraged or approximately enforced in two complementary ways. First, the design of atomic classes (choice of polynomial bases, interaction terms, and activation heads) can be chosen so that injectivity and linear independence are more plausible by construction. Second, model selection and post-processing (e.g., pruning atoms with near-zero coefficients, enforcing a deterministic tie-breaking rule for ordering) can be applied after training to realize minimality and canonical ordering. These practical measures make the theoretical assumptions operationally meaningful in empirical applications. B.3.2 Proof of Theorems Lemma B.1 (Identifiability of Linear Combinations). Let Z be a set. For each j=1,…,mj=1,…,m, let Φj _j be a parameter space and define atomic functions gψ:=f(⋅;ϕj),ψ=(j,ϕj)∈Ψ,g_ψ:=f(·; _j), ψ=(j, _j)∈ , where Ψ:=⨆j=1mΦj := _j=1^m _j is the disjoint union. Let Ψ¯ be the quotient atomic parameter space, and denote its elements by ψ¯∈Ψ¯ ψ∈ . Define the quotient parameter space of the model as Ξ¯:=∏j=1m((ℝ∖0)×Ψ¯),ξ¯=((a1,ψ¯1),…,(am,ψ¯m)). := _j=1^m ((R \0\)× ), ξ=((a_1, ψ_1),…,(a_m, ψ_m)). The associated linear combination model is Sξ¯:=∑j=1majgψ¯j.S_ ξ:= _j=1^ma_jg_ ψ_j. By virtue of Assumption B.1, the model is identifiable in the quotient parameter space Ξ¯ : if Sξ¯≡Sξ¯′S_ ξ≡ S_ ξ on Z, then ξ¯=ξ¯′ ξ= ξ . Proof. Suppose Sξ¯≡Sξ¯′S_ ξ≡ S_ ξ on Z, i.e., ∑j=1majg(j,ϕj)−∑j=1maj′g(j,ϕj′)≡0. _j=1^ma_j\,g_(j, _j)- _j=1^ma _j\,g_(j,φ _j)≡ 0. Let U be the set of distinct atoms in the quotient Ψ¯ that appear on either side, and for each ψ¯∈ ψ let β(ψ¯):=∑j:[j,ϕj]=ψ¯aj−∑j′:[j′,ϕj′]=ψ¯aj′β( ψ):= _ subarraycj:\,[j, _j]= ψ subarraya_j\;-\!\!\!\!\! _ subarraycj :\,[j ,φ _j ]= ψ subarraya _j be the net coefficient of gψ¯g_ ψ. Then ∑ψ¯∈β(ψ¯)gψ¯≡ 0. _ ψ β( ψ)\,g_ ψ\;≡\;0. By the linear independence condition (Assumption B.1:2) of pairwise distinct atoms in Ψ¯ , we must have β(ψ¯)=0β( ψ)=0 for all ψ¯∈ ψ . Furthermore, by the Minimality requirement (Assumption B.1:3), each ψ¯ ψ appears exactly once on each side and with nonzero coefficient. Thus the two sides must contain the exact same list of coefficient–atom pairs (aj,ψ¯j)j=1m\(a_j, ψ_j)\_j=1^m, and since a canonical ordering is imposed (Assumption B.1:4), it follows that ξ¯=ξ¯′. ξ= ξ . ∎ Theorem B.1 (Identifiability of IBL(Single)). The IBL(Single) architecture uses the atom set tanh(pu,i),softplus(pc,i),(pt,i)2:i=1,…,du;i=1,…,dc;i=1,…,dt. \\, (p_u,i),\ softplus(p_c,i),\ (p_t,i)^2:\ i=1,…,d_u;\ i=1,…,d_c;\ i=1,…,d_t\, \. Under Assumption B.1, the model is identifiable in the quotient space Θ¯ : if ℬθid≡ℬθ′idB^id_θ ^id_θ on ×X×Y, then θ=θ′θ=θ in Θ¯ . Proof. Write ℬθid=∑j=1majf(⋅;ϕj),m:=du+dc+dt,B^id_θ= _j=1^ma_j\,f(·; _j), m:=d_u+d_c+d_t, where each f(⋅;ϕj)f(·; _j) is one of the atoms tanh(pu,i) (p_u,i), softplus(pc,i)softplus(p_c,i), or (pt,i)2(p_t,i)^2, and aja_j is the corresponding entry in (λ0,λ1,λ2)( _0, _1, _2), with a fixed ordering over all indices. If ℬθid≡ℬθ′idB^id_θ ^id_θ on ×X×Y, then Lemma B.1 and Assumption B.1 imply that all atoms and coefficients must agree in the quotient atomic space Ψ¯ . Since the ordering is fixed, this implies θ=θ′θ=θ in Θ¯ . ∎ Theorem B.2 (Identifiability of IBL(Shallow)). The IBL(Shallow) architecture uses the atom set ℬθ1,jid(x,y)j=1d1, \\,B^id_ _1,j(x,y)\, \_j=1^d_1, where each ℬθ1,jid:×→ℝB^id_ _1,j:X×Y is a single-block IBL module parametrized by θ1,j∈Θ1 _1,j∈ _1. The full parameter is denoted θ:=((θ1,1,…,θ1,d1),1∘)∈Θ:=(Θ1)d1×ℝm×d1.θ:= (( _1,1,…, _1,d_1),\ W_1 )∈ :=( _1)^d_1×R^m× d_1. Under Assumption B.1, the mapping θ↦IBLShallowθ _Shallow is identifiable in the quotient space Θ¯ : if IBLShallow(x,y;θ)≡IBLShallow(x,y;θ′)on ×,IBL_Shallow(x,y;θ) _Shallow(x,y;θ ) X×Y, then θ=θ′in Θ¯.θ=θ . Proof. Write the k-th output component as a linear combination of atoms: sθ(k)(x,y)=∑j=1d1wj(k)ℬθ1,jid(x,y),k=1,…,m,s^(k)_θ(x,y)= _j=1^d_1w^(k)_j\,B^id_ _1,j(x,y), k=1,…,m, where wj(k)w^(k)_j denotes the (k,j)(k,j)-th entry of 1∘W_1 . Suppose two parameter tuples (1∘,θ1,jj=1d1)(W_1 ,\ _1,j\_j=1^d_1) and (1∘′,θ1,j′j=1d1)(W_1 \, ,\θ _1,j\_j=1^d_1) yield identical vector scores on ×X×Y. Then for each k, we have sθ(k)≡sθ′(k)s^(k)_θ≡ s^(k)_θ on ×X×Y. Fix any k. Under Assumption B.1, Lemma B.1 ensures that the coefficient–atom pairs (wj(k),ℬθ1,jid)j=1d1\(w^(k)_j,\,B^id_ _1,j)\_j=1^d_1 are uniquely determined (up to equivalence in the quotient Θ¯ ). In particular, for each j=1,…,d1j=1,…,d_1, we must have wj(k)=wj′(k),ℬθ1,jid≡ℬθ1,j′id.w^(k)_j=w (k)_j, ^id_ _1,j ^id_θ _1,j. Because this holds for all k=1,…,mk=1,…,m, it follows that 1∘=1∘′W_1 =W_1 \, and θ1,j=θ1,j′ _1,j=θ _1,j in the quotient parameter space for all j. Thus θ=θ′θ=θ in Θ¯ , establishing full identifiability under fixed ordering. ∎ Theorem B.3 (Identifiability of IBL(Deep)). Fix integers L>2L>2 and widths d1,…,dL−1d_1,…,d_L-1. The IBL(Deep) architecture uses the final-layer atom set ℬϑL,jid(x,y)j=1dL⊂ℝ×, \\,B^id_ _L,j(x,y)\, \_j=1^d_L ^X×Y, where each ℬϑL,jid:ℝdL−1→ℝB^id_ _L,j:R^d_L-1 is a scalar-valued block applied to the output of layer L−1L-1. Only the first-layer blocks (ℓ=1 =1) are IBL(Single) modules as in Theorem B.1. For architectures with skip connections, the final-layer atoms can be extended to include skipped features (e.g., from earlier layers), which are treated as elements of ℬϑL,jid(x,y)j=1dL. \\,B^id_ _L,j(x,y)\, \_j=1^d_L. The full parameter is θ:=(ϑℓ,jℓ=1,j=1L,dℓ,out)∈Θ:=∏ℓ=1L(Θ1)dℓ×ℝm×dL.θ:= (\ _ ,j\_ =1,j=1^L,d_ ,\ W_out )∈ := _ =1^L( _1)^d_ ×R^m× d_L. Under Assumption B.1, the mapping θ↦IBLDeep(x,y;θ)θ _Deep(x,y;θ) is identifiable in the quotient space Θ¯ . Proof. Under the given architecture, the IBL(Deep) model ultimately takes the form s(k)(x,y)=∑j=1dLwj(k)ℬϑL,jid(x,y),k=1,…,m,s^(k)(x,y)\;=\; _j=1^d_Lw^(k)_j\,B^id_ _L,j(x,y), k=1,…,m, where each ℬϑL,jidB^id_ _L,j is a scalar-valued function applied to the output of preceding layers. By treating the set ℬϑL,jid(x,y)j=1dL\B^id_ _L,j(x,y)\_j=1^d_L as the atom set, we reduce the model to an IBL(Shallow) form: (x,y)=outL(x,y).s(x,y)=W_out\,B_L(x,y). Under Assumption B.1, Theorem B.2 applies, implying that the full parameter θ=(ϑℓ,jℓ,j,out)θ=(\ _ ,j\_ ,j,\,W_out) is identifiable in the quotient space Θ¯ . ∎ Theorem 2.4 (Identifiability of IBL). Under Assumption B.1, the architectures IBL(Single), IBL (Shallow), and IBL(Deep) are all identifiable in the quotient space Θ¯ . Proof. Immediate from Theorems B.1, B.2, and B.3. ∎ Theorem 2.5 (Loss Identifiability of IBL). Let IBLθ(x,y)IBL_θ(x,y) denote an IBL model, and consider the conditional Gibbs distribution pθ(y∣x)=exp(IBLθ(x,y))∫exp(IBLθ(x,y′))y′.p_θ(y x)= \! (IBL_θ(x,y) ) _Y \! (IBL_θ(x,y ) )\,dy . Define the population risk on the symmetry quotient Θ¯ as in equation 38. Assume that the parameter space Θ is compact. Then, under Assumption B.1, the following holds: (i) If γc>0 _c>0, the risk functional ℳM admits a unique minimizer in Θ¯ . Moreover, ℳ(θ¯1)=ℳ(θ¯2)⟹θ¯1=θ¯2.M( θ_1)=M( θ_2)\;\; \;\; θ_1= θ_2. (i) If γc=0 _c=0, the risk functional ℳM admits a unique minimizer in the scale-invariant quotient Θ~ . Moreover, ℳ(θ~1)=ℳ(θ~2)⟹θ~1=θ~2.M( θ_1)=M( θ_2)\;\; \;\; θ_1= θ_2. Proof. Under Assumption B.1, the IBL architecture is identifiable modulo the symmetry group defined by Θ¯ , as established in Theorem 2.4. Let θ∙∈argminθ∈Θℳ(θ)θ ∈ _θ∈ M(θ) and set p⋆(⋅∣x):=pθ∙(⋅∣x)p (· x):=p_θ (· x). Since Θ is compact and the loss ℳM is continuous, a global minimizer exists. We show that it is unique in the stated quotient. Case γc>0 _c>0. At any minimizer we have both pθ(⋅∣x)=p⋆(⋅∣x)p_θ(· x)=p (· x) and ∇ylogpθ(⋅∣x)=∇ylogp⋆(⋅∣x) _y p_θ(· x)= _y p (· x) a.e. Since ∇ylogpθ(y∣x)=∇yIBLθ(x,y)−∇ylogZθ(x)=∇yIBLθ(x,y), _y p_θ(y x)= _yIBL_θ(x,y)- _y Z_θ(x)= _yIBL_θ(x,y), (the partition function Zθ(x)Z_θ(x) is y-independent), score equality yields ∇y(IBLθ−IBLθ∙)(y;x)=0 _y (IBL_θ-IBL_θ )(y;x)=0 a.e. IBL contains no y-independent terms. Therefore, IBLθ(x,y)=IBLθ∙(x,y)a.e.IBL_θ(x,y)=IBL_θ (x,y) .e. By Theorem 2.4 (identifiability in Θ¯ ), the minimizer is unique in Θ¯ ; in particular, ℳ(θ¯1)=ℳ(θ¯2)⟹θ¯1=θ¯2.M( θ_1)=M( θ_2)\ \ θ_1= θ_2. Case γc=0 _c=0. Here, ℳM reduces to the cross-entropy risk, which is minimized if and only if pθ(⋅∣x)=p⋆(⋅∣x)p_θ(· x)=p (· x) almost everywhere.The cross-entropy loss depends on IBLθ(x,y)IBL_θ(x,y) only through its relative values across y, and is invariant under additive shifts and positive rescalings of the compositional utility. Hence, the loss depends only on the equivalence class θ~∈Θ~ θ∈ \. As a result, ℳ(θ~1)=ℳ(θ~2)⟹θ~1=θ~2.M( θ_1)=M( θ_2) θ_1= θ_2. i.e., the minimizer is unique in Θ~ . Hence, the minimizer is unique in the stated quotient space. This completes the proof. ∎ Theorem B.4 (Uniform M-estimation consistency (Newey & McFadden, 1994, Theorem 2.1)). Let (,d)(A,d) be a compact metric space, and let L^n:→ℝ L_n:A be a sequence of random objective functions, with population objective L:→ℝL:A such that: 1. L(α)L(α) is uniquely minimized at α⋆∈α ; 2. A is compact; 3. L(α)L(α) is continuous; 4. L^n(α)→L(α) L_n(α) pL(α) uniformly in α∈α . Then any sequence α^n∈argminα∈L^n(α) α_n∈ _α L_n(α) satisfies α^n→α⋆ α_n pα . Theorem B.5 (Consistency of IBL). Let ℳM be the population risk defined in equation 38, and let ℳnM_n denote its empirical analogue. Suppose: 1. (Xi,Yi)i=1n\(X_i,Y_i)\_i=1^n are i.i.d. samples; 2. Θ is compact; 3. θ↦ℳ(θ)θ (θ) is continuous, and the loss class admits an integrable envelope such that supθ∈Θ|ℳn(θ)−ℳ(θ)|→ 0; _θ∈ |M_n(θ)-M(θ) |\; p\;0; Let Ξ denote the relevant quotient space (Θ¯ if γc>0 _c>0, Θ~ if γc=0 _c=0), and let θ^n∈argminθ∈Θℳn(θ) θ_n∈ _θ∈ M_n(θ) and θ∙∈argminθ∈Θℳ(θ)θ ∈ _θ∈ M(θ). Then θ^n→θ∙in Ξ,ℳ(θ^n)→ℳ(θ∙). θ_n\; p\;θ , ( θ_n)\; p\;M(θ ). If the model is correctly specified (the data law is realized by some θ⋆∈Θθ ∈ ), then θ∙=θ⋆θ =θ in Ξ , so θ^n→θ⋆ θ_n pθ . Proof. Let Ξ denote the relevant quotient space: Ξ=Θ¯ = if γc>0 _c>0 and Ξ=Θ~ = if γc=0 _c=0. Let π:Θ→Ξπ: → be the canonical quotient map. Since Θ is compact and π is continuous and onto, Ξ is compact. By assumption, ℳM and ℳnM_n are invariant under the corresponding symmetry, hence they factor through π: ℳ~(ξ):=ℳ(θ),ℳ~n(ξ):=ℳn(θ)(any θ∈π−1(ξ)). M(ξ):=M(θ), M_n(ξ):=M_n(θ) (any θ∈π^-1(ξ)). These are well-defined and continuous on Ξ because ℳM is continuous on Θ . Moreover, supξ∈Ξ|ℳ~n(ξ)−ℳ~(ξ)|≤supθ∈Θ|ℳn(θ)−ℳ(θ)|→ 0, _ξ∈ | M_n(ξ)- M(ξ) |\;≤\; _θ∈ |M_n(θ)-M(θ) |\; p\;0, so uniform convergence in probability holds on Ξ . By Loss Identifiability of IBL (Theorem 2.5), ℳ~ M has a unique minimizer ξ∙∈Ξξ ∈ . Let ξ^n∈argminξ∈Ξℳ~n(ξ) ξ_n∈ _ξ∈ M_n(ξ) (equivalently, choose θ^n∈argminθ∈Θℳn(θ) θ_n∈ _θ∈ M_n(θ) and set ξ^n=π(θ^n) ξ_n=π( θ_n)). Then the conditions of Theorem B.4 hold on the compact metric space (Ξ,d)( ,d), whence ξ^n→ξ∙. ξ_n pξ . Since ℳ~ M is continuous on Ξ and ℳ~(ξ^n)=ℳ(θ^n) M( ξ_n)=M( θ_n), ℳ~(ξ∙)=ℳ(θ∙) M(ξ )=M(θ ) for any representative θ∙∈π−1(ξ∙)θ ∈π^-1(ξ ), we also obtain ℳ(θ^n)=ℳ~(ξ^n)→ℳ~(ξ∙)=ℳ(θ∙).M( θ_n)= M( ξ_n) p M(ξ )=M(θ ). If the model is correctly specified (there exists θ⋆∈Θθ ∈ inducing the data law), the strict propriety of the CE/DSM terms implies that the unique minimizer in the quotient is the class of θ⋆θ ; hence θ^n θ_n converges in probability to θ⋆θ in the corresponding quotient space. ∎ Theorem B.6 (Universal Approximation of IBL). Let ⊂ℝdX ^d and ⊂ℝmY ^m be compact sets, and let p⋆(y∣x)p (y x) be any continuous conditional density such that p⋆(y∣x)>0p (y x)>0 for all (x,y)∈×(x,y) ×Y. Then for any τ>0τ>0 and ε>0 >0, there exists a finite IBL architecture (with some depth and width depending on ε ) and a parameter θ⋆θ such that the Gibbs distribution pτ(y∣x;θ⋆)=exp(IBLθ⋆(x,y)/τ)∫exp(IBLθ⋆(x,y′)/τ)dy′p_τ(y x;θ )= (IBL_θ (x,y)/τ ) _Y (IBL_θ (x,y )/τ )dy (42) satisfies supx∈KL(p⋆(⋅∣x)∥pτ(⋅∣x;θ⋆))<ε. _x KL (p (· x)\,\|\,p_τ(· x;θ ) )< . (43) Proof. The argument follows the same construction as in the proof of Theorem 2.3, with only notational modifications due to the IBL parameterization. For brevity, the details are omitted. ∎ Lemma B.2 (Sieve Approximation Lemma). Let :Θ→[0,∞)C: →[0,∞) be a complexity measure on the parameter space, and let (cn)n≥1(c_n)_n≥ 1 be a nondecreasing sequence with cn↑∞c_n ∞. Define the sieve Θn:=θ∈Θ:(θ)≤cn, _n:=\θ∈ :C(θ)≤ c_n\, and for a fixed data-generating distribution p†p , set δn(p†):=infθ∈Θnsupx∈KL(p†(⋅∣x)∥pθ(⋅∣x)). _n(p )\ :=\ _θ∈ _n\ _x KL\! (p (· x)\,\|\,p_θ(· x) ). Then the following are equivalent: 1. Sieve universal approximation: For every ε>0 >0 there exists a constant Cε<∞C_ <∞ such that infθ:(θ)≤Cεsupx∈KL(p†(⋅∣x)∥pθ(⋅∣x))<ε. _θ:\,C(θ)≤ C_ \ _x KL\! (p (· x)\,\|\,p_θ(· x) )\ <\ . 2. Vanishing approximation error: δn(p†)↓0\ _n(p ) 0 as n→∞n→∞. Moreover, if each Θn _n is compact and θ↦supxKL(p†∥pθ)θ _xKL(p \|p_θ) is continuous on Θn _n, then the infimum in δn(p†) _n(p ) is attained for every n. Proof. (1) ⇒ (2). Fix ε>0 >0 and let Cε(p†)C_ (p ) be as in (i). Since cn↑∞c_n ∞, there exists N such that cn≥Cε(p†)c_n≥ C_ (p ) for all n≥Nn≥ N. Hence Θn⊇θ:(θ)≤Cε(p†) _n \θ:C(θ)≤ C_ (p )\ for all n≥Nn≥ N, and therefore δn(p†)=infθ∈ΘnsupxKL(p†∥pθ)≤infθ:(θ)≤Cε(p†)supxKL(p†∥pθ)<ε, _n(p )\ =\ _θ∈ _n\ _xKL(p \|p_θ)\ ≤\ _θ:\,C(θ)≤ C_ (p )\ _xKL(p \|p_θ)\ <\ , for all n≥Nn≥ N. Since (δn)( _n) is nonincreasing in n (because Θn↑ _n ), it follows that δn(p†)↓0 _n(p ) 0. (2) ⇒ (1). Fix ε>0 >0. By (i) choose N such that δN(p†)<ε _N(p )< . Set Cε(p†):=cNC_ (p ):=c_N. Then infθ:(θ)≤Cε(p†)supxKL(p†∥pθ)≤infθ∈ΘNsupxKL(p†∥pθ)=δN(p†)<ε, _θ:\,C(θ)≤ C_ (p )\ _xKL(p \|p_θ)\ ≤\ _θ∈ _N\ _xKL(p \|p_θ)\ =\ _N(p )\ <\ , which is (i). The attainment statement follows immediately from compactness of Θn _n and continuity of θ↦supxKL(p†∥pθ)θ _xKL(p \|p_θ) on Θn _n. ∎ Theorem B.7 (Universal Consistency of IBL). Consider a parameter space Θ for a class of IBL models, and let :Θ→[0,∞)C: →[0,∞) be a lower semi-continuous complexity measure (e.g., network depth, width, or parameter norm). Let (cn)n≥1(c_n)_n≥ 1 be a nondecreasing sequence with cn↑∞c_n ∞, and define the sieve Θn:=θ∈Θ:(θ)≤cn. _n:=\θ∈ :C(θ)≤ c_n\. Assume: 1. The map θ↦supxKL(p†∥pθ)θ _xKL(p \|p_θ) is continuous on each compact Θn _n. 2. The sequence of empirical minimizers θ^n\ θ_n\ is relatively compact in ⋃nΘn _n _n, as ensured by the uniform LLN together with compactness and continuity. Then for any admissible data-generating distribution p†p satisfying the regularity assumptions of Theorem B.6, the IBL posterior sequence pθ^n\p_ θ_n\ satisfies supx∈KL(p†(⋅∣x)∥pθ^n(⋅∣x))→0, _x KL\! (p (· x)\,\|\,p_ θ_n(· x) ) p0, i.e. pθ^n\p_ θ_n\ converges to p†p uniformly in x (in KL). Proof. Fix an admissible data law p†p (satisfying the regularity of Theorem B.6). For θ∈⋃nΘnθ∈ _n _n define F(θ):=supx∈KL(p†(⋅∣x)∥pθ(⋅∣x)),δn:=infθ∈ΘnF(θ).F(θ)\;:=\; _x KL\! (p (· x)\,\|\,p_θ(· x) ), _n\;:=\; _θ∈ _nF(θ). Then Theorem B.6 and Lemma B.2 together imply that δn↓0 _n 0. By assumption 1, F is continuous on each compact Θn _n. Let θ^n∈argminθ∈Θnℳn(θ) θ_n∈ _θ∈ _nM_n(θ) be any sequence of ERM solutions. We show F(θ^n)→0F( θ_n) p0. Step 1 (subsequence reduction and precompactness). Take an arbitrary subsequence (θ^nk)k( θ_n_k)_k. By assumption 2 there exists a further subsequence, still denoted (θ^nk)k( θ_n_k)_k, and a (possibly k-dependent) index set Nk≤nkN_k≤ n_k with a parameter limit θ∞∈ΘN _∞∈ _N (for some finite N) such that θ^nk→θ∞ θ_n_k→ _∞ in probability. Passing to a further subsequence if needed, we may assume Nk≡N_k≡ N. Step 2 (risk domination against ΘN _N-approximants). For each k pick θk∈ΘN _k∈ _N with F(θk)≤δN+1/kF( _k)≤ _N+1/k (attainment follows from compactness and continuity of F on ΘN _N). By the ERM property and uniform LLN on ΘN _N, ℳ(θ^nk)≤ℳ(θk)+op(1)(k→∞).M( θ_n_k)\;≤\;M( _k)+o_p(1) (k→∞). Assume (w.l.o.g.) the CE component is present with a positive weight, so that the population risk decomposes as ℳ(θ)=const+γdX[KL(p†(⋅∣X)∥pθ(⋅∣X))]+γcℒDSM(θ),M(θ)\;=\;const+ _d\,E_X\! [KL\! (p (· X)\,\|\,p_θ(· X) ) ]\;+\; _c\,L DSM(θ), with γd>0 _d>0 (the DSM-only case is handled analogously by replacing KL with Fisher divergence). Using X[KL(⋅∥⋅)]≤F(⋅)E_X[KL(·\|·)]≤ F(·), we obtain lim supk→∞X[KL(p†(⋅∣X)∥pθ^nk(⋅∣X))]≤lim supk→∞F(θk)≤δN. _k→∞E_X\! [KL\! (p (· X)\,\|\,p_ θ_n_k(· X) ) ]\;≤\; _k→∞F( _k)\;≤\; _N. Hence, along the subsequence, X[KL(p†(⋅∣X)∥pθ^nk(⋅∣X))]→0.E_X\! [KL\! (p (· X)\,\|\,p_ θ_n_k(· X) ) ] p0. Step 3 (identification of the subsequential limit). By continuity of the model map θ↦pθ(⋅∣x)θ p_θ(· x) (from Theorem B.6 regularity) and bounded convergence, X[KL(p†(⋅∣X)∥pθ∞(⋅∣X))]=0.E_X\! [KL\! (p (· X)\,\|\,p_ _∞(· X) ) ]=0. Thus f(x):=KL(p†(⋅∣x)∥pθ∞(⋅∣x))f(x):=KL\! (p (· x)\,\|\,p_ _∞(· x) ) equals 0 for PXP_X-a.e. x. Since f is continuous on compact X (by the same regularity) and PXP_X has full support (admissible law), we conclude f(x)≡0f(x)≡ 0 on X, i.e. F(θ∞)=supx∈f(x)=0.F( _∞)= _x f(x)=0. Step 4 (conclude F(θ^nk)→0F( θ_n_k)→ 0 in probability, hence F(θ^n)→0F( θ_n)→ 0 in probability). By assumption 1 continuity of F on ΘN _N and θ^nk→θ∞ θ_n_k→ _∞ in probability, we have F(θ^nk)→F(θ∞)=0F( θ_n_k) pF( _∞)=0. Since the original subsequence was arbitrary and every subsequence admits a further subsequence with F(θ^nk)→0F( θ_n_k) p0, the full sequence satisfies F(θ^n)→0F( θ_n) p0. Therefore, supx∈KL(p†(⋅∣x)∥pθ^n(⋅∣x))→0, _x KL\! (p (· x)\,\|\,p_ θ_n(· x) ) p0, i.e. pθ^n(⋅∣x)→p†(⋅∣x)p_ θ_n(· x)→ p (· x) uniformly in x in KL. ∎ Theorem B.8 (Asymptotic normality of extremum estimators (Newey & McFadden, 1994, Theorem 3.1)). Suppose that the estimator θ^n θ_n satisfies θ^n→θ0 θ_n p _0, and: 1. θ0 _0 lies in the interior of the parameter space Θ ; 2. the criterion function Q^n(θ) Q_n(θ) is twice continuously differentiable in a neighborhood N of θ0 _0; 3. the score satisfies n∇θQ^n(θ0)→(0,Σ); n\, _θ Q_n( _0)\; d\;N(0, ); 4. there exists a function H(θ)H(θ), continuous at θ0 _0, such that supθ∈‖∇θ2Q^n(θ)−H(θ)‖→ 0; _θ \; \|\,∇^2_θ Q_n(θ)-H(θ)\, \|\; p\;0; 5. the limiting Hessian H:=H(θ0)H:=H( _0) is nonsingular. Then the estimator is asymptotically normal: n(θ^n−θ0)→(0,H−1ΣH−1). n\,( θ_n- _0)\; d\;N\! (0,\,H^-1 H^-1 ). Theorem B.9 (Asymptotic Normality of IBL). Consider the IBL family pθ(y∣x)∝exp(IBLθ(x,y))p_θ(y x) (IBL_θ(x,y)) with empirical criterion as in equation 41. Assume (Xi,Yi)i=1n(X_i,Y_i)_i=1^n are i.i.d. from an admissible data law, and that the true parameter θ0 _0 is an interior point of a locally identifiable chart. For each observation Z=(X,Y)Z=(X,Y), let ℓ(θ;Z) (θ;Z) denote the per-sample loss defined in equation 40, so that Q^n(θ)=1n∑i=1nℓ(θ;Zi) Q_n(θ)= 1n _i=1^n (θ;Z_i) and Q(θ):=[ℓ(θ;Z)]Q(θ):=E[ (θ;Z)]. Suppose, in addition: 1. Score moments. s(Z):=∇θℓ(θ0;Z)s(Z):= _θ ( _0;Z) satisfies [s(Z)]=0E[s(Z)]=0, Σ:=Var(s(Z))<∞ :=Var(s(Z))<∞, and 1n∑i=1ns(Zi)⇒(0,Σ) 1 n _i=1^ns(Z_i) (0, ). 2. Derivative envelopes. There exists a neighborhood N of θ0 _0 and envelopes G1,G2G_1,G_2 with supθ∈‖∇θℓ(θ;Z)‖≤G1(Z) _θ \| _θ (θ;Z)\|≤ G_1(Z), supθ∈‖∇θ2ℓ(θ;Z)‖≤G2(Z) _θ \| _θ^2 (θ;Z)\|≤ G_2(Z), [G12]+[G2]<∞\ E[G_1^2]+E[G_2]<∞. 3. Nondegenerate curvature. H:=∇θ2Q(θ0)H:= _θ^2Q( _0) exists, is continuous at θ0 _0, and is positive definite, where Q(θ):=[Q^n(θ)]Q(θ):=E[ Q_n(θ)]. Then, under conditions of Theorem 2.7, n(θ^n−θ0)⇒(0,H−1ΣH−1). n\,( θ_n- _0)\ \ N\! (0,\ H^-1 H^-1 ). Proof. We verify the hypotheses of Theorem B.8 with Q^n Q_n as above. (i) Interior & consistency. By quotient identifiability, fix a local chart in which the population minimizer admits a unique interior representative θ0 _0. Consistency θ^n→θ0 θ_n p _0 follows from uniform M-estimation consistency for IBL (Theorem B.5). (i) C2C^2 criterion. Since IBLθIBL_θ is C2C^2 in θ, the loss ℓ(θ;Z) (θ;Z) is twice continuously differentiable in a neighborhood N of θ0 _0, and so is Q^n Q_n. (i) Score CLT. By Score moments, n∇θQ^n(θ0)=1n∑i=1ns(Zi)⇒(0,Σ). n\, _θ Q_n( _0)\;=\; 1 n _i=1^ns(Z_i)\ \ N(0, ). (iv) Hessian limit. By Derivative envelopes and dominated convergence, supθ∈‖∇θ2Q^n(θ)−∇θ2Q(θ)‖→0, _θ \| _θ^2 Q_n(θ)- _θ^2Q(θ) \| p0, so Assumption 4 of Theorem B.8 holds with H(θ):=∇θ2Q(θ)H(θ):= _θ^2Q(θ), continuous at θ0 _0. (v) Nonsingularity. By Nondegenerate curvature, H:=H(θ0)H:=H( _0) is positive definite. All assumptions of Theorem B.8 are thus verified; consequently, n(θ^n−θ0)⇒(0,H−1ΣH−1) n( θ_n- _0) (0,H^-1 H^-1). ∎ Theorem B.10 (Efficiency of IBL Estimators). Under the regularity conditions of Theorem B.9, consider the estimating function associated with the per-sample loss equation 40: ψθ(Z):=∇θℓ(θ;Z),Z=(X,Y). _θ(Z)\;:=\; _θ (θ;Z), Z=(X,Y). At any population minimizer θ⋆θ , the moment condition [ψθ⋆(Z)]=0E[ _θ (Z)]=0 holds. Define the sensitivity and variability matrices J:=[∇θψθ(Z)]|θ=θ⋆,K:=Var(ψθ⋆(Z)).J\;:=\;E\! [ _θ _θ(Z) ] |_θ=θ , K\;:=\;Var\! ( _θ (Z) ). Then the asymptotic covariance of θ^n θ_n is given by the Godambe information matrix (sandwich form): n(θ^n−θ⋆)⇒(0,J−1KJ−1). n\,( θ_n-θ )\; \;N\! (0,\ J^-1KJ^-1 ). In particular: 1. CE-only. If γc=0 _c=0 (pure cross-entropy) and the model is correctly specified and regular, then ψθ(Z) _θ(Z) coincides (up to sign) with the log-likelihood score sθ(Z)s_θ(Z). Hence J=−I(θ⋆)J=-I(θ ) and K=I(θ⋆)K=I(θ ), where I(θ⋆)I(θ ) denotes the Fisher information matrix. It follows that n(θ^n−θ⋆)⇒(0,I(θ⋆)−1), n( θ_n-θ )\; \;N\! (0,\,I(θ )^-1 ), so the estimator is asymptotically efficient, attaining the Cramér–Rao lower bound. 2. CE+DSM or DSM-only. Suppose there exists a nonsingular matrix R (constant in a neighborhood of θ⋆θ ) such that ψθ⋆(Z)=Rsθ⋆(Z)a.s., _θ (Z)\;=\;R\,s_θ (Z) .s., where sθ(Z)=∇θlogpθ(Z)s_θ(Z)= _θ p_θ(Z) denotes the parametric score in a local chart. Then J=RI(θ⋆)R⊤J=RI(θ )R and K=RI(θ⋆)R⊤K=RI(θ )R , so the sandwich covariance again reduces to I(θ⋆)−1I(θ )^-1. Hence the estimator remains asymptotically efficient. Proof. The empirical first–order condition is 0=1n∑i=1nψθ^n(Zi),ψθ(Z):=∇θℓ(θ;Z).0\;=\; 1n _i=1^n _ θ_n(Z_i), _θ(Z):= _θ (θ;Z). A mean–value expansion around the population minimizer θ⋆θ yields 0=Sn+Gn(θ^n−θ⋆),0\;=\;S_n+G_n( θ_n-θ ), where Sn:=1n∑i=1nψθ⋆(Zi),Gn:=1n∑i=1n∇θψθ~(Zi),S_n:= 1n _i=1^n _θ (Z_i), G_n:= 1n _i=1^n _θ _ θ(Z_i), for some intermediate point θ~ θ lying on the line segment between θ^n θ_n and θ⋆θ . Under the regularity conditions of Theorem B.9, we have Gn→J:=[∇θψθ⋆(Z)],nSn⇒(0,K),K:=Var(ψθ⋆(Z)).G_n pJ:=E[ _θ _θ (Z)], n\,S_n\; \;N(0,K), K:=Var( _θ (Z)). Since J is nonsingular, GnG_n is invertible with probability tending to one, and hence n(θ^n−θ⋆)=−Gn−1nSn⇒(0,J−1K(J−1)⊤). n( θ_n-θ )=-\,G_n^-1\, n\,S_n\; \;N\! (0,\;J^-1K(J^-1) ). Because here ψθ=∇θℓ(θ;⋅) _θ= _θ (θ;·), the matrix J coincides with the expected Hessian of the loss, which is symmetric. Thus the asymptotic covariance may equivalently be written as J−1KJ−1J^-1KJ^-1. (i) CE-only. When γc=0 _c=0, the per-sample loss reduces to ℓ(θ;Z)=−logpθ(Z) (θ;Z)=- p_θ(Z), so that ψθ(Z)=−sθ(Z) _θ(Z)=-s_θ(Z), with sθ(Z)=∇θlogpθ(Z)s_θ(Z)= _θ p_θ(Z) denoting the likelihood score. Under correct specification and standard likelihood regularity conditions, the information identities hold: [sθ⋆(Z)]=0,Var(sθ⋆(Z))=I(θ⋆),−[∇θsθ⋆(Z)]=I(θ⋆).E[s_θ (Z)]=0, (s_θ (Z))=I(θ ), -E[ _θs_θ (Z)]=I(θ ). Therefore, K=Var(ψθ⋆(Z))=I(θ⋆),J=[∇θψθ⋆(Z)]=I(θ⋆),K=Var( _θ (Z))=I(θ ), J=E[ _θ _θ (Z)]=I(θ ), and the asymptotic covariance simplifies to I(θ⋆)−1I(θ )^-1. Thus the estimator is asymptotically efficient, attaining the Cramér–Rao lower bound (see also Van der Vaart, 2000, Theorem 5.39). (i) CE+DSM or DSM-only under score-span. Suppose there exists a nonsingular matrix R (constant in a neighborhood of θ⋆θ ) such that ψθ⋆(Z)=Rsθ⋆(Z)a.s., _θ (Z)\;=\;R\,s_θ (Z) .s., where sθ(Z)s_θ(Z) is again the parametric score. In this case, K=Var(ψθ⋆(Z))=RI(θ⋆)R⊤,J=[∇θψθ⋆(Z)]=−RI(θ⋆).K=Var( _θ (Z))=R\,I(θ )\,R , J=E[ _θ _θ (Z)]=-R\,I(θ ). Consequently, J−1K(J−1)⊤=(−RI(θ⋆))−1(RI(θ⋆)R⊤)(−RI(θ⋆))−⊤=I(θ⋆)−1.J^-1K(J^-1) = (-RI(θ ) )^-1 (RI(θ )R ) (-RI(θ ) )^- =I(θ )^-1. Hence the sandwich covariance reduces to the Fisher information bound, and the estimator is asymptotically efficient. This corresponds to the general efficiency condition for minimum-distance or GMM estimators (see Newey & McFadden, 1994, Section 5): the condition ψθ⋆=Rsθ⋆ _θ =R\,s_θ is equivalent to their moment–span condition G′W=CG′Ω−1G W=C\,G ^-1 (Newey & McFadden, 1994, Equation 5.4), under which the Godambe information collapses to the Fisher bound. The two claims are thereby established. ∎ Appendix C Experimental Details C.1 Hardware Most experiments are conducted on a single NVIDIA L40S GPU. A small number of runs are performed on a laptop equipped with an NVIDIA GeForce RTX 2050 GPU and an Intel Core i7–12700H CPU. C.2 Standard Prediction Tasks Datasets. In the Standard Prediction Task, we use 10 OpenML datasets across diverse application domains. Details are given in Table 4. Table 4: Standard OpenML datasets used in our task. #Features denotes the number of input variables (excluding the target and ID). Name Size #Features Task type Field German Credit 1,000 20 Binary cls. Finance Adult Income 48,842 14 Binary cls. Economics COMPAS (two-years) 5,278 13 Binary cls. Law & Society Bank Marketing 45,211 16 Binary cls. Marketing Planning Relax 182 12 Binary cls. Psychology EEG Eye State 14,980 14 Binary cls. Neuroscience MAGIC Gamma Telescope 19,020 10 Binary cls. Physics Electricity 45,312 8 Binary cls. Electrical Engineering Wine Quality (Red) 1,599 11 Multiclass Chemistry Steel Plates Faults 1,941 27 Multiclass Industrial Engineering Baseline Models. For comparison, we include the following baselines: MLP, Neural Additive Model (NAM) (Agarwal et al., 2020; Kayid et al., 2020), ElasticNet, Random Forest, Stochastic Variational Gaussian Process (SVGP) (Gardner et al., 2018), Logistic Regression, Decision Tree, TabNet (Arik & Pfister, 2021), Polynomial Logistic Regression, and LightGBM (Ke et al., 2017). Table 5: Overview of baseline models in the standard prediction task Methodological Family Model Name Neural networks Standard MLP Neural Additive Model (NAM) TabNet Linear regressors ElasticNet Logistic Regression Polynomial Logistic Regression Tree-based models Random Forest Decision Tree Gradient boosting methods LightGBM Bayesian methods Stochastic Variational Gaussian Process (SVGP) Data preprocessing. For all ten datasets, we apply a consistent preprocessing strategy. Ordinal categorical variables are mapped to integer levels to preserve their inherent order. Nominal categorical variables without natural ordering are transformed using one-hot encoding. Continuous variables are standardized to zero mean and unit variance. Each dataset is randomly partitioned into train/validation/test splits with a 7:1:2 ratio. Hyperparameter Tuning Protocol. We perform hyperparameter optimization for most models using the TPE sampler from the Optuna package (Akiba et al., 2019), with 50 trials per dataset. For each model and dataset, the tuned configuration is evaluated under 8 random seeds. BL Model Hyperparameter Space. For BL(Single) and BL(Shallow), we optimize cross-entropy loss for classification. Both Adam (Kingma, 2014) and AdamW (Loshchilov & Hutter, 2017) optimizers are considered, and the better-performing variant is reported for each dataset. No data augmentation is applied. Batch sizes are chosen in a dataset-specific manner. • BL(Single): A unified setting is reported across all experiments: degreeU=[2],degreeC=[2,2,2],degreeT=[2,2],degree_U=[2], _C=[2,2,2], _T=[2,2], σparams=0.01,σλ0=0.01,σλ1=0.01,σλ2=0.01. _params=0.01, _ _0=0.01, _ _1=0.01, _ _2=0.01. Here, degreeUdegree_U, degreeCdegree_C, and degreeTdegree_T denote the polynomial degrees of the blocks that parameterize U(x,y)U(x,y), C(x,y)C(x,y), and T(x,y)T(x,y), respectively. Lists indicate both the number of blocks and each block’s degree: degreeU=[2]degree_U=[2] means a single quadratic block for U, degreeC=[2,2,2]degree_C=[2,2,2] means three quadratic constraint blocks, and degreeT=[2,2]degree_T=[2,2] means two quadratic belief blocks. σparams _params initializes coefficients of all polynomial blocks, while σλ0 _ _0, σλ1 _ _1, and σλ2 _ _2 initialize the UMP weights (λ0,λ1,λ2)( _0, _1, _2). The search grid is reported in Table 6. • BL(Shallow): We use global gradient clipping of 1.0 and an early stopping patience of 20 epochs without validation improvement. Shallow architectures with depth L≤3L≤ 3 are considered. The search grid is reported in Table 7. Baseline Model Hyperparameter Spaces. For baseline models, we also consider both Adam and AdamW for the neural network–based variants, and report results with the better-performing optimizer on each dataset. Batch sizes are tuned separately for each dataset. The detailed hyperparameter search spaces are summarized in Table LABEL:tab:tuning_space. Table 6: Hyperparameter tuning space for BL(Single) Model Parameter Search space BL(Single) learning_rate 1e−3-3, 1e−1-1 batch_size 64, 128, 256, 512 max_grad_norm 1.0, 2.0, 5.0 Table 7: Hyperparameter tuning space for BL(Shallow) Model Parameter Search space BL(Shallow) learning_rate LogUniform5e−5-5, 5e−3-3 batch_size 64, 128, 256, 512 n_layers UniformInt1, 3 n_first_layer 24, 30, 36, 40 n_middle_layer 8, 6, 4 n_last_layer 2, 4, 6 weight_decay LogUniform1e−4-4, 1e−1-1 Table 8: The hyperparameter tuning space for baseline models used in the standard prediction tasks Model Parameter Search space MLP learning_rate LogUniform1e−5-5, 1e−1-1 batch_size 32, 64, 128, 256 n_layers UniformInt2, 4 hidden_size UniformInt32, 256 weight_decay LogUniform1e−6-6, 1e−2-2 NAM learning_rate LogUniform1e−3-3, 1e−1-1 batch_size 128, 256, 512, 1024 patience UniformInt10, 30 ElasticNet (SGD) alpha LogUniform1e−4-4, 1e+2+2 l1_ratio Uniform0.0, 1.0 max_iter UniformInt100, 2000 tol LogUniform1e−6-6, 1e−2-2 fit_intercept true, false learning_rate optimal, constant, invscaling, adaptive eta0 LogUniform1e−4-4, 1e−1-1 validation_fraction Uniform0.05, 0.30 n_iter_no_change UniformInt3, 20 PolyLogistic degree 2, 3 penalty ℓ2 _2, ℓ1 _1, ”elasticnet” C LogUniform1e−3-3, 1e+2+2 l1_ratio Uniform0.1, 0.9 solver ”liblinear”, ”lbfgs”, “newton-cg”, “saga” max_iter UniformInt500, 2000 tol LogUniform1e−5-5, 1e−3-3 Logistic (ElasticNet) C LogUniform1e−3-3, 1e+2+2 l1_ratio Uniform0.0, 1.0 max_iter UniformInt100, 2000 tol LogUniform1e−6-6, 1e−2-2 fit_intercept true, false LogisticRegression solver ”liblinear”, ”lbfgs”, ”sag” C LogUniform1e−4-4, 1e+2+2 max_iter UniformInt100, 2000 tol LogUniform1e−6-6, 1e−2-2 fit_intercept True, False intercept_scaling Uniform0.1, 10.0 TabNet learning_rate LogUniform1e−4-4, 3e−2-2 batch_size 128, 256, 512, 1024 virtual_batch_size 64, 128 n_d=n_a UniformInt16, 64 n_steps UniformInt3, 7 gamma Uniform1.2, 1.7 lambda_sparse LogUniform1e−6-6, 1e−3-3 DecisionTree criterion ”gini”, ”entropy”, ”log_loss” max_depth UniformInt3, 20 min_samples_split UniformInt2, 20 min_samples_leaf UniformInt1, 10 min_weight_fraction_leaf Uniform0.0, 0.5 max_features ”sqrt”, ”log2” max_leaf_nodes UniformInt10, 1000 min_impurity_decrease Uniform0.0, 0.1 ccp_alpha Uniform0.0, 0.1 GP (SVGP) kernel rbf, matern, rational_quadratic lengthscale LogUniform0.1, 10.0 rq_alpha LogUniform0.1, 5.0 num_inducing UniformInt100, 500 learning_rate LogUniform1e−2-2, 5e−1-1 training_iters UniformInt50, 200 RandomForest n_estimators UniformInt100, 500 max_depth UniformInt3, 30 max_features ”sqrt”, ”log2” min_samples_leaf UniformInt1, 10 min_samples_split UniformInt2, 20 C.3 Interpreting BL: A Case Study C.3.1 Interpreting BL(Deep): High-Level Overview Deeper variants of BL are constructed by stacking multiple BL(Single) modules into hierarchical layers, followed by a final affine transformation. This forms a system of interacting UMPs (each of which can be viewed as an agent), where each internal block ℬB represents a single interpretable UMP. As shown in Figure 4, first-layer modules correspond to individual UMPs, while the second-layer module performs optimal coordination by aggregating or allocating their outputs. This layered structure offers a compositional interpretation of deeper BL models as systems of interacting, interpretable UMPs. C.3.2 Case Study: Additional Details Table 9: Boston Housing dataset variables and descriptions. Variable Description CRIM Per-capita crime rate by town ZN Proportion of residential land zoned for lots over 25,000 sq.ft. INDUS Proportion of non-retail business acres per town CHAS Charles River dummy variable (=1 if tract bounds river) NOX Nitric oxide concentration (parts per 10 million) RM Average number of rooms per dwelling AGE Proportion of owner-occupied units built prior to 1940 DIS Weighted distances to five Boston employment centers RAD Index of accessibility to radial highways TAX Full-value property-tax rate per $10,000 PTRATIO Pupil–teacher ratio by town B 1000(Bk−0.63)21000(B_k-0.63)^2 where BkB_k is the proportion of Black residents by town LSTAT Percentage of lower-status population MEDV Median value of owner-occupied homes in $1000s Table 10: Semantic roles of blocks in the deep BL architecture. Layer Block Representative preference Layer 1 Location-Sensitive Buyer Values river access, transport accessibility, and neighborhood amenities. Risk-Sensitive Buyer Averse to local disamenities such as pollution and environmental risk. Economic-Sensitive Buyer Sensitive to school quality and neighborhood socio-economic composition. Zoning-Contrast Buyer Responds to zoning and land-use patterns that shape local housing supply. Affordability-Preferring Buyer Strongly prefers more affordable housing and dislikes high prices. Layer 2 Integrated Location–Economic Buyer Jointly evaluates location and socio-economic attributes in an integrated way. Budget-Conflict Buyer Exhibits strong preferences for desirable locations but faces binding budget constraints. Balanced Trade-off Buyer Jointly considers multiple housing attributes in a balanced manner. Layer 3 Representative Composite Buyer Aggregates all lower-level preference components into a representative household. Table 11: Each block in the deep BL architecture is aligned with a classic preference mechanism documented in the economics literature. Layer / Block Representative reference Layer 1: Location-Sensitive Buyer Gibbons & Machin (2005) Layer 1: Risk-Sensitive Buyer Chay & Greenstone (2005) Layer 1: Economic-Sensitive Buyer Black (1999) Layer 1: Zoning-Contrast Buyer Glaeser & Gyourko (2002) Layer 1: Affordability-Preferring Buyer McFadden (1977) Layer 2: Integrated Location–Economic Buyer Bayer et al. (2007) Layer 2: Budget-Conflict Buyer Balseiro et al. (2019) Layer 2: Balanced Trade-off Buyer Rosen (1974) C.4 Prediction on High-Dimensional Inputs Datasets Description and Preprocessing. For image datasets, we use the official train/test splits of MNIST and Fashion-MNIST: Inputs are converted to single-channel images scaled to [0,1][0,1] and standardized with dataset-specific statistics. No resizing or data augmentation is applied. Training uses shuffled mini-batches of size 6464. For text datasets, we apply the following procedures: 1 Data sources and official splits. We use the official training and test splits for AG News and Yelp Review Polarity without any custom re-partitioning. Both datasets are class-balanced across labels, and we do not perform any resampling. 2 Dataset sizes. AG News: 120,000 training / 7,600 test samples with four balanced classes. Yelp Review Polarity: 560,000 training / 38,000 test samples with two balanced classes. 3 Label mapping. AG News: labels 1–4 are mapped to 0–3. Yelp Review Polarity: labels 1–2 are mapped to 0–1. 4 Text preprocessing and feature representation. All texts are lowercased and tokenized at the word level. The vocabulary is built with unigrams and bigrams, discarding words that appear fewer than two times in the training corpus. The vocabulary size is capped (AG News: 200,000; Yelp: 100,000). We compute TF–IDF weights on the training split and apply the learned weights to the test split. Dimensionality is reduced to 128 latent components using truncated singular value decomposition (SVD). Features are standardized to zero mean and unit variance and finally ℓ2 _2-normalized. We fix the random seed for reproducibility and reuse the learned preprocessing components across runs. Additional OOD Detection Results. In addition to accuracy and AUROC, we also report AUPR and FPR@95 for both image and text datasets; the results are shown in Table 12. On image datasets, BL (depth=1) achieves the best overall balance: it ranks first on Fashion-MNIST AUPR and second on Fashion-MNIST FPR@95. On MNIST, it is second in AUPR but underperforms in FPR@95 compared with E-MLP (depth=2). These results suggest that BL yields separable score distributions, particularly on Fashion-MNIST, although its 95% FPR threshold admits more OOD samples than E-MLP at the same recall. On text datasets, OOD detection performance is dataset-dependent: E-MLP performs better on AG News, whereas BL achieves stronger OOD performance on Yelp. Table 12: OOD AUPR and FPR@95 (%) on image and text datasets. BL and E-MLP are evaluated at depths 1–3 with matched parameter counts, both without skip connections. Top-two per column are blue and red. Model MNIST Fashion-MNIST AUPR FPR@95 AUPR FPR@95 E-MLP (depth=1) 89.37±\;±\;1.52 35.57±\;±\;5.87 91.35±\;±\;1.25 28.24±\;±\;4.37 BL (depth=1) 91.57±\;±\;2.39 47.81±\;±\;11.29 91.79±\;±\;0.90 38.86±\;±\;2.57 E-MLP (depth=2) 91.52±\;±\;1.27 28.89±\;±\;2.85 86.19±\;±\;2.27 47.72±\;±\;4.79 BL (depth=2) 91.20±\;±\;1.22 52.71±\;±\;18.66 89.30±\;±\;2.47 42.65±\;±\;9.53 E-MLP (depth=3) 90.04±\;±\;1.89 31.92±\;±\;5.76 84.30±\;±\;1.50 54.49±\;±\;2.74 BL (depth=3) 92.36±\;±\;2.03 32.32±\;±\;5.76 88.41±\;±\;4.04 41.19±\;±\;13.36 Model AG News Yelp AUPR FPR@95 AUPR FPR@95 E-MLP (depth=1) 44.52±\;±\;15.10 33.82±\;±\;4.89 3.31±\;±\;1.60 54.21±\;±\;2.11 BL (depth=1) 18.68±\;±\;16.48 42.03±\;±\;5.99 12.70±\;±\;2.29 40.95±\;±\;1.56 E-MLP (depth=2) 31.48±\;±\;23.94 40.20±\;±\;9.82 1.47±\;±\;2.31 57.80±\;±\;4.88 BL (depth=2) 10.76±\;±\;15.94 53.71±\;±\;9.68 6.73±\;±\;2.61 46.54±\;±\;1.86 E-MLP (depth=3) 51.24±\;±\;9.13 32.96±\;±\;3.23 3.14±\;±\;2.22 55.24±\;±\;5.33 BL (depth=3) 16.99±\;±\;17.12 45.24±\;±\;6.11 10.96±\;±\;1.10 42.27±\;±\;1.94 Number of Parameters. To ensure a fair comparison between E-MLP and BL, we match the number of trainable parameters as closely as possible for models with the same depth (see Table 13). Running Time. To evaluate computational cost, we compare the training time of BL and Energy-based MLP across image and text datasets (Table 3). Under comparable parameter budgets, BL generally requires slightly higher training time than E-MLP across datasets. In particular, BL is moderately slower on image datasets and AG News, while exhibiting comparable running time on Yelp. Calibration We report ECE and NLL metrics to assess calibration quality, and the results are presented in Table 2. On image datasets, BL provides substantially better calibration, with BL models occupying the top two positions in each column. On text datasets, calibration performance is broadly comparable, with BL showing slightly lower NLL on Yelp. Overall, these results indicate that BL delivers strong predictive performance together with reliable probability estimates. Table 13: Number of trainable parameters for E-MLP and BL models across high-dimension datasets. Dataset Model # Parameters MNIST & FashionMNIST E-MLP (depth=1) 203,530 BL (depth=1) 208,384 E-MLP (depth=2) 235,146 BL (depth=2) 219,264 E-MLP (depth=3) 238,314 BL (depth=3) 221,684 AGNews E-MLP (depth=1) 136,196 BL (depth=1) 149,720 E-MLP (depth=2) 386,284 BL (depth=2) 397,568 E-MLP (depth=3) 230,788 BL (depth=3) 224,128 Yelp E-MLP (depth=1) 134,146 BL (depth=1) 148,960 E-MLP (depth=2) 385,770 BL (depth=2) 397,312 E-MLP (depth=3) 230,530 BL (depth=3) 224,000 C.5 Case Study: Estimation Results of BL on the Boston Housing Dataset Table 14: Estimated UMP block parameters learned by the BL model (layer = [2, 1]) on the Boston Housing dataset. For each block, U denotes the Utility component, C the Inequality-Constraint component, and T the Equality-Constraint component. Block 11 Block 12 Variable U11U_11 C11C_11 T11T_11 U12U_12 C12C_12 T12T_12 λ 1.003 0.997 0.999 0.997 1.003 1.000 per capita crime rate (CRIM) 0.21 0.14 0.03 0.12 0.09 0.25 residential land proportion (ZN) 0.23 -0.04 -0.27 0.25 0.00 0.09 non-retail business acreage (INDUS) -0.06 0.21 0.25 0.16 0.22 0.27 Charles River dummy (CHAS) 0.25 0.04 -0.24 -0.12 -0.20 -0.23 nitric oxide concentration (NOX) -0.06 -0.13 0.21 0.16 0.02 -0.28 average rooms per dwelling (RM) 0.06 0.07 0.05 0.05 -0.19 -0.22 proportion of older units (AGE) -0.13 -0.12 -0.09 0.14 0.08 -0.18 distance to employment centres (DIS) 0.16 -0.03 0.17 -0.17 -0.09 0.11 radial highway accessibility (RAD) 0.24 -0.11 0.04 -0.28 0.09 0.10 property tax rate (TAX) -0.20 0.18 0.22 -0.11 -0.06 0.23 low-income population (LSTAT) 0.05 -0.12 -0.09 0.23 -0.16 -0.19 median home value (MEDV) 0.21 -0.08 0.07 0.08 -0.17 0.15 Constant term (C) 0.03 -0.17 -0.07 0.11 -0.16 -0.12 Block 21 Variable U21U_21 C21C_21 T21T_21 λ 1.000 1.003 0.999 Block 11 output (b1,1b_1,1) 0.428 -0.551 0.147 Block 12 output (b1,2b_1,2) -0.168 -0.356 -0.178 Constant term (C) 0.406 0.219 0.421 Table 15: Estimated UMP parameters for the Layer 1 blocks of the BL model (layer = [5, 3, 1]) trained on the Boston Housing dataset. Here, U denotes the Utility component, C the Inequality-Constraint component, and T the Equality-Constraint component. Variable U11U_11 U12U_12 U13U_13 U14U_14 U15U_15 λ 1.000 0.998 1.003 1.002 1.000 per capita crime rate (CRIM) 0.21 0.12 0.17 -0.09 0.06 residential land proportion (ZN) 0.23 0.25 -0.07 -0.22 -0.16 non-retail business acreage (INDUS) -0.06 0.16 0.16 0.23 -0.14 Charles River dummy (CHAS) 0.25 -0.12 -0.22 -0.05 -0.01 nitric oxide concentration (NOX) -0.06 0.16 -0.14 0.24 0.16 average rooms per dwelling (RM) 0.05 0.05 0.08 0.09 -0.07 proportion of older units (AGE) -0.13 0.14 0.06 -0.23 -0.15 distance to employment centres (DIS) 0.16 -0.17 -0.07 0.19 -0.10 radial highway accessibility (RAD) 0.24 -0.27 0.17 -0.08 -0.20 property tax rate (TAX) -0.20 -0.11 0.19 -0.11 0.10 low-income population (LSTAT) 0.05 0.23 -0.15 -0.28 -0.26 median home value (MEDV) 0.21 0.08 0.25 0.08 0.06 Constant term (C) 0.03 0.12 -0.09 -0.06 0.15 C11C_11 C12C_12 C13C_13 C14C_14 C15C_15 λ 0.999 1.001 1.000 0.997 1.002 per capita crime rate (CRIM) 0.13 0.09 -0.10 0.11 0.05 residential land proportion (ZN) -0.04 -0.01 -0.27 -0.23 -0.10 non-retail business acreage (INDUS) 0.21 0.22 -0.16 0.20 0.15 Charles River dummy (CHAS) 0.04 -0.19 0.07 -0.20 0.15 nitric oxide concentration (NOX) -0.13 0.02 -0.04 -0.05 0.11 average rooms per dwelling (RM) 0.07 -0.19 -0.20 0.06 -0.05 proportion of older units (AGE) -0.13 0.08 0.01 0.14 -0.07 distance to employment centres (DIS) -0.03 -0.09 -0.19 0.22 0.03 radial highway accessibility (RAD) -0.11 0.08 -0.23 0.25 -0.05 property tax rate (TAX) 0.18 -0.06 -0.15 -0.22 -0.08 low-income population (LSTAT) -0.13 -0.17 -0.18 -0.12 0.24 median home value (MEDV) -0.08 -0.17 0.28 -0.03 -0.03 Constant term (C) -0.17 -0.16 0.05 -0.21 -0.06 T11T_11 T12T_12 T13T_13 T14T_14 T15T_15 λ 0.999 1.002 0.999 1.004 1.001 per capita crime rate (CRIM) 0.03 0.25 0.08 0.25 0.00 residential land proportion (ZN) -0.27 0.10 -0.26 -0.20 -0.02 non-retail business acreage (INDUS) 0.25 0.26 -0.18 0.15 0.07 Charles River dummy (CHAS) -0.23 -0.23 -0.09 0.10 0.08 nitric oxide concentration (NOX) 0.21 -0.28 0.04 0.09 -0.25 average rooms per dwelling (RM) 0.05 -0.21 -0.24 -0.15 -0.10 proportion of older units (AGE) -0.09 -0.19 -0.12 0.26 0.24 distance to employment centres (DIS) 0.17 0.12 -0.17 0.06 0.10 radial highway accessibility (RAD) 0.04 0.10 0.00 0.04 -0.01 property tax rate (TAX) 0.22 0.23 -0.10 -0.24 -0.17 low-income population (LSTAT) -0.09 -0.19 -0.19 -0.04 -0.25 median home value (MEDV) 0.08 0.15 -0.19 -0.13 -0.09 Constant term (C) -0.07 -0.11 -0.16 0.24 0.10 Table 16: Layer 2 and Layer 3 UMP parameters (U, C, T) for Blocks in the BL model (layer = [5, 3, 1]). Block 21 Block 22 Block 23 Variable U21U_21 C21C_21 T21T_21 U22U_22 C22C_22 T22T_22 U23U_23 C23C_23 T23T_23 λ 1.000 1.000 1.000 0.999 1.003 1.002 1.001 1.002 0.999 Block 11 output (b1,1b_1,1) 0.28 0.06 -0.20 -0.31 0.24 0.18 -0.29 -0.08 0.22 Block 12 output (b1,2b_1,2) 0.21 -0.11 -0.09 -0.44 0.12 -0.22 0.15 -0.22 0.20 Block 13 output (b1,3b_1,3) -0.40 0.18 -0.44 -0.36 -0.01 -0.09 -0.13 -0.14 0.32 Block 14 output (b1,4b_1,4) -0.27 -0.17 0.30 0.33 -0.34 -0.26 0.28 -0.42 -0.34 Block 15 output (b1,5b_1,5) -0.07 -0.29 0.34 0.22 -0.38 -0.08 -0.14 0.25 0.32 Constant term (C) 0.43 0.33 0.16 0.38 -0.42 -0.32 -0.33 -0.31 -0.21 Variable U31U_31 C31C_31 T31T_31 λ 1.002 0.998 1.000 Block 21 output (b2,1b_2,1) 0.21 -0.13 0.36 Block 22 output (b2,2b_2,2) 0.54 -0.48 0.43 Block 23 output (b2,3b_2,3) -0.08 0.28 0.55 Constant term (C) -0.01 -0.58 -0.14