A Complete Decomposition of KL Error using Refined Information and Mode Interaction Selection
Authors: James Enouen
Mahito Sugiyama
Venue: Accepted by TMLR ()

Abstract: The log-linear model has received a significant amount of theoretical attention in previous decades and remains the fundamental tool used for learning probability distributions over discrete variables. Despite its large popularity in statistical mechanics and high-dimensional statistics, the vast majority of related energy-based models only focus on the two-variable relationships, such as Boltzmann machines and Markov graphical models. Although these approaches have easier-to-solve structure learning problems and easier-to-optimize parametric distributions, they often ignore the rich structure which exists in the higher- order interactions between different variables. Using more recent tools from the field of information geometry, we revisit the classical formulation of the log-linear model with a focus on higher-order mode interactions, going beyond the 1-body modes of independent distributions and the 2-body modes of Boltzmann distributions. This perspective allows us to define a complete decomposition of the KL error. This then motivates the formulation of a sparse selection problem over the set of possible mode interactions. In the same way as sparse graph selection allows for better generalization, we find that our learned distributions are able to more efficiently use the finite amount of data which is available in practice. We develop an algorithm called MAHGenTa which leverages a novel Monte-Carlo sampling technique for energy-based models alongside a greedy heuristic for incorporating statistical robustness. On both synthetic and real-world datasets, we demonstrate our algorithm’s effectiveness in maximizing the log-likelihood for the generative task and also the ease of adaptability to the discriminative task of classification.

Published in Transactions on Machine Learning Research (04/2026)
A Complete Decomposition of KL Error using Refined Information and Mode Interaction Selection
James Enouen"
enouen@usc.eddv
Deportment of Computer Science University of Southern California Los Angeles, CA 90089-2905, USA
Mahito Sugiyama
mahito@nii..cc.p
Department of Computer Science
National Institute of Informatics, SOKENDAI Chiyoda Cily, Tokyo 101-8430, Japan
Reviewed on OpenReview: https: oppereevee.. net/ pormf id= DeC3I15C38
Abstract
The log-linear model has received a significant amount of theoretical attention in previous decades and remains the fundamental tool used for learning probability distributions over discrete variables. Despite its large popularity in statistical mechanics and high- dimensional statistics, the majority of related energy-based models only focus on the two-variable relationships, such as Boltzmann machines and Markov graphical models. Although these approaches have casier-to-solve structure learning problems and casier-to- optimize parametric distributions, they often ignore the rich structure which exists in the higher-order interactions between different variables. Using more recent tools from the field of information geometry, we revisit the classical formulation of the log-lincar model with a focus on higher-order mode interactions, going beyond the 1-body modes of independent distributions and the 2-body modes of Boltzmann distributions. This perspective allows us to define a complete decomposition of the KL error. This then motivates the formulation of sparse selection problem over the set of possible mode interactions. In the same way as sparse graph selection allows for better generalization, distributions are able to more efficiently use the finite amount of data which is available in practice. We develop an algorithm called MAHGen'Ta which leverages a novel Monte-Carlo sampling technique for energy-based models alongside a greedy heuristic for incorporating statistical robustness. On both synthetic and real-world datasets, we demonstrate our algorithm's effectiveness in maximizing the log- likeelikooood for the generative task and also the ease of adaptability to the discrininative task of classification.
find that our learned
1 Introduction
Distribution learning is a fundamental task in machine learning and statistics, with applications ranging from generative tasks like density estimation and generative modeling all the way to discriminative regression tasks and unsupervised clustering tasks. This fundamental problem remains at the heart. of many super- vised decision problems and as cornerstone of unsupervised learning and knowledge discovery. Probability distributions parametrized by exponential families are representative class widely chosen for distribution modeling. Although already covering wide variety of classical distributions for continuous variables (Gaus- sian, exponential, gamma, etc.), for finite and discrete variables, the hicrarchical log- lin-lincar model completely
Work done while at National Institute of Informatics.

---

Published in Transactions on Machine Learning Research (04/2026)
describes all positive distributions over the space. Also called energy-based model, it has remained the de facto choice of model for learning over discrete feature space for decades and has amassed considerable attention over the years (Ackley et al.]
1985; Sejnowski 1986 Lee et al., 2006; WWinwright et al., 2006
Shpitser et al. 2013 Van Haaren & Davis 2012: Lowd & Davis 2010 Nyman et al. 2014 Højsgaard
2004).
Despite this significant amount of work to date, the vast majority of existing approaches only deal with bivariate correlations or two-body interactions, prototypical examples being Boltzmann machines and Markov graphical models 4Ackley et al. 1985 Dempster 1972 Buhl 1993). Although this assumption is natural to fore onto continuous variables by making the simplifying Gaussian assumption, this restriction is too severe for most real-world data distributions and is not necessary for dealing with the case of discrete variables. Many existing amendments to these approaches like maximal cliques, chordal graphs, and stratified graphical models (Shpitser et al. 2013 Nyman et al., 2014 Hojsgaard 2004] are still only graph- based or two boo- structure approximations, remaining limited in their ability to describe the underlying higher-order structures which can exist within data. In this work, we instead offer a more unified perspective which further includes the hypergraphical structure encoded by higher-order interactions.
By replacing the (two-dimensional) edge graph between features with the higher-order hypergraph, we in- troduce a structure lcarning problem which is seemingly even more challenging than the usual graphical approaches. Despite this greater complexity a priori, using recent. developments in the field of information geometry (Ghalamkari et al. 2023: Sugiyama et al. 2018) allows us to construct a complete decomposition of a distribution's information content, instead actually providing a greater fine-grained understanding into the structure of a probability distribution. By associating the hypergraph of the hierarchical log-linear model of discrete distribution with the partially-ordered set (poset) of mode interactions in the corresponding probability tensor, are able to devise a higher-order definition of non-negative information which provides complete decomposition of the KL error for a given probability distribution.
Altogether, we demonstrate that this alternate perspective is theoretically well-supported, allowing us to de- fine more fine-grained measurements of the information between variables and opening up new opportunities for the study of higher-order structure. Additionally, we provide practically useful learning algorithms for both the combinatorial structure and the parametric value of the distribution. We summarize our contribu- tions as the following:
We first define the 'refined information' of a set of two or more variables, generalizing the mutual information of a set of two variables in a way which always returns positive quantity measuring the information content. We show that this yields a complete dccomposition of the KL error with applications in structure discovery.
We provide the first theoretical underpinnings for the better generalization properties of higher-order Boltzmann machines via the problem of 'mode interaction selection', showing how to yield better sample complexity for real-world scenarios with finite datasets. We then show how the combinatorially large space of all possible interaction hierarchies can be effectively tackled by greedy approach.
We finally develop our model called the (MAHGenTa) which implements GPU-based gradient descent algorithm to efficiently learn the hier- archical log-linear model on both synthetic and real-world datasets. We further demonstrate how such energy-based models trained to achieve good generative performance will have automatically emergent capabilities in discriminative tasks like classification, paralleling the wide success of generative pretraining.
Hierarchy for Generating Tabular data
2 Background
Before introducing the main task of distribution learning, we first review the modern approaches from feature selection which provide the high-dimensional intuition for our statistical arguments and the information geometric approaches to tensor decomposition which inspires the main distributionntensor correspondence.
2

---

Published in Transactions on Machine Learning Research (04/2026)
2.1 Feature Selection and Feature Interaction Selection
Feature selection (FS) has long been staple of machine learning for dealing with high-dimensional data, prescreening large number of features to remove both irrelevant and redundant features from the input before the training a predictive model. This provides many benefits like reducing overfitting, faster training, and better understanding of the data structure. Historical approaches determine relevant, irrelevant, and redundant features by understanding the mutual information between the inputs and the target variable (Shannon 1948 McGill 1954), whereas modern fcature selection approaches concern themselves with giving the proper credence to higher-order interactions and correlations during selection, generally called ure- interaction-awan feature selection' (Zeng et al. 2015 Nakariyakul 2018 2015
Chen et al. 2015 Bennasar et al.
In recent years, however, there has been a parallel interest in feature interactions via a more general problem than feature selection called feature interaction selection (FTS) (Fan et al.] 2016 Sugiyama & Borgwardi, 2019 Enouen & Liu, 2022. Lyu et al. 2023). This procedure further specifies how feature combinations are allowed to interact the final model. For example, in a random forest model, feature interaction selection would dictate how each different. tree can only use specific interaction subset,, rather than any possible subset out of the selected features. Although FIS raises the combinatorial complexity of the search problem from all possible subsets to all possible collections of subsets, the finer-grained structure further amplifies the same typical benefits of feature selection: reduced overfitting, faster training, and better understanding. In this work, we will apply this same mothodology to the generative task of distribution learning instead of to the discriminative task of classification.
2.2 Non-Negative Tensor Decomposition
Recent works in tensor decomposition have been able to avoid the optimization difficulties of typical low- rank decompositions by instead focusing on non-negative tensors and replacing the squared error loss with the KL divergence error [Aswani, 2016 Sugiyama et al.) 2018 Ghalamkari et al. 2023 Previous works assuriing that a full decomposition is leasible (Sugiyama et al. 2018) or that a Lensor's modes can be partitioned into independent components [Aswani] (0166) have recently been replaced by specific control over the 'mdny-body interactions' within the tensor (Ghalamkari al. 2023). Following this notation, using one- body interactions corresponds to independent distributions and using (wo-body interactions corresponds to Boltzmann machines.
Following this correspondence between the problems of approximating non-negative probability tensor and learning a discrete disuribution via minimizing the same KI-divergence objective, we adopt the language 'mode interaction' to describe variables interacting during generative modeling, instead the more popular 'feature interaction' of Section [2.1 which focuses on the discriminative cases. this work, we leverage the consequent connections to information geometry (Amari) 2016)) to develop procedure of mode interaction selection (MTS) which parallels the higher-order FTS selection problem, but applied to the variable interactions between tensor modes.
2.3 Distribution Learning
The Log-Linear Model The log-linear model has always been a fundamental tool of statistics, with its use dating as far back as the historical works of Fisher (Fisher) 1934). In continuous variables, focusing on exponential familics of distributions may limit us to certain types distributions (Ganssians, Poissons, etc.); however, in the case of finite variables, the log-linear model has no such limitations. Accordingly, the log- linear model has been staple of describing any categorical distribution, being the focus of many Bayesian optimal inference frameworks as well as central tool of statistical physics.
In the previous decade, this method received a significant amount of allention alongside the wave of spar- sity methods enabled by LASSO [Tibshiran] 1996). This primarily led to an abundance of graph-based approaches which perform selection over the pairs of 2-body correlators between fcatures, such as Markov
graphs with L1 regularization {Lee et al., (2006) Wainwright et al., 20066] Lowd & Davis 2010 Van Haaren & Davis 2012). This was later followed up with specific types of graphical assumptions which further
3

---

