Generalized Score Matching for Non-Negative Data
Shiqing Yu Mathias Drton Ali Shojaie
































































 

 
 


















-00pt 






 
 
 
	 
	 
	 
	
	 
	 
	
 
	 
	 
	
	 
	 


 
 

 
 
 
 
 
 
 






theoremTheorem


fthm
  
  
  
lemmaLemma

  
 flem
  
  
  
  corollaryCorollary


  
    remarkRemark


paragraphsection

definitionDefinition



§ 


§ 




 §.§ 


 §.§ 




  §.§.§ 


  §.§.§ 











1

2020191-LastPage5/18; Revised
4/194/1918-278Shiqing Yu, Mathias Drton, Ali Shojaie
Generalized Score MatchingYu, Drton and Shojaie







Aapo Hyvarinen

Generalized Score Matching for Non-Negative Data
    Shiqing Yu sqyu@uw.edu 

       Department of Statistics

       University of Washington, Seattle, WA, U.S.A. 
       
              Mathias Drton md5@uw.edu 

       Department of Mathematical Sciences

       University of Copenhagen, Copenhagen, Denmark

       
       and

       Department of Statistics

       University of Washington, Seattle, WA, U.S.A. 
       
       Ali Shojaie ashojaie@uw.edu 

       Department of Biostatistics

       University of Washington, Seattle, WA, U.S.A.
       
       
     
    December 2017 (Updated Feb 2018 and Dec 2018)
=============================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






A common challenge in estimating parameters of probability density
functions is the intractability of the normalizing constant. While
in such cases maximum likelihood estimation may be implemented using
numerical integration, the approach becomes computationally
intensive.  The score matching method of <cit.>
avoids direct calculation of the normalizing constant and yields
closed-form estimates for exponential families of continuous
distributions over ℝ^m. <cit.> extended the
approach to distributions supported on the non-negative orthant, 
ℝ_+^m. In this paper, we give a generalized form of score
matching for non-negative data that improves estimation
efficiency. As an example, we consider a general class of pairwise
interaction models. Addressing an overlooked inexistence problem, we
generalize the regularized score matching method of
<cit.> and improve its theoretical guarantees for non-negative Gaussian graphical models.



exponential family, graphical model, positive data, score matching, sparsity






§ INTRODUCTION



Score matching was first developed in <cit.> for continuous distributions supported on all of ℝ^m.  Consider such a distribution P_0, with density p_0 and support equal to ℝ^m.


Let 𝒫 be a family of distributions with twice continuously differentiable densities. The score matching estimator of p_0 using 𝒫 as a model is the minimizer of the expected squared ℓ_2 distance between the gradients of
log p_0 and a log-density from 𝒫.  So we minimize the
loss
∫_ℝ^mp_0(x)∇log
p(x)-∇log
p_0(x)_2^2x 
with respect to densities p from 𝒫.  The loss depends on p_0, but integration by parts can be used to rewrite it in a form that can be
approximated by averaging over the sample without knowing p_0.  A key feature of score matching is that normalizing constants cancel in gradients of log-densities, allowing for simple treatment of models with intractable normalizing constants.  
For exponential families, the loss is quadratic in the canonical parameter, making optimization straightforward. 

If the considered distributions are supported on a proper subset of ℝ^m, then the integration by parts arguments underlying the score matching estimator may fail due to discontinuities at the boundary of the support.  For data supported on the non-negative orthant ℝ_+^m, <cit.> addresses this problem by modifying the loss to 
∫_ℝ^mp_0(x)∇log
p(x)∘x-∇log
p_0(x)∘x_2^2x, where ∘ denotes entrywise multiplication. In this loss, boundary effects are dampened by multiplying gradients elementwise with the identity functions x_j. 

In this paper, we propose generalized score matching methods that are based on element­wise multiplication with functions other than x_j. As we show, this can lead to drastically improved estimation accuracy, both theoretically and empirically.  To demonstrate these advantages, we consider a family of 
graphical models on ℝ_+^m, which does not have tractable normalizing constants and hence serves as a practical example.


Graphical models specify conditional independence relations for a
random vector X=(X_i)_i∈ V indexed by the nodes
of a graph <cit.>.
For undirected graphs, variables X_i and X_j are required to be conditionally independent given  (X_k)_k≠i,j if there is no edge between i and j. The smallest undirected graph with this property is the conditional independence graph of X.  Estimation of this graph and associated interaction parameters has been a topic of continued research as reviewed by <cit.>.

Largely due to their tractability, Gaussian graphical models (GGMs)
have gained great popularity.  The conditional independence graph of
a multivariate normal vector
X∼𝒩(μ,Σ)
is determined by the inverse covariance matrix
𝐊≡Σ^-1, also termed
concentration or precision matrix.  Specifically, X_i and X_j are conditionally independent given all other variables  if and only if the (i,j)-th and the (j,i)-th entries of 𝐊 are both zero. This simple relation underlies a rich literature including <cit.>, <cit.>, <cit.> and <cit.>, among others.

More recent work has provided tractable procedures also for
non-Gaussian graphical models. This includes Gaussian copula models
<cit.>, Ising models <cit.>, other
exponential family models <cit.>, as well as semi- or
non-parametric estimation techniques <cit.>. In this
paper, we apply our method to a class of pairwise interaction models
that generalizes non-negative Gaussian random variables, as recently
considered by <cit.> and <cit.>, as well as square root
graphical models proposed by <cit.> when the sufficient
statistic function is a pure power. However, our main ideas can also
be applied for other classes of exponential families whose
support is restricted to a rectangular set.

Our focus will be on pairwise interaction power models with
probability distributions having (Lebesgue) densities proportional to

    exp{-1/2ax^a^⊤𝐊x^a+η^⊤x^b-1_m/b}

on ℝ_+^m≡[0,∞)^m. Here a>0 and b≥ 0 are known
constants, and 𝐊∈ℝ^m× m and
η∈ℝ^m are unknown parameters of
interest. When b=0 we define (x^b-1)/b≡log x and
ℝ_+^m≡(0,∞)^m.  This class of models is motivated
by the form of important univariate distributions for non-negative
data, including gamma and truncated normal distributions.  It provides a
framework for pairwise interaction that is concrete yet rich enough to
capture key differences in how densities may behave at the boundary of
the non-negative orthant, ℝ_+^m.  Moreover, the conditional
independence graph of a random vector X with
distribution as in (<ref>) is determined just as in the
Gaussian case: X_i and X_j are conditionally independent given all
other variables if and only if κ_ij=κ_ji=0 in the
interaction matrix 𝐊.  Section <ref> gives further details on these models. We will 
develop estimators of
(η,𝐊) in (<ref>) and the associated conditional independence graph using the proposed generalized score matching.



A special case of (<ref>) are truncated Gaussian
graphical models, with a=b=1.  Let
μ∈ℝ^m, and let 𝐊 be a positive
definite matrix.  Then a non-negative random vector X
follows a truncated normal distribution for mean parameter
μ and inverse covariance parameter 𝐊, in
symbols X∼TN(μ,𝐊),
if it has density proportional to

    exp{-1/2(x-μ)^⊤𝐊(x-μ)}

on ℝ_+^m. We refer to
Σ=𝐊^-1 as the covariance parameter of
the distribution, and note that the η parameter in (<ref>) is 𝐊μ.  Another special case of (<ref>) is the exponential square root graphical models in <cit.>, where a=b=1/2.

<cit.> estimate truncated GGMs based on Hyvärinen's modification, with an ℓ_1 penalty on the entries of 𝐊 added to the loss. However, the paper overlooks the fact that the loss can be unbounded from below in the high-dimensional setting even with an ℓ_1 penalty, such that no minimizer may exist. Since the unpenalized loss is quadratic in the parameter to be estimated, we propose modifying it by adding small positive values to the diagonals of the positive semi-definite matrix that defines the quadratic part, in order to ensure that the loss is bounded and strongly convex and admits a unique minimizer. We apply this to the estimator for GGMs considered in <cit.>, which uses score-matching on ℝ^m, and to the generalized score matching estimator for pairwise interaction power models on ℝ_+^m proposed in this paper. In these cases, we show, both empirically and theoretically, that the consistency results still hold (or even improve) if the positive values added are smaller than a threshold that is readily computable.

The rest of the paper is organized as follows. Section <ref> introduces score matching and our proposed
generalized score matching. In Section
<ref>, we apply generalized score matching to
exponential families, with univariate truncated normal distributions
as an example. Regularized generalized score matching for
graphical models is formulated in Section <ref>. The estimators for pairwise interaction
power models are shown in Section <ref>, while
theoretical consistency results are presented in Section <ref>, where we treat the probabilistically most
tractable case of truncated GGMs. Simulation results and applications
to RNAseq data are given in Section <ref>. Proofs for theorems in Sections <ref>–<ref> are presented in
Appendices <ref> and <ref>. 
                                
Additional experimental results are presented in Appendix <ref>.







 §.§ Notation


Constant scalars, vectors, and functions are written in lower-case
(e.g., a, a), random scalars and vectors in
upper-case (e.g., X, X).  Regular font is used for
scalars (e.g. a, X), and boldface for vectors
(e.g. a, X).  Matrices are in upright
bold, with constant matrices in upper-case (𝐊,
𝐌) and random matrices holding observations in lower-case
(𝐱, 𝐲). Subscripts refer to entries in vectors
and columns in matrices. Superscripts refer to rows in matrices.  So
X_j is the j-th component of a random vector X.
For a data matrix 𝐱∈ℝ^n× m, each row
comprising one observation of m variables/features, X_j^(i) is
the j-th feature for the i-th observation.  Stacking the columns
of a matrix 𝐊=[κ_ij]_i,j∈ℝ^q× r
gives its vectorization
vec(𝐊)=(κ_11,…,κ_q1,κ_12,…,κ_q2,…,κ_1r,…,κ_qr)^⊤.
For a matrix 𝐊∈ℝ^q× q,
diag(𝐊)∈ℝ^q denotes its diagonal, and
for a vector v∈ℝ^q,
diag(v) is the q× q diagonal matrix
with diagonals v_1,…,v_q.

For a≥ 1, the ℓ_a-norm of a vector v∈ℝ^q is denoted 

    v_a=(∑_j=1^q |v_j|^a)^1/a,
 with
v_∞=max_j=1,…,q|v_j|.  A
matrix 𝐊=[κ_ij]_i,j∈ℝ^q× r has Frobenius norm

    𝐊_F≡vec(𝐊)_2≡√(∑_i=1^q∑_j=1^rκ_ij^2),

and max norm
𝐊_∞≡vec(𝐊)_∞≡max_i,j|κ_ij|.
Its ℓ_a-ℓ_b operator norm  is

    𝐊_a,b≡max_x≠0𝐊x_b/x_a

with shorthand notation
𝐊_a≡𝐊_a,a; for instance,
𝐊_∞≡max_i=1.…,q∑_j=1^r|κ_ij|. 

For a function f:ℝ^m→ℝ, we define
∂_j f(x) as the partial derivative with respect
to 
x_j, and
∂_jjf(x)=∂_j∂_j
f(x).  For f:ℝ→ℝ^m,
f(x)=(f_1(x),…,f_m(x))^⊤, we let
f'(x)=(f_1'(x), …, f_m'(x))^⊤ be the vector
of derivatives.  Likewise f”(x) is used for second
derivatives.  The symbol 1_A(·) denotes the indicator
function of the set A, while 1_n∈ℝ^n is the
vector of all 1's. For a,
b∈ℝ^m,
a∘b≡(a_1b_1,…,a_mb_m)^⊤.
A density of a distribution is always a probability density function
with respect to Lebesgue measure.  When it is clear from the context,
𝔼_0 denotes the expectation under a true distribution P_0.





§ SCORE MATCHING


In this section, we review the original score matching and develop our generalized score matching estimators.


 §.§ Original Score Matching



Let X be a random vector taking values in ℝ^m with distribution P_0 and density p_0. Let 𝒫 be a family of distributions of interest with twice continuously differentiable densities supported on ℝ^m. Suppose P_0∈𝒫. 
The score matching loss for P∈𝒫, with density p, is given by

    J(P)=∫_ℝ^mp_0(x)∇log p(x)-∇log p_0(x)_2^2x.

The gradients in (<ref>) can be thought of as gradients with respect to a hypothetical location parameter, evaluated at the origin <cit.>.




The loss J(P) is minimized if and only if P=P_0, which forms the
basis for estimation of P_0.  Importantly, since the loss depends on
p only through its log-gradient, it suffices to know p up to a
normalizing constant. Under mild conditions, (<ref>) can be
rewritten as

    J(P)=∫_ℝ^mp_0(x)∑_j=1^m[∂_jjlog p(x)+(∂_jlog p(x))^2/2]x,

plus a constant independent of p.  The integral in (<ref>)
can be approximated by a sample average; this alleviates the need for knowing the true density p_0, and provides a way to estimate p_0.



 §.§ Generalized Score Matching for Non-Negative Data


When the true density p_0 is supported on a proper subset of
ℝ^m, the integration by parts underlying the equivalence of
(<ref>) and (<ref>) may fail due to discontinuity at the boundary.  For distributions supported on the
non-negative orthant, ℝ_+^m, <cit.> addressed this
issue by instead minimizing the non-negative score matching
  loss

    J_+(P)=∫_ℝ^m_+p_0(x)∇log p(x)∘x-∇log p_0(x)∘x_2^2x.

This loss can be motivated via gradients with respect to a
hypothetical scale parameter <cit.>. Under mild conditions,
J_+(P) can again be rewritten in terms of an expectation of a
function independent of p_0, thus allowing one to form a sample
loss.

In this work, we consider generalizing the non-negative score matching
loss as follows.


  Let 𝒫_+ be the family of distributions of interest, and
  assume every P∈𝒫_+ has a twice continuously
  differentiable density supported on ℝ_+^m. Suppose the
  m-variate random vector X has true distribution
  P_0∈𝒫_+, and let p_0 be its twice continuously
  differentiable density.  Let
  h_1,…,h_m:ℝ_+→ℝ_+ be a.s. positive
  functions that are absolutely continuous in every bounded
  sub-interval of ℝ_+,
  
  
  and set h(x)=(h_1(x_1),…,h_m(x_m))^⊤. For P∈𝒫_+ with density p, the generalized
    h-score matching loss is

    J_h(P)=∫_ℝ_+^m1/2p_0(x)∇log p(x)∘h(x)^1/2-∇log p_0(x)∘h(x)^1/2_2^2x,

where h^1/2(x)≡(h_1^1/2(x_1),…,h_m^1/2(x_m))^⊤.



  
  The distribution P_0 is the unique minimizer of
  J_h(P) for P∈𝒫_+.


 
  First, observe that J_h(P)≥ 0 and
  J_h(P_0)=0.  For uniqueness, suppose
  J_h(P_1)=0 for some P_1∈𝒫_+.  Let
  p_0 and p_1 be the respective densities.  By
  assumption p_0(x)>0 a.s. and
  h_j^1/2(x)>0 a.s. for all j=1,…,m.  Therefore,
  we must have
  ∇log p_1(x)=∇log
  p_0(x) a.s., or equivalently,
  p_1(x)=const×
  p_0(x) almost surely in
  ℝ_+^m. Since p_1
  and p_0 are continuous densities supported on ℝ_+^m, it follows that
  p_1(x)=p_0(x) for
  all x∈ℝ_+^m.


Choosing all h_j(x)=x^2 recovers the loss from (<ref>).  In
our generalization, we will focus on using functions h_j that are
increasing but are bounded or grow rather slowly.



This will alleviate the need to estimate higher moments,
leading to better practical performance and improved theoretical
guarantees. 

We will consider the following assumptions:

    (A1)        p_0(x)h_j(x_j)∂_jlog p(x)|^x_j↗+∞_x_j↘ 0^+=0.,  ∀x_-j∈ℝ_+^m-1,  ∀ p∈𝒫_+;
    
    (A2)       𝔼_p_0∇log p(X)∘h^1/2(X)_2^2<+∞,    𝔼_p_0(∇log p(X)∘h(X))'_1<+∞,   ∀ p∈𝒫_+,

where ∂_jlog p(x)≡.∂log p(y)/∂ y_j|_y=x, f(x)|_x_j↘ 0^+^x_j↗+∞.≡lim_x_j↗+∞f(x)-lim_x_j↘ 0f(x), “∀ p∈𝒫_+” is a shorthand for “for all p being the density of some P∈𝒫_+”, and the prime symbol denotes component-wise differentiation. While the second half of (A2) was not made explicit in  <cit.>, (A1)-(A2) were both required for integration by parts and Fubini-Tonelli to apply.

Once the forms of p_0 and p are given, sufficient conditions for
h for Assumptions (A1)-(A2) to hold are easy to find.
In particular, (A1) and (A2) are easily satisfied and verified for exponential families.












Integration by parts yields the following theorem which shows that
J_h from (<ref>) is an expectation (under
P_0) of a function that does not depend on p_0, similar to
(<ref>). The proof is given in Appendix <ref>.


Under (A1) and (A2), the loss from (<ref>) equals

    J_h(P)=∫_ℝ_+^mp_0(x)∑_j=1^m[h_j'(x_j)∂_j (log p(x))+h_j(x_j)∂_jj(log p(x))+1/2.
    .+1/2h_j(x_j)(∂_j(log p(x)))^2]x

plus a constant independent of p.




Given a data matrix 𝐱∈ℝ^n× m with rows
X^(i), we define the sample version of (<ref>)
as

    Ĵ_h(P)=1/n∑_i=1^n∑_j=1^m{h_j'(X_j^(i))∂_j (log p(X^(i)))[1/2(∂_j(log p(X^(i))))^2].
    .+h_j(X_j^(i))[∂_jj(log p(X^(i)))+1/2(∂_j(log p(X^(i))))^2]}.



Subsequently, for a distribution P with density p, we let
J_h(p)≡ J_h(P).  Similarly, when
a distribution P_θ with density
p_θ is associated to a parameter vector
θ, we write
J_h(θ)≡
J_h(p_θ) ≡
J_h(P_θ).  We apply similar
conventions to the sample version Ĵ_h(P).  We
note that this type of loss is also treated in slightly different
settings in <cit.> and <cit.>.

In the one-dimensional case, using the notation in <cit.>, J_h(P) and Ĵ_h(P) correspond to d(P_0,P) and S(x,P), respectively, and can be generated by ϕ(x,p,p_1)≡-h(x)p_1^2/(2p) (c.f. Equations (39), (51), (53) and Section 10.1 therein). Thus Theorem <ref> follows from this correspondence. While (A1) is equivalent to the condition implied by the boundary divergence d_b=0 in that paper, (A2), which we assume for invoking Fubini-Tonelli due to multi-dimensionality, is not present. On the other hand, while <cit.> treats the multivariate case, it does not cover the connection between our J_h and Ĵ_h. Since ϕ is concave but not strictly concave in (p,p_1), the results in <cit.> only imply that P_0 is a minimizer, a weaker conclusion than Proposition <ref>.












 



§ EXPONENTIAL FAMILIES


