Conditional Matrix Flows for Gaussian Graphical Models
Authors: Marcello Massimo Negri
Fabricio Arend Torres
Volker Roth
Venue: NeurIPS 2023 poster (2023)
Keywords: normalizing flow
variational inference
graphical lasso
gaussian graphical model
bayesian inference
TL;DR: General framework for variational inference with matrix-variate Normalizing Flow in Gaussian Graphical Models

Abstract: Studying conditional independence among many variables with few observations is a challenging task.
Gaussian Graphical Models (GGMs) tackle this problem by encouraging sparsity in the precision matrix through $l_q$ regularization with $q\leq1$.
However, most GMMs rely on the $l_1$ norm because the objective is highly non-convex for sub-$l_1$ pseudo-norms.
In the frequentist formulation, the $l_1$ norm relaxation provides the solution path as a function of the shrinkage parameter $\lambda$.
In the Bayesian formulation, sparsity is instead encouraged through a Laplace prior, but posterior inference for different $\lambda$ requires repeated runs of expensive Gibbs samplers.
Here we propose a general framework for variational inference with matrix-variate Normalizing Flow in GGMs, which unifies the benefits of frequentist and Bayesian frameworks.
As a key improvement on previous work, we train with one flow a continuum of sparse regression models jointly for all regularization parameters $\lambda$ and all $l_q$ norms, including non-convex sub-$l_1$ pseudo-norms.
Within one model we thus have access to (i) the evolution of the posterior for any $\lambda$ and any $l_q$ (pseudo-) norm, (ii) the marginal log-likelihood for model selection, and (iii) the frequentist solution paths through simulated annealing in the MAP limit.

