Verifying Global Optimality of Candidate Solutions to Polynomial
Optimization Problems using a Determinant Relaxation Hierarchy

Consider a polynomial optimization problem

where we use the multi-index notation $x^{\alpha}:=x_{1}^{\alpha_{1}}\cdots x_{n}^{\alpha_{n}}$ for $x\in{\mathbb{R}}^{n}$, $\alpha\in{\mathbb{N}}^{n}$, and where the data are polynomials $f,g_{1},\ldots,g_{m}\in{\mathbb{R}}[x]$ so that in the above sums only a finite number of coefficients $f_{\alpha}$ and $g_{i,\alpha}$ are nonzero. We will use the notation $|\alpha|:=\sum_{k=1}^{n}\alpha_{k}$.

Polynomial optimization problems (1) are generally non-convex and may therefore have multiple local and global optima. Local optima are feasible points for which no nearby point has a lower-cost objective value. Global optima are feasible points which have the lowest objective value among all feasible points. Many applications require (or at least benefit from) finding a global optimum to (1). Provably obtaining a global optimum is often accomplished by:

1. 1.

   Using a local solver to compute a locally optimal point that is a candidate for being a global optimum.

2. 2.

   Computing the optimality gap based on the difference between the objective value of the candidate point and the lower bound on the objective value provided by a convex relaxation of (1).

The candidate point is declared to be globally optimal if the optimality gap is sufficiently small. When the optimality gap is not sufficiently small, typical global optimization methods attempt to close the optimality gap by applying a variety of techniques for tightening the relaxation and obtaining a better local optimum [1, 2].

For certain polynomial optimization problems, especially those modeling physical systems, it is often straightforward to obtain high-quality local optima due to preexisting knowledge of appropriate initializations for local solvers. Consider, for instance, the optimal power flow (OPF) problem used in the operation of electric power systems [3, 4]. OPF problems are known to be NP-hard [5, 6]. Accordingly, there exist various test cases that challenge many solution algorithms [7, 8, 9, 10]. However, despite this challenging worst-case complexity, local solvers with physically justified initializations find the global optima for many practical OPF problems as verified by the small or zero optimality gaps obtained via a variety of convex relaxation techniques [11, 12, 13, 14, 15, 16]. (See [17] for a detailed survey.) Moreover, appropriately initialized local solvers can be significantly more computationally tractable than the convex relaxations used to certify global optimality [18, 19]. Thus, the computational effort for provably obtaining a globally optimal solution is often dominated by the time spent solving relaxations which show that a candidate feasible point from a local solver is, in fact, globally optimal.

Such problems motivate the development of methods for checking whether a given feasible point (e.g., the output of a local solver) is globally optimal. In contrast to convex relaxation techniques, whose objective values provide bounds on the potential suboptimality associated with a local solution, this paper solely focuses on methods for certifying whether a given feasible point is globally optimal. This less ambitious aim in comparison to existing convex relaxation techniques facilitates computational advantages for the proposed approach. Specifically, we reduce the problem of certifying global optimality of a candidate point to testing feasibility of a linear program.

Our proposed approach is particularly useful in settings that require the solution of many polynomial optimization problems. When applied as a screening step, our approach limits the need for explicitly solving computationally intensive relaxations to problems for which the candidate point and the selected relaxation result in a non-zero optimality gap. This allows us to leverage much of the extensive theoretical work regarding relaxations of polynomial optimization problems that has heretofore been practically limited by computational challenges associated with existing conic optimization solvers.

We close the introduction by reviewing two papers that aim to certify global optimality of polynomial optimization problems. Reviews of the specific relaxation techniques we employ will be presented in the following section.

The approach in [20] finds a objective function $\hat{f}$ for which a candidate point, $\hat{x}$, is globally optimal:

where $\mathcal{K}$ denotes the feasible space defined by the constraints in (1) and $\|\cdot\|$ is, for instance, the 2-norm of the coefficients. If the optimal objective value of (2) is equal to zero, then $\hat{x}$ is globally optimal. If not, then the solution of (2) provides insight into how far $\hat{x}$ is from the global optimum. Specifically, the solution to (2) yields an objective function, as close to the original objective function as possible, for which the candidate point $\hat{x}$ is globally optimal. In order to solve (2), [20] proposes to replace the non-negativity constraint using Putinar’s Positivstellensatz. The drawback of this approach is that it is a priori as computationally costly as actually computing a global solution.

