Non-asymptotic error controlled sparse high dimensional precision matrix estimation
Adam B Kashlak


Fermi surface topology and large magnetoresistance in the topological semimetal candidate PrBi
    Amit Vashist ^1,*, R. K. Gopal^1,*, Divya Srivastava^2,3, M. Karppinen^2, and Yogesh Singh ^1
    December 30, 2023
=================================================================================================




  Estimation of a high dimensional precision matrix 
  is a critical problem to many areas of statistics
  including Gaussian graphical models and inference on 
  high dimensional data.
  Working under the structural assumption of sparsity, we propose 
  a novel methodology for estimating such matrices 
  while controlling the false positive rate–percentage of 
  matrix entries incorrectly chosen to be non-zero.
  We specifically focus on false positive rates 
  tending towards zero with finite sample guarantees.
  This methodology is distribution free, but
  is particularly applicable to the problem of Gaussian network 
  recovery.  We also consider applications to 
  constructing gene networks in genomics data.






Biometrika

2019
00
000


Advance Access publication on ...


1 1 1999
1 1 1999


A. B. KashlakNon-asymptotic error controlled estimation





 Department of Mathematical & Statistical Sciences,

University of Alberta, Edmonton, T6G 2G1 Canadakashlak@ualberta.ca


  Non-asymptotic error controlled sparse high dimensional 
  precision matrix estimation

    A. B. Kashlak
    December 30, 2023
==========================================================================================




  Estimation of a high dimensional precision matrix 
  is a critical problem to many areas of statistics
  including Gaussian graphical models and inference on 
  high dimensional data.
  Working under the structural assumption of sparsity, we propose 
  a novel methodology for estimating such matrices 
  while controlling the false positive rate, percentage of 
  matrix entries incorrectly chosen to be non-zero.
  We specifically focus on false positive rates 
  tending towards zero with finite sample guarantees.
  This methodology is distribution free, but
  is particularly applicable to the problem of Gaussian network 
  recovery.  We also consider applications to 
  constructing gene networks in genomics data.



  debiased estimator; Gaussian graphical models;
  random matrix theory; thresholding









§ INTRODUCTION


Attempting to estimate the graphical structure of a high dimensional 
network is problemsome when edges are rare but the number of nodes
is large.  
In standard statistical classification problems, 
we fix a palatable false positive
rate and aim to recover as many true positives as possible.
Thus, we treat support recovery as a binary classification
problem.
For p nodes, the number of potential edges to consider is on 
the order of p^2.  With a sample size n<p, there are too 
many parameters to accurately estimate.
However, when trying to classify edges as
significant or not under the assumption of sparsity, 
the assumption that most edges are not significant,  
if we have even a small false positive rate, then
this will result in many erroneous connections potentially obscuring 
follow up research.
For example, a genomics study considering the conditional correlation
structure of, say, 2000 genes will have to consider over two million 
potential edges.  A standard 1% false positive rate will result in
tens of thousands of erroneous connections.

To address this issue, we consider the extreme estimation setting
of recovering the support of a precision matrix Ω in the high
dimensional setting, p≫ n, assuming Ω belongs to a 
class of sparse positive definite matrices, and for false positive
rates α→0^+.  Making use of debiased 
estimators <cit.>,
the non-asymptotic results of <cit.>, 
and a clever subsampling methodology, we can achieve 
finite sample guarantees in this extreme setting.

There has been much past work on covariance and precision
matrix estimation much of which is summarized in the 
survey <cit.>.  Most estimators for high dimensional
precision matrices are based on ℓ^1 penalization 
including the graphical lasso <cit.>,
the CLIME and ACLIME estimators 
<cit.>, and 
the debiased estimator of <cit.> used for 
constructing confidence sets and running hypothesis tests.
These articles all rely on high dimensional asymptotics, 
which is that estimation is successful in the limit as 
p and n grow to infinity together generally such that
(log p) /n, or some variant thereof, is o(1).
Our work differs as it considers guaranteed results for
controlling the support recovery of Ω as the 
false positive rate α→0 for
fixed finite p≫ n, which is asymptotic in α, 
a controllable tuning parameter, instead of in n and p,
which are generally fixed in the real world experimental
setting.

Our method can be deconstructed into two steps.  The first is to 
find a bias-corrected initial estimate Ω̂_0 for the 
precision matrix Ω.  The second is to construct 
a ball in the operator norm topology corresponding to the target 
false positive rate and then
search this ball for a sparse estimator.
In actuality, given a false positive rate of α,
we step down in a binary fashion constructing successive
false positive rates of γ^-1, γ^-2,… 
and a sequence of estimators tending towards one with the
desired α.
Depending on the sample size, we can 
randomly partition our data set into smaller sets,
apply this two-step procedure to each subsample in parallel,
and combine the results for improved performance.
In step two, we use a binary search procedure which 
converges rapidly minimizing the computational burden. 

