Efficient Regularization Parameter Selection for Latent Variable Graphical Models via Bi-Level Optimization

1 Introduction

Multivariate Gaussians N (µ, Σ) with mean µ and covariance matrix Σ are still among the most popular probabilistic models. The mean vector and the covariance matrix need to be estimated from data. The covariance matrix of an n-variate Gaussian has (n+1 2 ) free parameters which is fairly large already for moderately high dimension n. This large number of parameters may lead to overfitting, or even worse, to nonregular estimates of the covariance matrix if there are fewer than n data points.

∗Contact Author

Factor analysis, developed by [Spearman, 1904], can be used to address the problems of overfitting and non-regular covariance matrices of standard multivariate Gaussians. In a Gaussian factor model the dimensions are subdivided into observed and a small number of latent dimensions. The joint density for all dimensions is assumed to be a multivariate Gaussian. Hence, also the marginals for the observed and latent dimensions, respectively, are multivariate Gaussians. The marginal for the latent dimensions is assumed to be N(0, 1) and the marginal for the observed dimensions is assumed to be N(µ, Ψ + ΓΓT), where Ψ is a diagonal matrix with non-zero diagonal entries and Γ is a rectangular matrix that maps from the latent to the observed dimensions. The covariance matrix Ψ + ΓΓT has only n + kn free parameters, where n is now the number of observed and k the number of latent dimensions. The non-zero diagonal of Ψ ensures the regularity of the covariance matrix.

A different approach for dealing with the overfitting problem is the assumption of a graphical model structure. In the case of multivariate Gaussians a graphical model structure is defined by the zeros in the inverse covariance matrix (precision matrix). The structure, i.e., the precision matrix Λ, of a Gaussian graphical model can be learned consistently through the graphical Lasso [Yuan and Lin, 2007], which is given by the following optimization problem

where ℓ(Λ) = trace (Λ̂Σ) - log det(Λ) is the negative loglikelihood with empirical covariance matrix ̂Σ and regularization parameter λ > 0 for the sparsity promoting ℓ1regularization term. The learned graphical model structure of course depends on λ that has to be chosen carefully.

[Chandrasekaran et al., 2012] observed that it is possible to combine the ideas of Gaussian factor models and Gaussian graphical models. As in factor analysis, they assume that the observed dimensions along with a small number of latent dimensions are jointly multivariate Gaussian. The effect of the latent dimensions on the precision matrix of the observed dimensions can be computed from the Schur complement Λo - ΛolΛ-1 l Λlo, where Λ = ( Λo Λol Λlo Λl ) is the precision matrix of the joint distribution. Hence, the precision matrix for the observed dimensions is the difference of the two matrices S = Λo and L = ΛolΛ-1 l Λlo, where L has low rank, if there are much fewer latent than observed dimensions. In the spirit of Gaussian graphical models, the matrix S is assumed to be sparse. Thus, the precision matrix of the observed dimensions is the difference of a sparse and a low-rank matrix. [Chandrasekaran et al., 2012] show that these matrices can be learned consistently through the semidefinite program (latent variable graphical Lasso)

Two examples for such a sparse + low-rank decomposition of the precision matrix can be seen in Figure 1.

In a comment that was published together with the original work on latent variable graphical models it was pointed out by [Giraud and Tsybakov, 2012] that "Proposing a reasonable data-driven selector for [. . . the regularization parameters λ and ρ . . . ] would be very helpful for the practice". Choosing a good value for the regularization parameter is already an issue for the graphical Lasso since the statistical performance of the resulting model is highly dependent on this choice. The problem is aggravated for latent variable graphical models, where two instead of only one parameters need to be explored, compare Figure 1.

The problem of choosing a good value for a single regularization parameter has received a lot of attention, see for example [Giraud et al., 2012] who study the special case of the graphical Lasso. The case of two or more regularization parameters has received less attention and is addressed rigorously only in [Blechschmidt et al., 2015] and recently in [Giesen et al., 2019].

