Convex Minimization with Integer Minima in $\widetilde{O}(n^{4})$ Time^††thanks: A preliminary version of this paper appeared at the 35th Annual ACM-SIAM Symposium on Discrete Algorithms (SODA 2024).

The main result of this paper is a much more efficient algorithm for minimizing convex functions with integer minimizers:

The strong assumption that all extreme points of $K^{*}$ are integer vectors guarantees that the algorithm outputs an integer minimizer of $f$. In fact, this assumption is necessary for the algorithm to efficiently compute an integer minimizer (see Remark 1.4 in [32] for an example). Without such an assumption, the algorithm can still output a minimizer, though not neccessarily an integer one, with the same guarantee as in Theorem 1.1 (see Remark 1.5 in [32]).

Previously, under the same assumptions as in Theorem 1.1, Jiang [32] gave an algorithm that computes an integer minimizer of $f$ using $O(n(n\log\log n/\log n+\log R))$ separation oracle calls, which is better than the oracle complexity in Theorem 1.1. However, his algorithm uses an enormous $\widetilde{O}(n^{8})$ number of additional arithmetic operations because of the following reasons: to obtain a subquadratic oracle complexity, [32] resorts to the random walk based cutting plane method by Bertsimas and Vempala [6] to work directly with the polytope $K$ formed by all the separating hyperplanes output by $\operatorname{\mathsf{SO}}$. Unfortunately, the sampling approach in [6] requires $\Omega(n^{6})$ arithmetic operations to perform one iteration of the cutting plane method. In addition, to guarantee that the sampling can be performed efficiently, [32] has to maintain two “sandwiching” ellipsoids, contained in and containing the polytope $K$ respectively, which requires additional computational cost to preserve after reducing the dimension of the problem.

Another computational bottleneck in [32] is that the algorithm approximates the length of the shortest vector after every iteration of the cutting plane method. To do so, [32] leverages the sieve algorithm in [2] that computes a $2^{n\log\log n/\log n}$-approximation to the shortest vector using $\widetilde{O}(n^{4})$ arithmetic operations. Compounding with the overall $\widetilde{O}(n^{2})$ iterations of the algorithm, this step of approximating the shortest vector takes $\Omega(n^{6})$ arithmetic operations in total.

To get around the aforementioned computational bottlenecks, we utilize an intricate combination of a computationally more efficient cutting plane method based on the volumetric center of $K$ due to Vaidya [67] and a faster version of the LLL algorithm given in [51]. Instead of working directly with the polytope $K$ and performing volume reduction in a step-by-step manner, we run Vaidya’s cutting plane method [67] in blocks of $O(n\log n)$ consecutive steps, the volume decreases by a constant factor per step only in an amortized sense within each block. In particular, recent work by Jiang, Lee, Song and Wong [34] shows that $O(n\log n)$ steps of Vaidya’s method can be implemented using a total of $O(n^{3}\log n)$ arithmetic operations. Running Vaidya’s method in blocks also enables us to compute only $\widetilde{O}(n)$ approximate shortest vectors in total, while still being able to avoid the appearance of extremely short vectors. Harnessing the faster implementation of the LLL algorithm in [51] allows us to compute the $\widetilde{O}(n)$ approximate shortest vectors using a total of $\widetilde{O}(n^{4})$ arithmetic operations.

Submodular function minimization (SFM) has been one of the cornerstone problems in combinatorial optimization since the seminal work of Edmonds [19]. Many classic combinatorial problems can be abstracted as optimization over a submodular function, such as graph cut function, set coverage function and economic utility function. For more comprehensive reviews of SFM, we refer readers to [46, 30].

We briefly recall the standard setting of SFM: we are given a submodular function $f:2^{V}\rightarrow\operatorname{\mathbb{Z}}$ over an $n$-element ground set $V$ where the function $f$ is accessed via an evaluation oracle, which returns the value $f(S)$ for a query set $S$ using time $\mathsf{EO}$. The goal is to find the minimizer of $f$ using queries to the evaluation oracle and additional arithmetic operations.