Our search methodology is effectively a method for controlled
thesholding of the initial debiased estimator.  The literature
on thresholding for precision matrices is quite light when 
compared to thresholding covariance matrices 
<cit.>.  The challenge is that unlike
the unbiased empirical covariance estimator, there is no 
unbiased empirical precision estimator.  Debiasing such 
precision matrix estimators allows for thresholding as 
mentioned in <cit.>.  Furthermore, 
<cit.> proposed an entrywise thresholding method
given some false positive rate α, which while 
working well in some high dimensional asymptotic sense,
proved to not achieve the target false positive rate
on simulated data. 




§ METHODOLOGY




 §.§ Notation



In this article, we primarily make use of the family of 
Schatten norms.  For Ω∈^p× p with 
i,jth entry Ω_i,j and eigenvalues _1,…,_p, 
we define the trace norm 
Ω_1 = ∑_i=1^p _i, 
the Hilbert-Schmidt norm
Ω_2 = (∑_i=1^p _i^2)^1/2 = 
  (∑_i,j=1^p Ω_i,j^2)^1/2, 
and the operator norm 
Ω_∞ = max_i=1,…,p_i.
We will also use the standard ℓ^p norms applied to vectors
in ^p
and entrywise to matrices in ^p× p.

The algorithm in Section <ref> makes use of 
hard thresholding.  For Ω∈^p× p, we
define the hard thresholding operator 
φ:^p× p×^+ →^p× p,
which returns a matrix with i,jth entry

    φ(Ω; t)_i,j = 
      {[           0 if Ω_i,j< t;       Ω_i,j if Ω_i,j≥ t ].

which simply removes entries from Ω with magnitude less 
than t.



 §.§ Initial Estimator



Let X_1,…,X_n∈^p be an collection of mean zero
random vectors with common unknown positive definite 
covariance matrix Σ and corresponding precision matrix 
Ω = Σ^-1.  Of course, we require the covariance
and its inverse 
to exist, but put no further assumptions on the distribution 
of the X_i at the moment except to assume that Ω is sparse.  
Specifically, 

    Ω∈𝒰(κ,δ) = { 
        A ∈^p× p : .
    . 
        A=A, ∑_i=1^p 1(A_i,j0)≤κ for j=1,…,p,
         if A_i,j0 then A_i,j≥δ}

where 𝒰(κ,δ) is the class of sparse 
^p× p
matrices with no more than κ entries in each row or column
non-zero and with non-zero entries bounded away from zero
allowing them to be detectable.  
In practice, we attempt to recover the support of the normalized 
matrix Ω with diagonal entries of 1 to avoid scale issues.

Our approach is similar to <cit.> who attempt to 
recover the support of a covariance matrix Σ by starting 
with an initial estimator, constructing a confidence ball, and
searching this ball for a sparser estimator.  However,
whereas this preceding work can use the empirical estimator,
Σ̂=n^-1∑_i=1^nX_iX_i, as an unbiased
initial estimator for controlled thresholding, 
we cannot construct an unbiased estimator
for the precision matrix Ω in the p>n setting. 

A standard estimator for sparse precision matrices is the
graphical lasso <cit.>, which is based on an ℓ^1
penalized maximum likelihood under the Gaussian distribution:

    Ω̂^GL = 
      Θ∈^p× p{Σ̂Θ - log(Θ)
        + ∑_i=1^p∑_j iΘ_i,j}

for some tuning parameter >0.  The graphical lasso is 
debiased in <cit.>, which uses a correction
factor based on the subgradient of the above optimization
extending the work of <cit.> that debiases 
the classic lasso estimator for linear regression.  
The resulting debiased estimator is

    Ω̂^deGL = 
      2Ω̂^GL - 
      Ω̂^GLΣ̂Ω̂^GL.


Alternatively, the CLIME method <cit.> solves
a constrained ℓ^1 optimization problem to find Ω̂^CL.

    min_Θ∈^p× pΘ_ℓ^1 such that max_i,jΣ̂Θ - I_p _i,j≤

with a more sophisticated version ACLIME <cit.> adapting
to individual entries.

A different regularized estimator is the ridge estimator, 



    Ω̂^RD =
      Θ∈^p× p(
        Σ̂Θ - I _2^2 + Θ_2^2
      ),

whose closed form solution is 
Ω̂^RD=(Σ+ I_p)^-1.
Consider the singular value decomposition X = UDV 
where U∈^n× r, D∈^r× r, and
V∈^p× r for r=min( n,p ).
In <cit.>, a bias corrected estimator for 
ridge regression is proposed.  This is achieved by finding
some other estimator, which is used to correct for the
projection bias.  We can apply the same method to the ridge 
precision matrix estimator to get

    Ω̂^deRD = 
      Ω̂^RD + P_0Ω̃

where Ω̃ is some other estimate for Ω
and P_0 = VV - diag(VV) where
VV is
the projection in ^p onto an n dimensional subspace
spanned by the columns of X, and P_0 is that projection 
with the diagonal entries set to zero.  
<cit.> uses the graphical lasso 
estimator for
Ω̃.








 §.§ Binary Search



Given an initial estimator from the previous setting
henceforth denoted Ω̂_0
for simplicity of notation and some γ>1 such that
α = γ^-m for some positive integer m corresponding
to the number of iterates of the below algorithm, 
we aim to construct 
an estimator
Ω̂_m with false positive rate α=γ^-m 
by carefully thresholding the 
entries in Ω̂_0 so that 

    *{
          (i,j)  :  
          (Ω̂_m)_i,j0 and 
          (Ω)_i,j=0, i j 
        }/*{
          (i,j)  :  (Ω)_i,j=0, i j
        }≈α,

which is that the desired false positive rate is achieved.
The two extreme estimators are the initial estimator 
Ω̂_0 and the diagonal matrix 
Ω̂_∞ with diagonal entries coinciding with 
Ω̂_0 and off-diagonal entries set to zero.
These correspond to the 100% and 0% false positive 
cases, respectively.  For the remainder, we normalize
Ω̂_0 to have unit diagonal thus making 
Ω̂_∞=I_p.

For values of α tending towards zero, we consider the 
operator norm balls centred at I_p being
B_α = {Π∈^p× p : Π - I_p _∞≤Ω̂_m - I_p _∞} 
where ·_∞ refers to the operator norm 
for ℓ^2(^p)→ℓ^2(^p), which is
the principal eigenvalue.  This motivates the following 
algorithm:


Constructing an error controlled estimator 

      Begin with an initial estimator from Section <ref>
         denoted Ω̂_0.
 
      Given the sth iterated estimator Ω̂_s,
    we construct the (s+1)th estimator

            Compute r_s = Ω̂_s-I_p_∞.

            Compute r_s'= r_s γ^-1/2.

            Find t_s+1=min(t) such that 
          φ(Ω̂_s;t)-I_p_∞≤ r_s'

             
          where φ is the hard threshold operator.

          
          
          
          
            Set 
      Ω̂_s+1 = φ(Ω̂_s;t_s+1)

      Repeat step 2 with Ω̂_s+1 until
                  s=m




  
  
  * Begin with an initial estimator from Section <ref>
    denoted Ω̂_2^0. 
    
    
    
  
  * Given the sth iterated estimator Ω̂_2^-s,
    we construct the (s+1)th estimator
  
    
  * Compute r_s = Ω̂_2^-s-I_p_∞.
    
  * Compute r_s'= r_s 2^-1/2.
    
  * Apply hard thresholding to Ω̂_2^-s
          to find the densest Ω̂_2^-s-1
          such that 
          Ω̂_2^-s-1-I_p_∞≤ r_s'.
  
  
  * Repeat step 2 with Ω̂_2^-s-1.
  
  * Stop once 2^-s-1 is small enough.



In this algorithm, we quickly locate the densest estimator close 
to I_p, alternatively being the sparsest estimator close to 
Ω̂_0, by using a binary search procedure.
Given the sth iterated matrix as a starting point,
we set the smallest half of the non-zero entries in magnitude in 
Ω̂_s to zero, compute the distance to I_p
and then if the distance is greater than r_s', we remove half 
of the remaining entries whereas if the distance is less than r_s',
we reintroduce half of the removed entries.


  
  For choice of γ>1 and number of steps m, values of 
  γ closer to 1 result in smaller steps adding stability
  but requiring a larger m to achieve a small false positive rate.
  In the simulations of Section <ref>, we use γ=2. 
  
  For the choice of initial estimator from Section <ref>, 
  the best performance in
  simulated data experiments on multivariate normal data
  was achieved by
  beginning the above procedure with the debiased
  graphical lasso estimator of <cit.>. 
  This is mainly because their estimator follows the requirements
  of the theorems in the subsequent section as long as the 
  data under analysis is sub-Gaussian.
  However, good performance is still observed when the 
  data is sub-exponential as can be seen in the supplementary
  material.



  A similar search algorithm is proposed in <cit.>.
  The main differences are that in the cited work, the operator
  norm balls are confidence sets centred about the empirical
  covariance estimator rather than centred about the identity
  matrix I_p.  Furthermore, the radii are reduced by a factor
  of γ^-1 in that previous work due to the low rank structure 
  of the empirical covariance estimator.  Specifically, 
  Σ̂ = O(p) rather than O(p^1/2).




 §.§ Theoretical Guarantees


The reason for shrinking the radius by γ^-1/2   
in the above algorithm comes from applying tools from random matrix
theory <cit.> to sparse matrices—see the supplementary
material for proofs.
The following result states that 
when Ω is sufficiently sparse and p≫ n, we can reduce 
the false positives by γ^-1 by shrinking the radius of an operator
norm ball centred around the p-dimensional identity matrix I_p
by γ^-1/2.


  
  Let Ω∈𝒰(κ,δ)  
  with κ=O(p^ν) for ν<1/2, and
  Ω_∞=o(p^1/2). 
  For some false positive rate α=γ^-s with s∈^+, 
  let the bias of the initial estimator be
  bias(Ω̂)_∞ =
    Ω̂-Ω_∞ = o(γ^sp^1/2).
  Then, 
  
    (a)     Ω̂_s+1-Ω_∞/Ω̂_s-Ω_∞ = 1 + o(1)/γ^-1/2 + o(1)   (b)     Ω̂_s+1-I_p_∞/Ω̂_s-I_p_∞ = 1 + o(1)/γ^-1/2 + o(1)

 


  To control the false positive rate, we require a few assumptions
  in the theorem statement.  Namely, the number of non-zero entries
  per row κ and the operator norm of Ω cannot grow
  at a rate faster than p^1/2 as p increases.  In the 
  simulations of Section <ref>, we have the much 
  nicer setting where these quantities remain bounded as p increases.


  The accuracy of the conclusion of Theorem <ref>
  is improved as p increases for a fixed sample size n.
  This fact motivates the subsampling methodology in 
  Section <ref>.  
  Secondly, as the false positive rate decreases, the
  bias is allowed to be larger.  Hence, we can choose γ^-s
  based on p.  This is because as s increases, we threshold
  more aggressively effectively fighting against increases in the
  bias.


As noted in Remark <ref>, we chose the method of 
<cit.> for our simulations in Section <ref>.
It is shown in their work that this debiased estimator has the following 
maximal entrywise bias,





b_max = 
max_i,j{bias(Ω̂)} = 
O(κlog(p)n^-1/2)
or =O(κ^3/2log(p)n^-1/2) depending on specific assumptions
on the sub-Gaussian nature of the data.
Hence, 
bias(Ω̂)_∞=o(b_max p^1/2γ^-s).
As we have control over the false positive rate γ^-s, we can 
choose this to satisfy the conditions of the above theorem.
Namely, considering false positive rates less than p^-1/2
or p^-3/4, which are generally of more interest than large
false positive rates.


Further considering the debiased estimator described in 
<cit.> and given the assumptions made in that
article, we can recover the support of the matrix asymptotically
as n,p increase and α→0.
Indeed, Lemma 9 from <cit.> and equivalently
Theorem 1 from <cit.> assume an 
irrepresentability condition common in the lasso literature
and get convergence rates of the graphical lasso estimator 
depending on the tail behaviour of the random vectors X_1,…,X_n.
Thus, the graphical lasso estimator asymptotically recovers the support.
Further assuming sub-Gaussian tails for the X_i, 
<cit.> shows that the remainder term in the
debiased graphical lasso estimator is asymptotically negligible.
Thus, thresholding the debiased estimator will in turn re-recover
the support.  
Similarly, the literature on sparse covariance matrix estimation 
generally requires sub-Gaussian tails for asymptotic support 
recovery 
<cit.>.

Beyond distributional assumptions, a quick calculation can
demonstrate that the true positive probability is 
necessarily greater than the false positive probability.
For some threshold t_s∈[0,1] corresponding to a false
positive rate of α = γ^-s, we have assuming
symmetry of the distribution of Ω̂_i,j that

    P(true positive) 
         = 
      P(
        Ω̂_i,j > t_s | Ω_i,j0
      ) 
      
       =
      2P(
        Ω̂_i,j-Ω_i,j > t_s-Ω_i,j | Ω_i,j0
      ) 
      
       =
      P(
        Π̂_i,j > t_s'
      )
      = P(false positive for t_s')

which is that the probability of a true positive 
at threshold t_s
corresponds to the probability of a false positive
at threshold t_s'>t_s.  Thus, this method is 
guaranteed to return at least as many true positives
proportionally as false positives and generally,
as will be seen in Section <ref>,
performs much better.

More interesting is to consider the finite 
sample performance of our estimator.
Under significantly fewer assumptions than the results in those 
previous articles, 
we can extend the previous theorem to give us a rough guarantee
on the true positive recoveries by effectively treating
them as false positives for a different estimator.


  
  Given the assumptions of Theorem <ref>,
  let t_α>0 be the maximal threshold such that for the
  initial estimator Ω̂, we have 
  P(Ω̂_i,j>t_α | Ω_i,j=0) ≤α.
  Then, for any i,j,
  
    P(Ω̂_i,j≥ t_α | Ω_i,j0
        ) ≥α(t_α-Ω_i,j)^-2

  and
  
    inf_i,jP(Ω̂_i,j≥ t_α | Ω_i,j0
        ) ≥α(1-δ)^-2



This implies that the worst case probability of a true positive
is directly related to the false positive rate α 
and δ, the lower bound on nonzero entries in Ω.
In the simulations of Section <ref>, we see this
occurring in the case of block diagonal matrices where, even
though the true positive rate is quite low, it is still greater
than 2.25α as can be seen in Figure <ref>.





 §.§ Subsampling



As Theorem <ref>, becomes more accurate for 
large p, we can run the above methodology
in parallel by randomly partitioning the sample of size n
into subsamples of size n/k.
As a result, we run our method in parallel k times returning
estimators Ω̂^(1),…,Ω̂^(k).
For a false positive rate of 0<α≪1 and the assumption 
from Theorem <ref>
that Ω∈𝒰(κ,δ) with 
κ=O(p^ν) and ν<1/2, then the expected number of 
false positive recoveries is 
α[p(p+1)/2-Cp^1+ν] = O(α p^2).  
After splitting the sample into k disjoint pieces, 
we can run the above algorithm in parallel to construct
k independent estimators Ω̂^(1),…,Ω̂^(k).
The probability
of recovering the same false positive entry in d or more of the 
Ω̂^(i) is a binomial tail sum: 
∑_i=d^k k iα^i(1-α)^n-i.
When α < 0.001, then choosing d=2 is 
generally sufficient.
Indeed, we can bound this binomial tail sum by 
the Kullback-Leibler divergence. 
For Bkα, 

    P(B≥ d) 
      ≤exp{
          -dlog( d/k/α)  
          -(k-d)log(1-d/k/1-α) 
      }
    ≤ k^k(α/d)^d
           (1-α/k-d)^k-d≤k^k/d^d(k-d)^k-dα^d
      ≤ 2^kα^d.

If we want a target false positive rate of 2^-10, 
we can choose α such that 
2^kα^d = 2^-10 or α = 2^-(10+k)/d
to achieve the same false positive rate as when not
subsampling.  For example, if we take d=2 and keep k<10, then we 
can relax the false positive rate for each individual 
subsampled estimator.



This addition of subsampling to the methodology allows for 
faster runtimes, as the k estimators can be computed in 
parallel with fewer iterations, and also increases accuracy
as the target false positive rate decreases below p^-1
as will be seen in Section <ref>.




§ APPLICATIONS




 §.§ Simulated Data



We test our methodology for three different graphical structures.
In the first case,
Ω tridiagonal–i.e. Ω_i,i=1, Ω_i,j=1/3
if i-j=1, Ω_i,j=0 if i-j>1.  This is 
the sparsest connected graph structure with p-1 edges on p
nodes.
In the second case,
Ω represents a binary tree graph with depth d
resulting in d(d+1)/2 nodes and d(d-1) edges, which is 
roughly 2p edges for p nodes.  Similar to case one,
diagonal entries are set to 1 and non-zero off-diagonal entries
are set to 1/3. 
In the third case,
Ω represents a disjoint collection of k complete graphs
on d nodes.  Thus, there are kd nodes and kd(d+1)/2 edges,
and Ω is a block diagonal matrix with blocks 
1/31_d + 2/3I_d with 1_d being
the d× d matrix of all 1's and I_d, the d× d
identity matrix.

For simulations, the sample size was fixed at n=50.
The dimensions considered are p=496,1128,2080 as they
correspond to binary trees of depth 31,47,64, respectively,
and correspond to block diagonal matrices with 
62,141,260 (8×8)-blocks, respectively. In all nine 
cases considered, the non-zero off-diagonal entries are all 
set to 1/3 with the diagonal entries set to 1.

In Figure <ref>, we plot the empirical 
true and false positive rates on the vertical axis against
the target false positive rate on the horizontal axis.  The
target false positive rate acts as a tunable parameter.
Three types of matrices, three dimensions, and three values
of the graphical lasso regularization parameter = 0.5,1,2 are 
considered.
Ideally, the observed false positive rates will be close
to the solid black line, which is where the observed and 
target rates coincide.  This occurs as the dimension
increases in all three cases.  For the sparser models–tridiagonal
and binary tree structures–we are able to maintain a true
positive rate above 20% as the false positive rate is taken to zero.
Performance is much worse when trying to recover the 
block diagonal matrix.  For the lowest dimension, p=496,
setting the penalization parameter =0.5 gives the best
results whereas in the higher dimensions =1 gives the 
best performance albeit only slightly.  Hence, it appears that
for choosing an initial estimator, the graphical lasso penalization parameter 
can be set to 1 regardless of dimension and true Ω 
for this methodology.
On this note, we also considered =0.25
but the performance was terrible and thus it is not included
in the figures.  For an alternative look at the simulations,
receiver operating characteristic curves are included in the 
supplementary material.
The same simulations were also run on multivariate Laplace data
with details in the supplementary material.  In short, 
choosing =1 allowed for the empirical false positive 
rate to stay close to the target rate whereas 
=2 did not perform as well.


In Figure <ref>, we consider the tridiagonal and 
binary tree models estimated by the standard approach as before
and by subsampling as in Section <ref> with k=2.
By subsampling, we can extend the false positive rate from 
around α=p^-1, where the standard method begins to 
breakdown, all the way to around α=p^-2.  Hence, 
we can have a controlled reduction in the empirical false positive 
rate to effectively zero, as the precision matrix only has p^2 
entries, using the subsampling method.

In Table <ref>, we tabulate empirically achieved true
and false positive rates from the cross validated graphical lasso 
<cit.> and 
ACLIME estimators <cit.>. These estimators were 
computed by splitting the n=50 sample in half and optimizing
the tuning parameters with respect to the operator norm distance.
The graphical lasso estimator's achieved false positive was roughly 
halved as the dimension doubled.
The ACLIME estimator much more aggressively penalized the 
matrix entries.  
The threshold method of <cit.> was also tested, which
for a false positive rate α, sets the individual entries
in the debiased estimator Ω̂ to zero if

    Ω̂_i,j < 
      Φ^-1(1-α/p(p-1))
      Σ̂_i,j/√(n)

where Σ̂ is the empirical covariance matrix and Φ(·) is
the cumulative distribution function for the standard normal distribution.  
A figure in line with Figure <ref> for this method 
is included in the supplementary 
material.  In short, this approach should work asymptotically as 
n,p→∞.  In
our simulations, the empirically achieved false positive rate was not 
close to the target false positive rate.










 §.§ Geonomics data


We apply our methodology to the dataset of gene expressions
for a small round blue cell tumours mircoarray experiment from 
<cit.>, which is also analyzed in other works 
<cit.>.
The data set consists of
a training set of 64 vectors containing 2308 gene expressions.
The data contains four types of tumours.
Considering this X∈^64×2308, we construct
error controlled precision matrix estimators using the 
algorithm from Section <ref> choosing the
debiased estimator of <cit.> as the initial 
estimator and for graphical lasso penalization parameter
=0.5,1,2.

In Figure <ref>, we plot the cardinality of the 
support of the estimators for the three penalization 
parameters as the false positive rate ranges from 
2^-1 to 2^-10.  For =0.5, the method
plateaus quickly returning a small support of 74 nonzero
off-diagonal entries.  The other two lines continue to decay
towards zero off-diagonal entries.

In Figure <ref>, we see kernel density plots of
the number of non-zero entries per row aggregated over all 
2308 rows.  In these plots, 
= 0.5,1,2 and α=0.5,0.001 are considered.
Comparing the top and bottom row of plots, we see that 
there are many rows with approx 500 nonzero entries 
when α=0.5, which quickly drops to fewer than 
10 nonzero entries per row once α=0.001.








§ SUPPLEMENTARY MATERIAL


Supplementary material available at Biometrika online includes
proofs for the main results and auxiliary results. 
Receiver operating characteristic curves for the data in 
Figure <ref> are included, and
additional simulations for exponential data and the method from 
<cit.> are also included.





biometrika





§ PROOFS


The theory behind the methodology in this article comes from 
redoing the results from random matrix theory, which generally
achieves an operator norm of O(p^1/2),  under the assumption
that only α∈(0,1) of the entries are non-zero.  This
gives similar results but with an operator norm 
of O(α^1/2p^1/2).

We first prove a lazy version of Hoeffding's concentration
inequality.

  
  For some M>0,
  let Z_1,…,Z_p∈[-M,M] be independent mean zero random 
  variables not necessarily equally distributed,
  and let b_1,…,b_pα also independent
  of the Z_i.  Then,
  
    P(∑_i=1^p b_iZ_i≥ tM(α p)^1/2) ≤ 2^-t^2/4.

 

  The proof is more or less the same as that of the standard
  Hoeffding inequality.  Without loss of generality, 
  let M=1.
  For any i=1,…,p, we can bound the variance of b_iZ_i
  by 
  
    b_iZ_i = b_i^2Z_i^2≤α.

  Thus, we have a bound on the cumulant generating function
  of b_iZ_i that
  log^ b_iZ_i≤α^2 and 
  log^∑_i=1^p b_iZ_i≤α p^2.
  Therefore, via Chernoff's inequality, we have for any >0 that
  
  
    P(∑_i=1^p b_iZ_i ≥ t)
        ≤P(^∑_i=1^p b_iZ_i≥^ t)
        ≤^- texp(
          ∑_i=1^p b_iZ_i
        )
        ≤exp(
          - t + pα^2
        )

  
  Minimizing the righthand side 
  over  gives = t/(2pα) and finally
  that 
  P(∑_i=1^p b_iZ_i ≥ t) ≤^ -t^2/(4α p) .
  Running the proof for ∑_i=1^p b_iZ_i ≤ -t
  gives the reverse inequality.  Combining those gives
  the desired result.
  Lastly, to adjust for M1, replace Z_i with Z_i/M
  in the above result.


Using the lazy version of Hoeffding's inequality, we can 
get concentration of the operator norm for a 
sparse random matrix using a modification of the 
standard -net argument as presented in <cit.> Section 2.3.


  
  Let A∈^p× p be a symmetric random matrix
  with entries A_i,i=0, A_i,j<1, and 
  A_i,j=A_j,i.  Let B∈{0,1}^p× p be 
  a symmetric Bernoulli random matrix with iid entries
  such that P(B_i,j=1)=α∈[0,1].  If the lower 
  triangular entries of A are independent then,
  
    A∘ B_∞ = O{ (α p)^1/2}

  where A∘ B denotes the entrywise or Hadamard product
  of A and B.


  For now, we consider A and B as iid ensembles—i.e. remove the 
  A_i,j=A_j,i and B_i,j=B_j,i condition—and adjust for the
  symmetry at the end of the proof.
  Let Π = A∘ B.
  For any vector v∈^p such that v_2=1,
  we have that 
  Π v_2 = 
    (∑_i=1^p(∑_j=1^pA_i,jB_i,jv_j)^2)^1/2.
  Therefore, from the proof of Lemma <ref>, independence 
  of entries, and v_2=1, 
  
    log^Π v_2 =  
        ∑_i=1^plog^∑_j=1^p A_i,jB_i,jv_j≤α p^2.

  Thus, for any arbitrary unit vector v, we can 
  apply Lemma <ref> to get
  
    P(Π v_2 > t(α p)^1/2)
        ≤ 2^-t^2/4.

  To extend to Π_∞, we construct a maximal 
  (1/2)-net of the unit sphere S = {v : v_2=1}.
  That is, let V be the maximal set of points in S such that 
  for any u,v∈ V, u-v_2≥0.5.  Let v^*∈ S be
  the vector such that Π_∞ = Π v^*_2,
  which exists via compactness.
  Thus, there exists a v_0∈ V such that v^*-v_0_2<0.5
  as otherwise, the set V would not be maximal.
  Furthermore, Π(v^*-v_0)_2≤Π_∞/2
  and thus Π v_0_2≥Π_∞/2.
  Therefore, by the union bound,
  
    P(Π_∞ > t(α p)^1/2) ≤P(⋃_v∈ V{Π v_2 > t/2(α p)^1/2}) ≤
        2V^-t^2/16.

  From Lemma 2.3.4 of <cit.>, V = (2C)^p for
  some absolute constant C>0. Thus, for t large enough,
  the righthand side becomes negligible.
  
  Now, considering A and B symmetric with zero diagonal
  as in the theorem statement, we have that 
  Π_∞≤Π_low_∞ + 
    Π_up_∞
  for Π_low and Π_up the strict lower and upper 
  triangular portions of Π.  Hence,
  
    P(Π_∞ > t(α p)^1/2) ≤P(Π_low_∞ > t(α p)^1/2/2
      ) +
      P(Π_up_∞ > t(α p)^1/2/2
      )
      ≤
      4V^-t^2/32.



  The constants in the above proof are not necessarily
  optimal, but sufficient to justify our approach to controlled 
  precision matrix estimation.



  For (a), Ω̂_s-Ω
  has zero diagonal and off-diagonal entries bounded by two.
  Further, if 
  Ω∈𝒰(κ,δ), then 
  Ω_∞≤max_i∑_j Ω_i,j≤κ.
  Hence, if κ = O(p^ν), then so is Ω_∞.
  For the bias, let B_s be a Bernoulli random 
  matrix with probability γ^-s that entry B_i,j is 1.
  Then,
  
    bias(Ω̂_s)
        = (Ω̂_s)-Ω
        = ( B_s∘Ω̂)-Ω
        = γ^-sbias(Ω̂) + 
          (1-γ^-s)Ω = o(p^1/2).

  Thus, we have by assumption and by applying 
  theorem <ref> as well as the Bai-Yin 
  theorems in <cit.> Section 2.3 that
  
    Ω̂_s-Ω_∞≤Ω̂_s-Ω-bias(Ω̂_s)
        _∞ + bias(Ω̂_s)_∞≤
    ≤ 
        (2+o(1))p^1/2γ^(s-1)/2 + o(p^1/2)

  and
  
    Ω̂_s-Ω_∞≥Ω̂_s-Ω-bias(Ω̂_s)
        _∞ - bias(Ω̂_s)_∞≥
    ≥
        (2+o(1))p^1/2γ^(s-1)/2 - o(p^1/2).

  Therefore,
  
    Ω̂_s+1-Ω_∞/Ω̂_s-Ω_∞ = 
       (2+o(1))p^1/2γ^(s-1)/2 + o(p^1/2) / (2+o(1))p^1/2γ^s/2 + o(p^1/2) 
      = 1 + o(1)/γ^-1/2 + o(1).
  
  
  For (b),
  let K = Ω-I_p_∞ which is 
  K≤max_i∑_j=1^pΩ_i,j≤κ = o(p^1/2)
  by the sparsity assumption.
  
  
  
    Ω̂_s+1-I_p_∞/Ω̂_s-I_p_∞≤Ω̂_s+1-Ω_∞
      	+ K 
      /Ω̂_s-Ω_∞
      	- K
       
      = 
     =
       (2+o(1))γ^(s-1)/2  + o(1)+ Kp^-1/2/ 
        (2+o(1))γ^s/2 + o(1) - Kp^-1/2
      = 1 + o(1)/γ^-1/2 + o(1).

  and similarly
  
    Ω̂_s+1-I_p_∞/Ω̂_s-I_p_∞≥Ω̂_s+1-Ω_∞
      	- K 
      /Ω̂_s-Ω_∞
      	+ K
       
      = 
     =
       (2+o(1))γ^(s-1)/2  + o(1)- Kp^-1/2/ 
        (2+o(1))γ^s/2  + o(1)+ Kp^-1/2
      = 1 + o(1)/γ^-1/2 + o(1).




  Let β_i,j = Ω̂_i,j≥ t_α
  be the probability of rejecting Ω_i,j=0, which
  is the power associated with the threshold t_α.
  Thus, with a false positive rate α, we have a 
  true negative rate 1-α and a false negative probability
  1-β_i,j.

  We begin by noting that the false negative probability for some 
  threshold t_α   
  can be related to the false positive 
  rate for some other threshold τ_i,j.
  Indeed, 
  
    1-β_i,j 
        =
        Ω̂_i,j<t_α | Ω_i,j0
        ≥Ω̂_i,j - Ω_i,j<t_α-Ω_i,j | Ω_i,j0
        
    
        =
        Ω̂_i,j^⋆<τ_i,j | Ω_i,j^⋆=0
         = 1-α^⋆

  for some τ_i,j<t_α and α^⋆>α.
  Thus, 
  we are left with the task of determining α^⋆ based on 
  this new threshold.

  From the previous results, we have that 
  Ω̂_α - I_p_∞ = (2+o(1))(α p)^1/2
  and that 
  Ω̂_α - Ω̂_∞ = 
    (2+o(1))t_α(1-α)^1/2 p^1/2
  as Ω̂_α - Ω̂ has approximately
  1-α of its entries nonzero with bound t_α on 
  the entries.  Hence, 
  
    (2+o(1))p^1/2 = Ω̂-I_p_∞≤ 
        (2+o(1))(α p)^1/2 + 
        (2+o(1))t_α(1-α)^1/2 p^1/2.

  Noting trivially that 
  Ω̂-I_p_∞/Ω̂-I_p_∞=1,
  and for p large enough to make the o(1) negligible, we have
  that 
  
    1    ≤{α^1/2 + t_α(1-α)^1/2}{(α^⋆)^1/2 + τ(1-α^⋆)^1/2}
    
        (α^⋆)^1/2   ≥{α^1/2 + t_α(1-α)^1/2}^-1 - 
        τ(1-α^⋆)^1/2






§ ADDITIONAL SIMULATIONS




 §.§ Receiver operating characteristic curves


In Figure <ref>, we consider plots of the observed
false positives against the observed true positives–that is,
Receiver operating characteristic curves.  
Similarly to the analysis of Figure 1 from the main article,
We see better performance for the tridiagonal and binary tree 
matrices than for the block diagonal matrix.  Also, the 
graphical lasso penalization parameter does not have a large effect 
on the methodology.  Though, =1,2 perform marginally better
than =0.5.  




 §.§ Sub-Exponential Data


The same simulations as in Section 3.1
were 
rerun replacing the multivariate Gaussian distribution with
the multivariate Laplace distribution, and are 
displayed in Figure <ref>.  As expected,
the true positive rate is not as large as in the Gaussian
setting.  When the penalization parameter for the 
graphical lasso is set to =2, we see that our 
methodology does not maintain the desired false positive
rate as α→0.  However, for =1, the
empirical false positive rate does approximately track 
with the desired false positive rate in the three simulation
settings.





 §.§ Asymptotic Threshold


In the work of <cit.>, a thresholding method is proposed 
for the debiased graphical lasso estimator making use of the normal distribution
function.  Namely, 

    Ω̂_i,j < 
    Φ^-1(1-α/p(p-1))
    Σ̂_i,j/√(n)

where Σ̂ is the empirical covariance matrix and Φ(·) is
the cumulative distribution function for the standard normal distribution.  
The same simulations as in Section 3.1

were run on this method
for graphical lasso penalization parameters of = 2^-4,2^-3,2^-2,2^-1,1,2.
The larger values of  resulted in all of the off-diagonal entries being
set to zero.  Hence, Figure <ref> displays the false and true
positive rates for only =1/16,1/8.  This method should work asymptotically
as n,p→∞.  However, for our specific choices of n and p,
this method failed to achieve anything close to the target false positive rate. 