Conditional Matrix Flows for Gaussian Graphical Models
Marcello Massimo Negri University of Basel
Fabricio Arend Torres
University of Basel
fabricio. fabricio. wwt..aartttttreeessaaa....
Volker Roth
University of Basel
volker www.rroottsss..oommm
Abstract
Studying conditional independence among many variables with few observations is a challenging task. Gaussian Graphical Models (GGMs) tackle this problem by encouraging sparsity in the precision matrix through l regularization with q S 1. However, most GMMs rely on the norm because the objective is highly non-convex for sub-l, pseudo-norms. In the frequentist formulation, the 1 norm relaxation provides the solution path as a function of the shrinkage parameter A. the Bayesian formulation, sparsity is instead encouraged through Laplace prior, but posterior inference for different À requires repeated runs of expensive Gibbs samplers. Here propose general framework for variational inference with matrix-variats Normalizing Flow GGMs, which unifies the benefits frequentist and Bayesian frameworks. As key improvement on previous work, we train with one flow a continuum parameters À and all lg norms, including non-convex sub-l1 pseudo-norms. Within one model we thus have access to (i) the evolution of the posterior for any À and any l (pseudo-) norm, (ii) the marginal log-likelihood for model selection, and (iii) the frequentist solution paths through simulated annealing in the MAP limit.
sparse regression models jointly for all regularization
Introduction
Estimating complex relationships between random variables central problem in science. When only few observations are available, the problem becomes even more challenging. Examples include functional connectivity in fMRI data (Smith al. 2013], networks of interactions from microarray data [Castelo and Roverato 2006), or correlation patterns in longitudinal studies (Diggle 2002). In such cases, it is particularly challening to infer the conditional independence among random variables. Graphical models are commonly used to represent such an independence structure in the form of network graphs. In this work, we focus specifically on Gaussian Graphical Models (GGMs). GGMs assume observations X E [Rdd to be generated from multivariate Gaussian distribution X N (ja, Σ) with X,, being realizations of X' for E (I,..., n). Here μa E Rdd denotes the mean and Σ Rdxd the covariance matrix, which symmetric positive definite Σ > 0. GGMs provide a simple interpretation of conditional independence through the precision matrix Ω Σ whenever Σ is non-singular. Specifically, 12 = implies that the pair variables (i.j) is conditionally
independent given all remaining variables. That is to say, there is no edge between nodes i and ý in the underlying undirected graph.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Penalized likelihood formulation Given the centered observations X and the associated sample covariance matrix S x"XX, X. we are interested in reconstructing the precision matrix 52. hence the underlying graph structure. This is particularly challenging when d > n. because the sample covariance matrix becomes singular and the precision matrix can no longer be obtained by simply inverting the MLE of the covariance matrix. One way to overcome this problem is to consider saclihood--pnaliiee models that encourage sparsity in the precision matrix. In other words, we trade off the likelihood with the number of zeros in S2:
Ω arg max (log det Ω Tr(().s) 211211a) 12>D
where log det 52 - Tr( SS2) is the log- log-likelihood term and ||2||o Σiej luw # 0J is the tn norm, which counts the number of non-zero dff-diagonal elements of S2. The trade-off between the two terms is controlled by the parameter À Z 0. Note that the optimization must be performed over the space of symmetrie positive definite matrices Z > 0. In practice, the objective in Eq. TT is highly non-convex and cannot be optimized easily. Even in the simpler linear regression setting. lo regularization can be computed exactly only for few features [Hastic et al. 2015]. Similar problems arise for all non-convex sub-11 pseudo-norms.
Lasso relaxation: Frequentist and Bayesian approaches In order to simplify Eq. (I. most approaches replace the highly non-convex in norm with the closest convex norm, the T norm: ||a||] = Exej [tyl- Meinshausen and Bühlmann (2066) first proposed using l regularization in the
Graphical Lasso model, which sparked the development of several likelihood penalized algorithms
2007 Yuan and Lin 2007 Banerjee et al.] 2007]. These frequentist approaches (Friedman et al.
make it possible to elegantly compute the solution path as function of the shrinkage parameter À
[Mazumder and Hastie 2012). The -reelaxed problem can also be easily extended to a Bayesian
formulation. Waang] [2012] assumed a Laplace prior on the off-diagonal elements of the precision matrix and showed that the frequentist solution can be recovered as the maximum a posteriori estimate (MAP). In the Bayesian framework, it is possible to explore the full posterior distribution, to formulate the posterior predictive, and to use the marginal likelihood for model selection. But in sigg--iimnniooal settings, MCMC samplers suffer from poor mixing behaviour, which becomes increasingly difficult to diagnose higher dimensions. Furthermore, different hyperparameters values, e.g., X, require independent expensive Markov chains. In the literature, alternative priors have also been proposed (Li et al. 2019]. However, the dependence of Gibbs samplers on tractable posterior conditionals typically restricts the choice of priors significantly.
Variational Inference Variational inference approaches [Bei et al. 2017] approximate intractable densities with tractable parameterized distributions. Compared to MCMC samplers, variational infer- ence turns a sampling problem into an optimization one, which generally requires less computational time. Let be the observed variables and z the unobserved ones. Given family Q (go(2) of parameterized distributions, variational inference involves finding the distribution qo- (z) that best approximates the posterior p( (z).. The distance between the two distributions is usually measured as the Kullback-Leibler divergence, which defines the following optimization problem:
(z(2) log
go (z) = arg min (L(g = arg min Ez 8E0
p( (x)
The goodness of go- (2) the expressive power of the variational family Q and the tractability of the optimization in Eq. (2) Common approaches rely on the so-called mean field approximation, which assumes mutual independence among variables: (z)z)= = II qo, (2.). But this is not viable in GGMs as we intend to model precisely the dependence structure. Several approaches have been proposed for Bayesian Lasso Alves et al.] [202]] and for its group-sparse variant [Babacan et al. 2014) Yet, to the best of our knowledge, variational approaches have not been studied in the context of Bayesian GGMs.
approximation of the posterior p(z(æ:) creates trade-off between
Contribution We present a unified approach to sparse regression for the whole family of all 1, (pseudo-) norms with 0 < < xo, which allows for a fully Bayesian and frequentist-type interpretation. Specifically, we propose a very general framework for variational inference in Bayesian GGMs through a Normalizing Flow defined over the space of symmetric positive definite matrices. By conditioning the flow À and on q. we simultancously train a continuum of sparse regression

---

models for all choices of shrinkage parameters 0 < À < Amax and all ig (pseudo-) norms with 0 <<g 9 < 9max: including the highly non-convex sub 4 pseudo-norms, We use as prior the generalized Normal distribution, which for q recovers the Laplace prior and the Lasso regularization in the MAP limit. On the one hand, our approach inherits the advantages of the Bayesian framework while altogether avoiding problems arising from Gibbs sampling strategies. On the other hand, we can still recover the frequentist solution path in the MAP limit by training through simulated annealing.
In summary, our main contributions are the following:
1. We propose general framework for variational inference in GGMs through normalizing flow defined directly over the space of symmetric positive definite matrices. We condition such flow on the shrinkage parameter À and on g >> 0 to model L, (pseudo-) norms.
2. We combine the advantages of the frequentist and Bayesian frameworks: in single model we have access to posterior inference and to the marginal likelihood as function and g. We can further recover the frequentist solution path as the MAP by annealing the system.
3. To the best of our knowledge, we are the first to propose a unified framework for sub-l1 pseudo-norms that does not require surrogate penalties and that enables consistent explo- ration of the full solution paths in terms of À and q.
2 Related Work
In this section, we provide an overview of the Bayesian Graphical Lasso and illustrate the limitations of current inference methods. We then outline the limitations of Lasso regularization and briefly review existing approaches for sparsity with sub-l2 pseudo-norms. Lastly, as an alternative MCMC approaches, we briefly review variational inference with normalizing flows.
Bayesian Graphical Lasso Similarly the Bayesian Lasso Park and Casella 2008 Wang 2012
provided a Bayesian interpretation of the Graphical Lasso (BGL) for posterior inference. Specifically, they showed that the 1 relleee optimization of Eq. (I) equivalent to the maximum a posteriori estimate (MAP) of the model defined by p(. S, A) x p(S(() p(1 A), with
P(S( x Wa(n,
((SI))) x Wald + 11 II DE > 0).
where the indicator function [[ > 0) imposes positive definiteness on 52. The likelihood term reflects
that the sample covariance matrix S = X X is distributed according to the Wishart distribution exp Tr(-An.s). The prior is instead composed of
Wa(n, 1-1) x d2t dct(S) two terms. The Wishart term Wa(d + 1, Ala) imposes symmetry and positive definiteness on Ω. In addition, the double exponential (or Laplace) prior DE(wa)) (A) = 2/2 (X encourages sparsity on the off-diagonal elements of f.
Inference in the BGL, and more generally for Bayesian Lassos, is performed through Gibbs samplers. However, these MCMC strategies are computationally expensive and can suffer from high rejection rates, especially in high dimensions [Mohammadi and Wil 2015]. Furthermore, to recover the frequentist solution path in the MAP limit, the Markov Chain must be restarted for each À value. Its derivation also requires to expand the Laplace prior as an infinite mixture of Gaussians [Andrews) and Mallows, 1974 West 1987] and to define an ad hoc mixing density. Lastly, the obtained Gibbs sampler does not generalize other priors.
Sparsity with sub-l1 pseudo-norms The to norm is particularly suited to enforce sparsity because directly penalizes the non-zero entries. In other words, it is equivalent to best subset selection of the variables. Convex relaxation with the l norm makes the problem tractable, but comes at the cost of encouraging shrinkage also on the retained variables [Hastie et al.) 2015). This motivated interest in sub-l, pseudo-normy which in the limit reduce to the original to norm formulation. Due to their non-convexity and combinatorial complexity, multiple surrogate objectives have been proposed. Notably, with the proposed approach we do not have to resort to alternative objectives or penalties, as we can directly and exactly enforce sub-l, pseudo-norm regularization with suitable prior. In the
context of linear regression, several algorithms have been proposed. The popular SCAD Fan and Li
3

---

2001] guarantees unbiasedness, sparsity, and continuity while reducing the overshrinking behaviour. Zhang 2010] proposed the nearly unbiased MC+ method, which consists of concave penalty and a selection algorithm that bridges the gap between 1 and l-- In the literature, other penalties that approximate the la penalty have been proposed, such as the SICA penalty [Lv and Fan 2009], the log-penalty [Mazumder et al., [01T], the seamless-lo penalty [Dicker et al. 2013] and the atan penalty [Wang and Zhu [016]. In the Bayesian framework, Ishwaran and Rao 20055 used a rescaled spike
and slab prior to encourage variable selection and drew connections to Ridge regularization.
Normalizing Flows for variational inference Normalizing Flows (NFs) are flexible models that can perform accurate density estimation For this reason, NFs represent a very attractive solution to perform efficient and accurate variational inference [Rezende and Mohamed, 2013 van den Berg et al. 2019). Suppose X' is a continuous d-dimensional random variable with unknown distribution Pa. and let 2 be any continuous di d-dimensional random variable with some known base distribution p2. The key idea of NFs is to construct a diffeomorphism T supp(X() U supp(Z). i.e., a differentiable bijection, in order to rewrite px through the change of variable formula as
(x(=) Pz (T((T)))) det (x))
(4)
where JT is the Jacobian of the transformation T. Assuming that we have successfully learned the transformation T, we can then evaluate the target density through Eq. (4) We can further sample from the target distribution by simply transforming samples of the base distribution through T. namely æ = T(z) with z being realizations of pz. In practice, the crucial part of designing NFs is to construct arbitrarily complicated bijections. To do NFs exploit that compositions T = o of bijections } remain bijections. The determinant of the resulting transformation decomposes into the determinants of each T,
det JT(x)
II del JT. (uC-1) with u:-1 = TT( ...... Ti- 1;=1(u;-2)
(5)
Therefore, by composing computationally tractable non-linear bijections, it is possible to define arbitrarily expressive bijections and hence transform the base distribution pz into arbitrarily complicated distributions px. Among others, |Huang et al. [2018] proved that autoregressive flows are universal density approximators of continuous random variables. For a comprehensive review of NFs and their different architectures, see Papamakarios et al. [2021] and Kod/eee et al. [2020). As first proposed by Atanov et al. [2020). is also possible to condition the transformation 7 on some parameter E * The resulting Conditional NF (Kobyzev et al. 2020) models with one flow the family of conditional distributions px (x (c)c for continuous values of c.
3 Proposed approach
In this section, we describe how to infer the posterior of Ω in GGMs with conditional NFs and how to do as function of the shrinkage parameter À and of the l (pseudo-) norm regularization, all within a single model. We illustrate how to define the conditional flow directly over the space of positive definite matrices and argue why training and posterior inference is particularly efficient. Furthermore, we show that NFs provide direct access to the marginal log- ligliiked for model selection and show how to recover the frequentist solution path by training through simulated annealing. We term the resulting model Conditional Matrix Flow (CMF).
3.1 Generalized Normal distribution for lg (pseudo-) norms
Unlike previous Bayesian approaches to GGMs, we can flexibly specify any prior (and likelihood) in Eq. without worrying about the existence of suitable Gibbs sampler. We exploit such flexibility to extend the model beyond the standard I relaxation. Wang [2012) showed that in the MAP limit 4 regularization is recovered with Laplace prior on the off-diagonals of J2. We generalize this idea to any lg (pseudo-) norm with the generalized Normal distribution as prior, which reduces to Laplace and Normal distributions for q and g 2, respectively. Its probability density is defined as:
B exp {
x 8
f(xla, 3) = (0((1/3)
(6)

---

where r() is the Gamma function, B > 0 are the scale and shape parameters, respectively. Here, we assumed the distribution to be centered around zero, i.c. μl == 0. In order link the prior the l norm, we rename its parameters as À = a- > and q = B. Our full model can now be conditioned through the prior on both À and and is defined as p(22) S, X, g) x p(S[2) p(2) X q) with
((S(S) Sa(
g) & Wa(d+1,,14 II i<j
9A 2T(1/)
exp(- [|[[ 0
(7)
If now consider the log- peg-robabilityy and drop the normalization constant, which does not affect the KL optimization, we recover the l (pseudo-) norm. For = 1 the prior on the off-diagonals reduces to the Laplace distribution and we recover the Bayesian Graphical Lasso in Eq. (3) For
= 2 we recover instead the normal distribution, i.c., Ridge regularization. Overall, when employing effectively model solutions corresponding to lg (pseudo-) norm
the distribution in Eq. (7), regularization for > 0, including the non-convex sub-11 aseudo-norms
3.2 Variational inference with conditional flows for GGMs
As want to model the posterior SS over the precision matrix Z design flows directly over the space of symmetric positive definite matrices. Specifically, we train a NF conditioned on x E [λI, X2] and on q E [41, qa] with X1, Ag. 21. > in order to study the evolution of the posterior as function of À and From Eq. (4) we see that the resulting flow To(>.g) implicitly
defines the probability 9o(x.g) (S2) = Pbase (To(..g) (S2)) det (To().g) (52)), (S2)), which we want to train to approximate the posterior p(52) S, X, 4) x (()) p(2]), q). We minimize the KL divergence Eq.
(2) and approximate the expectation with Monte Carlo samples:
C(0; KL (1o(A.) (S))
M
(012.a) (S2;)
(8)
M Σ log ((S),,, p(
+ (N)(),))
Note that this is particularly efficient in NFs because it only requires sampling from the base distribution Phose( ( (To(A.g) (I2), which is computationally cheap. We further need to evaluate the unnormalized posterior SS X, q) & ((S p(S2) q), i.c. the product of the likelihood and the prior, which is also trivial. Relevantly, in Eq. 88) we do not need compute the normalization constant Nc(A, = log p( S X, 2) because it does not influence the optimization.
Unlike standard Bayesian approaches, the proposed conditional flow provides access directly the marginal log-likelihood Nc(A, q) log p( X, g) as function of and without further calcula- tions. This is particularly interesting because the marginal ligeliill extremely expensive to compute with classical approaches. If we assume that the flow is expressive enough and that the optimization reached the global minimum, then KL (9o:)...)) (S2))I((S])S X, 4)) ss In this case, the loss function at the end of training gives us access directly to the negative marginal log- liglliieidoood for continuous values of (and q). We can then perform model selection by simply choosing X = max log (SS)) 4).
3.3 Conditional Matrix Flow
We now describe the proposed Conditional Matrix Flow (CMF) and how we define over the space of symmetric positive definite matrices by construction. We exploit the well-known Cholesky decomposition [Horn and Johnson 1985). which states that any symmetric positive definite matrix s
M Rdxdd can be decomposed as M LLT, where L E R is a lower triangular matrix. The decomposition is unique if L has positive diagonal elements. As a result, the symmetric positive definite M can be fully and uniquely identified by d(d + 1)/2 numbers, d of which are constrained to be positive. We can then use d(d + 1)/2-dimensional normalizing flow to model probability densities over symmetric- positive matrices. We now illustrate how to transform the initial vector [d()+1) into a symmetric positive definite matrix Z E dRdxd To the best of our knowledge, this is the first attempt to define a flow T(2) on the space of symmetric positive definite matrices by construction.
5

---

Hyper- NeterrN 2,
Sum-of -iggmoods
Fill- Triangular
Positive-Diagonal
Cholesky-Produc
did+1)
#id+1) 2
d
m
x(#)
p(D)A,))
L(=
Figure 1: Architecture of the proposed Conditional Matrix Flow (CMF) model. Light grey is used to denote unconstrained values, dark grey for positive values, and white for zeros.
A high-level visualization of the proposed architecture is shown in Figure ! The first n layers of the (x)) network consist of arbitrary NF transformations and map the initial vector to =x.g We propose to use class of flexible transformations that we called Sum-of-Sigmoids which worked particularly
well for our applications. We provide further insight on the implementation Appendix [A.]] Furthermore, condition the flow on À and q via hypemeeworr. The hypemnetwrr takes À and g (n) as inputs and returns the parameters of the Matrix Flow. We transform Zx.q into a symmetric positive
definite matrix in three steps. First, the vector is raveled into lower triangular matrix L E [Rdxd which we call Fill- Fill-Trinnular. This transformation has unit Jacobian determinant because it just reshapes the vector into matrix. In second step, we bijectively map the diagonal of the resulting lower triangular matrix to positive values with a softplus activation Softplus(2 (a) log(1 + exp x). This transformation, which we call Jacobian (log) determinant. Lastly, we compute the Cholesky product Chol L) + LL", which again has an inexpensive Jacobian (log) determinant [Gupta and Nagar 1999
acts element-wise and hence admits cheap
det JChol (L) 2d
II (L.) d-141 1=1
(9)
which depends only on the d diagonal elements of L. Note that by enforcing the triangular matrix to be positive, we make sure that the Cholesky decomposition is unique, and hence bijection.
We open-source the implementation the conditional bijective layers We based our library, which was first used in Torres et al. [2023] on the high-level structure provided by the nffows library [Durkan
al. 2020
3.4 Training through simulated annealing
With our CMF, we can perform posterior inference as a function of À and and to perform model selection. These are key advantages of the Bayesian perspective. In addition, frequentist solution path as function of without modifying the model. Specifically, we employ optimization through Simulated Annealing (Kirkpatrick al. [983] approximately sample from the global maxima of the distribution, which for posteriors is the maximum a posteriori estimate (MAP). The idea comes from statistical mechanics, where slowly cooling processes used to study the ground (optimal) state of the system. When applied to more general optimization tasks, simulated annealing consists accepting iterative improvements they lead lower cost function. Meanwhile, the temperature is slowly decreased from the initial value To to T, s 0, where the system frozen and no further changes occur. In our setting, this corresponds to introducing an artificial temperature T; for the target posterior in Eq. 7
also recover the
T p(2 S.A.))
(10)
where T; is the temperature at the i-th iteration. If the initial temperature is high enough, ( will likely be very flat, allowing for a better exploration of the support of the distribution. As the temperature decreases, the distribution becomes more peaked. In the limit T, + 0 the distribution
FlowConductor: (Conditional) Normalizing Flows and bijective Layers for Pytorch http://ggtttt.tooooo//ccccceetet

---

p()(11 should concentrate on the global maximum, hence on the MAP solution. Geman and Geman 1984 formally showed that convergence the set of global minima is guaranteed for logarithmic cooling schedules. Unfortunately, logarithmic cooling schedules are not viable in practice, so alternative schemes have been explored [Abramson et al. 1999). In this paper, we use the popular geometric cooling schedule [Andrieu and Doucet 2000, Yuan et al., 2004), which works well in practice: T; = Tyat/" for some a > 0 and i E {O..... n). In practice, we train the objective in Eq. (8) while slowly cooling down the system, i.e. decreasing T, in Eq. (IO), until we reach a sufficiently low temperature Tū By saving the model at T, = 1 and Tn we have access to both the Bayesian and frequentist solutions, respectively.
3.5 Limitations
The proposed approach provides general framework for posterior inference in GGMs and generalizes to a large family of likelihood and priors, which requires an efficient evaluation of their analytic expression. In some cases, e.g. for posterior predictive distributions, this would require an additional integration step, which could become computationally expensive. The proposed approach is also limited by the expressive power of the bijective layers of the flow. Even though sute-of-the-an layers are extremely powerful in modeling high- digh-dimensiona distributions, each layer might still be limited in terms of the number of modes that can be modeled [Liao and He, 2021]. Lastly, the proposed model assumes that we can model a family of posterior distributions as function of the conditioning parameters g and X, which ultimately depends on the flexibility of the hypernetwork used.
4 Experiments
In this section, we showcase the effectiveness of the proposed CMF first on artificial data and then on real application. In particular, study the evolution of the variational posterior a function of X and q. We then perform model selection on À through marginal likelihood maximization. We further illustrate the effect of training through simulated annealing and show that we recover the frequentist solution path through the MAP. Lastly, we show that the proposed method can be readily applied to real data in higher-dimennional settings. Results show that sub-/1 pseudo-norms provide sparser solutions and contrast the well-known overshrinking effect of Lasso relaxation. We provide the code reproduce the experiments at https://github. http:/aww.ccoo..aaaaelleeee/////
4.1 Toy example: synthetic data
We illustrate how the proposed CMF works on artificially generated sparse precision matrices. The data generation process consists of sampling a sparse precision matrix [Pedregosa et al. 2011] with d features and sparsity level and then generating Gaussian samples accordingly. For illustrative purposes, we show results for one precision matrix with a = 90% and n = d = 15 (n = d ensures model for 10'000 epochs through
the invertibility of the empirical covariance). We trained simulated annealing with an initial temperature of To = 5 to TT = 0.01 and performed 100 geometric cooling steps with T; = S/" for a = T../To-
Figure 2 shows the effect of simulated annealing on the posterior. At temperature T 1 the conditional flow approximates the true posterior in Eq. 7 which for q = coincides with the Bayesian Graphical Lasso (BGL) model. As we decrease the temperature, the (unnormalized) target posterior Eq. (TOp becomes more peaked and we observe shrinking credible intervals. At the final temperature Tn. the distribution converges to the MAP and we recover the frequentist path. As shown in Figure 2 (right), our model accurately reconstructs the solution path over À (MSE 0.052). Furthermore, we perform model selection on À for the model in Eq. 7 which is obtained at temperature T = 1. In the Appendix A..4 in Figure 7 we show the (approximate) marginal function of À and the resulting optimal XCMP = 3.52 that maximizes it. We
log- likelihhood compare the result with the frequentist estimate obtained through MLE with 5-fold cross validation GLaaso 3.36. Again, the model agrees with the frequentist solution.
One of the most significant drawbacks of relaxing the lo-morm regularization with I is that sparsity is achieved through shrinkage. This means that the selected relevant features are over-shrunk. With the proposed CMF we show that by exploring sub-l, pseudo-norms we can reduce the shrinkage effect and virtually overcome it in the q + 0 limit. In Figure 3 (left) we show 95% posterior
7

---

1.0
0.5
10
10a
A
Figure 2: 75% posterior credible intervals function of for q 1. Left: at T 1.0 the CMF reduces to the BGL. Right: at T = 0.01 the CMF reduces to the frequentist solution path (dashed).
CMF (g=0.25)
CM Iq=0.50)
lQ=0.75)
CMF (q-1.00)
BGL
Glasso
0.3
1.0
3.2
10.0
10P
10A
Figure 3: Left: 95% posterior credibility intervals of the proposed CMF for one entry of the precision matrix as function the pseudo-norm = (1, 0.75, 0.5, 0.25]. Results are compared with BGL and the frequentist solution (GLasso). Right: MAP estimate of the CMF posterior function of À for 4 = 0.5 (T 0.01). The dashed line the frequentist solution path for =
credibility intervals for entry of the reconstructed precision matrix. We observe that, the l pseudo-norm gets closer to l1. the posterior median is progressively shrunk towards zero. With the proposed CMF sub-l, pseudo-norms. In particular, note that for À = 0.3 the posterior 95% credibility interval does not include the value 0 for g (0.25, 0.5, 0.75), (right) we show the posterior solution path of the CMF for 0.5 in the MAP limit (T' = 0.01). Compared the solution path, sub-/1 solution paths show less shrinkage effect (higher values for the selected features) and steepest decrease towards zero, This behaviour becomes more evident as decreases, as we show more extensively in the Appendix A.4 Figure for {1, 0.75, 0.5, 0.25). These results support the interest in sub-l1 pseudo-norms norms,,,, and in lo norm in the limit.
prevent this shrinkage effect by exploring the posterior solution path for
opposed to the BGL (or q 1). In Figure
4.2 Edge recovery with sub-l, pseudo-norms
We now study the behaviour of the proposed CMF in terms of Fl score for edge recovery as a function of the number of samples. We show that the proposed CMF outperforms competing methods, especially in the low sample regime. We compare the CMF with the BGL, which is the only alternative Bayesian model. Then, compare the CMF in the q 0 limit against frequentist approaches that use approximate In norm penalties: Atan ("atan") (Wang and Zhu 2016), Seamless Io ("selo") Diikker et al. 2013), log-penalty ("log") [Mazumder al. 201 I]. SICA ("sica"" [Ly and Fan 2009
use ground truth precision matrices of dimension 30 and generate n Gaussian samples accordingly. We show results for the relevant < d regime and around d, namely for samples

---

{15, 25, 35, 45). For Bayesian approaches (the proposed CMF and the BGL) we draw 1000 samples from the (approximate) posteriors and consider 90% credibility intervals. All results are averaged over the 10 precision matrices. The proposed CMF is trained for 5000 epochs to a final temperature = 1, which corresponds to the Bayesian model in Eq. 7 For the BGL. we run the Gibbs sampler for 4000 iterations with burn-in of 1000 and keep every fourth sample. The frequentist approaches were run with the specific syper-parameters suggested their original papers.
0.7
0.7
0.6
0.5
0.5
(p=0.255
CMF (p-0.25)
CMF (p=0.5)
stan
CMF (p-0.75)
CMF (j=1)
selo
BGL
0.3
0.3
15
25
35
45
15
35 # samples
45
# samples
Figure F1 score for edge recovery. Left: the proposed CMF for g = {0.25, 0.5, 0.75, 1} against the Bayesian Graphical I Lasso (BGL.). Right: the proposed CMF for 0.25 against various frequentist approaches with penalties that approximate the to norm. Results are averaged over 10 precision matrices with d 30.
As the proposed method is inherently Bayesian we first compare it against the BGL. We show that for q = obtain results compatible with the BGL.. Furthermore, with the proposed method can additionally infer the posterior for q < 1. The results Figure 4 (left) show that in the low sample regime (n < d) sub-l1 pseudo-norms are beneficial and result in higher F1 score. The effect stronger as - 0, while in the n > regime sub-l, pseudo-norms do not provide a significant advantage. We also compare the proposed CMF with classical frequentist approaches with surrogate penalties that approximate the lo norm. The results in Figure 4 (right) show that the proposed CMF with q 0.25 outperforms all competing methods across all sample regimes, but especially for very low sample sizes (n 15). Note that in the < regime frequentist algorithms require an ad hoc initialization of the precision matrix, which we provided through the Scdoit-Wol: shrinkage estimator.
4.3 Real data application
We now consider a real-world application to showcase that the proposed CMF can be easily used in high-dimensional settings. We consider setting in which are interested in studying only a subset of query variables. In particular, we show that we can avoid inferring the posterior on the full precision matrix, which is expensive, by simply redefining the target posterior, i.e. the prior and likelihood. We consider a colorectal cancer dataset (Sheffer et al.) 2009), which contains measurements variables together with 312 gene measurements from biopsies for 260 cancer patients. As the dataset contains many missing values, we drop the p53 mutation status clinical variable and only consider 190 fully measured patients. We study the connections between the s 6 clinical variables and the f 312 gene expression measurements. Like Kaulmann et al. [2015] we consider the partition
clinical
Q Ω
Ω12 222 1
8
S
S
Ω
S = ST.
S22
(11)
where < t. Instead of inferring the expensive (s + f) x (s + f) precision matrix Ω, we study the Q2 and 2212 sub-matrices, which overall require only s x (s + t) dimensions. Following Kaufmann et al. (2055). we show in Appendix A..pendixAA.21thie that we can infer the posterior p(l11, S12/S, A) independently of the large Ω22211 1. But contrast to Kaufmann et al. 2015), we (i) can enforce the correct double-exponential prior on Q1i we (ii) do not need to invert the st x sl matrix, which takes (ss13) operations, and we (iii) do not require Gibbs sampler, which can still suffer from poor mixing behaviour. In our framework, we only need define a different posterior as the target

---

distribution. In practice, we define the flow jointly over the sub-muttren motr-matiiess U211 which must be positive definite, and 12 which unconstrained. It is then sufficient define the bijective layers on the s(s + 1)/2 plus x / dimensions and to perform the Cholesky product only on the dimensions encoding S211- Note that the determinant of the Jacobian of the full transformation still the same as in Eq. (9) because the x dimensions encoding L2 are just raveled into a matrix.
Figure 5: Inferred network structure with the proposed CMF for q 1.0 (left) and 0.6 (right). Positive edges are shown in black and negative ones in grey. Clinical variables are highlighted in red.
In Figure 5 we show the inferred network structure with the proposed CMF for 1.0 (left) and for q = 0.6 (right). Significant edges are obtained by considering 80% credible intervals on the posterior of Pun and 12 The clinical variables in the data are age, sex, cancer group stage (GS), and TNM classification for colorectal cancer, which measures its size (T), whether it spreads to lymph nodes (N) and if metastases develop (M). Note that, in contrast to [Kaufmann et al.. 2015].
infer the posterior on the 121 block as well, which allows to study the conditional independence structure among clinical variables. For instance, results suggest strong connections between the cancer group stage (GS) and each variable of the TNM classification (T, N, M). More interestingly, with the proposed CMF we can infer the posterior for different sub-/1 pseudo-norms, As expected, Figure shows that with lower values we obtain sparser solutions. Additional quantitative comparisons for q = [1.0, 0.9, 0.8, 0.7, 0.6] are included in the Appendix A.5 Figure
Conclusions
We propose very general framework for variational inference in Gaussian Graphical Models through conditional normalizing flows. Our model unifies the benefits of Bayesian and frequentist approaches while avoiding most of their specific problems. Compared to existing approaches, the most important advantage of our method is that can jointly train a continuum of sparse regression models for all regularization parameters and all l (pseudo-) norms. All of these models can be analyzed both in a Bayesian fashion (at temperature T 1) or, alternatively. in the frequentist limit (i.e. penalized likelihood) as T approaches 0. To the best of our knowledge, this is the first Gaussian Graphical Model that can continuously explore all sparsity-induein inducing priors from the norm family. Moreover, thanks our variational formalism, marginal likelihood for model selection without any additional computational costs.
integrate out the model parameters and the (approximate)
10

---

References
Stephen M Smith, Christian F Beckmann, Jesper Andersson, Edward J Auerbach, Janine Bijsterbosch, Gwenaëlle Douaud, Eugene Duff, David A Feinberg, Ludovica Griffanti, Michael P Harms, et al. Resting-state fmri in the human connectome project. Neuroimage, 00:144-168, 2013.
Robert Castelo and Alberto Roverato. A robust procedure for gaussian graphical model search from microarray data with larger than n. Journal of Machine Learning Research, 7(94))2611-2650, 2006.
P. Diggle. Analysis of Longitudinal Data. OUP Oxford, 2002. ISBN 9780198524847.
Trevor Hastie, Robert Tibshirani, and Martin Wainwright. Statistical Learning with Sparsity: The Lasso and Generalizations. Chapman &; Hall/CRC, 2015. ISBN 1498712169.
Nicolai Meinshausen and Peter Bühhmann. High-dimensionall graphs and variable selection with the Lasso. The Annals of Statistics, 34(3):1436 1462, 2006. doi: 10. 00.2114/099555560000001111
Jerome Friedman, Trevor Hastie, and Robert Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 93):432-441,, 12 2007. ISSN 1465-4644, doi: 00.t0///iiostaiistics// kxm045.
Ming Yuan and Yi Lin. Model selection and estimation in the gaussian graphical model. Biometrika, 94(1):19-35, 2007. ISSN 00063444, 14643510.
Onureena Banerjee, Laurent El Ghaoui, and Alexandre d' Aspremont. Model selection through sparse maximum likelihood estimation, 2007.
Rahul Mazumder and Trevor Hastie. The graphical lasso: New insights and alternatives. Electronic Journal of Statistics, (inone):2125 2149, 2012. doi: 10.1214/12-EJS444
Hao Wang. Bayesian Graphical Lasso Models and Efficient Posterior Computation. Bayesian Analysis, 7(4):867 886, 2012. doi: 10.1214/12-BA72).
Yunfan Li, Bruce A. Craig, and Anindya Bhadra. The graphical horseshoe estimator for inverse covariance matrices, 2019.
David M. Blei, Alp Kucukelbir, and Jon D. McAuliffe. Variational inference: A review for statisticians. Journal of the American Statistical Association, 112(518):859-877, apr 2017. doi: 1080/ 01621459.2017.1282
Larissa Alves, Ronaldo Dias, and Helio S. Migon. Variational full bayes lasso: Knots selection in regression splines, 2021.
S. Derin Babacan, Shinichi Nakajima, and Minh N. Do. Bayesian group-sparse modeling and variational inference. IEEE Transactions on Signal Processing, 62211)2996, 2014. doi: 10.1109/TSP2201::331777
Trevor Park and George Casella, The bayesian lasso. Journal of the American Statistical Association, 103(482):8816686 2008. doi: 10. 011119000661155000000000333
A. Mohammadi and F. C. Wit. Bayesian Structure Learning in Sparse Gaussian Graphical Models. Bayesian Analysis, 10(1):109 - 138, 2015, doi: 10.1214/14-BA89)
D. F. Andrews and C. L. Mallows. Scale mixtures of normal distributions. Journal of the Royal Staristical Society. Series B (Methodological), 36(1):99-102, 1974. ISSN 00359246.
Mike West. On scale mixtures of normal distributions. Biometrika, 74(3):6466443 09 1987. ISSN 0006-3444. doi: 10. (0.1099/biomet/74..6666
Jianqing Fan and Runze Li. Variable selection via nonconcave penalized likelihood and its oracle properties. Journal of the American Statistical Association, 96:1 1111148-11 1360, 02 2001.
Cun-Hui Zhang. Nearly unbiased variable selection under minimax concave penalty. The Annals of Statistics, 38(2):894 942, 2010. doi: 4/011

---

Jinchi Lv and Yingying Fan. A unified approach to model selection and sparse recovery using regularized least squares. The Annals of Statistics, 37(6A), dec 2009, doi: 10.1214/09-aos88..
Rahul Mazumder, Jerome H. Friedman, and Trevor Hastie. Sparsenet: Coordinate descent with nonconvex penalties. Journal of the American Statistical Association, 106(495): 125-1 138, 2011.
Lee Dicker, Baosheng Huang, and Xihong Lin. Variable selection and estimation with the seamless-I 0 penalty. Statistica Sinica, 23, 04 2013. doi: 10.5705/s.2011.072.
Yanxin Wang and Li Zhu. Variable selection and parameter estimation with the atan regularization method. Journal of Probability and Statistics, 2016:1 01 2016. doi: 10.1 15//0066665455417
Hemant Ishwaran and Sunil Rao. Spike and slab variable selection: Frequentist and bayesian strategies. The Annals of Statistics, 33(2), apr 2005. doi: 10. 4009000000000000000000001 147.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning. volume 37, pages 1530-1538. PMLR, 07-09 Jul 2015.
Rianne van den Berg, Leonard Hasenclever, Jakub M. Tomezak, and Max Welling. Sylvester normalizing flows for variational inference, 2019.
Chin-Wad Huang, David Krueger, Alexandre Lacoste, and Aaron Courville. Neural autoregressive flows, 2018
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshninarayanan. Normalizing flows for probabilistic modeling and inference. J. Mach. Learn. Res., 22(1), jan 2021. ISSN 1532-4435.
Ivan Kobyzev, Simon Prince, and Marcus Brubaker. Normalizing flows: An introduction and review of current methods. IEEE Transactions on Pattern Analysis and Machine Intelligence, PP:1-1, 05 2020. doi: 10.1 LO99/TTAA...00002999992.9.
Andrei Atanov, Alexandra Volokhova, Arsenii Ashukha, Ivan Sosnovik, and Dmitry Vetrov. Semi- conditional normalizing flows for semi-supervised learning, 2020.
Roger A. Horn and Charles R. Johnson. Matrix analysis. Cambridge University Press, repr. with corr edition, 1985.
Arjun K Gupta and Daya K Nagar. Matrix variate distributions, volume 104. CRC Press, 1999.
Fabricio Arend Torres, Marcello Massimo Negri, Marco Inversi, Jonathan Aellen, and Volker Roth. Lagrangian flow networks for conservation laws, 2023.
Conor Durkan, Artur Bekasov, Iain Murray, and George Papamakarios. nflows: normalizing flows in pytorch, 2020.
S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi. Optimization by simulated annealing. Science, 220 (4598):671-688) 1983. doi: 10. 211226scienee.2224559.....
Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions Pattern Analysis and Machine Intelligence, PAMI-6 (6):721-741, 1984. doi: 10.1119/TPAMI.199 1L109/TPAMI.19444777555009
David Abramson, Mohan Krishnamoorthy, and Henry Dang- Simulated annealing cooling schedules for the school timetabling problem. Asia-Pacific Journal of Operational Research, 1999.
Christophe Andrieu and Arnaud Doucet. Simulated annealing for maximum a posteriori parameter estimation of hidden markov models. Information Theory, IEEE Transactions on, 46:994 - 1004, 06 2000. doi: 10.1 109/18.841176
Changhe Yuan, Tsai-Ching Lu, and Marek Druzdzel. Annealed map. In Proceedings of the 20th Conference on Uncertainty in Artificial Intelligence, page 628-635. AUAI Press, 2004.
Huadong Liao and Jiawei He. Jacobian determinant of normalizing flows, 2021.
12

---

F. Pedregosa. G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel. M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and H. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 22:2825-2830 2011.
Michal Sheffer, Manny D. Bacolod, Or Zuk, Sarah F. Giardina, Hanna Pincas, Francis Barany, Philip B. Paty, William L. Gerald, Daniel A. Notterman, and Eytan Domany. Association of survival and disease progression with chromosomal instability: A genomic exploration of colorectal cancer. Proceedings of the National Academy of Sciences, 106(17::1111-1336 2009. doi: 000073//paas.9902222166....
Dinu Kaufmann, Sonali Parbhoo, Aleksander Wieczorek, Sebastian Keller, David Adametz, and Volker Roth. Bayesian markov blanket estimation, 2015.
Durk Kingma, Tim Salimans, Rafal Jozefowicz. Xi Chen, Ilya Sutskever, and Max Welling. Improved variational inference with inverse autoregressive flow. In Advances in Neural Information Processing Systems, volume 29, Curran Associates, Inc., 2016.
George Papamakarios, Theo Pavlakou, and lain Murray. Masked autoregressive flow for density estimation, 2018.
Fabricio Arend Torres. Sampling and annealing for dependency subnetwork estimation. Master's thesis, University of Basel, 2018.
13

---

A Appendix
A.1 Sum-of-Sigmoids layers
We implement the proposed model with a conditional flow architecture that relies on element wise transformations through monotonic functions, which term of--mo--s We further combine these element-wise monotonic transformations in an autoregressive fashion using a MADE-like approach [Kingma al. 2016 Huung et al.] 2018). The proposed layers are light and very flexible and work particularly well for our purposes already with only 4 layers. Note that, even though it is needed in our framework, the inverse can be computed numerically, which relatively cheap since the transformation is element-wise monotonic. Specifically, we implement flexible monotonic transformations by combining shifted and scaled sigmoid activations. Differently from Huang et al. (2118), we also add shifted (and flipped) softplus functions to the element-wine activations. This leads to a linear behaviour outside specific range [-s, . s). >> alleviating the effect of inputs that lic outside the seen training data. The strictly monotonic element-wist transformation is given by
Osos (2(10)
By a (10z= b)) + In (((00-0))) - In (1+e(-z()) el-
(12)
monotonic functions
upprax. lineur [z(() [ and ((r]
<
where is the sigmoid function and {2y, wy, by وزر and a are learnable parameters such that Wk, u/, a > 0 and Σド=1 V 1. The advantage of using an comenntwwist transformation is that the associated Jacobian is diagonal and therefore its log-determinant is straightforward to compute. The proposed Sum-of-Sigmoide layers provide very flexible transformation, but due to the element-wise nature, different dimensions are not mixed together. We overcome this by predicting the parameters {uz, Wy, 1ー b,ز and of the tum-of-Sigmoide with masked autoregressive hypernetworks, similar to other masked autoregressive flows. When implemented through masking [Papamakarios et al.) 2018), autoregressive flows allow for an elegant extension to conditional settings as well. Relevantly, this autoregressive structure admits a simple log-determinant Jacobian since the Jacobian is lower- triangular by construction.
A.2 Joint posterior for Ωn and P12
We show that given the partition in Eq. III and the prior in Eq. 7 the joint posterior factorizes as
p( 122 S222.1 ()) - p(2n S S2
S, q) p(02..11S,A,.0)
(13)
We further provide the analytic expression for p(S211, Q22 S, A, g). The latter is used target distribution in the training loss in Eq. (8) for the real data application in Section For this proof we followed the approach of Torres 2018
Let X E x be the design matrix containing independent observations. We are interested in estimating the connections between query variables with respect to the remaining ones, typically with & t. Given the partition of Ω and S in Eq. (II. we define the full posterior p(22, 5)),g) p(121 D12, S22 S),,q) as the product of the Wishart likelihood p( (2) and suitable prior p( A, g):
p(S[S2) & Wo+e(m 2- det(s)) exp (xp(IT(33 p(lA) g) & Wste(s
(14)
A1,+:) p( |A. q) p(2)2)))
Specifically, select Wishart prior J2 to ensure positive definiteness on the full matrix. In order to encourage sparsity we choose a generalized Normal distribution prior over the off-diagonal elements of Qn and over the full Ω12:
(SlulAg x II exp ( Al(2)))19) p(aia|A ②) x II exp ( (1((2)))11) (15)
[.
As first step explicitly write down the expression for the likelihood ((S p( P12, 22221)
(w+t(n,(t-1) x det( 122 (et( (1-(()))-- exp (- (- (r((ss))
(16)
14

---

and for the prior p(22/,,9) = p (f 8212 S22211 g) x Wetr(s+ W++ 1, (1)+1) ((121111) p(112)))
(17)
exp - II i
exp ( (1(5111) II exp ((((12))))
Overall, the posterior reads p( S212 222.11 SIA, 9) det( m/2 det (1-
1 Tr Tr + =lsI
x exp
(18)
II exp ( (1
II exp (S11())))))))))
i<j
Following Kaufmann al. [015].
re- wrrit the posterior through the change of variables
2 52 Ω22 (2 S22, P22) - (Ω 12 122211) involving the Schur component 1222.1 S222
Note that the transformation has unit Jacobian: (( S212 9222) - (S21 S22, 122211) = 1.
We can now explicitly re-write the posterior as p(2211 22 8222.1 S|,, 4) by using the following substitutions, which are straightforward show:
det(S2) det(S2n1) det 22 012) = det (ct
Tr(S2)) Tr T1 S125 + 2 122.1 S2 2 1 S22S22
(19)
Tr(At) Tr [>iı1 A0222. + 2 812]
By re-arranging the terms the posterior simplifies to:
p(l11 012, 1 det(S21))" )²/² dct(S) --(-)))-1
x det(S2221)
n1/2 exp
Tr 122. Tr (S22 A1t)
1 exp - Tr Ω (S
(L 2(2)) + 2 (2 + x1,)1)
II exp (- (1((1)))1)) Il exp ( (((()))))) iej
1.j
If we condition S, and specifically S22, we can clearly that the posterior factorizes as p(2n, D12, S22.1SSA X, 4) x p(2n, 1 SS X, p(22211 x) In particular, we are interested in estimating the joint posterior p( 12 S, X, q), which reads
p( 012 SS X, g) det(SLn1)m12
Tr (S + (L + S12) + 2 112 + xl,)1)
X exp
2
exp ( (1(()))))))
II
exp ( (1(22)))
(20)
A.3 Run-time comparison with Gibbs sampler
The proposed CMF provides a significant speed-up in terms of sampling with respect to standard Gibbs sampling algorithms for posterior inference in Gaussian Graphical Models. For realistic and fair comparison, we measure run-time on our real data experiment. For this purpose, we trained the Conditional Matrix Flow for g and compare its sampling speed against the Gibbs sampler introduced in Kaufmann et al. [2015). The setting presented here is the same one used to obtain the results shown the Experiment section. The Gibbs sampler takes 89 seconds to generate 500 samples, which translates to about 5.6 samples per second. This result was obtained on Intel(R) Xeon(R) CPU ES- 1660 v3 @ 3.00GHz. On the ohter hand, on the consumer consumer-grade GPU NVIDIA TITAN X (12GB VRAM) sampling is extremely efficient: each second we can generate 2000 independent samples from the approximate posterior. Relevantly, order to retrieve the full posterior path as a function of X, we would need independent Markov Chains for each À value, rendering the approach infeasible. In contrast, with the proposed framework we can draw independent samples for different A values at the same computational cost.
15

---

A.4 Artificial data
1.0
0.5
-0.5
10a
10a
A
(a) 4 1.00
(b) q 0.75
C
10a
10:
10a
10a
À
(c) = 0.50
(d) 0.25
Figure 6: MAP estimate function of X for different sub-l1 pseudo-norms. The dashed line the frequentist solution path for = 1.
GLasso CV
log pixxa
-200
-250
10a
10!
A
Figure 7: Evolution of the (approximate) marginal log- logeliieldhood of the CMF function of X for T 1.00 and g = 1. We select the optimal À as its maximum (ΛeMF 3.52, in blue) and compare with the frequentist estimate obtained through cross validation (AGLasso 3.36, in red).
16

---

A.5
Real data
a
CMF (q)0.6)
CMF (q-0.)
CMF (q-0.8)
(q90.9)
CMF (q-1.0)
(M,GS)
(N.GS)
(N,M)
(SEX.CTNNAL)
Figure 8: 95% posterior credibility intervals of the proposed CMF for the entries with highest absolute median posterior as a function of the pseudo-norm {1, 0.9, 0.8, 0.7, 0.6).
17


Reviews:
[{"review_id": "8DQ2f2lWKF", "reviewer": "Reviewer_CnQu", "summary": "This paper concerns the estimation of precision matrix under $l_p$ norm sparcity penal. The solution is a variational inference through normalizing flow, which is a function of shrinkage parameter $\\lambda$ and non-negative norm parameter $p$. It allows for straightforward computation of solution paths for the intervals of $\\lambda$ and $p$, and was empirically evaluated on two relatively small data sets.", "questions": "In synthetic data example, why did you choose to have more samples than dimensions ( $n>d$ )? Since in that case GGM can be obtained with matrix inverse, and no need for penalized objective.", "limitations": "Limitations were not discussed.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Framework for GGM estimation based on conditional normalizing flows, indeed appears novel. Supporting math seems solid. \n\nUsing simulated annealing algorithm to recover a path of solutions for varying $\\lambda$ and $p$ is useful, in particular for the case of $p$, as in case of $\\lambda$ it was fairly straightforward to perform it with other methods too. I am just wondering how costly and scalable it is under the new framework, an empirical/theoretical analysis would be appreciated.", "weaknesses": "Empirical evaluation appears limited. It does not contain comparison with other (e.g. frequentist) approaches to derive the solution paths. Both in terms of estimation accuracy and in terms of computational cost.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for acknowledging the novelty and soundness of our work. We hope to provide satisfying answers to the concerns raised.\n\n- \"Using simulated annealing algorithm to recover a path of solutions for varying $\\lambda$ and $p$ is useful, in particular for the case of $p$, as in case of $\\lambda$ it was fairly straightforward to perform it with other methods too. I am just wondering how costly and scalable it is under the new framework, an empirical/theoretical analysis would be appreciated.\"\n\nThe proposed approach is particularly suitable to study the evolution of the posterior both as a function of $\\lambda$ and $p$ (and of the two together). Bayesian Graphical Lasso allows to infer the posterior only for $p=1$ and requires to re-start the Markov Chain for every different $\\lambda$ value. Clearly, this becomes unfeasible very quickly. In contrast, posterior inference with the proposed CMF can be performed in few seconds on consumer-grade GPUs. The only computational cost is thus training, which however can be performed in around 10 minutes. Frequentist approaches are indeed faster and can compute solution paths in less than a minute. However, this is not surprising since they only recover a single point-wise path, while the presented CMF models a (infinite) family of posterior distributions, conditioned on both $\\lambda$ and $p$. It is however relevant that as a particular instance of our model we can recover the frequentist solution paths.\n    \n- \"Empirical evaluation appears limited. It does not contain comparison with other (e.g. frequentist) approaches to derive the solution paths. Both in terms of estimation accuracy and in terms of computational cost\"\n\nFor a more detailed comparison with frequentist (and Bayesian) approaches in terms of estimation accuracy, we provide further experiments in the global response. In terms of computational cost we refer to the previous answer.\n\n- \"In synthetic data example, why did you choose to have more samples than dimensions ($n>d$)? Since in that case GGM can be obtained with matrix inverse, and no need for penalized objective\"\n\nWe illustrated the presented model for $n>d$ just as a toy experiment. This is similar to the experiments in the original BGL paper by Wang et al (2012), where they consider $(d=30, n=50)$ and $(d=100, n=200)$. However, we agree that the relevant range is indeed $n<d$. In the global answer we included further experiments to specifically address this regime.\n    \n- \"Limitations were not discussed\"\n\nFor a discussion on the limitations of the proposed model we refer to the global response."}, {"review_id": "4gZSINDkER", "reviewer": "Reviewer_pTVu", "summary": "This paper targets the structure learning problem in Gaussian Graphical Models via (Normalising) Flow-based Variational approximation of the elements of weight metrics that correspond to the Gaussian Bayesian network. \nThey use sub-l1 pseudo norms to penalize dense precision metrics (which correspond to graphs with numerous links) without imposing an extra high penalty for large non-zero values (which typically occurs if $l_{\\geq1}$ is used).", "questions": "* In line 141, what do you mean by \"contradiction\"?", "limitations": "The authors should compare their method with the relevant structure learning lierature and reveal its points of strength as well as its limitations. \n\nThis work is theoretical/methodological and does not have any positive or negative social/ethical impact on its own.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "1. Up to my knowledge, this is the first time flows are applied to the space of positive definite matrices. \n2. The proposed approach is flexible meaning the class of applicable prior and likelihood functions is quite large.\n3. Using sub-l1 norm is suitable for structure learning. \n4. The proposed algorithm is mathematically sound (as far as I can follow) and is quite interesting. \n5. The paper is well-written, and the relevant work is sufficiently discussed.    \n6. Due to its flexibility, the proposed method has the potential of having a large impact.", "weaknesses": "Due to the factors mentioned in the previous section, I find this work impressive and beautiful. However, unfortunately, the carried out experiments are minimal. Most notably, the algorithm is compared to no alternative work (neither in the main paper nor in the supplementary material). With no quantitative comparisons, it is impossible to evaluate the performance of the proposed algorithm compared to the existing methods. \n\nNOTE: In the Rebuttal, some experiments are carried out (though the code is still not accessible).    \n\nMinor suggestion: \n1. Though it is clear in the context, I suggest that the authors do not use the same letter \"p\" (with the same font) for both probability density and norm parameter.  \n2. Fix minor typos e.g. the end sentence period in line 214.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the very encouraging comments and for acknowledging the potentially large impact of the proposed model. We hope to provide satisfying answers to the concerns raised.\n\n- \"Most notably, the algorithm is compared to no alternative work (neither in the main paper nor in the supplementary material). With no quantitative comparisons, it is impossible to evaluate the performance of the proposed algorithm compared to the existing methods. Due to weak evaluation, I do not think this work is ready for publication except if the authors can provide sufficient quantitative experimental results (along with code) in the rebuttal.\"\n\nFor a more detailed comparison with both frequentist and Bayesian approaches, we provide further experiments in the global response. We provide quantitative comparison in terms of accuracy in graphical structure learning and we compare the presented approach with the Bayesian Graphical Lasso and several frequentist method. In the Bayesian setting we recover BGL results and obtain higher F1 scores in the sub $l_1$ regime ($p<1$). The presented approach outperforms also classical frequentist algorithms across the sub-$l_1$ pseudo-norm and in the $l_0$ limit as well. We provide the code to reproduce our results as a comment to the AC.\n    \n- \"Though it is clear in the context, I suggest that the authors do not use the same letter 'p' (with the same font) for both probability density and norm parameter.\"\n\nOur notation was indeed confusing. We now changed it to $l_q$ norm and kept $p$ for probability density.\n    \n- \"In line 141, what do you mean by 'contradiction'?\"\n\nKaufmann et al. (2015) showed how to efficiently infer the posterior on the blocks $\\Omega_{11}$ and $\\Omega_{22}$ independently of the large block $\\Omega_{22}$. However, in order to derive the algorithm for the Gibbs sampler, they drop the sparsity prior on the $\\Omega_{11}$ block. In our opinion, this creates a contradiction in the model because the resulting model fails to enforce sparsity on the $\\Omega_{11}$ block.\n    \n- \"The authors should compare their method with the relevant structure learning literature and reveal its points of strength as well as its limitations.\"\n\nIn the global response, we provided further experiments specifically to evaluate structure learning and compare the proposed approach with both the Bayesian Graphical Lasso and several frequentist approaches. For a discussion on the limitations of the proposed model we refer to the global response."}, {"review_id": "7JKoCjh5YL", "reviewer": "Reviewer_3sWQ", "summary": "This paper proposes a method that can be used to infer conditional independencies in a Gaussian model. These conditional independencies are related to zeros in the precision matrix. Typically, sparse enforcing norms are used to estimate the precision matrix while enforcing zeros in the elements outside of the diagonal. In this paper a Bayesian approach is considered. For this a pseudo-distribution for the data is considered by taking the exponential to the p-norm. The method is trained via variational inference combined with normalizing flows to increase the accuracy of the posterior approximation. The variational distribution is tuned via simulated annealing and a temperature parameter allows to interpolate between the Bayesian and the Map solution.", "questions": "None", "limitations": "The authors have not commented on the limitations of their approach.", "rating": 5, "confidence": 3, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "- Well written paper.\n\n        - Illustrative toy experiments.", "weaknesses": "- The proposed method is a combination of already known techniques.\n\n        - The experimental section is weak as only a single real problem is considered.\n\n        - Although the proposed method is a generalization of several known techniques, I have found in the experimental section a lack of comparisons with other related methods.\n\n        My main point of criticism is the weak experimental section which only considers a single real problem and no comparisons with other related methods are carried out in real problems.\n\n        Another point of criticism is that, for some particular values of the p parameter one does not actually observe sparsity in the Bayesian solution. For example, when sampling from the Laplace distribution one never observes zeros in practice. Spike and slab priors (a mix between a Gaussian and a point of mass center at zero) are the ones that actually lead to zeros.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the useful feedback on the experiment section and for acknowledging the clarity in the presentation of the proposed approach. We hope to provide satisfying answers to the concerns raised.\n\n- \"The proposed method is a combination of already known techniques.\"\n\nWe believe that one of the benefits of the proposed approach is indeed to provide a general framework that unifies the two known approaches, namely the Bayesian and frequentist ones. We do so by combining their advantages: we have access to the marginal likelihood for model selection, we can compute credible intervals and we can recover frequentist solution paths. However, we also improve upon existing methods. As a key advantage to current approaches, we (i) are able to model the posterior in the sub-$l_1$ norm regime and we (ii) have access to a continuum of Bayesian models as a function of both $\\lambda$ and $p$. In the global response we (iii) show quantitatively that the proposed approach outperforms both the Bayesian Graphical Lasso and several frequentist approaches in terms of graphical structure recovery.\n\n- \"The experimental section is weak as only a single real problem is considered\"\n\nIn the paper we considered one real problem to showcase that the proposed method can be readily applied to high dimensional real world settings. As part of the rebuttal, we decided to focus on further synthetic experiments in order to have access to the ground truth for a quantitative comparison with competing methods.\n    \n- \"Although the proposed method is a generalization of several known techniques, I have found in the experimental section a lack of comparisons with other related methods\"\n\nFor a more detailed comparison with Bayesian and frequentist approaches, we provide further experiments in the global response.\n    \n- \"for some particular values of the p parameter one does not actually observe sparsity in the Bayesian solution. For example, when sampling from the Laplace distribution one never observes zeros in practice. Spike and slab priors (a mix between a Gaussian and a point of mass center at zero) are the ones that actually lead to zeros\"\n\nContrary to frequentist approaches, Bayesian approaches provide credible intervals on the solution. As correctly pointed out, this implies that individual samples may not contain exact zeros. Indeed, this is precisely what motivates our interest in studying solutions in the sub-$l_1$ pseudo-norms, which lead to exact zeros in the $p\\rightarrow 0$ limit. \n    \n- \"The authors have not commented on the limitations of their approach\"\n\nFor a discussion on the limitations of the proposed model we refer to the global response."}, {"review_id": "ImVCECqCY5", "reviewer": "Reviewer_vqBu", "summary": "This work proposed a framework for performing inference on Gaussian Graphical Models by approximating the posterior with a normalizing flow over PSD matrices. In this way, the authors can investigate $l_p$-norm regularized GGMs for any value of $p$ in an efficient way.", "questions": "- in the frequentist case, how does the proposed framework compares against more classical techniques to obtain the solution paths, e.g., iterative reweighted l1-norm? ", "limitations": "The authors adequately addressed the limitations of their proposed framework. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The idea of using normalizing flows for GGM inference definitely brings in advantages of both Bayesian and frequentist worlds; to me, that's an innovative idea.", "weaknesses": "The main weakness that I identified is the lack of comparison between the proposed framework and the well-studied graphical lasso with concave approximations of the $l_0$-norm. More precisely, the authors show that their framework obtains frequentist solution paths through simulated annealing, therefore, it'd be of great interest to see a comparison between these solution paths and those obtained by iterative algorithms such as iterative reweighted l1-norm for graphical lasso.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for acknowledging the novelty of the proposed framework and for suggesting relevant experiments, which we have included in the global response. We hope to provide satisfying answers to the concerns raised.\n\n- \"The main weakness that I identified is the lack of comparison between the proposed framework and the well-studied graphical lasso with concave approximations of the $l_0$ norm.\"\n\nIn the global response we provide further experiments to compare the proposed method against several frequentist approaches with penalties that approximate the $l_0$ norm. We include Atan by Wang and Zhu (2016), Seamless L0 by Dicker et al. (2013), Log by Mazumder et al. (2011) and SICA by Lv and Fan (2009). Overall, results show that the proposed method outperforms all the considered frequentist approaches in terms of F1 score for structural recovery. \n    \n- \"More precisely, the authors show that their framework obtains frequentist solution paths through simulated annealing, therefore, it'd be of great interest to see a comparison between these solution paths and those obtained by iterative algorithms such as iterative reweighted l1 norm for graphical lasso. Questions: in the frequentist case, how does the proposed framework compares against more classical techniques to obtain the solution paths, e.g., iterative reweighted l1-norm?\"\n\nIn the global answer we provide an experimental comparison specifically for an iterative re-weighted lasso approach, namely the Adaptive Lasso by Zou (2006). Results show that the proposed method consistently outperforms Adaptive lasso over a range of $l_p$ norms and across all number of samples. Furthermore, we compare our model against two other frequentist approaches with penalties that allow to interpolate between $l_1$ and $l_0$ norms: SCAD by Fan and Li (2001) and MC+ by Zhang (2010). Also in this case we outperform competing methods. For lack of space we have not included these results in the 1-page PDF. However, they can be reproduced by running the code linked in the rebuttal."}]

Author rebuttal:
## Global response
We would like to thank all reviewers for the very encouraging comments and for acknowledging the novelty and the potentially large impact of the proposed approach. As all reviewers expressed one main concern, we try to address it in the global response. The main point of criticism concerns the experimental section and in particular the lack of comparisons with classical frequentist approaches for Gaussian Graphical Models (GGMs). As pointed out by reviewers vqBu and 3sWQ, the relevance and novelty of the proposed method is to unify the benefits of the Bayesian and frequentist approaches. However, as a key advantage on current Bayesian methods, the proposed method enables us to additionally explore all sparsity-inducing priors from the sub $l_1$ norm family. Even though the proposed method is inherently Bayesian, we agree that it is indeed interesting to compare its behaviour with the frequentist approaches. As a sanity check, in the paper we showcased that the proposed model recovers the frequentist solution path in the annealed limit. In the hope of resolving the concerns raised, we provide three additional experiments to compare the proposed approach with Bayesian and frequentist methods. As suggested, we study the performance in graphical structure learning (reviewer pTVu) in the $n<d$ regime (reviewer CnQu) and we focus on penalties that approximate the $l_0$ norm and sub-$l_1$ pseudo-norms (reviewer vqBu).

### Further experiments
We show that in the in the sub-$l_1$ regime the proposed Conditional Matrix Flow approach (CMF) outperforms both Bayesian Graphical Lasso (BGL) and the classical frequentist approaches.  In particular, we measure the F1 score for edge recovery as a function of the number of samples.  We first compare the CMF with the BGL, which is the only alternative Bayesian model. Then, we compare the CMF in the $p\rightarrow0$ limit against the frequentist approaches with approximate $l_0$ norm penalties: Atan ("atan") by Wang and Zhu (2016), Seamless L0 ("selo") by Dicker et al. (2013), Log ("log") by Mazumder et al. (2011), SICA ("sica") by Lv and Fan (2009).  Finally, we compare the CMF solution as a function of $p$ against frequentist approaches with penalties that allow to interpolate between $l_1$ norm and $l_0$ norm: SCAD ("scad") by Fan and Li (2001), MC+ ("mcp") by Zhang (2010) and Adaptive lasso ("adapt") by Zou (2006).

We use 10 ground truth precision matrices of dimension $d=30$ and generate $n$ Gaussian samples accordingly. For illustrative purposes we report results for the relevant regime $n<d$ and around $n=d$, namely for samples $n=\{15, 25, 35, 45\}$. More exhaustive results can be obtained by running the provided code. For Bayesian approaches (the proposed CMF and BGL) we draw 1000 samples from the (approximate) posteriors and consider $90\%$ credibility intervals. All results are averaged over the 10 precision matrices. The proposed CMF is trained for 5000 epochs to a final temperature $T=1$, which corresponds to the Bayesian model in Eq. (9). For the BGL we run the Gibbs sampler for 4000 iterations with a burn-in of 1000. We thin the samples and keep every fourth. The frequentist approaches were run each with the specific hyper-parameters suggested in their original papers. 

1.  As the proposed method is inherently Bayesian we first compare it against the BGL. We show that for $p=1$ we recover the BGL results. Furthermore, with the proposed method we can additionally infer the posterior for $p<1$. The results in Figure 1 show that in the low sample regime ($n<d$) sub $l_1$-pseudo norms are beneficial and result in higher F1 score. The effect is stronger as $p \rightarrow 0$, while in the $n>d$ regime sub-$l_1$ pseudo-norms do not provide a significant advantage.
2. We compare the proposed CMF also with classical frequentist approaches with surrogate penalties that approximate the $l_0$ norm. Results in Figure 2 show that the proposed CMF with $p=0.25$ outperforms all competing methods across all number of samples, especially for very low sample sizes ($n=15$). Note that in the $n<d$ regime frequentist algorithms require an ad-hoc initialization of the precision matrix, which we provided through the Ledoit-Wolf shrinkage estimator.
3. Lastly, we compare the CMF against frequentist approaches with penalties that allow to smoothly interpolate between the $l_1$ and $l_0$ norms as a function of the parameter $\gamma$. We evaluate the proposed CMF for $p=\{0.25,0.5,0.75,1\}$. Results show that we outperform all baseline methods. In Figure 3 we illustrate the results for the best baseline approach, namely Adaptive Lasso. Both the CMF and the Adaptive Lasso show increased performance with more regularization (for $p \rightarrow 0$ and $\gamma \rightarrow 0$) but the CMF still provides a higher F1 score across all number of samples. Results for SCAD and MC+ can be obtained by running the provided code.

The link to the code for the experiments can be found in a separate comment to the AC.

### Limitations
The proposed approach provides a general framework for posterior inference in GGMs and generalizes to a large family of likelihood and priors, which requires to efficiently evaluate their analytic expression. In some cases, e.g. for posterior predictive distributions, this would require an additional integration step, which could become computationally expensive. The proposed approach is also limited by the expressive power of the bijective layers of the flow. Even though state-of-the-art layers are extremely powerful in modelling high-dimensional distributions, each layer might still be limited in terms of the number of modes that can be modeled (see Liao et al (2021)). Lastly, the proposed model assumes that we can model a family of posterior distributions as a function of the conditioning parameters $\lambda$ and $p$, which ultimately depends on the flexibility of the hyper-network used.

Meta-review: This paper introduces a novel approach for inferring conditional independencies within a Gaussian model. These independencies correspond to the presence of zeros in the precision matrix. While conventional techniques involve employing norms that encourage sparsity to estimate the precision matrix while maintaining zeros off the diagonal, this study adopts a Bayesian perspective.

In this context, a pseudo-distribution is constructed for the data by exponentiating the p-norm. The method leverages variational inference along with normalizing flows to enhance the accuracy of the posterior approximation. The variational distribution is fine-tuned using simulated annealing, and a temperature parameter facilitates interpolation between the Bayesian and Maximum A Posteriori (MAP) solutions.

A particularly innovative aspect is the framework for estimating Gaussian Graphical Model (GGM) using conditional normalizing flows. The utilization of a simulated annealing algorithm allows for the exploration of a solution path distribution, contrasting the single path typical of frequentist approaches. Although an initial concern about limited experimental validation was raised, the authors addressed this issue during the rebuttal phase.

Given these insights, I recommend the acceptance of this paper.
Decision: Accept (poster)