Especially in the higher-dimensional case, i.e., for two or more regularization parameters, we found it very convenient to model the parameter-selection task as a bi-level optimization problem (Stackelberg game [von Stackelberg, 2010]). The two levels in such an approach are mostly referred to as upper and lower level, but are also known as leader and follower in game theoretic settings. Any solution or decision taken by the leaders on the upper level to optimize their goals is affected by the response of the followers on the lower level, who will seek to optimize their own outcomes. Here, on the upper level a measure of generalization error is minimized, for instance a cross validation error in a supervised setting or a likelihood value in an unsupervised setting. On the lower level the set of feasible solutions for all combinations of regularization parameters, the so called solution gamut, is described.

Typically, it is not feasible to provide a description of the full solution gamut and thus schemes for approximating it have been devised. [Blechschmidt et al., 2015] provided the definition of an approximate solution gamut, devised an algorithm for computing it, and proved that the algorithmic complexity of their algorithm is asymptotically optimal. Recently, [Giesen et al., 2019] revisited the problem and showed that Benson's vector optimization algorithm [Benson, 1998] can be adapted for approximating the solution gamut. So far, Benson's algorithm has been used only for regularized problems on standard vector domains but not for semidefinite problems.

Benson's algorithm avoids many of the practical issues of the original solution gamut approximation algorithm. These issues include choosing a sufficiently fine grid, evaluating the objective function on this grid, and solving dual problems next to the primal ones. The derivation of the dual problem is necessary for providing a stopping criterion in form of a small duality gap, in terms of a prescribed accuracy ε > 0, at every point of the solution gamut. Benson's algorithm uses a slightly different stopping criterion, but nevertheless comes with strong approximation guarantees. Here, we adapt and extend Benson's algorithm to the semidefinite programming setting and apply it to the regularization parameter selection problem for latent variable graphical models. It should be noted, however, that Benson's algorithm, in contrast to the traditional approach, can be used outof-the-box for regularization parameter optimization for any other machine learning problem that is cast as a regularized semidefinite program like for instance [Cand`es et al., 2011; Cand`es and Recht, 2009; Chen et al., 2017; Chen et al., 2011; Song et al., 2007].

In Section 2, we provide the details on our bi-level optimization approach and our extension of Benson's algorithm. We also discuss how to implement the traditional solution gamut approximation algorithm for the latent variable graphical model. That includes deriving the dual of Problem (LVGL)-a maximum entropy problem that is interesting in itself. In Section 3, we experimentally compare Benson's algorithm to the original solution gamut approximation algorithm on a range of data sets from various domains. It turns out that Benson's algorithm that, in contrast to the original algorithm, can be used out-of-the-box is also more efficient in practice.

2 Regularization Parameter Selection by Bi-Level Optimization

By rescaling the objective of Problem (LVGL), one can see that the latent variable graphical Lasso is an optimization problem of type

with convex feasible set C ⊆ Rn×n . The functions ℓ and ri, i = 1, . . . , q, are assumed to be convex. The objective fw is weighted by the regularization parameters w = (w0, w1, . . . , wq) ∈ S ⊆ Rq+1 , where

defines the q-dimensional standard simplex. By Xw , we denote a global optimal solution to Problem (Pw) for some given weight combination w ∈ S.

Remark 2.1. Let us emphasize that next to constraints of type g(X) ≤ 0 with some convex function g, in particular, also positive semidefinite matrix constraints of type X ⪰ 0 are allowed for defining the feasible set C. ⋄

Choosing good weight parameters w is typically done in a two-step procedure:

1. Solve (Pw) for "every" w ≥ 0 on training data.
2. Choose "the best" w, such that Xw minimizes some type of generalization error GE(w) on validation data.

From the viewpoint of optimization theory, Steps (1) and (2) can be combined and also be understood as a bi-level optimization problem

