Accelerating operator Sinkhorn iteration with overrelaxation

Technical Summary: Accelerating Operator Sinkhorn Iteration with Overrelaxation

Problem Statement
The paper addresses the operator scaling problem, which seeks to find invertible square matrices $L \in GL_m(\mathbb{R})$ and $R \in GL_n(\mathbb{R})$ that jointly transform a given collection of matrices $A_1, \dots, A_k \in \mathbb{R}^{m \times n}$ into a "scaled" form $\bar{A}_i = L A_i R^\top$. The goal is to satisfy two simultaneous conditions:

1. $\sum_{i=1}^k \bar{A}_i \bar{A}_i^\top = \frac{1}{m} I_m$
2. $\sum_{i=1}^k \bar{A}_i^\top \bar{A}_i = \frac{1}{n} I_n$

This problem generalizes classical matrix scaling and arises in applications such as non-commutative polynomial identity testing, invariant theory, and signal processing. The standard solution method is the operator Sinkhorn iteration (OSI), an alternating fixed-point algorithm. While OSI is globally convergent under certain conditions, its convergence can be slow. The authors aim to accelerate this process using Successive Overrelaxation (SOR), a technique well-established for linear systems but less explored for this specific nonlinear operator scaling context.

Methodology
The authors propose three distinct variants of overrelaxation for the operator Sinkhorn iteration, all based on the interpretation of the problem as an alternating minimization on cones of positive definite (PD) matrices.

1. Euclidean Overrelaxation (PD-OR): This approach applies standard affine combinations to the iterates $(X_t, Y_t)$, where $X_t = L_t^\top L_t$ and $Y_t = R_t^\top R_t$. The update rule is $(1-\omega)X_t + \omega S_1(Y_t)$, where $S_1$ is the scaling map. While computationally simple, this method lacks theoretical guarantees that the resulting matrices remain positive definite when the relaxation parameter $\omega > 1$.
2. Coordinate Transformation Overrelaxation (Cholesky-OR): To ensure the iterates remain in the PD cone, the authors suggest applying overrelaxation in a transformed coordinate space. They specifically utilize the Cholesky factors ($L$ and $R$) of the PD matrices. Since the Cholesky factors belong to linear spaces (triangular matrices), affine combinations are well-defined and preserve the structure required for the algorithm. This variant is presented as the most numerically efficient.
3. Geodesic Overrelaxation (Geodesic-OR): This method respects the hyperbolic geometry of the PD cones induced by the Hilbert metric. Instead of affine combinations, the authors define updates along geodesics using the operator geometric mean $X \#_\omega \tilde{X}$. This approach is theoretically robust, preserving the orbital equivariance of the problem and ensuring iterates remain within the PD cone.

Key Contributions and Theoretical Results

• Local Convergence Analysis: The authors demonstrate that near a fixed point, all three overrelaxation variants are equivalent to the first order. By linearizing the iterations and utilizing the orbital equivariance (invariance under scaling of the solution orbit), they derive a unified formula for the asymptotic local convergence rate.

  • They establish that the spectral radius of the linearized SOR iteration follows the classic Young's SOR theorem pattern.
  • They provide a formula for the asymptotically optimal relaxation parameter $\omega_{opt}$, which depends on the spectral radius $\beta$ of the standard (unrelaxed) operator Sinkhorn iteration: $\omega_{opt} = \frac{2}{1 + \sqrt{1-\beta^2}}$.
  • This analysis relies on the fact that the Hessian of the underlying geodesically convex objective function has a one-dimensional null space corresponding to the scaling orbit, a property they rigorously prove.

• Global Convergence for Geodesic SOR: While global convergence for nonlinear SOR is generally difficult, the authors prove global convergence for the Geodesic-OR variant within a specific range of $\omega$.

  • They utilize the Hilbert metric and the contractivity of the standard Sinkhorn iteration.
  • A novel contribution is a metric inequality for the Hilbert metric on PD matrices (Lemma 3.4), which bounds the distance of a geodesic point from a third point. This inequality allows them to extend contractivity results to the overrelaxed case, proving convergence for $\omega \in (0, \frac{2}{1+\sqrt{\Lambda_1 \Lambda_2}})$, where $\Lambda_1, \Lambda_2$ are Lipschitz constants of the scaling maps.

• Adaptive Parameter Selection: Recognizing that the optimal $\omega$ depends on the unknown convergence rate of the specific instance, the authors adopt an adaptive strategy. They estimate the convergence rate $\beta$ from the initial iterations of the standard method and dynamically update $\omega$ to $\omega_{opt}$, a technique adapted from recent work on matrix scaling.

Numerical Results
The authors conducted experiments on two types of instances:

1. Frame Scaling: In this application, the SOR methods (particularly Cholesky-OR and Geodesic-OR) demonstrated significantly faster convergence than the standard OSI. The methods achieved errors of order $10^{-14}$ in roughly half the iterations required by OSI to reach $10^{-8}$. The observed convergence rates matched the theoretical predictions closely.
2. Ill-Conditioned Operators: Using Hilbert matrices as inputs, the authors observed that while SOR accelerated convergence initially, it stagnated at a lower accuracy ($10^{-6}$) compared to the standard OSI ($10^{-11}$). They attribute this to the fact that SOR variants based on "half-scaled" intermediate matrices do not fully compensate for the initial ill-conditioning as effectively as the full alternating updates in the standard OSI.

Significance and Claims
The paper claims to generalize existing results on matrix scaling (specifically works by Thibault et al. and Lehmann et al.) to the more complex setting of operator scaling.

• Theoretical Insight: The authors emphasize that their rigorous local convergence analysis, which handles the non-positive definite nature of the Hessian due to group equivariance, provides new insight into the structure of operator scaling problems.
• Practical Utility: The numerical experiments suggest that overrelaxation can accelerate operator scaling by several orders of magnitude in well-conditioned scenarios (like frame scaling) with negligible computational overhead, particularly when using the Cholesky-based approach.
• Limitations: The authors modestly acknowledge that their global convergence proof for the geodesic version covers a range of $\omega$ that is likely too pessimistic (often close to 1), and that the current SOR implementations may suffer from reduced numerical accuracy on ill-conditioned inputs. They identify establishing global convergence for a wider range of $\omega$ and improving accuracy on ill-conditioned problems as directions for future work.