Beyond Procrustes distances: a multilinear Gromov-Wasserstein distance capturing chirality

Motivation

• Standard Procrustes distance cannot distinguish chiral shapes, although molecular chirality can determine whether compounds are therapeutic or toxic.D-penicillamine treats rheumatoid arthritis, whereas its reflected form, L-penicillamine, is toxic.Source: §1 · p. 1

Shape representation

• Shapes are represented as probability measures in P2(X), with distances made invariant to spatial transformations by quotienting Wasserstein distance over a transformation group G.The framework considers X ⊂ R^d with d = 2 or 3 and groups including SO(d), E(d), and SL(d).Source: §1 · pp. 1–2

Computational challenge

• Procrustes-Wasserstein optimization is non-convex, locally initialization-sensitive, and lacks proved global convergence guarantees, with chirality requiring SO(d) rather than O(d).Fixed-coupling or fixed-transformation cases are tractable, but global optimization remains computationally challenging.Source: §1 · p. 2

Related distances

• Gromov-Wasserstein distances provide direct invariance to group actions but are generally NP-hard to compute, motivating multilinear extensions for transform-sensitive shape comparison.GW variants have been applied to brain images, proteins, graphs, and structured objects; Inner Gromov-Wasserstein uses inner-product similarity [2] [42].Source: §1 · p. 2

Contributions

• The paper introduces multilinear Gromov-Wasserstein discrepancies that can become proper shape distances, with stability, low-dimensional structure, and optimal correspondence guarantees.The framework targets distances sensitive to specified transforms, including SE(d), which yields the Chiral Gromov-Wasserstein distance.Source: §1 · p. 2

Algorithms

• A low-dimensional reformulation enables a fully polynomial-time approximation scheme for globally solving the multilinear Gromov-Wasserstein problem.The authors present this as the first estimation of complexity for globally solving a Gromov-Wasserstein-type problem.Source: §1 · p. 2

2 Multilinear Gromov-Wasserstein distances and their fundamental properties

• Shapes are probability distributions modulo a specified matrix Lie-group action, generalizing classical shape invariance but not providing scaling invariance automatically.The framework denotes these equivalence classes as G-shapes and specializes G as needed.Source: §2 · p. 3

2.1 A distance between shapes

• Multilinear forms define invariant discrepancies for transformation groups preserving those forms, including IGW for O(d), DGW for SL(d), and chirality-sensitive CGW for SO(d).CGW combines inner-product and determinant invariances, with a fixed parameter 0 < t < 1.Source: §2.1 · pp. 3–4
• Under non-degeneracy, GWm separates compactly supported absolutely continuous measures and positive-volume point clouds exactly up to the corresponding group action.Centering the marginals yields translation invariance, while scaling transports optimal couplings by pushing them forward on the scaled marginal.Source: §2.1 · pp. 4–5

2.2 Characterization of the multilinear Gromov-Wasserstein problem with the covariance matrix

• GWm optimization depends only on the coupling’s d^2-dimensional cross-covariance projection, reducing computation to a polynomial over the compact convex projected coupling set.This reformulation exposes the low-dimensional structure used for approximate global optimization.Source: §2.2 · p. 6, Figure 1 · p. 7

2.3 Properties of the optimal coupling

• For W2-convergent marginals, the polynomial GWm optimum is attained, optimal values converge, and minimizing couplings are precompact with limit points that remain optimal.The distance also admits quantitative marginal-stability bounds under uniformly bounded second moments.Source: §2.3 · p. 8
• If the first marginal has compact support and a density, an optimal GWm plan exists as a Monge map; general optimal plans need not be maps and converge only up to subsequences.These guarantees support shape alignment and registration, including controlled approximation when continuous marginals are discretized.Source: §2.3 · p. 9

3.1 A Sandwich algorithm to approximate polytopes

• The method approximates the projected coupling polytope PΠ by iteratively refining inner and outer bounding boxes that converge at a known rate.The inner box is formed from boundary points, whereas the outer box is an intersection of supporting half-spaces; under Assumption 1 and controlled initialization, convergence has rate 1/d^2−1.Source: §3 · p. 9, §3.1 · p. 10

3.2 Optimization of the cost

• After bounding-box refinement, polynomial optimization supplies lower and upper cost bounds that certify global optimality, while an optimal transport problem recovers an approximate coupling.For general costs, the polynomial problem can use simulated annealing, branch and bound, or low-level Lasserre semidefinite relaxations; local optimization can then refine the coupling.Source: §3 · p. 9, §3.2 · p. 12