In this section, we study the case where 𝒫_+≡{p_θ:θ∈Θ} is an
exponential family comprising continuous distributions with support
ℝ_+^m.  More specifically, we consider densities that are indexed by the canonical parameter
θ∈ℝ^r and have the form

    log p_θ(x)=θ^⊤t(x)-ψ(θ)+b(x),  x∈ℝ_+^m,

where t(x)∈ℝ_+^r comprises the sufficient statistics, ψ(θ) is a normalizing constant depending on θ only, and b(x) is the base measure, with t and b a.s. differentiable with respect to each component. Define t_j'(x)≡(∂_j t_1(x),…,∂_j t_r(x))^⊤ and b_j'(x)≡∂_j b(x).
 

 In the 1-D case h is a scalar function and (A1)-(A2) become


  * lim_x↗+∞h(x)exp(θ_0^⊤t(x)-ψ(θ)+b(x))(θ^⊤t'(x)+b'(x))=0, ∀θ∈Θ,

  * lim_x↘ 0^+h(x)exp(θ_0^⊤t(x)-ψ(θ)+b(x))(θ^⊤t'(x)+b'(x))=0, ∀θ∈Θ,

  * ∫_ℝ_+h(x)exp(θ_0^⊤t(x)-ψ(θ)+b(x))(θ^⊤t'(x)+b'(x))^2<+∞.




Under Assumptions (A1)-(A2) from Section <ref>, the empirical generalized h-score matching loss (<ref>) can be rewritten as a quadratic function in θ∈ℝ^r:

    Ĵ_h(p_θ)   =1/2θ^⊤Γ(𝐱)θ-g(𝐱)^⊤θ+const,  where
    Γ(𝐱)   =1/n∑_i=1^n∑_j=1^mh_j(X_j^(i))t_j'(X^(i))
    t_j'(X^(i))^⊤  and
    g(𝐱)   =-1/n∑_i=1^n∑_j=1^m[h_j(X_j^(i))b_j'(X^(i))
    t'_j(X^(i))+h_j(X_j^(i))t_j”(X^(i))+
    h_j'(X_j^(i))t_j'(X_i)]

are sample averages of functions of the data matrix 𝐱 only. 


Define Γ_0≡𝔼_p_0Γ(𝐱), g_0≡𝔼_p_0g(𝐱), and Σ_0≡𝔼_p_0[(Γ(𝐱)θ_0-g(𝐱))(Γ(𝐱)θ_0-g(𝐱))^⊤].


Suppose that
 
 
   (C1)  Γ is a.s. invertible, and  

  
   (C2)  Γ_0, Γ_0^-1, g_0
    and Σ_0 exist and are entry-wise finite.
  
  Then the minimizer of (<ref>) is a.s. unique with closed-form solution θ̂≡Γ(𝐱)^-1g(𝐱). Moreover, 

    θ̂→_a.s.θ_0
      and  √(n)(θ̂-θ_0)→_d𝒩_r(0,Γ_0^-1Σ_0Γ_0^-1)  as n→∞.



Theorems <ref> and <ref>
are proved in Appendix <ref>.
Theorem <ref> clarifies the quadratic nature
of the loss, and Theorem <ref> provides a basis for
asymptotically valid tests and confidence intervals for the parameter
θ.  Note that Condition (C1) holds if and only if
h_j(X_j)>0 a.s. and [t_j'(X^(1)),…,t_j'(X^(n))]∈ℝ^r× n has rank r a.s. for some j=1,…,m.

The conclusion in Theorem <ref> indicates that, similar to the estimator in <cit.> with h_j(x)=x^2, the closed-form solution for our generalized θ̂ allows one to consistently estimate the canonical parameter in an exponential family distribution without needing to calculate the often complicated normalizing constant ψ(θ) or resort to numerical methods. Computational details are explicated in Section <ref>.

Below we illustrate the estimator θ̂ in the case of univariate truncated normal distributions. We assume (A1)-(A2) and (C1)-(C2) throughout.



Univariate (m=r=1) truncated normal distributions for
mean parameter μ and variance parameter σ^2 have density 

    p_μ,σ^2(x)∝exp{-(x-μ)^2/2σ^2},   x∈ℝ_+.

If σ^2 is known but μ unknown, then writing the density in canonical form as
in (<ref>) yields

    p_θ(x)∝exp{θ t(x)+b(x)},  θ≡μ/σ^2,   t(x)≡ x,   b(x)=-x^2/2σ^2.


Given an i.i.d. sample X_1,…,X_n∼ p_μ_0,σ^2, the generalized h-score matching estimator of μ is

    μ̂_h≡∑_i=1^nh(X_i)X_i-σ^2h'(X_i)/∑_i=1^n h(X_i).