An SFM algorithm is called strongly-polynomial time if its runtime depends polynomially only on $\mathsf{EO}$ and the dimension $n$, but not on the range of the function $f$. In their seminal work, Grötschel, Lovász, and Schrijver [23, 24] gave the first strongly-polynomial time algorithm for SFM based on the ellipsoid method. Since then, there has been a long history of efforts in designing faster strongly-polynomial time algorithms for SFM (see Table 1).

The state-of-the-art strongly-polynomial time algorithms for SFM, in terms of the number of arithmetic operations used, were given by Lee, Sidford, and Wong [44] and Dadush, Végh, and Zambelli [17]. Both algorithms have runtime $O(n^{3}\log^{2}n\cdot\mathsf{EO}+n^{4}\log^{O(1)}n)$. The oracle complexity of these algorithms were later improved to sub-cubic in terms of $n$ by Jiang [31, 32], where he gave a strongly-polynomial time algorithm with runtime $O(n^{3}\log\log n/\log n\cdot\mathsf{EO}+n^{8}\log^{O(1)}n)$ alongside an exponential time algorithm with a nearly-quadratic $O(n^{2}\log n)$ oracle complexity. It remains a major open problem in the area of SFM whether there exists a strongly-polynomial time algorithm with truly sub-cubic, i.e. $O(n^{3-c})$ for some absolute constant $c>0$, oracle complexity. When the problem has certain structures, such as the minimizer is $k$-sparse, Graur, Jiang and Sidford [21] have developed an algorithm with $\widetilde{O}(|V|\cdot\operatorname{poly}(k))$ queries. SFM has also been heavily studied in the parallel setting, see e.g., [8].

Unfortunately, while the algorithm in [32] answers major open questions in [44] and makes significant progress on the oracle complexity for SFM, the number of additional arithmetic operations used by Jiang’s algorithm is a factor of $\widetilde{\Theta}(n^{4})$ larger than the algorithms in [44, 17]. In this paper, we complement Jiang’s algorithm by giving the following strongly-polynomial time SFM algorithm that improves both the oracle complexity and the number of arithmetic operations used by the algorithms in [44, 17].

Corollary 1.2 almost immediately follows from Theorem 1.1 together with the standard fact that a separation oracle for the Lovász extension of the submodular function $f$ can be implemented by making $n$ queries to the evaluation oracle of $f$ and $O(n\log n)$ additional arithmetic operations (e.g., Theorem 6.4 in [32]). The proof of Corollary 1.2 is essentially identical to the proof of Theorem 1.7 in [32], and is therefore omitted.

Compared to the previous best runtime algorithms in [44] and [17], our algorithm improves their oracle complexity from $O(n^{3}\log^{2}n)$ to $O(n^{3}\log n)$ while also improving the number of arithmetic operations from $O(n^{4}\log^{O(1)}n)$ to $O(n^{4}\log n)$. We highlight that in [44], the authors manage to achieve an $O(n^{3}\log n)$ oracle complexity, but at the expense of an exponential runtime. It is also important to note that our improvements are achieved via a much more general algorithm, whereas the algorithms in [44, 17] work specifically for SFM.

In comparison to the current best oracle complexity algorithm in [32], our algorithm has a slightly worse oracle complexity, but we significantly improve the $\widetilde{O}(n^{8})$ additional arithmetic operations in his algorithm down to $O(n^{4}\log n)$.

Finally, note that the current best implementation of a separation oracle for the Lovász extension requires $n$ queries to $\mathsf{EO}$, and the current fastest cutting plane method requires $O(n^{2})$ arithmetic operations per step. So for any cutting plane algorithm for SFM that uses $T$ iterations, the current best runtime we can hope for such a method is $O(Tn\cdot\mathsf{EO}+Tn^{2})$ using state-of-the-art techniques. Our algorithm, in fact, matches such a runtime bound. In particular, we use $O(n^{2}\log n)$ iterations of the cutting plane method with a total runtime of $O(n^{3}\log n\cdot\mathsf{EO}+n^{4}\log n)$. So in some sense, our result matches what can be achieved using the current best known algorithmic techniques for cutting plane methods. In contrast, the algorithm in [32] does not have such a feature, with the $\Omega(n^{8})$ additional arithmetic operations being much larger than its oracle complexity.

