We consider the MAP-MRF inference task, that is, minimizing a function of discrete variables represented as a sum of unary and pairwise terms. A prominent approach for tackling this NP-hard problem in practice is to solve its natural LP relaxation and then iteratively tighten the relaxation by adding clusters. Based on some theoretical observations, we propose a new technique for identifying such clusters. It works by running the Singleton Arc Consistency algorithm in a certain CSP instance. Experimental results indicate that the new tightening technique outperforms the previous approach by [Sontag et al. UAI 2012] that searches for frustrated cycles. Our code will be made available at https://github.com/vnk-ist/MAP-MRF/.

We present a family of fast pseudo-approximation algorithms for the minimum balanced vertex separator problem in a graph. Given a graph G=(V,E) with n vertices and m edges, and a (constant) balance parameter c∈(0,1/2), where G has some (unknown) c-balanced vertex separator of size OPT_c, we give a (Monte-Carlo randomized) algorithm running in O(n^O(ε)m^1+o(1)) time that produces a Θ(1)-balanced vertex separator of size O( OPT_c·√((log n)/ε)) for any value ε∈[Θ(1/log(n)),Θ(1)]. In particular, for any function f(n)=ω(1) (including f(n)=loglog n, for instance), we can produce a vertex separator of size O( OPT_c·√(log n)· f(n)) in time O(m^1+o(1)). Moreover, for an arbitrarily small constant ε=Θ(1), our algorithm also achieves the best-known approximation ratio for this problem in O(m^1+Θ(ε)) time. The algorithms are based on a semidefinite programming (SDP) relaxation of the problem, which we solve using the Matrix Multiplicative Weight Update (MMWU) framework of Arora and Kale. Our oracle for MMWU uses O(n^O(ε)polylog(n)) almost-linear time maximum-flow computations, and would be sped up if the time complexity of maximum-flow improves.

We consider several problems related to packing forests in graphs. The first one is to find k edge-disjoint forests in a directed graph G of maximal size such that the indegree of each vertex in these forests is at most k. We describe a min-max characterization for this problem and show that it can be solved in almost linear time for fixed k, extending the algorithm of [Gabow, 1995]. Specifically, the complexity is O(k δm log n), where n, m are the number of vertices and edges in G respectively, and δ= max{1, k - k_G}, where k_G is the edge connectivity of the graph. Using our solution to this problem, we improve complexities for two existing applications: (1) k-forest problem: find k forests in an undirected graph G maximizing the number of edges in their union. We show how to solve this problem in O(k^3 min{kn, m}log^2 n + k · MAXFLOW(m, m) log n) time, breaking the O_k(n^3/2) complexity barrier of previously known approaches. (2) Directed edge-connectivity augmentation problem: find a smallest set of directed edges whose addition to the given directed graph makes it strongly k-connected. We improve the deterministic complexity for this problem from O(k δ(m+δn)log n) [Gabow, STOC 1994] to O(k δm log n). A similar approach with the same complexity also works for the undirected version of the problem.

In Numerical Algebraic Geometry (NAG) isolated solutions of polynomial systems are usually computed by tracking a solution curve defined by a homotopy equation. The tracking problem becomes especially challenging close to a singular root (the “endgame” regime). Existing approaches include power series endgames, Cauchy endgames, and various methods that regularize the system via dual-space-based deflation. We make the following contributions. (1) For corank-1 systems we introduce a new “Arclength Endgame” which combines the idea of the classical pseudo-arclength continuation method with the estimation of the Puiseux series of the curve. We formally prove that it has a superlinear rate of convergence in some neighborhood of the root. The method uses only evaluations of the system and its Jacobian, whereas previous techniques with proven superlinear convergence (such as deflation) require computing additional derivatives of the system. (2) For systems with a larger corank we propose a heuristic “Lifted Arclength Endgame”, which shows promising experimental results. (3) A key step in our approach (as well as in the standard power series endgame) is estimating the Puiseux series of the curve, which is characterized by fractional exponents k_i/c for i≥ 1 together with associated coefficients. Previous work addressed only estimating the ratio k_1/c. We present a new method for that which empirically appears to be more stable than previous methods, and also show how to estimate k_i/c for i≥ 2.

