  In this work, we construct a novel numerical method for solving the
multi-marginal optimal transport problems with Coulomb cost. This type of
optimal transport problems arises in quantum physics and plays an important
role in understanding the strongly correlated quantum systems. With a
Monge-like ansatz, the orginal high-dimensional problems are transferred into
mathematical programmings with generalized complementarity constraints, and
thus the curse of dimensionality is surmounted. However, the latter ones are
themselves hard to deal with from both theoretical and practical perspective.
Moreover in the presence of nonconvexity, brute-force searching for global
solutions becomes prohibitive as the problem size grows large. To this end, we
propose a global optimization approach for solving the nonconvex optimization
problems, by exploiting an efficient proximal block coordinate descent local
solver and an initialization subroutine based on hierarchical grid refinements.
We provide numerical simulations on some typical physical systems to show the
efficiency of our approach. The results match well with both theoretical
predictions and physical intuitions, and give the first visualization of
optimal transport maps for some two dimensional systems.
