Non-negative Alternate Scaling (NASA)

• Non-negative Alternate Scaling (NASA) is an iterative method that transforms non-negative matrices into doubly stochastic matrices via alternating row and column normalizations.
• The algorithm generalizes classical matrix scaling to incorporate Bregman divergences and non-negativity constraints, enabling applications in regularized optimal transport and statistical estimation.
• Supported by the Sinkhorn–Knopp theorem, NASA exhibits geometric convergence with O(n²) operations per iteration, making it efficient for high-dimensional and structured data problems.

The Non-negative Alternate Scaling Algorithm (NASA), also known as the Sinkhorn algorithm or the Sinkhorn-Knopp iteration, is a foundational class of iterative methods for transforming a non-negative (or strictly positive) matrix into a doubly stochastic matrix via alternating normalizations. Generalizations of NASA play a central role in scaling problems, optimal transport, matrix analysis, and algorithmic applications, including polynomial-time permanent approximation, design matrix analysis, statistical estimation, and regularized optimal transport. NASA encompasses both classical matrix scaling and modern extensions to general Bregman divergences and non-automatic non-negativity constraints.

1. Problem Formulation and Definitions

Given a non-negative matrix $A \in \mathbb{R}^{n \times n}_{\ge0}$, doubly stochasticity requires that all row and column sums equal 1: $$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$ A key functional to measure deviation from a doubly stochastic matrix is

$$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$

where $r_i$ and $c_j$ denote the $i$-th row sum and $j$-th column sum, respectively.

The $\varepsilon$-scaling problem asks: for given $\varepsilon > 0$, find diagonal $B, C \succ 0$ such that $$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$0 (where possible). Scalability—existence of such scaling for every $$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$1—is equivalent to positivity of the matrix permanent and to the bipartite support graph having a perfect matching (Garg et al., 2018).

For regularized optimal transport (ROT), the NASA is applied to solve problems of the form

$$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$2

where $$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$3 is the set of couplings (joint distributions) with prescribed marginals $$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$4, $$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$5 is a cost matrix, and $$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$6 is a smooth convex regularizer, possibly not restricting the domain to the non-negative orthant (Dessein et al., 2016).

2. Iterative Algorithmic Structure

The classical NASA iterates between affine projections that alternately normalize rows and columns:

• Row-scaling: For a matrix $$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$7, define

$$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$8

• Column-scaling: Likewise,

$$\sum_{j=1}^n A_{ij} = 1 \quad \forall i,\qquad \sum_{i=1}^n A_{ij} = 1 \quad \forall j$$9

The update sequence is then: $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$0 This sequence converges to a unique doubly stochastic scaling $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$1, where $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$2 are positive diagonals (Nathanson, 2019, Nathanson, 2019).

For convex regularized transport, NASA generalizes to compute Bregman projections in dual coordinates. When the regularizer’s domain does not guarantee non-negativity, explicit projections onto the non-negative orthant are required. Corrections (Dykstra’s algorithm) are maintained to ensure convergence when projecting cyclically onto the constraints:

• $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$3 (non-negativity)
• $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$4 (row marginals)
• $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$5 (column marginals) NASA cycles as $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$6, with Newton-type inner steps for affine projections in the case of separable $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$7 (Dessein et al., 2016).

3. Convergence Theory and Complexity

The classical Sinkhorn–Knopp theorem asserts that for every strictly positive $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$8, the NASA produces positive diagonal matrices $$\operatorname{ds}(A) = \sum_{i=1}^n (r_i - 1)^2 + \sum_{j=1}^n (c_j - 1)^2$$9 such that $r_i$0 is uniquely doubly stochastic. The iterative map is non-expansive in the Hilbert metric on the positive cone, and the process exhibits geometric (linear) convergence (Nathanson, 2019, Nathanson, 2019).

In the entropy-regularized setting, convergence is seen as monotonic minimization of the Kullback–Leibler divergence $r_i$1 over diagonal $r_i$2 scalings, and each iteration strictly decreases a potential. The number of NASA iterations needed to achieve $r_i$3 is $r_i$4, assuming the matrix entries have bit complexity $r_i$5. Each iteration incurs $r_i$6 arithmetic operations (Garg et al., 2018).