The computational equivalence between approximate counting and sampling is well established for polynomial-time algorithms. The most efficient general reduction from counting to sampling is achieved via simulated annealing, where the counting problem is formulated in terms of estimating the ratio Q=Z(β_max)/Z(β_min) between partition functions Z(β)=∑_x∈ Ωexp(βH(x)) of Gibbs distributions μ_β over Ω with Hamiltonian H, given access to a sampling oracle that produces samples from μ_β for β∈ [β_min, β_max]. The best bound achieved by known annealing algorithms with relative error ε is O(q log h / ε^2), where q, h are parameters which respectively bound ln Q and H. However, all known algorithms attaining this near-optimal complexity are inherently sequential, or *adaptive*: the queried parameters β depend on previous samples. We develop a simple non-adaptive algorithm for approximate counting using O(q log^2 h / ε^2) samples, as well as an algorithm that achieves O(q log h / ε^2) samples with just two rounds of adaptivity, matching the best sample complexity of sequential algorithms. These algorithms naturally give rise to work-efficient parallel (RNC) counting algorithms. We discuss applications to RNC counting algorithms for several classic models, including the anti-ferromagnetic 2-spin, monomer-dimer and ferromagnetic Ising models.

We implement an algorithm for solving the minimum weight perfect matching problem. Our code significantly outperforms the current state-of-the-art Blossom V algorithm on those families of instances where Blossom V takes superlinear time. In practice, our implementation shows almost-linear runtime on every family of instances on which we have tested it. Our algorithm relies on solving the maximum-cardinality unweighted matching problems during its primal phase. Following the state-of-the-art cherry blossom algorithm, we use cherry trees instead of traditional alternating trees and cherry blossoms instead of traditional blossoms. We shrink cherry blossoms rather than traditional blossoms into supernodes. This strategy allows us to deal with much shallower supernodes.

The Lovasz Local Lemma (LLL) is a powerful tool in probabilistic combinatorics which can be used to establish the existence of objects that satisfy certain properties. The breakthrough paper of Moser and Tardos and follow-up works revealed that the LLL has intimate connections with a class of stochastic local search algorithms for finding such desirable objects. In particular, it can be seen as a sufficient condition for this type of algorithms to converge fast. Besides conditions for existence of and fast convergence to desirable objects, one may naturally ask further questions regarding properties of these algorithms. For instance, are they parallelizable?, how many solutions can they output?, what is the expected weight of a solution?, etc. These questions and more have been answered for a class of LLL-inspired algorithms called commutative. In this paper we introduce a new, very natural and more general notion of commutativity (essentially matrix commutativity) which allows us to show a number of new refined properties of LLL-inspired local search algorithms with significantly simpler proofs.

Seymour's decomposition theorem is a hallmark result in matroid theory presenting a structural characterization of the class of regular matroids. Formalization of matroid theory faces many challenges, most importantly that only a limited number of notions and results have been implemented so far. In this work, we formalize the proof of the forward (composition) direction of Seymour's theorem for regular matroids. To this end, we develop a library in Lean 4 that implements definitions and results about totally unimodular matrices, vector matroids, their standard representations, regular matroids, and 1-, 2-, and 3-sums of matrices and binary matroids given by their standard representations. Using this framework, we formally state Seymour's decomposition theorem and implement a formally verified proof of the composition direction in the setting where the matroids have finite rank and may have infinite ground sets.