Published in Transactions on Machine Learning Research (04/2026)
Massam et al, 2009). These graphical models have allowed for more fine-grained control over the structure of the learned distribution compared with previous approaches, receiving the benefits of sparsity which allows for easier generalization with fewer data samples. However, in the same way that FS does not allow for the full, higher-order control of FIS, graphical selection does not give the full, higher-order control of MIS.
simplify the learning problem (Shpitser et al. 2013 Nyman et al. 2014 Højsgaard 2004
Higher Order Boltzmann The 'fully-visible higher-order Boltzmann machine', extending 2-body corre- lations to all higher order interactions, was formulated long before it could be made practical (Sejnowski, Early works focus on binary variables with very small event spaces (Amari 2001 Nakahara & Amari, 2002 Ganmor et al.) 2011), mainly interested in biological applications like neuronal activity and protein interactions. These older works mainly ignored the computational issues associated with scaling to more serious sizes and accordingly remained limited in their application. The most closely related works to ours are the two works that have extended the sparse graphical modeling formulation to higher-order interac- tions, (Schmidt & Murphy 2010) and (Min et al.) 2014). Although not extending beyond binary variables, these works attempt scalability by formulating the hierarchical structure learning problem and trying to overcome the computational challenges of scaling higher-order Boltzmann machines to learning real-world data distributions.
1986
Although these two works make strides towards defining the problem and attempting to scale beyond syn- thetic datasets, they still apply only to binary variables with many challenges remaining even for medium- scale datasets. Moreover, a deeper theoretical understanding of the statistical benefits had yet to be de- veloped. We push further what is possible in practice by leveraging a theoretical grounding from tools in information geometry and efficient GPU-based training methods for modern applications. This progress, however, must be contrasted against the significant gap in capability when comparing to likelihood-free neural approaches.
SOTA Generative Models It is reminded that most recent work in distribution learning has shifted entirely away from having direct control over the distribution at all. Instead, enabling and later enabled by deep learning, recent models have been developed to instead reshape large set of hidden or latent variables towards the distribution of the data, following the trail blazed by RBMs and DBMs (Hinton & Salakhutdinov 2006 Salakhuudiinov & Hinton 2009) and extending into the modern day VAEs, GANs, and dillusion models (Kingma & Welling 2014 Goodfellow et al.) 2014 Ho el al. 2020
In contrast to these methods, this work learns the hierarchical linelr model on the data features, directly providing predictions of likelihood. For tabular datasets with interpretable variables, there is the great benefit of having interpretable insights into the data structure compared to what is available from latent approaches. Furthermore, is imagined that this work's revisiting of the theoretical foundations for the statistical generalization properties of 'visible-only' generative modeling may ultimately have downstream effects on bettering our understanding of the generalization of generative modeling in more general cases including latent variables.
3 Refined Information
Notation Consider distributions over d E N variable dimensions, with each feature k E [d] = (1,..., d
having 4 E N discrete (and disjoint) possibilities called events. We write indices as i E [41] x ... x [1a] and index distributions as p(i). When we deal with subsets of the features S C [d]] and their conjugates -S : [d]- - - (where we use + and - as shorthand for set union and set difference), we write the subindices is E Is = @RES( and the marginal distributions (is) == Ej_seLs p(is.j-s). We write the powerset as P((d) = (S _ [d] = 10 and collection of interaction subsets as I S P([d),, or equally I P(P((d))) To reiterate, we will write: A E [ E P((4), I E P(P((d)), hence also k E S.S
Every distribution will be considered simultaneously as a discrete distribution and as a finite tensor, meaning that we interpret the tensor product on distributions as (p^ & p (i == := ((A)) -p (in). We write u to represent the uniform distribution and we will later write Pern and Pval represent the empirical training and validation distributions.
4

---

Published in Transactions on Machine Learning Research (04/2026)
(a)
S
[1,2,3}
(b) (() bits
(c) I(s)
-1
bits
(d)
J(S)
1 bit
{1.2)
{1,3}
{
2 bits
bits
2 bits
bits
0
bits
bits
bits
bits
0 bits
X
X
(1)
(2}
[3]
bit
bit
1 bit
1 bit
bit
bit
0 bits
bits
0 bits
0
0 bits
0 bits
0 bits
Figure 1: Three types of higher-order information (Hs))), I(ps), (ps)) defined for each S S [3].
3.1 Information Theory
Information theory was born out of the fundamental contributions of Shannon (Shannon 1948) defining the entropy of variable and the mutual information (MI) between two variables. Shortly after, extension to three or more variables was constructed with the multiple mutual information (MMT) (McGill 1954). We write the definitions of entropy, (ps), and MI/MMI, (pss) as:
H(ps) = - Σ (og(p))))sssss)
(1)
I(ps) = Σ
(-1) HPT)),
(2)
TCS
J(ps) = Σ D
(3)
TES
We additionally write the definition of J(ps), which is closely related to the original MI and MMI via Jps)) = (-1))s)((ps) for |S| 2 2. As will later see in the definitions to follow, this simple rephrasing is born out of our information geometric viewpoint rather than the original communication theory viewpoint.
In Figure D we can see their different properties on a simple distribution over three binary variables. In particular, take X1 and X₂ to be Bernoulli and the third X3 to be the XOR of the first two, X3 = X1 @ X2. All three variables are symmetric and any one/ two of them is indistinguishable from a single/ pair of random bits. Accordingly, there is no mutual information between any of the pairs. However, despite one or two of the three variables seeming totally random, knowing all three variables simultaneously completely different from the case of three independent variables. Accordingly, this distribution exemplifies the need to discuss 'purely third-order' types of information.
Unfortunately, (ps)) and (ps) may return positive or negative values whenever |S] 2 3, diminishing the ability to interpret MMI as 'information content' as is possible for the mutual information ion between two variables. This inspires our new definition in Equation 3 of Section 3.3 which generalizes mutual information but still returns a non-negative quantity for higher-order interaction information.
3.2 Information Geometry
It is fairly well known that we can model any discrete distribution by an exponential family:
go()) = exp Σ sgi 05
for some continuous parameters 8 Rs It is lesser known that when equipped with KL divergence, this results in a Riemannian manifold called a statistical manifold. Due to the exponential structure, this further results in a dually flat or Bregman flat manifold, allowing for an even richer theory to be developed over the space (Rao) 1992 Amari & Nagaoka]
20000 Amari 2016 Niclsen et al. 2017). The natural parameters and the dual expectation 17 parameters mi = Ej>ge() [[ig js) % (is) are 'orthogonal' to each other
in particular sense, obeying the Legendre transform and Bregman duality. Details in Appendix A
5

---

Published in Transactions on Machine Learning Research (04/2026)
Of particular importance for work are the Pythagorean theorem and the projection theorem Nagaoka & Amari 1982), which imply the uniqueness of the distribution projected onto a flat submanifold of the distribution space as well as the convexity of the corresponding optimization problem. See Theorems 2 and
in Appendix A for details. In particular, for a fixed p and a fixed collection of interactions I _ P((d), the best forward DKL approximation of p within the submanifold Mz allows for unique projection pr:
Mz = {00 go exp { SET }
PI = (p) = argmin { (KLL
(4)
qEMı
where MMz denotes the projection onto that submanifold. particular, the solution to our optimization problem hence guaranteed have a unique solution, and this solution respects the dually flat manifold.
3.3 Definition of Refined Information
In order to achieve our goal of defining fine-grained and high-order definition of information content, we use this information geometry lens to define a sequence of projected distributions, and then use the distance (divergence) between each of the distributions in that sequence of projections the definition of information content. Because mutual information can also be defined in this way, our definition gives a natural extension to higher-order and non-negative information content which we call refined information.
Let us say that a collection I is hierarchical if it downwards closed with respect to subsets (S E I and T S S T I We can define chain of collections such that:
(0) C To I & ... I P((d)).
We will say that a chain complete whenever Io = {0} and Zr = P(dd) and is hierarchical whenever each I, is hierarchical. We will hereafter restrict our attention to complete and hierarchical chains. Moreover, we will almost always considered mazimally-reffied chains, which equivalent to saying IN (I+-1 + S,) for some subset S,, for all choices of t. Any shorter chain can be extended into maximal chain by further refining the sequence.
We use our chain of hierarchical collections in conjunction with the information geometry projection to construct in parallel a chain of distributions P......,Pi. -P,,- These can be seen as the repeated projection onto submanifolds which slowly approach our target distribution: IT, (p),..., TTy (p). Accordingly, we know that each drop in divergence defines unique and positive quantity which will call the refined information from I to J: RIT-JJ (p) DKL. (p)iPT) (P:( (KI(P:P)). From this definition, is clear that:
DEL u) -'r 11z-1-1((1))
where we recall that the uniform distribution u is the null model the space of finite distributions (all 0 coordinates are zero).
Definition 1. The refined information of S at I is:
RII,s(p) RII (I+)s())) DKL (PL4s:::)))
(5)
This leads to a full decomposition of the KL error as:
DKL u) = E RII.,...
(6)
Accordingly, after fixing a chain, this formula altributes each positive drop in KL error to a single interaction set S. Since the goal of distribution learning is to reduce the KL divergence to zero, this decomposition allows for extremely fino-grained control by directly corresponding cach effective parameter @ we may choose inchude with a decrease error. We discuss the implications of this for generalization performance in the coming Section 4.1

---

Published in Transactions on Machine Learning Research (04/2026)
4 MAHGenTa
Here we introduce the Mode-Attributing Hicrarchy for Generating Tabular data (MAHGenTa) to tackle this donbly-exponential combinatorial problem and efficiently learn an arbitrary probability distribution. Our procedure consists of two major components: (1) a mode interaction selection algorithm in conjunction with an early stopping procedure to guarantee a low gap between the train and test performances; and (2) an efficient gradient descent training algorithm which overcomes the challenges of the normalizing constant with energy-based modeling and a GPU-enabled pytorch implementation which extends existing Gibbs samplers to higher-order tensors.
argmin
(KII(POaI: ) where % argmin (DKz (Prri ())mmm))))) I,0z
(7)
W.MMz
We write our learning objective as bilevel optimization problem in Equation Further details justifying this choice under the thcoretical framework decomposing the KL error are provided in the Appendix.
4.1 Mode Interaction Selection
Even when ignoring the difficulties associated with learning the continuous 0 parameters in Section 4.2 there are major practical challenges in finding good collection I from the combinatorially explosive set of available choices. Particularly, the question of how to select good collection of mode interactions which accurately describe the distribution without overfitting to the training set. We must leverage an appropriate hcuristic for selecting interacting modes amongst the 2 choices of interaction and 2 choices of final collection.
We follow similar greedy heuristics as have been explored in previous literature based on the strong or weak 'heredity' assumption for choosing pairs (Peixoto 1987 Bien al.] 2013). Generalizing this to higher order interactions allows to only consider a polynomial number of candidate interactions. In particular, we will explore starting from the smallest S, only considering S if its heredity score, (SS) is greater than the threshold = 30%. Further discussion of heredity in Appendix B.4
w(S) = [T: S = U (t} for some i) n (TT: T has already been selected]]
(8)
Based on our thcoretical developments Equation 6 we would like to add cach Ø parameter which corresponds to the greatest amount of refined information, continuing until our validation KL error stops decreasing alongside our training KL error. Early stopping in this way is just because every parameter of the log-linear model is an effective parameter, and sequential projections along the statistical manifold curve.
will cause our model to obey the classical mnderfitiig-overffti
Unfortunately, exactly computing the refined information is difficult because for degree three and higher as there is no closed form available and one must resort to the continuous optimization approaches leveraged herein. Instcad, we must a priori choose some heuristic measurement which corresponds with high refined information gain from including a specific mode interaction within the log-linea, model. Accordingly, we use the absolute value of Jg as introduced in Section 3.1
easy-to-compute alternative to RIs. We present our search and learning algorithm in Algorithu | which continuously alternates between adding new mode interactions to the model and using gradient descent to train with the new parameters. Each subroutine is presented in full detail in the appendix.
4.2 Gradient Descent Learning
As mentioned in previous sections, the optimization of {0S) is always a convex problem. Accordingly, for a small enough learning rate, we can always guarantee the convergence of the gradient descent algorithm. Nevertheless, we find there are still multiple challenges to overcome for fast. training of log-linear models when attempting to scale to real-world datasets.
We first recall the gradient of a linelinear model when optimizing for forward KL divergence is 8s DKL. (pumi go) u bh when evaluated on the empirical training distribution Pim- However,
7

