Efficient Learning of Discrete Graphical Models

TL;DR

This paper introduces a sample-efficient method based on Interaction Screening for learning discrete graphical models with complex interactions, providing rigorous guarantees and improving upon previous approaches in terms of sample complexity.

Contribution

The work presents the first provably sample-efficient method for learning general discrete factor models with multi-body interactions using the Interaction Screening framework.

Findings

Provides rigorous guarantees on model structure and parameter recovery.

Includes all previously studied models as special cases with improved sample complexity.

Distinguishes between model parameters and prior inputs explicitly.

Abstract

Graphical models are useful tools for describing structured high-dimensional probability distributions. Development of efficient algorithms for learning graphical models with least amount of data remains an active research topic. Reconstruction of graphical models that describe the statistics of discrete variables is a particularly challenging problem, for which the maximum likelihood approach is intractable. In this work, we provide the first sample-efficient method based on the Interaction Screening framework that allows one to provably learn fully general discrete factor models with node-specific discrete alphabets and multi-body interactions, specified in an arbitrary basis. We identify a single condition related to model parametrization that leads to rigorous guarantees on the recovery of model structure and parameters in any error norm, and is readily verifiable for a large class of models. Importantly, our bounds make explicit distinction between parameters that are proper to the model and priors used as an input to the algorithm. Finally, we show that the Interaction Screening framework includes all models previously considered in the literature as special cases, and for which our analysis shows a systematic improvement in sample complexity.

Tables

Table 1. Table 1: Sample complexity and computational complexity of Suprise over special cases.

Model name	Inter. order	Alphabet size	Recovery type	Sample complexity	Algo. complexity
Ising	2	2	structure	$O\left({\alpha }^{-4}{e}^{8\gamma }\log p\right)$	$\tilde{O}\left(p^2\right)$
Ising	2	2	${\ell }_2$-parameter	$O\left({\chi }^2{\alpha }^{-4}{e}^{8\gamma }\log p\right)$	$\tilde{O}\left(p^2\right)$
Binary	$L$	2	structure	$O\left({\alpha }^{-4}{4}^L{e}^{4\gamma L}L\log p\right)$	$\tilde{O}\left(p^L\right)$
Pairwise	2	$q$	structure	$O\left({\alpha }^{-4}{q}^4{e}^{12\gamma }\log \left(pq\right)\right)$	$\tilde{O}\left(p^2\right)$
Pairwise	2	$q$	${\ell }_2$-parameter	$O\left({\chi }^2{\alpha }^{-4}{q}^4{e}^{12\gamma }\log \left(pq\right)\right)$	$\tilde{O}\left(p^2\right)$
General	$L$	$q$	structure	$O\left({\alpha }^{-4}{q}^{2L}{e}^{4\gamma \left(L+1\right)}L\log \left(pq\right)\right)$	$\tilde{O}\left(p^L\right)$

Equations

$$\mu (\underline{\sigma })=\frac{1}{Z}\exp \left(\sum _{k\in \mathcal{K}}{\theta }_k^{\ast }{f}_k({\underline{\sigma }}_k)\right),$$

$${\varphi }_{ik}({\underline{\sigma }}_{k\setminus i}):=\frac{1}{|{\mathcal{A}}_i|}\sum _{{\sigma }_i\in {\mathcal{A}}_i}{f}_k({\underline{\sigma }}_k),$$

$$g_{ik}({\underline{\sigma }}_k)=f_k({\underline{\sigma }}_k)-{\varphi }_{ik}({\underline{\sigma }}_{k\setminus i}).$$

$$\underset{{\underline{\sigma }}_k}{\max }|g_{ik}({\underline{\sigma }}_k)|\le 1,$$

$$\Vert {\underline{\theta }}_i^{\ast }{\Vert }_1\le \hat{\gamma },$$

$${\underline{\theta }}_i^{\ast }\in {\mathcal{Y}}_i,$$

$$\Vert \hat{}{\underline{\theta }}_{{\mathcal{T}}_i}-{\underline{\theta }}_{{\mathcal{T}}_i}^{\ast }\Vert \le \frac{\alpha }{2},\forall i\in [p],$$

$${\mathcal{X}}_i=\{\underline{x}={\underline{y}}_1-{\underline{y}}_2\mid {\underline{y}}_1,{\underline{y}}_2\in {\mathcal{Y}}_i,\Vert {\underline{y}}_1{\Vert }_1\le \hat{\gamma },\Vert {\underline{y}}_2{\Vert }_1\le \hat{\gamma }\},$$

$$\text{mathbbm}E\left[{\left(\sum _{k\in {\mathcal{K}}_i}{x}_k{g}_{ik}({\underline{\sigma }}_k)\right)}^2\right]\ge {\rho }_i\Vert {\underline{x}}_{{\mathcal{T}}_i}{\Vert }^2,$$

$$\gamma =\underset{i\in \mathcal{V}}{\max }|\underset{\underline{\sigma }}{\max }\sum _{k\in {\mathcal{K}}_i}{\theta }_k^{\ast }{g}_{ik}({\underline{\sigma }}_k)|<\infty .$$

$$I_{k,k^{\prime }}$$

$${\mathcal{S}}_n({\underline{\theta }}_u)=\frac{1}{n}\sum _{t=1}^n\exp \left(-\sum _{k\in {\mathcal{K}}_u}{\theta }_k{g}_{uk}({\underline{\sigma }}_k^{(t)})\right),$$

$${\hat{\underline{\theta }}}_u=\underset{{\underline{\theta }}_u\in {\mathcal{Y}}_u:\Vert {\underline{\theta }}_u{\Vert }_1\le \hat{\gamma }}{\mathrm{argmin}}{\mathcal{S}}_n({\underline{\theta }}_u),$$

$${\mathcal{S}}_n({\hat{\underline{\theta }}}_u)\le \underset{{\underline{\theta }}_u\in {\mathcal{Y}}_u:\Vert {\underline{\theta }}_u{\Vert }_1\le \hat{\gamma }}{\min }{\mathcal{S}}_n({\underline{\theta }}_u)+ϵ.$$

$$n\ge 2^{14}\frac{{\hat{\gamma }}^2(1+\hat{\gamma }{)}^2{e}^{4\gamma }}{{\alpha }^4{\rho }_u^2}\log (\frac{4{\mathbf{K}}_u^2}{\delta }),$$

$$\sum _{k\in {\mathcal{K}}_u}{\underline{\theta }}_k{g}_{uk}({\underline{\sigma }}_k)=\sum _{k\in {\mathcal{K}}_u}\hat{}{\underline{\theta }}_k{g}_{uk}({\underline{\sigma }}_k).$$

$$\hat{}{\underline{\theta }}_u\in {\mathcal{P}}_{{\mathcal{Y}}_u}({\underline{\theta }}_u).$$

$$\mathcal{E}=\left\{(i,k)\subseteq \mathcal{V}\times \mathcal{K}\mid {\sigma }_i\in {\underline{\sigma }}_k\right\}.$$

$${\mathcal{M}}_{\mathrm{fac}}\left(\mathcal{G}\right)=\left\{k\in \mathcal{K}\mid \nexists k^{\prime }\in \mathcal{K}\text{s.t}\partial k\subset \partial k^{\prime }\right\}.$$

$${\mathcal{M}}_{\mathrm{cli}}\left(\mathcal{G}\right)=\left\{c\in P(\mathcal{V})\mid \exists k\in {\mathcal{M}}_{\mathrm{fac}}\left(\mathcal{G}\right)\text{ s.t. }c=\partial k\right\}.$$

$$\text{mathbbm}S({\mathcal{G}}^{\ast })={\mathcal{M}}_{\mathrm{cli}}\left({\mathcal{G}}^{\ast }\right).$$

$$h_k\left({\underline{\sigma }}_k\right)=f_k({\underline{\sigma }}_k)+\sum _{r\in P(\partial k)\setminus \varnothing }\frac{(-1{)}^{|r|}}{|{\mathcal{A}}_r|}\sum _{{\underline{\sigma }}_r}{f}_k({\underline{\sigma }}_k),$$

$${\rho }_i^{\mathrm{NPC}}=\underset{\begin{array}{c}c\in {\mathcal{M}}_{\mathrm{cli}}\left(\mathcal{G}\right)\\ c\ni i\end{array}}{\min }\underset{\begin{array}{c}\Vert {\underline{x}}_c{\Vert }_2=1\\ {\underline{x}}_c\in {\mathcal{X}}_i^c\end{array}}{\min }{\text{mathbbm}E}_{\left({\sigma }_i\right)}\left[\sum _{{\underline{\sigma }}_{c\setminus i}\in {\mathcal{A}}_{c\setminus i}}{\left(\sum _{k\in {\left[c\right]}_{\mathrm{sp}}}{x}_k{h}_k({\underline{\sigma }}_k)\right)}^2\right],$$

$$\text{mathbbm}E\left[{\left(\sum _{k\in {\mathcal{K}}_i}{x}_k{g}_{ik}({\underline{\sigma }}_k)\right)}^2\right]\ge {\rho }_i^{\mathrm{NPC}}{\left(\frac{\exp (-2\gamma )}{q}\right)}^{L-1}\Vert {\underline{x}}_{{\mathcal{T}}_i}{\Vert }_{\infty ,2}^2.$$

$$\text{mathbbm}E\left[{\left(\sum _{k\in {\mathcal{K}}_i}{x}_k{g}_{ik}({\underline{\sigma }}_k)\right)}^2\right]\ge \frac{{\rho }_i^{\mathrm{NPC}}}{\chi }\frac{\exp (-2\gamma )}{q}\Vert {\underline{x}}_{{\mathcal{T}}_i}{\Vert }_2^2.$$

$$n\ge 2^{14}{q}^{2\left(L-1\right)}\frac{{\hat{\gamma }}^2(1+\hat{\gamma }{)}^2{e}^{4\gamma L}}{{\alpha }^4{\rho }_{\mathrm{NPC}}^2}\log \left(\frac{4pL{\mathbf{K}}^2}{\delta }\right),$$

$$\exp \left(\sum _{k\in {\mathcal{K}}_i}{\theta }_k^{\ast }{f}_k({\underline{\sigma }}_k)\right)$$

$${\mu }_{\underline{\theta }}(\sigma ,s)\propto \exp \left({\theta }_1\sigma (s-1)+{\theta }_{2}s(\sigma -1)+{\theta }_3(\sigma +s)\right),$$

$${\mu }_{\underline{\eta }}(\sigma ,s)\propto \exp \left({\eta }_1\sigma s+{\eta }_2\sigma +{\eta }_{3}s\right).$$

$$\text{mathbbm}P(\sigma \mid s)$$