Currently, the best known tradeoff between approximation ratio and complexity for the Sparsest Cut problem is achieved by the algorithm in Sherman [FOCS, 2009]: It computes \(O(\sqrt{(\log n)/\varepsilon})\) -approximation using \(O(n^{\varepsilon}\log^{O(1)}n)\) maxflows for any \(\varepsilon\in[\Theta(1/\log n),\Theta(1)]\) . It works by solving the SDP relaxation of Arora et al. [STOC, 2004] using the Multiplicative Weights (MW) Update algorithm of Arora and Kale [JACM, 2016]. To implement one MW step, Sherman approximately solves a multicommodity flow problem using another application of MW. Nested MW steps are solved via a certain “chaining” algorithm that combines results of multiple calls to the maxflow algorithm. We present an alternative approach that avoids solving the multicommodity flow problem and instead computes “violating paths.” This simplifies Sherman’s algorithm by removing a need for a nested application of MW and also allows parallelization: We show how to compute \(O(\sqrt{(\log n)/\varepsilon})\) -approximation via \(O(\log^{O(1)}n)\) maxflows using \(O(n^{\varepsilon})\) processors. We also revisit Sherman’s chaining algorithm and present a simpler version together with a new analysis.

Currently, the best known tradeoff between approximation ratio and complexity for the Sparsest Cut problem is achieved by the algorithm in [Sherman, FOCS 2009]: it computes O(√((log n)/ε))-approximation using O(n^εlog^O(1)n) maxflows for any ε∈[Θ(1/log n),Θ(1)]. It works by solving the SDP relaxation of [Arora-Rao-Vazirani, STOC 2004] using the Multiplicative Weights Update algorithm (MW) of [Arora-Kale, JACM 2016]. To implement one MW step, Sherman approximately solves a multicommodity flow problem using another application of MW. Nested MW steps are solved via a certain “chaining” algorithm that combines results of multiple calls to the maxflow algorithm. We present an alternative approach that avoids solving the multicommodity flow problem and instead computes “violating paths”. This simplifies Sherman's algorithm by removing a need for a nested application of MW, and also allows parallelization: we show how to compute O(√((log n)/ε))-approximation via O(log^O(1)n) maxflows using O(n^ε) processors. We also revisit Sherman's chaining algorithm, and present a simpler version together with a new analysis.