If lim_x↘ 0^+h(x)=0, lim_x↗+∞h^2(x)(x-μ_0)p_μ_0,σ^2(x)=0 and the expectations are finite (for example, when h(x)=o(exp(Mx^2)) for M<1/4σ^2), then 

    √(n)(μ̂_h-μ_0)→_d𝒩(0,𝔼_0[σ^2 h^2(X)+σ^4 h'^2(X)]/𝔼_0^2[h(X)]).

We recall that the Cramér-Rao lower bound (i.e. the lower bound on the variance of any unbiased estimator) for estimating μ_ is 
    σ^4/(X-μ_0).





  Consider the univariate truncated normal distributions from
  (<ref>) in the setting where the mean parameter μ
  is known but the variance parameter σ^2>0 is unknown.

In canonical form as in (<ref>), we write

    p_θ(x)∝exp{θ t(x)+b(x)},  θ≡1/σ^2,   t(x)≡ -(x-μ)^2/2,   b(x)=0.

Given an i.i.d. sample X_1,…,X_n∼ p_μ,σ_0^2, the generalized h-score matching estimator of σ^2 is

    σ̂_h^2≡∑_i=1^nh(X_i)(X_i-μ)^2/∑_i=1^n
      h(X_i)+h'(X_i)(X_i-μ).

If, in addition to the assumptions in Example <ref>, lim_x↗+∞h^2(x)(x-μ)^3p_μ,σ_0^2(x)=0, then 

    √(n)(σ̂_h^2-σ_0^2)→_d𝒩(0,2σ^6_0𝔼_0[h^2(X)(X-μ)^2]+σ^8_0𝔼_0[h'^2(X)(X-μ)^2]/𝔼^2_0[h(X)(X-μ)^2]).

Moreover, the Cramér-Rao lower bound for estimating σ^2_ is 
    4σ_0^8/(X-μ)^2.




  In Example <ref>, if μ_0=0, then h(x)≡ 1 also satisfies (A1)-(A2) and
  (C1)-(C2) and one recovers the sample variance
  1n∑_i X_i^2, which obtains the Cramér-Rao lower
  bound.


In these examples, there is a benefit in using a bounded function h, which can be explained as follows. When μ≫σ, there is effectively no truncation to the Gaussian distribution, and our method adapts to using low moments in (<ref>), since a bounded and increasing h(x) becomes almost constant as it reaches its asymptote for x large.  Hence, we effectively revert to the original score matching (recall Section <ref>). In the other cases, the truncation effect is significant and 
our estimator uses higher moments accordingly.





Figure <ref> plots the asymptotic
variance of μ̂_h from Example
<ref>, with σ=1
known. Efficiency as measured by the Cramér-Rao lower bound divided
by the asymptotic variance is also shown. We see that two truncated
versions of log(1+x) have asymptotic variance close to the
Cramér-Rao
bound. 
This asymptotic variance is also reflective of the variance for smaller
finite samples.


Figure <ref> is the analog of Figure
<ref> for σ̂^2_h from Example <ref> with
μ=0.5 known.



While the specifics are a bit different the benefits of using bounded
or slowly growing h are again clear.  We note that when
σ is small, the effect of truncation to the positive part of
the real line is small.

In both plots we order/color the curves based on their overall efficiency, so they have different colors in one from the other, although the same functions are presented. For all functions presented here (A1)–(A2) and (C1)–(C2) are satisfied.









§ REGULARIZED GENERALIZED SCORE MATCHING


In high-dimensional settings, when the number r of parameters to
estimate may be larger than the sample size n, it is hard, if not
impossible, to estimate the parameters consistently without turning
to some form of 
regularization. More specifically, for exponential families, condition
(C1) in Section <ref> fails when r>n. A
popular approach is then the use of ℓ_1 regularization to exploit
possible sparsity.






Let the data matrix 𝐱∈ℝ^n× m comprise n
i.i.d. samples from distribution P_0.  Assume P_0 has density
p_0 belonging to an exponential family
𝒫_+≡{p_θ:θ∈Θ},
where Θ⊆ℝ^r.  Adding an
ℓ_1 penalty to (<ref>), we obtain the
regularized generalized score matching loss 

    1/2θ^⊤Γ(𝐱)θ-g(𝐱)^⊤θ+λθ_1

as in <cit.>.  The loss in (<ref>) involves a
quadratic smooth part as in the familiar lasso loss for linear
regression.  However, although the matrix Γ is
positive semidefinite, the regularized loss in (<ref>) is not
guaranteed to be bounded unless the tuning parameter λ is
sufficiently large—a problem that does not occur in lasso.  We note
that here, and throughout, we suppress the dependence on the data
𝐱 for Γ(𝐱),
g(𝐱) and derived quantities.

For a more detailed explanation, note that that by (<ref>),

Γ=𝐇^⊤𝐇 for some
𝐇∈ℝ^nm× r.  In the high-dimensional case,
the rank of Γ, or equivalently 𝐇, is at
most nm<r.  Hence, Γ is not invertible and
g does not necessarily lie in the column span of
Γ.  Let Ker(Γ) be the
kernel of Γ.  Then there may exist
ν∈Ker(Γ) with
g^⊤ν≠0.  In this case, if

    0≤λ<sup_ν∈Ker(Γ)|
      g^⊤ν|/ν_1,
 
there exists ν∈Ker(Γ)
with
1/2ν^⊤Γν=0
and
-g^⊤ν+λν_1<0.
Evaluating at θ(a)=a·ν for
scalar a>0, the loss becomes
a(-g^⊤ν+λν_1),
which is negative and linear in a, and thus unbounded below. In this
case no minimizer of (<ref>) exists for small values of
λ. This issue also exists for the estimators from
<cit.> and <cit.>, which correspond to score
matching for GGMs.  We note that in the context of estimating the interaction matrix in pairwise models, r=m^2; thus, the condition nm<r reduces to n<m, or n<m+1 when both 𝐊 and η are estimated.

To circumvent the unboundedness problem, we add small values
γ_ℓ>0 to the diagonal entries of Γ,
which become Γ_ℓ,ℓ+γ_ℓ,
ℓ=1,…,r. This is in the spirit of work such as <cit.>
and corresponds to an elastic net-type penalty <cit.> with
weighted ℓ_2 penalty
∑_ℓ=1^rγ_ℓθ_ℓ^2.  After this
modification, Γ is positive definite, our
regularized loss is strongly convex in θ, and a
unique minimizer exists for all λ≥ 0.  For the special case
of truncated GGMs, we will show that a result on consistent estimation holds if we choose γ_ℓ=δ_0Γ_ℓ,ℓ for a suitably small constant δ_0>0, for which we propose a particular choice to avoid tuning. This choice of γ_ℓ depends on the data through Γ_ℓ,ℓ.


  For
  γ∈ℝ_+^r\{0},
  let
  Γ_γ≡Γ+diag(γ).
  The regularized generalized h-score matching
    estimator with tuning parameter λ≥ 0 and amplifier
  γ is the estimator

    θ̂ ∈_θ∈ΘĴ_h,λ,γ(θ)≡_θ∈Θ1/2θ^⊤Γ_γ(𝐱)θ-g(𝐱)^⊤θ+λθ_1.







In the case where γ=(δ-1)diag(Γ) for some δ>1, we also call δ the multiplier. We note that θ̂ from (<ref>) is a piecewise linear function of λ <cit.>.











§ SCORE MATCHING FOR GRAPHICAL MODELS FOR NON-NEGATIVE DATA



In this section we apply our generalized score matching estimator to a general class of graphical models for non-negative data.



 §.§ A General Framework of Pairwise Interaction Models





We consider the class of pairwise interaction power models with
density introduced in (<ref>).  We recall the form of
the density:

    p_η,𝐊(x)∝exp(-1/2ax^a^⊤𝐊x^a+η^⊤x^b-1_m/b)1_ℝ_+^m(x),

where a and b are known constants, and the interaction matrix
𝐊 and the vector η are parameters.  When
b=0, we use the convention that x^0-1/0≡log x and
apply the logarithm element-wise.  Our focus will be on the
interaction matrix 𝐊 that determines the conditional
independence graph through its support
S(𝐊)≡{(i,j):κ_ij≠ 0}.  However, unless
η is known or assumed to be zero, we also need to
estimate η as a nuisance parameter. In the case where we assume η≡0 is known (i.e. the linear part (x^b-1_m)/b is not present), we call the distribution (and the corresponding estimator) a centered distribution (estimator), in contrast to the general case termed non-centered when we assume η≠0 or unknown.



We first give a set of sufficient conditions for the density to be
valid, i.e., the right-hand side of (<ref>) to be
integrable.  The proof is given in Appendix <ref>.


Define conditions 


  * 𝐊 is strictly co-positive, i.e., v^⊤𝐊v>0 for all v∈ℝ_+^m\{0};

  * 2a>b>0;

  * a>0, b=0, and η_j>-1 for j=1,…,m (η≻-1_m).

In the non-centered case, if (CC1) and one of (CC2) and (CC3) holds,
then the function on the right-hand side of (<ref>) is
integrable over ℝ_+^m. In the centered case, (CC1) and
a>0 are sufficient.


We emphasize that (CC1) is a weaker condition than positive
definiteness.  Criteria for strict co-positivity are discussed in
<cit.>.




 §.§ Implementation for Different Models

In this section we give some implementation details for the
regularized generalized h-score matching estimator
defined in (<ref>) applied to the pairwise
interaction models from (<ref>).
We again let
Ψ≡(𝐊^⊤,
  η)^⊤∈ℝ^(m+1)× m.  The
unregularized loss is then

    Ĵ_h(P)   =1/2vec(Ψ)^⊤Γ(x)vec(Ψ)-g(x)^⊤vec(Ψ).




The general form of the matrix Γ and the vector
g in the loss were given in equations
(<ref>)–(<ref>).



Here Γ∈ℝ^(m+1)m× (m+1)m is block-diagonal, with the j-th ℝ^(m+1)× (m+1) block

    Γ_j(𝐱)   ≡[   Γ_11,j   Γ_12,j; Γ_12,j^⊤   Γ_22,j ]
       ≡1/n∑_i=1^n
    [ h_j(X_j^(i))X_j^(i)^2a-2X^(i)^aX^(i)^a^⊤        -h_j(X_j^(i))X_j^(i)^a+b-2X^(i)^a;      -h_j(X_j^(i))X_j^(i)^a+b-2X^(i)^a^⊤                 h_j(X_j^(i))X_j^(i)^2b-2 ]
       =1/n𝐲^⊤𝐲,  𝐲≡[ -(√(h_j(X_j))∘X_j^a-1)∘𝐱^a        √(h_j(X_j))∘X_j^b-1 ]∈ℝ^n,m+1,

where the ∘ product between a vector and a matrix means an
elementwise multiplication of the vector with each column of
the matrix, and
h_j(X_j)≡[h_j(X_j^(1)),…,h_j(X_j^(n))^⊤]∈ℝ^m.
Furthermore, g≡[ vec(𝐠_1);      g_2 ]∈ℝ^(m+1)m, where 𝐠_1 and g_2 correspond to each entry of 𝐊 and η, respectively. The j-th column of 𝐠_1∈ℝ^m× m, written as g_1,j(𝐱), is

    1/n∑_i=1^n (h_j'(X_j^(i))X_j^(i)^a-1+(a-1)h_j(X_j^(i))X_j^(i)^a-2)X^(i)^a+ah_j(X_j^(i))X_j^(i)^2a-2e_j,m,

where e_j,m is the m-vector with 1 at the j-th position and 0 elsewhere, and the j-th entry of g_2∈ℝ^m is

    g_2,j=1/n∑_i=1^n-h_j'(X_j^(i))X_j^(i)^b-1-(b-1)h_j(X_j^(i))X_j^(i)^b-2.




These formulae also hold for b=0 since Γ and g only depend on the gradient of the log density, and (̣x^b-1)/b/x̣=x^b-1 also holds for b=0. In the centered case where we know η_0≡0, we only estimate 𝐊∈ℝ^m× m, and Γ∈ℝ^m^2× m^2 is still block-diagonal, with the j-th block being the Γ_11,j submatrix in (<ref>), while g is just vec(𝐠_1). Since b only appears in the η part of the density, the formulae only depend on a in the centered case.


We emphasize that it is indeed necessary to introduce amplifiers
γ≻0 or a multiplier δ>1 in
addition to the ℓ_1 penalty.  It is clear
from (<ref>) that
rank(Γ_j)≤min{n,m+1} (or min{n,m} if centered).
Thus, Γ is non-invertible when n≤ m (or n<m if centered) and
g need not lie in its column span.  

We claim that including amplifiers/multipliers for the submatrices
Γ_11,j only is sufficient for unique existence of a solution for all penalty parameters λ≥ 0.  To see this,
consider any nonzero vector ν∈ℝ^m+1.
Partition it as
ν≡(ν_1,ν_2) with
ν_1∈ℝ^m.  Let
Γ_j,γ be our amplified version of the
matrix Γ_j from (<ref>), so

    Γ_j,γ = [ Γ_11,j+diag(γ_1,…,γ_m)                 Γ_12,j;               Γ_12,j^⊤                 Γ_22,j ].

As Γ_j itself is positive semidefinite, we find that if at least
one of the first m entries of ν is nonzero then

    ν^⊤Γ_j,γν ≥ ν^⊤Γ_jν + ∑_k=1^m
      ν_k^2γ_k  ≥ ∑_k=1^m
      ν_k^2γ_k  > 0.

If only the last entry of  ν is nonzero then

    ν^⊤Γ_j,γν = ν_m+1^2Γ_22,j > 0

almost surely; recall that
Γ_22,j=1n∑_i=1^nh_j(X_j^(i))X_j^2b-2.
We conclude that Γ_j,γ (and thus the entire amplified Γ) is a.s. positive
definite, which ensures unique existence of the loss minimizer.




Given the formulae for Γ and g, one adds the ℓ_1 penalty on Ψ to get the regularized loss (<ref>). Our methodology readily accommodates two different choices of the penalty parameter λ for 𝐊 and η. This is also theoretically supported for truncated GGMs, since if the ratio of the respective values λ_𝐊 and λ_η is fixed, the proof of the theorems in Section <ref> can be easily modified by replacing η by (λ_η/λ_𝐊)η.
  
  
  
  To avoid picking two tuning parameters, one may also
  choose to remove the penalty on η altogether by
  profiling out η and solve for
  η̂≡Γ_22^-1(g_2-Γ_12^⊤vec(𝐊̂)),
  with 𝐊̂ the minimizer of the profiled loss

    Ĵ_h,λ,γ,profile(𝐊)≡1/2vec(𝐊)^⊤Γ_γ,11.2vec(𝐊)
    -(g_1-Γ_12Γ_22^-1g_2)^⊤vec(𝐊)+λ𝐊_1,

where the Schur complement Γ_γ,11.2≡Γ_γ,11-
Γ_12Γ_22^-1Γ_12^⊤ is a.s. positive definite such that the profiled estimator exists a.s. for all λ≥ 0.  This profiled approach corresponds to choosing λ_η/λ_𝐊=0. A detailed theoretical analysis of the profiled estimator is beyond the scope of
this paper, however.  We note that in the other extreme, with λ_η/λ_𝐊=+∞, the non-centered estimator reduces to the estimator from the centered case.






The truncated normal
model comprises the density

    p_μ,𝐊(x) ∝ exp{-1/2(x-μ)^⊤𝐊(x-μ)}1_[0,∞)^m(x).













This corresponds to (<ref>) with a=b=1, and η=𝐊μ.
The 
j-th (m+1)× (m+1) block of Γ(𝐱) is

    1/n[ 𝐱^⊤diag(h_j(X_j))𝐱       -𝐱^⊤h_j(X_j);       -h_j(X_j)^⊤𝐱      h_j(X_j)^⊤1_n ].

Partitioning the vector g(𝐱) into m subvectors
g_j(𝐱)∈ℝ^m+1, where the entries of
g_j(𝐱) correspond to column
Ψ_j, the k-th entry of
g_j(𝐱) is

    g_jk(𝐱)≡1n∑_i=1^nh'_j(X_j^(i))X_k^(i)
            if  k≤ m,  k≠ j,
    1n∑_i=1^nh'_j(X_j^(i))X_k^(i)
        +h_j(X_j^(i))
            if  k= j,
    
        -1n∑_i=1^nh_j'(X_j^(i))
            if  k=m+1.






The exponential square-root graphical model in <cit.> has

    p_η,𝐊(x)∝exp(-√(x)^⊤𝐊√(x)+2η^⊤√(x))1_[0,∞)^m(x),

which corresponds to (<ref>) with a=b=1/2. We refer to this as the exponential model. In this case, the j-th ℝ^(m+1)× (m+1) block of Γ is

    Γ_j(𝐱)≡1/n∑_i=1^nh_j(X_j^(i))/X_j^(i)(-√(X^(i))
     1)(-√(X^(i))^⊤,1)

and g=vec(𝐠_0), where the j-th column of 𝐠_0∈ℝ^(m+1)× m is

    g_j(𝐱)≡1/n∑_i=1^n2h'_j(X_j^(i))X_j^(i)-h_j(X_j^(i))/2X_j^(i)^3/2(√(X^(i))
     -1)+h_j(X_j^(i))/2X_j^(i)e_j,m+1.




If a=1/2 and b=0, then (<ref>) becomes

    p_η,𝐊(x)∝exp(-√(x)^⊤𝐊√(x)+η^⊤log(x))1_(0,∞)^m(x).

If 𝐊 is diagonal in this case, then
X∼ p_η,𝐊 has independent
entries with X_j following the gamma distribution with rate
κ_jj and shape η_j+1, which gives an intuition for
condition (CC3) η_j>-1 in Theorem <ref>. We
can thus view (<ref>) as a multivariate gamma
distribution with pairwise interactions among the covariates, and call this the gamma model.  For this model, the j-th block
of Γ is

    Γ_j(𝐱)≡1/n∑_i=1^nh_j(X_j^(i))/X_j^(i)^2(-√(X_j^(i)X^(i))
     1)(-√(X_j^(i)X^(i))^⊤,1)

and the part of g corresponding to 𝐊_j is 

    g_1,j(𝐱)≡1/n∑_i=1^n2h'_j(X_j^(i))X_j^(i)-h_j(X_j^(i))/2X_j^(i)^3/2√(X^(i))+h_j(X_j^(i))/2X_j^(i)e_j,m,

while the part for η_j is 

    g_2.j(𝐱)=1/n∑_i=1^nh_j(X_j^(i))/X_j^(i)^2-h_j'(X_j^(i))/X_j^(i).



We note that the Γ_11,j sub-matrix of Γ_j and the g_1,j sub-vector of g_j for the gamma model are the same as those for the exponential model, since a=1/2 in both cases and the parts involving 𝐊 in the densities are the same.



 §.§ Computational Details


In the most general exponential family setting, as in Eq. (<ref>)–(<ref>) in Theorem <ref>, the time complexity for forming Γ∈ℝ^r× r and g∈ℝ^r is 𝒪(nm(f_b'(m)+r^2+r(f_t'(m)+f_t”(m)))). Here f_b'(m) is the average time complexity for calculating ∂_j b(x) over j=1,…,m, and similarly f_t'(m) for ∂_j t_ℓ(x) and f_t”(m) for ∂_jj t_ℓ(x) over j=1,…,m and ℓ=1…,r. In many applications, however, these three functions would be constant in m, thus giving an 𝒪(nmr^2) computational complexity, with the dominating term coming from the operations for t_j't_j'^⊤ in Γ since Γ is of dimension r× r. 

For pairwise interaction power models, r=m^2 and the formula above becomes 𝒪(nm^5). However, since Γ is block-diagonal with only m^3 nonzero entries and by the special form of t(x)=x^ax^a^⊤, the true complexity is in fact 𝒪(nm^3).

While the introduction of the ℓ_1 penalty inevitably precludes the estimator from having a closed-form solution and introduces non-differentiability, state-of-art numerical optimization algorithms, such as coordinate-descent <cit.>, can be applied for fast estimation. To speed up estimation, one can usually use warm starts using the solution from the previous λ's, as well as lasso-type strong screening rules <cit.> to eliminate components of θ̂ that are known a priori to have zero estimates.

In our implementation for pairwise interaction models of Section <ref> (that will become available in an R package), we optimize our loss functions with respect to a symmetric matrix 𝐊̂; in the non-centered case the vector η̂ is also included. We use a coordinate-descent method analogous to Algorithm 2 in <cit.>, where in each step we update each element of 𝐊̂ and η̂ based on the other entries from the previous steps, while maintaining symmetry. In our simulations in Section <ref> we always scale the data matrix by column ℓ_2 norms before proceeding to estimation. Note that estimation of 𝐊̂ without symmetry can be parallelized as the loss can be decomposed into a sum over the columns.





 §.§ Choice of the Function h

In this subsection we discuss the requirements on the function h as well as some reasonable choices of h.


  §.§.§ Requirements on h

In Section <ref>, we presented two assumptions (A1) and (A2) under which the generalized score-matching loss is valid, i.e., the integration by parts is justified and Theorem <ref> holds. In this section, we present some sufficient (and nearly necessary) requirements on h such that (A1) and (A2) are satisfied.


Suppose h:ℝ_+^m→ℝ_+^m with h(x)=(h_1(x_1),…,h_m(x_m))^⊤. We write that h∈ℋ_a,b (for simplicity we omit the dependency on m) if for all j=1,…,m:



  * h_j is absolutely continuous in every bounded sub-interval of ℝ_+, and thus has derivative h_j' a.s.;

  * h_j(x)>0 a.s. on ℝ_+;

  * h_j and h_j' are both bounded by some piecewise powers of x a.s. on ℝ_+; 

  *  





  Assume every P in the family of distribution 𝒫_+ satisfies (CC1)–(CC3) and thus has finite normalizing constants. If h∈ℋ_a,b, then (A1) and (A2) are satisfied.



In centered models, where η≡0, we can assume b=2a and iv) in the definition of ℋ_a,2a has q=1-a. For truncated GGMs, a=b=1, so <ref> in Definition <ref> is simply lim_x_j↘ 0^+h_j(x_j)=0.

In the case of b=0, η is an unknown parameter, and (CC3) requires each of its component to be greater than -1. If one has prior information on η or restricts the parameter space for η, the requirement reduces to h_j(x_j)=o(x_j^1-η_0,j) as x_j↘ 0^+.



Otherwise, it suffices to require h_j(x_j)=o(x_j^2).  Note that this is only a condition for x_j↘ 0^+, and the globally quadratic behavior of h_j(x_j)=x_j^2 from the original score matching is not needed on the entire ℝ_+, leaving
opportunities for improvements.




  §.§.§ Reasonable Choices of h

Assume a common univariate h for all components in
h. Inspired by Theorem <ref>, we consider h that behaves like a power of x both as x↗+∞ and as x↘ 0^+. Since the requirements on the two tails are separate, we can choose h to be a piecewise defined function that joins two
powers with possibly different degrees.  In other words, h(x)=min(x^p_1,cx^p_2) for some powers
p_1≥ p_2≥ 0 and constant c>0.


Only one constant c is
required since generalized score matching is invariant to scaling of
h.  In determining the exact power of p_1 we have the following
considerations:




  * In the centered case:


  * (A1) and (A2): Theorem <ref> requires that p_1≥ 1-a.

  * “Controlled Γ and g  for
  x^a”:  We propose avoiding poles at the origin for the entries of
  Γ and g.  The formula for
  Γ_11 in (<ref>) shows that to this end
  √(h(x))x^a-1 needs to have a non-negative degree  . This
  requires p_1≥ 2-2a.  
  The formula for g_1 similarly shows that h'(x)x^a-1, h(x)x^a-2
  and h(x)x^2a-2 all need to have a non-negative degree for small
  x.  This requires p_1≥ 2-a.


  * In the non-centered case, in addition to (i) and (ii),



  * (A1) and (A2): Theorem <ref> requires p_1≥max{1-a,1-b} for b>0, or 1-min_jη_0,j for b=0.

  * “Controlled Γ and g for x^b”: From the definition of Γ_22 and g_2 and by the same reasoning as above, √(h(x))x^b-1, h'(x)x^b-1 and h(x)x^b-2 need to be non-negative powers of x, thus requiring p_1≥max{2-b,2-2b}=2-b.




The choice of p_2, is only relevant for large data
points.  Our main consideration is then merely how well Γ and g concentrate on their true population values (Theorem <ref>).  From this perspective, our intuition is that p_2 should be chosen small so that the tails of the distributions of the entries of Γ and g are well-behaved.  Thus, we can choose p_2=0, in
which case h(x)=min(x^p_1,c) is a truncated power.



  * Centered case:


  * min(x^2-a,c), which satisfies both (1) and (2).

  * min(x^1-a,c), which satisfies (1) only ((A1) and (A2)).


  * Non-centered case:


  * min(x^max{2-a,2-b},c), which satisfies (1)—(4);

  * min(x^max{2-a,1-min_jη_0,j},c) when b=0, which satisfies (1)—(3);

  * min(x^2-a,c), which satisfies (1) and (2) and possibly (3) and (4);

  * min(x^2-b,c), which satisfies (3) and (4) and possibly (1) and (2);

  * min(x^max{1-a,1-b},c) when b>0, which satisfies (1) and (3) ((A1) and (A2));

  * min(x^1-min{a,min_jη_0,j},c) when b=0, which satisfies (1) and (3) ((A1) and (A2)).







 §.§ Tuning Parameter Selection

By treating the unpenalized loss (i.e., λ=0,
γ=0) as a negative log-likelihood, we may use the extended Bayesian Information Criterion (eBIC) to choose the tuning parameter <cit.>. Consider the centered case as an example.
Let Ŝ^λ≡{(i,j):κ̂_ij^λ≠ 0,i<j}, where 𝐊̂^λ be the estimate associated with tuning parameter λ. The eBIC is then

    eBIC(λ)=   -nvec(𝐊̂)^⊤Γ(𝐱)vec(𝐊̂)+2ng(𝐱)^⊤vec(𝐊̂)+|Ŝ^λ|log n+2log(p(p-1)/2 
     |Ŝ^λ|),

where 𝐊̂ can be either the original estimate associated with λ, or a refitted solution obtained by restricting the support to Ŝ^λ.




We use the eBIC instead of the ordinary BIC (Bayesian Information Criterion) since the BIC tends to choose an overly complex model when the model space is large, as encountered in the high-dimensional setting.  The extension in eBIC comes from the last term in the above display which can be motivated by a prior distribution under which the number of edges in the conditional independence graph is uniformly distributed; see also <cit.> and <cit.>.




§ THEORY FOR GRAPHICAL MODELS



In our regularized generalized score matching framework, we introduced the amplifiers/ multipliers to address the inexistence problem. We also proposed using a general function h in place of x^2 as a means to improve estimation accuracy. This section provides a theoretical analysis of these two aspects.

In Section <ref>, we present the theory for our regularized generalized score matching estimators for general pairwise interaction models before going into the details for the special cases of (truncated) GGMs.
Next, we show that a specific choice of amplifiers/multipliers yields consistent estimation without the need for tuning. This point is important even in the case of Gaussian models on all of ℝ^m. Therefore, in Section <ref> we digress from non-negative data and consider 

the original score matching of <cit.> for centered Gaussian distributions.



Finally, in Section <ref>, we derive probabilistic results for Ψ̂ based on Theorem <ref>, justifying the benefits of using a general bounded h over x^2 in the non-negative setting. As the most important models from the class of pairwise interaction power models over ℝ_+^m, we only treat truncated GGMs since they have the most tractable concentration bounds; this case  also provides a comparison to Corollary 2 in <cit.>, which uses x^2.



 §.§ Theory for Pairwise Interaction Models


The graphical models we treat are parametrized by the interaction matrix 𝐊 and the coefficients η on (x^b-1_m)/b. 
It is convenient to accommodate this setting with a matrix-valued parameter
Ψ∈ℝ^r_1× r_2 (in place of
θ) and specify our regularized
h-score matching loss as

    Ĵ_h,λ,γ(Ψ)≡_Ψ∈ℝ^r_1× r_21/2vec(Ψ)^⊤Γ_γ(𝐱)vec(Ψ)-g(𝐱)^⊤vec(Ψ)+λΨ_1.




In the non-centered case we thus take Ψ=[𝐊,η]^⊤∈ℝ^m(m+1)× m. In the centered case, Ψ is simply the m× m interaction matrix 𝐊. Following related prior
work such as <cit.>, for ease of proof we allow the matrix 𝐊 to be
nonsymmetric, which allows us to decouple optimization over the
different columns of 𝐊 or Ψ, while in our implementations we ensure that 𝐊 is symmetric.


  Let
  Γ_0≡𝔼_0Γ(𝐱)
  and g_0≡𝔼_0g(𝐱)
  be the population versions of Γ(𝐱) and
  g(𝐱) under the distribution given by a true
  parameter matrix Ψ_0.  The support of a matrix
  Ψ is
  S(Ψ)≡{(i,j):ψ_ij≠ 0}, and we let
  S_0=S(Ψ_0).  For a matrix Ψ_0, we
  define d_Ψ_0 to be the maximum number of non-zero
  entries in any column, and
  c_Ψ_0≡Ψ_0_∞,∞.
  Writing Γ_0,AB for the A× B submatrix of
  Γ_0, we define
  
    c_Γ_0≡(Γ_0,S_0S_0)^-1_∞,∞.

  Finally, Γ_0 satisfies the
    irrepresentability condition with incoherence parameter
    α∈(0,1] and edge set S_0 if

    Γ_0,S_0^cS_0(Γ_0,S_0S_0)^-1_∞,∞≤ (1-α).








Our analysis of the regularized generalized h-score
matching estimator builds on the following theorem taken from
<cit.>.  









 
  
  Suppose Γ_0 has
  Γ_0,S_0S_0 invertible and satisfies the
  irrepresentability condition (<ref>)
  with
  incoherence parameter α∈(0,1]. Assume that

    Γ_γ(𝐱)-Γ_0_∞<ϵ_1,        g(𝐱)-g_0_∞<ϵ_2,

with d_Ψ_0ϵ_1≤α/(6c_Γ_0). If 

    λ>3(2-α)/αmax{c_Ψ_0ϵ_1,ϵ_2},

then the following holds:


  * The regularized generalized h-score matching estimator Ψ̂ minimizing (<ref>) is unique, with support Ŝ≡ S(Ψ̂)⊆ S_0, and satisfies

    Ψ̂-Ψ_0_∞   ≤c_Γ_0/2-αλ.


  * If 

    min_1≤ j<k≤ m|Ψ_0,jk|>c_Γ_0/2-αλ,

then Ŝ=S_0 and sign(Ψ̂_jk)=sign(Ψ_0.jk) for all (j,k)∈ S_0.


This result is deterministic, and the improvement of our generalized
estimator over the one in <cit.> is in its probabilistic guarantees, as shown for truncated GGMs in Theorems <ref> and <ref> in Section <ref>.  Before going into these examples, we
state a general corollary.  






  Under the assumptions of Theorem <ref>, the matrix
  Ψ̂ minimizing (<ref>)
  satisfies

    Ψ̂-Ψ_0_F     ≤ c_Γ_0/2-αλ√(|S_0|) ≤ c_Γ_0/2-αλ√(d_Ψ_0m),
    Ψ̂-Ψ_0_2    ≤ c_Γ_0/2-αλmin(√(|S_0|),d_Ψ_0).


 




 §.§ Revisiting Gaussian Score Matching



In this section we consider estimating the inverse covariance matrix
𝐊 of a centered Gaussian distribution N(0,𝐊), which of course has
density proportional to (<ref>) on all of
ℝ^m.   As shown, e.g., in Example 1 of <cit.>,  the
ℓ_1-regularized score matching loss then takes the form

    1/2tr(𝐊𝐊𝐱𝐱^⊤)-tr(𝐊)+λ𝐊_1,

which can be written as (<ref>) with θ=vec(𝐊), Γ=diag(𝐱𝐱^⊤,…,𝐱𝐱^⊤) and g=vec(𝐈_m). Thus, in general, the kernel of Γ need not be orthogonal to g, and for λ small the loss can be unbounded below as discussed above. Hence, an amplifier/multiplier on the diagonals of Γ is needed. We have the following theorem on the estimator using the amplification.


Suppose the data matrix 𝐱 holds n i.i.d. copies of
X∼N(0,𝐊_0).
Adopt the amplifying in Section <ref> and redefine the loss in (<ref>) as

    1/2tr(𝐊𝐊𝐆)-tr(𝐊)+λ𝐊_1,  𝐆_jk=(𝐱𝐱^⊤)_jk(1_{j≠ k}+(δ-1)1_{j=k}),

where 1<δ<2-(1+80√(log m/n))^-1.  Let 𝐊̂ be the resulting estimator.  Let c^*≡ 12800(max_jΣ_0,jj)^2 and c_1=4c_Γ_0/α.  If for some τ>2, the regularization parameter and the sample size satisfy

    λ   >(2c_𝐊_0(2-α)√(c^*(τlog m+log
      4)/n))/α,
    
      n   >max(c^*c_1^2d_𝐊_0^2,2)(τlog m+log 4),

then
     𝐊̂-𝐊_0_∞≤c_Γ_0/2-αλ
     with probability 1-m^2-τ.

In Corollary 1 of <cit.> the same results were shown with c^*≡ 3200(max_jΣ_0,jj)^2 when a unique minimizer exists, but the existence was not guaranteed.



 §.§ Generalized Score Matching for Truncated GGMs

Next, we provide theory for the regularized generalized h-score matching estimator Ψ̂ in the special case of truncated GGMs. Again, assume a common h for all components in
h.

  Suppose the data matrix 𝐱 holds n i.i.d. copies of
  X∼TN(0,𝐊_0), where the mean parameter is known to be zero.
 Assume that h∈ℋ_1,1 and that 0≤ h≤ M, 0≤ h'≤ M' a.s. for constants M,M', and choose γ=(δ-1)diag(Γ) with
  
    1<δ<C(n,m)≡ 2-(1+4emax{6log m/n,√(6log
            m/n)})^-1.

  Suppose that the Γ_0,S_0S_0 block of Γ_0 is invertible and Γ_0 satisfies the irrepresentability condition (<ref>) with α∈(0,1] and true edge set S_0. Define
  c_X≡
  2max_j(2√((𝐊_0^-1)_jj)+√(e) 𝔼_0
  X_j).  If
  for τ>3 the sample size and the regularization
  parameter satisfy

    n   >𝒪(τlog mmax{M^2c_Γ_0^2c_X^4d_𝐊_0^2/α^2,Mc_Γ_0c_X^2d_𝐊_0/α}),
    λ   >𝒪[(Mc_𝐊_0c_X^2+M'c_X+M)(
    √(τlog m/n)+τlog m/n)],

 then the following statements hold with probability 1-m^3-τ:


  * The regularized generalized h-score matching estimator 𝐊̂ that minimizes (<ref>) is unique, has its support included in the true support, Ŝ≡ S(𝐊̂)⊆ S_0, and satisfies

    𝐊̂-𝐊_0_∞   ≤c_Γ_0/2-αλ,
      
    𝐊̂-𝐊_0_F   ≤c_Γ_0/2-αλ√(|S_0|),
    𝐊̂-𝐊_0_2   ≤c_Γ_0/2-αλmin(√(|S_0|),d_𝐊_0),

where c_Γ_0 is defined in (<ref>).

  * Moreover, if

    min_j,k:(j,k)∈ S_0|κ_0,jk|>c_Γ_0/2-αλ,

then Ŝ=S_0 and sign(κ̂_jk)=sign(κ_0,jk) for all (j,k)∈ S_0.



The theorem is proved in Appendix <ref>, where details on the
dependencies on constants are provided.

A key ingredient of the proof is a tail bound on
Γ_γ-Γ_0_∞, which features products of the X_j^(i)'s. In <cit.>, the products are up to fourth order.  Using bounded h, our products automatically calibrates to a quadratic polynomial when the observed values are large, and resort to higher moments only when they are small.  This leads to improved bounds and convergence rates, underscored in the new requirement on the sample size n, which should be compared to n≥𝒪(d_𝐊_0^2(log m^τ)^8) in <cit.>.


For the non-centered case, by definition,
c_Ψ_0≡Ψ_0^⊤_∞,∞≤
c_𝐊_0+η_0_∞,
d_Ψ_0≤ d_𝐊_0+1.  The proof given for
Theorem <ref> goes through again here, and we have the
following consistency results.  


  Suppose the data matrix holds n i.i.d. copies of
  X∼TN(μ_0,𝐊_0).
  Assume that h∈ℋ_1,1 and that 0≤ h≤ M, 0≤ h'≤ M' a.s. for constants M,M'.  Let γ be a vector of amplifiers that are
  non-zero only for the diagonal entries of the matrices
  Γ_11,j, amplifying those by
  (δ-1)diag(Γ_11,j) with
  
    1<δ<C(n,m)≡ 2-(1+4emax{6log m/n,√(6log
            m/n)})^-1.

  Suppose further that Γ_0,S_0S_0 is invertible
  and satisfies the irrepresentability condition
  (<ref>) with α∈(0,1].  Define
  c_X≡2max_j(2√((𝐊_0^-1)_jj)+√(e) 𝔼_0
  X_j).  Suppose for τ>3 the sample size and the regularization
  parameter satisfy

    n   >𝒪(τlog mmax{M^2c_Γ_0,Ψ_0^2c_X^4d_Ψ_0^2/α^2,Mc_Γ_0,Ψ_0c_X^2d_Ψ_0/α}),
    λ   >𝒪[(Mc_Ψ_0c_X^2+M'c_X+M)(
    √(τlog m/n)+τlog m/n)],

where c_Γ_0,Ψ_0 is c_Γ_0 as in (<ref>)  but with notation Ψ_0 to differentiate it from the centered case.
Then the following statements hold with probability 1-m^3-τ:


  * The regularized generalized h-score matching estimator Ψ̂ that minimizes (<ref>) is unique, has its support included in the true support, Ŝ≡ S(Ψ̂)⊆ S_0, and satisfies

𝐊̂-𝐊_0_∞   ≤c_Γ_0,Ψ_0/2-αλ,      η̂-η_0_∞≤c_Γ_0,Ψ_0/2-αλ,

𝐊̂-𝐊_0_F   ≤c_Γ_0, Ψ_0/2-αλ√(|S_0|),      η̂-η_0_F≤c_Γ_0, Ψ_0/2-αλ√(|S_0|),

𝐊̂-𝐊_0_2   ≤c_Γ_0,Ψ_0/2-αλmin(√(|S_0|),d_Ψ_0),        η̂-η_0_2≤c_Γ_0,Ψ_0/2-αλmin(√(|S_0|),d_Ψ_0).


  * Moreover, if

    min_j,k:(j,k)∈ S_0|κ_0,jk|>c_Γ_0/2-αλ    and    min_j:(m+1,j)∈ S_0|η_0,j|>c_Γ_0/2-αλ,

then Ŝ=S_0 and sign(κ̂_jk)=sign(κ_0,jk) for all (j,k)∈ S_0 and sign(η̂_j)=sign(η_0j) for (m+1,j)∈ S_0.




  The quantity c_X in Theorem <ref>
  depends on 𝔼_0X_j, which in turn depends on the
  structure of both μ_0 and 𝐊_0.  If
  μ_0,j is large compared to (𝐊_0)^-1_jj, then
  c_X seems to scale as μ_0, which
  negatively impacts the guarantees stated in
  Theorem <ref>.  However, as in the one-dimensional case
  for estimation of μ_0 (Example <ref>), our
  estimator should automatically adapt to the large mean parameter.
  This suggests that it might be possible to improve our analysis
  involving c_X.
  
  
  
  
  
  
  
  





§ NUMERICAL EXPERIMENTS


In this section, we compare the performance of our estimator with different choices of h to the existing approaches for pairwise interaction power models. In our simulation experiments, we consider m=100 variables and n=80 and n=1000 samples, corresponding to high- and low-dimensional settings. We also tried intermediate sample sizes between these two extremes, but found no interesting result worth reporting. For n=80, amplification is necessary. Except in Section <ref>, the amplifier is set based on Theorem <ref> to δ=C(n,m)=1.8647 for truncated GGMs. The same amplifier is also used for settings with other a and b. For n=1000, we consider δ=1, i.e., no amplification, and δ=C(n,m)=1.6438 (again, based on Theorem <ref>). Throughout, we assume a common univariate h for all components in h.












 §.§ Structure of K

The underlying interaction matrices are selected as follows:
Proceeding as in Section 4.2 of <cit.>, the graph is chosen to have 10 disconnected subgraphs, each containing m/10 nodes. Thus, 𝐊_0 is block-diagonal. In each block, each lower-triangular element is set to 0 with probability 1-π for some π∈(0,1), and is otherwise drawn from Uniform[0.5,1]. The upper triangular elements are determined by symmetry.  The diagonal elements of 𝐊_0 are chosen as a common positive value such that the minimum eigenvalue of 𝐊_0 is 0.1.





We generate 5 different true precision matrices 𝐊_0, and run 10 trials with each of these precision matrices.  For n=1000, we choose π=0.8, 
which is in accordance with <cit.>. 
For n=80, we set π=0.2. 
This way n/(d_𝐊_0^2log m) is roughly constant; recall Theorems <ref> and <ref> for truncated GGMs.

In Appendix <ref>, we report results on Erdös-Rényi
  graphs, which lead to similar conclusions.




 §.§ Truncated GGMs

Given our focus on truncated GGMs and their relevance in graphical modeling applications, we start with experiments for these models.




  §.§.§ Choice of h


Our estimator requires choosing a function
h:ℝ_+^m→ℝ_+^m. For simplicity, we will always specify h(x)=(h(x_1),…,h(x_m)) for a single non-decreasing univariate function
h:ℝ_+→ℝ_+, i.e. all coordinates share the same h function.

As previously explained, h∈ℋ_a,b is a sufficient condition for assumptions (A1)-(A2), as well as (C1)-(C2) in the case of unregularized estimators. Only in the proofs of our theoretical guarantees in Section <ref> for truncated GGMs, did we require h to be bounded and to have bounded derivatives. As motivated by the discussion in Section <ref>, we consider truncated and untruncated powers, min(x,c) and x (since 2-a=2-b=1); we evaluate this choice by contrasting them with powers x^1.5 and x^2. 





We also explore functions like log(1+x) that seem natural and are linear near 0. In particular, we make a further comparison to functions linear near 0 with a finite asymptote as x↗+∞ but differentiable everywhere: MCP- <cit.> and SCAD-like <cit.> functions defined below. The results we report are based on selections of best performing choices of h.

    SCAD(x;λ,γ)≡λ x    if 0≤ x≤λ,
    2γλ x-x^2-λ^2/2(γ-1)   if λ<x<γλ,
    λ^2(γ+1)/2   if x≥γλ;MCP(x;λ,γ)≡λ x-x^2/2γ   if 0≤ x≤γλ, 
    1/2γλ^2    if x>γλ.



We do not observe any clear relationship between features such as convexity, differentiability or the slope of h at 0, and performance of the estimator. Nonetheless, for many choices of rather simple functions h, our estimator provides a significant improvement over existing methods. In particular, most h functions that behave linearly for small x, namely log(1+x) and x and their truncations, and additionally MCP and SCAD, always perform better than x^1.5 and x^2. This agrees with our discussion in  Section <ref>, where 2-a=1 is a reasonable choice of the power for small x; also see Section <ref>. However, we conclude that there is no real gain from making the function smoother by using MCP or SCAD.















Truncated Centered GGMs:For data from a truncated centered Gaussian distribution, we compare our generalized score matching estimator 
with various choices of h, to SpaCE JAM (SJ,
), which estimates graphs using additive models for conditional means, a pseudo-likelihood method SPACE
<cit.> in the reformulation of <cit.>, graphical lasso (GLASSO, ), the neighborhood selection estimator (NS) of <cit.>, and nonparanormal SKEPTIC <cit.> with Kendall's τ. Recall that the choice of h(x)=x^2 corresponds to the estimator from <cit.>. 







The ROC (receiver operating characteristic) curves for different estimators are shown in Figure <ref> on Page plot_chooseh_centered. Each plotted curve corresponds to the average of 50 ROC curves, where the averaging is based on the vertical averaging from Algorithm 3 in <cit.>, and is mean AUC-preserving. The x and y axes of each ROC curve represent the false positive and true positive rates at varying levels of penalty parameter λ, defined as


    FPR≡|Ŝ_off\ S_0,off|/m(m-1)-|S_0,off|    and    TPR≡|Ŝ_off∩ S_0,off|/|S_0,off|,

where S_0,off≡{(i,j):i≠ j∧κ_0,ij≠ 0}, and Ŝ_off≡{(i,j):i≠ j∧κ̂_ij≠ 0}.

To reduce clutter, we only report the results for the top performing competing methods. In particular, results for nonparanormal SKEPTIC are omitted, as the method always performs the worst in our experiments. The corresponding means and standard deviations of AUCs (areas under the curves) over 50 curves are given in Table <ref>. 

Looking at the mean AUCs, with the standard deviations in mind, all choices of h considered here perform better than h(x)=x^2 from <cit.> and <cit.> and the competing methods. The results for n=1000 in Table <ref> also show that the multiplier does help improve the AUCs, a matter to be discussed in Section <ref>.



Truncated Non-Centered GGMs: We generate data from a truncated noncentered Gaussian distribution with both parameters μ and 𝐊
unknown. In each trial, we form the true 𝐊_0 as in
Section <ref>, and generate each component of
μ_0 independently from the normal distribution with
mean 0 and standard deviation 0.5.



We compare the performance of our profiled estimator based on (<ref>), with different h functions, but with no penalty on η≡𝐊μ, to SPACE, SpaCE JAM (SJ), GLASSO, and neighborhood selection (NS). As before, we consider 50 trials. Representative ROC curves are plotted in Figure <ref>, and the corresponding AUCs are summarized in Table <ref>.








Even without tuning the extra penalty parameter on η≡𝐊μ, our profiled estimator beats the competing methods by a large margin when n=80. 
With multipliers 1 and n=1000, our estimators still do better than Space JAM and GLASSO, and have performance comparable to other competing methods. It might appear that the performance of our estimators deteriorate with a multiplier larger than 1; however, as we will see, there can be significant improvement in AUCs if we tune an additional parameter for the multiplier.
As in the centered case, the leading h functions in each category perform similarly, and the exact choice is not crucial.  Subsequently, we will simply use h(x)=min(x,3).









  §.§.§ Choice of multiplier

Truncated Centered GGMs:In Figure <ref>, the ROC curves for GLASSO, SPACE, and our estimator with h(x)=min(x,3), but with different levels of amplification, via different choices of multipliers δ, are compared for the centered case of Section <ref>. 

While Theorem <ref> guarantees consistency only for δ<C(n,m), we observe that there can be a gain from going beyond the
upper-bound multiplier C(n,m), which is 1.8647 for n=80 and 1.6438 for n=1000 (when n=1000, C(n,m) turns out to be the best-performing multiplier). However, the effect deteriorates fast as the multiplier grows larger. The figure suggests that while some additional gains are possible by tuning over the choice of multiplier, the upper-bound multiplier is a good default.











Truncated Non-Centered GGMs
:In Figure <ref>, we consider the non-centered case of  Section <ref>, and use the non-profiled estimator; that is, the non-centered estimator with ℓ_1 penalty on both 𝐊 and η≡𝐊μ. The ROC curves are compared to competing methods GLASSO and SPACE. For the choice of amplification in our estimator, we consider the upper-bound multiplier C(n,m) from Theorem <ref> as the default. We refer to this as high amplification.  We also consider lower amplification, with δ=2-(1+24elog m/n)^-1, referred to as medium.  For n=1000, we also consider a low multiplier 1, which corresponds to no amplification.  We compare these possible defaults to a finer grid of multipliers of which we show some representatives in the plots.






We see that among our defaults, the upper-bound choice C(n,m)
performs best. Some additional gains are possible by tuning
the multiplier over a grid of values containing this choice.
Moreover, we see that it can be beneficial to tune over both λ_𝐊 and λ_𝐊/λ_η.



We remark that while for each run, the best model picked by BIC falls on the ROC curve, a few squares are off the curve in Figure <ref> (c). This is because these squares correspond to the average of the true and false positive rates of the chosen BIC models over 50 runs, potentially due to multimodality of the distribution of the models. Nonetheless, in all cases, the average of the models picked by BIC tuned over both λ_𝐊 and λ_𝐊/λ_η looks reasonable.







 §.§ Other a/b Models

We now turn to the non-Gaussian (a≠ 1 or b≠ 1) setting. Based on the observations in Section <ref>, we focus on functions of type min(x^p,c) for some power p>0 and truncation point c>0. For simplicity, for the non-centered models we use the profiled estimator (<ref>) (i.e., λ_η=0) and use the multiplier C(n,m) in Theorem <ref> for truncated GGMs as a guidance. We note that tuning over the λ_η parameter and the multiplier can potentially give a significant improvement as seen in Section <ref>.

These simulations suggest that among the class of functions of the form min(x^p,c), x^2-a or min(x^2-a,c) with a moderately large c can be used as the default choice of h(x). This agrees with our findings in Section <ref>. We note that bounded h functions were only used in the proof for truncated GGMs, and picking a moderately large truncation point can correspond to having an untruncated power.



  §.§.§ Exponential Setting

For the exponential models, a=b=1/2. Since a=b, for both centered and non-centered settings, based on the principle in Section <ref>, choosing h(x)=min(x^3/2,c) satisfies (A1) and (A2) and also ensures that entries in Γ and g are bounded (for small x), while choosing h(x)=min(√(x),c) only guarantees (A1) and (A2).

In Figure <ref>, we present the AUCs for the ROC curves of edge recovery with different choices of h(x)=min(x^pow,c). As before, we set n=80 or 1000 and m=100, but we use an η_0 with each component uniformly equal to -0.5, 0 or 0.5; for η_0≡0, we assume this information is known and use the centered estimator. The results suggest that pow=3/2=2-a is the best choice of power. For this optimal choice, the performance improves with larger c, so x^2-a gives the best results. For sub-optimal powers, including truncation gives better results.




  §.§.§ Gamma Setting

The centered gamma models reduce to the centered exponential models. Thus, in this section, we only consider the non-centered settings, with  a=1/2, b=0. From Section <ref>, we have the following choices:


  * min(x^2,c) both satisfies (A1)–(A2) and ensures Γ and g are bounded;

  * min(x^max{3/2,1-min_jη_0,j},c) ensures (A1)–(A2) and bounds Γ_11 and g_1; by default without prior information on η_0 this is min(x^2,c);

  * min(x^3/2,c) satisfies both conditions on the interaction part only (x^a), but does not guarantee (A1)–(A2);

  * min(x^1/2,c) satisfies the sufficient conditions for (A1)–(A2) on the interaction only.

The results are shown in Figure <ref>, where we consider n=80, 1000, and η=± 0.51_100. They suggest that pow=2-a=1.5 works consistently well, although slightly outperformed by 1 and 1.25 in one case. As in the exponential case, with the optimal power it is beneficial to choose a large truncation point, or work with an untruncated power x^1.5. We conclude that the performance is likely only dependent on the (2-a) power requirement for the x^a^⊤𝐊x^a part or 2-min_jη_0,j; simulations in the next section rule out the possibility of the latter.







  §.§.§ Other Choices of a and b

In this section, we consider other choices of a and b. Specifically, a=3/2 and b=1/2 or 0. These combinations are chosen to confirm, in a more extreme setting, that the performance is mainly determined by requirements on the power based on a, which correspond to choosing a power of 1-a or 2-a, but not those on b (or on η when b=0) that correspond to 1-b and 2-b. The relationship between these  two settings is analogous to that between the exponential and gamma models (same a, b nonzero/zero). 

The results are shown in Figures <ref> and <ref>, and indeed confirm that x^2-a=x^0.5 consistently gives the optimal results, even though η^⊤x^b is in favor of x^2-b=x^1.5 for b=0.5, and η^⊤log(x) is in favor of x^2 or at least x^1-min_jη_0,j when b=0. 
There are two possible explanations for the optimality of 2-a over max{2-a,2-b} or max{2-a,1-min_jη_0,j}: (1) The AUC metric is measured only on our interest, edge recovery for the interaction matrix, which only depends on x^a; (2) using the profiled estimator weakens the effect of b.











 §.§ RNAseq Data


In this section we apply our regularized generalized h-score matching estimator for truncated non-centered GGMs to RNAseq data also studied in <cit.>, since the same model is considered therein. The data consists of n=487 prostate adenocarcinoma samples from The Cancer Genome Atlas (TCGA) data set.  Following <cit.>, we focus on m=333 genes that belong to the known cancer pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and that have no more than 10% missing values.
Missing values are set to 0.  We choose h(x)=min(x,3) and use the
upper-bound multiplier (high), as discussed in Section
<ref>.  For simplicity, we use the profiled
estimator, and choose the regularization parameter λ so that
the estimated graph has exactly m=333 edges, all these choices being
as in <cit.>. 

We compare our graph to the one in <cit.>, which corresponds to h(x)=x^2 with no multiplier. Shown in Figure <ref> are the estimated graphs, with their intersection in the middle. To improve visualization, isolated nodes are removed and the layouts are optimized for each plot. Red-colored points are the “hub nodes”, namely nodes with degree at least 10. In Figure <ref>, we plot the same graphs in a fixed layout optimized for the graph corresponding to h=min(x,3), and include the isolated nodes.

Out of 333 edges, the two estimated graphs share 117 edges in
common. Assuming that edges are placed at random between nodes and the
two graphs are independent, the distribution of the number R of
common edges follow a hypergeometric distribution, so
P(R=r)=([ m; r ])([ m(m-1)/2-m;        m-r ])/([ m(m-1)/2;        m ]).   For m=333 the probability of at
least 117 common edges is essentially zero. The large number of shared edges between the two methods can be explained by the fact that they both minimize the same underlying score-matching loss.






The graph using h(x)=min(x,3) has much more isolated nodes (204) than the other (108), and has a slightly smaller max degree (16 versus 19). 

Table <ref> provides another way of comparing between the two graphs by listing the genes with the highest node degrees.

























































In Table <ref> we list the top ten genes in terms of node degree for both estimated graphs. Due to ties, 13 genes are listed for h(x)=min(x,3) and 11 for <cit.>. As noted in <cit.>, genes with high node-degrees  are known to be important in biological networks <cit.>. Among these top genes, six are common in both graphs, and are discussed in <cit.>. We next elaborate on the evidence supporting the first four of the newly discovered genes.



  * MMP2 (Matrix metalloproteinase 2): According to <cit.>,
  increased MMP-2 expression is an independent predictor of decreased
  prostate cancer disease-free survival. <cit.> state that
  activity of MMP-2 can be useful in diagnosis, therapy, and
  assessment of malignant progression in prostate cancer.
  

  * GLI2 (GLI family zinc finger 2): GLI2 is a primary mediator of the hedgehog signaling pathway, which has been reported in prostate cancer, and plays a critical role in the malignant phenotype of prostate cancer cells <cit.>. Its increased level of expression is also related to AI prostate cancer, and may be a therapeutic target in castrate-resistant prostate cancer <cit.>.

  * LAMA4 (Laminin subunit alpha 4): LAMA4 is consistently upregulated in benign prostatic hyperplasia when compared to normal prostate tissues <cit.>.

  * RASSF5 (RAS association domain family member 5): The combination of RASSF5 along with four other DNA methylation markers can effectively differentiate between benign prostate biopsy cores from non-cancer patients and cancer cores, and can be used to identify patients at risk without repeat biopsies <cit.>.


We note that the two methods indeed use different estimators (different h functions and multipliers), and it is thus not surprising to see that some of the top genes by one method are not among those for the other. In particular, CCNE2, BRCA2, SKP2 and STAT5B, while previously reported as newly discovered in <cit.>, are dropped by our new analysis. Testing and inference (potentially using bootstrapping) is an important problem but is beyond the scope of this paper.








§ DISCUSSION


In this paper, we proposed a generalization of the score matching
estimator of <cit.>, based on scaling the log-gradients to be
matched with a suitably chosen function h.  The
generalization retains the advantages of Hyvärinen's method:
Estimates can be computed without knowledge of normalizing constants,
and for canonical parameters of exponential families, the estimation
loss is a quadratic function.

For high-dimensional exponential family graphical models, following
 <cit.>, we add an ℓ_1 penalty to regularize
the generalized score matching loss.



One practical issue that is overlooked in <cit.> is the fact
that the score matching loss can be unbounded below for a small tuning
parameter, when the dimension m exceeds the sample size
n. We fix this issue by amplifying the diagonal entries in the
quadratic matrix in the definition of the generalized score matching
loss by a factor/multiplier, and we give an upper bound on that
multiplier that guarantees consistency.



As examples we consider pairwise interaction power models on
the non-negative orthant ℝ_+^m.  Specifically, the
considered models are exponential families in which the log density is
the sum of pairwise interactions between entries in of powers
x^a plus linearly weighted effects x^b,
or log(x) when b=0.  Our main interest is in the
matrix of interaction parameters whose support determines the
distributions' conditional independence graph.  The considered
framework covers truncated normal distributions (a=b=1),
exponential square root graphical models (a=b=1/2) from
<cit.>, as well as a class of multivariate gamma distributions
(a=1/2, b=0).  

In the case of multivariate truncated normal distributions, where
the conditional independence graph is given by the underlying Gaussian
inverse covariance matrix, the sample size required for the
consistency of our method using bounded h is
Ω(d^2log m), where d is the degree of the graph.  This
matches the rates for Gaussian graphical models in <cit.> and
<cit.>.  In contrast, the sample complexity for truncated
Gaussian models given in <cit.> is Ω(d^2log^8 m).

For the considered class of pairwise interaction models, we recommend
using the function h with coordinates
h_j(x)=min(x^2-a,c) for some moderately large c, or
simply h_j(x)=x^2-a.  While this choice is effective, it would be an
interesting problem for future work to develop a method that
adaptively chooses an optimized function h from data.







This work was partially supported by grant DMS/NIGMS-1561814 from the National Science Foundation (NSF). AS also gratefully acknowledges funding by grant R01-GM114029 from the National Institute of Health (NIH).






§ PROOFS




 §.§ Proof of Theorem <ref>

The following integration by parts lemma is used in the proof of Theorem <ref>.









Let f,g:ℝ_+→ℝ be functions that are absolutely continuous in every bounded sub-interval of ℝ_+. Then

    lim_x↗+∞f(x)g(x)-lim_x↘0^+f(x)g(x)
    =∫_0^∞f(x)g̣(x)/x̣x̣+∫_0^∞g(x)f̣(x)/x̣x̣.





This is an analog of Lemma 4 from <cit.> that can be proved by integrating the partial derivatives, and follows from the fundamental theorem of calculus for absolutely continuous functions and the product rule. In particular, we work on integrals in bounded [0,c], where the product of two absolutely continuous functions in a bounded interval is again absolutely continuous, and the result is then obtained by letting c↗+∞.



Recall the following assumptions from Section <ref>:

!

    (A1)        p_0(x)h_j(x_j)∂_jlog p(x)|^x_j↗+∞_x_j↘ 0^+=0.,   ∀x_-j∈ℝ_+^m-1,  ∀ p∈𝒫_+;
    
    (A2)       𝔼_p_0∇log p(X)∘h^1/2(X)_2^2<+∞,    𝔼_p_0(∇log p(X)∘h(X))'_1<+∞,   ∀ p∈𝒫_+.




Without explicitly writing the domains ℝ_+ or ℝ_+^m in all integrals, by (<ref>) we have

J_h(p)   =1/2∫p_0(x)[∇logp(x)∘h^1/2(x)_2^2.

       .-2(∇logp(x)∘h^1/2(x))^⊤(∇logp_0(x)∘h^1/2(x))
+∇logp_0(x)∘h^1/2(x)_2^2]x

   =1/2∫p_0(x)∑_j=1^mh_j(x_j)(∂logp(x)/∂x_j)^2x_≡ A+1/2∫p_0(x)∑_j=1^mh_j(x_j)(∂logp_0(x)/∂x_j)^2x_≡ C

           -∫p_0(x)∑_j=1^mh_j(x_j)∂logp(x)/∂x_j∂logp_0(x)/∂x_jx_≡ B,

where A will simply appear in the final display as is, C is a constant as it only involves the true pdf p_0, and we wish to simplify B by integration by parts. We can split the integral into these three parts since A and C are assumed finite in the first part of (A2), and the integrand in B is integrable since |2ab|≤ a^2+b^2. Thus, by linearity and Fubini's theorem, we can write












    B   =-∑_j=1^m∫ p_0(x)h_j(x_j)∂log p(x)/∂ x_j∂log p_0(x)/∂ x_jx
       =-∑_j=1^m∫[∫ p_0(x)h_j(x_j)∂log p(x)/∂ x_j∂log p_0(x)/∂ x_jx̣_j]x_-j.

By the fact that ∂log p_0(x)/∂ x_j=1/p_0(x)∂ p_0(x)/∂ x_j, this can be simplified to

    B=-∑_j=1^m∫[∫∂ p_0(x)/∂ x_jh_j(x_j)∂log p(x)/∂ x_jx̣_j]x_-j.


But, we assume p_0 and p are twice continuously differentiable, for every j=1,…,m and fixed x_-j∈ℝ_+^m-1. Hence, in every bounded sub-interval of ℝ_+, p_0(x_-j;x_j) is an absolutely continuous function of x_j, ∂_jlog p(x_-j,x_j)=∂_j p(x_-j,x_j)/p(x_-j,x_j) is a continuously differentiable (and hence absolutely continuous) function of x_j by the quotient rule. Thus h_j(x_j)∂_jlog p(x_-j;x_j) is also absolutely continuous by the absolute continuity assumption on h_j. Then, by Lemma <ref>, where we take f≡ p_0(x_-j;x_j) and g≡ h_j(x_j)∂_jlog p(x_-j;x_j) as functions of x_j, followed by assumption (A1),

    B   =-∑_j=1^m∫[lim_a↗+∞,b↘ 0^+[p_0(x_-j;a)h_j(a)∂_jlog p(x_-j,a)-p_0(x_-j;b)h_j(b)∂_jlog p(x_-j,b)].
                   -.∫ p_0(x)∂(h_j(x_j)∂_jlog p(x))/∂ x_jx̣_j]x_-j
       =∑_j=1^m∫[∫ p_0(x)∂(h_j(x_j)∂_j log p(x))/∂ x_jx̣_j]x_-j.

Justified by the second half of (A2), by Fubini-Tonelli and linearity again

    B   =∑_j=1^m∫ p_0(x)∂(h_j(x_j)∂_j log p(x))/∂ x_jx,
       =∑_j=1^m∫ h'_j(x_j) ∂log p(x)/∂ x_jp_0(x)x+∑_j=1^m∫ h_j(x_j)∂^2log p(x)/∂ x_j^2p_0(x)x.

Thus,

    J_h(p)
    
    =    B+A+C
    
    =    ∫_ℝ_+^mp_0(x)∑_j=1^m[h_j'(x_j)∂log p(x)/∂ x_j+h_j(x_j)∂^2log p(x)/∂ x_j^2+1/2h_j(x_j)(∂log p(x)/∂ x_j)^2]x+C,

where C is a constant that does not depend on p.




 §.§ Proof of Theorems and Examples in Section <ref>



For 
exponential families and under the assumptions, the empirical loss
Ĵ_h(p_θ) in
(<ref>) becomes (up to an additive constant)

    Ĵ_h(p_θ)
    
    =    1/n∑_i=1^n∑_j=1^m[h_j'(X_j^(i))∂log p_θ(X^(i))/∂ X_j^(i)+h_j(X_j^(i))∂^2 log p_θ(X^(i))/∂ (X_j^(i))^2.
       .+
    	1/2h_j(X_j^(i))(∂log p_θ(X^(i))/∂ X_j^(i))^2]
    
    =    1/n∑_i=1^n∑_j=1^m[h_j'(X_j^(i))(θ^⊤t_j'(X^(i))+b_j'(X^(i)))+h_j(X_j^(i))(θ^⊤t_j”(X^(i))+b_j”(X^(i))).
                  +.1/2h_j(X_j^(i))(θ^⊤t'_j(X^(i))+b'_j(X^(i)))^2]
    
    =    1/n{1/2θ^⊤[∑_i=1^n∑_j=1^m
    h_j(X_j^(i))t_j'(X^(i))
    t_j'(X^(i))^⊤]θ+.
       .[∑_i=1^n∑_j=1^mh_j(X_j^(i))b_j'(X^(i))t_j'(X^(i))+h_j(X_j^(i))
    t_j”(X^(i))+h_j'(X_j^(i))t_j'
    (X^(i))]^⊤θ}+const,

which is quadratic in θ. Let

    Γ(𝐱)   =1/n∑_i=1^n∑_j=1^mh_j(X_j^(i))t_j'(X^(i))
    t_j'(X^(i))^⊤,
    g(𝐱)   =-1/n∑_i=1^n∑_j=1^m[h_j(X_j^(i))b_j'(X^(i))
    t'_j(X^(i))+h_j(X_j^(i))t_j”(X^(i))+
    h_j'(X_j^(i))t_j'(X_i)].

Then we can write Ĵ_h(p_θ)=1/2θ^⊤Γ(𝐱)θ-g(𝐱)^⊤θ+const.



By Theorem <ref>, Ĵ_h(p_θ)=1/2θ^⊤Γθ-g^⊤θ+const. The minimizer of Ĵ_h(p_θ) is thus available in the unique closed form θ̂≡Γ(𝐱)^-1g(𝐱) as long as Γ is invertible (C1).

Since Γ and g are sample averages, 
the weak law of large numbers yields that 

Γ→_p𝔼_p_0Γ≡Γ_0 and g→_p𝔼_p_0g≡g_0,

where existence of Γ_0 and g_0 is assumed in (C2). Since J_h(p_θ)=𝔼[Ĵ_h(p_θ)]=𝔼[1/2θ^⊤Γ(𝐱)θ-g(𝐱)^⊤θ]=1/2θ^⊤Γ_0θ-g_0θ and we know θ_0 minimizes J_h(p_θ) by definition, by the first-order condition we must have Γ_0θ_0=g_0.  Then by the Lindeberg-Lévy central limit theorem, 

    √(n)(g(𝐱)-Γ(𝐱)θ_0)→_d𝒩_m(0,Σ_0),

where Σ_0≡𝔼_p_0[(Γ(𝐱)θ_0-g(𝐱))(Γ(𝐱)θ_0-g(𝐱))^⊤], as long as Σ_0 exists (C2).
Thus, by Slutsky's theorem, 

    √(n)(θ̂-θ_0)≡√(n)(Γ(𝐱)^-1(g(𝐱)-Γ(𝐱)θ_0))→_d𝒩_r(0,Γ_0^-1ΣΓ_0^-1),

as long as Γ_0 is invertible (C2).



For the second half of the theorem, (C2) 𝔼_p_0Γ(𝐱)<∞ and 𝔼_p_0g(𝐱)<∞ implies 𝔼_p_0|Γ(𝐱)|<∞ and 𝔼_p_0|g(𝐱)|<∞, so by strong law of large numbers (and a union bound on at most k^2 null sets)

    Γ(𝐱)→_a.s.Γ_0,  g(𝐱)→_a.s.g_0.

Then outside a null set,

    θ̂≡Γ(𝐱)^-1g(𝐱)→_a.s.Γ_0^-1g_0=θ_0.





We choose to estimate θ≡μ/σ^2. Then by (<ref>) and (<ref>), 

    μ̂_h   =σ^2θ̂≡σ^2Γ(x)^-1g(x)
       =-σ^2[∑_i=1^nh(X_i)t'(X_i)^2]^-1[∑_i=1^nh(X_i)b'(X_i)t'(X_i)+h(X_i)t”(X_i)+h'(X_i)t'(X_i)]
       =-σ^2[∑_i=1^nh(X_i)]^-1[∑_i=1^n-h(X_i)X_i/σ^2+h'(X_i)].

By Theorem <ref>,

    √(n)(μ̂_h-μ_0)   →_d 𝒩(0,σ^4𝔼_0[-h(X)X-μ_0/σ^2+h'(X)]^2/𝔼_0^2[h(X)])
       ∼𝒩(0,𝔼_0[-h(X)(X-μ_0)+σ^2 h'(X)]^2/𝔼_0^2[h(X)]).

By integration by parts, (suppressing the dependence of p_μ_0 on μ_0)

    𝔼_0[h(X)h'(X)(X-μ_0)]
    
    =    ∫_0^∞h'(x)h(x)(x-μ_0)p(x)x̣
    
    =   ∫_0^∞h(x)(x-μ_0)p(x)ḥ(x)
    
    =    .h^2(x)(x-μ_0)p(x)|_0^∞-∫ h(x)ḥ(x)(x-μ_0)p(x)
    
    =    -∫ h^2(x)p(x)x̣-∫ h(x)h'(x)(x-μ_0)p(x)x̣+∫ h^2(x)(x-μ_0)^2/σ^2p(x)x̣,

where the last step follows from the assumptions lim_x↘ 0^+h(x)=0 and lim_x↗+∞h^2(x)(x-μ_0)p_μ_0(x)=0. So

    𝔼_0[h(X)h'(X)(X-μ_0)]   =𝔼[h^2(X)((X-μ_0)^2/σ^2-1)]/2.

The asymptotic variance is thus 

    𝔼_0[-h(X)(X-μ_0)+σ^2 h'(X)]^2/𝔼_0^2[h(X)]
    
    =     𝔼_0[h^2(X)(X-μ_0)^2-2σ^2h^2(X)((X-μ_0)^2/σ^2-1)/2+σ^4 h'^2(X)]/𝔼_0^2[h(X)]
    
    =     𝔼_0[σ^2 h^2(X)+σ^4 h'^2(X)]/𝔼_0^2[h(X)].

The Cramér-Rao lower bound follows from taking the second derivative of log p_μ_0 with respect to μ_0.



We estimate θ≡ 1/σ^2. By (<ref>) and (<ref>), 

    θ̂   ≡Γ(x)^-1g(x)
       =-[∑_i=1^nh(X_i)t'(X_i)^2]^-1[∑_i=1^nh(X_i)b'(X_i)t'(X_i)+h(X_i)t”(X_i)+h'(X_i)t'(X_i)]
       =[∑_i=1^nh(X_i)(X_i-μ)^2]^-1[∑_i=1^nh(X_i)+h'(X_i)(X_i-μ)].

By Theorem <ref>, √(n)(θ̂-θ)→_d 𝒩(0,ς^2), where

    ς^2   ≡𝔼_0[h(X)((X-μ)^2/σ_0^2-1)-h'(X)(X-μ)]^2/𝔼_0^2[h(X)(X-μ)^2]
       =
    1/𝔼_0^2[h(X)(X-μ)^2](𝔼_0[h^2(X)(X-μ)^4/σ_0^4-2h^2(X)(X-μ)^2/σ_0^2+h^2(X)
    
                     +h'^2(X)(X-μ)^2-2h(X)h'(X)(X-μ)^3/σ_0^2+2h(X)h'(X)(X-μ)).

By integration by parts, (suppressing the dependence of p_σ_0^2 on σ_0^2)

    𝔼_0[h(X)h'(X)(X-μ)^3]
    
    =    ∫_0^∞h'(x)h(x)(x-μ)^3p(x)x̣
    
    =   ∫_0^∞h(x)(x-μ)^3p(x)ḥ(x)
    
    =    .h^2(x)(x-μ)^3p(x)|_0^∞-∫ h(x)ḥ(x)(x-μ)^3p(x)
    
    =    -∫ h(x)h'(x)(x-μ)^3p(x)x̣-3∫ h^2(x)(x-μ)^2p(x)x̣+∫ h^2(x)(x-μ)^4/σ_0^2p(x)x̣,

where the last step follows from the assumptions lim_x↘ 0^+h(x)=0 and lim_x↗+∞h^2(x)(x-μ)^3p_σ_0^2(x)=0. Combining this with (<ref>) we get

    √(n)(θ̂-θ)→_d𝒩(0,ς^2)∼𝒩(0,2𝔼_0[h^2(X)(X-μ)^2/σ_0^2]+𝔼_0[h'^2(X-μ)^2]/𝔼_0^2[h(X)(X-μ)^2]),

and so by the delta method, for σ̂_k^2≡θ̂^-1,

    √(n)(σ̂_h^2-σ_0^2)→_d𝒩(0,2σ_0^6𝔼_0[h^2(X)(X-μ)^2]+σ_0^8𝔼_0[h'^2(X-μ)^2]/𝔼_0^2[h(X)(X-μ)^2]).

The Cramér-Rao lower bound follows from taking the second derivative of log p_σ_0^2 with respect to σ_0^2.






The proof of Theorem <ref> from <cit.> does not
rely on the fact that the original Γ is an
unbiased estimator for the population Γ_0, but
instead only requires one to bound
Γ-Γ_0_∞. Thus,
writing the new matrix as
Γ+diag(γ), by Theorem
<ref> it suffices to prove that for any τ>3, we can
bound
Γ(𝐱)+diag(γ(𝐱))-Γ_0_∞
by some ϵ_1 and
g(𝐱)-g_0_∞ by some
ϵ_2, uniformly with probability 1-1/m^τ-3.
Recall from (<ref>) that the j^th block of Γ+diag(γ)∈ℝ^m^2× m^2 has (k,ℓ)-th entry 

    1/n∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))+γ_kj·1_{k=ℓ}.

The entry in g∈ℝ^m^2 (obtained by linearizing a m× m matrix) corresponding to (j,k) is 

    1/n∑_i=1^nX^(i)_kh_j'(X^(i)_j)+1_{j=k}·1/n∑_i=1^nh_j(X^(i)_j).

Denote M≡max_jsup_x>0h_j(x) and M'≡max_jsup_x>0h_j'(x), and let c_X≡ 2max_j(2√(Σ_jj)+√(e) 𝔼_0 X_j). Using results for sub-Gaussian random variables from Lemma <ref>.2 in Appendix <ref>, we have for any t_1>0,

    (|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))-𝔼_0 X_kX_ℓh_j(X_j)|>t_1)   ≤2exp(-min(nt_1^2/2M^2c_X^4,nt_1/2Mc_X^2)).

Thus choosing ϵ_1≡ 2Mc_X^2max{2(log m^τ+log 12)/n,√(2(log m^τ+log 12)/n)}, for γ_kj≤ϵ_1/2, we have

    (|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))+γ_kj1_{k=ℓ}-𝔼_0 X_kX_ℓh_j(X_j)|>ϵ_1)
    ≤    (|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))-𝔼_0 X_kX_ℓh_j(X_j)|>ϵ_1/2)
    ≤    2exp(-min(nϵ_1^2/8M^2c_X^4,nϵ_1/4Mc_X^2))≤1/6m^τ.

On the other hand, for c_n,m≡max{2(log m^τ-2+log 6)/n,√(2(log m^τ-2+log 6)/n)}, by Lemma <ref>.2 again, for any k=1,…,m,

    (n^-1∑_i=1^n (X_k^(i))^2-𝔼_0X_k^2>c_n,mc_X^2)≤exp(-min(nc_n,m^2/2,nc_n,m/2))≤1/6m^τ-2.

By definition c_X^2=4max_k(4Σ_kk+4√(e)√(Σ_kk) 𝔼_0X_k+e(𝔼_0X_k)^2)≥ 4emax_k(Σ_kk+(𝔼_0X_k)^2). By Lemmas <ref>.2 and <ref>.1, (X_k)≤Σ_kk, so c_X^2≥ 4emax_k 𝔼_0 X_k^2. Thus, (<ref>) implies that 1/n(c_n,m+1/4e)∑_i=1^n (X_k^(i))^2≤c_n,mc_X^2+𝔼_0X_j^2/c_n,m+1/4e≤ c_X^2 with probability at least 1/(6m^τ-2). Then by our definition of ϵ_1 and by a union bound over m events,

    Mmax{2(log m^τ+log 12)/n,√(2(log m^τ+log 12)/n)}/1/4e+max{2(log m^τ-2+log 6)/n,√(2(log m^τ-2+log 6)/n)}n^-1∑_i=1^n(X_k^(i))^2≤ϵ_1/2

holds for all k=1,…,m, with probability ≥ 1/(6m^τ-3). Note that the max in the numerator of the previous display is always larger than the max in the denominator, and that 2(log m^τ-2+log 6)/n≥2(log m+log 6)/n red(or simply 2log 6/n to make γ smaller? too small!) by the assumption that τ>3. Since x/1/4e+x is increasing in x and h_j(X_j^(i))≤ M, we define

    γ_kj≡(1-1/1+4emax{2(log m+log 6)/n,√(2(log m+log 6)/n)})n^-1∑_i=1^n(X_k^(i))^2h_j(X_j^(i))

red(or simply √(2log 6/n)/1/4e+1X_k^2h_j(X_j) assuming n>2log 6≈ 3.58 to make γ smaller – try!!), which is just a constant multiple of the (k,k)-th entry of Γ_j itself, so that max_k,jγ_kj≤ϵ_1/2 with probability 1/(6m^τ-3). Then by a union bound for (<ref>) over m^3 events,

    (max_j,k,ℓ|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))+γ_kj1_{k=ℓ}-𝔼_0 X_kX_ℓh_j(X_j)|≥ϵ_1)≤1/3m^τ-3.

Now on the other hand, Lemma <ref>.1 and Hoeffding's inequality give for any t_2,1,t_2,2>0 that,

    (|n^-1∑_i=1^nX^(i)_kh_j'(X^(i)_j)-𝔼_0 X_kh_j'(X_j)|≥ t_2,1)   ≤2exp(-nt_2,1^2/2M'^2c_X^2),
    (|n^-1∑_i=1^nh_j(X^(i)_j)-𝔼_0h_j(X_j)|≥ t_2,2)   ≤2exp(-2nt_2,2^2/M^2).

Choosing ϵ_2,1≡√(2)M'c_X√(log m^τ-1+log 6/n), ϵ_2,2≡ M√(log m^τ-2+log 6/2n) and taking union bounds over m^2, and m events, respectively, we have

    (max_j,k|n^-1∑_i=1^nX^(i)_kh_j'(X^(i)_j)-𝔼_0 X_kh_j'(X_j)|≥ϵ_2,1)   ≤1/3m^τ-3,
    (max_j|n^-1∑_i=1^nh_j(X^(i)_j)-𝔼_0h_j(X_j)|≥ϵ_2,2)   ≤1/3m^τ-3.

Hence, by (<ref>) (<ref>) (<ref>), with probability at least 1-m^3-τ, Γ(𝐱)-Γ_0_∞<ϵ_1 and g(𝐱)-g_0_∞<ϵ_2≡ϵ_2,1+ϵ_2,2. Then by Theorem <ref>, for any τ>3, let

    c_2   ≡6/αc_Γ_0,         n≥max{2M^2c_X^4c_2^2d_𝐊_0^2(τlog m+log 12), 2Mc_X^2c_2d_𝐊_0(τlog m+log 12)},
    λ   >3(2-α)/αmax{c_𝐊_0ϵ_1,ϵ_2}
       ≡3(2-α)/αmax{4Mc_𝐊_0c_X^2(log m^τ+log 12)/n,.
           .2Mc_𝐊_0c_X^2√(2(log m^τ+log 12)/n),√(2)M'c_X√(log m^τ-1+log 6/n)+M√(log m^τ-2+log 6/2n)},

then dϵ_1≤α/(6c_Γ_0) and the results follow.





 §.§ Proof of Theorems in Section <ref>


 

Case b≠ 0: We use a strategy similar to that of <cit.>.  Let 𝒱_1={v:v_1=1,v∈ℝ_+^m}. Then by Fubini-Tonelli the normalizing constant is,

    =   ∫_ℝ_+^mexp(η^⊤x^b-1_m/b-1/2ax^a^⊤𝐊x^a)x
    
    =   ∫_𝒱_1∫_0^∞exp(η^⊤z^bv^b-1_m/b-1/2az^2av^a^⊤𝐊v^a)ẓv
    ∝   ∫_𝒱_1∫_0^∞exp(z^b(η^⊤v^b)/b-z^2a(v^a^⊤𝐊v^a)/(2a))ẓv.


Here 𝒱_1 is compact and the inner integral, if finite, is continuous in v. It thus suffices to show that the inner integral is finite at every single v∈𝒱_1.

Fixing v∈𝒱_1, write A≡ A(v)≡v^a^⊤𝐊v^a/(2a) and B≡ B(v)≡(η^⊤v^b)/b. We need to show that

    N(A,B,a,b)≡∫_0^∞exp(-Az^2a+Bz^b)ẓ<+∞.

Recall that (CC1) v^⊤𝐊v>0 for all v∈ℝ_+^m\{0}, so A>0.



















  * Suppose B≤ 0. If a>0, then N(A,B,a,b)≤∫_0^∞exp(-N_𝐊,1z^2a)ẓ, a finite constant since a>0.

  * Suppose B>0.  First note that if a>0 and b>0, exp(-Az^2a+Bz^b)≤exp(-N_𝐊,1z^2a+N_𝐊,2z^b). We then want to bound exp(-N_𝐊,1z^2a+N_𝐊,2z^b)≤ N_0exp(-N_𝐊,1z^2a/2) by some finite constant N_0>0, so that N(A,B,a,b)≤ N_0∫_0^∞exp(-N_𝐊,1z^2a/2)ẓ, a finite constant if a>0. Thus, it remains to give conditions so that exp(-N_𝐊,1z^2a/2+Bz^b) is bounded by some finite constant N_0, which by continuity only requires a finite limit as z↘ 0 and as z↗+∞. As z↗+∞, while B>0, since N_𝐊,1>0 and a>0 we have -N_𝐊,1z^2a/2↘ -∞. We thus need b<2a so that the sum of the two does not go to positive infinity. On the other hand, as z↘ 0, -N_𝐊,1z^2a/2↗ 0, so b>0, otherwise z^b↗+∞. In conclusion, we require that 2a>b>0.

It thus suffices to require (CC1) and that 2a>b>0 to eliminate restrictions on B, and hence on η. That is, η can take value in the entirety of ℝ^m.




  * Suppose B≤ 0. Then N(A,B,a,b)≤∫_0^∞exp(-Az^2a)ẓ=A^-a/2Γ(1+1/(2a)), a finite constant since A>0 and a>0.

  * Suppose B>0.  We first want  to bound exp(-Az^2a+Bz^b)≤ N_0exp(-Az^2a/2) by some finite constant N_0>0, so that N(A,B,a,b)≤ N_0∫_0^∞exp(-Az^2a/2)ẓ, a finite constant for a>0. Thus, it remains to give conditions so that exp(-Az^2a/2+Bz^b) is bounded by some finite constant N_0, which by continuity only requires a finite limit as z↘ 0 and as z↗+∞. As z↗+∞, Bz^b↗+∞, while -Az^2a/2↘ -∞. We thus need b<2a so that the sum of the two does not go to positive infinity. On the other hand, as z↘ 0, -Az^2a/2↗ 0, so we need b>0, otherwise z^b↗+∞. In conclusion, we require that 2a>b>0.

It thus suffices to require (CC1) and (CC2) 2a>b>0 to eliminate restrictions on B, and hence on η. That is, η can take value in the entirety of ℝ^m.



Case b=0: Again in (CC1) we assume v^⊤𝐊v>0 for all v∈ℝ_+^m\{0}. Since 𝒱_2≡{v:v_2=1,v∈ℝ_+^m} is compact and v^⊤𝐊v is continuous in v and strictly positive on 𝒱_2, the image of 𝒱_2 under v^⊤𝐊v is a compact subset of (0,∞), i.e. N_𝐊≡min_v∈ℝ_+^m\{0}v^⊤𝐊v/v^⊤v≡min_v∈𝒱_2v^⊤𝐊v>0. We thus have





    =   ∫_ℝ_+^mexp(η^⊤log(x)-1/2ax^a^⊤𝐊x^a)x
    ≤   ∫_ℝ_+^mexp(η^⊤log(x)-N_𝐊/2ax^a^⊤x^a)x
    
    =   ∏_j=1^m∫_0^∞exp(η_jlog(x_j)-N_𝐊/2ax_j^2a)x̣_j
    
    =   ∏_j=1^m[Γ(η_j+1/2a)(N_𝐊/2a)^-η_j+1/2a/2a],

where the integration follows by change of variable and requires a>0. Assuming a>0, the last quantity is finite if and only if η≻-1_m, by definition of the gamma function.


In conclusion, given conditions (CC1) min_v∈ℝ_+^m\{0}v^⊤𝐊v>0, (CC2) 2a>b>0, and (CC3) a>0, b=0 and η≻-1_m, the unnormalized density (<ref>) has a finite normalizing constant when (CC1) and (CC2) both hold, or (CC1) and (CC3) both hold.

The centered settings, where the term involving x^b is excluded, can be considered as a special case of both (1) and (2) with η≡0, and thus (CC1) and a>0 are sufficient.



Recall assumptions (A1) and (A2):


  * p_0(x)h_j(x_j)∂_jlog p(x)|^x_j↗+∞_x_j↘ 0^+=0.,  ∀x_-j∈ℝ_+^m-1,  ∀ p∈𝒫_+;

  * 𝔼_p_0∇log p(X)∘h^1/2(X)_2^2<+∞,    𝔼_p_0(∇log p(X)∘h(X))'_1<+∞,   ∀ p∈𝒫_+.

Let 𝐊_0 and η_0 be the true parameters so that p_0∈𝒫_+, with 𝒫_+ corresponding to a parameter space in which all parameters satisfy the conditions for a finite normalizing constant. We now give sufficient conditions for h to satisfy (A1) and (A2).


Conditions for (A1):
Fix j=1,…,m and x_-j∈ℝ_+^m-1. We show that the conditions on h_j imply that the limits go to 0 as x_j↗+∞ and as x_j↘ 0^+, which is stronger than (A1); in fact, from (<ref>) below, the limits cannot go to a nonzero finite constant assuming an h with polynomial tail, since a>0 and B_1≡κ_0,jj>0 for all j. Now,

    ∝= p_0(x)h_j(x_j)∂_jlog p(x)
       ∝ h_j(x_j)exp(-1/2ax^a^⊤𝐊_0x^a+η_0^⊤x^b-1_m/b)∂_j(-1/2ax^a^⊤𝐊x^a+η^⊤x^b-1_m/b)
       ∝ h_j(x_j)exp(-1/a(k_0,j,-j^⊤x^a_-j)x_j^a-κ_0,jj/2ax_j^2a+η_0,jx_j^b-1/b)×
       (-k_j,-j^⊤x^a_-jx_j^a-1-κ_jjx_j^2a-1+η_j x_j^b-1)
       ≡ h_j(x_j)exp(A_1x_j^a/a+B_1x_j^2a/2a+C_1x_j^b-1/b)(A_2x_j^a-1+B_2x_j^2a-1+C_2x_j^b-1),

where A_1≡-k_0,j,-j^⊤x^a_-j, A_2≡-k_j,-j^⊤x^a_-j, B_1≡-κ_0,jj<0 and B_2≡-κ_jj<0 by condition (CC1). Finally C_1≡η_0,j, C_2≡η_j.


  * Let x_j↗+∞. If b>0, since 2a>b>0 and B_1<0, the exponential term in (<ref>) decreases to 0 exponentially and its reciprocal dominates any polynomial functions. Thus, the entire product goes to 0 if h_j(x_j) grows no faster than polynomially as x_j↗+∞. If b=0, the C_1log x_j term is again dominated by B_1x_j^2a/(2a), and the same conclusion holds.

  * Let x_j↘ 0.


  * Let b>0. Then the exponential term in (<ref>) goes to constant exp(-C_1/b), and we only need 

    lim_x_j↘ 0^+h_j(x_j)(A_2x_j^a-1+B_2x_j^2a-1+C_2x_j^b-1)=0.



  * If a>1 and b>1, the second term in (<ref>) is a polynomial with three terms having powers ≥min{a-1,b-1}. The product goes to zero if and only if h_j(x_j)=o(x_j^max{1-a,1-b}) as x_j↘ 0. Note that this is satisfied by any h_j that has a finite right limit at 0.

  * If a=1 and b≥ 1, or a≥ 1 and b=1, then the second term in (<ref>) is a polynomial of non-negative power plus a potentially nonzero constant. A 
sufficient condition for (<ref>) is thus lim_x_j↘ 0h_j(x_j)=0.

  * If a<1 or b<1, then the second part in (<ref>) is a polynomial having terms with negative degree ≥min{a-1,b-1}. To counteract this a 
sufficient condition is h_j(x_j)=o(x_j^max{1-a,1-b}).

In conclusion, lim_x_j↘ 0^+p_0(x)h_j(x_j)∂_jlog p(x)=0 if and only if 
    lim_x_j↘ 0^+h_j(x_j)/x_j^max{1-a,1-b}=0.


  * Now assume b=0. Then, (<ref>) now becomes

    h_j(x_j)exp(A_1x_j^a/a+B_1x_j^2a/(2a)+C_1log x_j)(A_2x_j^a-1+B_2x_j^2a-1+C_2/x_j).

With C_1log x_j dominating, the exponential part scales as x_j^C_1. We thus require

    lim_x_j↘ 0^+h_j(x_j)(A_2x_j^a-1+C_1+B_2x_j^2a-1+C_1+C_2x_j^C_1-1)=0,

which by the previous discussion on (<ref>) holds if and only if

    lim_x_j↘ 0^+h_j(x_j)/x_j^1-C_1=0

since 1-a-C_1<1-C_1. 


In summary, (A1) is satisfied if h_j(x_j) grows at most polynomially as x_j↗ +∞, and lim_x_j↘ 0^+h_j(x_j)/x_j^max{1-a,1-b}=0 if b>0, or lim_x_j↘ 0^+h_j(x_j)/x_j^1-η_0,j=0 if b=0.


Conditions for (A2): For (A2), we consider powers of x as the h functions for simplicity; conclusions for other functions that have the same tail behavior (big-O scaling) as x↘0 and x↗+∞ follow similarly. Sufficiency results for piecewise power functions follow by partitioning, and similarly for other functions h whose function values and derivatives can be bounded by those of some piecewise power function (e.g. truncated powers), since (A2) is on integrability of products involving positive powers of h and h'.

Let 𝐊_0 and η_0 be the true parameters from the parameter space that satisfies the conditions for finite normalizing constant. By part (2) of the proof of Theorem <ref>, the assumption that 𝐊_0 satisfies (CC1) implies that min_v∈ℝ_+^m\{0}v^⊤𝐊_0v/v^⊤v≡ N_𝐊_0>0. Then we have the following decomposition

    p_𝐊_0,η_0(x)   ≡exp(-1/2ax^a^⊤𝐊_0x^a+η_0^⊤x^b-1_m/b)
       ≤∏_j=1^mexp(-N_𝐊_0/2ax^2a_j+η_0,jx_j^b-1/b).


Then for any other 𝐊 and η in the parameter space, for the first part of (A2) it suffices to show for any j=1,…,m that D<∞, where

    D   ≡∫_ℝ_+^m∏_j=1^mexp(-N_𝐊_0/2ax^2a_j+η_0,jx_j^b-1/b)h_j(x_j)×
       (-κ_jjx^2a-1_j-∑_i≠ jκ_jix_i^ax_j^a-1+η_j x_j^b-1)^2x
       ≥∫_ℝ_+^mp_0(x)h_j(x_j)(∂_jlog p(x))^2x.

Note that 

    =(-κ_jjx^2a-1_j-∑_i≠ jκ_jix_i^ax_j^a-1+η_j x_j^b-1)^2
       =κ_jj^2x_j^4a-2+∑_i≠ j,ℓ≠ jκ_jiκ_jℓx_i^ax_ℓ^ax_j^2a-2+η_j^2x_j^2b-2+2∑_i≠ jκ_jjκ_jix_i^ax_j^3a-2
       -2∑_i≠ jκ_jiη_j x_i^ax_j^a+b-2-2κ_jjη_jx_j^2a+b-2.

Thus, plugging this back in the definition of D, we can split D into a sum of six terms D_1 through D_6, each of which is a sum of terms of the form

    =   ∫_ℝ_+^m∏_k=1^mexp(-N_𝐊_0/2ax_k^2a+η_0,kx_k^b-1/b)h_j(x_j)x_i^pow_ix_ℓ^pow_ℓx_j^pow_jx
    
    =   ∏_k≠ j∫_0^∞exp(-N_𝐊_0/2ax_k^2a+η_0,kx_k^b-1/b)x_k^pow_kx̣_k
    
       ×∫_0^∞exp(-N_𝐊_0/2ax_j^2a+η_0,jx_j^b-1/b)h_j(x_j)x_j^pow_jx̣_j

times a constant involving 𝐊 and η_j, where pow_k≥ 0 for each k≠ j. We have thus decomposed the integral into a product of univariate integrals. Note that

    ∫_0^∞exp(-N_𝐊_0/2ax_i^2a+η_0,ix_i^b-1/b)x_i^pow_ix̣_i

is finite for all pow_i≥ 0 regardless of whether b is nonzero, since we assumed 𝐊_0 and η_0 to lie in the parameter space with a finite normalizing constant. Indeed, if b>0 then the terms in the exponential is a regular polynomial with positive degree and a negative leading term; if b=0 then integrability follows from η_0,i+pow_i≥η_0,i>-1. Thus, we only need to consider the univariate integral that involve the x_j terms, namely

    ∫_0^∞exp(-N_𝐊_0/2ax_j^2a+η_0,jx_j^b-1/b)h_j(x_j)x_j^pow_jx̣_j,

where pow_j takes value in {4a-2,2a-2,2b-2,3a-2,a+b-2,2a+b-2}⊆[2min{a,b}-2,4a-2]. 
We split the integral into two parts over [0,1] and [1,∞], respectively. 


  * If b>0, on [0,1] the exponential part is bounded above and below by positive constants, and for (A1) we require h_j(x)=o(x^1-min{a,b}) as x↘ 0^+, so the integrand is o(x^min{a,b}-1)=o(x^-1) and is thus integrable on [0,1]. The integrand on [1,∞) is integrable as  in (A1) we assume h to grow at most polynomially. 

  * If b=0, pow_j∈[-2,4a-2] and the integrand becomes 
    exp(-N_𝐊_0x_j^2a/(2a))h_j(x_j)x_j^pow_j+η_0,j.
 On [0,1], (A1) requires h_j(x)=o(x^1-min_jη_0,j), so h_j(x_j)x_j^pow_j+η_0,j=o(x^-1) and the integrand is again integrable. Integrability on [1,∞) follows similarly to the case with b>0.


Now consider the second part of (A2). By definition 𝔼_p_0(∇log p(X)∘h(X))'_1 equals

    ∫_ℝ_+^mp_0(x)∑_j=1^m|h_j'(X_j)∂_jlog p(X)+h_j(X_j)∂_j^2log p(X)|x
    ≤    ∑_j=1^m∫_ℝ_+^mexp(-N_𝐊_0/2ax^2a_j+η_0,jx_j^b-1/b)|h_j'(x_j)(-κ_jjx^2a-1_j-∑_i≠ jκ_jix_i^ax_j^a-1+η_j x_j^b-1).
       .+h_j(x_j)(-κ_jj(2a-1)x^2a-2_j-∑_i≠ jκ_ji(a-1)x_i^ax_j^a-2+(b-1)η_j x_j^b-2)|x.

By the triangle inequality and the fact that h_j≥ 0 and h_j'≥ 0, similar to the proof for the first part, for each j the integral can be bounded by a sum of six integrals, each of the form

    const×∫_ℝ_+^m∏_k=1^mexp(-N_𝐊_0/2ax_k^2a+η_0,kx_k^b-1/b)h_j(x_j)x_i^pow_ix_j^pow_jx,

or with h_j replaced by h_j'. Finiteness thus follows from the same type of discussion by noting that h_j(x)=o(x^1-min{a,b}) and h'_j(x)=o(x^-min{a,b}).


We conclude that if the true and the proposed parameters give densities with finite normalizing constants, and if h satisfies assumption (A1), then (A2) is automatically satisfied.

In the centered case where we assume η≡0, we only need lim_x_j↘ 0^+h_j(x_j)/x_j^1-a=0 as it is a special case with b=2a.




 §.§ Proof of Theorems in Section <ref>



By Theorem <ref>, under assumptions in that theorem, the support of Ψ̂ is a subset of the true support of Ψ_0, and Ψ̂-Ψ_0_∞≤c_Γ_0/2-αλ. Since Ψ_0 has |S_0| nonzero entries, 
    Ψ̂-Ψ_0_F=[∑_Ψ_0,jk≠ 0(Ψ̂_jk-Ψ_0,jk)^2]^1/2≤√(|S_0|)Ψ̂-Ψ_0_∞≤c_Γ_0/2-αλ√(|S_0|).

Similarly, by the definition of matrix ℓ_∞-ℓ_∞ norm, 

    Ψ̂-Ψ_0_2≤Ψ̂-Ψ_0_∞=max_j=1,…,m∑_k=1^m|Ψ̂_jk-Ψ_0,jk|≤c_Γ_0/2-αλ d_Ψ_0.

The result follows by also noting that Ψ̂-Ψ_0_2≤Ψ̂-Ψ_0_F.



The proof is based on Theorem <ref> and a probabilistic bound on Γ_γ-Γ_0_∞, where in the case of centered Gaussian Γ=diag(𝐱𝐱^⊤,…,𝐱𝐱^⊤). Denote Σ_0=𝐊_0^-1. In particular, given τ>2 we wish to show that for ϵ=80√(2)c_0max_j(Σ_0,jj), assuming c_0≡√((τlog m+log 4)/n)<1/√(2),

    ℙ(|n^-1∑_i=1^nX_j^(i)X_k^(i)+γ_1j1_{j=k}-𝔼X_jX_k|>ϵ)≤ m^2-τ,

and so the results follow from Theorem <ref>.

By Lemma 1 of <cit.>, since X_j/√(Σ_0,jj) is Gaussian with mean 0 and standard deviation 1,

    (|n^-1∑_i=1^nX_j^(i)X_k^(i)-𝔼X_jX_k|>t)≤ 4exp(-nt^2/3200max_j(Σ_0,jj)^2)

for t∈ (0,40max_j(Σ_0,jj)). Denote the event as ℰ_j,k(t). Note that X_j^2≤max_jΣ_0,jj=ϵ/(80√(2)c_0). Then letting t=ϵ/2 and conditioning on the complement of ℰ_j,j(ϵ/2), we have

    n^-1∑_i=1^n X_j^(i)^2≤𝔼X_j^2+ϵ/2≤ϵ/2(1+1/40√(2)c_0).

Thus, choosing γ_ℓ j=(δ-1)∑_i=1^n X_j^(i)^2/n for ℓ=1,…,m (Γ has m identical blocks) with 1<δ<1+(1+1/(40√(2)c_0))^-1, by the triangle inequality and a union bound we have

    ℙ(max_j,k|n^-1∑_i=1^nX_j^(i)X_k^(i)+γ_1j1_{j=k}-𝔼X_jX_k|>ϵ)≤(ℰ_j,k(ϵ/2))=m^2-τ.

Since τ>2,  it holds that 1+(1+1/(40√(2)c_0))^-1=2-(1+40√(2)c_0)^-1 is larger than 2-(1+80√(log m/n))^-1≡ C(n,m), so it is safe to choose any δ∈(1,C(n,m)). Thus by the requirement on ϵ, the theorem statement holds when n>max(c^*c_1^2d_𝐊^2,2)(τlog m+log 4) with c^*=12800max_j(Σ_0,jj)^2.



The proof of Theorem <ref> from <cit.> does not rely on the fact that the original Γ is an unbiased estimator for the population Γ_0, but instead only requires one to bound Γ-Γ_0_∞. Thus, for Γ_γ=Γ+diag(γ), by Theorem <ref> it suffices to prove that for any τ>3, we can bound Γ(𝐱)+diag(γ(𝐱))-Γ_0_∞ by some ϵ_1 and g(𝐱)-g_0_∞ by some ϵ_2, uniformly with probability 1-m^3-τ. 
Recall from (<ref>) that the j^th block of Γ_γ∈ℝ^m^2× m^2 has (k,ℓ)-th entry 

    n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))+γ_kj·1_{k=ℓ}.

The entry in g∈ℝ^m^2 (obtained by linearizing a m× m matrix) corresponding to (j,k) is 

    n^-1∑_i=1^nX^(i)_kh_j'(X^(i)_j)+n^-11_{j=k}∑_i=1^nh_j(X^(i)_j).

Denote M≡max_jsup_x>0h_j(x), M'≡max_jsup_x>0h_j'(x), and c_X≡ 2max_j(2√(Σ_jj)+√(e) 𝔼_0 X_j). Using results for sub-Gaussian random variables from Lemma <ref>.2 in Appendix <ref>, we have for any t_1>0,


    (|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))-𝔼_0 X_kX_ℓh_j(X_j)|>t_1)   ≤2exp(-min(nt_1^2/2M^2c_X^4,nt_1/2Mc_X^2)).


Thus, choosing ϵ_1≡ 2Mc_X^2c_n,m, where c_n,m≡max{2(log m^τ+log 6)/n,√(2(log m^τ+log 6)/n)}, for γ_kj≤ϵ_1/2, we have

    (|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))+γ_kj1_{k=ℓ}-𝔼_0 X_kX_ℓh_j(X_j)|>ϵ_1)
    ≤    (|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))-𝔼_0 X_kX_ℓh_j(X_j)|>ϵ_1/2)
    ≤    2exp(-min(nϵ_1^2/8M^2c_X^4,nϵ_1/4Mc_X^2))≤1/3m^τ.

Denote the event inside the probability in (<ref>) as ℰ_k,ℓ,j(ϵ_1/2).

By definition, 

    c_X^2=4max_k(4Σ_kk+4√(e)√(Σ_kk) 𝔼_0X_k+e(𝔼_0X_k)^2)≥ 4emax_k(Σ_kk+(𝔼_0X_k)^2).
 
By Lemmas <ref>.2 and <ref>.1 from Appendix <ref>, (X_k)≤Σ_kk, so c_X^2≥ 4emax_k 𝔼_0 X_k^2≥ 4e𝔼_0X_k^2 h_j(X_j)/M. Thus, setting ϵ_1=2Mc_X^2c_n,m, on the complement of ℰ_k,k,j(ϵ_1/2) we have

    n^-1∑_i=1^nX_k^(i)^2h_j(X_j^(i))≤𝔼_0X_k^2h_j(X_j)+ϵ_1/2≤ϵ_1/2(1+1/4ec_n,m).


Then 

    1/1+1/(4ec_n,m)1/n∑_i=1^nX_k^(i)^2h_j(X_j^(i))≤ϵ_1/2

on the complement of ℰ_k,k,j(ϵ_1/2), again with c_n,m≡max{2(log m^τ+log 6)/n,√(2(log m^τ+log 6)/n)}.
Note that the multiplier on the left of (<ref>) is increasing in c_n,m, and that 2(log m^τ+log 6)>6log m by the assumption that τ>3. 
Thus, if we let 

    γ_kj≡1/1+1/(4emax{6log m/n,√(6log m/n)})1/n∑_i=1^nX_k^(i)^2h_j(X_j^(i)),

which is just a constant multiple of the (k,k)-th entry of Γ_j itself, with the constant explicitly calculable and a function of p and n only, then for k=ℓ

    (|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))+γ_kj-𝔼_0 X_kX_ℓh_j(X_j)|≥ϵ_1)≤(ℰ_k,ℓ,j(ϵ_1/2))≤1/3m^τ.

Since this also holds for k≠ℓ without the γ_kj term, by a union bound over m^3 events,

    (max_j,k,ℓ|n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))+γ_kj1_{k=ℓ}-𝔼_0 X_kX_ℓh_j(X_j)|≥ϵ_1)≤1/3m^τ-3.

Now, on the other hand, Lemma <ref>.1 and Hoeffding's inequality give for any t_2,1,t_2,2>0 that

    (|n^-1∑_i=1^nX^(i)_kh_j'(X^(i)_j)-𝔼_0 X_kh_j'(X_j)|≥ t_2,1)   ≤2exp(-nt_2,1^2/2M'^2c_X^2),
    (|n^-1∑_i=1^nh_j(X^(i)_j)-𝔼_0h_j(X_j)|≥ t_2,2)   ≤2exp(-2nt_2,2^2/M^2).

Choosing ϵ_2,1≡√(2)M'c_X√(log m^τ-1+log 6/n), ϵ_2,2≡ M√(log m^τ-2+log 6/2n) and taking union bounds over m^2, and m events, respectively, we have

    (max_j,k|n^-1∑_i=1^nX^(i)_kh_j'(X^(i)_j)-𝔼_0 X_kh_j'(X_j)|≥ϵ_2,1)   ≤1/3m^τ-3,
    (max_j|n^-1∑_i=1^nh_j(X^(i)_j)-𝔼_0h_j(X_j)|≥ϵ_2,2)   ≤1/3m^τ-3.

Hence, by (<ref>) (<ref>) (<ref>), with probability at least 1-m^3-τ, Γ_γ(𝐱)-Γ_0_∞<ϵ_1 and g(𝐱)-g_0_∞<ϵ_2≡ϵ_2,1+ϵ_2,2. Consider any τ>3, and let

    c_2   ≡6/αc_Γ_0,
    
    
    n   >max{2M^2c_X^4c_2^2d_𝐊_0^2(τlog m+log 6), 2Mc_X^2c_2d_𝐊_0(τlog m+log 6)},
    λ   >3(2-α)/αmax{c_𝐊_0ϵ_1,ϵ_2}
       ≡3(2-α)/αmax{4Mc_𝐊_0c_X^2(log m^τ+log 6)/n,.
         .2Mc_𝐊_0c_X^2√(2(log m^τ+log 6)/n),√(2)M'c_X√(log m^τ-1+log 6/n)+M√(log m^τ-2+log 6/2n)}.

Then d_𝐊_0ϵ_1≤α/(6c_Γ_0) and the results follow from Theorem <ref>.



Similar to the proof of Theorem <ref>, by Theorem <ref> it suffices to prove that for any τ>3, we can bound Γ_γ(𝐱)-Γ_0_∞ by some ϵ_1 and g(𝐱)-g_0_∞ by some ϵ_2, uniformly with probability 1-m^3-τ. 
Recall that Γ∈ℝ^(m^2+m)× (m^2+m) is a rearrangement of Γ^(*), which is in turn formed by Γ_11∈ℝ^m^2× m^2, Γ_12∈ℝ^m^2× m, Γ_12^⊤ and Γ_22∈ℝ^m× m, all of which are block-diagonal with m blocks.

The j^th block of Γ_11∈ℝ^m^2× m^2 has (k,ℓ)-th entry 

    n^-1∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i)),

the k^th entry in the j^th block of Γ_12 is 

    -n^-1∑_i=1^nX_k^(i)h_j(X_j^(i)),

the j^th diagonal entry of Γ_22 is 

    n^-1∑_i=1^nh_j(X_j^(i)).

On the other hand, g∈ℝ^(m^2+m) is a rearrangement of g^(*)≡ [g_1^⊤,g_2^⊤]^⊤, where the entry in g_1∈ℝ^m^2 (obtained by linearizing a m× m matrix) corresponding to (j,k), is 

    n^-1∑_i=1^nX^(i)_kh_j'(X^(i)_j)+n^-11_{j=k}∑_i=1^nh_j(X_j^(i)),

while the j-th component of g_2∈ℝ^m is 

    -n^-1∑_i=1^nh_j'(X_j^(i)).

Recalling that the bounds in Lemma <ref> also hold when μ≠ 0, we may then use bounds similar to those in the proof of Theorem <ref>, and use union bounds to arrive at analogous consistency results, modulus different constants. The amplifiers γ can be incorporated analogously.


WRONG PROOF

Recall from the proof of <ref> that defining M≡max_jsup_x>0h_j(x),M'≡max_jsup_x>0h_j'(x), and c_X≡ 2max_j(2√(Σ_jj)+√(e) 𝔼_0 X_j), we have for any t_1>0,

    (|1/n∑_i=1^n X_k^(i)X_ℓ^(i)h_j(X_j^(i))-𝔼_0 X_kX_ℓh_j(X_j)|>t_1)   ≤2exp(-min(nt_1^2/2M^2c_X^4,nt_1/2Mc_X^2)).

On the other hand, by <ref>.<ref> and <ref> we know ℙ(|X_j-𝔼_0X_j|≥ t_2)≤ 2exp(-ϵ^2/2Σ_jj) for all ϵ>0, from which we obtain

    ℙ(∃ i: X_j^(i)>𝔼_0X_j+√(Σ_jj2log n)+t_2)   ≤∑_i=1^nℙ(|X_j^(i)-𝔼_0X_j|>√(Σ_jj2log n)+t_2)
       ≤ 2nexp(-t_2^2+Σ_jj2log n+2t_2√(Σ_jj2log n)/2Σ_jj)
       ≤ 2exp(-t_2^2/(2Σ_jj)).

On the other hand,

    𝔼_0X_j+√(Σ_jj2log n)≤max_j√(e) 𝔼_0X_j+2√(log n)√(Σ_jj)/√(2)≤c_X√(log n)/2√(2).

Combining the previous two equations we thus have

    ℙ(∃ i: X_j^(i)>c_X√(log n)/2√(2)+t_2)≤ 2exp(-t_2^2/2Σ_jj)≤ 2exp(-8t_2^2/c_X^2).

We wish to use a union bound on (<ref>) and (<ref>) to derive a tail bound for 1/n∑_i X_k^(i)X_ℓ^(i)(X_j^(i))^a for a=1,2,…. Notice that the M in (<ref>) would correspond to the c_X√(log n)/2√(2)+t_2 on the left-hand side of (<ref>), and we wish to match the right-hand sides so that they are comparable. Thus, choosing t_2=c_Xn^1/2/(2√(2)),

    (|1/n∑_i=1^n X_k^(i)X_ℓ^(i)(X_j^(i))^a-𝔼_0 X_kX_ℓ(X_j)^a|>t_1)
    ≤    (|1/n∑ X_k^(i)X_ℓ^(i)|)












§ AUXILIARY LEMMAS AND DEFINITIONS



In this appendix, to simplify notation, when it is clear from the context, the operator 𝔼 is defined as the expectation under the true distribution, unless otherwise noted.








































 

The sub-Gaussian (r=2) and sub-exponential (r=1) norms of a random variable are 

    X_ψ_r≡sup_q≥ 1q^-1/r(|X|^rq)^1/(rq)≡sup_q≥ 1q^-1/rX_rq.

If X_ψ_2<∞ we say X is sub-Gaussian; if X_ψ_1<∞ we call X sub-exponential.
For a zero-mean sub-Gaussian random variable X also define the sub-Gaussian parameter

    τ(X)=inf{τ≥ 0:𝔼exp(tX)≤exp(τ^2 t^2/2), ∀ t∈}.



The definition of sub-Gaussian norm here allows for a non-centered
variable and differs from the one in <cit.>, which uses
X_q.  Instead, it coincides with θ_2 in
<cit.>. The sub-Gaussian parameter is defined  as in <cit.> and the sub-exponential norm as in <cit.>.




  *  For any X and r=1,2, X-𝔼 X_ψ_r≤ 2X_ψ_r and X_ψ_r≤X-𝔼 X_ψ_r+|𝔼 X|, as long as the expectation and norms are finite. 

  *  <cit.> τ(X) is
  a norm on the space of all zero-mean sub-Gaussian variables; so
  τ(X+Y)≤τ(X)+τ(Y).
  
If X is zero-mean sub-Gaussian, then var(X)≤τ^2(X), X_ψ_2≤2τ(X)/√(e), τ(X)≤√(e)X_ψ_2.
If X_1,…,X_n are i.i.d. zero-mean sub-Gaussian, τ(n^-1∑_i=1^nX_i)≤ n^-1/2τ(X_i).


  *  If 
X_1 and X_2 are sub-Gaussian (not necessarily independent) with X_1_ψ_2≤ K_1 and X_2_ψ_2≤ K_2, then X_1X_2 is sub-exponential with X_1X_2_ψ_1≤ K_1K_2.

  *  <cit.> If X is
  zero-mean sub-Gaussian and q>0, then

    |X|^q≤ 2(q/e)^q/2τ^q(X).



  *  <cit.> If X_1,…,X_n are independent zero-mean, sub-Gaussian variables, then for any ϵ>0,

    (|X_1|≥ϵ)   ≤ 2exp(-ϵ^2/2τ^2(X_1)),
    (|n^-1∑_i=1^nX_i|>ϵ)   ≤ 2exp(-nϵ^2/2max_iτ^2(X_i)).



  *  <cit.> If X_1,…,X_n are independent zero-mean sub-exponential random variables with K≥max_iX_i_ψ_1, 

then for any ϵ>0,

    (|X_1|≥ϵ)   ≤ 2exp(-min(ϵ^2/8e^2K^2,ϵ/4eK)),
    (|n^-1∑_i=1^nX_i|≥ϵ)   ≤ 2exp(-min(nϵ^2/8e^2K^2,nϵ/4e K)).


  *  <cit.> If for X_i i.i.d. there exists some B>0 such that

    sup_q≥ 2(|X|^q/q!)^1/q≤ B/2

then for all ϵ>0,

    (|n^-1∑_i=1^n(X_i-𝔼 X_i)|≥ϵ)≤ 2exp(-min(nϵ^2/2B^2,nϵ/2B)).

 






  1)
  For r=1,2, by the triangle inequality, X-𝔼
  X_ψ_r≤X_ψ_r+𝔼
  X_ψ_r=X_ψ_r+|𝔼
  X|≤X_ψ_r+|X|≤ 2X_ψ_r, where in the last
  step we used the definition of ·_ψ_r with q=1 for
  r=1 and |X|≤(|X|^2)^1/2 with q=2 for r=2. On the
  other hand, X_ψ_r≤X-𝔼
  X_ψ_r+𝔼 X_ψ_r=X-𝔼
  X_ψ_r+|𝔼 X|.



  2)
  
  These follow from Theorems 1.2 and 1.3 and Lemmas 1.2 and 1.7 from
  <cit.>, and √(3.1)e^9/16/√(2)≈
  1.6467≤1.6487≈√(e).


  3)
  
  By Hölder's inequality (or Cauchy-Schwarz),

    X_1X_2_ψ_1   =sup_q≥ 1q^-1(|X_1X_2|^q)^1/q=sup_q≥ 1q^-1(|X_1^qX_2^q|)^1/q
       ≤sup_q≥ 1q^-1[(|X_1|^2q)^1/2(|X_2|^2q)^1/2]^1/q
       ≤sup_q≥ 1[q^-1/2(|X_1|^2q)^1/2q]sup_q≥ 1[q^-1/2(|X_2|^2q)^1/2q]
       =X_1_ψ_2X_2_ψ_2≤ K_1K_2.



4-6) These are Lemma 1.4 and Theorem 1.5 in <cit.>,
and a consequence of Corollary 5.17 in <cit.>.






7) By Theorem 2.10 of <cit.> wherein we let v≡ nB^2/2 and c≡ B/2, we have 

    (|n^-1∑_i=1^n(X_i-𝔼 X_i)|≥ϵ)≤ 2exp(-nϵ^2/B^2+Bϵ)