In this hierarchical optimization problem, the upper level aims for optimizing the objective GE with variable w. On the lower level, for each choice of w one has to find a solution Xw of Problem (Pw). It is known (cf. [Dempe, 2002]) that, in general, problems of type (BiP) are non-convexeven if (Pw) is a convex problem and GE is a convex objective. This can be observed practically in Figure 3, where several local minima for GE exist. Indeed, bi-level programming is known to be NP-hard [Hansen et al., 1992].

For solving Problem (BiP), one has to deal with its feasible set, which is given by the set of solutions of the lower level problem (Pw), depending on the upper level variable w. Since, in most cases, it is not possible to give a closed formula for Xw , we use an approximation based on the definition of the solution gamut for Problem (Pw), cf. [Blechschmidt et al., 2015].

Definition 2.2. Let ε > 0 be given. We call some function ̂X: S → Rn×n an ε-approximative solution gamut of Problem (Pw), if for all w ∈ S

In the following, we will briefly present two methods for the computation of such an ε-approximative solution gamut.

The first one is the algorithm presented in [Blechschmidt et al., 2015], which we just call solution gamut method. The second one is a variant of Benson's dual algorithm, an established method from the area of vector optimization.

Since both algorithms yield a finite representation of the εapproximative solution gamut ̂X, minimization of the upper level objective GE just boils down to function evaluations.

2.1 Solution Gamut Method

The algorithm described in [Blechschmidt et al., 2015, Sec. 4] computes a piecewise linear ε-approximative solution gamut for Problem (Pw). For that, the regularization parameter domain S is covered by a sufficiently fine grid. The spacing of this grid depends on the objective and the problem data, and thus is difficult to choose in practice. For the initialization, the method computes the optimal solution Xw to Problem (Pw) for w = 0q+1, as well as the optimal solution to the dual problem of (Pw) which has to be derived first. Then, for every grid point, i.e., every weight w ∈ S, the duality gap of this primal-dual solution pair is computed. The duality gap serves as upper bound to the optimality gap used in Def. 2.2. In the iterative process, the algorithm now chooses the grid point ˆw ∈ S with the maximal stored duality gap. Then, the primal and dual solutions for Problem (Pw) with w = ˆw are computed, and the duality gaps are updated at all grid points, where the stored duality gap is larger than ε. The algorithm stops as soon as the duality gap is smaller than ε for every grid point w ∈ S.

2.2 Adaptive Benson Algorithm

In a recent work [Giesen et al., 2019], it turned out that the so-called Benson algorithm is well-suited for the computation of an ε-approximative solution gamut of Problem (Pw).

Vector Optimization

To apply Benson's algorithm, Problem (Pw) has to be studied in the context of vector optimization. Thereby, Problem (BiP) becomes a semi-vectorial bi-level optimization problem, where the scalar objective GE has to be minimized over the Pareto set of the following convex vector optimization problem

The objective function F : Rn → Rq+1 is vector-valued and minimized w.r.t. the component-wise partial ordering

where Rq+1 + = {y ∈ Rq+1 | yi ≥ 0, i = 1, . . . , q + 1}. Problem (Pw) is the weighted sum scalarization of (VP).

In [Giesen et al., 2019, Sec. 2] it has been observed that the full solution gamut (Def. 2.2 with ε = 0) coincides with the set of weak minimizers to Problem (VP).

Definition 2.3. A point X∗ ∈ C is called a weak minimizer of Problem (VP) if ({F(X∗ )} - intRq+1 + ) ∩ F(C) = ∅, where F(C) = {F(X) ∈ Rq+1 | X ∈ C} is the image of the feasible set. ⋄

Since, in general, the set of weak minimizers is infinite, we have to consider finite ε-solutions. Therefore, it is important to understand the geometry of Problem (VP)-especially its upper image defined by

The practical computation of an ε-approximation of P is motivated by the following definition.