Before discussing our main technical insights, we first review the approach of [32] that achieves subquadratic oracle calls for minimizing convex functions with integer minimizers. The key ingredient in Jiang’s algorithm is an interplay between the polytope, formed by the separating hyperplanes returned by the separation oracle, which contains the integer minimizer, and a lattice that captures the structures of integer points.

In particular, in each iteration of the cutting plane method, Jiang’s algorithm examines the length of the (approximate) shortest vector in the lattice under the ${\rm Cov}(K)$-norm, where ${\rm Cov}(K)$ is the covariance matrix of the polytope $K$. The length of this vector under ${\rm Cov}(K)$-norm provides a measurement for the width of the outer ellipsoid induced by ${\rm Cov}(K)^{-1}$. Thus, when this norm is small, it implies that the outer ellipsoid has a small width and we can safely proceed by reducing dimension and projecting the lattice onto the orthogonal complement of the shortest vector. On the other hand, if the shortest vector still has a relatively large norm, then the polytope $K$ can be further refined by using cutting plane method. To achieve the best possible oracle complexity, the algorithm in [32] performs the width measurement in a step-by-step fashion, meaning that it actively checks the length of the shortest vector after every single cutting plane step. This ensures that the algorithm can enter the dimension reduction phase as soon as a short lattice vector can be obtained, which happens when the volume of $K$ is small enough.

Unfortunately, such a step-by-step strategy and careful width measurement come at a price – the approximate shortest vector subprocedure will be called in every iteration of the cutting plane method. Moreover, an expensive cutting plane method, such as the random walk based center of gravity method [6], that guarantees the volume shrinks in every step needs to be used.

Now we discuss our approach for overcoming the major computational bottleneck of [32], which is the step-by-step measurement of width mentioned above. For simplicity, we assume the integer minimizer $x^{*}\in\{0,1\}^{n}$ (i.e. $R=1$) in the subsequent discussion. By delaying the width measurement, i.e. find an approximate shortest vector after a block of cutting plane steps, we can utilize the much more efficient cutting plane method due to Vaidya [67]. However, this comes at the cost of a possibly much larger number of oracle calls, since the length of the shortest vector in the lattice might become very small when we actually perform the measurement, which will lead to a large increase in the volume of $K$ after dimension reduction. Nevertheless, we show that by carefully balancing the block size of cutting plane steps and the loss incurred during dimension reduction due to short lattice vectors, we can still achieve an $O(n^{2}\log n)$ oracle complexity.

To explain our approach, we first provide a brief introduction to Vaidya’s cutting plane method. Unlike classical cutting plane methods, such as the ellipsoid method and the center of gravity method, where each step is guaranteed to shrink the volume, Vaidya’s algorithm and analysis rely on controlling the volume of the Dikin ellipsoid induced by the log barrier on polytope $K$. The algorithm iteratively finds the point $\omega_{K}$ whose corresponding Dikin ellipsoid has the largest volume, which is called the volumetric center. Subsuquently, the algorithm uses a separation oracle on the volumetric center $\omega_{K}$. If $\omega_{K}$ is a minimizer, then we are done; otherwise, the algorithm computes the leverage score of each constraint with respect to the Hessian matrix, which measures the relative importance of each constraint. If all constraints are relatively important, then a new constraint is added based on the separating hyperplane returned by the separation oracle at $\omega_{K}$. If one of the constraints has leverage score smaller than some tiny constant, then it is dropped to ensure that the polytope $K$ always has at most $O(n)$ constraints. Due to the increase in the volume of the polytope $K$ when constraints are dropped, the volume of $K$ only shrinks in an amortized sense. In particular, only after $O(n\log n)$ steps, the volume of $K$ is guaranteed to decrease by a multiplicative factor of $2^{-O(n\log n)}$ from the initial volume. It is crucial to note that no guarantee on the volume shrinkage of $K$ is known if $o(n\log n)$ iterations of Vaidya’s cutting plane method is run.