We consider the problem of estimating the partition function Z(β)=∑_x exp(β(H(x)) of a Gibbs distribution with the Hamiltonian H:Ω→{0}∪[1,n]. As shown in [Harris Kolmogorov 2024], the log-ratio q=ln (Z(β_max)/Z(β_min)) can be estimated with accuracy ε using O(q log n/ε^2) calls to an oracle that produces a sample from the Gibbs distribution for parameter β∈[β_min,β_max]. That algorithm is inherently sequential, or adaptive: the queried values of β depend on previous samples. Recently, [Liu, Yin Zhang 2024] developed a non-adaptive version that needs O( q (log^2 n) (log q + loglog n + ε^-2) ) samples. We improve the number of samples to O(q log^2 n/ε^2) for a non-adaptive algorithm, and to O(q log n/ε^2) for an algorithm that uses just two rounds of adaptivity (matching the complexity of the sequential version). Furthermore, our algorithm simplifies previous techniques. In particular, we use just a single estimator, whereas methods in [Harris Kolmogorov 2024, Liu, Yin Zhang 2024] employ two different estimators for different regimes.

This paper deals with the algorithmic aspects of solving feasibility problems of semidefinite programming (SDP), aka linear matrix inequalities (LMI). Since in some SDP instances all feasible solutions have irrational entries, numerical solvers that work with rational numbers can only find an approximate solution. We study the following question: is it possible to certify feasibility of a given SDP using an approximate solution that is sufficiently close to some exact solution? Existing approaches make the assumption that there exist rational feasible solutions (and use techniques such as rounding and lattice reduction algorithms). We propose an alternative approach that does not need this assumption. More specifically, we show how to construct a system of polynomial equations whose set of real solutions is guaranteed to have an isolated correct solution (assuming that the target exact solution is maximum-rank). This allows, in particular, to use algorithms from real algebraic geometry for solving systems of polynomial equations, yielding a hybrid (or symbolic-numerical) method for SDPs. We experimentally compare it with a pure symbolic method; the hybrid method was able to certify feasibility of many SDP instances on which the exact method failed. Our approach may have further applications, such as refining an approximate solution using methods of numerical algebraic geometry for systems of polynomial equations.

Greedy minimum weight spanning tree packings have proven to be useful in connectivity-related problems. We study the process of greedy minimum weight base packings in general matroids and explore its algorithmic applications. When specialized to bicircular matroids, our results yield an algorithm for the approximate fully-dynamic densest subgraph density ρ. We maintain a (1+ε)-approximation of the density with a worst-case update time O((ρε^-2+ε^-4)ρlog^3 m). It improves the dependency on ε from the current state-of-the-art worst-case update time complexity O(ε^-6log^3 n) [Chekuri, Christiansen, Holm, van der Hoog, Quanrud, Rotenberg, Schwiegelshohn, SODA'24]. We also can maintain an implicit fractional out-orientation with a guarantee that all out-degrees are at most (1+ε)ρ. Our algorithms above work by greedily packing pseudoforests, and require maintenance of a minimum-weight pseudoforest in a dynamically changing graph. We show that this problem can be solved in O(log n) worst-case time per edge insertion or deletion. For general matroids, we observe two characterizations of the limit of the base packings (“the vector of ideal loads”), which imply the characterizations from [Cen, Fleischmann, Li, Li, Panigrahi, FOCS'25], namely, their entropy-minimization theorem and their bottom-up cut hierarchy. Finally, we give combinatorial results on the greedy tree packings. We show that a tree packing of O(λ^5log m) trees contains a tree crossing some min-cut once, which improves the bound O(λ^7log^3 m) from [Thorup, Combinatorica'07]. We also strengthen the lower bound on the edge load convergence rate from [de Vos, Christiansen, SODA'25], showing that Thorup's upper bound is tight up to a logarithmic factor.

We consider two problems for a directed graph G, which we show to be closely related. The first one is to find k edge-disjoint forests in G of maximal size such that the indegree of each vertex in these forests is at most k. We describe a min-max characterization for this problem and show that it can be solved in O(k δm log n) time, where (n,m) is the size of G and δ is the difference between k and the edge connectivity of the graph. The second problem is the directed edge-connectivity augmentation problem, which has been extensively studied before: find a smallest set of directed edges whose addition to the graph makes it strongly k-connected. We improve the complexity for this problem from O(k δ(m+δn)log n) [Gabow, STOC 1994] to O(k δm log n), by exploiting our solution for the first problem. A similar approach with the same complexity also works for the undirected version of the problem.

Given a fixed finite metric space (V,μ ) , the minimum 0-extension problem, denoted as 0-𝙴𝚡𝚝[μ] , is equivalent to the following optimization problem: minimize function of the form min _x∈ V^n∑ _i f_i(x_i) + ∑ _ij c_ijμ (x_i,x_j) where f_i:V→ℝ are functions given by f_i(x_i)=∑ _v∈ V c_viμ (x_i,v) and c_ij,c_vi are given nonnegative costs. The computational complexity of 0-𝙴𝚡𝚝[μ] has been recently established by Karzanov and by Hirai: if metric μ is orientable modular then 0-𝙴𝚡𝚝[μ] can be solved in polynomial time, otherwise 0-𝙴𝚡𝚝[μ] is NP-hard. To prove the tractability part, Hirai developed a theory of discrete convex functions on orientable modular graphs generalizing several known classes of functions in discrete convex analysis, such as L^♮ -convex functions. We consider a more general version of the problem in which unary functions f_i(x_i) can additionally have terms of the form c_uv;iμ (x_i,{u,v}) for {u,v}∈ F , where set F⊆( [ V; 2 ]) is fixed. We extend the complexity classification above by providing an explicit condition on (μ ,F) for the problem to be tractable. In order to prove the tractability part, we generalize Hirai’s theory and define a larger class of discrete convex functions. It covers, in particular, another well-known class of functions, namely submodular functions on an integer lattice. Finally, we improve the complexity of Hirai’s algorithm for solving 0-𝙴𝚡𝚝[μ] on orientable modular graphs.

Farkas established that a system of linear inequalities has a solution if and only if we cannot obtain a contradiction by taking a linear combination of the inequalities. We state and formally prove several Farkas-like theorems over linearly ordered fields in Lean 4. Furthermore, we extend duality theory to the case when some coefficients are allowed to take “infinite values”.

We consider the problem of solving LP relaxations of MAP-MRF inference problems, and in particular the method proposed recently in [16], [35]. As a key computational subroutine, it uses a variant of the Frank-Wolfe (FW) method to minimize a smooth convex function over a combinatorial polytope. We propose an efficient implementation of this subroutine based on in-face Frank-Wolfe directions, introduced in [4] in a different context. More generally, we define an abstract data structure for a combinatorial subproblem that enables in-face FW directions, and describe its specialization for tree-structured MAP-MRF inference subproblems. Experimental results indicate that the resulting method is the current state-of-art LP solver for some classes of problems. Our code is available at pub.ist.ac.at/~vnk/papers/IN-FACE-FW.html.

A central problem in computational statistics is to convert a procedure for sampling combinatorial objects into a procedure for counting those objects, and vice versa. We consider sampling problems coming from Gibbs distributions , which are families of probability distributions over a discrete space \(\Omega\) with probability mass function of the form \(\mu^{\Omega}_{\beta}(\omega)\propto e^{\beta H(\omega)}\) for \(\beta\) in an interval \([\beta_{\min},\beta_{\max}]\) and \(H(\omega)\in\{0\}\cup[1,n]\) . Two important parameters are the partition function , which is the normalization factor \(Z(\beta)=\sum_{\omega\in\Omega}e^{\beta H(\omega)}\) , and the vector of preimage counts \(c_{x}=|H^{-1}(x)|\) . We develop black-box sampling algorithms to estimate the counts using roughly \(\tilde{O}(\frac{n^{2}}{\varepsilon^{2}})\) samples for integer-valued distributions and \(\tilde{O}(\frac{q}{\varepsilon^{2}})\) samples for general distributions, where \(q=\log\frac{Z(\beta_{\max})}{Z(\beta_{\min})}\) (ignoring some second-order terms and parameters). We show this is optimal up to logarithmic factors. We illustrate with improved algorithms for counting connected subgraphs, independent sets, and perfect matchings. As a key subroutine, we estimate all values of the partition function using \(\tilde{O}(\frac{n^{2}}{\varepsilon^{2}})\) samples for integer-valued distributions and \(\tilde{O}(\frac{q}{\varepsilon^{2}})\) samples for general distributions. This improves over a prior algorithm of Huber (2015) which computes a single point estimate \(Z(\beta_{\max})\) and which uses a slightly larger amount of samples. We show matching lower bounds, demonstrating this complexity is optimal as a function of \(n\) and \(q\) up to logarithmic terms.

This paper studies algorithms for computing a Gomory-Hu tree, which is a classical data structure that compactly stores all minimum $s$-$t$ cuts of an undirected weighted graph. We consider two classes of algorithms: the original method by Gomory and Hu and the method based on "OrderedCuts" that we recently proposed. We describe practical implementations of these methods, and compare them experimentally with the algorithms from the previous experimental studies by Goldberg and Tsioutsiouliklis (2001) and by Akibo et al. (2016) (designed for unweighted simple graphs). Results indicate that the method based on OrderedCuts is the most robust, and often outperforms other implementations by a large factor.

We consider two models for the sequence labeling (tagging) problem. The first one is a Pattern-Based Conditional Random Field (PB), in which the energy of a string (chain labeling) x=x1⁢…⁢xn∈Dn is a sum of terms over intervals [i,j] where each term is non-zero only if the substring xi⁢…⁢xj equals a prespecified word w∈Λ. The second model is a Weighted Context-Free Grammar (WCFG) frequently used for natural language processing. PB and WCFG encode local and non-local interactions respectively, and thus can be viewed as complementary. We propose a Grammatical Pattern-Based CRF model (GPB) that combines the two in a natural way. We argue that it has certain advantages over existing approaches such as the Hybrid model of Benedí and Sanchez that combines N -grams and WCFGs. The focus of this paper is to analyze the complexity of inference tasks in a GPB such as computing MAP. We present a polynomial-time algorithm for general GPBs and a faster version for a special case that we call Interaction Grammars .