Full text

1 Introduction

Representing and understanding the structure of direct correlations between distinct random variables with graphical models is a fundamental task that is essential to scientific and engineering endeavors. It is the first step towards an understanding of interactions between interleaved constituents of elaborated systems [10]; it is key for developing causal theories [6]; and it is at the core of automated decision making [8], cybersecurity [5] and artificial intelligence [19].

The problem of reconstruction of graphical models from samples traces back to the seminal work of [7] for tree-structured graphical models, and as of today is still at the center of attention of the learning community. For factor models defined over general hypergraphs, the learning problem is particularly challenging in graphical models over discrete variables, for which the maximum likelihood estimator is in general computationally intractable. One of the earlier tractable algorithms that has been suggested to provably reconstruct the structure of a subset of pairwise binary graphical models is based on inferring the sparsity pattern of the so-called regularized pseudo-likelihood estimator, equivalent to regularized logistic regression in the binary case [16]. However, additional assumptions required for this algorithm to succeed severely limit the set of pairwise binary models that can be learned [15]. After it was proven that reconstruction of any discrete graphical models with bounded degree can be done in polynomial time in the system size [4], Bresler showed that it is possible to bring the computational complexity down to quasi-quadratic in the number of variables for Ising models (pairwise graphical models over binary variables); however, the resulting algorithm has non-optimal sample requirements that are double-exponential in other model parameters [3]. The first computationally efficient reconstruction algorithm for sparse pairwise binary graphical models with a near-optimal sample complexity with respect to the information theoretic lower bound [17], called Rise, was designed and analyzed in [18]. The algorithm Rise suggested in this work is based on the minimization of a novel local convex loss function, called the Interaction Screening objective, supplemented with an $\ell_{1}$ penalty to promote sparsity. Even though it has been later shown in [13] that regularized pseudo-likelihood supplemented with a crucial post-processing step also leads to a structure estimator for pairwise binary models, strong numerical and theoretical evidence provided in that work demonstrated that Rise is superior in terms of worst-case sample complexity.

Algorithms for learning discrete graphical models beyond pairwise and binary alphabets have been proposed only recently in [9] and [12]. The method in [9] works for arbitrary models with bounded degrees, but being a generalization of Bresler’s algorithm for Ising models [3], it suffers from similar prohibitive sample requirements growing double-exponentially in the strength of model parameters. The so-called Sparsitron algorithm in [12] has the flavor of a stochastic first order method with multiplicative updates. It has a low computational complexity and is sample-efficient for structure recovery of two subclasses of discrete graphical models: multiwise graphical models over binary variables or pairwise models with general alphabets. A recent follow-up work [20] considered an $\ell_{2,1}$ constrained logistic regression, and showed that it provides a slight improvement of the sample complexity compared to [12] in the case of pairwise models over non-binary variables.

In this work, we propose a general framework for learning general discrete factor models expressed in an arbitrary parametric form. Our estimator termed Grise is based on a significant generalization of the Interaction Screening method of [18, 13], previously introduced for pairwise binary models. Our primary insight lies in the identification of a single general condition related to model parameterization that is sufficient to obtain bounds on sample complexity. We show that this condition can be reduced to a set of local identifiability conditions that only depend on the size of the maximal clique of the factor graph and can be explicitly verified in an efficient way. We propose an iterative algorithm called Suprise which is based on Grise and show that it can efficiently perform structure and parameter estimation for arbitrary graphical models. Existing results in the literature on this topic [18, 9, 12, 20] can be obtained as special cases of our general reconstruction results, which noticeably includes the challenging case of multi-body interactions defined over general discrete alphabets. Our theoretical guarantees can be expressed in any error norm, and explicitly includes distinction between bounds on the parameters of the underlying model and the prior parameters used in the optimization; as a result prior information that is not tight only has moderate effect on the sample complexity bounds. Finally, we also provide a fully parallelizable algorithmic formulation for the Grise estimator and Suprise algorithm, and show that they have efficient run times of $\widetilde{O}(p^{L})$ for a model of size $p$ with $L$-order interactions, that includes the best-known $\widetilde{O}(p^{2})$ scaling for pairwise models.

2 Problem formulation

In this Section, we formulate the general discrete graphical model selection problem that we consider and describe conditions that makes this problem well-posed.

2.1 Parameterized family of models

We consider positive joint probability distributions over $p$ variables $\sigma_{i}\in\mathcal{A}_{i}$ for $i=1,\ldots,p$. The set of variable indices $i$ is referred to as vertices $\mathcal{V}=1,\ldots,p$. Node-dependent alphabets $\mathcal{A}_{i}$ are assumed to be discrete and of size bounded by $q>0$. Without loss of generality, the positive probability distribution over the $p$-dimensional vector $\underline{\sigma}$ can be expressed as

[TABLE]

where $\{f_{k},\ k\in\mathcal{K}\}$ is a set of basis functions acting upon subsets of variables $\underline{\sigma}_{k}\subseteq\underline{\sigma}$ that specify a family of distributions and $\theta^{*}_{k}$ are parameters that specify a model within this family. The quantity $Z$ denotes the partition function and serves as a normalization constant that enforces that the $\mu$ in (1) is a probability distribution. For $i\in\{1,\ldots,p\}$, let $\mathcal{K}_{i}\subseteq\mathcal{K}$ denote the set of factors corresponding to basis functions acting upon subsets $\underline{\sigma}_{k}$ that contain the variable $\sigma_{i}$ and $|\mathcal{K}_{i}|=\mathbf{K}_{i}$.

Given any set of basis functions, we can locally center them by first defining for a given $i\in[p]$, the local centering functions

[TABLE]

where $\underline{\sigma}_{k\setminus i}$ denotes the vector $\underline{\sigma}_{k}$ without $\sigma_{i}$, and define the locally centered basis functions,

[TABLE]

As their name suggests, the locally centered basis functions sum to zero $\sum_{\sigma_{i}\in\mathcal{A}_{i}}g_{ik}(\underline{\sigma}_{k})=0$. To ensure the scales of the parameters are well defined, we assume that $\theta^{*}_{k}$ are chosen or rescaled such that all locally centered basis functions are normalized in the following sense:

[TABLE]

for all vertices $i\in\mathcal{V}$ and basis factor $k\in\mathcal{K}_{i}$. This normalization can always be achieved by choosing bounded basis functions $|f_{k}(\underline{\sigma}_{k})|\leq 1/2$. An important special case is when the basis functions are already centered, i.e. $g_{ik}(\underline{\sigma}_{k})=f_{k}(\underline{\sigma}_{k})$. In this case the basis functions are directly normalized $\max_{\underline{\sigma}_{k}}|f_{k}(\underline{\sigma}_{k})|=1$. Note that one of the reasons to define the normalization in (4) in terms of the centered functions $g_{k}$ instead of $f_{k}$ is to avoid spurious cases where the functions $f_{k}$ have inflated magnitudes due to addition of constants $f_{k}\leftarrow f_{k}+C$. In Appendix A, we show that the other important reason to employ centered functions is that degeneracy of the local parameterization with these functions translates to degeneracy of the global distribution in Eq. (1).

2.2 Model selection problem

For each $i\in[p]$, let $\mathcal{T}_{i}\subseteq\mathcal{K}_{i}$ denote the set of target factors that we aim at reconstructing accurately and let $\mathcal{R}_{i}=\mathcal{K}_{i}\setminus\mathcal{T}_{i}$ be the set of residual factors for which we do not need learning guarantees. The target and residual parameters are defined similarly as $\underline{\theta}^{*}_{\mathcal{T}_{i}}=\{\theta_{k}^{*}\mid k\in\mathcal{T}_{i}\}$ and $\underline{\theta}^{*}_{\mathcal{R}_{i}}=\{\theta_{k}^{*}\mid k\in\mathcal{R}_{i}\}$ respectively. Given independent samples from a model in the family in Section 2.1, the goal of the model selection problem is to reconstruct the target parameters of the model.

Definition 1** (Model Selection Problem).**

Given $n$ i.i.d. samples $\underline{\sigma}^{(1)},\ldots,\underline{\sigma}^{(n)}$ drawn from some distribution $\mu(\underline{\sigma})$ in Eq. (1) defined by $\underline{\theta}^{*}$, and prior information on $\underline{\theta}^{*}$ given in form of an upper bound on the $\ell_{1}$-norm of the local sub-components

[TABLE]

and a local constraint set $\mathcal{Y}_{i}\subseteq\mathbbm{R}^{\mathbf{K}_{i}}$ for each $i\in[p]$ such that

[TABLE]

compute estimates $\widehat{}\underline{\theta}$ of $\underline{\theta}^{*}$ such that the estimates of the target parameters satisfy

[TABLE]

where $\|\cdot\|$ denotes some norm of interest with respect to which the error is measured.

The bound on the $\ell_{1}$-norm in (5) is a natural generalization of the sparse case where $\underline{\theta}^{*}$ only has a small number of non-zero components; in the context of parameter estimation in graphical models, the setting of parameters bounded in the $\ell_{1}$-norm has been previously considered in [12]. The constraint sets $\mathcal{Y}_{i}$ are used to encode any other side information that may be known about the model.

2.3 Sufficient conditions for well-posedness

We describe some conditions on the model in (1) that makes the model selection problem in Definition 1 well-posed. We first state the conditions formally.

Condition 1**.**

The model from which the samples are drawn in the model selection problem in Definition 1 satisfies the following:

• (C1)

Local Learnability Condition for Graphical Models:* There exists a constant $\rho_{i}>0$ such that for every vertex $i$ and any vector in the perturbation set $\underline{x}\in\mathcal{X}_{i}\subseteq\mathbbm{R}^{\mathbf{K}_{i}}$ defined as*

[TABLE]

the following holds:

[TABLE]

where $\underline{x}_{\mathcal{T}_{i}}$ denotes the components $k\in\mathcal{T}_{i}$ of $x$, and $\|\cdot\|$ is the norm used in Definition 1.

• (C2)

Finite Maximum Interaction Strength:* The following quantity $\gamma$ is finite,*

[TABLE]

Condition (C1) consists in satisfying the inequality in Eq. (9) involving a quadratic form $\underline{x}^{\top}\widetilde{I}\underline{x}$ where the matrix $\widetilde{I}$ has indices $k,k^{\prime}\in\mathcal{K}_{i}$ and is explicitly defined as $\widetilde{I}_{k,k^{\prime}}={\mathbbm{E}}\left[g_{ik}(\underline{\sigma}_{k})g_{ik^{\prime}}(\underline{\sigma}_{k^{\prime}})\right]$. This matrix $\widetilde{I}$ is in fact related to the conditional Fisher information matrix.

