Entropy-Penalized Semidefinite Programming

1 Introduction

Semidefinite programming (SDP) has become a key tool in solving numerous problems across operations research, machine learning, and artificial intelligence. While there are too many applications of SDP to present even a representative sample, inference in graphical models [Wainwright and Jordan, 2008; Erdogdu et al., 2017], multi-camera computer vision [Torr, 2003], and applications of polynomial optimization [Parrilo, 2003; Lasserre, 2015] in power systems [Ghaddar et al., 2016] stand out. Under some assumptions [Madani et al., 2014], the rank at the optimum of the SDP relaxation is bounded from above by the tree-width of a certain hypergraph, plus one. When a rank-one solution is not available, it is often not needed [Mareˇcek and Takáˇc, 2017], as one should like to construct a stronger SDP relaxation.

Penalization of the objective is a popular approach for obtaining low-rank solutions, at least in theory [Recht et al., 2010; Lemon et al., 2016; Zhou, 2019; Fawzi et al., 2019]. Notice that without a further penalization, an interior-point method for SDP provides a solution on the boundary of the feasible set, where SDP corresponds to the optimum of the highest rank, whenever there are optima of multiple ranks available. The use of a penalization provides a counter-balance in this respect. In practice, however, the penalties are often ignored, as it is believed that their computation is too demanding for large-scale problems and does not guarantee low-rank solutions in general.

∗Contact Author

1 See https://github.com/mkrechetov/epsdp for the implementation details.

An alternative approach develops numerical optimization methods that seek a priori low-rank solutions. This approach, widely attributed to [Burer and Monteiro, 2003], considers a factorization of a semidefinite matricial variable X = V · V ⊤ with V ∈ Rn×k for increasing 1 ≤ k ≪ n. In general, the resulting problems are non-convex. Early analyses required determinant-based penalty terms [Burer and Monteiro, 2005], although no efficient implementations were known. Under mild assumptions, for large-enough k, there is a unique optimum over such a factorization even without a penalization and it recovers the optimum of the initial SDP problem [Boumal et al., 2016]. For smaller values of k, it is known that the low-rank relaxation achieves O(1/k) relative error [Mei et al., 2017]. Much more elaborate analyses [Erdogdu et al., 2018] are now available. Especially when combined with efficient gradient computation, e.g. within low-rank coordinate descent [Mareˇcek and Takáˇc, 2017, e.g.], this approach can tackle sufficiently large instances and is increasingly popular.

In this paper, we aim to develop a method combining both approaches, i.e., utilize an efficient low-rank-promoting penalty in the Burer-Monteiro approach. We present efficient first-order numerical algorithms for solving the resulting penalized problem, with (almost) linear-time per-iteration complexity. This makes the combined approach applicable to a wide range of practical problems.

In a case study, we focus on certain combinatorial optimization problems and inference in graphical models. We show that despite the non-convexity of the penalized problem, our approach successfully recovers rank-one solutions in practice. We compare our solutions against non-penalized SDP, belief propagation, and state-of-the-art branch-and-bound solvers of [Krislock et al., 2014].

Contribution. Our contributions can be summarized as follows. We show

1. convergence properties of optimization methods employing a wide class of penalty functions that promote low-

rank solutions;

1. linear-time algorithms for computing the gradient of these penalty functions;
2. computational results on the penalized SDP relaxation of maximum a posteriori (MAP) estimates in Markov random fields (MRF), which considerably improve upon the results obtained by interior-point methods and randomized rounding.

This allows for both well-performing and easy-to-analyze low-rank methods for SDPs coming from graphical models, combinatorial optimization, and machine learning.

Paper structure. This paper is organized as follows. First, we define the conic optimization problem together with a penalized form with a list of suitable penalization functions. Next, we present theoretical guarantees for solution recovery. These extend known results for solution recovery to the penalty case. Then, we consider the MAP problem in Markov Random Fields (MRF) and introduce an iterative procedure for it, together with a first-order method for solving a subproblem at each step; we also show how to compute the gradients efficiently. Finally, we provide computational experiments for different inference problems in MRF.