For general Bregman divergences, under mild technical conditions (Legendre type, cofinite, full domain for conjugate, etc.), Dykstra’s algorithm ensures unique Bregman projection convergence at a linear rate (Theorem 2.1 of (Dessein et al., 2016)). For separable $r_i$7, each outer iteration has $r_i$8 complexity, dominated by vectorized Newton updates. For non-separable regularizers (e.g., Mahalanobis), complexity grows to $r_i$9. Empirical studies confirm practical quadratic scaling for moderate problem sizes (Dessein et al., 2016).

4. Algorithmic Pseudocode and Procedures

The concrete NASA scheme for matrix scaling (Sinkhorn’s algorithm) is as follows (Garg et al., 2018, Nathanson, 2019):

1. If any row or column of $c_j$0 is all-zero, terminate and report “not scalable.”
2. Initialize $c_j$1, $c_j$2.
3. For $c_j$3:
   • Compute $c_j$4 and $c_j$5.
   • If $c_j$6, return $c_j$7.
   • If $c_j$8, set $c_j$9; else $i$0.
4. Output “not scalable” if no feasible scaling is found.

For regularized optimal transport:

• NASA performs in the dual parameterization $i$1, maintaining correction variables for each constraint.
• Each cycle consists of projections (generally using Newton–Raphson for affine marginals), interleaved with non-negativity projections applied elementwise (Dessein et al., 2016).

Pseudocode details for both cases are explicitly provided in (Garg et al., 2018) and (Dessein et al., 2016).

5. Applications and Illustrative Use Cases

The NASA has a spectrum of applications across mathematics, statistics, and machine learning:

• Permanent Approximation: Any $i$2 with $i$3 can be nearly doubly-stochastically scaled. Since $i$4 for a doubly stochastic $i$5 satisfies $i$6, NASA delivers a deterministic $i$7-approximation to $i$8 in polynomial time (Garg et al., 2018).
• Rank Bounds in Design Matrices: Scaling a design matrix to (near) doubly stochasticity provides lower bounds on matrix rank, crucial for quantitative Sylvester–Gallai theorems and improved bounds in combinatorial geometry (Garg et al., 2018).
• Maximum-Entropy Estimation: The scaling problem is equivalent to entropy-regularized convex optimization for contingency tables and joint distribution estimation, with NASA iterations realized as block coordinate minimizations (Garg et al., 2018, Dessein et al., 2016).
• Regularized Optimal Transport (ROT): In problems where regularization is not entropic, NASA enables efficient solution via explicit alternated projections, supporting sparsity, additional constraints, and richer regularization (Dessein et al., 2016). Experiments include synthetic 1D distributions (Gaussian and bimodal), runtime comparisons for various regularizers, and real-world audio scene classification with improved classification rates using the NASA-optimized kernel (Dessein et al., 2016).

6. Special Cases, Finite-Termination Phenomena, and Limitations

Classically, NASA achieves convergence only asymptotically. Recent work constructs explicit families of $i$9 positive matrices that terminate after finitely many steps: for each $j$0, the family of non-invertible, row-stochastic but not column-stochastic matrices constructed by Nathanson become doubly stochastic after a single column-scaling (yet no earlier) (Nathanson, 2019). All such examples have determinant zero, and it is open whether a strictly positive, invertible $j$1 can have finite NASA termination.

Further, explicit Sinkhorn limits have been calculated for all symmetric $j$2 matrices with two values, providing links to Diophantine approximation and sequences of rational numbers converging to algebraic irrationals via NASA iterates (Nathanson, 2019).

7. Practical Implementations and Advanced Considerations

NASA is well-suited for high-dimensional, structured, or sparse data due to its operation via vectorized row/column-wise updates, which are GPU/parallel-friendly (Dessein et al., 2016). For problems with infinite cost entries, NASA naturally restricts computation to active submatrices. Choice of regularization controls solution smoothness and sparsity, with larger $j$3 promoting faster convergence and increased smoothing. Termination criteria are application-dependent, involving either direct error in margin deviation or iterations on dual variables.

In the separable case, practical performance matches the theoretical $j$4 scaling, outperforming classical, combinatorial Earth Mover’s Distance solvers (which are super-cubic). For non-separable $j$5, NASA remains applicable but computationally more intensive.

NASA provides a unifying approach to matrix scaling and regularized transport, extending the reach of classical iterative scaling to general convex frameworks and broad application domains (Garg et al., 2018, Dessein et al., 2016, Nathanson, 2019, Nathanson, 2019).