The conditional Fisher information matrix $I$ with indices $k,k^{\prime}\in\mathcal{K}_{i}$ is derived from the conditional distribution of $\sigma_{i}$ given the remaining variables and reads,

[TABLE]

We immediately see that the matrix $\widetilde{I}$ dominates the conditional Fisher information matrix in the positive semi-definite sense, that is $\underline{x}^{\top}\widetilde{I}(\theta^{*})\underline{x}\geq\underline{x}^{\top}I(\theta^{*})\underline{x}$ for all $\underline{x}\in\mathbbm{R}^{\mathbf{K}_{i}}$. Therefore, Condition (C1) is satisfied whenever the conditional Fisher information matrix is non-singular in the parameter subspace $\underline{x}_{\mathcal{T}_{i}}$ that we care to reconstruct and which is compatible with our priors, i.e. for $\underline{x}\in\mathcal{X}_{i}$. We would like to add that the conditional Fisher information matrix is a natural quantity to consider in this problem as we deliberately focus on using conditional statistics rather than global ones in order to bypass the intractability of the global log-likelihood approach. We are strongly convinced that it should appear in any analysis that entails conditional statistics.

Condition (C2) is required to ensure that the model can be recovered with finitely many samples. For many special cases, such as the Ising model, the minimum number of samples required to estimate the parameters must grow exponentially with the maximum interaction strength [17]. A more detailed discussion about well-posedness and Conditions (C1) and (C2) can be found in Appendix A.

Conditions (C1) and (C2) differ from the concepts in [18] called restricted strong convexity property and bound on the interaction strength, respectively, in a subtle but critical manner. Conditions (C1) can be identified with restricted strong convexity only when the $\ell_{2}$-norm is used in Eq. (9). We will see later that the notion of restricted strong convexity is not required for the $\ell_{\infty}$-norm that appears to be a natural metric for which the local learnability condition can be verified for general models. Moreover, for general models it remains unclear whether the restricted strong convexity holds for values of $\rho_{i}$ that are independent of the problem dimension $p$. Condition (C2) is a weaker assumption than the bound on the interaction strength from [17] for it does not require an extra assumption on the maximum degree of the graphical model.

3 Generalized interaction screening

In this Section, we introduce the algorithm that efficiently solves the model selection problem in Definition 1 and provides rigorous guarantees on its reconstruction error and computational complexity.

3.1 Generalized regularized interaction screening estimator

We propose a generalization of the estimator Rise, first introduced in [18] for pairwise binary graphical models, in order to reconstruct general discrete graphical models defined in (1). The generalized interaction screening objective (GISO) is defined for each vertex $u$ separately and is given by

[TABLE]

where $\underline{\sigma}^{(1)},\ldots,\underline{\sigma}^{(n)}$ are $n$ i.i.d samples drawn from $\mu(\underline{\sigma})$ in Eq. (1), $\underline{\theta}_{u}:=(\underline{\theta}_{k})_{k\in\mathcal{K}_{u}}$ is the vector of parameters associated with the factors in $\mathcal{K}_{u}$ and the locally centered basis functions $g_{uk}$ are defined as in Eq. (3). The GISO retains the main feature of the interaction screening objective (ISO) in [18]: it is proportional to the inverse of the factor in $\mu(\underline{\sigma})$, except for the additional centering terms $\phi_{uk}$. The GISO is a convex function of $\underline{\theta}_{u}$ and retains the “screening” property of the original ISO. The GISO is used to define the generalized regularized interaction screening estimator (Grise) for the parameters given by

[TABLE]

where $\widehat{\gamma}$ and $\mathcal{Y}_{u}$ are the prior information available on $\underline{\theta}_{u}^{*}$ as defined in (5) and (6).

3.2 Error bound on parameter estimation with Grise

We now state our first main result regarding the theoretical guarantees on the parameters reconstructed by Grise. We call $\widehat{}\underline{\theta}_{u}$ an $\epsilon$-optimal solution of (13) if

[TABLE]

Theorem 1** (Error Bound on Grise Estimates).**

Let $\underline{\sigma}^{(1)},\ldots,\underline{\sigma}^{(n)}$ be i.i.d. samples drawn according to $\mu(\underline{\sigma})$ in (1). For some node $u\in\mathcal{V}$, assume that the model satisfies Condition 1 for some norm $\|\cdot\|$ and some constraint set $\mathcal{Y}_{u}$, and let $\alpha>0$ be the prescribed accuracy level. If the number of samples satisfies

[TABLE]

then, with probability at least $1-\delta$, any estimate that is an $\epsilon$-minimizer of Grise, with $\epsilon\leq(\rho_{u}\alpha^{2}e^{-\gamma})/(20(1+\widehat{\gamma}))$, satisfies $\|\widehat{\underline{\theta}}_{\mathcal{T}_{u}}-\underline{\theta}^{*}_{\mathcal{T}_{u}}\|\leq\frac{\alpha}{2}$.

The proof of Theorem 1 can be found in Appendix B.

The computational complexity of finding an $\epsilon$-optimal solution of Grise for a trivial constraint set $\mathcal{Y}_{u}=\mathbbm{R}^{\mathbf{K}_{u}}$ is $C\frac{c_{g}n\mathbf{K}_{u}}{\epsilon^{2}}\ln(1+\mathbf{K}_{u})$, where $c_{g}$ is an upper-bound on the computational complexity of evaluating any $g_{ik}(\underline{\sigma}_{k})$ for $k\in\mathcal{K}_{i}$, and $C$ is a universal constant independent of all the parameters of the problem, see Proposition 5 in Appendix C. For a certain class of constraint sets $\mathcal{Y}_{u}$, which we term parametrically complete, the problem can be solved in two steps: first, finding a solution to an unconstrained problem, and then projecting onto this set. Note, however, that in general the problem of finding $\epsilon$-optimal solutions to constrained Grise can still be difficult since the constraint set $\mathcal{Y}_{u}$ can be arbitrarily complicated.

Definition 2**.**

The constraint set $\mathcal{Y}_{u}$ is called a parametrically complete set if for all $\underline{\theta}_{u}\in\mathbbm{R}^{|\mathbf{K}_{u}|}$, there exists $\widehat{}\underline{\theta}_{u}\in\mathcal{Y}_{u}$ such that for all $\underline{\sigma}_{u}$, we have

[TABLE]

Any $\widehat{}\underline{\theta}_{k}\in\mathcal{Y}_{u}$ satisfying (16) is called an equi-cost projection of $\underline{\theta}_{u}$ onto $\mathcal{Y}_{u}$ and is denoted by

[TABLE]

The computational complexity of finding of an $\epsilon$-optimal solution of GRISE with parametrically complete set is $C\frac{c_{g}n\mathbf{K}_{u}}{\epsilon^{2}}\ln(1+\mathbf{K}_{u})+\mathcal{C}(\mathcal{P}_{\mathcal{Y}_{u}}(\widehat{}\underline{\theta}_{u}^{\text{unc}}))$, where $\mathcal{C}(\mathcal{P}_{\mathcal{Y}_{u}}(\widehat{}\underline{\theta}_{u}^{\text{unc}}))$ denotes the computational complexity of the projection step, see Theorem 3 in Appendix C.

As we will see, for many graphical models it is often possible to explicitly construct parametrically complete sets for which the computational complexity of the projection step $\mathcal{C}(\mathcal{P}_{\mathcal{Y}_{u}}(\widehat{}\underline{\theta}_{u}^{\text{unc}}))$ is insignificant compared to the computational complexity of unconstrained Grise.

4 Structure identification and parameter estimation

In this Section we show that the structure of graphical models, which is the collection of maximal subsets of variables that are associated through basis functions, as well as the associated parameters, can be efficiently recovered. The key elements are twofold. First, we prove that for maximal cliques, the Local Learnability Condition (LLC) in (C1) can be easily verified and yields a LLC constant independent of the system size. Second, we leverage this property to design an efficient structure and parameter learning algorithm coined Suprise that requires iterative calls of Grise.

4.1 The structure of graphical models

The structure plays a central role in graphical model learning for it contains all the information about the conditional independence or Markov property of the distribution $\mu(\underline{\sigma})$ from Eq (1). In order to reach the definition of the structure presented in Eq. (21), we have to introduce graph theoretic concepts specific to graphical models.

The factor graph associated with the model family is a bipartite graph $\mathcal{G}=\left(\mathcal{V},\mathcal{K},\mathcal{E}\right)$ with vertex set $\mathcal{V}$, factor set $\mathcal{K}$ and edges connecting factors and vertices,

[TABLE]

We see from Eq. (18) that the edge $(i,k)$ exists when the variable $\sigma_{i}$ associated with the vertex $i$ is an argument of the basis function $f_{k}(\underline{\sigma}_{k})$ associated with the factor $k$. Note that this definition only depends on the set of basis functions $\mathcal{K}$ and does not refer to a particular choice of model within the family. The factor graph $\mathcal{G}^{*}=\left(\mathcal{V},\mathcal{K}^{*},\mathcal{E}^{*}\right)$ associated with a model, as defined in Eq. (1), is the induced subgraph of $\mathcal{G}$ obtained from the vertex set $\mathcal{V}$ and factor subset $\mathcal{K}^{*}=\left\{k\in\mathcal{K}\mid\theta_{k}^{*}\neq 0\right\}$. We also use the shorthand notation $\mathcal{G}^{*}=\mathcal{G}\left[(\mathcal{V},\mathcal{K}^{*})\right]$ to denote an induced subgraph of $\mathcal{G}$.

We define the neighbors of a factor $k$ as the set of vertices linked by an edge to $k$ and denote it by $\partial k=\left\{i\in\mathcal{V}\mid(i,k)\in\mathcal{E}\right\}$. The largest factor neighborhood size $L=\max_{k\in\mathcal{K}}\left|\partial k\right|$ is called the interaction order. Families of graphical models with $L\leq 2$ are referred to as pairwise models as opposed to the generic denomination of $L$-wise models when $L$ is expected to be arbitrary.

The set of maximal factors of a graph is the set of factors whose neighborhood is not strictly contained in the neighborhood of another factor,

[TABLE]

Notice that multiple maximal factors may have the same neighborhood. This motivates the definition of the set of maximal cliques which is contained in the powerset $P(\mathcal{V})$ and consists of all neighborhoods of maximal factors,

[TABLE]

The set of factors whose neighborhoods are the same maximal clique $c\in\mathcal{M}_{\rm{cli}}\left(\mathcal{G}\right)$ is called the span of the clique defined as $\left[c\right]_{\rm{sp}}=\left\{k\in\mathcal{M}_{\rm{fac}}\left(\mathcal{G}\right)\mid c=\partial k\right\}$. Finally, the structure $\mathbbm{S}$ of a graphical model is the set of maximal cliques associated with the factor graph of the model,