2 Background

SDP is the following conic optimization problem:

where X ∈ Sn+ denotes that the n × n matrix variable W is symmetric positive semidefinite, I and J are finite index sets, each fi and gj are convex functions Rn×n → R, and Cj ∈ Rn×n and Si ∈ Rn×n are constant matrices.

In the context of combinatorial optimization, one may also consider even more powerful methods such as Sum-of-Squares hierarchies of [Parrilo, 2003]. However, even an SDP relaxation, which is, in fact, the first step of this hierarchy, may be too computationally challenging. It is usually solved by interior-point methods in the dimension that is quadratic in the number of variables and thus becomes intractable even for medium-scale problems (with a few thousand variables). The problem becomes even less scalable for higher orders of the hierarchy since it requires one to solve SDP with nΘ(d) variables, where d is the level of the hierarchy.

2.1 Low-rank Relaxations and Penalized Problem

First, let us formally define our notion of a penalty function and explain related work on first-order methods for SDP.

Assumption 1. Eq. (SDP) has an optimum solution with rank r.

Let us consider the following proxy problem:

Where q > r and Rq,λ(V ) satisfies the following:

Definition 1 (Strict penalty function). A function Rq,λ(V ) : Rn×n → R is a penalty function that promotes low-rank solutions if for some integers q′ ≤ q and a multiplier λ ∈ R+:

Moreover, if q′ = q and the rank (X) < q, then Rq,λ(V ) = 0, ∀λ, R(V ) is a strict penalty function.

We use the word penalty instead of penalty function that promotes low-rank solutions, where there is no risk of confusion. This notion of a penalty is rather wide. When multiplied by λ, a determinant is a prime example. One may also consider functions of the following quasi-norms.

1. The nuclear norm:

where σi is the i-th singular value, cf. [Lemon et al., 2016]. The norm is also known as a trace norm, Schatten 1-norm, and Ky Fan norm. As shown by [Srebro et al., 2004], in the method of [Burer and Monteiro, 2003], one can benefit from a bi-Frobenius reformulation:

There are also truncated [Hu et al., 2013] and capped variants [Sun et al., 2013].

1. Schatten-p quasi-norm for p > 0:

where σi(X) denotes the i-th singular value of X.

1. A smoothed variant of Schatten-p quasi-norm by [Pogodin et al., 2017] for p, ϵ > 0:

which is also the Schatten-1/3 quasi-norm.

1. Bi-nuclear (BiN) quasi-norm of [Shang et al., 2018]:

which is also the Schatten-1/2 quasi-norm.

1. Frobenius/nuclear quasi-norm of [Shang et al., 2018]:

which is also the Schatten-2/3 quasi-norm.

We also note there has been considerable interest in the analysis of low-rank approaches without penalization, especially in matrix-completion applications. Much of the analysis goes back to the work of [Keshavan et al., 2010]. For further important contributions, see [Arora et al., 2012].

2.2 Entropy Viewpoint

One could see the penalty functions introduced above from the entropy-penalization perspective. This is useful not only from a methodological standpoint, but also from a computational one. To this end, we consider the Tsallis entropy:

The Tsallis entropy is crucial in our study because it generalizes many popular penalties considered earlier. The Schatten p-norm coincides with the Tsallis entropy S Tp over a set of matrices with a fixed trace norm, so that the tri-quasi norm and bi-nuclear norm (2-6) are covered as well. The Log-Det function, - log detX, which is also used in low-rank SDP, is up to an additive constant factor relative (Shannon) entropy taken concerning a unit matrix, while Renyi (S R) and von Neumann (S N ) entropies,

respectively, can also be used as penalties to promote a lowrank solution. To the best of our knowledge, neither Renyi, von Neumann, nor Tsallis entropies have been studied in the context of low-rank SDP.