Definition 2.4. Let c ∈ intRq+1 + be arbitrary but fixed and assume that Problem (VP) is bounded, i.e., P ⊆ {y} +Rq+1 + holds for some y ∈ Rq+1 . Then a non-empty finite set C ∗ ⊆ C is called an ε-infimizer if

An ε-infimizer C ∗ of (VP) is called a weak ε-solution to (VP) if it only consists of weak minimizers. ⋄

By setting Pε = conv F (C ∗ )+Rq+1 + , such an ε-solution C ∗ provides an inner polyhedral approximation of the upper image P by finitely many minimizers, see Figure 2. Following the arguments in [Giesen et al., 2019], a weak ε-solution to Problem (VP) coincides with an ε-approximative solution gamut for Problem (Pw).

Adaptive Benson Algorithm

For the computation of a weak ε-solution for Problem (VP), we develop an adaptive variant of Benson's dual algorithm. The class of dual Benson algorithms aims for approximating the upper image P by iteratively generating a growing sequence of inner approximating polyhedra until some prescribed accuracy is reached. Since polyhedra play a crucial role in this method, remember that they are either given in H-representation (intersection of half spaces) or in Vrepresentation (set of vertices and directions). The conversion between both representations is done by vertex and facet enumeration, resp. [Bremner et al., 1998].

While the standard Benson algorithm requires the choice of an approximation accuracy ε beforehand, we devise a variant which works adaptively, since the choice of ε cannot be done generically in practice. Our Algorithm 1 starts with a coarse accuracy ε0 > 0 and calculates an ε0-approximation I ε0 of the upper image P. Then the accuracy is successively refined by setting εi+1 = εi 2 until the resulting approximation of P satisfies some stopping criterion. This leads to a growing sequence of inner approximation polyhedra given by

To calculate an (intermediate) inner approximation I εi , a hyperplane wTx = b from the H-representation of the current inner approximation is moved outwards until it becomes a supporting hyperplane of the upper image P. A contact point of this supporting hyperplane and P is given by F(Xw ), where Xw is the solution to the scalarized Problem (Pw). The distance the hyperplane is moved outwards in a given direction c is tracked and calculated by

If dc(Xw ) ≥ εi, the vertex F(Xw ) is added to the working V-representation followed by an update of the working H-representation, done by facet enumeration. If otherwise dc(Xw ) < εi, we continue by checking the next hyperplane. The procedure is repeated until the εi-approximation of P is complete. For the next iteration, we refine the accuracy εi+1 and initialize the working V-representation with all vertices F(Xw ) where dc(Xw ) ≥ εi+1, since the inner approximation in the neighborhood of those vertices may need further improvement.

The initialization of Algorithm 1 itself can be done with a V-representation that consists of a single point F(Xw0 ) for some initial weight w0. For instance, one can choose w0 = ( 1 q+1 , . . . , 1 q+1 ) ∈ Rq+1 and set Ipoi = {F(Xw0 )}.

Stopping criteria can be (i) attaining a target ε-approximation of P, (ii) completing a maximum number of iterations, or (iii) improving the generalization error by a certain degree.

Algorithm 1 Adaptive Dual Benson Algorithm

Input: Problem data (F, C), initialization Ipoi, chosen weights T = {w0}, direction c, initial accuracy ε = ε0 Output: V-representation Ipoi

1: function ADAPTIVEBENSONALGORITHM

2: repeat 3: I ε poi ←{F(Xw ) ∈ Ipoi | dc(Xw ) ≥ ε} 4: compute H-representation I ε of I ε poi 5: repeat 6: choose w ∈ I ε \ T 7: Xw ← arg min(Pw) 8: if dc(Xw ) > ε then 9: I ε poi ← I ε poi ∪ {F(Xw )} 10: update H-representation I ε of I ε poi 11: end if 12: T ← T ∪ {w} 13: until I \ T = ∅ 14: Ipoi ← I ε poi ∪ Ipoi 15: ε ← ε/2 16: until Stopping Criterion 17: end function