[TABLE]

We would like to stress that the structure of a model is different from the set of maximal cliques of the family of graphical models $\mathcal{M}_{\rm{cli}}\left(\mathcal{G}\right)$ as the former is constructed with the set of factors associated with non-zero parameters while the latter consists of all potential maximal factors.

4.2 From local learnability condition to nonsingular parametrization of cliques

We show that the learning problem of reconstructing maximal cliques is well-posed in general and especially for non-degenerate globally centered basis functions. To this end, we demonstrate that the LLC in (C1) is automatically satisfied whenever the target sets $\mathcal{T}_{i}$ consist of factors corresponding to maximal cliques of the graphical model family. Importantly, we prove that the LLC constant $\rho_{i}$ does not depend on the dimension of the model for the $\ell_{\infty,2}$-norm but rather relies on the Nonsingular Parametrization of Clique (NPC) by the basis functions. Similarly, we also guarantee that the LLC holds for the $\ell_{2}$-norm in the case of pairwise colorable models.

We introduce globally centered basis functions defined for any factor $k\in\mathcal{K}$ through the inclusion–exclusion formula,

[TABLE]

where $\mathcal{A}_{r}=\bigotimes_{j\in r}\mathcal{A}_{j}$. It is straightforward to see that globally centered functions sum partially to zero for any variables, i.e. $\sum_{\sigma_{i}\in\mathcal{A}_{i}}h_{k}(\underline{\sigma}_{k})=0$ for all $i\in\partial k$. It is worth noticing that when the functions $f_{k}$ are already globally centered, we have $f=g=h$. We would also like to point out that unlike locally centered functions $g_{ik}$, globally centered functions cannot in general be interchanged with functions $f_{k}$ without modifying conditional distributions. However they play an important role in determining the independence of basis functions around cliques through the Nonsingular Parametrization of Cliques (NPC) constant introduced below. Given a perturbation set $\mathcal{X}_{i}$, as defined in Eq. (8), the NPC constant is defined through the following minimization,

[TABLE]

where the vector $\underline{x}_{c}\in\mathbbm{R}^{|\left[c\right]_{\rm{sp}}|}$ belongs to $\mathcal{X}^{c}_{i}$, the projection of the constraint $\mathcal{X}_{i}\subseteq\mathbbm{R}^{\mathbf{K}_{i}}$ to the components $k\in\left[c\right]_{\rm{sp}}$ and the expectation is with respect to the marginal distribution of $\sigma_{i}$. Note that NPC constant only depends on $L$ and not on the size of the system, and can be explicitly computed in time $O(\mathbf{K})$. A detailed discussion can be found in Appendix D. The importance of the NPC constant is highlighted by the following proposition that guarantees that the LLC is satisfied for maximal factors in $\ell_{\infty,2}$-norm as long as $\rho^{\rm{NPC}}_{i}>0$.

Proposition 1** (LLC in $\ell_{\infty,2}$-norm).**

For a specific vertex $i\in\mathcal{V}$, let the target set be maximal factors with $i$ as neighbor, i.e. $\mathcal{T}_{i}=\left\{k\in\mathcal{M}_{\rm{fac}}\left(\mathcal{G}\right)\mid\partial k\ni i\right\}$. For vectors in the perturbation set $\underline{x}\in\mathcal{X}_{i}\subseteq\mathbbm{R}^{\mathbf{K}_{i}}$, define the $\ell_{\infty,2}$-norm over components that are maximal factors as $\|\underline{x}_{\mathcal{T}_{i}}\|_{\infty,2}=\max_{\begin{subarray}{c}c\in\mathcal{M}_{\rm{cli}}\left(\mathcal{G}\right)\\ c\ni i\end{subarray}}\sqrt{\sum_{k\in\left[c\right]_{\rm{sp}}}x^{2}_{k}}$. Then for discrete graphical models with maximum alphabet size $q$, interaction order $L$ and models with finite maximum interaction strength $\gamma$ as defined in Eq. (10), the Local Learnability Condition (C1) is satisfied whenever the Nonsingular Parameterization of Cliques constant $\rho^{\rm{NPC}}_{i}$ is nonzero and we have,

[TABLE]

Proposition 1 guarantees that the LLC can be satisfied uniformly in the size $p$ of the model whenever $\rho_{i}^{\rm{NPC}}>0$. The proof of Proposition 1 can be found in Appendix D.

For family of models whose factors involve at most $L=2$ variables, the so-called pairwise models, we can show that the LLC conditions for maximal factors also holds for the $\ell_{2}$-norm. This LLC conditions depends on the vertex chromatic number $\chi$ of the model factor graph. We recall that a vertex coloring of a graph $\mathcal{G}^{*}=\left(\mathcal{V},\mathcal{K}^{*},\mathcal{E}^{*}\right)$ is a partition $\{S_{r}\}_{r\in\mathbb{N}}\in P(\mathcal{V})$ of the vertex set such that no two vertices with the same color are connected to the same factor node, i.e. $i,j\in S_{r}\Rightarrow\nexists k\in\mathcal{K}^{*}$ s.t. $i,j\in\partial k$. The chromatic number is the cardinality of the smallest graph coloring.

Proposition 2** (LLC in $\ell_{2}$-norm for pairwise models).**

For a specific vertex $i\in\mathcal{V}$, let the target set be maximal factors with $i$ as neighbor, i.e. $\mathcal{T}_{i}=\left\{k\in\mathcal{M}_{\rm{fac}}\left(\mathcal{G}\right)\mid\partial k\ni i\right\}$. For $\underline{x}\in\mathcal{X}_{i}\subseteq\mathbbm{R}^{\mathbf{K}_{i}}$, define the $\ell_{2}$-norm over components that are maximal factors $\|\underline{x}_{\mathcal{T}_{i}}\|_{2}=\sqrt{\sum_{k\in\mathcal{T}_{i}}x^{2}_{k}}$. Then for discrete pairwise graphical models with maximum alphabet size $q$ and models with chromatic number $\chi$ and finite maximum interaction strength $\gamma$ as defined in Eq. (10), the Local Learnability Condition (C1) is satisfied whenever the NPC constant $\rho^{\rm{NPC}}_{i}$ is nonzero and we have,

[TABLE]

The reader will find the proof of Proposition 2 in Appendix D.

4.3 Structure unveiling and parameter reconstruction with interaction screening estimation

Suppose that we know $\alpha>0$, a lower-bound on the minimal intensity of the parameters associated with the structure in the sense that $\alpha\leq\min_{c\in\mathbbm{S}(\mathcal{G}^{*})}\sqrt{\sum_{k\in\left[c\right]_{\rm{sp}}}{\theta^{*}_{k}}^{2}}$. Then we can recover the structure and parameters associated with maximal factors of any graphical models using Algorithm 1, coined Suprise for Structure Unveiling and Parameter Reconstruction with Interaction Screening Estimation. Suprise that implements an iterative use of Grise is shown to have a sample complexity logarithmic in the system size for models with non-zero NPC constants. Our second main result is the following Theorem 2, proved in Appendix D, which provides guarantees on Suprise.

Theorem 2** (Reconstruction and Estimation Guarantees for Suprise).**

Let $\mu(\underline{\sigma})$ in (1) be the probability distribution of a discrete graphical model with maximum alphabet size $q$, interaction order $L$, finite maximum interaction strength $\gamma$ and smallest Nonsingular Parameterization of Cliques constant greater than zero, i.e. $\rho_{\rm{NPC}}=\min_{u\in\mathcal{V}}\rho^{\rm{NPC}}_{u}>0$. Let $\underline{\sigma}^{(1)},\ldots,\underline{\sigma}^{(n)}$ be i.i.d. samples drawn according to $\mu(\underline{\sigma})$ and assume that

[TABLE]

where $\mathbf{K}=\max_{u\in\mathcal{V}}\mathbf{K}_{u}$ is the maximal number of basis functions in which a variable can appear and $\widehat{\gamma}\geq\max_{u\in\mathcal{V}}\|\underline{\theta}^{*}_{u}\|_{1}$ is our $\ell_{1}$-prior on the components of the parameters. Then the structure of the general graphical model is perfectly recovered using Algorithm 1, i.e. $\widehat{\mathbbm{S}}=\mathbbm{S}$, with probability $1-\delta$. In addition, the parameters associated with maximal factors are reconstructed with precision $\max_{c\in\mathbbm{S}}\sum_{k\in\left[c\right]_{\rm{sp}}}\left(\widehat{\theta}_{k}-\theta^{*}_{k}\right)^{2}\leq\alpha^{2}/4$ for general models and with $\sum_{c\in\mathbbm{S}}\sum_{k\in\left[c\right]_{\rm{sp}}}\left(\widehat{\theta}_{k}-\theta^{*}_{k}\right)^{2}\leq\chi^{2}\alpha^{2}/4$ for pairwise models with chromatic number $\chi$.

The total computational complexity scales as $\widetilde{\mathcal{O}}(p\mathbf{K})$, for fixed $L$, $\alpha$, $\gamma$, $\widehat{\gamma}$ and $\delta$, if the constraint sets $\mathcal{Y}_{u}$ are parametrically complete.

As an application of Theorem 2, we show the sample and computational complexity of recovering parameter values and the structure of some well-known special cases in Table 1. The parameter $\alpha$ appearing in Table 1 is the precision to which parameters are recovered in the considered norm, $\chi$ is the chromatic number of the graph, $L$ is the interaction order, $q$ is the maximum alphabet size, $\gamma$ is the maximum interaction strength and $p$ is the number of variables. At this point, it is instructive to compare our sample complexity requirements to existing results. A direct application of bounds of [12] and [20] to the case of pairwise multi-alphabet models that we consider below yields $O(\exp(14\gamma))$ dependence, whereas Suprise has a complexity that scales as $O(\exp(12\gamma))$. In the case of binary $L$-wise models, while [12] shows the $O(\exp(O(\gamma L)))$ scaling, Suprise enjoys a sample complexity $O(\exp(4\gamma L))$. The algorithm of [9] recovers a subclass of general graphical models with bounded degree, but has a sub-optimal double-exponential scaling in $\gamma$, while Suprise leads to recovery of arbitrary discrete graphical models with a single-exponential dependence in $\gamma$ and needs no bound on the degree. In terms of the computational complexity, Suprise achieves the efficient scaling $\widetilde{O}(p^{L})$ for models with the maximum interaction order $L$, which matches the best-known $\widetilde{O}(p^{2})$ scaling for pairwise models [12, 20]. In summary, Suprise generalizes, improves and extends the existing results in the literature. The proofs for special cases can be found in Section E of the Supplementary Material.