for all ϵ>0. (Theorem 2.10 gives an one-sided bound; bound for the other side is obtained by taking X_i=-X_i). The inequality follows by splitting into cases ϵ≤ B and ϵ>B.




Suppose X follows a truncated normal distribution on ℝ_+^m with parameters μ and Σ=𝐊^-1≻0. Let X^(1),…,X^(n) be i.i.d. copies of X, with j-th component of the i-th copy being X^(i)_j. Then 


  * For j=1,…,p, τ(X_j-𝔼 X_j)≤√(Σ_jj). That is, the sub-Gaussian parameter of any marginal distribution of X, after centering, is bounded by the square root of its corresponding diagonal entry in the covariance parameter Σ. Then, for any ϵ>0,

    (|n^-1∑_i=1^nX^(i)_j-𝔼 X_j|>ϵ)≤2exp(-nϵ^2/2Σ_jj).

In particular, if h_0 is a function bounded by M_0, then for any ϵ>0,

    (|n^-1∑_i=1^nX^(i)_jh_0(X^(i)_k)-𝔼 X_jh_0(X_k)|≥ϵ)   ≤2exp(-nϵ^2/8M_0^2(2√(Σ_jj)+√(e) 𝔼 X_j)^2),
    τ(n^-1∑_i=1^n X^(i)_jh_0(X^(i)_k)-𝔼 X_jh_0(X_k))   ≤2M_0/√(n)(2√(Σ_jj)+√(e) 𝔼 X_j),
    n^-1∑_i=1^n X^(i)_jh_0(X_k^(i))-𝔼 X_jh_0(X_k)_ψ_2   ≤4M_0/√(en)(2√(Σ_jj)+√(e) 𝔼 X_j).


  * For j,k,ℓ∈{1,…,p}, if h_0 is a function bounded by M_0, then 

    X_jX_kh_0(X_ℓ)-𝔼 X_jX_kh_0(X_ℓ)_ψ_1≤M_0/2ec_X^2,
 
