Genetic column generation: Fast computation of high-dimensional multi-marginal optimal transport problems

Abstract

We introduce a simple, accurate, and extremely efficient method for numerically solving the multi-marginal optimal transport (MMOT) problems arising in density functional theory. The method relies on (i) the sparsity of optimal plans [for N marginals discretized by gridpoints each, general Kantorovich plans require N gridpoints but the support of optimizers is of size O( • N ) [FV18]], (ii) the method of column generation (CG) from discrete optimization which to our knowledge has not hitherto been used in MMOT, and (iii) ideas from machine learning. The well-known bottleneck in CG consists in generating new candidate columns efficiently; we prove that in our context, finding the best new column is an NP-complete problem. To overcome this bottleneck we use a genetic learning method tailormade for MMOT in which the dual state within CG plays the role of an "adversary", in loose similarity to Wasserstein GANs. On a sequence of benchmark problems with up to 120 gridpoints and up to 30 marginals, our method always found the exact optimizers. Moreover, empirically the number of computational steps needed to find them appears to scale only polynomially when both N and are simultaneously increased (while keeping their ratio fixed to mimic a thermodynamic limit of the particle system).

The strong-interaction limit of density functional theory

“…This method was proposed recently by Friesecke, Schulz, and Vögler [51]. It directly solves the discretized SIL problem (108), by combining the sparse but exact quasi-SCE or quasi-Monge ansatz (see Theorem 2.23), the method of column generation from discrete optimization, and basic ideas from machine learning.…”

“…N ) are arbitrary N -point configurations in X N and β > 1 is a hyperparameter (taken to be 5 in [51]) which limits the number of N -point configurations to O( ) instead of the naively required O( N ). To achieve a unique correspondence between symmetrized Diracs and N -point configurations one restricts the r (ν) to the sector X N sym = {(a i 1 , ..., a i N ) ∈ X N : i 1 ≤ ... ≤ i N }, making the expansion coefficients α ν in (120) unique.…”

“…Adding this configuration to the set Ω "cuts off" u new from the optimization domain of the dual problem, yielding a new dual solution and an energy decrease. For a rigorous justification see [51].…”

“…Figure 10, taken from [51], shows the solution of the SIL problem (25) computed by the GenCol algorithm for 10 electrons in a 1D interval discretized by 100 gridpoints. In this example, the grid spacing is normalized to 1, the density is taken to be ρ(x) = const(0.2 + sin 2 ( x +1 )), and the interaction is the soft Coulomb potential (104) with a = 0.1.…”

“…Tests reported in [51] on larger 1D systems with up to N = 30 electrons on 120 grid points (corresponding to a space of N -point densities of dimension N ≈ 2.4 × 10 62 ) show only a slow polynomial growth in N of the number of iterations required to find the exact solution to machine precision, with the average number of samples needed per iteration to satisfy the acceptance criterion remaining approximately constant.…”