{'fraction_non_alphanumeric': 0.03853181427343078, 'fraction_numerical': 0.01633705932932072, 'mean_word_length': 4.457038123167155, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thresholded Adaptive Validation: Tuning the Graphical Lasso for Graph Recovery

The graphical lasso is the most popular estimator in Gaussian graphical models, but its performance hinges on a regularization parameter that needs to be calibrated to each application at hand. In this paper, we propose a novel calibration scheme for this parameter. The scheme is equipped with theoretical guarantees and motivates a thresholding pipeline that can improve graph recovery. Moreover, requiring at most one line search over the regularization path of the graphical lasso, the calibration scheme is computationally more efficient than competing schemes that are based on resampling. Finally, we show in simulations that our approach can improve on the graph recovery of other approaches considerably.


Introduction
Over the last decades full of technical achievements, we experienced a revolution of the field of Data Science, confronting us with large-scale data sets. In order to handle, and to infer new insights from the appearing wealth of data, it presupposes us to put effort into the development of new, scaleable procedures. One approach to address this problem are graphical models, which proved to serve as an intuitive, easy-understanding visualization of the underlying interaction network that can then be further analyzed. Typical applications for graphical models occure in several modern sciences, including genetics (Dobra et al., 2004), the analysis of brain connectivity networks (Bu and Lederer, 2017), or the investigation of complex financial networks (Denev, 2015). In all of these cases, graphical models can reduce the network of interactions into its relevant parts to lighten the challenge of high-dimensionality. Similar to classical dimension reduction algorithms like Principal Component Analysis (Pearson, 1901;Hotelling, 1936;Jolliffe, 1986), our goal is to find these relevant features, structures and clusters in the underlying data. In the context of graphical modelling, we call this task graph recovery. A popular approach to face this challenge is to consider undirected graphical models, which structure the underlying data by representing the conditional dependence structure of the measured parameters.
An important class of undirected graphical models are Gaussian graphical models. There are numerous estimators for Gaussian graphical models, including such that account for high dimensionality: the graphical lasso (Yuan and Lin, 2007;Banerjee et al., 2007;Friedman et al., 2008), SCAD (Fan et al., 2009), and MCP (Zhang, 2010), which are based on the idea of regularized Maximum-Likelihood estimation. But there are also completely different approaches such as MB (Meinshausen and Bühlmann, 2006), TIGER (Liu and Wang, 2012), or scaled lasso (Sun and Zhang, 2012) that rely on piecewise estimation of the graph. These estimators reduce the effective dimensionality of the model through a regularization term that is adjusted to the setting at hand with a regularization parameter. The regularization parameter should not only satisfy theoretical finite sample properties, but should also be efficiently computable as the complexity of a graph scales quadratically in the amount of variables. This paper tackles the problem of regularization parameter selection of the graphical lasso for graph recovery tasks under consideration of these 2 aforementioned points: finite sample theory and computational efficency. Although, arXiv:2005.00466v1 [stat.ML] 1 May 2020 the graphical lasso being one of the simplest among them, we can make use of its outstanding theoretical properties (Ravikumar et al., 2008) to use a technique, originally introduced by Chichignoud et al. (2014) for regularized linear regression, to calibrate the graphical lasso. Important features of the proposed data-driven calibration scheme are that it comes with a finite sample result, and is computationally efficient as it requires only 1 graphical lasso path.
The rest of the paper is organized as follows. In Section 2 we give a short introduction to the graphical lasso. Section 3 reviews the calibration scheme including its finite sample theory that allows us to motivate a thresholded graphical lasso approach, which we call thresholded adaptive validation graphical lasso (thAV). Running multiple simulations in Section 4, we verify that the proposed method is indeed suited for graph recovery tasks. We apply the thAV on real-world data to recover biological networks in Section 5. Finally, we conclude with a discussion in Section 6 to present generalizations of the proposed method and ideas for future research. The software and experiments are provided in our repository https://github.com/MikeLasz/thav.glasso . All proofs are deferred to the Appendix.

Brief Review of Gaussian Graphical Models
The theory of undirected graphical modelling seeks to unravel the relationship structure between components of a highdimensional random variable z given some samples. To be precise, it aims at identifying the conditional dependence structure of a real-valued random variable z ∈ R d , that is, to find out whether where z i , z j are the ith and jth component of z, respectively, and z \{i,j} are the remaining d − 2 components. This leads us to the definition of an undirected graphical model as a pair (z, G), where G := (V, E) is a graph with vertices V := {1, ..., d} and edge set E := {(i, j) ∈ V × V : z i ⊥ z j |z \{i,j} }. In regard to the large theoretical background of the Gaussian distribution, we want to consider z ∼ N d (0 d , Σ), where Σ ∈ R d×d is a positive definite covariance matrix. It is well-known that in the Gaussian case we can find an elegant characterization of the conditional independence structure: for any i = j ∈ V, where Θ := Σ −1 is the so called precision matrix. The equivalence (1) is also referred to as a special case of the Hammersley-Clifford Theorem (Grimmett, 1973;Besag, 1974;Lauritzen, 1996). Hence, in order to estimate the conditional dependence graph G, one aims at estimating the precision matrix, which in turn allows to estimatê whereΘ is an estimate of Θ. Let z (1) , ..., z (n) be n samples drawn from z ∼ N d (0 d , Σ). Being a well-studied workhorse in statistical estimation, we apply Maximum-Likelihood estimation (MLE) what results in where S + d is the set of positive definite and symmetric matrices in R d×d , see chapter 7 in Giraud (2015). Even though there is a large theoretical background on MLE (Lehmann and Casella, 2003), most results are of asymptotic nature. Moreover, in the high-dimensional setting where n ≈ d or n < d, it is well-known that MLE (if it even exists) performs badly if one aims to recover the relation structure. A typical approach to overcome the burdens that come with high-dimensionality is to assume a sparsity structure on the target, that is, we assume Θ to have many zero-entries. It does not only improve theoretical results of the estimation, but it also makes the graphical representation of the conditional dependence graph more interpretable. Further, imposing a sparsity structure is in accordance with the scientific beliefs in typical areas in which graphical models are being used (Thieffry et al., 1998;Jeong et al., 2001). Applying a 1 -regularization in (3), we end up with the probably most-frequently used estimation procedure for Gaussian graphical models: the graphical lasso (Yuan and Lin, 2007) Θ r = argmin where r is a problem-dependent regularization parameter, and Ω 1,off := i =j |Ω ij | denotes the 1 -norm of Ω ∈ S + d on its off-diagonal. And while there is much theory on the properties of the graphical lasso (Rothman et al., 2008;Ravikumar et al., 2008;Jankova and van de Geer, 2018), as well as general theory for graphical models (Zhuang and Lederer, 2018), there is no theory that includes the practical choice of r.

Thresholded Adaptive Validation
In this section we derive a novel scheme for calibrating the regularization parameter of the graphical lasso. We call it the thresholded Adaptive Validation Regularization Parameter Calibration (thAV). It is based on the adaptive validation scheme for calibrating the lasso regularization parameter (Chichignoud et al., 2014), which we generalize to our setting in Section 3.1. The proposed approach along with its theoretical guarantees is introduced in Section 3.2. Its computational efficiency is then discussed in Section 3.3.

Adaptive Validation Regularization Parameter Calibration
In the following, we review the adaptive validation scheme (Chichignoud et al., 2014) and apply it to calibrate the graphical lasso. We show that this leads to ∞ -bounds on the graph's off-diagonal, as described by Theorem 4 in the end of this section. These bounds will also be the basis for the thresholding approach described in the following section.
Let us fix a finite set of possible regularization parameters R. Definition 1 (AV). The Adaptive Validation calibration scheme (AV) selects the regularization parameter according tô where C is a constant and ∞,off is the elementwise ∞ -norm on the off-diagonals We callΘr the AV graphical lasso, or simply AV.
Note that C is a constant that we specify in the following. For the ease of notation we will write := ∞,off and r :=r C henceforth. Further, to obtain theoretical results of the AV regularization parameter, we need to impose a setting which can be generalized by the following theoretical assumption. Assumption 2. Let T r r∈R be a class of events indexed by r ∈ R. Assume that T r is monotone in r in the sense that T r ⊂ T r for r ≤ r .
Suppose further that, conditioned on T r , it holds that whereΘ r is the graphical lasso (4), Θ is the true precision matrix and C is the constant used to define the AV.
The validity of Assumption 2 under standard assumptions is shown in the Appendix. Nevertheless, in practice, we can not refer to these theoretical properties as the proposed C depends on quantities that are unknown in practice. One might think that we simply swap one parameter, the regularization parameter r, with another parameter, C. However, the big difference is that the final thresholded estimator (Section 3.2), is not decisively dependent on the particular choice of C but rather returns a stable estimate for a range of C (Section 4.3).
In regard of Assumption 2, it seems natural to seek for the smallest regularization parameter such that we can apply (5) with some fixed probability 1 − δ. This regularization parameter is denoted as the oracle regularization parameter: Definition 3 (Oracle Regularization Parameter). Consider a constant δ ∈ (0, 1). Then, the oracle regularization parameter is given by In practice, the oracle estimator is inaccessible as we cannot measure P(T r ), yet, using the AV construction, we obtain a graphical lasso estimateΘr that performs almost as well as the oracle graphical lassoΘ r * δ . In fact, it is a very common strategy in high-dimensional statistics to define an oracle tuning parameter based on an oracle inequality as in (5) to motivate an estimation target. For instance, a well-known oracle inequality in lasso regression is conditioned on the set T r := {r ≥ 2 X ε ∞ /n} (Bühlmann and Geer, 2011), where the quantity 2 X ε ∞ is also known as the effective noise. Using (7), one might calibrate the lasso estimator for linear regression by estimating the quantiles of the effective noise, or equivalently speaking by estimating argmin r∈R {P(T r ≥ 1 − δ} such as it has been done in Lederer and Vogt (2020). As we will see in the Appendix, Theorem 12, the equivalent to T r in graphical lasso estimation can be defined as where κ Σ ,Σ emp , deg(Θ), and α are defined in Section A.1 in the Appendix. Hence, the quantity 8 max{ ∞ (Σ emp , Σ), κ Σ /(6 deg(Θ))} can be interpreted as the effective noise for the graphical lasso. This effective noise can then be used to define an oracle regularization parameter as it has been done in (6). Theorem 4 (Finite Sample Bound for the AV). Suppose that Assumption 2 holds and thatr is the regularization parameter selected by the AV method. Then, for any δ ∈ (0, 1), it holds that with probability at least 1 − δ.
The proof follows the same steps as in Chichignoud et al. (2014) and can be found in the Appendix.

Thresholding the AV for Graph Recovery
With Theorem 4, we have an outstanding result to address the estimation of Θ in -loss. However, it turns out that, being an overly dense estimation,Θr is unsuited to recover a sparse dependence structure. Nonetheless, we can still make use of Theorem 4 by pursuing a thresholded graphical lasso approach to obtain graph recovery properties.
For simplicity of notation, let us denote the AV estimator byΘ :=Θr from here on.
Definition 5 (thAV). LetΘ be the AV estimator. Then, we define the thresholded adaptive validation graphical lasso (thAV) by where t := λCr is the threshold, λ ∈ (0, 3], and 1 A is the indicator function over the set A.
As far as we know, there is no theoretical finite sample result that justifies a specific choice for a threshold in a thresholded version of the graphical lasso. However, our proposed method comes with a threshold for the AV which is based on the finite sample guarantees that are stated in Theorem 4.
1. Consider the set of zero entries of Θ, denoted by S : 2. Consider the set of significant true entries of Θ, denoted by S sig (λ) : Hence, part 1 of Corollary 6 can be interpreted as a false positive control result. In fact, with probability 1 − δ, thAV only includes false positives if the corresponding graphical lasso estimate is significantly wrong, that is, |Θ ij | > λCr even though Θ ij = 0. Proposing λ = 3, the rationale is first, that we suspect the graphical lasso to perform that good, that it assigns false positives a small absolute value, and secondly, thatr is a reasonably good approximation of r * δ which would result in a small interval (3Cr, 3Cr * δ ] in which thAV missclassifies zero-entries. Furthermore, part 2 of Corollary 6 is a result that ensures false negative control: with probabilty 1 − δ, we include every edge that is significantly large. Hence, we can specify the set of false negatives as S\S sig (λ), where S := {(i, j) ∈ V × V : Θ ij = 0} is the set of true edges, which is also called the support of the graphical model. According to our needs, we can regulate λ to decrease the false negative rate (part 2, as S sig (λ) is decreasing in λ) at the cost of increasing the interval in which thAV missclassifies true negatives (part 1), what increases the false positive rate. The proof of Corollary 6 is a direct consequence of (8) and can be found in the Appendix.
Motivated by these results, we propose the thAV estimator (9) using t := 3Cr, which we use in the following simulations in Section 4. In summary, we can expect the thAV to have an attractive performance in terms of false positive rate and to be able to recover significant edges. The simulations verify these results and expectations on the thAV.

Computational efficiency of the thAV
In addition to its finite sample guarantees, thAV also comes with notable computational benefits, which make it possible to efficiently compute it even if d is in the thousands. In comparison to calibration schemes, which are based on data splitting and resampling methods like StARS or other CV-like approaches, thAV only requires at most one solution path (Θ r ) r∈R . We can readily implement thAV using Algorithm 1. We observe that using Algorithm 1, we 1. calculate at most 1 graphical lasso path since we can store and reuse every graphical lasso solution; 2. start to calculate graphical lasso solutions for large regularization parameters and decrease the regularization until we break. Therefore, we most likely avoid to calculate the graphical lasso for very small regularization parameters, which are computationally most demanding; 3. can use warm starts, which are implemented in the glasso R package (Friedman et al., 2008;Witten et al., 2011) to calculateΘ r , decreasing the computation time once more.
Empirically, we can see that the AV estimator does not only overtop other calibration schemes that are based on data splitting and resampling methods like CV and StARS in runtime, but also in computationally stability and memory usage.

Simulation study
In this section we want to compare our proposed method with 2 other commonly used regularization parameter calibration schemes for the graphical lasso, which are the StARS (Liu et al., 2010) and the RIC (Zhao et al., 2012). We sample synthetic data from z ∼ N d (0 d , Θ −1 ), whereby we adopt the precision matrix generation procedure from (Caballe et al., 2015): • random graph: First, we generate a Gilbert graph (Gilbert, 1959), that is, we independently connect two nodes i = j with some fixed probability p. The default value is p = 3/d. Then, we set the non-zero pattern of Θ according to the adjacency matrix of the Gilbert graph. The second step is to assign all non-zero entries of the precision matrix a value. To do so, we set where U (I) is the uniform distribution over I = [−0.9, −0.5] ∪ [0.5, 0.9] and Adj ∈ {0, 1} d×d is the adjacency matrix corresponding to Θ. To guarantee the positive definiteness of Θ, we use the fact that where λ min (Θ) is the smallest eigenvalue of Θ, µ is some constant, and I d×d ∈ R d×d is the identity matrix. Hence, (10) suggests to add a diagonal matrix with µ := min{ε, − λ min (Θ) 1 } where · 1 denotes the floor function that returns the downwards rounded to 1 decimal place, and ε is some small positive constant. Finally, we scale the precision matrix to obtain ones at the diagonal: • scale-free graph: In this case, the adjacency matrix is generated according to a Barabasi-Albert Algorithm (Barabasi, 2016) that generates the corresponding graph subsequently. It starts with 2 connected nodes and adds sequentially a node that will be connected with one of the existent nodes in the network. The probability to be connceted to a node i is proportional to its degree, that is, where deg(i), deg(j) denote the degree of the nodes i and j respectively. We repeat this step, till we have a graph with d nodes. From here on, we proceed as in the above case and generate values forΘ ij , ensure positive definiteness, and finally scale the precision matrix. Graph 1(a) represents a random graph with p = 3/100, and graph 1(b) is a scale-free graph according to the described Barabasi-Albert Algorithm.
The generation of both types of precision matrices is implemented in our attached repository. Figure 1 shows an example of both graph structures. In contrast to usual precision matrix generation methods, a Θ generated according to this 6 (a) (b) (c) Figure 2: Histogram of the absolute valued non-zero weights |Θ ij | of 50 scale-free and 50 random graphs according to our generation procedure described in Section 4. We considered the cases d ∈ {200, 300, 500} in Figure 2(a), 2(b), and 2(c), respectively.
method includes values Θ ij that can be either positive or negative, and vary in absolute value. However, a disadvantage of this generation method is that the entries of a scale-free graph tend to be differently distributed than the ones of a random graph. Especially, the absolute value of the signals of a scale-free graph is biased towards lower magnitudes, see Figure 2. This observation suggests that a scale-free graph might be, not only due to its topology, but also because of its smaller signal strengths, harder to estimate.
In our simulations, we use the pre-implemented functions in the R package huge (Zhao et al., 2012) to generate the adjacency matrices. In the following, if not stated differently, we use t = 3Cr and C = 0.3 as suggested in Section 4.3. The set of possible regularization parameters is defined as the grid R := {0.2 · 0.95 i : i ∈ {1, ..., 40}}. We set the probability to connect two nodes in a Gilbert graph (used to generate a random graph) to p = 3/d. RIC and StARS are both with their default settings in huge. The results of all simulations are averaged over 50 iterations and standard deviations are shown in brackets. The software, containing all simulations and definitions, can be found in our provided repository.

Performance in F 1 -score
To evaluate the graph recovery performance we consider the F 1 -score, which is the harmonic mean of the precision and the recall, which are defined as where E,Ê are the edge sets obtained from the non-zero patterns of Θ and some estimateΘ, respectively. Hence, precision measures the ratio of correctly assigned edges and the number of selected edges, recall on the other hand measures the ratio of correctly assigned edges and the number of all true edges.The F 1 -score puts both measures into In the best case, where we correctly assign all edges, that is, ifÊ = E, we get a F 1 -score equal to 1. The worst F 1 -score is 0, if either precision or recall is equal to 0.
As precision and recall behave very similar to the false positive rate and false negative rate, respectively, these quantities serve as a witness for the validity of Corollary 6, whereas the F 1 -score acts as a general performance measure. Starting with the results for a random graph, we can see from Table 1 that our procedure not only clearly outperforms the state-of-the-art methods, but also has a big advantage over the oracle estimator, which is the (non-thresholded) graphical lasso estimator that achieves maximal F 1 -score among all regularization parameters. Remarkably, we observe that the thAV always achieves a precision very close to 1.00 what means that we include just a few false positives. This goes in line with the theoretical result in part 1 of Corollary 6. Figure 3 shows one exemplary graph recovery of a random graph and further recovered graphs can be found in Figure 4. Moreover, it can be seen that the scale-free graph is in general much harder to estimate but thAV still remains superior to the RIC and StARS, and even to the oracle estimator. Nevertheless, we observe that the estimator suffers to maintain a good recall, but still, the precision of thAV remains close to 1. This means that the resulting estimate comes with a 9 reasonable amount of false negatives but only a few false positives. The results of the simulation in Table 2 demonstrate that we can substantially increase the F 1 -score by thresholding by 2rC instead of 3rC in the case of a scale-free graph. In addition, we make some investigations on the performance of the thAV on large-scale examples in Table 3, that is, we dramatically increased the dimensionality d. Surprisingly, in the case of a random graph, we observe that thAV has an outstanding performance for every d ∈ {600, 700, ..., 1000} that does not even decay if we increase d. On the other hand, thAV quite suffers to recover the scale-free graph for this setting. But, remarkably we can still achieve a precision of almost 1 in every investigated case, which comes at the costs of a low recall . Once more, this affirms that thAV confidently recovers a sub-network of the true network, what goes in line with Corollary 6.

Recovery of significant edges
This section aims to verify the second part of Corollary 6, which states that we can reliably recover significant edges. To make investigations on the recovery of these edges, we consider E 50 , which is the edge set containing the 50 most significant edges (i, j) of Θ, that is, we consider the pairs of distinct nodes (i, j) such that |Θ ij | are the 50 largest values among all distinct nodes (i, j). Recall, part 2 of Corollary 6 states that we can recover all edges |Θ ij | > 6Cr * δ with probability 1 − δ. Although r * δ being unknown, we know with Theorem 4 thatr ≤ r * δ with high probability. In this case we can suspect the thAV to recover any entry that is notably larger than the reference quantity 6Cr. Table 4 demonstrates an exemplary run for various settings and compares |Θ ij | for (i, j) ∈ E 50 with 6Cr. Denote the 50th largest off-diagonal entry of Θ as top 50 (Θ). Then, we observe that the difference between top 50 (Θ) and 6Cr of a random graph is notably larger than the same quantity in a scale-free graph. As can be seen from Table 4 and Table 5, the larger difference goes in line with a better recovery of the top-50 edges. This, once more, is in accordance to our theory.

Dependence on C
Despite of having a theoretical result for a suitable C that respects Assumption 5, we must ask if we can properly determine a suitable C with practical relevance.
To do so, let us first note that the AV regularization parameter,r C as a quantity of C, is monotonely decreasing in C (Proposition 5 in the Appendix). Therefore, the suggested threshold t = 3Cr C is a product of an increasing part in C, namely 3C, and a decreasing part in C,r C . Further, we know that an increment of the regularization parameter of the graphical lasso (4) biases its values towards zero, and hence it is inappreciable to choose a threshold independently of the regularization. And now, here comes the trick: the threshold seems to automatically adopt to the chosen regularization parameter and results, even if we have some inaccuracies or some error in C, in a stable graph that does not depend decisively on the specific choice of C. Figure 5 and Figure 6 show examples of how the proposed threshold adapts on C and thresholds the AV such that we almost maximize the F 1 -score. In all cases, even though we choose different C, we can obtain approximately the same performance. To further investigate the stability of the thAV in C, we calculate the difference in F 1 -score between 2 thAV estimators F 1 (Θ t C ,Θ t C ) using different values for C (Note, that the threshold t is dependent on the particular C. However, to   avoid notational overhead, we leave out an index in t). Table 6 demonstrates the stability in C in for various different settings. Estimating a random graph, we can observe that we do not only achieve a good performance F 1 (Θ,Θ C ) for any C (Table 7), but also that the different estimates are all very similar: F 1 (Θ C ,Θ C ) is always above 0.90. The estimation of a scale-free graph on the other hand is more dependent on C, but thAV still remains superior to StARS and RIC for any choice of C.

Stability in the graphs density
The last simulation study is about the estimation results if we vary the graphs density, that is, the number of edges in relation to the graphs size. Recall, the density of a random graph, which was determined by a Gilbert graph, is controlled via the probability p of each pair of nodes to be connected. A variation in the density of a scale-free graph is not implemented in the huge package, so we leave it out for now. However, as pointed out in Barabasi (2016), we could use a general Barabasi-Albert algorithm and connect each node with m > 1 other nodes, or impose a random distribution on the number of new edges in each step to vary the graphs density.
To account for different densities of the random graph, we investigate the cases p ∈ {2/d, 3/d, 4/d}. Figure 7 demonstrates the effect of p on a graph with d = 200 nodes. It is conceivable that methods that prefer to overestimate a graph like StARS (what can be seen by a high recall but in comparison low precision , and by the visualizations of the graph) tends to perform better when the density increases. We can see in Table 8 that the F 1 -score of the thAV decreases as p increases what is not a surprising result as thAV tends to rather underestimate than overestimate the graph. We see that both other methods, StARS and RIC become worse if we violate the standard setting p = 3/d. Again, in all considered cases thAV clearly outperforms the competing estimation procedures.   Table 6: Similarity F 1 (Θ t C ,Θ t C ) for different choices of C in various settings for n and d.

Applications
Graphical model recovery plays a big role in understanding biological networks. In this Section, we apply our procedure on 2 open-source data sets and compare them with published results. Note that in practice, it can be reasonable to regulate λ in the threshold, if one is rather interested in finding many possible relations instead of finding true positives with reasonable confidence.

Recovering a Microbial Network
It is believed that the human microbiome plays a fundamental role in human health. Thus, the American Gut Project (McDonald et al., 2018) was launched to pave the way to find associations among the microbiome, but also to confirm associations between the microbiome and other aspects of human health, like psychiatric stability. Furthermore, it is believed that the bacterial diversity is of great relevance for the human lifestyle or can even be used to inspect for malnutrition. In this case, we apply the thAV to recover the microbial network to enhance the understanding of the roles and the relations between the microbes. Since microbial datasets come with some technical problems, it is vital to preprocess the data. We refer to Kurtz et al. (2015); Yoon et al. (2019) for details about the problems, but also about a suitable preprocessing routine for microbial data. We apply the procedure on the amgut2.filt.phy data, which is included in the SpiecEasi (Kurtz et al., 2015) R package. The dataset measures the abundance of microbial operational taxonomic units (OTUs) and consists of n = 296 samples and d = 138 different OTUs (Kurtz et al., 2015). The proposed method in Kurtz et al. (2015), SPIEC-EASI, applies a clr-transformation (Aitchison, 1981) and uses a MB-method based on StARS regularization parameter calibration to infer the graph. Figure 8 shows the results of the SPIEC-EASI and the thAV. We observe that thAV returns a much sparser graph than the procedure in Kurtz et al. (2015). More precisely, the thAV has 112 edges, the graph based on the MB-method has 243 edges, and the graph that consists of similar edges, that is, of edges that are in Figure 8(a) but also in Figure 8(b), consists of 97 edges. The fact that the graph of similar edges has almost as many edges as the thAV graph suggests once more that the thAV most likely achieves a high precision. Based on our observations in Section 4, we also applied the thAV approach using λ = 2, which result is depicted in Figure 9.

Recovering a Gene Network
In pharmacology, riboflavin is industrially produced using diverse microorganisms. In the age of technological revolutions in genomics, biologists found a genetically modified version of Bacillus subtilis that has an increased riboflavin production. Being able to fully understand the bacterias genome, biologists promise to further optimize the riboflavin production. We apply the thAV on the riboflavin data set, which is provided by the DSM in Switzerland and contains n = 71 samples and d = 4088 gen expressions. The R package hdi (Dezeure et al., 2015) provides this data set in its implementation.
We compare our results with these in Bühlmann et al. (2014), who analyzed the same data set using a MB approach. Similar to Bühlmann et al. (2014), we shrunk the data set by only considering the 100 genes with the highest empirical variance and scaled the data using the nonparanormal transformation (Liu et al., 2009). We observe that the thAV (Figure 10(a)) returns a sparser graph with more cluster-like structures than the MB estimate ( Figure 10(b)). This does not only increase the interpretability of the graph but also imposes some tight connections between several genes within these clusters. In regard of Corollary 6 and our simulation studies in Section 4, we suspect a the thAV to have a high true-positive rate.

Discussion
In this paper, we extend the ideas of Chichignoud et al. (2014) to calibrating the graphical lasso, the most popular estimator for Gaussian graphical models. The resulting calibration method, that we call thresholded adaptive validation (thAV), outmatches existing calibration methods empirically and provides a simple and fast implementation and a finite sample guarantee. The latter allows to derive a corollary suggesting how to optimally choose the threshold-parameter to obtain finite sample results on the recovery performance.
Using thAV with the proposed threshold in the empirical simulations, we see that the method convinces with a good F 1 -score and excellent results in precision. Sometimes the high precision comes at the cost of a too sparse estimate. If the user is interested in a denser graph, one could contemplate to lower the threshold as also suggested by our theory and experiments.
Despite our focus on graphical lasso in this paper, one could readily adapt our scheme to other graph estimation methods, like TIGER (Liu and Wang, 2012) or methods which aim to recover a specific graph topology. For instance, there are methods aiming to recover the support of a particularly scale-free graph, see Liu and Ihler (2011) and references 15  Figure 7(a), 7(b), 7(c), respectively to generate the adjacency matrix.
therein. Extending the theory for such estimators, one might eventually be able to improve the recovery performance of scale-free graphs by apply the thresholded adaptive validation technique for such kind of estimators. Furthermore, one could replace ∞,off in Assumption 2 by any loss function that satisfies the triangle inequality and is symmetric. This includes, for example, the Wasserstein distance (Kolouri et al., 2017). Bounds in Wasserstein distance could be, for instance, useful for PCA, where one tries to recover the eigenstructure of Θ, or Linear Discriminant Analysis (LDA) (Fisher and R. A., 1936), where one is interested in Θ itself rather than in its conditional dependence structure. Both methods are known to suffer in a high-dimensional setting, PCA does so because the eigenstructure ofΘ emp is a poor estimator for the true eigenstructure (Wainwright, 2019), and LDA does so becauseΘ emp is generally a bad estimator of Θ. The details of the estimation of Θ using another loss could be the topic of further research.
In addition, it would be interesting to investigate the robustness of thAV against model miss-assumptions, for instance if z ∼ N d (0 d , Σ). Especially in view of Section 5.1, one could do some research on the performance of thAV based on synthetically generated data for microbial networks (Kurtz et al., 2015;Yoon et al., 2019).   Figure 8: Based on the clr-transformed American Gut Dataset provided by the R package SpiecEasi, we applied the proposed method in Kurtz et al. (2015) (Figure 8(a)) and the thAV (Figure 8(b)) on the non-paranormal transformed data (Liu et al., 2009) to relax the normality assumption. To avoid too large graphics, we excluded isolated vertices. The color of a node implies the biological cluster (taxonomic rank 3) of each OTU, which is also provided by the SpiecEasi package. 18 Figure 9: Based on the clr-transformed American Gut Dataset provided by the R package SpiecEasi, we applied the thAV using threshold t := 2Cr on the non-paranormal transformed data (Liu et al., 2009) to relax the normality assumption. To avoid too large graphics, we excluded isolated vertices. The color of a node implies the biological cluster (taxonomic rank 3) of each OTU, which is also provided by the SpiecEasi package.
(a) (b) Figure 10: Estimated graphs using the riboflavin data. Figure 10(a) depicts the graph obtained by thAV. Figure 10(b) is obtained by applying a MB method where each regularization parameter is calibrated using StARS (Bühlmann et al., 2014) 20

A.1.2 Incoherence and bounded norm of the Primal
Now we are ready to impose a well behaviour of Θ that ensures the PDW to succeed and implies the validity of Assumption 2. Assumption 8 (Θ well-behaved). There exists an α ∈ (0, 1] such that −1 for all r ∈ R, where deg(Θ) denotes the maximum degree of the corresponding conditional dependency graph induced by Θ.
The first assumption is widely used and is often referred to as incoherence or irrepresentability condition. To shed some light onto this condition Ravikumar et al. (2008) consider for i, j ∈ V the centered random variable y ij := z i z j − E[z i z j ] , which can be interpreted as the interaction between components z i and z j . One can show that the incoherence condition requires y ij for (i, j) ∈ E to be not highly correlated with y lm with (l, m) ∈ E. This interpretation goes in line with the incoherence condition for the lasso, see Ravikumar et al. (2008) and references therein. Hence, the incoherence assumption bounds the correlation between irrelevant interactions (y ij with (i, j) ∈ E) and relevant interactions (y ij with (i, j) ∈ E). The second assumption states that the constrained graphical lasso,Θ r , is good enough in the sense that it satisfies the ∞ -bound. Note, that the bound scales with the maximal degree of the graph. Hence, the higher the degree of the true graph, the better the primal solution must be to satisfy Assumption 8. This indicates a reason why scale-free graphs, which have in general a larger maximal degree than random graphs, tend to be harder to be estimated by the graphical lasso, as also indicated by the results of our experiments.