where c_X≡ 2max_j(2√(Σ_jj)+√(e) 𝔼X_j). In particular, for any ϵ>0,

    (|n^-1∑_i=1^n X^(i)_jX^(i)_kh_0(X_ℓ^(i))-𝔼 X_jX_kh_0(X_ℓ)|>ϵ)
    ≤2exp(-min(nϵ^2/2M_0^2c_X^4,nϵ/2M_0c_X^2)).







  1. Without loss of generality choose j=1. By the definition of sub-Gaussian parameters, we need to show that for all t∈,

    𝔼exp(t X_1)≤exp(t^2Σ_11/2+t𝔼 X_1),

which is equivalent to

    t^2Σ_11/2+t𝔼 X_1-log𝔼exp(tX_1)≥ 0  ∀ t∈.

Since the left-hand side of (<ref>) equals 0 at t=0, it suffices to show that its derivative,

    tΣ_11+𝔼 X_1-ḷọg̣𝔼exp(t X_1)/ṭ=tΣ_11+𝔼 X_1-Eexp(t X_1)/ṭ/𝔼exp(t X_1),

is non-negative on (0,∞) and non-positive on (-∞,0). By properties of moment-generating functions, /ṭ𝔼exp(t X_1) evaluated at t=0 equals 𝔼 X_1, so (<ref>) equals 0 at t=0. It in turn suffices to show the derivative of (<ref>), namely

    Σ_11-^̣2log𝔼exp(tX_1)/ṭ^2