5 Conclusion and future work

A key result of our paper is the existence of a computationally efficient algorithm that is able to recover arbitrary discrete graphical models with multi-body interactions. This result is a particular case of the general framework that we have introduced, which considers arbitrary model parametrization and makes distinction between the bounds on the parameters of the underlying model and the prior parameters. The computational complexity $\widetilde{O}(p^{L})$ that we achieve is believed to be efficient for this problem [12]. In terms of sample complexity, the information-theoretic bounds for recovery of general discrete graphical models are unknown. In the case of binary pairwise models, the sample complexity bounds resulting from our general analysis are near-optimal with respect to known information-theoretic lower bounds [17]. It would be interesting to see if the $1/\alpha^{4}$ factor in our sample complexity bounds can be improved to $1/\alpha^{2}$ using an $\ell_{1}$-norm penalty rather than an $\ell_{1}$-norm constraint, as it has been shown for the particular case of Ising models [13, 18].

Other open questions left for future exploration include the possibility to extend the analysis to the case of graphical models with nonlinear parametrizations like in [11], and to graphical models with continuous variables. It is particularly interesting to see whether the computationally efficient and nearly sample-optimal method introduced in the present work could be useful for designing efficient learning algorithms that can improve the state-of-the-art in the well-studied case of Gaussian graphical models, for which it has been recently shown that the information-theoretic lower bound on sample complexity is tight [14].

Acknowledgments

Research presented in this article was supported by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project numbers 20190059DR, 20190195ER, 20190351ER, and 20210078DR.

Appendices

Appendix A contains a detailed discussion about Conditions (C1) and (C2). In Appendix B, the reader can find the proofs for the error bound on Grise estimates. Appendix C contains the proofs relatives to the computational complexity of Grise and its efficient implementation. In Appendix D, the reader can find the proofs in connection with the NPC constant and relative to structure and parameter estimation with Suprise. Finally, Section E contains the proofs for the applications of Suprise for structure and parameter estimation of iconic special cases.

Appendix A About well-posedness and local learnability conditions

A.1 About condition (C1)

To illustrate further why Condition (C1) is required, we look at a case where the local constraint set is trivial, i.e. $\mathcal{Y}_{i}=\mathbbm{R}^{\mathbf{K}_{i}}$ and we consider a model that violates Condition (C1). This implies that there exists a sequence $\underline{x}_{n}\in\mathcal{X}_{i}$ such that $\underline{x}_{n}^{\top}\widetilde{I}(\theta^{*})\underline{x}_{n}/\|\underline{x}_{\mathcal{T}_{i}}\|^{2}<\rho_{n}$ with $\rho_{n}\rightarrow 0$. In the limit, we can find a vector $\underline{x}$ such that $\underline{x}^{\top}\widetilde{I}(\theta^{*})\underline{x}=0$ and $\|\underline{x}_{\mathcal{T}_{i}}\|=1$. In other words, it implies that for this model there exists a vertex $i$ and a perturbation vector $\underline{x}\in\mathcal{X}_{i}$ such that $\underline{x}_{\mathcal{T}_{i}}\neq 0$ and for which ${\mathbbm{E}}\left[\left(\sum_{k\in\mathcal{K}_{i}}x_{k}g_{ik}(\underline{\sigma}_{k})\right)^{2}\right]=0$. Since the probability distribution in Eq. (1) is positive, it further implies that for all configurations $\underline{\sigma}$ we have the functional equality $\sum_{k\in\mathcal{K}_{i}}x_{k}g_{ik}(\underline{\sigma}_{k})=0$. This enables us to locally reparameterize the distribution:

[TABLE]

where $c(\underline{x},\underline{\sigma}_{\setminus i})=\sum_{k\in\mathcal{K}_{i}}x_{k}\phi_{ik}(\underline{\sigma}_{k\setminus i})$ is a sum of locally centered functions that does not involve the variable $\sigma_{i}$. At this point, we should distinguish between the two cases when the basis functions $f_{k}$ are centered or not. When the basis functions are centered, i.e. $f_{k}=g_{k}$, the residual in Eq. (27) is identically zero, $c(\underline{x},\underline{\sigma}_{\setminus i})=0$. Therefore, the probability distribution of the model in (1) can be reparameterized entirely with $\theta^{*}_{k}\rightarrow\theta^{*}_{k}-x_{k}$ for $k\in\mathcal{K}_{i}$. It implies, as $\underline{x}_{\mathcal{T}_{i}}\neq 0$, that there exists two parameterization of the same models with different target parameters and the model selection problem as stated in Definition 1 is ill-posed. In the case when the basis functions are not centered, i.e. $f_{k}\neq g_{k}$, it may not be possible to reparameterized the whole distribution of the model. However, the conditional probability distribution $\mathbbm{P}(\sigma_{i}\mid\underline{\sigma}_{\setminus i})$ can be reparameterized as it is proportional to $\exp\left(\sum_{k\in\mathcal{K}_{i}}\theta^{*}_{k}f_{k}(\underline{\sigma}_{k})\right)$ and $\exp\left(\sum_{k\in\mathcal{K}_{i}}\left(\theta^{*}_{k}-x_{k}\right)f_{k}(\underline{\sigma}_{k})\right)$ thanks to Eq. (27). Thereby, even if the model is uniquely parameterized with $\underline{\theta}_{\mathcal{T}_{i}}$, local methods based on independent neighborhood reconstructions using conditional distributions will fail at selecting a unique model as shown in the following example.

Example 1**.**

Consider two family of models over two binary variables $\sigma,s\in\{-1,1\}$, parameterized by $\underline{\theta}$ and $\underline{\eta}$,

[TABLE]

and

[TABLE]

Both models are equivalent through the invertible mapping $\eta_{1}=\theta_{1}+\theta_{2}$, $\eta_{2}=\theta_{3}-\theta_{1}$ and $\eta_{3}=\theta_{3}-\theta_{2}$. However the model in Eq. (29) that has centered basis functions satisfies (C1) from Condition 1, while the model in Eq. (28) that has non-centered basis functions does not. This implies that the parameters $\underline{\theta}$ cannot be recovered by looking independently at conditional distributions as they are degenerate in this basis,

[TABLE]

Indeed the change of parameters $\theta_{1}\rightarrow\theta_{1}+\epsilon$, $\theta_{2}\rightarrow\theta_{2}-\epsilon$ and $\theta_{3}\rightarrow\theta_{3}+\epsilon$ leaves the conditional distribution (30) unchanged while the change of parameters $\theta_{1}\rightarrow\theta_{1}-\epsilon$, $\theta_{2}\rightarrow\theta_{2}+\epsilon$ and $\theta_{3}\rightarrow\theta_{3}-\epsilon$ leaves the conditional distribution (31) unaffected. Note that there does not exist a change of parameters that leaves both (30) and (31) unchanged. This is in agreement with the fact that the model in Eq. (28) is uniquely parameterized and can be in principle recovered by looking jointly at both conditional distributions.

For the specific models that we considered in Section E, the basis functions are always centered, which implies that failure to satisfy (C1) means that the model selection problem is ill-posed.

A.2 About condition (C2)

The bound on the interaction strength in (C2) translates directly into a uniform bound on the conditional probabilities of the models as shown in the following lemma.

Lemma 1** (Lower-Bounded Conditional Probabilities).**

Consider a graphical model with bounded maximum interaction strength of $\gamma=\max_{i\in\mathcal{V}}|\max_{\underline{\sigma}}\sum_{k\in\mathcal{K}_{i}}\theta^{*}_{k}g_{ik}(\underline{\sigma}_{k})|$. Then for any two disjoint subsets of vertices $A,B\subseteq\mathcal{V}$ the conditional probability of $\underline{\sigma}_{A}$ given $\underline{\sigma}_{B}$ is bounded from below,

[TABLE]

where $|\mathcal{A}_{i}|$ is the alphabet size of $\sigma_{i}$.

Proof of Lemma 1: Lower-bound on conditional probabilities.

We start by bounding the conditional probability of one variable $\sigma_{i}$ given the rest $\underline{\sigma}_{\setminus{i}}$. This is given by the following expression,

[TABLE]

as the centering functions $\phi_{ik}(\sigma_{k\setminus i})$ are independent of $\sigma_{i}$. The last expression can be simply bounded away from zero,

[TABLE]

using $\gamma\geq|\sum_{k\in\mathcal{K}_{i}}\theta^{*}_{k}g_{ik}(\underline{\sigma}_{k})|$.

Now we consider the conditional probability of one variable $\sigma_{i}$ given a subset of variable $\underline{\sigma}_{B}$ where $B\subseteq\mathcal{V}$ and $i\notin B$. Denote the complementary set of ${i}$ and $B$ by $S=\mathcal{V}\setminus{\left(\{i\}\cup B\right)}$. Then using the chain rule and the inequality from Eq. (36) we find the following lower-bound,

[TABLE]

Finally we consider the conditional probability of a set of variable $\underline{\sigma}_{A}$ given another set $\underline{\sigma}_{B}$ where $A\cap B=\emptyset$. Denote the vertices in $A$ by $\left\{1,2,\ldots,|A|\right\}$. Then using the chain rule and the inequality from Eq. (36), we obtain the desired result,

[TABLE]

∎

Appendix B Proofs of Grise estimation error bound

Proposition 3** (Gradient Concentration for Grise).**

For some node $u\in V$, let $n>\frac{2e^{2\gamma}}{\epsilon^{2}_{1}}\log(\frac{2\mathbf{K}_{u}}{\delta_{1}})$, then with probability at least $1-\delta_{1}$ the components of the gradient of the GISO are bounded from above as

[TABLE]

Define the residual of the first order Taylor expansion as

[TABLE]

Proposition 4** (Restricted Strong Convexity for Grise).**

For some node $u\in\mathcal{V}$, let $n>\frac{2}{\epsilon_{2}^{2}}\log\left(\frac{2\mathbf{K}_{u}^{2}}{\delta_{2}}\right)$ and assume that Condition 1 holds for some norm $\|\cdot\|$. Then, with probability at least $1-\delta_{2}$ the error of the first order Taylor expansion of the GISO satisfies

[TABLE]

for all $\Delta\in\mathcal{X}_{u}\subseteq\mathbbm{R}^{\mathbf{K}_{u}}$.

We first prove Theorem 1 before proving the propositions.

Proof of Theorem 1: Error Bound on Grise.

