10.2.13 Earth Mover’s Distance Analysis

Let $𝐫$ and $𝐜$ denote two discrete nonnegative distributions defined on a real-space grid, with equal total weight. In the EMD framework, these are interpreted as the source and target distributions, respectively. For a distance matrix $𝐌$, whose elements $M_{ij}$ represent the distance between grid points $i$ and $j$, the Earth Mover’s Distance is defined as

where $𝐏$ is the transport plan and

The admissible transport plans are therefore constrained to reproduce the prescribed source and target marginals. In other words, each row sum of $𝐏$ must equal the source distribution, and each column sum must equal the target distribution.

In the present application, the source and target distributions are constructed from the positive and negative parts of the real-space density difference associated with an excitation. Let

where ${\rho }^{\mathrm{GS}}(𝐫)$ and ${\rho }^{\mathrm{ES}}(𝐫)$ are the ground-state and excited-state densities, respectively. The positive and negative components of the density difference define the two charge distributions to be compared,

After discretization on a grid with points $\{{𝐫}_i\}$ and quadrature weights $\{w_i\}$, the effective charge associated with each grid point is

and the two nonnegative discrete distributions entering the EMD problem are obtained from the positive and negative components of the density difference. The total transferred charge is then

and the corresponding transport cost,

has units of charge times distance. A corresponding length-like quantity may be defined as

The exact EMD may be obtained by solving this constrained optimization problem using a transportation simplex algorithm. While this approach yields the exact optimal transport cost, its computational cost grows rapidly with grid size and can become prohibitive for large systems or dense real-space grids.

A direct EMD calculation on the full quadrature grid used in a DFT or TDDFT calculation is generally impractical because the number of grid points may be very large. Since the optimal transport problem scales unfavorably with the size of the discrete distributions, the present implementation employs a reduced grid for the EMD evaluation.

In this procedure, the electronic density is first evaluated on the underlying quadrature grid used in the electronic-structure calculation. A smaller atom-centered grid is then introduced for the EMD analysis. Each grid point from the original quadrature grid is assigned to a point on the smaller EMD grid according to spatial proximity, and the effective charges from the original grid are accumulated onto the selected EMD grid. In this way, the fine-grid density is projected onto a smaller discrete representation that preserves the overall distribution while substantially reducing the number of points entering the transport problem.

Let $\{{𝐫}_k^{\mathrm{big}},w_k^{\mathrm{big}},{\rho }_k^{\mathrm{big}}\}$ denote the grid points, weights, and densities on the original quadrature grid, and let $\{{𝐫}_i^{\mathrm{small}}\}$ denote the reduced grid used for the EMD evaluation. The effective charge assigned to point $i$ on the reduced grid is constructed by accumulating the charges from all original grid points associated with that reduced point,

where ${\mathrm{\Omega }}_i$ denotes the set of original grid points assigned to reduced grid point $i$. The resulting reduced-grid charges are then used to define the source and target distributions entering the EMD calculation.

This reduced-grid procedure greatly lowers the cost of both exact and regularized EMD calculations while retaining a physically meaningful representation of the density difference. The reduced grid is chosen from standard atom-centered quadrature grids in Q-Chem. The quality of the EMD result therefore depends on the quality of the reduced-grid representation of the density difference, and convergence with respect to the EMD grid should be checked whenever high quantitative accuracy is required.

To further reduce the cost of the EMD evaluation, one may replace the exact optimal transport problem by an entropically regularized one,

where $h(𝐏)$ is the entropy of the transport plan and $\lambda$ is the regularization parameter. The corresponding regularized transport cost is

The entropy term favors smoother transport plans and removes the need to solve the original linear programming problem directly. The parameter $\lambda$ controls the balance between smoothness and fidelity to the exact EMD. Larger values of $\lambda$ reduce the influence of the entropy term and drive the solution closer to the exact EMD, whereas smaller values of $\lambda$ produce a smoother transport plan and generally improve numerical stability.

The regularized problem admits a factorized form for the optimal transport plan,

where

The vectors $𝐮$ and $𝐯$ are positive scaling factors chosen such that the transport plan satisfies the marginal constraints,

Substituting the factorized form into these constraints yields the coupled fixed-point equations

These equations form the basis of the Sinkhorn–Knopp algorithm, which alternates updates of $𝐮$ and $𝐯$ until the marginal constraints are satisfied to within a specified threshold.

The computational steps of the method are therefore as follows. First, the density difference is evaluated on the underlying quadrature grid. Second, the associated charges are projected onto a smaller EMD grid. Third, the distance matrix $𝐌$ is constructed on that reduced grid, and the kernel matrix $𝐊$ is formed according to Eq. (10.49). Finally, the Sinkhorn iterations are performed to obtain the regularized transport plan and the corresponding transport cost. The dominant cost of the regularized method is associated with repeated matrix-vector products involving $𝐊$ and ${𝐊}^T$, which is substantially less expensive than solving the exact transportation simplex problem for large grids.

In practice, the principal numerical parameter is the regularization strength $\lambda$. This parameter must be chosen large enough that the regularized transport cost remains close to the exact EMD, but not so large that the kernel matrix becomes severely ill-conditioned. If $\lambda M_{ij}$ is too large, the exponential kernel may underflow, leading to unstable scaling vectors and poor convergence behavior. The present implementation therefore exposes user-adjustable regularization and convergence parameters, where The user-specified regularization parameter is scaled by the maximum element of the distance matrix. In other words, the effective regularization strength used in the calculation is ${\lambda }_{\mathrm{eff}}=\frac{{\lambda }_{\mathrm{in}}}{M_{ij}^{\max }}$, where $M_{ij}^{\max }={\max }_{ij}{M}_{ij}$. In general, increasing ${\lambda }_{\mathrm{eff}}$ improves agreement with the exact EMD, whereas decreasing ${\lambda }_{\mathrm{eff}}$ improves numerical stability and convergence behavior.

In summary, this feature evaluates the Earth Mover’s Distance between two discrete charge distributions on a real-space grid either using exact transport algorithm or using an entropically regularized optimal transport formulation solved by the Sinkhorn–Knopp algorithm, with a reduced atom-centered grid for the EMD calculation rather than the full quadrature grid employed in the underlying electronic-structure calculation. The quality of the result therefore depends both on the EMD grid and on the regularization parameters. The principal user-controlled options are the activation of EMD analysis itself EMD_ANALYSIS, the grid used for constructing the reduced real-space charge distributions EMD_FIT_GRID, the use of the Sinkhorn solver SINKHORN, the regularization strength for Sinkhorn solver SINKHORN_REG, the convergence threshold for Sinkhorn solver SINKHORN_THRESHOLD, and the maximum cycles MAX_SINKHORN_CYCLES. Larger regularization parameters generally improve agreement with the exact EMD but may lead to numerical instability if the kernel matrix becomes too sharply peaked. Smaller regularization parameters improve numerical robustness, but introduce a larger approximation error relative to the exact EMD.