is non-negative in t∈.

Given any vector v∈ℝ^p, define ℝ^p_+-v≡{u-v:u∈ℝ^p_+}. By <cit.>, denoting the first column of Σ as Σ_1, the moment-generating function of the marginal distribution of X_1 is 

    ∫_ℝ_+^p-(μ+tΣ_1)exp(-1/2x^⊤Σ^-1x)x/∫_ℝ_+^p-μexp(-1/2x^⊤Σ^-1x)xexp(tμ_1+1/2t^2Σ_11).

(<ref>) thus becomes

    -^̣2/ṭ^2log∫_ℝ_+^p-(μ+tΣ_1)exp(-1/2x^⊤Σ^-1x)x.

Showing this is non-negative in t∈ is equivalent to showing that the integral itself is log-concave in t. But

    ∫_ℝ_+^p-(μ+tΣ_1)exp(-1/2x^⊤Σ^-1x)x=
    ∫_ℝ^pexp(-1/2x^⊤Σ^-1x)1_ℝ_+^p-μ(x+tΣ_1)x

with exp(-1/2x^⊤Σ^-1x) log-concave in x and 1_ℝ_+^p-μ(x+tΣ_1) log-concave in (x,t) since ℝ_+^p-μ is a convex set. Since log-concavity is closed under multiplication and integration over ℝ^p, the integral is indeed log-concave, and our proof of the bound on the sub-Gaussian parameter of X_j-𝔼X_j is complete. The tail bound follows from <ref> of Lemma <ref>.