For some node $u\in\mathcal{V}$, let $n\geq\frac{2^{14}\widehat{\gamma}^{2}(1+\widehat{\gamma})^{2}e^{4\gamma}}{\alpha^{4}\rho_{u}^{2}}\log(\frac{4\mathbf{K}_{u}^{2}}{\delta})$. As the estimate $\widehat{\underline{\theta}}_{u}$ is an $\epsilon$-optimal point of the GISO and $\underline{\theta}^{*}_{u}$ lies in the constraint set from Eq. (6), we find that for $\Delta=\widehat{\underline{\theta}}_{u}-\underline{\theta}^{*}_{u}$

[TABLE]

Using the union bound on Proposition 3 and Proposition 4 with $\delta_{1}=\delta_{2}=\frac{\delta}{2}$ and

[TABLE]

we can express the inequality as

[TABLE]

Since by assumptions $\|\underline{\theta}^{*}_{u}\|_{1}\leq\gamma$ and $\|\widehat{\underline{\theta}}_{u}\|_{1}\leq\widehat{\gamma}$ for $\gamma\leq\widehat{\gamma}$ as the estimate is an $\epsilon$-optimal point of the $\ell_{1}$-constrained GISO, the error $\|\Delta\|_{1}$ is bounded by $2\widehat{\gamma}$. By choosing

and after some algebra, we obtain that

[TABLE]

∎

B.1 Gradient concentration

The components of the gradient of the GISO is given by

[TABLE]

Each term in the summation above is distributed as the random variable

[TABLE]

Lemma 2**.**

For any $u\in\mathcal{V}$ and $k\in\mathcal{K}_{u}$, we have

[TABLE]

Proof.

Simple computation. ∎

Proof of Proposition 3: Gradient Concentration for Grise.

The random variable $X_{uk}$ is bounded as

[TABLE]

Using Lemma 2 and the Hoeffding inequality, we get

[TABLE]

The proof follows by using (56) and the union bound over all $k\in\mathcal{K}_{u}$. ∎

B.2 Restricted strong convexity

We make use of the following deterministic functional inequality derived in [18].

Lemma 3**.**

The following inequality holds for all $z\in\mathbbm{R}$.

[TABLE]

Proof of Lemma 3.

Note that the inequality is true for $z=0$ and the first derivative of the difference is positive for $z>0$ and negative for $z<0$. ∎

Let $H_{k_{1}k_{2}}$ denote the correlation between $g_{k_{1}}$ and $g_{k_{2}}$ defined as

[TABLE]

and let $H=[H_{k_{1}k_{2}}]\in\mathbbm{R}^{|\mathcal{K}_{u}|\times|\mathcal{K}_{u}|}$ be the corresponding matrix. We define $\hat{H}$ similarly based on the empirical estimates of the correlation $\hat{H}_{k_{1}k_{2}}=\frac{1}{n}\sum_{t=1}^{n}g_{uk_{1}}(\underline{\sigma}_{k_{1}}^{(t)})g_{uk_{2}}(\underline{\sigma}_{k_{2}}^{(t)})$. The following lemma bounds the deviation between the above two quantities.

Lemma 4**.**

Choose some node $u\in\mathcal{V}$. With probability at least $1-2\mathbf{K}_{u}^{2}\exp\left(-\frac{n\epsilon_{2}^{2}}{2}\right)$, we have

[TABLE]

for all $k_{1},k_{2}\in\mathcal{K}_{u}$.

Proof of Lemma 4.

Fix $k_{1},k_{2}\in\mathcal{K}_{u}$. Then the random variable defined as $Y_{k_{1}k_{2}}=g_{uk_{1}}(\underline{\sigma}_{k_{1}})g_{uk_{2}}(\underline{\sigma}_{k_{2}})$ satisfies $|Y_{k_{1}k_{2}}|\leq 1$. Using the Hoeffding inequality we get

[TABLE]

The proof follows by using the union bound over $k_{1},k_{2}\in\mathcal{K}_{u}$. ∎

Lemma 5**.**

The residual of the first order Taylor expansion of the GISO satisfies

[TABLE]

Proof of Lemma 5.

Using Lemma 3 we have

[TABLE]

The proof follows by observing that $|\sum_{k\in\mathcal{K}_{u}}\Delta_{k}g_{uk}(\underline{\sigma}_{k}^{(t)})|\leq\|\Delta\|_{1}$. ∎

We are now in a position to complete the proof of Proposition 4.

Proof of Proposition 4: Restricted Strong Convexity for Grise.

Using Lemma 5 we have

[TABLE]

where $(a)$ follows from Lemma 4 and $(b)$ follows from Condition 1 as

[TABLE]

∎

Appendix C Efficient implementation of Grise and its computational complexity

The iterative Algorithm 2 takes as input a number of steps $T$ and output an $\epsilon$-optimal solution of Grise without constraints in Eq. (74). This algorithm is an application of the Entropic Descent Algorithm introduced by [1] to reformulation of Eq. (13) as a minimization over the probability simplex. Note that there exist other efficient iterative methods for minimizing the GISO, such as the mirror gradient descent of [2]. The following proposition provides guarantees on the computational complexity of unconstrained Grise.

Proposition 5** (Computational Complexity for Unconstrained Grise).**

Let $1\geq\epsilon>0$ be the optimality gap and $T\geq 6\epsilon^{-2}\ln{(2\mathbf{K}_{u}+1)}$ be the maximum number of iterations. Then Algorithm 2 is guaranteed to produce an $\epsilon$-optimal solution of Grise without a constraint set $\mathcal{Y}_{u}$ with a number of operation less than $C\frac{c_{g}n\mathbf{K}_{u}}{\epsilon^{2}}\ln(1+\mathbf{K}_{u})$, where $c_{g}$ is an upper bound on the computational complexity of evaluating any $g_{ik}(\underline{\sigma}_{k})$ for $k\in\mathcal{K}_{i}$ and $C$ is a universal constant that is independent of all parameters of the problem.

Proof of Proposition 5: Computational complexity of unconstrained Grise.

We start by showing that the minimization of Grise in Eq. (13) in the unconstrained case where $\mathcal{Y}_{u}=\mathbb{R}^{\mathbf{K}_{u}}$ is equivalent to the following lifted minimization on the logarithm of Grise,

[TABLE]

We first show that for all $\underline{\theta}_{u}\in\mathbb{R}^{\mathbf{K}_{u}}$ such that $\|\underline{\theta}_{u}\|_{1}\leq\widehat{\gamma}$, there exists $\underline{x}^{+},\underline{x}^{-},y$ satisfying constraints (71), (72), (73). This is easily done by choosing $x^{+}_{k}=\max(\theta_{k}/\widehat{\gamma},0)$ , $x^{-}_{k}=\max(-\theta_{k}/\widehat{\gamma},0)$ and $y=1-\|\underline{\theta}_{u}\|_{1}/\widehat{\gamma}$. Second, we trivially see that for all $\underline{\theta}_{u},\underline{x}^{+},\underline{x}^{-},y$ satisfying constraints (71), (72), (73), it implies that $\underline{\theta}_{u}$ also satisfies $\|\underline{\theta}_{u}\|_{1}\leq\widehat{\gamma}$. Therefore, any $\underline{\theta}^{\rm min}_{u}$ that is an argmin of Eq. (70) is also an argmin of Eq. (13) without constraint set $\mathcal{Y}_{u}$. Moreover, if we find $\underline{\theta}^{\rm\epsilon}_{u}$ such that $\log\mathcal{S}_{n}(\underline{\theta}^{\rm\epsilon}_{u})-\log\mathcal{S}_{n}(\underline{\theta}^{\min}_{u})\leq\epsilon/\sqrt{3}$, we obtain an $\epsilon$-minimizer of Eq. (13) without constraint set $\mathcal{Y}_{u}$. Indeed, since $\epsilon/\sqrt{3}\leq\log(1+\epsilon)$ for $1\geq\epsilon>0$, we have that $\mathcal{S}_{n}(\underline{\theta}^{\rm\epsilon}_{u})-\mathcal{S}_{n}(\underline{\theta}^{\min}_{u})\leq\mathcal{S}_{n}(\underline{\theta}^{\min\epsilon}_{u})\epsilon\leq\epsilon$ as $\mathcal{S}_{n}(\underline{\theta}^{\min}_{u})\leq\mathcal{S}_{n}(0)=1$. The remainder of the proof is a straightforward application of the analysis of the Entropic Descent Algorithm in [1, Th. 5.1] to the above minimization where $\underline{\theta}_{u}$ has been replaced by $\underline{x}^{+},\underline{x}^{-},y$ using Eq. (71). In this analysis we use the fact that the logarithm of Grise remains a convex function as it is a sum of exponential functions and also that the gradient of our objective function is bounded uniformly by $\|\nabla\log\mathcal{S}_{n}(\underline{\theta}_{u})\|_{\infty}=\|\nabla\mathcal{S}_{n}(\underline{\theta}_{u})/\mathcal{S}_{n}(\underline{\theta}_{u})\|_{\infty}\leq 1$ as $|g_{k}(\underline{\sigma}_{k})|\leq 1$. Note that the computational complexity of the gradient evaluation is proportional to $n\mathbf{K}_{u}c_{g}$. This is because for each sample, one has to first compute an exponential containing $\mathbf{K}_{u}$ terms $g_{k}(\underline{\sigma}_{k})$ with an evaluation cost of $c_{g}$ and then multiply the exponential by the factor $-g_{k}(\underline{\sigma}_{k})$ corresponding to each of the $\mathbf{K}_{u}$ components of the gradient. ∎

When the constraint set $\mathcal{Y}_{u}$ is parametrically complete, an $\epsilon$-optimal solution to (13) can be found by first solving the unconstrained version of Grise and then performing an equi-cost projection onto $\mathcal{Y}_{u}$. We define an $\epsilon$-optimal solution to the unconstrained Grise problem as

[TABLE]

Lemma 6**.**

Let $\widehat{\underline{\theta}}_{u}^{unc}$ be an $\epsilon$-optimal solution of the unconstrained Grise problem in Eq. (74). Then an equi-cost projection of $\widehat{\underline{\theta}}_{u}^{unc}$ is an $\epsilon$-optimal solution of the constrained Grise problem,

[TABLE]

Proof of Lemma 6.

Since unconstrained Grise is a relaxation of Grise with the constraints $\underline{\theta}_{u}\in\mathcal{Y}_{u}$, we must have

[TABLE]

Since, $\mathcal{Y}_{u}$ is parametrically complete, by definition,

[TABLE]

The estimates $\mathcal{P}_{\mathcal{Y}_{u}}(\widehat{}\underline{\theta}_{u}^{\text{unc}})$ are feasible for the constrained Grise problem, completing the proof. ∎

Lemma 6 implies that the computational complexity of Grise for parametrically complete cases is the sum of the computational complexity of the unconstrained Grise and the projection step.