---

Published in Transactions on Machine Learning Research (04/2026)
Algorithm 1: MAHGenTa Algorithm
MARGENTA(T, a, K,T)
Ertbest oc, I 4- (0), O (0}
3
while Error < ErTbest do
4
Errbest Error(e)
5
J -AAATAAALLLLLLLLTBAAATTTOSSSSS
I 1 for S E K do
0$ - RER's
8
θ θ ou{
9
10
θ
GRADIENDEESCET((O...].
return
because the parameter space of the hicrarchical model is invariant under constant shifts along any paramcter tensor (so long another parameter absorbs the negative shift), will restrict each 0 such that its sum across every mode/file /fiber is equal to zero. Practically, this leads to the use of the purified gradient:
as (Pura 90)] = (-)5s(-
T 2 S-T
(9)
To facilitate the modern applicability of our algorithm, we implement GPU-based gradient descent. training algorithm in Pytorch (Paszke al. 2019). The major challenge of any gradient implementation for energy- based models is the calculation of the partition function normalizing constant. Even after implementation tricks and virtualization of the tensor, exact computation is simply too slow to handle the billions of events which are possible even medium-dimensionall tabular datasets, and we must resort to a new variant of the classical Gibbs sampling approach |Geman & Geman 1984 Gelfand 2000) in conjunction with the annealed importance sampling technique (Neal, 2001]. We find that the compound use of several tricks like this is critical for achieving the fastest implementation of gradient descent for higher-order energy based models, provide detailed explanation of each trick in Appendix B.5
5
Experiments
To first get a glimpse into the theoretical properties of refined information and the sample complexity of MAHGenTa we generate a suite synthetic distributions by choosing random Ø and sampling data from the induced distribution, details in the appendix. We demonstrato the impact of structure learning and the value of refined information in this setting where we have full control over the structure in the ground truth, We next apply method to three real-world distributions from UCT machine learning datasets. The three datasets used are shown in Table with their numbers of samples, features, and total possible events.
5.1 Synthetic Results
In Figure 2 we show the sample complexity of training when the underlying nou-dimenssioall distribution has low complexity, medium complexity, or high complexity (top to bottom). For each of the three data distribution, we then train models of three different complexities and evaluate their train-time and test-time KL error. In the bottom row, which peaks at subpar performance. In the top row, leads to overfitting which learns more slowly and less efficiently than the low-complexity model (even with multiple thousands of samples). In addition to showing the importance of matching the correct structure to achieve optimal performance, these experiments also show how achieving good generation performance automatically generalizes to classification performance.
how the underspecified, low-complexity model leads underfituing how the overspecified, high-complexity model
In Figure 3 we plot all different values of refined information for our high complexity distribution. This gives preliminary insights into the properties of refined information and the range of values single
8

---

Published in Transactions on Machine Learning Research (04/2026)
KL Error v. Sample Size (log)
Clax Sample Size (lgg)
low drain)
meduum traam)
high terain)
law 0:5t
(see))
High (tea)]
Figure Model Performance vs. Number Training Samples. Left three panels: KL error as a function of sample size. Right twelve panels: Normalized classification error a function of sample size. Performance is evaluated across three different model complexities. Each corresponds the complexity of the underlying data distribution which from top to bottom has low complexity, medium complexity, and high complexity. The top shows the high complexity model overfitting and slowly fitting. The bottom row shows the low complexity model underfitting. Error bars are with respect to 5 different resampling of the synthetic training dataset.
Refined Information (log) v. Mode Interaction
Retined Information
2D Refined Information
3D Refined Infotmation
Refined Information
-1.5
Mutual Information
-2.5
-3.5
of
-4.5
-5.5
Figure Synthetic four-dimensional data with all possible values of refined information plotted. Each column represents an interaction [d]. Each point represents a particular value of RIz,s for some possible pair (I.S). Random horizontal spread is used for visual clarity.
interaction S can take depending on the context Z. We additionally plot the 2D marginal refined information which corresponds to the classical definition of mutual information. Further discussion of its properties and applications to structure learning and generalization are saved for Appendices B and
5.2 RorldWorld Results
For the real-world datasets, first demonstrate the tuning of our MIS hyperparameters (heredity strength, heuristic norm, and loss function) on small subset of the real-world dataset. By using only the first ten dimensions of the mushroom dataset, we work in regime where the exact gradients can be readily calculated and we perform our algorithm with collections of up to size 300.
Figure 4 shows the performance terms of the capacity curves, plotting the training and validation errors as function of the size of the interaction collection, which is a measure of the log-linear model's capacity. We train significantly beyond the point of early stopping to help fully illustrate the underfitiing-ovrfftingg

---

Published in Transactions on Machine Learning Research (04/2026)
Table 1: Statistics for real-worke datasets.
71
d
|44
mushroom
8,124
23
2.4e14
adults
32,561
14
6.5e11
breast cancer
286
10
6.0e05
behavior of the log-linear model. This provides empirical support for our theoretically principled approach of early stopping as soon as the validation error stops improving alongside the training error. Further hyperparameters are provided in the appendix. Overall, we find that using the weak hierarchy of 30% strength was the most effective choice for achicving minimal validation crror for our MAHGen Ta algorithm and we keep this choice consistent throughout.
We then apply these hyperparameters our two large-scale datasets where we cannot directly calculate the KL divergence and resort to the AIS approximation discussed in Section 4.2 For our third real-world dataset, exact KL gradients are still calculable with an event space smaller than one million. We compare against Boltzmann machine and an independent distribution also trained with gradient descent on the same objective. In Table [ we compare our approach which has the capacity to learn sparse and higher-order structures against both the 1-body and 2-body log- linear models. We find that our MAHGenTa approach is able to consistently deliver improvements generation performance terms of the KT. divergence or log-likes li-liieeihoodd.
30% Weak Heredity with Normalization
Strong Heredity with Normalization
2-body
2-body
3-bodly
3-body
KL Error (log)
4-body
KL Error (log)
4-body
body
body
Callection (log)
Collection (log)
Figure 4: Comparison of the training (dark) and validation (light) error as we continue to add interaction subsets to our collection. On the log-log plot, we see clear distinction of the underfitting phase and the overfitting phase. We also see how assuming strong heredity can lead to 'discontinuttes" in the performance, caused by important higher-order interactions being blocked in the greedy algorithm.
Table 2: KL divergence across all real-world datasets.
mushroom
adults
breast
independent (1D)
15.477 1 0.056
8.692 1 0.047
5.991 1 0.210
boltzmanu (2D)
4.472 1 0.069
6.444
0.042
5.652 1 0.105
mahgenta (3D+)
2.212 + 0.062
5.832
0.012
5.176 = 0.052
5.3 Discussion
In Table 3 we see how the training of generative model automatically leads to emergent capabilties in classification via the mode interactions simplifying into feature interactions. In particular, the encrgy-based MAHGenTa and Boltzmann machine are able to simultaneously predict across multiple classes, unlike the discriminative approaches which must be retrained for each task. Although the discriminative approaches have the advantage of reusing the dataset learn only one of the conditional distributions the generative approaches nevertheless yield a comparable accuracy performance across a variety of tasks simultaneously.
time,
10

---

Published in Transactions on Machine Learning Research (04/2026)
Table 3: Class-wise accuracy across all real-world datasets.
mushroom odor
adults race
breast cancer
poison
habitat
income
gender
recurrence
malignance
mahgenta
99.7
79.7
66.0
85.2
86.5
84.5
80.2
51.6
boltzmann
98.2
78.5
63.4
84.2
84.9
83.6
72.1
50.4
logistic
100.0
85.6
80.6
88.0
71.3
regression
65.8
84.4
42.7
naive
94.8
78.6
81.6
85.3
72.0
bayes
63.3
82.1
44.1
In the adults dataset, we can clearly see how single generative model trained to adcquately model the data easily obtains good accuracy not only for the original target of income level, but also sensitive features like race and gender. the classification setting, may be unclear that a model biased using sensitive features to predict income; however, our energy-based model working directly on the observed variables, the learned connections between variables are made explicit. This could have implications for algorithm fairness approaches, where removing sensitive feature labels from the training data is not sufficient remove the fundamental bias which exists within the dataset. In contrast, biased energy terms in the log-lincar model could be directly inspected, analyzed, and removed.
6 Conclusion
Overall, we find that refined information opens up many directions for further exploration of higher-order information and that mode-interactionn-elctiio for hierarchical log-linear modeling is an effective tool in reducing the number of parameters to be learned in a principled way. Theoretical developments allow for a complete decomposition of the KL error in terms of the refined information content. The regularizing effect of choosing simpler structure is made clear on both synthetic and wealdworld datasets, with an easy-to-use carly stopping hcuristic achieve optimal performance. The benefits of gencrative distribution lcarning general pretraining objective for multiple downstream tasks also reinforced.
11

---

Published in Transactions on Machine Learning Research (04/2026)
References
David H. Ackley, Geoffrey E. Hinton, and Terrence J. Sejnowski. A learning algorithm for boltz-
mann machines. Cogmitive Science, 9(1):147 169, 1985. ISSN 0364-0213. doi: https:/ /di:
URL
50364021385800124
S. Amari and H. Nagaoka. Methods of Information Geometry. Translations of mathematical monographs. American Mathematical Society, 2000. ISBN 9780821843024.
URL httpss//books: http://bbok...oogl...... co. jp/
boots/1d=vc22W507/.00C
S.-I. Amari. Information geometry on hicrarchy of probability distributions. IEEE Transactions on Infor- mation Theory, 47(5))1701 1711, 2001. doi: 10. 11109/18.930911
Shun-Ichi Amari. Information Geometry and Its Applications. Springer Publishing Company, Incorporated, 1st edition, 2016. ISBN 4431559779.
Anil Aswani. Low-rank approximation and completion of positive tensors. SIAM Journal Matric Analysis and Applications, (7(33)133741364, 2016. doi: 10.11371161107818.. URL https://doi.org/10..11
16M1078318
Mohamed Bennasar, Yulia Hicks, and Rossitza Setchi. Feature selection using joint mutual information
maximisation. Expert Systems with Applications, 42(22) 8532, 2015. ISSN 0957-4174. doi: https:
10..1101.1.......055...000000000 URL https://wwwt wwtp///caaaeeet....////////////a
50957417415004674
Jacob Bien, Jonathan Taylor, and Robert Tibshirani. A lasso for hierarchical interactions. slalistics, 413):1111,, 2013.
Annals of
Søren L. Buhl. On the existence of maximum likelihood estimators for graphical gaussian models. Scandi- navian Journal of Statistics, 20(3):263 270, 1993. ISSN 03036898, 14679469. URL http::/www.jsttee
org/stable/4616281
Zhijun Chen, Chaozhong Wu, Yishi Zhang, Zhen Huang, Bin Ran, Ming Zhong, and Nengchao Lyu. Feature selection with rohuddncg--ommlemntrriesss dispersion. Knowledge- Based Systems, 09:203-217, 2015. ISSN 0950-7051. doi:
http:::/wwwwcccenndddreer
A. P. Dempster. Covariance selection. Biometrics, 28()):1177 175, 1972. ISSN 0006341X, 15410420. URL http://www.p5
orggsttable/228966
James Enouen and Yan Liu. Sparse interaction additive networks via feature interaction detec- tion and sparse selection. In S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and
A. Oh (eds.), Advances in Neural Information Processing Systems, volume 35, pp. 13908 13920. Cur- ran Associales, Inc., 2022. URL https:
cttrippscc/paperrr http://///////ppppaaaaa////////p ffiles//ppperr/222//ieee 63307888006613330009383332283333 Paper Conf-rCooeeeeeeeee pdf
Yingying Fan, Yinfei Kong, Daoji Li, and Jinchi Lv. Interaction pursuit with feature screening and selection, 2016.
R. A. Fisher. Two new properties of mathematical likelihood. Proceedings of the Royal Sociely of London. Series A, Containing Papers of Mathematical and Physical Character, 144(852):285 307, 1934. ISSN 09501207. URL http:://ww, jstorr org/stablee/935559
Elad Ganmor, Ronen Segev, and Elad Schneidman. Sparse low-order interaction network underlies highly correlated and learnable neural population code. Proc. Notl. Acad. Sci. U. S. A., 108((2):967999666)) June 2011.
Alan E. Gelfand. Gibbs sampling. Journal of the American Statistical Association, 95(452): 1304, 2000. ISSN 01621459. URL http://www jiw.rrstoooooooo org//staalee26667755
12