Now by <ref> and <ref> of Lemma <ref>, 


    X_j_ψ_2≤ 2√(Σ_jj/e)+𝔼 X_j.
 

If h_0 is a function bounded by M_0, then by definition 


    X_jh_0(X_k)_ψ_2≤ M_0(2√(Σ_jj/e)+𝔼 X_j).


By <ref> and <ref> of Lemma <ref> again, 


    τ(X_jh_0(X_k)-𝔼 X_jh_0(X_k))   ≤√(e)X_jh_0(X_k)-𝔼 X_jh_0(X_k)_ψ_2
       ≤ 2√(e)X_jh_0(X_k)_ψ_2
       ≤ 2M_0(2√(Σ_jj)+√(e) 𝔼 X_j).

The tail bound thus follows from the first inequality using <ref> of Lemma <ref>. By <ref> of the Lemma <ref>,

    τ(n^-1∑_i=1^n X^(i)_jh_0(X^(i)_k)-𝔼 X_jh_0(X_k))   ≤2M_0/√(n)(2√(Σ_jj)+√(e) 𝔼 X_j),
    n^-1∑_i=1^n X^(i)_jh_0(X^(i)_k)-𝔼 X_jh_0(X_k)_ψ_2   ≤4M_0/√(en)(2√(Σ_jj)+√(e) 𝔼 X_j).