3.4 Local optimization

• The low-dimensional Frank-Wolfe reformulation computes local optima using only the coupling cross-covariance matrix, with objective convergence at rate O(1/√k).Each iteration solves an optimal transport problem, and the line search is one-dimensional polynomial optimization that is analytically implemented for IGW and for DGW and CGW in dimensions 2 or 3.Source: §3 · p. 9, §3.4 · p. 15

3.3 Global optimization for a concave cost

• For concave costs, global optimization reduces to evaluating the cost at bounding-box vertices and solving classical optimal transport problems.This applies to IGW in all dimensions and to CGW for 1/2 ≤ t < 1 in dimension d = 2, while in d = 3 CGW concavity depends on the marginal supports and parameter.Source: §3 · pp. 9–10, §3.3 · pp. 12–13
• For discrete marginals with at most N points, the concave-cost algorithm is an FPTAS whose iteration count is independent of N, with data-size dependence arising through optimal transport.The algorithm uses O(N^3) optimal transport per iteration, or entropic regularization with complexity O(N^2 ln(N)ϵ^-3), although Sinkhorn noise leaves the convergence rate unknown.Source: §3.3 · p. 14

• CGW assigns a positive distance to mirrored hands and larger distances to reflected than non-reflected penicillamine conformations, capturing chirality.Classical, inner, and 2-Gromov-Wasserstein comparisons fail to detect the hand reflection; for penicillamine, CGW and RMSD are linearly correlated without reflection, and enumeration confirms the 9-atom optimization reached the global solution.Source: §4 · pp. 16–17, Figure 2 · p. 16
• The global optimizer scales better empirically than its theoretical bound, with linear vertices and constraints in 2D and O(k^3) growth in 3D versus O(k^4).The 2D behavior enabled several thousand iterations without memory issues on a standard laptop.Source: §4 · p. 17
• Bounding boxes with approximately 10^5 vertices or constraints create the main memory-time bottleneck for GWm, even when optimal transport over the shape points is fast.The limitation was observed in the penicillamine case study despite its nine points.Source: §4 · p. 17
• Despite a slow theoretical convergence rate, Algorithm 3 converged in fewer than 50 iterations in the numerical experiments.This suggests convergence rate was not the primary limitation in globally solving the problem.Source: §4 · p. 17
• The proposed algorithm performed similarly to Ryner and colleagues’ method [35], with differences in iteration counts and marginal size.The comparison builds on prior benchmarking showing that Ryner and colleagues’ method [35] outperformed generic nonlinear global optimization algorithms.Source: §4 · p. 17

• GWm compares shapes up to matrix Lie-group transformations when an associated non-degenerate multilinear form exists, with robust marginals and optimal correspondence maps.These properties support shape registration and alignment.Source: §5 · p. 17
• The proposed local scheme and fully polynomial-time global approximation make GWm tractable when the cost is concave.For sufficiently large t, CGW_t is concave and computable by the global algorithm.Source: §5 · pp. 17–18
• Tests in dimensions 2 and 3 confirm that CGW captures changes in shape chirality, while further optimization could enable more iterations in three dimensions.Global-optimality certificates can control per-iteration precision, and large datasets remain feasible when their optimal transport problems can be solved.Source: §5 · p. 18
• Entropy-regularized or low-rank optimal transport could scale the method, but the stability impact of noise from these approximations remains open.The implementation, tests, datasets, and numerical experiments are publicly provided in the Python repository; the penicillamine structure comes from PubChem.Source: §5 · p. 18

Notation

• Vectors x ∈ X are represented by coordinates (x_1, …, x_d), and families of k vectors are identified with elements of R^(dk).Source: Notation · p. 22
• P(X) denotes probability measures on X, while P2(X) restricts to measures with finite second moment.Source: Notation · p. 22
• G denotes a matrix Lie group over R, including O(d), SO(d), SL(d), E(d), and SE(d).These examples cover orthogonal, special orthogonal, special linear, Euclidean, and special Euclidean transformations.Source: Notation · p. 22
• Π(µ, ν) is the set of couplings between µ ∈ P(X) and ν ∈ P(Y), contained in P(X × Y).Product measures µ × ν and µ^k are defined for repeated-product constructions.Source: Notation · p. 22
• For π ∈ Π(µ, ν), π⊗k is the k-fold product coupling defined by multiplying π(x_i, y_i) across i.The notation also includes coordinatewise map extensions g×n and measure pushforwards f♯µ.Source: Notation · p. 22
• For ξ ∈ P2(X × Y), Σ_ξ denotes its cross-covariance matrix, while Spec denotes the spectrum of a symmetric matrix.Expectation is written E[·], and H2 denotes Hausdorff distance between polytopes under the 2-norm.Source: Notation · p. 22