In [21], another related approach provides a sufficient (but not necessary) condition for certifying global optimality of a candidate point by evaluating whether the point satisfies the Karush-Kuhn-Tucker (KKT) conditions [22, 23] for a relaxation. Assuming certain constraint qualification conditions (e.g., Slater’s condition), a local solution necessarily satisfies the KKT conditions for the non-convex polynomial optimization problem (1). If the local solution additionally satisfies the KKT conditions for a convex relaxation of (1), then the local solution is, in fact, a global optimum. The approach in [21] specifically tests whether a local solution to an OPF problem satisfies the KKT conditions for the Shor relaxation [24], which is equivalent, in the case of OPF problems, to the first-order moment/sum-of-squares relaxation in the Lasserre hierarchy to be discussed below. While explicitly computing the Shor relaxation requires the solution of a semidefinite program, evaluating whether the primal/dual variables obtained from a local solver satisfy the Shor relaxation’s KKT conditions only requires a single Cholesky factorization of a sparse matrix and a vector dot product. Thus, certifying global optimality by evaluating the condition in [21] is up to two orders of magnitude faster than explicitly solving the Shor relaxation for various OPF test cases, even after applying advanced techniques for improving the relaxation’s computational tractability [25, 26]. The approach proposed in the remainder of this paper can be viewed as an extension of the condition in [21] to exploit certain tighter relaxations, specifically, higher-order moment/sum-of-squares relaxations in the Lasserre hierarchy.

The remainder of paper is organized as follows. Section II presents the proposed approach. Section III applies the approach to several polynomial optimization problems including an instance of the optimal power flow problem. Section IV concludes the paper and proposes avenues for future research.

In 2001, the Lasserre hierarchy [27, 28] (see also [29, 30]) was proposed to find global solutions to polynomial optimization problems. It is also known as moment/sum-of-squares hierarchy in reference to the primal moment hierarchy and the dual sum-of-squares hierarchy. Its global convergence is guaranteed by Putinar’s Positivstellensatz [31] proven in 1993. Typically, if one of the constraints is a ball $x_{1}^{2}+\ldots+x_{n}^{2}\leqslant 1$, then the sequence of lower bounds provided by the hierarchy converges to the global infimum of the polynomial optimization problem. In addition, there is zero duality at all relaxation orders [32].

The moment problem of order $d$ is defined as

where $\succcurlyeq 0$ denotes positive semidefiniteness, and where the Riesz functional, the moment matrix, and the localizing matrices are respectively defined by

Above, $\lceil\cdot\rceil$ denotes the ceiling of a real number. Note that the relaxation order $d$ must be greater than or equal to $d^{\text{min}}:=\max_{i}k_{i}$. See Section III for an illustrative example of a second-order moment/sum-of-squares relaxation of a quadratically constrained quadratic program corresponding to a two-bus OPF problem.

Our approach extends a relaxation technique proposed in [33] that leverages theory developed in [34, 35] regarding the importance of the convexity of the feasible space rather than its algebraic representation.¹¹1In particular, under a mild degeneracy condition, log-barrier solution methods are guaranteed to converge to a global minimizer of a problem with a convex feasible space even if the algebraic representation of the feasible space is non-convex. The technique in [33] further relaxes the Shor relaxation [24] using necessary conditions for positive semidefiniteness based on matrix determinants. Specifically, the approach in [33] exploits the fact that a symmetric matrix is positive semidefinite if and only if all its principal minors (i.e., the determinants of its principal submatrices) are non-negative [36]. Extending this approach to the Lasserre hierarchy (3) yields

where $n_{0}:=\frac{\left(n+d\right)!}{n!\,d!}$ is the size of the matrix $M_{d}(y)$, $n_{i}:=\frac{\left(n+d-k_{i}\right)!}{n!\,\left(d-k_{i}\right)!}$ is the size of the matrix $M_{d-k_{i}}(y)$, and $\det_{I}$ denotes the principal minor corresponding to the index $I$ (i.e., the determinant of the submatrix whose rows and columns are indexed by $I$), likewise for $\det_{J_{i}}$.