Our idea is then to view the $O(n\log n)$ steps as a “block”, meaning that each time we measure the shortest vector and realize that it still has a large norm, we execute Vaidya’s cutting plane steps for $O(n\log n)$ steps and then re-examine the norm. This naturally induces a strategy that measures the width of the ellipsoid in a lazy fashion. The volume decrease of $K$ via this strategy is similar to that when running the cutting plane method in [6] for $O(n\log n)$ steps. Therefore, the only issue is the significant decrease in the norm of the shortest vector during a block of cutting plane steps, which can be as small as $2^{-O(n\log n)}$. Fortunately, this shrinkage of the norm is acceptable and it incurs only an additional $O(\log n)$ factor in the number of oracle calls.

One major advantage of using Vaidya’s cutting plane method is that the iterations can be performed very efficiently. Using the fast implementation in [34], a block of $O(n\log n)$ cutting plane steps can be done in $O(n^{3}\log n)$ time. Lazy width measurement also reduces the total number of approximate shortest vectors we need to compute from $\widetilde{O}(n^{2})$ to $\widetilde{O}(n)$, which opens up the gate of using faster approximate shortest vector algorithms in [51].

Although the efficiency issue of the cutting plane step and width measurement is addressed, we still need to carefully control the complexity of the dimension reduction step. This step can be particularly expensive due to the loss of structural information of the polytope $K$ after collapsing it onto a proper subspace $P$ on which the algorithm recurs, as discussed in [32]. In the standard setting of Vaidya’s cutting plane method, the polytope evolves in a “slow-changing” manner that makes it easy to maintain and update the volumetric center and its corresponding Hessian matrix. However, after reducing a dimension, we no longer have such a “slow-changing” property, and thus the new volumetric center and Hessian matrix have to be computed from scratch.

A natural idea would be to formulate the problem of recomputing the volumetric center and Hessian as a convex optimization task constrained to the polytope $K\cap P$ we have access to. On the surface, this problem can be solved straightforwardly using the same cutting plane procedure efficiently. The caveat is that evaluating the volumetric function or its gradient requires $O(n^{\omega})$ time⁴⁴4$\omega$ is the exponent of matrix multiplication [18, 69, 40]. Currently, $\omega\approx 2.373$.. As the cutting plane method evolves for $O(n\log n)$ iterations, this will lead to a total of $O(n^{\omega+1}\log n)$ arithmetic operations. To circumvent this issue, we start from a simpler convex body containing $K\cap P$ whose volumetric center can be easily computed. Specifically, we choose the hyperrectangle centered at the center of the outer ellipsoid that contains $K\cap P$. We show that this method only affects the number of oracle calls by a factor of $O(\log n)$, despite causing the volume to blow up by a factor of $n^{O(n)}$. A similar approach is also taken in [32] to ease the computation of centroid and covariance matrix of $K\cap P$. However, his method requires iterativelty refining the hyperrectangle using constraints of $K\cap P$ until it coincides with $K\cap P$, at which point the algorithm relearns the collapsed polytope $K\cap P$. This approach can take as many as $\widetilde{O}(n^{2})$ steps (without calling the separation oracle) with $\Omega(n^{2})$ operations per step, and is mainly for achieving the best possible oracle complexity. Our approach is arguably simpler and much more efficient.

Another challenge is to prove that using volumetric centers and Dikin ellipsoids [15] is sufficient to progress the algorithm. In Jiang’s algorithm, the use of centroid and covariance matrix makes the analysis straightforward due to the standard nice geometric property that the polytope $K$ is sandwiched between ellipsoids induced by the covariance matrix at the centroid. However, in our analysis, we use the volumetric center and Dikin ellipsoid of the log barrier, which requires a different approach. Past works that exclusively analyze the performance of the cutting plane method based on this approach [67] take a functional value point of view, measuring the progress of the algorithm via the change of the volumetric function. Instead, our analysis extracts key geometric information of Vaidya’s cutting plane method, showing that the volumetric center and Dikin ellipsoid evolve in a similar spirit as the centroid and covariance matrix.