3 Exact Recovery

Let us now present a unified view of the penalties and their properties:

Lemma 1. Any of:

1. λdet(X);
2. λσq(X), where σi(X) denotes the i-th singular value of X;
3. Tsallis, Renyi, and von Neumann entropies defined on the last n - q + 1 singular values;

is a penalty function that promotes low-rank solutions. Moreover, penalties 1-3 are strict.

Proof. Sketch. (1.) The proof is by simple algebra. (2.) If σq(X) is 0, we know the rank is q - 1 or less. Otherwise, for large values of λ, the value of the penalties goes to infinity, and hence q′ = q. (3.) The definition of entropy assumes that S(0, ..., 0) = 0, thus all entropies are strict penalty functions by definition. (4.) First, consider the case where all non-zero singular values are equal. In that case, ‖X‖∗/σmin(X) = rank (X), and subtracting q results either in a non-positive number when the rank is less than q or a positive number otherwise. If the singular values are non-equal, ‖X‖∗/σmin(X) provides an upper bound on the rank of X, which can be improved as suggested. The use of the upper bound results in the value of the penalty tending to infinity for values between q′ and q in the large limit of λ.

Crucially, under mild assumptions, any penalty allows for the recovery of the optimum of a feasible instance of (SDP) from the iterations of an algorithm on the non-convex problem in variable V ∈ Rn×r , such as in the methods of [Frostig et al., 2014] or [Burer and Monteiro, 2003]. In contrast to the traditional results of [Burer and Monteiro, 2005], who consider the det penalty, we allow for the use of any strict penalty function.

Theorem 2. Assume that we solve the proxy problem (P-SDP) iteratively and Rq,λ(V ) is a strict penalty function that promotes low-rank solutions. In each iteration, if Rq,λ(V ) ̸= 0, we increase λ (e.g., set λt+1 = γλt, with u > 1 as some fixed parameter). Furthermore, let us assume that the solution we found is denoted by ˜ Vq with rank ( ˜ Vq) = q′ < q. Let us also denote ˜ Vq′ ∈ Rn×q′ some factorization of ˜ Vq ˜ V ⊤ q (such factorization exists because rank ( ˜ Vq) = q′ ). Also assume that we have an optimal solution of (SDP), X∗ with a rank (X∗ ) = r. If

is a local minimum of Pq′+r+1,λ, then (˜ Vq′ ) ˜ V ⊤ q′ is a global solution of SDP.

Proof. Let us define a family of matrices for τ ∈ [0, 1] as follows:

where (V∗)⊤(V∗) is some factorization of X∗ with V∗ ∈ Rn×r .

Note that ∀τ , we have rank (X(τ )) < r+q′ +1, and hence ∀λ, τ : Rq′+r+1,λ(V (τ)) = 0. Now, assume the contradiction, that is, Vq′+r+1 is a local optimum solution but ˜ Vq′ is not a global solution.

We show that ∀τ ∈ [0, 1], V (τ) is a feasible solution. Indeed, for any j ∈ J we have

We just showed that for each V (τ), τ ∈ [0, 1] is a feasible point. Now, let us compute the objective value at this point. For all τ ∈ [0, 1], we have

which is a contradiction under the assumption that ˜ Vq is a local optimum.

4 Efficient Implementation for MAP Inference

A pairwise Markov Random Field (MRF) is defined for an arbitrary graph G = (V, E) with n vertices. We associate a binary variable xi ∈ {-1, +1} with each vertex i ∈ V . Let θi : {±1} → R and θij : {±1}2 → R defined for each vertex and edge of the graph be vertex and pairwise potential, respectively. Thus, a posteriori distribution of x follows the Gibbs distribution:p(x|θ) = 1

with U(x; θ) = ∑i∈V θi(xi) + ∑(i,j)∈E θij(xi, xj). The maximum a posteriori (MAP) estimate is then

which is its turn an NP-hard binary quadratic optimization problem,