Pseudo-distance properties

• GW_m is finite, non-negative, symmetric, vanishes on identical measures, and satisfies the triangle inequality, so it is a pseudo-distance on P_2(X).The triangle inequality follows by gluing optimal transport plans.Source: §A.1 · pp. 23–24

Invariance properties

• For centered measures, translating the arguments does not change the objective, and optimal couplings correspond bijectively after transforming one measure by U_α.Multilinearity and centering eliminate the mixed translation terms, while pushforward transports map optimal solutions between the two problems.Source: §A.1 · pp. 24–25

Structural lemma

• If a non-degenerate n-form agrees on a basis under a map T, its unique linear extension is invertible and preserves the multilinear form.The proof establishes linearity and surjectivity of T, then concludes m ◦ T^×n = m.Source: §A.1 · pp. 25–26

Continuous measures

• For compact X and absolutely continuous μ, GW_m(μ,ν)=0 implies that μ and ν are related by a transformation T ∈ G.Optimal transport becomes a Monge map, and measure-one basis arguments combined with the structural lemma show T ∈ GL(d) ∩ G.Source: §A.1 · p. 26

Point clouds

• For distinct equal-weight point clouds with non-zero volume, zero GW_m forces the optimal coupling to be a scaled permutation, induced by a map in G.Non-degeneracy rules out splitting one source point across multiple targets; the resulting permutation preserves m and belongs to GL(d) ∩ G.Source: §A.1 · p. 27

Determinant form

• For the determinant multilinear form, the preserving transformation group is G_det = SL(d).The group is characterized by det = det ◦ T^×d.Source: §A.1 · p. 27

A.2 Proofs of Section 2.2

• The polynomial Q_m is non-negative on all symmetric positive semidefinite matrices and has d^2 variables with degree at most n.This follows because any symmetric positive semidefinite matrix can be realized by a suitable marginal.Source: §A.2 · p. 28
• The proofs show that the projected coupling set is compact and convex because the projection is continuous and the original coupling set is compact and convex.The construction uses coordinate-product functions and, under compactness, an orthogonal projection onto their span.Source: §A.2 · pp. 28–29
• The proofs linearize the multilinear costs, including the determinant-based DGW and CGW terms, for numerical schemes such as Frank-Wolfe.The determinant expression is expanded using first minors and cofactors, while the CGW expression follows by linearity.Source: §A.2 · pp. 29–30

A.3 Low-dimensional structure of the 2-Gromov-Wasserstein

• The coupling-dependent information in the 2-Gromov-Wasserstein problem is represented through second-order coordinate products and squared-norm terms.For centered marginals, the relevant feature space is generated by g(x,y)=||x||^2||y||^2 and f_i,j(x,y)=x_i y_j.Source: §A.3 · p. 31
• The classical 2-Gromov-Wasserstein cost depends polynomially on the coupling projection onto a vector space of dimension d^2 + 1.The projection image of the coupling set remains continuous, convex, and bounded under the stated boundedness and centering assumptions.Source: §A.3 · p. 31
• Because the 2-Gromov-Wasserstein cost is concave in the coupling, a variant of Algorithm 1 computes its global optimum.The low-dimensional reformulation is the basis for the global optimization procedure.Source: §A.3 · p. 31
• The same reformulation enables Frank-Wolfe optimization for a local solution of the 2-Gromov-Wasserstein problem.The coupling-dependent cost is expressed in the projected variables used by the algorithm.Source: §A.3 · p. 31

Auxiliary results