Algorithm 3 is an implementation of Grise for parametrically complete cases. Its computational complexity is obtained easily by combining Lemma 6 and Proposition 5.

Theorem 3** (Computational Complexity for Grise with P.C. Constraints).**

Let $\mathcal{Y}_{u}$ be a parametrically complete set and let $1\geq\epsilon>0$ be given. Then Algorithm 3 computes an $\epsilon$-optimal solution to Grise with a number of operations bounded by $C\frac{c_{g}n\mathbf{K}_{u}}{\epsilon^{2}}\ln(1+\mathbf{K}_{u})+\mathcal{C}(\mathcal{P}_{\mathcal{Y}_{u}}(\widehat{}\underline{\theta}_{u}^{\text{unc}}))$, where $c_{g}$ is an upper bound on the computational complexity of evaluating any $g_{ik}(\underline{\sigma}_{k})$ for $k\in\mathcal{K}_{i}$ and where $\mathcal{C}(\mathcal{P}_{\mathcal{Y}_{u}}(\widehat{}\underline{\theta}_{u}^{\text{unc}}))$ denotes the computational complexity of the projection step.

Appendix D Proofs & algorithms for structure and parameter estimation

D.1 Dimension independence and easier computation of NPC constants

We recall the definition of the NPC constant,

[TABLE]

In order to give an intuition for the intricate formula in Eq. (78), let us define for maximal cliques $c\in\mathcal{M}_{\rm{cli}}\left(\mathcal{G}\right)$, their clique parameterization matrix $G^{c}$, with indices being maximal factors $k,k^{\prime}\in\left[c\right]_{\rm{sp}}$. The clique parameterization matrix is obtained by summing over variable $\underline{\sigma}_{c}$ the globally centered basis functions,

[TABLE]

Note that the clique parametertization matrices are positive semi-definite matrices by construction. Bounding the expectation over $\sigma_{i}$ using Lemma 1, we see that the NPC constant is linked to the smallest eigenvalue of the clique parameterization matrix,

[TABLE]

Clique parametertization matrices have a typical size of $\mathcal{O}(q^{L}\times q^{L})$ since variables in a clique can take up to $\mathcal{O}(q^{L})$ different configurations. Therefore, Eq. (80) emphasizes that that the NPC constant does not depend on the dimension of the model $p$ but rather on local properties of the parameterization of the family.

D.2 Proofs for local learnability condition from nonsingular parametrization of cliques

Proof of Proposition 1: LLC in $\ell_{\infty,2}$-norm.

For a given vertex $i\in\mathcal{V}$, let $\underline{x}\in\mathcal{X}_{i}\subseteq\mathbbm{R}^{\mathbf{K}_{i}}$ be a vector in the perturbation set. First, suppose that $\{i\}$ is not a maximal clique and choose any maximal clique $c\in\mathcal{M}_{\rm{cli}}\left(\mathcal{G}\right)$ that contains the vertex $i$ and let the set $S=c\setminus\{i\}$ be the set of nodes in the clique without $i$. The expression characterizing the LLC can be evaluated conditioning the expectation over $S$. Denoting the marginal and the conditional probability distribution used to compute expectation by subscripts, we find,

[TABLE]

where in the last line, we bounded the probability of $\mathbbm{P}(\underline{\sigma}_{S}\mid\underline{\sigma}_{\setminus S})$ using Lemma 1. We want to rewrite the sum over $\underline{\sigma}_{S}$ in Eq. (82) using globally centered functions $h_{k}$ for factors $k\in\left[c\right]_{\rm{sp}}$ instead of locally centered functions $g_{ik}$. Using definitions of locally centered functions in Eq. (3) and globally centered functions in Eq. (22), we see that $g_{ik}(\underline{\sigma}_{k})=h_{k}\left(\underline{\sigma}_{k}\right)+R_{ik}(\underline{\sigma}_{k})$, where

[TABLE]

The sum in Eq. (82) can be expanded into the four contributions,

[TABLE]

We start by evaluating the contribution from terms in Eq. (85). For $k\in\left[c\right]_{\rm{sp}}$ and $l\in\mathcal{K}_{i}\setminus{\left[c\right]_{\rm{sp}}}$, there exists at least one node $u\in c$ such that $u\neq i$ and $u\notin\partial l$. Summing over the variable $\sigma_{u}$ cancels the expression,

[TABLE]

as $h_{k}(\underline{\sigma}_{c})$ is globally centered and $g_{il}(\underline{\sigma}_{l})$ does not depend on $\sigma_{u}$.

The contribution from Eq. (86) is also null. To see this we expand the sums using the formula for the reminder in Eq. (83),

[TABLE]

where the sum over $\underline{\sigma}_{S}$ vanishes as $h_{k}(\underline{\sigma}_{c})$ depends on $\underline{\sigma}_{r}\neq\sigma_{i}$ while $\sum_{\underline{\sigma}_{r}}f_{l}(\underline{\sigma}_{c})$ does not. As the contribution from Eq.(87) is non-negative, we can lower-bound Eq. (82) by the following expression,

[TABLE]

where in the last line we have recognized the definition of the NPC constant from Eq. (23). Since (91) holds for any $c\in\mathcal{M}_{\rm{cli}}\left(\mathcal{G}\right)$ that contains the vertex $i$, the Local Learnability Condition is satisfied for a weighted $\ell_{\infty,2}$-norm with LLC constant equal to $\rho^{\rm{NPC}}$,

[TABLE]

where the weighted $\ell_{\infty,2}$-norm is defined as follows,

[TABLE]

As the weighted $\ell_{\infty,2}$-norm in Eq. (93) is lower-bounded by the $\ell_{\infty,2}$-norm,

[TABLE]

we have that the LLC is also satisfied for the $\ell_{\infty,2}$-norm with LLC constant equal to $\rho^{\rm{NPC}}_{i}\left(\frac{\exp(-2\gamma)}{q}\right)^{L-1}$.

When $\{i\}$ is a maximal clique, then $\mathcal{K}_{i}=\left[\{i\}\right]_{\rm{sp}}$ and it straightforward to see that

[TABLE]

∎

Lemma 7**.**

Let $\sigma\in\mathcal{A}$, be a discrete random variable with probability distribution $p(\sigma)$. Consider $x_{\sigma}\in\mathbb{R}$, a function defined over $\sigma$ that is centered, i.e. $\sum_{\sigma\in\mathcal{A}}x_{\sigma}=0$. The variance of the function $x_{\sigma}$ is lower-bounded by,

[TABLE]

where $p_{\text{min}}=\min_{\sigma\in\mathcal{A}}p(\sigma)$.

Proof.

The proof goes as follows,

[TABLE]

where in the last line we used that $\sum_{\sigma\in\mathcal{A}}x_{\sigma}=0$ and $\left(\sum_{\sigma^{\prime}\in\mathcal{A}}p(\sigma^{\prime})x_{\sigma^{\prime}}\right)^{2}\geq 0$. ∎

Proof of Proposition 2: LLC in $\ell_{2}$-norm for pairwise models.

For a given vertex $i\in\mathcal{V}$, let $\underline{x}\in\mathcal{X}_{i}\subseteq\mathbbm{R}^{\mathbf{K}_{i}}$ be a vector in the perturbation set. When $\{i\}$ is a maximal clique, then $\mathcal{K}_{i}=\left[\{i\}\right]_{\rm{sp}}$ and we immediately see that

[TABLE]

Now suppose that $\{i\}$ is not a maximal clique, i.e. there exists $j\in\mathcal{V}$ such that $\{i,j\}\in\mathcal{M}_{\rm{cli}}\left(\mathcal{G}^{*}\right)$. The expectation that arises in the LLC is lower-bounded by its variance,

[TABLE]

Let $\{S_{r}\}_{r=1,\ldots,\chi}$ be a minimal coloring of the graph $\mathcal{G}^{*}$. For a given color $r$, define the set $C_{r}=S_{r}\setminus\{i\}$ and apply the law of total variance on the right-hand side of Eq. (104), conditioning on $\underline{\sigma}_{\setminus C_{r}}$,

[TABLE]

where the marginal and the conditional probability distribution used to compute expectation and variance respectively are indicated by subscripts. As the variance on the right-hand side of Eq. (105) is conditioned on $\underline{\sigma}_{\setminus C_{r}}$, only basis functions involving a pair $(\sigma_{i},\sigma_{j})$ with $j\in C_{r}$ are giving a non-zero contribution to the conditional variance,

[TABLE]

We can rewrite the locally centered functions with respect to globally centered functions using their definitions found in Eq. (3) and in Eq. (22),

[TABLE]

We see from Eq. (107) that the difference between locally and globally centered functions only depends on the variable $\sigma_{i}$. This means that we can interchange locally centered functions with globally centered functions in the right-hand side of Eq. (106) as the variance is conditioned on $\sigma_{i}$,

[TABLE]

Since $\{S_{r}\}_{r=1,\ldots,\chi}$ is a vertex coloring, by definition all nodes $j\in C_{r}$ having the same color are not sharing a factor node, i.e. $\forall j_{1},j_{2}\in C_{r}$, $\nexists k\in\mathcal{K}^{*}$ such that $j_{1},j_{2}\in\partial k$. This implies that variables $\sigma_{j}$ with $j\in C_{r}$ are independent conditioned on the remaining variables $\underline{\sigma}_{\setminus C_{r}}$ and the variance in Eq. (108) can be rewritten,

[TABLE]

The right-hand side of Eq. (109) is centered with respect to $\sigma_{j}$ and we can apply Lemma 7 and Lemma 1 to find a lower-bound that is only dependant on the random variable $\sigma_{i}$,

[TABLE]

Plugging back the results derived in Eq. (108), Eq. (109), and Eq. (110) into the initial inequality in Eq. (105), we find,

[TABLE]

where in Eq. (112) we used the definition of the NPC constant in Eq. (23) to bound the quadratic form involving $\underline{x}$.

Finally, we average the inequality described by Eq. (112) over the different colors and hence possible conditioning sets $C_{r}$ to conclude the proof,

[TABLE]

∎

D.3 Proofs of estimation guarantees for the Suprise algorithm

Proof of Theorem 2: Reconstruction and Estimation Guarantees for Suprise.