Complexity. Surprisingly, it seems that Benson's algorithm has never been used for semidefinite programs before. In [L¨ohne et al., 2014], the authors only allow convex constraints for Problem (VP) defined by polyhedral cones. By adding semidefinite constraints, their proof of convergence still works, since by [L¨ohne et al., 2014, Rem. 3 (Sec. 4.3)] one only has to assume that (i) int C ̸= ∅ (Slater's condition) and (ii) Problem (Pw) has a solution for all w ∈ S. Thereby, Algorithm 1 returns a weak ε-solution on termination. It is known from [Blechschmidt et al., 2015] that the lower bound for the number of optimization problems that need to be solved for an ε-approximative solution gamut is Ω(ε-q/2 ). An upper bound for the complexity of Bensontype algorithms can be found in [Kamenev, 1994, Theorems 3 and 4]. This bound is given by O(ε-q ) and O(ε-q/2 ), resp., where for the sharp second bound a twice continuously differentiable boundary of the upper image P is assumed.

3 Experiments

We perform regularization parameter selection for Problem (LVGL) from the introduction using the following data sets [Tsanas et al., 2014; Higuera et al., 2015; Zhou et al., 2014; Dua and Karra Taniskidou, 2017]

• GENE1 with n = 100 features and m = 255 samples,
• TCGA with n = 500 features and m = 801 samples,
• MICE with n = 81 features and m = 552 samples,
• ROSETTA with n = 100 features and m = 301 samples,
• SONAR with n = 60 features and m = 208 samples,
• OR70 with n = 70 features and m = 1059 samples,
• LSVT with n = 310 features and m = 126 samples,
• S&P500 with n = 471 features and m = 60 samples.

These data sets cover a wide range of applications. GENE1, TCGA, and MICE are biological data sets, ROSETTA and SONAR are geological data sets, OR70 was recorded for investigating the geographical origins of music, LSVT is about voice rehabilitation in psychology, and S&P500 includes monthly stock return data from major US companies over the course of 5 years. From S&P500 we removed 29 companies because their data was incomplete. From the original data sets GENE1, ROSETTA, and TCGA we selected the n features with the highest variance, similarly as [Chandrasekaran et al., 2012] who also used only subsets of GENE1 and ROSETTA.

The Dual of Problem (LVGL). The solution gamut method requires the dual of Problem (LVGL) that we derive here. [Dud´ık et al., 2004] have developed a fairly general duality theory of discrete maximum likelihood and maximum entropy problems. In this theory, regularization terms on the maximum likelihood side translate themselves into relaxations of moment-matching constraints on the maximum entropy side. The theory can be extended to the continuous case that we need here, if Shannon entropy is replaced by differential entropy. In our case, it turns out that the maximum entropy solution is a multivariate Gaussian with zero mean and covariance matrix Σ. It follows that the objective function of the dual of Problem (LVGL) is logdet(Σ) since the entropy of a multivariate Gaussian is proportional to the log determinant of its covariance matrix. The regularization terms of Problem (LVGL) become relaxations of the momentmatching constraint Σ = ̂Σ in the dual maximum entropy problem, which reads as

Experimental Setup

Any algorithm for approximating the solution gamut needs to guarantee a good solution in terms of the generalization error GE while being computationally efficient. In our experiments we study both aspects, the quality of the solutions and the computational efficiency of both approximation algorithms, the solution gamut method and our adaptive Benson algorithm. For solving the optimization problems, we adapted the ADMM-based algorithm discussed in [Ma et al., 2013] and used CDD [Bremner et al., 1998] for facet enumeration. All experiments were run on a Linux machine with an Intel Core i5-2500K (4 × 3.30 GHz) CPU and 16 GB RAM.