• Couplings remain stable under marginal perturbations: any coupling of (µ1,ν1) can be matched to one of (µ2,ν2) within W2(µ1,µ2) + W2(ν1,ν2).This stability is used to transfer limiting couplings and compare optimization problems with changing marginals.Source: Auxiliary results · pp. 32–33
• Minimizing couplings exist because the coupling set is W2-precompact, remains closed under marginal constraints, and the multilinear cost is continuous.A convergent minimizing subsequence retains its prescribed marginals and achieves the infimum.Source: Auxiliary results · pp. 32–33
• The displayed total second-moment quantity equals m2(π1) + m2(π2), which also equals m2(µ0) + m2(ν0) + m2(µ) + m2(ν) and is at most 4α.Source: Auxiliary results · p. 34
• The multilinear Gromov-Wasserstein objective is robust because the polynomial Qm is Lipschitz on the bounded moment-image ball, with constants determined by its coefficients and α.The resulting objective bounds depend on the polynomial degree and coefficients through the Lipschitz estimate.Source: Auxiliary results · pp. 34–35
• Under sufficiently high q-th moments, empirical multilinear Gromov-Wasserstein convergence follows with a constant depending only on marginal q-th moments, m, and dimension.The required condition is q > 4 when d ≤2 or q > 2d/(d − 2) otherwise.Source: Auxiliary results · p. 35
• An optimal Monge map exists for the multilinear Gromov-Wasserstein problem by reducing its cost through orthogonal transformations and projection to a lower-dimensional transport problem.The reduction uses the rank h of the relevant matrix and a projected cost satisfying the twist condition.Source: Auxiliary results · p. 36
• The map-existence theorem requires compact support, an injective measurable coordinate representation, absolute continuity of the projected marginal, and conditional absolute continuity on fibers.Under these assumptions, an optimal map is induced for the original transport cost.Source: Auxiliary results · pp. 36–37
• For projected costs of the form −⟨ψ1(x),ψ2(y)⟩, absolute continuity of the source yields a unique optimal plan induced by a transport map.The result applies when the measures have compact support and ψ1, ψ2 are diffeomorphisms.Source: Auxiliary results · p. 37

Projection geometry

• Under Assumption 1, Lemma 9 establishes that PΠ contains a Euclidean ball of positive radius r0 depending only on λ, d, and R.The construction takes r0 = r2/rR after establishing r2 > 0 for centered η, ξ ∈ P([−R, R]) with ˜Ση, ˜Σξ ≥ λ.Source: Auxiliary results · pp. 37–38

Bounding-box initialization

• Proposition 8 initializes bounding boxes around PΠ with inner and outer radii whose resulting guarantees depend only on d and the radii r0, r1.The construction produces a rectangular bounding box and a translated inner ball with radius bounded below by a positive quantity depending only on these parameters.Source: Auxiliary results · p. 39

B.2 Proofs of Section 3.3

• For discrete marginals with at most N points, each optimal transport solve costs O(N3), while the iteration count k_f is bounded for precision ϵ.The iteration bound has constants depending on d, R, β, and λ, but not on the marginals.Source: §B.2 · p. 40
• At iteration k, direction selection costs O(k2⌊d2/2⌋), transport solving costs O(N3), and polytope representation updates cost O(k2⌊d2/2⌋+1).The new optimal cost is computed by enumerating polytope vertices with complexity O(k⌊d2/2⌋).Source: §B.2 · pp. 40–41
• The final algorithmic complexity is O(N 3k2⌊d2/2⌋+2 f ), with memory O(N 2 + k⌊d2/2⌋+1).The per-iteration complexity is O(N3k2⌊d2/2⌋+1).Source: §B.2 · p. 41

Complexity for solving the 2-Gromov-Wasserstein distance

• Classical GW2 optimization may fail the volume assumption because PΠ can have zero volume, so the method first computes its effective dimension using supporting hyperplanes and extreme points.For centered unit-sphere supports, gdπ = 1 for every coupling π and PΠ has zero volume.Source: Complexity for solving the… · p. 42
• Algorithm 4 constructs an orthonormal basis of a subspace ˆV containing PΠ, where PΠ has non-zero volume, enabling the approximation algorithm to run in the reduced space.The reduced dimension changes the complexity analysis from d2 to d2 + 1.Source: Complexity for solving the… · p. 42
• Theorem 6 states that Algorithm 5 approximates GW2(µ, ν)2 to precision ϵ for discrete measures with bounded support, with constants depending on d and µ, ν.The theorem provides both a complexity bound and a memory bound, although the supplied passage omits their displayed expressions.Source: Complexity for solving the… · p. 42
• Algorithm 5 computes the reduced basis, initializes primal and dual bounding boxes, iteratively closes their cost gap, and finally returns an optimal coupling from the convex subproblem.The loop continues while |c+ − c−| > ϵ, after which the returned value is C0 + 2c−.Source: Complexity for solving the… · p. 43

Appendix C Convexity of the different costs