---

Published in Transactions on Machine Learning Research (04/2026)
Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEFF Transactions on Pollern Analysis and Machine Intelligence, PAMMI-(6):721--77 -411, 1984.
doi:
Kazu Ghalamkari, Mahito Sugiyama, and Yoshinobu Kawahara. Many-body approximation for non- negative lensors. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (ods.): Advances in Neural Information Processing Systems, volume 36, pp. 74077-74102. Curran Associates, Inc., 2023. URL https://prooeeddnng..m nww.gg.....ri...... httpar/issccppppee..a//////////p pteee///e/////////////1e eB2357760D01aat11155556668648883
Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Parley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. In Z. Ghahramani, M. Welling, C. Cortes, N. Lawrence, and K.Q. Weinberger (eds.), Advances in Newral Information Processing Systems, vol- ume 27. Curran Associates, Inc., 2014. URL https://prooceedings www.igs...ui.s... htcppss.cc/pppppee ftaperrffile//pperrr 2004:///hnnccacce2222266666666e2
Te Sun Han. Multiple mutual informations and multiple interactions in frequency data. Information and Control, 46(1):26-45, 1980. ISSN 0019-9958. doi:
URL
https://www. wtt.://wwcec.ettccccoococcaaaaaa
G. E. Hinton and R. R. Salakhuldinov. Reducing the dimensionality of data with neural networks. Science 313(5786):504 507, 2006. doi: (11126/science.1276644 URL bttps::/www.sciemeee wwtpaa/ww..///////////aa/a/a////
1126/science 1127647
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusson probabilistic models. In H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (eds.), Advances in Neural Information Processing
Systems, volume 33, pp. 6840 6851. Curran Associates, Inc., 2020. http:///tee/pprreeenn....o..
Søren Højsgaard. Statistical inference in context specific interaction models for contingency tables. Scandi- navian Journal of Statistics, 31(1):143-158, 2004. ISSN 03036898, 14679469. URL http://www.jstaee
Diederik Kingma and Max Welling. Auto-Encoding Variational Bayes. In 2nd International Confer- ence on Learning Representations, ICLR 2014, BanjJ, AB, Canada, April 14-16, 2014, Conference Track Proceedings, 2014.
Su-in Lee, Varun Ganapathi, and Daphne Koller. Efficient structure learning of markov networks using 1_1 regularization. In B. Schölkopf, J. Platt, and T. Hoffman (eds.), Advances in Neural Information Processing Systems, volume 19. MIT Press, 2006. URL httppp/pp//poooooio.g...oo.opppp
pdf
Daniel Lowd and Jesse Davis. Learning markov network structure with decision trees. In 2010 TEEE International Conference on Data Mining, pp. 334 343, 2010. doi: 10.1109/ICDM.22100.222
Fuyuan Lyu, Xing Tang, Dugang Liu, Chen Ma, Weihong Luo, Liang Chen, xiuqiang Hc, and Xue (Steve) Liu. Towards hybrid-grained feature interaction selection for deep sparse nelwork. In A. Oh, T. Naumann, A. Globerson, K. Saenko, M. Hardt, and S. Levine (eds.), Advances in Newral Information Processing Sys- tems, volume 36, pp. 49325 49340. Curran Associates, Inc., 2023. URL URL/http:::/pproceednnng neurips. .dd
Hélène Massam, Jinnan Liu, and Adrian Dobra. A conjugate prior for discrete hierarchical log-linear models. The Annals of Statistícs, 37(6A):3431 3467, 2009. doi: 10.1214/08-AAO AOS669, URL https: //do:. org/ 10. 1094//00--A
W. McGill. Multivariate information transmission. Transactions of the IRE Professional Group on Infor- mation Theory, 4(4):93-11) 1954. doi: 10.1109/TTT. 1T51//TTTT.994..7..655
13

---

Published in Transactions on Machine Learning Research (04/2026)
Martin Renqiang Min, Xia Ning, Chao Cheng, and Mark Gerstein. Interpretable Sparse High- Boltz- mann Machines. In Proceedings of the Seventeenth International Conference on Artificial Intelligence and Statistics, volume 33 of Proceedings of Machine Learning Research, pp. 614 622, Reykjavik, Iceland, 22 25 Apr 2014. PMLR. URL https:://proceedings .ml. www.s/..e./e............e...eeee aattass/33//mini/ html
Hiroshi Nagaoka and Shun-ichi Amari. Differential geometry of smooth families of probability distributions. Technical report, University of Tokyo, 1982.
Hiroyuki Nakahara and Shun-Ichi Amari. Information--oometric mcasure for neural spikes. Neural Comput.. 14(10):2269 2316, October 2002.
Songyot. Nakariyakul. High-dimensional hybrid feature selection using interaction informationngguideed search. Knowledge- Based Systems, 145:59 66, 2018. ISSN 0950-7051. doi: http:://ooo.orr///000006..... 100o 01.002. URL
Radford M. Neal. Anncaled importance sampling. Statistics and Computing, 11:125-139,
Frank Nielsen, Frank Critchley, and Christopher T. Dodson. Computational Information Geometry. Signals and Communication Technology. Springer Publishing Company, Incorporated, 2017. ISBN 9783319470566. URL
wwt.cooo/bbooecoo 1007 cotbo//tc.oeooooooooo/////////// /ttp:///////0011000000000------0 1007/978---319--70000
Henrik Nyman, Johan Pensar, Timo Koski, and Jukka Corander. Stratified Graphical Models Context- Specific Independence in Graphical Models. Bayesian Analysis, 9(4):883 908, 2014. doi: 10.1214/
14-BA882. URL https://doi. www.ioorr.1... 10. 1214/14-B8882
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. Pytorch: an imperative style, high-performance deep learning library. In Proceedings of the 33rd International Conference on Newral Information Processing Systems, Red Hook, NY, USA, 2019. Curran Associales Inc.
Julio L. Peixoto. Hierarchical variable selection in polynomial regression models. The American Statistician, 41(4):311 313, 1987. ISSN 00031305. URL hnttp://ww. http://www.jttrr.oo/../oo jstor org/stabble/664755
C. Radhakrishna Rao. Information and the Accuracy Attainable in the Estimation of Statistical Pa- rameters, pp. 235-247. Springer New York, New York, NY, 1992. ISBN 978 978-1-4612- doi: 16. URL https://dii.orrr
10.1007/978--446122000
Ruslan Salakhutdinov and Geoffrey Hinton. Deep boltzmann machines. In David van Dyk and Max Welling (eds.), Proceedings of the Twelth International Conference on Artificial Intelligence and Statistics, vol- ume of Proceedings of Machine Learning Research, pp. 448-455, Hilton Clearwater Beach Resort, Clearwater Beach, Florida USA, 16-18 Apr 2009. PMTR. URL https://proceedings mlr. f.l...pres/vv// salakhutdinorOSa.fP
G. Schay. Constrained differentiation. Mathematical and Computer Modelling, 21 88, 1995. ISSN 0895-7177. doi:
URL https://wwwt /www www.sceeiiiicet....
Mark Schmidt and Kevin Murphy. Convex structure learning in log-lincar models: Beyond pairwise potentials. In Yee Whye Teh and Mike Titterington (eds.), Proceedings of the Thirleenth Interna- lional Conference on Artificial Intelligence and Stalistics, volume 9 of Proceedings of Machine Learn- ing Research, pp. 709 716, Chia Laguna Resort, Sardinia, Italy, 13-15 May 2010. PMLR. URL httts::/prrocedding mlr. wwr..press/vee http://////ccmmmi1........t
Terrence J. Sejnowski. Higher-order Boltzmann machines. In Neural Networks for Computing. volume 151 of American Instâtute of Physics Conference Series, pp. 398 403.4 ATP, August 1986. doi: 10.1063/1 00.1063/1.3624 11.6446.
14

---

Published in Transactions on Machine Learning Research (04/2026)
C. E. Shannon. A mathematical theory of communication. 379-423, 1948. doi: 10.1002/j-.. 110002//j-13--1111111 1538- 194...994...013333.
The Bell System Technicol Journal, 27(3):
Iya Shpitser, Robin J. Evans, Thomas S. Richardson, and James M. Robins. Sparse nested markov models with log-linear parameters. In Proceedings of the Twenty-Nindl Conference on Uncertainty in Artificial Intelligence, 1Ar13, pp. 576-585, Arlington, Virginia, USA, 2013. AUAI Press.
Mahito Sugiyama and Karsten Borgwardt. Finding statistically significant interactions between continuous features. In Proceedings of the Twenty-Eightt International Joint Conference on Artificial Intelligence, IJCAT-19, pp. 3490 3498. International Joint Conferences Artificial Intelligence Organization, 7 2019. doi: 00.24963/ijcai2200 2019 484. URL https://dol www.orog//........oo.ooo 24963/ 2019//84
Mahito Sugiyama, Hiroyuki Nakahara, and Koji Tsuda. Legendre decomposition for tensors. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (eds.), Advances in Neural Informalion Processing Systems, volume 31. Curran Associates, Inc., 2018 URL https://proceddddng.. neurips. http:///w.....atoooo///////////a
Robert Tibshirani. Regression shrinkage and selection via the lasso. Journal of the Royal Statistical Society. Series B (Methodological), 58(1))267-283) 1996. ISSN 00359246. URL http::/www.jstor.orr/aaalll
2346178
Jan Van Haaren and Jesse Davis. Markov network structure learning: A randomized feature generation approach. Proceedings of the AAAI Conference on Artificial Intelligence, 26(1):1148-1154, Sep. 2012. doi:
10.1609/aaai.v2 .8315. URL. https://ojjs aaai www./odex../oeeeeooooooooooee..t http://wwa..pppphhh//aaaaattatta
Martin Wainwright, John Lafferty, and Pradeep Ravikumar. High- nigh-dimensional graphical model selection using f1-regularized logistic regression. In B. Schölkopf, J. Platt, and T. Hoffman (eds.), Advances in Newral Information Processing Systems, volume 19. MIT Press, 2006. URL httpss://pooceeddingg www.:gg...uuriirr.om cc/
paper_
Zilin Zeng, Hongjun Zhang, Rui Zhang, and Chengxiang Yin. A novel feature selection method consid- ering feature interaction. Pattern Recognition, 48(8):2656-2666) 2015. ISSN 0031-3203. doi: https:// doi.org/
500311220115000350
15