2. By the proof of 1) of this lemma, X_j_ψ_2≤ 2√(Σ_jj/e)+𝔼 X_j, and by <ref> of Lemma <ref>, 

    X_jX_k_ψ_1≤(2√(Σ_jj/e)+𝔼 X_j)(2√(Σ_kk/e)+𝔼 X_k)≤max_j(2√(Σ_jj/e)+𝔼 X_j)^2.
 Since h_0 is a function bounded by M_0, by definition 

    X_jX_k h_0(X_ℓ)_ψ_1≤ M_0
    max_j(2√(Σ_jj/e)+𝔼 X_j)^2.
 Then by <ref> of Lemma <ref> again, 
    X_jX_kh_0(X_ℓ)-𝔼 X_jX_kh_0(X_ℓ)_ψ_1≤ 2M_0max_j(2√(Σ_jj/e)+𝔼 X_j)^2.
 The tail bound then follows from <ref> of Lemma <ref>.

Although not used for our consistency results, in the special case of h_0≡ 1, we also have the following lemma. The notable difference between bounds (<ref>) below and (<ref>) from Lemma <ref>.2 is in the constants and dependency on 𝔼X_j: The constants in the denominator in the right-hand side of  (<ref>) is smaller and thus gives a tighter bound, but (<ref>) is preferred when 𝔼X_j is notably large compared to √(Σ_jj), since the constant is only linear in 𝔼X_j.

Consider the setting in Lemma <ref>. Then for j,k∈{1,…,p}, for any ϵ>0,

    (n^-1|∑_i=1^n X^(i)_jX^(i)_k-𝔼 X_jX_k|≥ϵ)≤ 4exp(-min(2nϵ^2/C_1^2,nϵ/C_1)),

where C_1≡ 91max_jΣ_jj+72max_j𝔼 X_jmax_j√(Σ_jj).



We use a proof similar to Lemma 1 in <cit.> (note that 
𝔼X_j may be nonzero in our case). Define

    U^(i)_jk≡ X^(i)_j+X^(i)_k,   U_jk≡ X_j+X_k,   V^(i)_jk≡ X^(i)_j-X^(i)_k,   V_jk≡ X_j-X_k.

Since X^(i)_jX^(i)_k=1/4(U^(i)_jk^2-V^(i)_jk^2), by union bound we have

    (|n^-1∑_i=1^nX^(i)_jX^(i)_k-𝔼 X_jX_k|≥ϵ)
    ≤    (|n^-1∑_i=1^nU^(i)_jk^2-𝔼 U_jk^2|≥ 2ϵ)+(|n^-1∑_i=1^nV^(i)_jk^2-𝔼 V_jk^2|≥ 2ϵ).


We next define

Z^(i)_jk   
≡U^(i)_jk^2-𝔼 U_jk^2=A^(i)_jk^2+B_jk^(i)+C_jk,        X^(i)_j≡X^(i)_j-𝔼 X_j,

A^(i)_jk   ≡X^(i)_j+X^(i)_k,   B^(i)_jk≡2(𝔼 X_j+𝔼 X_k)(X^(i)_j+X^(i)_k),         C_jk≡-(X^(i)_j+X^(i)_k)^2.


Then since τ is a norm by <ref> of Lemma <ref>, A_jk is sub-Gaussian with parameter ≤√(Σ_jj)+√(Σ_kk), and B_jk is sub-Gaussian with parameter ≤2(𝔼 X_j+𝔼 X_k)(√(Σ_jj)+√(Σ_kk)). Using <ref> of Lemma <ref> together with  the inequality (a+b+c)^q≤ (3max{a,b,c})^q≤ 3^q(a^q+b^q+c^q) for all a,b,c≥ 0 and q>0, we have for any q≥ 2

    (|Z_jk|^q)^1/q   ≤(3^q(𝔼 |A_jk|^2q+𝔼 |B_jk|^q+|C_jk|^q))^1/q
       ≤ 3^1+1/q((|A_jk|^2q)^1/q + (|B_jk|^q)^1/q+|C_jk|)
       ≤ 3^1+1/q(2^1/q(2q/e)(√(Σ_jj)+√(Σ_kk))^2.
                  .+2^1/q√(q/e)2(𝔼 X_j+𝔼 X_k)(√(Σ_jj)+√(Σ_kk))+(X_j+X_k)).


Using (X+Y)≤2((X)+(Y)) and the fact that
(X_j)=(X_j-𝔼X_j)≤τ^2(X_j-𝔼X_j)≤Σ_jj
(by <ref> of Lemma
<ref> and 1) of Lemma <ref>, we then have

    (|Z_jk|^q/q!)^1/q
    ≤3^1+1/q2^3+1/q(q/e)max_jΣ_jj+2^3+1/q√(q/e)max_j𝔼 X_j·max_j√(Σ_jj)+4max_jΣ_jj/(q!)^1/q.

Since all three coefficients involving q are decreasing in q≥ 2, we have

    sup_q≥ 2(|Z_jk|^q/q!)^1/q≤(48√(3)/e+6√(6))max_jΣ_jj+24√(6/e)max_j𝔼 X_jmax_j√(Σ_jj).

Thus by <ref> of Lemma <ref>, letting B≡(91max_jΣ_jj+72max_j𝔼 X_jmax_j√(Σ_jj)), we have for all ϵ>0:

    (n^-1|∑_i=1^n Z_jk^(i)|≥ 2ϵ)≤ 2exp(-min(2nϵ^2/B^2,nϵ/B)).

A tail bound for the sample average of V_jk^2 can be similarly derived, and the result follows.







§ SIMULATION RESULTS FOR ERDÖS-RÉNYI GRAPHS


We revisit the simulations from Section <ref> but use Erdös-Rényi (ER) graphs in which
each possible edge is independently included with probability π.  Independent
uniform draws from [0.5,1] are used to fill the
non-zero off-diagonal entries
 of the symmetric matrix 𝐊_0. The
  diagonal elements are set such
that 𝐊_0 has minimum eigenvalue 0.1.  We choose π=0.08 for n=1000, and π=0.02 for n=80.



 §.§ Truncated GGMs

In this section we present the results for truncated GGMs.


  §.§.§ Choice of h


The results for truncated centered GGMs are reported in
Table <ref> and Figure <ref>.  Those for truncated
non-centered GGMs using the profiled estimator are in Table <ref> and
Figure <ref>.


























  §.§.§ Choice of multiplier

The results for truncated centered GGMs where each curve represents a different multiplier are shown in Figure <ref>, and those for truncated non-centered GGMs are in Figure <ref>, where each curve corresponds to a different ratio λ_𝐊/λ_η.














 §.§ Other a/b Models



Figure <ref> exhibits the results for the exponential models.



Figure <ref> displays the results for the gamma models.



Figures <ref> and <ref> demonstrate the results for a=3/2 and b=1/2 or b=0, respectively.