with indefinite matrix S. The SDP relaxation for this problem is given by [Goemans and Williamson, 1995; Nesterov, 1998]:

which also covers the Ising model in statistical physics and a number of combinatorial optimization problems. We believe that the approach can be extended to a general setup given by Eq. (SDP).

An entropy-penalized SDP relaxation of (8) has the form

where V i is the i-th column of matrix V ∈ Rn×k , X = V V ⊤.

4.1 Numerical Method.

To solve Problem (EP-SDP), we use the Augmented Lagrangian method starting from a sufficiently small value of the penalty parameter λ > 0 and increasing it in geometric progression, λk+1 = λkγ, with γ > 1, as summarized in Algorithm ??. The efficiency of the method is due to the efficient computability of gradients of Tsallis, Renyi, and von Neumann entropies:

Lemma 3. For any matrix V ∈ Rn×k with k = O(1), let X(V ) = V ⊤V . Then, gradients of S Tα (X), S R(X), and S N (X) can be computed in O(n) time. Moreover, if the number of non-zero elements in matrix A is O(n), then the iteration complexity of Algorithm ?? is O(n).

Proof. We start our analysis with Tsallis entropy. First, compute the gradient of S Tα in V :

Similarly for Renyi, S R(X), and von Neumann, S N (X), entropies we have

Algorithm 1: Entropy-Penalized SDP.

Data: Quadratic matrix S of the MAP inference problem, staring point λ0, γ > 1, step size policy {ηk}k≥1 accuracy parameters ε, ϵ

Result: Solution V∗ as a local minimum of (EP-SDP) of

unit rank

begin

Return: first singular vector of Vt.;

and

Following [Holmes et al., 2007], the singular-value decomposition of matrix V = U1DU2 with U1 ∈ Rn×n , D ∈ Rn×k , and U2 ∈ Rk×k can be performed in O(min (nk2 , n2k)) = O(nk2 ) time.

For any α > 1, the product Xα-1 · V = U1D2α-1U2 can be computed in time O(n) together with trXα = trD2α and trX = trD2 . Thus, for a fixed k, the gradient ∂ S T (X) ∂V computation time is linear in its dimension. (Here, for any α ∈ (0, 1), we use the identity ∂λi = vi⊤∂Xvi.) To finish the proof of the statement, it remains to note that matrix-vector multiplication takes O(n) time for any matrix with O(n) nonzero entries.

5 Case Study

In this section, we compare our penalized algorithm with other conventional approaches to MAP problems. We fix the width of factorization to k = 10, since there is no significant gain in practice for larger values of k, cf. [Mei et al., 2017]. We choose 2ηkβ = 1, where β is the Lipschitz constant of the gradient in ℓ2 norm and γ = 3/2. Parameters λ0 and γ of Algorithm ?? are usually chosen by a few iterations of random search. It is usually enough to have about 35 iterations for penalty updates and a few hundred iterations to find a local minimum using Algorithm ??. We emphasize that matrices we obtain by solving EP-SDP are rank-one solutions on all MAP instances presented. Thus, we do not need any further rounding procedure.

First, in Table 1, we show the performance of our algorithm on selected hard MAP inference problems from the BiqMac collection 2 . We selected a few of the hardest instances ("gkaif" among them)-dense quadratic binary problems of 500 variables.

We compared our algorithm (EP-SDP with Tsallis entropy and α = 2) with the plain-vanilla semidefinite programming instance solved by the interior-point method, possibly with rounding using the best of one thousand roundings of [Goemans and Williamson, 1995] and also with Gurobi solver for Mixed-Integer Problems. To avoid any confusion, we solve the corresponding maximization problems; by the objective value, we mean the value at a feasible solution produced by the method (e.g., rounded solution of SDP relaxation), which is a lower bound for the corresponding problem. Because these problems are of the same size (but varying density), the running time of each method is almost constant. It took around 10 minutes for CVXPY to solve the SDP relaxation, and it runs out of memory for the two problems with higher density. Within five seconds, EP-SDP obtains results that are better than what Gurobi can produce in 70 seconds.