In particular, we prove structural result for using Dikin ellipsoids [15] to sandwich of an $n$-dimensional polytope defined by $m$ constraints. Let us parametrize $K=\{x\in\mathbb{R}^{n}:Ax\geq b\}$ where $A\in\mathbb{R}^{m\times n}$ and $b\in\mathbb{R}^{m}$, define the Hessian of the log-barrier as $H(x)=A^{\top}S_{x}^{-2}A$ where $S_{x}\in\mathbb{R}^{m\times m}$ is the diagonal matrix with $i$-th diagonal being $a_{i}^{\top}x-b_{i}$. The volumetric function is defined as $F(x)=\frac{1}{2}\log\det(H(x))$ and we let $\operatorname{vc}(K)$ denote the minimizer of $F$, i.e. $\operatorname{vc}(K)$ is the volumetric center of $K$. The key structural result we prove is that $E(\operatorname{vc}(K),H(\operatorname{vc}(K)))\subseteq K\subseteq\operatorname{poly}(mn)\cdot E(\operatorname{vc}(K),H(\operatorname{vc}(K)))$, where $E(\operatorname{vc}(K),H(\operatorname{vc}(K)))$ is the ellipsoid centered at $\operatorname{vc}(K)$ and defined by $H(\operatorname{vc}(K))$. While it is standard for $K$ to be sandwiched by Dikin ellipsoids induced by a self-concordant barrier function and centered at the minimizer of that function, we note that $\operatorname{vc}(K)$ is not the minimizer of the log-barrier function (but rather the minimizer of $F$), while the Dikin ellipsoid is defined w.r.t. the Hessian $H$ of the log-barrier. In fact, for the Dikin ellipsoid centered at the analytic center denoted by ${\rm ac}(K)$ (minimizer of the log-barrier function) and the induced Hessian at ${\rm ac}(K)$, we have

On the other hand, to progress Vaidya’s cutting plane method, we have to work with the Dikin ellipsoid centered and induced by $\operatorname{vc}(K)$ instead of $\operatorname{ac}(K)$. We first note that the scaled volumetric function indeed defines a self-concordant barrier function and we thus have

For the left containment, we can safely replace $\nabla^{2}F(\operatorname{vc}(K))$ as for any $x\in K$, it is true that $E(x,H(x))\subseteq K$. For the right containment, we utilize the fact that $H(x)$ and $\nabla^{2}F(x)$ are close spectral approximation of each other, thus we have $E(\operatorname{vc}(K),\nabla^{2}F(\operatorname{vc}(K)))\subseteq O(\sqrt{m})\cdot E(\operatorname{vc}(K),H(\operatorname{vc}(K)))$, and we conclude

We believe this sandwiching condition might be of independent interest for other geometric applications of Vaidya’s method.

To carry out the analysis, we use the potential function developed in [32], which captures the volume of the polytope and the density of the lattice simultaneously. We show that the block of cutting plane steps also leads to rapid decrement of the potential. Unlike [32], the block cutting plane steps and Dikin ellipsoids cause extra losses on the volume of $K$ after dimension reduction, due to the possible appearance of very short vectors. We observe that such blowups are always at most $n^{O(n)}$ for each dimension reduction step. Thus, the potential increases by at most $O(n^{2}\log n)$ in total, and we need to use a total of $O(n^{2}\log n)$ oracle calls to counter this increment. To summarize, by trading for a slightly-worse number of oracle calls, we make more room for the algorithm to gain extra speedup through leverage score maintenance, faster approximate shortest vector, and crude estimation of the convex body after reducing dimension.