As the NPC constant is non-zero $\rho^{\rm{NPC}}_{u}>0$ for all nodes $u\in\mathcal{V}$, we apply Proposition 1 in conjunction with Theorem 1 to find that for each step $t\in\{0,\ldots,L-1\}$ and with probability at least $1-\delta/(pL)$, Grise around a node $u\in\mathcal{V}$ recovers the parameters in each maximal clique $c\in\mathcal{M}_{\rm{cli}}\left(\mathcal{G}\left[(\mathcal{V},\mathcal{K}^{t})\right]\right)$ that contains $u$ with precision $\sum_{k\in\left[c\right]_{\rm{sp}}}(\underline{\theta}_{k}-\widehat{}\underline{\theta}_{k})^{2}\leq(\alpha/2)^{2}$. Therefore, at each step $t\in\{0,\ldots,L-1\}$ and with probability at least $1-\delta/L$, the factor removal procedure is guaranteed to remove all factors of size $L-t$ that are not present in the graph if all factors of size bigger than $L-t$ were correctly removed in the previous steps. Since there are at most $L$ removal steps, it implies that the overall procedure discovers all maximal cliques with probability at least $1-\delta$.

∎

Appendix E Application to special cases

In this section, we show how to apply Theorem 2 in order to derive the sample and computational complexity of reconstructing graphical models for some common basis functions.

E.1 Binary models on the monomial basis

In this subsection, we consider general models on binary alphabet $\mathcal{A}_{i}=\{-1,1\}$. Let the factors be all nonempty subsets of $\{1,\dots,p\}=\mathcal{V}$ of size at most $L$,

[TABLE]

The set $\mathcal{K}$ contains all potential subsets of variable of size at most $L$. The parameterization uses the monomial basis given by $f_{k}(\underline{\sigma}_{k})=\prod_{j\in k}\sigma_{j}$ with $k\in\mathcal{K}$. Note that the monomial basis functions are already globally centered $f_{k}\equiv g_{k}\equiv h_{k}$. The probability distribution for this model is expressed as

[TABLE]

When $L\leq 2$, the model in Eq. (116) is pairwise and it is referred as the Ising Model.

For each maximal clique there exists exactly one maximal factors in its span. Therefore, the NPC constant as defined in Eq. (23) is $\rho_{\rm{NPC}}=1$ since for any clique $c$ we have,

[TABLE]

and the minimum is achieved for cliques of size one. As every node is involved in at most $\mathbf{K}\leq p^{L-1}$ factor functions, the structure of binary models can be recovered as a corollary of Theorem 2.

Corollary 1** (Structure recovery for binary graphical models).**

Let $\underline{\sigma}^{(1)},\ldots,\underline{\sigma}^{(n)}$ be i.i.d. samples drawn according to $\mu(\underline{\sigma})$ in (116) and let $\alpha\leq\min_{k\in\mathbbm{S}(\mathcal{G}^{*})}|\theta^{*}_{k}|$ be the intensity of the smallest non-zero parameter. If

[TABLE]

then the structure of the binary graphical model is perfectly recovered using Algorithm 1, i.e. $\widehat{\mathbbm{S}}=\mathbbm{S}(\underline{\theta}^{*})$, with probability $1-\delta$. Moreover the total computational complexity scales as $\widetilde{\mathcal{O}}(p^{L})$, for fixed $L$, $\alpha$, $\gamma$, $\widehat{\gamma}$ and $\delta$.

For pairwise Ising models that are $\chi$ colorable, we have also guarantees on the $\ell_{2}$-norm reconstruction by Suprise of pairwise parameters.

Corollary 2** ($\ell_{2}$-parameter estimation for Ising models).**

Let $\underline{\sigma}^{(1)},\ldots,\underline{\sigma}^{(n)}$ be i.i.d. samples drawn according to $\mu(\underline{\sigma})$ in (116) for $L=2$ and let $\alpha>0$ be the prescribed estimation accuracy. If

[TABLE]

then, with probability at least $1-\delta$, the parameters are estimated by Algorithm 1 with the error

[TABLE]

The computational complexity of obtaining these estimates is $\widetilde{\mathcal{O}}(p^{2})$ for fixed $\chi$, $\alpha$, $\gamma$, $\widehat{\gamma}$ and $\delta$.

As graphs with bounded degree $d$ have a chromatic number at most $d+1\geq\chi$, Corollary 2 recovers the $\ell_{2}$-guarantees for sparse graphs recovery of [18] albeit with slightly worse dependence with respect to $\gamma$ and $\alpha$. The worse $\gamma$ dependence is an artifact of the general analysis presented in this paper. For models over binary variables one can improve the $e^{8\gamma}$ dependence to $e^{6\gamma}$ using Berstein’s inequality in Proposition 3 instead of Hoeffding’s inequality. However, the worse $\alpha$ dependence seems to be more fundamental. It is caused by the replacement of the $\ell_{1}$-penalty used in [18] by an $\ell_{1}$-constraint.

For graphs with unbounded vertex degree but low chromatic number, such as star graphs or bipartite graphs, Corollary 2 shows that the parameters of the corresponding Ising model can be fully recovered with a bounded $\ell_{2}$-error using a number of samples that is only logarithmic in the model size $p$.

E.2 L-wise models with arbitrary alphabets on the indicator basis

In this subsection, we consider $L$-wise graphical models over variables taking values in arbitrary alphabet $\mathcal{A}_{i}$ of size $q_{i}$, parametrized with indicator-type functions. The building block of the set of basis functions is the centered univariate indicator function defined as

[TABLE]

where $s_{i},\sigma_{i}\in\mathcal{A}_{i}$ are prescribed letters of the alphabet. The univariate indicator functions in Eq. (121) are centered Kronecker delta functions and possess similar properties such as symmetry $\Phi_{s_{i},\sigma_{i}}=\Phi_{\sigma_{i},s_{i}}$ and contraction under a summation,

[TABLE]

The set of factors $\mathcal{K}$ are pairs associating elements of $R=\{r\in P(\mathcal{V})\mid|r|\leq L\}$ which are subsets of variable of size at most $L$ with an alphabet configuration in $\mathcal{A}_{r}=\bigotimes_{i\in r}\mathcal{A}_{i}$,

[TABLE]

In what follows, we slightly abuse the notation of factors and parameters by shortening $(r,\underline{s}_{r})\equiv\underline{s}_{r}$. With these notations, the indicator basis functions are constructed as $f_{\underline{s}_{r}}(\underline{\sigma}_{r})=\prod_{i\in r}\Phi_{s_{i},\sigma_{i}}$. Note that the indicator basis functions are globally centered i.e. $f_{\underline{s}_{r}}\equiv g_{\underline{s}_{r}}\equiv h_{\underline{s}_{r}}$. The probability distribution of an $L$-wise graphical model with arbitrary alphabet is defined as follows,

[TABLE]

The family of distribution in Eq. (124) is not uniquely parameterized by the parameters $\underline{\theta}^{*}$. To see this, we introduce the linear application $\mathcal{P}_{r}$ acting on arrays $\theta_{\underline{s}_{r}}$ as follows,

[TABLE]

Using the contraction property from Eq. (122), it is easy to see that $\mathcal{P}_{r}$ is a projector, i.e $\mathcal{P}^{2}_{r}=\mathcal{P}_{r}$. It is also straightforward to verify that $\mathcal{P}_{r}\underline{\theta}$ is always a globally centered array and if $\underline{\theta}$ is already globally centered then $\mathcal{P}_{r}\underline{\theta}=\underline{\theta}$. Therefore, the applications $\mathcal{P}_{r}$ are projectors on the space of array $\theta_{\underline{s}_{r}}$ which are globally centered, i.e. $\sum_{s_{i}}\theta_{\underline{s}_{r}}=0$ for all $i\in r$. We lift the parametrization degeneracy in Eq. (124) by imposing that parameters $\underline{\theta}^{*}$ are in the range of the projector $\mathcal{P}_{r}$. We thus require that the parameters satisfy the following linear constraints at each vertex $u\in\mathcal{V}$,

[TABLE]

The constraint set in Eq. (126) is parametrically complete according to Definition 2 as we explicitly exhibited the equi-cost projection $\{\mathcal{P}_{r}\}_{r\in R}$ onto it. The computational complexity of this projection is no more than $\mathcal{O}(p^{L-1}q^{L})$.

As the constraint set in Eq. (126) forms a linear subspace, the perturbation set is simply $\mathcal{X}_{u}=\mathcal{Y}_{u}\cap B_{1}(2\widehat{\gamma})$, the intersection of the constraint set with the $\ell_{1}$-ball of radius $2\widehat{\gamma}$. Maximal cliques are subset of vertices and hence are also elements of $R$. Therefore, the NPC constant as defined in Eq. (23) is bounded by $\rho_{\rm{NPC}}\geq\exp(-2\gamma)/q$ since for each clique we have,

[TABLE]

as $\underline{x}\in\mathcal{X}_{u}$ is globally centered and thus is in the range of the projector $\mathcal{P}_{c}$. Every node is involved in at most $\mathbf{K}\leq p^{L-1}q^{L}$ factor functions and the structure of L-wise models with arbitrary alphabets can be recovered as a corollary of Theorem 2.

Corollary 3** (Structure recovery for L-wise graphical models).**

Let $\underline{\sigma}^{(1)},\ldots,\underline{\sigma}^{(n)}$ be i.i.d. samples drawn according to $\mu(\underline{\sigma})$ in (124) and let $\alpha\leq\min_{c\in\mathbbm{S}(\mathcal{G}^{*})}\sqrt{\sum_{\underline{s}_{c}\in\mathcal{A}_{c}}{\theta^{*}_{\underline{s}_{c}}}^{2}}$ be the intensity of the smallest non-zero parameter. If

[TABLE]

then the structure of the L-wise graphical model with arbitrary alphabets is perfectly recovered using Algorithm 1, i.e. $\widehat{\mathbbm{S}}=\mathbbm{S}(\underline{\theta}^{*})$, with probability $1-\delta$. Moreover the total computational complexity scales as $\widetilde{\mathcal{O}}(p^{L})$, for fixed $L$, $q$, $\alpha$, $\gamma$, $\widehat{\gamma}$ and $\delta$.

For pairwise models with arbitrary alphabet that are $\chi$ colorable, we have also guarantees on the $\ell_{2}$-norm reconstruction by Suprise of pairwise parameters.

Corollary 4** ($\ell_{2}$-parameter estimation for pairwise models).**

Let $\underline{\sigma}^{(1)},\ldots,\underline{\sigma}^{(n)}$ be i.i.d. samples drawn according to $\mu(\underline{\sigma})$ in (124) for $L=2$ and let $\alpha>0$ be the prescribed estimation accuracy. If

[TABLE]

then, with probability at least $1-\delta$, the parameters are estimated by Algorithm 1 with the error

[TABLE]

The computational complexity of obtaining these estimates is $\widetilde{\mathcal{O}}(p^{2})$ for fixed $\chi$, $q$, $\alpha$, $\gamma$, $\widehat{\gamma}$ and $\delta$.