• The polynomial costs differ substantially: IGW uses a convex squared-Frobenius term, DGW’s determinant is neither convex nor concave, and CGW combines both terms through parameter t.For CGW, the determinant contribution depends on dimension, requiring separate convexity analysis.Source: Appendix C · p. 44
• For d = 1, the CGW cost is convex for 0 < t ≤1, while d = 2 uses Hessian eigenvalues to characterize convexity.The supplied d = 2 passage states the eigenvalue structure but truncates the resulting condition.Source: Appendix C · p. 44
• For d = 3, determinant-Hessian eigenvalues are bounded using either ||x||∞≤r or numerical bounds based on ||x||2≤r.The coordinate-wise bound is [−4r, 4r], while experiments report α approximately 1.2 for the Euclidean-norm bound.Source: Appendix C · p. 44
• Choosing t from bounding-box radius bounds ensures the CGW cost has the desired curvature on Pπ, with sufficient conditions based on R or an upper bound r.The supplied passages give t ≥ 24R^2/(2+24R^2), an experimental bound t ≥ 8R^2/(2+8R^2), and the analogous r-based condition.Source: Appendix C · p. 44

Appendix D Algorithm, concave cost

• For concave costs, the appendix initializes bounding boxes by solving directional problems, selecting an optimal direction and cost, then iteratively refining bounds until convergence.The main algorithm returns an optimal coupling after solving the convex problem for the selected direction, with an optional Frank–Wolfe refinement.Source: Appendix D · pp. 45, 47, 48
• An alternative bounding-box initialization is practical and appears effective, but neither it nor the practical direction choice guarantees the convergence rate of Theorem 3.The appendix explicitly distinguishes practical performance from the theorem’s convergence guarantee.Source: Appendix D · p. 47

Appendix E Frank-Wolfe line search for the IGW, DGW and CGW problems

• The IGW Frank–Wolfe line search reduces to endpoint selection because its one-dimensional objective is convex, yielding solutions at τ = 0 or τ = 1.The appendix provides a dedicated IGW line-search algorithm.Source: Appendix E · p. 49
• For DGW and CGW, Frank–Wolfe line searches optimize polynomial objectives: quadratic in two dimensions and cubic for the determinant term in three dimensions.The 2D procedure compares endpoints and an interior critical point when applicable; the 3D procedure evaluates the cubic’s candidate maximizer and endpoints.Source: Appendix E · pp. 49, 51

Experimental setup

• Experiments use equally weighted marginal points randomly drawn from a standard Gaussian distribution.This sampling setup applies to Figures 4, 5, and 6.Source: Appendix F · p. 53

Figure 4

• In dimensions 2 and 3, bounding-box complexities remain below the theoretical complexity from Theorem 3.Figure 4 evaluates vertex and constraint counts in the upper and lower bounding boxes together with the global-optimality certificate |c+ −c−|.Source: Appendix F · p. 53

Figure 5

• Algorithm 3 converges in less than 30 iterations for both IGW and DGW costs, measured by squared cost error and transport-plan error.The experiments assess convergence in both objective cost and transport-plan accuracy.Source: Appendix F · p. 53

Figure 6

• Our cut-selection method converges continuously, with linear outer-bounding-box vertex growth versus quadratic growth for Ryner et al.’s method [35].Ryner’s method can stagnate indefinitely from machine imprecision, whereas our method converges exactly in finite (exponential) time for finite discrete measures.Source: Appendix F · p. 53

• The lower bound c− converges faster than the upper bound c+ when estimating the CGW cost, a recurrent behavior in simulations.The gap |c+ − c−| bounds the precision of the approximate global solution.Source: Figure 3 · p. 54, Appendix G · p. 55
• For IGW global optimization, the vertex bound grows as O(k^2) in dimension 2 and O(k^4) in dimension 3, with practical dimension-3 growth near O(k^3).In dimension 2, the algorithm runs several thousand iterations on an ordinary laptop without memory issues.Source: Appendix G · p. 55
• The Frank–Wolfe scheme converges quickly for both concave IGW and DGW costs in dimension 3, typically requiring fewer than 50 iterations.Plan error is measured by Euclidean 2-norm, while cost error compares final and per-iteration squared costs.Source: Figure 5 · p. 56
• For IGW global optimization, the supplementary comparison evaluates the proposed facet-based search direction against Ryner et al.’s gradient-based strategy [35].The proposed strategy selects the facet of P−Π facing x+, whereas Ryner et al.’s strategy uses the gradient of Qm at x+.Source: Figure 6 · p. 57, Appendix G · p. 57