To derive our approach, we exploit the KKT conditions of the formulation (5). To express these conditions, let $\lambda$ denote the KKT multiplier associated with the constraint $y_{0}=1$. Let $\lambda_{0,I}$ denote the KKT multiplier associated with the constraint $\det_{I}(M_{d}(y))\geqslant 0$. Let $\lambda_{i,J_{i}}$ denote the KKT multiplier associated with the constraint $\det_{J_{i}}(M_{d-k_{i}}(g_{i}y))\geqslant 0$. The KKT conditions of (5) are

We have that $\nabla\det(M)=\text{co}(M)$ where the comatrix $\text{co}(M)$ is defined by $\text{co}(M)_{i,j}:=(-1)^{i+j}\det(M^{(i,j)})$; $M^{(i,j)}$ is the matrix $M$ where row $i$ and column $j$ have been removed. Let $\text{co}_{I}(\cdot)$ denote the comatrix of the submatrix whose rows and columns are indexed by a subset of indices $I$. Below, we will consider matrices $B_{\alpha}$ and $C_{i,\alpha}$ such that the moment and localizing matrices are as follows:

and

The KKT conditions are then²²2We use the fact that for all sequences of multi-indexed real numbers $(h_{\alpha})_{|\alpha|\leqslant 2d}$, we have $d(\det\circ M_{d})(y).h=d(\det)(M_{d}(y))(d(M_{d})(y).h)=d(\det)(M_{d}(y)).M_{d}(h)=\text{trace}[\text{co}(M_{d}(y))M_{d}(h)]=\text{trace}[\text{co}(M_{d}(y))\sum_{\alpha}B_{\alpha}h_{\alpha}]=\sum\limits_{\alpha}\text{trace}[\text{co}(M_{d}(y))B_{\alpha}]h_{\alpha}$.

If we substitute $y$ using the feasible point $\hat{x}$, i.e., by taking $\hat{y}:=(\hat{x}^{\alpha})_{|\alpha|\leqslant 2d}$, the KKT conditions reduce to a linear feasibility problem, that is

Observe that all the determinants are equal to zero, except when the index sets $I$ and $J_{i}$ are singletons (i.e., $I=J_{i}=\{\alpha\}$ where $\alpha\in\mathbb{N}^{n}$). In that case, the determinant is simply equal to a diagonal term of the matrix.

Any point that satisfies (10) provides values for dual variables that, in combination with the primal variables $\hat{y}$ obtained from the candidate solution $\hat{x}$, satisfy the KKT conditions corresponding to a convex relaxation of (1). Hence, any feasible point for (10) certifies global optimality of the candidate feasible point $\hat{x}$ for the polynomial optimization problem (1).³³3Observe that neither [34, Assumption 2.1] nor Slater’s condition are required in order to show that every KKT point is a minimizer in [34, Theorem 2.3]. In other words, all assumptions needed for the theory in [34] to be applicable for our purposes are satisfied in our setting.

Since all values $\hat{y}$ are fixed, we emphasize that (10) is linear in the KKT multipliers $\lambda$, $\lambda_{0,I}$, and $\lambda_{0,J_{i}}$. Thus, the KKT conditions (10) can be verified using a linear programming solver, which is significantly more computationally tractable than solving the KKT conditions for the relaxation (9) with both primal variables $y$ and dual variables $\lambda$, $\lambda_{0,I}$, and $\lambda_{0,J_{i}}$ allowed to vary (a complementarity problem).

We implemented the algorithms with Matlab R2020a. We used Gloptipoly 3 to obtain the Lasserre hierarchy formulation. As for solving the KKT conditions, we used the LINPROG and LSQLIN solver provided by Matlab. The optimality tolerance of both solvers are set to $10^{-8}$.

We have proposed a method for verifying global optimality by checking for the feasibility of a linear system of equations. We have demonstrated the proof of concept on several small examples of polynomial optimization problems. A topic of further investigation is to exploit sparsity in order to make the approach tractable on larger optimal power flow instances. One could use techniques involving chordal graphs as was done in [15], which enabled the application of the Lasserre hierarchy to large-scale instances. Since global solutions were extracted in many instances within minutes, it is reasonable to believe that checking for global optimality should be possible with significantly faster runtimes. Another possible avenue for future research is to find another equivalent formulation of the Lasserre hierarchy whose KKT conditions would also lead to a tractable system for the dual multipliers (perhaps a second-order conic program).