---

Published in Transactions on Machine Learning Research (04/2026)
A Information Geometry
A.1 Necessary Background
Information geometry is the field which treats the set of parameters E IRP as a manifold with a Riemannian metric, called the statistical manifold. In our case, we will focus on the structure introduced by KL diver- gence. Accordingly, we will write the divergence as forward KL divergence D(p:)) g) (KL q). the dual divergence reverse KT. divergence D" (1:
DRL p), the Bregman potential as negative entropy s(po) E, po( log po)), and the dual Bregman potential as free energy (p(p)) log E exp ())g()))) For now, write Ej the one hot basis for the discrete space ei(i) == 1(i j), further discussion on alternate coordinate systems will be had in the next subsection.
We write the distribution corresponding to some 0 parametrization as:
(i( exp { -0e
(10)
where it can be seen that (pa) already chosen such that
Σvl(i) Σ{ exp E [oses(i) exp { (log Σ exp Σ [eee(())3
Σ exp Σ (eee((i)
(Σ exp Σ (rei(H)) -1 1.0
(11)
Using the one hot basis Gj, this means that we may write po(i) = (xpp1))))),.
exp(0,1).
However, to ensure there are no redundancies in the manifold, we are required to choose to drop one coordinate. Otherwise, we could rescale by shifting by a constant. so (1) = 1/(,, exp(0,)). Again, we discuss our alternative later in Section A.3
convention drop the first coordinate,
This choice of D, sp, and &b allows for a duality between the 0 coordinates and the r coordinates corre- sponding to the expectations of the basis functions πソ Eg [ei], which has the following nice consequences:
Theorem 1. (Legendre Transform Formula, Theorem 1.1 of (Amari & Nagaoka 2000
DKL(p:g) = p(p) ())) m(p))
(12)
Proof. DKL is a Bregman divergence of the form Dip, meaning that DKL (KI(P: s(p) - -(g)- (p g).
Because we know from the Bregman dually flat structure that 0 = 11" V)((n) and n = 0 = (v) and we also have the Legendre duality by our choice of and vi, giving v(g) + vb(g") = n(()) 0((), we have that:
(p)(p: = =(p)) =p()) - O(q) (n(p) - (()))
(13)
= ((p) (w(g) v(g)) -(()) 0(g) (n(p) n(4))
(14)
(p)) +x)) 0(g) ((p)
(15)
Theorem 2. (Pythagorean Theorem, Theorem 1.2 of (Amari & Nagaoka 2000) Given a dually flat manifold, with distribution p,q. such that p and g are connected by an ?) geodesic and g and r are connected by a 0 geodesic, which are orthogonal, then the generalized Pythagorcan theorem holds:
(p:(p:)) (RI(p: g) + DRI(grr).
(16)
16

---

Published in Transactions on Machine Learning Research (04/2026)
Proof. We may write that:
4) + D( = (v(p) + +(a) 7(p) -0()) + (p(() + (((r -(()))))()))))))))))))))))))))>) 0(r)]
(17)
(s(()) + ((r) - m(p) -
(p(2) +(x)) -(p) -(()) -7()) 0(r))
(18)
D(p;r) - (n(p) - --() (0(g) - 8(r)
(19)
= D(prr)
(20)
where the last step follows by the orthogonality of the geodesics. For each coordinate i, we have that either (0(p) (g)]a is zero [m(4)= -f]]] zero, meaning their dot product is zero overall.
Theorem 3. (Projection Theorem, Theorem 1.4 of (Amari & Nagaoka, 2000) Given a dually flat
manifold and considering a submanifold M, for example one defined by a subset of the coordinates, the point that minimizes the divergence is the dual geodesic projection and the point that minimizes the dual divergence is the geodesic projection.
Proof. Proof is completed by considering the geodesic project and dual geodesic projections a small neighborhood where the Pythaogrean theorem then shows that any small deviation can only create positive increase in divergence or dual divergence (due to the non-negativity of divergences). This confirms that the dual projection and projection are indeed critical point.
A.2 Hierarchical Coordinates
In this work, we follow the hierarchical coordinate scheme of Ghalamnkari el al. 2023 in order to adequately respect the many-body structure the mode interactions between the different variables.
qo(i) exp { 1 0
It should now be understood how naively using these coordinates results in an overspecified number of dimensions compared the actual manifold (which is (I -imdeeesii,nal) Accordingly, we first consider the approach of (Ghalamkari et al. 2023) which zeroes out the first coordinates at cach hicrarchical level.
See the later examples in Section A.3
These hierarchical coordinates are already sufficient for defining refined information, so we first confirm this hierarchical representation remains valid under information geometry.
First writing the free energy constant 0 (b) as:
g0_ -oo Σ exp o (Eocecia
{E
0 is
11
As mentioned, we can envision the parameter 0 either as function of the other 0 or as a parameter which is then constrained by the other parameters. Although the function perpsective needed to nicely align with the theory of information geometry, we also find it useful take the latter perspective during the practical fitting of these energy- based model. For simplicity, we will for now take the former perspective and look at the derivative of 0 or the free energy taken function of the other parameters (but this
17

---

Published in Transactions on Machine Learning Research (04/2026)
also be thought of as the derivative of the constraining equation).
{
(00
log
Σ exp
Σ
0/5
{}}}
00
007
yr
Σ exp Σ DESEEd
0 i8
$}}}
0
Σ exp
E
0%0 {%
00T
exp {
80)
0
Σ acscia
05+ }}}}
O07
exp
-eep
exp
) {
exp Σ @CSCLd
0 is)
Σ Sr exp {Σ SC 0
}}
Zómnt 20(1) = Edrn.r beritr
--1 Ser fr my - Th
Eder
Str.sr Thy 755
Now it will be easy for us to take the derivative of the objective function, the forward KL-divergence:
Ptrn(i) log Ptrn (i)
DKI(Ptrni go) = -Perrnl
Let us write:
DRL. (perni qo)
Σ
ptrn (i) (og(ptrn(1)) Ptra (1) log(ga(i))
T07 jr
0
Σ (og((((()))) = Σ
θ
pm (i) (1)
O0T
dhT 2T
JT
trn
00T JT
log(g (i)) = Σ
p
trn
20T
Σ sclu
JT
Σ
p
tr
{ 8 + p% }
trm
a
000) Σ p(r(())
{0f Jr
00T
000
20T
jr
fr
{
000
p tra S.r
007
0.1 M5r
0.1 trn.T prnn (ir) = =ījr
This completes the calculation of the gradient for the paramctrization of the probability distribution. Wo reiterate that this is the derivative with the free energy constraint, but without any additional constraints to enable identifiability of the parameters. In particular, this means that multiple 0 parameters can correspond to the same probability distribution. This is alleviated by the zeroing out to restrict the manifold done in (Ghalamkari et al., 2023). In the next section, we introduce our final version of the information-geomen coordinates, which instead use coordinates which are both hierarchical and centered, which we find to be critical for the practical deployment of these methods.
A.3 Centered and Hierarchical Coordinates
In this section, we first revisit the zeroing out constraints of Sections A.1 and A.2 using explicit examples before introducing the centered coordinates which we introduce as an alternative with computationally nicer properties needed for practical implementation.
Textbook Constraints The textbook method of providing additional constraints is to zero out all of the unnecessary parameters for the probability distribution. This is typically done by zeroing out the 'corners'
18

---

Published in Transactions on Machine Learning Research (04/2026)
of each 8 parameter via:
00 -0 if i = for any 8 8 S.
We omit the exact details of showing this will represent each probability distribution with a unique set of parameters; however, we note that this leaves each @ with I -1,1--) parameters out of the original [ees(( Tsl- Accordingly, all of the 0 parameters together have |T4a degrees of freedom, and after the sum-lo-one constraint. from the last section on 0% there are the required |/4/ - degrees of freedom.
In the textbook formulation, these zeroed out coefficients are usually completely dropped from consideration. After additionally considering the 8 a function of the other parameters rather than another parameter itself, we are left with parametrization of the manifold which has the same number of parameters are there are intrinsic dimensions in the manifold.
We could come to an equivalent such parametrization with any set of choices for the 'base events' i E [4,] instead of simply choosing the first event. However, some choices can be arbitrarily worse than others and searching over all possible choices of base events is much more work than it is worth. We will briefly introduce the concerns of numerical stability in some small, low-dimensional distributions. On one-dimensional distribution with four outcomes, this would become:
p(1) e
p(2) 6 +01
p(3) e0 +03
p(4) e
In this case, can write the closed form solution
0 = log( ((2)
0g log( p(3) p(1)
0 = log( p(4) p(1)
0e log(p(1)),,
p(1)
Fortunately, this seems to mcan that as long as p(1) is not extremely large or extremoly small compared to the other probabilities, there will not be a huge amount of issues with our arbitrary choice of base event for the exponential distribution. However, the risks associated with this arbitrary choice will only continue to grow and compound as we increase the dimensionality and the number of events.
In two dimensions, we can see:
p(1,1) = 0² p(2, 1) = +0
2(1,2) = e -00+
p(1,3) = 00000
((2 = +00
1(2,3) = e00 +00 +02+
Again, we may write closed form solution as:
03 = log pp(1,2 2)
03 log (P(1,3)
0 = log(p(1, 1));
pp(1,1)
p(1,
8- log (P(1,2) (p(1,
012 log (p(2,2) p(1, 1) (((,)) p(2,
(2(2,3) p(1, ((1,3) p(2, 1)
012 log
It can be imagined how increasingly high dimensional distributions would only exacerbate the issues of potentially imbalanced events within the space, especially the high dependence on the probability p(1...., .,, 1). It is important also remember that for our main application, there does not exist a closed form solution for the 0 parameters and we must use gradient-based learning approaches, necessitating adequate handling of these potential numerical issucs.
Balanced Parametrization Given these concerns of significant and unnecessary challenges in the gradient-based based learning process for the 0 parameters, we instead leverage an alternative identifiability con- straint which is more compatible with initialization at the all zeroes vector. In particular, we balance the 0 tensor around zero by assuming that each fiber sums to zero:
Σ
S i/ 13- =0
for any & E S for any i(s-0) E ((S-s)-
(21)
J.
19

---

Published in Transactions on Machine Learning Research (04/2026)
It is again relatively straightforward to verify that this gives unique parametrization for every probability distribution in the manifold. We can also see that this condition can be written in a more tensorial form as:
E =0 d.
for any E
where 0 refers to the zero tensor 0 R Let also recall the r parameters defined as ng := Bj-g(t) l[ig = js)] = (is) = En g(ss,-ss). We may immediately see that the parameters auto-
matically obey similar tensorial constraint as:
Eng, n s-=
for any S.
Both the 0 and the y sets parameters form a tower structure made of tensors. We may immediately try to purify the yj tower a similar way to our O's via the principle of inclusion-exclusion (related to the mobius inversion of the zeta function), namely let write:
S
Σ (-1)551-17)) "7" 8u S-T
(22)
TCS
It is fairly straightforward to see that these new parameters obey the centralized condition:
Σ = ō
for any s SES
a.
and that moreover the variables are recoverable in the opposite direction via:
17 =
Σ
π 2
S-T.
TCS
More importantly, may write the purified gradient of Equation in terms of which further corresponds to the constrained gradient of the KL divergence:
V DRL (Peni [o)] - m
Although this already sufficient for the purposes of our work, but us proceed slightly to round out the discussions about our parametrization as it relates to the typical topics of information geometry. First, we recall that the typical Bregman divergences energy and total entropy. We have already computed the 0 derivatives for free energy when looking al the derivative of forward KL, but let us also take look at the backwards case with entropy.
forward and backwards KL with Bregman functions of free
First recall that
H(p) = -E(O log (p(i)
and also that
Ox flo(ff) = =If -1f -lg(loo(f)) ff f (1gg
20

---

Published in Transactions on Machine Learning Research (04/2026)
0 H (qo)
a J
Σ
Σ
Tig
S
uusss) log
Tig
2 u
8
1 Σ
log Σ
8
Σ
8
& u s
iis
@
dt!
Ris
8
Σ |+ log (90(i)) |
ΣΣ
mig
& 4
5
Jr
055
" T
9[dd
Sir
[1+
0%5
scId
TT
1
0
r
Jrnsis
sCid
|T
1
+
SCT jirns
Eser
Σ
0 sgr irnsis r
Σ
+ Σ SgT
1
g Jrns.s T
SCT
Jms
[T
Σ SCT
0
+
Contimuing on,
H ((g))
Σ
S Tig
8 ur-s) log (prr((i)
mm T
= [ Σ SCT
0
T -- (Brrir u
Σ Escial
grra.s)
00
scld
Σ
SCT 05 - [1+ Σ SCT
Using again the purification of the tower of I's leaves us only with the top of the sum of f's, meaning that the purified derivative of the free energy is only left with the corresponding Of, parameter. Finally, we can write:
V (2:L. = =0S_ - 0
21

---

Published in Transactions on Machine Learning Research (04/2026)
B Experimental Details
B.1 Additional Results
In Figures
we show the capacity curves across all sets of MIS hyperparameters chosen for the experiments with the 10-dimensional subset of the 23-dimensional mushroom dalasets.
5 and B
True hood 35% Weak:
Wf Narm
Thue hood: 50% Weak:
of Nor
Tre Libslbozd: Sumii Shang: ai Verm
True Listhood Stron:
term
Secill Colerio the
She Coledionnnn
Calerie (he)
Figure 5: All hyperparameters of heredity strength and parameter count renormalization. Top 8: Full likelihood training.
5000
2000
Penusto Liet hhedt 32% Weal:
Narm
Punudo (tei hood: 52% Wait::
wnto Norm
Proust: the head: Sami Strngg of Verm
Pesudo- Hatihheddddd Streng:
W/ Nerm
Soa Necion teg
Ste Coledin
Ste Collletee Jeg)
Sin Colect 19g)
Figure 6: All hyperparameters of heredity strength and parameter count renormalization. Bottom 8: Pseu- dolikelihood training (masking turned off to avoid -cc).
Table 4: Best validation KL error from all MIS hyperparameters and all many-body solutions.
Heredity
Likelihood Type
Param Count Renorm
30% Weak
50% Weak
Semi
100% Strong
True Likelihood
with
0.2359
0.2466
0.2445
0.2559
True Likelihood
without
0.2372
0.2489
0.2472
0.2555
Pseudolikelihood
with
0.3659
0.3533
0.3659
0.3817
Pseudolikelihood
without
0.3589
0.3-198
0.3519
0.3583
Many- (no sparsity)
Likelihood Type
1D
2D
3D
True Likelihood
4.6062
0.8281
0.2644
I'seudolikelihodd
4.6062
1.5005
0.6579
The search over these hyperparameters can be seen as a comparison to previous works (Schmidt & Murphy 2010 Min al. 2014) bocause of their use of other types of stronger hicrarchical assumptions. Morcover, the stagewise-selection procedures can also be seen as a special case of this mode interaction selection framework. Accordingly, the only missing component of full comparison to these previous works would be tuning over
22

---

Published in Transactions on Machine Learning Research (04/2026)
the L1 regularization parameter. We find unlike those works, tuning an L1 parameter is not necessary in our work due to our L.0 selection with the MIS algorithm and our theoretically supported early-soppping procedure. It is moreover emphasized that both previous works have no available code and regardless were only designed for binary variables, making them inapplicable to any of our datasets used herein.
B.2 Synthetic Datasets
We generate small synthetic distributions by first drawing 0 from a Gaussian distribution with unit covari- ance. We then center each 0 such that the sum across any mode is equal to zero. A Afterwards, there are any S we have not yet zeroed out from the model, then we do this now. Finally, we compute the probability distribution as the exponential of the sum of the B paramcters and compute the renormalization constant if necessary. Synthetic training dalasets are then drawn iid from this final distribution.
In experiments, we use d = 4 dimensional distributions inside 4d [5] medium, and high complexity distributions, we use the following sparsity patterns:
[3]
5
For the low,
low complexity
To = 4, 13, 14, 23, 123] 12,
medium complexity high complexity
Imed = 1,,
Thigh {0, 1,, 12, 14, 23, 24, 34, 123, 124, 134, 234, 1234)
We use sample sizes of (10,20,40,88) 16,, 320, 640, 1280,, 2260, ,,2,,,512,,,,,,, ,,,,,,1,,,,,,, 10240] and plot from 10 to 10,000 on a logarithmic scale.
KL Error Sample Size (log)
Class Score v. Sample Size (log)
(tran)
low medium (sae)) high (laak)
necke
Figure Model Performance vs. Number of Training Samples. Performance evaluated across three different model complexities. Each row correspond to an underlying data distribution which from top bottom has low com- plexity, medium complexity, and high complexity. The top demonstrates the high complexity model overfitting. The bottom row demonstrates the low complexity model underfitting Left hand side is the KL error objective de- creasing; right hand side is the class-wise performance (which is automatically gained from generative performance). Error bars are with respect to 5 different resampling of the synthetic training dataset.
In Figure ? we show the sample complexity of training when the underlying four dine-imennional distribution has low complexity, medium complexity, or high complexity (top to bottom). On the left-hand side, we see the KL error optimized during training, whereas on the right-hand side, we see the calibrated classification score for each of the four dimensions (predicted using the other three features), automatically rising alongside improved generative performance. In the bottom row, we see how the underspecified, low-complexity model leads to underfitting which poaks at subpar performance. In the top row, we see how the overspecified, high-complexitty model leads overfitting which makes less efficient use of the finite dataset (even with multiple thousands of samples). In addition to showing the importance of matching the correct structure
23

---

Published in Transactions on Machine Learning Research (04/2026)
to achieve optimal performance, these experiments also show how achieving good generation performance automatically generalizes to classification performance.
KL Error v. Sample Size (log)
Score (lg)
High (baks)
hebut
"elux
matgeenia (tren)
low Itest
medium Noad1
Nigh (eas))
mratigenia (swe))
Figure & Model Performance vs. Number of Training Samples. Also with Mahgenta automatic selection performance (black) compared against oracle structure performance. As discussed in our theoretical section about overfitring, for low sample sizes, Mahgenta actually outperforms the oracle structure. In most cases, mahgenta very nearly matches oracle performance for larger sample sizes.
B.3 world Datasets
We next apply our method to the real-world distribution of three different UCI machine learning datasets. Testing split is generated from 50% of the data and kept fixed throughout. The remaining data is split 70%/30% into the training and the validation set. Experiments are run on a Tesla V100 32GB GPU, The three datasets used are shown in Table I with their nunbers of samples, features, and total possible events.
B.4 Algorithm Details
In Algorithm
2
we give full description of the subroutines used by main algorithm.
B.4.1 Heredity Styles
Recall that we define the heredity score as:
== 5 = {i} for some i (TT: T' has already been selected]]
(23)
Thus, writing the count for a particular S as ni = w(S) and its final score as = what percentage of subsets T C S, of size one less than S, are alrcady selocted by the algorithm. In the case of pairs, there are the two options of both singles already selected (strong heredity) and only one single already selected (weak heredity). Generalizing to larger S provides much greater fiexibility. We focus on class of heredities defined by single +% parameter, which asks that the percentage of 'oncless' subsets are at least T, i.e. SSS)/S|| > т. In the experiments in Section B.I we consider this for values of T = 30%, 50%, 100%.
we can ask
We also consider a heredity style which we call 'semi-strong' which instead asks that w(S) 2 |S| - 1, meaning that all but one of the 'oncless' subscts have already been selected. For higher-order interactions, more general heredity styles can also be used which consider, for instance, subsets of size less than S. We do not consider these herein.
24

---

Published in Transactions on Machine Learning Research (04/2026)
Algorithm 2: Full Details of the Mode Interaction Selection Algorithm
AAUVAAAILBALLNVBRRACCOOOCCCooooO I, heredity strength r)
J 01
3
for S E P([d)) do
4
vy +- E +=
5
if ns/ S) > then
L J JU
return J TOPINTERACTIONS(J K) scores +
10
for 9 E I do
11
SCoRss()))) return (a K]
scores]S] + ScoRE(S)
GRADIENTDISSCET
θ learning rale a, epochs T)
14
for t= 1 to do
for 0 E θ do
16
8 - 8
al f + n0)
17
return +
18
NOBELVEERACTTONSEEECOOOO = 30%, 0.50, 10, = 10)
Errbest
20
I < (0}
θ (0}
22
while Error(@) < ErTbest do
Errbest Error(O)
24
25
I. TOTOLTTERACTIONS(]
28
for S E K do
Ø$ - R'S
27
θ OU θ GRADENNTDESCEET((......
28
29
return +
B.5
Gradient Descent Details
Here we provide the additional necessary implementation details of the gradient descent algorithm. Although the gradient descent algorithm is itself very simple, the bag of practical tricks required to enable scaling large event spaces quickly becomes very large. Our major innovation is the usage of higher-order Gibbs sampling which allows for much more rapid convergence of MCMC sampling, but we also find the standard usage of annealed importance sampling to provide significant gains. The upsampling of new interaction terms and caching of GPU-computed energy terms provide additional improvements
Higher-Order Block Sampling Typical Gibbs sampling constitutes resampling a single coordinate at a time by calculating the conditional distribution of a single variable given all others, (o( This has been found be extremely slow for learning higher-order energy models (Min et al.] (014)), mainly due to the inability to simultaneously activate all entries of higher-order energy parameter 6 Accordingly, instead resample according to the conditional distribution of a particular subset of the variables go(isli-s), decidedly choosing 5 which are already included in our model's collection I. We find that this provides significant speedups over the coordinate-basee approach which requires repeatedly sampling from the separate & S coordinates before becoming close to the distribution (ii)(iss).
Annealed Importance Sampling Another critical component of our framework is the usage of annealed importance sampling (Neal 2001). This allows for the approximation of the normalizing constant Og without
25

---

Published in Transactions on Machine Learning Research (04/2026)
requiring the sum over billions of elements in the distribution space. This approach also significantly benefits from our development of higher-order Gibbs sampling used during the intermediary steps. We find that this is a critical method for being able to adequately track the progress of the model over time which allows the use of our carly stopping procedure during interaction selection.
addition to sampling uniformly across the set of included interac- tions in the model, we additionally upweight those 0 parameters which have only recently been included into the model. This is inspired by the fact that the more recently included parameters are likely to move more quickly during training whereas the older parameters will have already been mostly trained and less mobile. In practice, we use 50/50 split between the old and new 0 parameters, with cach round of Gibbs sampling updating the full set of new parameters once and then updating an equal amount of the old, existing parameters.
Upsampling Active Interactions
Energy Caching In conjunction with our technique of sampling multi-dimensional conditional distribu- tions, find it efficient to avoid redundant computation of the energy functionals. In particular, suppose are updating the (0(,, (2/i212/112111 then we will split up collection follows: I Ip L In LI2 I E 1 $ 8,2 $ S} ......((( E II: 1 E 8}. In order to calculate the energies for all possible values of i,, we only need to calculate the energies once for S E Tg and only for all values of i for S E IT We find that caching these values improves the performance while summing over S T.
Practical Warnings We close with some practical warnings for reproducibility the results. It is reminded that the key challenge of doing gradient descent for an energy- based model is the need calculate the log- partition function 0 but how this translates to the gradient descent. algorithm is through poor estimates of the ne values, lcading to divergence of the training. Fortunately, although it is difficult to know a priori how much Gibbs sampling and importance sampling is required for efficient training, the quick divergence of a training run does easily identify insufficient MCMC sampling (because the convexity of the KL minimization problem ensures the good behavior of gradient descent otherwise). Lastly, it is mentioned that even after these many optimizations, the final MAHGen'Ta runs on the mushroom and adults datasets took several days on a single GPU to complete.
B.6 Verification of Heuristic
Due to the usage of J(S) as a heuristic for measuring the higher-order information in a distribution, we include additional experiments verifying that the MMI and RI values loosely track one another on most distributions. In Figure g we sample 100 random distributions for dimensions 3, 4, and 5, computing the refined information of the corresponding dimension, alongside the multiple mutual information for that distribution. Finally, scatter plot of the 100 distributions is provided in Figure Saaureraaaa and the Pearson correlation is calculated to bc: 0.8741, 0.7563, and 0.6376.
3D Distributions
Distributions
Distributions
0.36
0.25
0.14
0.08
abs(MMD
0.20
AbxMMN0
0.32
3bs(MM))
0.07
0.15
0.30
0.06
0.10
0.08
0.05
0.00
0.00
6.05
0:15
0.20
0.30
0.04
0.06
0.08
0.12 0:14
0.18
0.03
0.08
0.30
0.12 0.13
Figure 9: For 100 random distributions of sizes = 3, 4,5 (with = 5), we plot the absolute value of the multiple mutual information against the marginal refined information, finding Pearson correlation coefficients of 0.8741, 0.7563, and 0.6376.
26

---

Published in Transactions on Machine Learning Research (04/2026)
C Further Theoretical Discussion
C.1 Order Invariant Values
Since there is an abundance of complete, supported chains on the lattice of P(dd), we may still be interested in some canonical information of a mode-interact ion, without too much consideration of its supporting information set. Accordingly, let us define two canonical values of the information of a set S not associated with a specific chain or interaction set. That will be the 'marginal refined information', given by using the minimal supporting information set and the 'conditional purified information', given by using the maximal supporting information set.
Rls R = --RIlgg Tmii .s Rigond:.. RI TUggssss
I P == P([d)) 7
We note that the mutual information corresponds the marginal refined information in the case that |S] 2. In other words, MI,X,,X) =: I((t,j)) = RI (i.s) = (T -
C.2 Relation to Causal Structure Learning
We further make clear the relationship structure learning with a simple example in causal structure learning. Suppose we have the distribution induced by the causal graph depicted in Figure TO Although there is mutual information between the variables B and C, there is no direct causal information between them. In fact, they are both controlled by the variable A and the correlations between them are thereafter induced.
In the case of = and |S| = 2, the refined information only takes two values (which correspond to the marginal and conditional values). In Figure 10 we have these values calculated for all three pairs of variables. If we take look the mutual information between B and C, we can indeed see there is a positive amount of information; however, in the presence of conditioning on A, there is no refined information between these two variables.
RII0AAB,C).. 0.1308
p(a,b,
exp [oog(3) 8a,5 1og(3)))))))))) da,c
B
32
R.O.A.DC..C..C....... 0.0992
O. 1.1308
A
0
A w.p.
A w-p.
RIBAARCABBBC)..
0.0992
w.p.
C
w.p. - w.p.
B
C=
A w.p.
RIq.ABBCC)
0.0316
R.B.AB.C....A..... 0.0000
Figure 10: Simple causal graph to help illustrate refined information.
Although the traditional tools of causal structure learning, namely a set of conditional independence tests, are sufficient to identify the causal structure this case (or at least the Markov equivalence class), the tool of refined information is imagined order interactions between variables and various hypergraphical extensions to existing graphical models.
more powerful tool which may distinguish additional higher
27

---

Published in Transactions on Machine Learning Research (04/2026)
C.3
Structure of Mode Interactions
Here provide some additional figures to more quickly provide intuition about the algebraic structures introduce. In Figure 11
we depict a possible chain which is maximally refined for d = 3.
To
I1
I2
Σ3
(1,2,))
(1,233)
(1,2,3)
{1,
(1,2,3)
(1
{1,3}
{
(1,2)
{1.3}
(2.3)
{1,2)
f1.3}
{2,3)
{1,2]
{1.3]
[2,3)
X
I
(1)
(2)
{3}
(1)
(2)
{3}
(1)
{2}
{3}
(1)
(2}
[3}
I4
I5
I6
I7
(1,2,3)
(1,2,3)
(1,2,3)
(1,2,3)
(1,2)
(1,3)
(2,3)
(1,2)
(1,3)
[2,3)
(1,2)
(1.3)
[2,3]
(1,2)
[1,3]
(2,3)
(1)
(2)
(3)
(1)
(2)
3
(1)
22)
3
(1)
2
3
Figure 11: Example of a possible chain of mode interaction collections.
we depict the algebraic structure representing all possible hierarchical chains which could be selected. Each mode interaction S [d] is represented by a different color and all arrows are drawn which are possible add while still obeying the hierarchical condition. Horizontally sorted by the size of each collection I, although redundancies are suppressed (e.g. {0, (1), (2), [1,23} {0, 1, 12) is written {12}).
In Figure 12
(1)
(1
(0)
{2)
(1,3)
(1,2,3)
123
(3)
(2,3)
{12)
{3, 12)
{12, 13)
{13)
{2, 13)
{12, 23)
{12, 13, 23)
$ (1)
8 (1,2)
{23)
(1,23)
{13, 23)
S = (2)
8 = (1,3)
S = (3)
= (2,3)
Figure 12: Algebraic structure of all hierarchical collections of mode interactions.
28

---

Published in Transactions on Machine Learning Research (04/2026)
C.4 Learning Problem Formulation
In its simplest form, distribution learning is about modeling a distribution g which accurately matches the true distribution p². In this work, we use the objective of forward KL divergence which is equivalent to the maximum likelihood approach.
(((1) ((7)
(1) Dki(p"; - ( log
Evr(i) log (p) E1'()) log (a(1)),
min { ( max Er((1) log (1(1))
Importantly, as discussed throughout the paper, we reframe this objective via the choice of a sparse set of mode interactions in order to achieve better generalization properties from the learned distribution. In particular, we may write that:
DK = (KI(P:P)) + DRI(Pr:TT)
This means that, after a choice of interactions Z, our KL error decomposes orthogonally into an estimable part and an inestimable part. This means that a choice of small I will have a large inestimable error but will very accurately predict the estimable part of the distribution, whereas a choice of large I will have a small incstimable residual but will have much more challenging distribution to estimate directly.
Accordingly, may write our complete objective as the following bilevel optimization problem with outer combinatorial search over the space of I and inner continuous optimization over the Oz {0$}sez parameters:
min I
min (oxi()reee))) = min (KL((IPE) + min ( 07
0z
As discussed, the inner optimization is handled via gradient descent algorithm which leverages multiple nec- essary Monte Carlo approaches in order to learn the distribution parameters while handling the intractability of the normalizing constant. The outer combinatorial optimization is handled with simple greedy heuristic. The major benefit of using a greedy heuristic is the fact that our search across the space of I will proceed along a single hierarchical chain. It follows that the inestimable portion will be monotonically decreasing along our chain and the remaining error from the estimable part will only increase as we continue to increase the complexity of the learned distribution with our fixed and finite number of samples. This allows fit snugly within the 'classical' regime of overfiuting where have our disposal the simple rule of stopping soon validation error no longer improves.
This results in our final approach of using
argmin (DrI.Pwat ( 2,0z
where
o =rgmin (D (perm; (3))
to validate the learned distribution subset of data which is different from the training samples which are used to fit the continuous f parameters.
29


Reviews:
[{"review_id": "t4jdKuXcgT", "reviewer": "Reviewer_ZRNL", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Summary of Contributions\n\nThis manuscript introduces a novel information-geometric framework to decompose the KL error of discrete probability distributions into a sum of non-negative contrbutions termed 'refined information'. This theoretical formulation enables a procedure for progressive selection of interaction sets for hierarchical log-linear models, allowing them to capture higher-order variable interactions that are out of the scope for (fully-visible) Boltzmann machines or Markov models. The authors develop MAHGenTa, a hierarchical log-linear algorithm that utilizes a greedy heuristic based on Multiple Mutual Information to capture higher-order interactions between discrete variables. The work demonstrates the feasibility of this framework on some synthetic datasets and simple tabular real-world dataset, showing that this higher-order generative models can be used for multi-target classification tasks and perform comparably to classical binary classification baselines such as logistic regression.\n\nStrengths:\n- provides a mathematically rigorous framework that projects the data distributions into flat manifolds, ensuring that the optimization problem of minimizing the KL divergence is convex.\n- \"refined Information\" is a positive and additive metric that provides a novel decomposition of KL divergence, and allows for a more interpretable selection of interaction subsets.\n- proposes an efficient sampling strategy as an alternative to Gibbs sampling\n- the GPU implementation allows to scale the approach to to real-world tabular data, and evaluations show that the generative training procedure allows a single model to be used for multi-target classification without retraining.\n\nWeaknesses\n- the presentation of visual data is lacking, with main text figures missing axis labels and descriptive captions.\n- heavy reliance on the Appendix and reference for fundamental theoretical and implementation details.\n- several sections suffer from unclear phrasing and grammatical inconsistencies that hinder the readability of the theoretical arguments.\n\nOverall the paper is certainly interesting, but given the density of information and the page limit of TMLR, I am not quite sure if this is the right venue for this work.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The theoretical claims claims regarding the non-negativity and additivity of Refined Information, as well as the convex nature of the optimization procedure are well-supported by the mathematical definitions provided. \nThe experimental evidence on synthetic datasets effectively demonstrates the expected underfitting-overfitting behavior associated with different model complexities. \n\nHowever, the criterium of clear evidence is somewhat undermined by the lacking presentation of the theoretical framework, as well as the poor presentation of some figures, particularly figures 2 and 3, which lack any labeling of the axes, which paired with the uninformative captions, which is not self-sufficient at all, make it unnecessarily difficult to understand and interpret the figures. \nFurthermore, the comparison to existing sparse graphical models makes sense within the context of this paper, but makes it difficult to place the relevance of this work within the context of contemporary machine learning research. Maybe some benchmarks about computation time can give a sense of how feasible is it to apply aspects of this work to more large scale datasets.", "audience_yn": "Yes", "audience_explanation": "Yes, the paper is likely to interest researchers working on problems related to information geometry, statistical mechanics, and tabular data modeling. For example, the bridge between many-body expansions in physics and hierarchical log-linear models in machine learning provides a compelling perspective for those working on interpretable generative models or the modelling physical systems using machine learning. Additionally, the demonstration of using generative energy-based models for simultaneous classification tasks and bias inspection in sensitive features (e.g., race/gender) aligns with current interests in explainable AI and algorithmic fairness.", "requested_changes": "- The lack of clear labeling and captioning of the figures is a big detriment to the paper. Add clear labels to the x and y axes for Figures 2, 3, and 4. Ensure that captions provide a self-contained explanation of what is being shown in the figure (critical).\n\n- Some more detailed and rigorous evaluations are needed to contextualize the utility of the model compared to other approaches, as well as empirically demonstrate the actual utility of some of the design choices of the framework. Provide a computation time benchmarking of the model for the different datasets, and compare it to the competing approaches. Also include some ablations with regards to the more critical components such as a comparison of the benefit (whether for accuracy or computation time) of the block sampling procedure, or a comparison on some toy examples of the differences in selection of interaction collections between the theoretical RI and mutual information (critical).\n\n- clearly refer to the the definition and explanation of \"strong\" and \"weak\" heredity from the Appendix into Section 4.1 to ensure the selection algorithm is reproducible from the main text, or ideally explain it in plain words in the main text since it's currently lacking a textual explanation.\n\n- correct grammatical errors and awkward phrasing throughout the manuscript. Some examples:\n* Page 2, Section 2.1: \"modern feature selection approaches mainly concern themselves with paying the proper credence \", paying credence does not work, it is usually phrased as giving credence.\n* Page 3, Section 2.3: \"The most related work to ours is the few works that have extended the sparse graphical modeling formulation...\", awkward because of mixup of plural and singular\n* Page 10, Section 5.3: \"however, in our energy model working directly on the observed variables, it is made explicit the learned connections between variables.\", an 'in' seems to be missing.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "MpUi8om4OB", "reviewer": "Reviewer_eaNw", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "When learning probability distributions over discrete variables, log-linear models such as Boltzmann machines and Markov random fields are widely used. However, most of these models are limited to pairwise (two-body) interactions. Incorporating higher-order interactions (among three or more variables) has been avoided due to increased model complexity and the difficulty of structure learning.\nThis paper employs the framework of information geometry to quantify how much purely new information a given set of interactions contains, i.e., information that cannot be explained by lower-order interactions. The authors claim that the decrease in KL divergence (error) resulting from adding a new interaction to the model can be completely and non-redundantly decomposed into a sum of refined information contributed by each interaction.\nBased on this principle, they propose a sparse structure learning algorithm, which sequentially selects interactions that most efficiently reduce the KL error.\nBecause it is based on a theoretically guaranteed error decomposition, it provides a clear criterion for which interactions to add. Consequently, the method can efficiently represent complex data structures that conventional second-order models (e.g., Boltzmann machines) fail to capture. Experiments on both synthetic and real datasets demonstrate that the proposed method achieves higher performance with fewer samples compared to existing approaches.", "claims_and_evidence_yn": "No", "claims_explanation": "First, regarding the information geometry that forms the core of the proposed method, although the procedures themselves are described in a self-contained manner, it is difficult to conclude that the paper sufficiently explains why this approach works.\n\nIn particular, to gain geometric intuition for why the KL divergence can be decomposed into a simple sum, analogous to the Pythagorean theorem, one needs external knowledge about the dual flat structure and foliation of statistical manifolds. Without such background, the exposition feels like a conceptual leap. Especially, the idea that $\\mathcal{D}(p | \\Pi_{\\mathcal{I}_k}(p))$ can be decomposed as a sequence of orthogonal projections—the geometric concept underlying Figure 1 and Figure 2—is likely to be very difficult for first-time readers to grasp. In this sense, the current presentation seems accessible only to a very limited audience, and I cannot judge that the claims are supported by \"convincing and clear evidence\".\nI will describe possible improvements to address this issue later.\n\nIn addition, there are several other points where the exposition is unclear. These will also be discussed below.", "audience_yn": "Yes", "audience_explanation": "Inference methods for energy-based models with higher-order interactions are important, and the proposed approach is likely to be of interest to a substantial portion of the TMLR readership. In particular, the method clearly quantifies the contribution of each interaction to the error (in terms of KL divergence), and it also discusses potential applications of this perspective, which strengthens its appeal.\nEspecially, the presentation in Figure 3 is excellent. In addition, the application to fairness discussed at the end of Section 5.3 is highly interesting.", "requested_changes": "* In the paper, for computational convenience, $R_I$ (refined information) is heuristically replaced by $J_S$ (a measure of higher-order interactions). However, an explanation bridging this gap is necessary. In particular, the paper lacks a justification of this approximation: it should be explained how $J_S$ relates to local coordinate changes in information geometry, or in what sense $J_S$ provides a good approximation to $R_I$ under small perturbations. A detailed explanation of the rationale and validity of this approximation is required.\n\n* In Section 4.2, the phrase “We first reiterate the gradient of a log-linear model...” is used. Since this is the first appearance of the gradient in the paper, “reiterate” is inappropriate and should be revised.\n\n* Please provide a derivation of Equation (7).\n\n* Regarding the clarity of the paper, the description of “our hyperparameters” in Section 5.2 is insufficient. Even including the appendix, it is unclear which specific hyperparameters were swept (e.g., learning rates, regularization terms, etc.). These details should be stated explicitly.\n\n---\n\nNote: If addressing the requested changes leads to page limit issues, it may not be necessary to devote a full paragraph each to Higher-Order Block Sampling, Upsampling Active Interactions, and Annealed IS in Section 4.2. It would be sufficient to briefly state that these techniques are used, with detailed explanations deferred to the appendix.", "broader_impact_concerns": "N.A.", "additional_comments": "Regarding the non–self-consistent explanation of information geometry, the following improvements could be made. These would help increase the transparency of the theory and assist readers in developing an intuitive understanding of why the KL error can be decomposed so cleanly.\n\nWhile I do not claim that these revisions are strictly mandatory, the authors should be aware that the current presentation is not sufficiently convincing in its present form.\n\n* At the beginning of Section 3.2 or 3.3 (not in the Appendix in its current form), the Pythagorean theorem for the KL divergence should be introduced explicitly, together with an illustrative figure. In doing so, the notion of orthogonality should be clearly defined: namely, that e-flat manifolds (exponential families) and m-flat manifolds (mixture families) are orthogonal, and that this orthogonality is precisely the condition under which the Pythagorean theorem holds. This is partially discussed in Appendix C.5, but a visual and intuitive explanation in the main text would be highly desirable.\n\n* In the context of the hierarchical model proposed in this paper, it should be further explained why each projection step satisfies this orthogonality condition, based on the inclusion relationship between manifolds (i.e., $M_{I_{t-1}} \\subset M_{I_{t}}$).\n\n* The “chain of collections” introduced in Section 3.3 could be visualized not merely as a sequence of sets, but as a path in the space of probability distributions. This would allow the procedure of e-/m-projections to be explained more concretely. Beyond the algebraic structure shown in Figure 10 and 11, one could depict how the distribution moves geometrically in probability space—from the uniform distribution to the target distribution—as specific interactions are added step by step, represented as a sequence of points and projections on manifolds. In this visualization, it would also be helpful to explain that different choices of ordering (i.e., different paths) can lead to different values of refined information at each step, interpreted geometrically as differences in the “route” taken through the space.\n\n* Finally, regarding the purified gradient introduced in Section 4.2, it would strengthen the presentation to explain it not merely as a computational trick (invariance under constant shifts), but as an operation that removes redundancy in the parameter space and corresponds to a well-defined operation on the statistical manifold. This would further improve the consistency between the optimization procedure and the information-geometric framework developed in the paper.", "extra_scores": {}, "extra_text": {}}, {"review_id": "WtfMvo1Znw", "reviewer": "Reviewer_Y1SM", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper revisits hierarchical log-linear (energy-based) models for discrete/tabular distribution learning through an information-geometric lens. Specifically, the authors focus on higher-order interactions beyond standard pairwise Boltzmann machines and Markov graphical models. The main contributions of this work include: \n\n* Refined Information (RI): Introduces a notion of \"refined information\" for interaction sets, which is defined via KL differences between information-geometric projections onto nested hierarchical submanifolds. This yields a complete additive decomposition of $D\\_{KL}(p;u)$ (relative to uniform) along a chosen maximally refined hierarchical chain.\n\n* Mode Interaction Selection: Uses the KL-decomposition viewpoint to motivate a sparse interaction selection problem for hierarchical log-linear models, with an early-stopping argument to control overfitting along a monotone capacity growth process.\n\n* MAHGenTa algorithm: Proposes an end-to-end practical system combining (1) a greedy Mode Interaction Selection heuristic that uses $|J(S)|$ as a computable proxy for RI gains, and (2) GPU-based training of hierarchical log-linear models with purified gradients plus higher-order block Gibbs sampling and annealed importance sampling to handle normalization and likelihood tracking.\n\n* Empirical evaluation: Demonstrates improved KL/log-likelihood vs 1-body and 2-body baselines on synthetic data and several UCI datasets, and shows that good generative modeling yields useful downstream conditional prediction/classification performance.\n\n\nStrengths:\n\n* Clear motivation for going beyond pairwise structure; the XOR example highlights why pairwise MI can miss genuine high-order dependence.\n* The projection/KL-drop definition makes the “information attribution” nonnegative by construction, which is conceptually appealing.\n* The paper attempts to bridge theory and practice for higher-order energy-based models. Empirical results indicate consistent gains in generative fit from adding selected higher-order interactions.\n\n\nWeaknesses/limitations:\n\n* The algorithm uses $|J(S)|$ as a proxy for refined information, but the paper appears to offer limited direct evidence that this proxy reliably tracks the true RI gain in realistic settings, while the success of Mode Interaction Selection may depend critically on this correlation.\n* Comparisons are mainly against 1-body/2-body log-linear models; it is unclear how MAHGenTa compares to other scalable discrete generative baselines (e.g., autoregressive models for tabular data, modern neural density estimators, or other higher-order sparse energy-based methods).", "claims_and_evidence_yn": "Yes", "claims_explanation": "Theoretical definitions (refined information via KL differences between successive projections) are internally consistent. The conceptual claims about KL decomposition and attribution are supported. Empirical results on synthetic data and several real-world datasets show the improvement made by MAHGenTa with KL (or log-likelihood proxy) over independent and pairwise Boltzmann baselines, supporting the claim that selected higher-order interactions can improve generative performance.", "audience_yn": "Yes", "audience_explanation": "This paper studies the theoretical frameworks and practical methods for learning and understanding distributions, especially for discrete data and energy-based models, which are relevant to the TMLR community.", "requested_changes": "* Validate the MIS scoring proxy vs true refined information: It would benefit if the authors could provide experiments showing that $|J(S)|$ correlates with the actual RI gain or with true KL improvement. \n\n* Comparisons against other scalable discrete generative baselines, e.g., autoregressive models for tabular data, modern neural density estimators, or other higher-order sparse energy-based methods.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)