Solution quality. We prepared the data sets by splitting them into training and validation data in a 2:1 ratio. We also centralized and standardized the data using empirical means and standard deviations of the training data. The generalization error, here the negative log-likelihood function value, is computed on the validation data.

For the experiments, we reformulated Problem (LVGL) such that it conforms to the theory in Section 2, where weights are chosen from a standard simplex.

As a baseline, we used grid search on a fine grid with more than 5000 points. At the grid points we solved Problem (LVGL) and computed the corresponding generalization error, see Figure 3. The interval between the worst and best solution, in terms of generalization error, serves as our estimate of the range of possible generalization errors. We consider a solution as good, if it is close to the best solution. We measure closeness as the normalized distance to the optimum, where the normalization factor is the length of the interval of possible generalization errors. Then, we let both the solution gamut method and the adaptive Benson algorithm run until their solutions were within 1% distance to the optimum from the grid search. On all data sets, we started Benson's algorithm with ε0 = 210 and fixed direction c = (1, 1, 1)T .

Computational efficiency. It is known that the solution gamut method asymptotically matches the theoretical lower bound of Ω(1/ε) optimization problems that need to be solved for an ε-approximation. The pessimistic known theoretical upper bound for Benson's algorithm is only in O(1/ε2 ). In practice, however, one also has to consider the overhead for function evaluations on a fine grid for the solution gamut method and the overhead incurred by the adaptive Benson algorithm for facet enumerations. Hence, in our experiments we counted not only the number of solved optimization problems, but also measured the elapsed CPU time.

Results

Solution quality. On the LSVT data set, we stopped Benson's algorithm at the fairly coarse accuracy of ε = 26 . This is already sufficient for finding good solutions with different algebraic properties (sparsity and rank), see Figure 3. Although the search region for the upper level problem, i.e., minimizing GE, is non-convex, Benson's algorithm finds solutions close to all local minima. Hence, in practice the best solutions returned by Benson's algorithm should be taken into account, because they may provide alternatives in terms of their algebraic properties, see Figure 1.

Computational efficiency. Experimental results for the comparison of the adaptive Benson algorithm and the solution gamut method are shown in Table 1. The solution gamut method requires solving between two and seven times more optimization problems than Benson's algorithm. The difference becomes even more pronounced when looking at the actual CPU times-it is a factor between 6 and 34 (median 13), i.e., about one order of magnitude.

Furthermore, it turns out that the adaptive Benson algorithm experimentally matches the theoretical lower bound of Ω(1/ε), see Figure 4. This suggests that the optimistic upper bound for Benson's algorithm is realistic in our case.

4 Conclusion

Latent variable graphical models are a sparse + low-rank variant of multivariate Gaussians that address some of the shortcomings of the latter. Like Gaussian factor models or the graphical Lasso, this model is less prone to overfitting and does not suffer from the problem of non-regular covariance matrix estimates. It is more flexible than both factor models and the graphical Lasso, which reflects itself in a higher likelihood value on test data at its optimal solution. A good solution, though, can only be found by choosing the regularization parameters carefully. In practice, this is mostly done by either grid, manual, or random search which do not come with any performance guarantees.

Here, we have modeled regularization parameter selection as a bi-level optimization problem. The lower level of this approach entails describing the solution gamut of all feasible solutions for all combinations of regularization parameters. For approximating the solution gamut, we have adapted and extended Benson's vector optimization algorithm. On data sets from different domains, our bi-level approach with Benson's algorithm for the lower level works better than the standard approach for solution gamut approximation for the latent variable graphical Lasso. Using Benson's algorithm on the lower level is a generic approach that works out-ofthe-box for any other problem that is given in the form of a regularized semidefinite program. Also, Benson's algorithm matches the known lower complexity bound for the solution gamut approximation problem.

Acknowledgements

This work has been supported by Carl-Zeiss-Stiftung and the DFG grant GI-711/3-2.