In our study, parameter α of entropies S T α , S N α , and S R α is chosen on an exponential grid from 1 to 10 with a step 1.1. After experimentation, we note that α = 1.1 and α = 5.0 seem to improve the results the best for the low-rank SDP with Tsallis and Renyi entropies, respectively, although the difference between different α ∈ (1, 10) is not very significant for either of the (Tsallis and Renyi) entropies.

Table 2 summarizes the results of solving the Max-Cut problem over a GSET collection of sparse graphs 3 . As we see from the experiments, the results of applying suitable entropy often outperform both the plain-vanilla SDP with the classical Goermans-Williamson rounding, the mean-field approximation, as well as the results of UGM solver 4 for loopy belief propagation and mean-field inference. It is worth noting that for several instances of the GSET graph collection, loopy belief propagation provides rather weak results. Usually, strong results of the loopy belief propagation are complementary to those of the mean-field approximation, which is supported by our empirical results. Results of both loopy belief propagation and mean-field approximation can be substantially improved using the linear-programming belief-propagation approach (LP-BP).

2 http://biqmac.uni-klu.ac.at/biqmaclib.html

3 https://sparse.tamu.edu/Gset

4 https://www.cs.ubc.ca/ schmidtm/Software/UGM.html

We also want to point out that our iterative algorithm successfully decreases the rank of the solution. The higher the penalization parameter, the lower the rank. We illustrate this in Figure 1, where for Tsallis entropies with α = 2 and α = 5, we plot the second singular value of matrix V . For this plot, we considered the Max-Cut problem for the first graph from the GSET collection and the Tsallis entropy as the penalization function. In Figure 2, we illustrate the same concept for fixed penalization (Tsallis entropy with α = 2) and different initial values of the multiplier λ. We observe that for different penalization functions and update schemes, the rank of the solution decreases gradually with each step. In practice, our iterative algorithm could be seen as a universal rounding procedure for SDP relaxations. Indeed, if we choose a large-enough penalization update (e.g., γ = 2 as in Figure 1), we easily obtain a rank-one solution that is not worse and often is substantially better than solutions obtained by randomized rounding.

Overall, we would like to stress that Algorithm ?? is very fast. This is shown in Figure 3 and Table 3, where we compare run times of EP-SDP, low-rank Burer-Monteiro approach (LR-SDP), and interior-point method solvers (SDP) for various Erdos-Renyi random graphs. From the data, we see that (assuming the fixed width of factorization k = 10) EP-SDP run time increases linearly with the number of vertices. Indeed, throughout the benchmark instances tested, the run time does not exceed a few seconds per each of the test cases. At the same time, the bound is often almost as good as that of the Branch and Bound Biq-Mac Solver of [Krislock et al., 2014], which requires a significant amount of time.

6 Conclusions

This paper presented a unified view of the penalty functions used in low-rank semidefinite programming using entropy as a penalty. This makes it possible to find a low-rank optimum, where there are optima of multiple ranks. Semidefinite programs with an entropy penalty can be solved efficiently using first-order optimization methods with linear-time per-iteration complexity, which makes them applicable to large-scale problems that appear in machine learning and polynomial optimization. Our case study illustrated the practical efficiency on binary MAP inference problems. The next step in this direction is to consider the structure of the SDP, which seems to be crucial for further scalability.

Acknowledgements

The work of Martin Takáˇc was partially supported by the U.S. National Science Foundation, under award numbers NSF:CCF:1618717, NSF:CMMI:1663256, and NSF:CCF:1740796. The work at LANL was supported by the U.S. Department of Energy through the Los Alamos National Laboratory as part of the GMLC "Emergency Monitoring and Control through new Technologies and Analytics" project. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of U.S. Department of Energy (Contract No. 89233218CNA000001).