A natural question that arises from our result is whether it is possible to obtain a strongly-polynomial time algorithm with quadratic or subquadratic oracle complexity, matching the one in [32], while achieving the same improved runtime as ours. Our algorithm has only one $\log n$ factor in the oracle complexity, but we contend that this is inherent for Vaidya’s approach since it can only guarantee a volume decrease after $O(n\log n)$ cutting plane steps. Additionally, in the dimension reduction phase, we use a hyperrectangle to approximate the convex body, which introduces a volume increase of $n^{O(n)}$. To resolve this issue, an algorithm that approximately computes the volumetric center in $O(n^{3}\log n)$ time would be necessary. However, designing such an algorithm is an interesting and nontrivial data structure task, as discussed in the preceding subsection, since it requires maintaining and updating the log determinant of a Gram matrix under slow diagonal changes.

To achieve a subquadratic oracle complexity, a crucial ingredient in [32] is an approximate shortest vector algorithm with a sub-exponential approximation factor, first given in [2]. The main reason for the sub-exponential approximation ratio is a block-reduction scheme introduced in [56] that computes a more general notion of reduced lattice basis. Recent improvements on basis reduction algorithms make use of this block-reduction idea to obtain more refined recursive structures. Therefore, it is of interest to design approximate shortest vector algorithms using $\widetilde{O}(n^{3})$ arithmetic operations while achieving a sub-exponential approximation factor.

Given a lattice $\Lambda$ and a corresponding basis $B\in\mathbb{R}^{n\times k}$, it is natural to seek an algorithm that finds the vector with norm $\lambda_{1}(L)$, or at least approximately finds it. The famous Lenstra-Lenstra-Lovász algorithm serves such a purpose:

Lately, [51] improves the runtime of the LLL algorithm by leveraging the block reduction technique introduced in [56]. This is a key component in our $O(n^{4}\log n)$ time algorithm.

In this section, we define notions about centroid and covariance. The work [32] heavily exploits the structure of objects to design a subqudratic oracle complexity algorithm for minimizing convex functions. Our approach also relates to these notions.

For a convex body $K$, we use $\mathrm{vol}(K)$ to denote its volume, i.e., $\mathrm{vol}(K):=\int_{x\in K}\mathrm{d}x$. We first define the centroid of convex body,

Then, we define the covariance of convex body,

It is well-known that any isotropic⁵⁵5A convex body is isotropic if the uniform distribution over the body has zero mean and covariance matrix being the identity. convex body is enclosed by two balls.

We collect some standard facts of lattice projection that is directly implied by Gram-Schmidt process. We use $\Pi_{W}(\cdot)$ to denote the orthogonal projection onto the subspace $W$.

We present a variant of Lemma 3.2 in [31] in our Lemma 4.17. In particular, we replace the norm with respect to $\mathrm{Cov}(K)$ to the norm with respect to $H_{K}^{-1}$. Before stating the our new lemma (Lemma 4.17), we first recall the high dimension slicing lemma of [31].

The key building block of both [31] and our algorithm is a lattice-dimension reduction step. The following result from [31] shows that all integral points are preserved after a dimension reduction step.

We start with defining the log-barrier, and it has been widely used in convex optimization [54, 50, 33, 26, 35, 45, 63, 25].

Let $K$ be the bounded full-dimensional polytope $L=\{x:Ax\geq b\}$ where $A\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$ and $x\in\mathbb{R}^{n}$.

$H(x)$ is the Hessian of the logarithmic barrier function $\sum_{i=1}^{m}-\ln(a_{i}^{\top}x-b_{i})$ and is positive definite for all $x$ in the interior of $K$.

Observe that $F$ is a convex function, hence one can run Newton-type algorithm to approximate $\operatorname{vc}(K)$ very fast.

We start with some definitions.

Using $\sigma_{i}$, we can write the gradient of $F$ in an convenient form:

It is well-known that $Q(x)$ is a good approximation of $\nabla^{2}F(x)$:

Finally, we define $\mu(x)$ which quantifies:

The following lemma provides an upper bound on $\mu(x)$:

The following result (Lemma 4.9) bounds the progress of adding or deleting a plane of Vaidya’s CPM.

Next, in Lemma 4.10, we show that after $T=O(n\log m)$ iterations, the volume of the resulting convex body is only $(\frac{1}{m})^{n}$ fraction of the original convex body.