On Sinkhorn's Algorithm and Choice Modeling

On Sinkhorn's Algorithm and Choice Modeling

\RUNAUTHOR{Qu, Galichon, Gao, Ugander}

\RUNTITLE{Sinkhorn's Algorithm and Choice Modeling}

\TITLE{On Sinkhorn's Algorithm and Choice Modeling}

\ARTICLEAUTHORS{ \AUTHOR{Zhaonan Qu} \AFF{Data Science Institute, Columbia University, New York, USA \EMAIL{[email removed]}} \AUTHOR{Alfred Galichon} \AFF{Department of Mathematics and Department of Economics, New York University, New York, USA\\ Department of Economics, Sciences Po, Paris, France \EMAIL{[email removed]}} \AUTHOR{Wenzhi Gao} \AFF{Institute for Computational and Mathematical Engineering, Stanford University, California, USA \EMAIL{[email removed]}} \AUTHOR{Johan Ugander} \AFF{Department of Management Science & Engineering, Stanford University, California, USA \EMAIL{[email removed]}} }

\ABSTRACT{For a broad class of models widely used in practice for choice and ranking data based on Luce's choice axiom, including the Bradley--Terry--Luce and Plackett--Luce models, we show that the associated maximum likelihood estimation problems are equivalent to a classic matrix balancing problem with target row and column sums. This perspective opens doors between two seemingly unrelated research areas, and allows us to unify existing algorithms in the choice modeling literature as special instances or analogs of Sinkhorn's celebrated algorithm for matrix balancing. We draw inspirations from these connections and resolve some open problems on the study of Sinkhorn's algorithm. We establish the global linear convergence of Sinkhorn's algorithm for non-negative matrices whenever finite scaling matrices exist, and characterize its linear convergence rate in terms of the algebraic connectivity of a weighted bipartite graph. We further derive the sharp asymptotic rate of linear convergence, which generalizes a classic result of Knight (2008). To our knowledge, these are the first quantitative linear convergence results for Sinkhorn's algorithm for general non-negative matrices and positive marginals. Our results highlight the importance of connectivity and orthogonality structures in matrix balancing and Sinkhorn's algorithm, which could be of independent interest. More broadly, the connections we establish in this paper between matrix balancing and choice modeling could also help motivate further transmission of ideas and lead to interesting results in both disciplines. } \KEYWORDS{Luce Choice Models, Matrix Balancing, Algebraic Connectivity, Linear Convergence}

Introduction

The modeling of choice and ranking data is an important topic across many disciplines. Given a collection of $m$ objects, a universal problem is to aggregate choice or partial ranking data over them to arrive at a reasonable description of either the behavior of decision makers, the intrinsic qualities of the objects, or both. Work on such problems dates back over a century at least to the work of Landau, who considered $m$ chess players and a record of their match results against one another, aiming to aggregate the pairwise comparisons to arrive at a global ranking of all players landau1895relativen,elo1978rating. More generally, comparison data can result from choices from subsets of varying sizes, from partial or complete rankings of objects, or from mixtures of different data types.

The modern rigorous study of comparisons primarily builds on the foundational works of thurstone1927method and zermelo1929berechnung. Both proposed models based on a numerical “score” for each item (e.g., chess player), but with different specifications of choice probabilities. zermelo1929berechnung builds on the intuition that choice probability should be proportional to the score, and proposes a iterative algorithm to estimate the scores from pairwise comparison data. As one of the foundational works in this direction, luce2012individual formalized the multinomial logit model of discrete choice starting from the axiom of independence of irrelevant alternatives (IIA). It states that the relative likelihood of choosing an item $j$ over another item $k$ is independent of the presence of other alternatives. In other words, if $S$ and $S'$ are two subsets of the $m$ alternatives, both containing $j$ and $k$, and $\Pr(j, S)$ denotes the probability of choosing item $j$ from $S$, then for $\Pr(k, S),\Pr(k, S')>0$,

This invariance property, together with a natural condition for zero probability alternatives, are often referred to as Luce's choice axioms. They guarantee that each alternative can be summarized by a non-negative score $s_{j}$ such that the probability of choice can be parameterized by

for any set $S$ that contains $j$. The parameters $s_{j}$ reflect the “intrinsic” value of item $j$, and are unique up to a normalization, which can be set to $\sum_{j}s_{j}=1$. This general choice model includes as a special case the Bradley--Terry--Luce model (BTL) for pairwise comparisons bradley1952rank, and also applies to ranking data when each $k$-way ranking is broken down into $k-1$ choice observations, where an item is chosen over the set of items ranked lower critchlow1991probability,plackett1975analysis,hausman1987specifying. The many subsequent works that build on Luce's choice axioms speak to its fundamental importance in choice modeling. Other works have also sought to address the limitations of Luce choice models. Prominent among them are probit models thurstone1927method,berkson1944application, random utility (RUM) models mcfadden2000mixed, context-dependent (CDM) models batsell1985new,seshadri2020learning, and behavioral models from psychology tversky1972elimination.

Matrix balancing (or scaling), meanwhile, is a seemingly unrelated mathematical problem with an equally long history. In its most common form that we study in this paper, the problem seeks positive diagonal matrices $D^0,D^1$ of a given (entry-wise) non-negative matrix $A\geq0$, such that the scaled matrix $D^0AD^1$ has row and column sums equal to some prescribed positive marginals $p,q$:

Over the years, numerous applications and problems across different domains, including statistics yule1912methods,deming1940least,ireland1968contingency, economics stone1962multiple,bacharach1970biproportional,galichon2021matching, transportation networks kruithof1937telefoonverkeersrekening,lamond1981bregman,chang2024inferring, optimization bregman1967relaxation,ruiz2001scaling, and machine learning cuturi2013sinkhorn,peyre2019computational, have found themselves essentially solving a new incarnation of the old matrix balancing problem, which attests to its universality and importance.

A major appeal of the matrix balancing problem lies in the simplicity and elegance of its popular solution method, widely known as Sinkhorn's algorithm sinkhorn1964relationship. Observe that it is easy to scale the rows or columns of $A$ such that the resulting matrix satisfies one of the two marginal constraints in (ref). However, it is more difficult to construct scalings $D^0,D^1$ that simultaneously satisfy both constraints. Sinkhorn's algorithm ((ref)) simply alternates between updating the scalings $D^0$ and $D^1$ to satisfy one of the two marginal conditions in the hope of converging to a solution, leading to lightweight implementations that have proved effective for practical problems of massive size. In particular, Sinkhorn's algorithm has gained much popularity in the recent decade thanks to its empirical success at approximating optimal transport (OT) distances cuturi2013sinkhorn,altschuler2017near, which are bedrocks of important recent topics in operations research, such as Wasserstein distributionally robust optimization mohajerin2018data,gao2023distributionally,kuhn2019wasserstein,blanchet2019robust,blanchet2022distributionally.

Despite the widespread popularity of Sinkhorn's algorithm, its convergence behavior is yet to be fully understood. In particular, while there have been extensive studies of convergence, many focus on the setting when the matrix $A>0$, i.e., entry-wise strictly positive, which includes most OT problems. In contrast, other applications of matrix balancing, particularly those with network structures, have $A\geq0$ with zero elements, and is therefore potentially sparse. In this setting, quantitative analyses are less common and more fragmented, employing different assumptions whose connections and distinctions remain less clear. On one hand, works such as kalantari2008complexity, chakrabarty2021better, and leger2021gradient have established global (that is, true for all iterations $t\geq1$) sub-linear convergence results, i.e., convergence to an $\varepsilon$ accuracy solution requires a total number of iterations that is polynomial in $1/\varepsilon$. On the other hand, knight2008sinkhorn establishes local (and more specifically asymptotic) linear convergence for square matrix $A\geq 0$ and uniform marginals $p,q$. In other words, as $t\rightarrow \infty$, solution accuracy at iteration $t+1$ improves over that at iteration $t$ with a constant factor. Furthermore, a general result in luo1992convergence implies global linear convergence of Sinkhorn's algorithm, i.e., convergence to an $\varepsilon$ accuracy solution requires iterations polynomial in $\log(1/\varepsilon)$. However, their result is implicit and does not characterize the dependence on problem parameter and structure as those in the sub-linear results.

These results leave open several questions on the convergence of Sinkhorn's algorithm. First, when does a quantitative global linear convergence result exist for $A\geq 0$, and how to characterize the global convergence rate in terms of the problem primitives? Second, how to characterize the sharp, i.e., best possible, asymptotic linear convergence rate $\lambda$ that is applicable to general non-negative $A$ and non-uniform $p,q$? Third, how to reconcile and clarify the linear vs. sub-linear convergence results under different assumptions on the problem structure? Given that many applications of matrix balancing with network structures correspond to the setting with sparse $A\geq0$, such as transportation and trade, it is therefore important to better understand the convergence of Sinkhorn's algorithm in this setting.

In this paper, we provide answers to these open questions in matrix balancing. Surprisingly, the inspirations for our solutions come from results in the seemingly unrelated topic of choice modeling. Our main contributions are summarized below.

Our first set of contributions, which we detail in (ref), is recognizing Luce choice models as yet another instance where a central problem reduces to that of matrix balancing. More precisely, we formally establish the equivalence between the maximum likelihood estimation of Luce choice models and matrix balancing problems with an $A\geq 0$ with binary elements ((ref)). We also clarify the relations and distinctions between problem assumptions in the two literatures ((ref)). More importantly, we demonstrate that classic and new algorithms from the choice literature, including those of zermelo1929berechnung,dykstra1956note,ford1957solution,hunter2004mm,maystre2017choicerank,agarwal2018accelerated, can be viewed as special cases or analogs of Sinkhorn's algorithm when applied to various problems in the choice setting ((ref)). These intimate mathematical and algorithmic connections allow us to provide a unifying perspective on works from both areas. More broadly, they enable researchers to import insights and tools from one domain to the other. In particular, recent works on choice modeling shah2015estimation,seshadri2020learning,vojnovic2020convergence have highlighted the importance of algebraic connectivity fiedler1973algebraic,spielman2007spectral of the data structure for efficient parameter learning, which motivates us to also consider this quantity in the convergence analysis of Sinkhorn's algorithm.

Our next set of contributions, detailed in (ref), is establishing novel convergence bounds on Sinkhorn's algorithm, drawing from the connections to choice modeling that we establish. First, we provide a global linear convergence bound for Sinkhorn's algorithm, whenever the matrix balancing problem has a finite solution pair $D^0,D^1$ ((ref)). We characterize the global convergence rate in terms of the algebraic connectivity of the weighted bipartite graph whose biadjacency matrix is precisely $A$. To our knowledge, this result is the first to highlight the fundamental role of algebraic connectivity in the study of matrix balancing with sparse matrices. In addition, we characterize the asymptotic linear rate of convergence in terms of the scaled matrix $D^0AD^1$ with target marginals $p,q$, generalizing a result of knight2008sinkhorn for uniform marginals and square matrices ((ref)). This result employs a more explicit analysis that exploits an intrinsic orthogonality structure of Sinkhorn's algorithm. We also clarify the convergence behavior of Sinkhorn's algorithm under two regimes: when a finite scaling pair $D^0,D^1$ exists, Sinkhorn's algorithm converges linearly; otherwise, it only converges sub-linearly under the minimal conditions required for convergence ((ref)).

{Besides the contributions above, we further discuss connections between Sinkhorn's algorithm and topics in optimization and choice modeling in (ref). For example, interpreting Sinkhorn's algorithm as a distributed optimization algorithm on a bipartite graph helps explain the importance of the spectral properties of the graph on its convergence. Inspired by Bayesian regularization of Luce choice models using gamma priors, we also design a regularized Sinkhorn's algorithm in (ref), that is guaranteed to converge even when the standard algorithm does not, which is not uncommon when the data is very sparse and there are measurement errors.}

We believe that the connections we establish in this paper between choice modeling and matrix balancing can lead to further interesting results in both disciplines, and are therefore relevant to researchers working on related topics. In particular, the fundamental role of algebraic connectivity in the study of matrix balancing for sparse matrices goes beyond quantifying the algorithmic efficiency of Sinkhorn's algorithm. See, for example, chang2024inferring, which quantifies the statistical efficiency of a network traffic model using algebraic connectivity.

Related Work

This section includes an extensive review of related works in choice modeling and matrix balancing. Well-versed readers may skip ahead to the mathematical preliminaries ((ref)) and our core results ((ref)).

Choice Modeling

Methods for aggregating choice and comparison data usually take one of two closely related approaches: maximum likelihood estimation of a statistical model or ranking according to the stationary distributions of a random walk on a Markov chain. Recent connections between maximum likelihood and spectral methods have put these two classes of approaches in increasingly close conversation with each other.

Spectral Methods. The most well-known spectral method for rank aggregation is perhaps the PageRank algorithm page1999pagerank, which ranks web pages based on the stationary distribution of a random walk on a hyperlink graph. The use of stationary distributions also features in the work of dwork2001rank, the Rank Centrality (RC) algorithm negahban2012iterative,negahban2016rank, which generates consistent estimates for the Bradley--Terry--Luce pairwise comparison model under assumptions on the sampling frame, and the Luce Spectral Ranking (LSR) and iterative LSR (I-LSR) algorithms of maystre2015fast for choices from pairs as well as larger sets. Following that work, agarwal2018accelerated proposed the Accelerated Spectral Ranking (ASR) algorithm with provably faster mixing times than RC and LSR, and better sample complexity bounds than negahban2016rank. knight2008sinkhorn is an intriguing work partially motivated by page1999pagerank that applies Sinkhorn's algorithm, which is central to the current work, to compute authority and hub scores similar to those proposed by kleinberg1999authoritative and tomlin2003new, although the focus in knight2008sinkhorn is on Markov chains rather than maximum likelihood estimation of choice models. For ranking data, soufiani2013generalized decompose rankings into pairwise comparisons and develop consistent estimators for Plackett--Luce models based on a generalized method of moments. Other notable works that make connections between Markov chains and choice modeling include blanchet2016markov and ragain2016pairwise.

Maximum Likelihood Methods. Maximum likelihood estimation of the Bradley--Terry--Luce model dates back to zermelo1929berechnung, dykstra1956note, and ford1957solution, which all give variants of the same iterative algorithm and prove its convergence to the MLE when the directed comparison graph is strongly connected. Much later, hunter2004mm observed that their algorithms are instances of a class of minorization-maximization or majorization-minimization (MM) algorithms, and develops MM algorithms for the Plackett--Luce model for ranking data, among others. vojnovic2020convergence,vojnovic2023accelerated further investigate the convergence rate of MM algorithm for choice models, quantifying it in terms of the algebraic connectivity of the comparison graph. newman2023efficient proposes an alternative to the classical iterative algorithm for pairwise comparisons based on a reformulated moment condition, achieving impressive empirical speedups. negahban2012iterative is arguably the first work that connects maximum likelihood estimation to Markov chains, followed by maystre2015fast, whose spectral method is based on a balance equation interpretation of the optimality condition. kumar2015inverting consider the problem of inverting the stationary distribution of a Markov chain, and embed the maximum likelihood problem of the Luce choice model into this framework, where the MLEs parameterize the desired transition matrix. maystre2017choicerank consider the estimation of a network choice model with similarly parameterized random walks. Lastly, a vast literature in econometrics on discrete choice also considers different aspects of the ML estimation problem. In particular, the present paper is related to the Berry--Levinsohn--Pakes (BLP) framework of berry1995automobile, well-known in econometrics. The matrix balancing interpretation of maximum likelihood estimation of choice models that we develop in this paper connects many of the aforementioned works.

Besides optimization problems related to maximum likelihood estimation, there have also been extensive studies on the statistical properties of maximum likelihood estimates themselves hajek2014minimax,rajkumar2014statistical. In particular, a line of recent works have highlighted the importance of {algebraic connectivity}---as quantified by the Fiedler eigenvalue fiedler1973algebraic,spielman2007spectral---on the statistical efficiency of the MLEs. shah2015estimation is the first to recognize this significance of data structure for the statistical efficiency of the BTL model, which they refer to as “topology-dependence”. As a byproduct of analysis for a context-dependent generalization of the Luce choice model, seshadri2020learning obtain tight expected risk and tail risk bounds for the MLEs of Luce choice models (which they call MNL) and Plackett--Luce ranking models in terms of the algebraic connectivity, extending and improving upon previous works by hajek2014minimax,shah2015estimation,vojnovic2016parameter. Other works with tight risk bounds on the BTL model include hendrickx2020minimax and bong2022generalized, who also provide the first high probability guarantees for the existence of finite MLEs of the BTL model under conditions on a Fiedler eigenvalue. Interestingly, the statistical significance of algebraic connectivity has also been highlighted in models of networks in econometrics and machine learning by the works of de2017econometrics,jochmans2019fixed,chang2024inferring, among others. Our present work is primarily concerned with the optimization aspects of the maximum likelihood estimation of choice models. Nevertheless, the statistical importance of algebraic connectivity in the aforementioned works also provide motivations for us to leverage it in our convergence analysis of Sinkhorn's algorithm for matrix balancing.

Lastly, a short note on terminology. Even though a choice model based on (ref) is technically a “multinomial logit model” with only intercept terms mcfadden1973conditional, there are subtle differences. When (ref) is applied to model ranking and choice data with distinct items, each observation $i$ usually consists of a possibly different subset $S_i$ of the universe of all alternatives, so that there is a large number of different configurations of the choice menu in the dataset. On the other hand, common applications of multinomial logit models, such as classification models in statistics and machine learning bishop2006pattern and discrete choice models in econometrics mcfadden1973conditional, often deal with repeated observations consisting of the same number of alternatives. However, these alternatives now possess “characteristics” that vary across observations, which are often mapped {parametrically} to the scores in (ref). In this paper, we primarily use the term Luce choice model to refer to the model (ref), although it is also called MNL (for multinomial logit) models in some works. We refrain from using the term MNL to avoid confusion with parametric models for featurized items used in ML and econometrics.

Matrix Balancing

Matrix balancing (or scaling) is an important topic in optimization and numerical linear algebra that underlies a diverse range of applications. The particular question of scaling rows and columns of a matrix $A$ so that the resulting matrix has target row and column norms $p,q$ has been studied as early as the 1930s, and continue to interest researchers from different disciplines today. The present paper only contains a partial survey of the vast literature on this topic. (ref) provides a summary of some popular applications to illustrate the ubiquity of the matrix balancing problem. schneider1990comparative and idel2016review also provide excellent discussions of many applications.

The standard iterative algorithm for the matrix balancing problem we study in this paper has been rediscovered independently quite a few times. As a result, it has domain-dependent names, including the iterative proportional fitting (IPF) procedure deming1940least, biproportional fitting bacharach1965estimating and the RAS algorithm stone1962multiple, but is perhaps most widely known as Sinkhorn's algorithm sinkhorn1964relationship. A precise description can be found in (ref). Sinkhorn's algorithm is also closely related to relaxation and coordinate descent type methods for solving the dual of entropy optimization problems bregman1967relaxation,cottle1986lagrangean,tseng1987relaxation,luo1992convergence, as well as message passing and belief propagation algorithms in distributed optimization balakrishnan2004polynomial,agarwal2018accelerated.

The convergence behavior of Sinkhorn's algorithm in different settings has been extensively studied by sinkhorn1964relationship,bregman1967proof,lamond1981bregman,franklin1989scaling,ruschendorf1995convergence,kalantari2008complexity,knight2008sinkhorn,pukelsheim2009iterative,altschuler2017near,dvurechensky2018computational,marino2020optimal,chakrabarty2021better,leger2021gradient,carlier2022linear, among many others. For $A$ with strictly positive entries, franklin1989scaling establish the global linear convergence of Sinkhorn's algorithm in the Hilbert projective metric $d$ bushell1973hilbert. More precisely, if $r^{(t)}$ denotes the row sum of the scaled matrix after $t$ iterations of Sinkhorn's algorithm that enforce column constraints, then

for some $\lambda \in (0,1)$ dependent on $A$. On the other hand, works such as kalantari1993rate, altschuler2017near, and dvurechensky2018computational develop complexity bounds on the number of iterations required for the $\ell^1$ distance $\|r^{(t)}- p\|_1 \leq \varepsilon$ for a given $\varepsilon>0$. Although these bounds imply a convergence that is sub-linear, i.e., $\|r^{(t)}- p\|_1 = \mathcal{O}(1/t)$, their focus is on optimal dependence on problem size and dimension. An important class of problems with $A>0$ is entropy regularized optimal transport cuturi2013sinkhorn, where $A$ is of the form $A=\exp(-c/\gamma)$ with a finite cost function (matrix) $c$, i.e., $A$ strictly positive everywhere. In this setting, convergence of Sinkhorn's algorithm in discrete and continuous problems has been studied by altschuler2017near,marino2020optimal,leger2021gradient,ghosal2022convergence, among others. The linear convergence of Sinkhorn's algorithm for $A>0$ has also been extended to the multi-marginal continuous setting by carlier2022linear, building on the work of marino2020optimal.

However, the matter of convergence is more delicate when the matrix contains zero entries, and additional assumptions on the problem structure are required to guarantee the existence of scalings $D^0,D^1$ and the convergence of Sinkhorn's algorithm. For non-negative $A$, convergence is first established by sinkhorn1967concerning in the special case of {square} $A\geq 0$ and uniform $p=q=\mathbf{1}_n=\mathbf{1}_m$. Their necessary and sufficient condition is that $A$ has support, i.e., there exists a permutation $\sigma$ such that the “diagonal” $(A_{1\sigma(1)},A_{2\sigma(2)},\dots,A_{n\sigma(n)})$ is strictly positive. soules1991rate and achilles1993implications further show that the convergence is linear if and only if the stronger condition of total support holds, i.e., any non-zero entry of $A$ must be in $(A_{1\sigma(1)},A_{2\sigma(2)},\dots,A_{n\sigma(n)})$ for some permutation $\sigma$. knight2008sinkhorn provides a tight asymptotic linear convergence rate in terms of the sub-dominant (second largest) singular value of the scaled doubly stochastic matrix $D^0AD^1$. The convergence in knight2008sinkhorn is measured by some implicit distance to the target marginals $p,q$. However, no asymptotic linear convergence rate is previously known for non-square $A\geq0$ and non-uniform marginals.

For general non-negative matrices and non-uniform marginals, the necessary and sufficient conditions on $A$ in the matrix balancing problem that generalize that of sinkhorn1967concerning have been studied by thionet1964note,bacharach1965estimating,brualdi1968convex,menon1968matrix,djokovic1970note,sinkhorn1974diagonal,balakrishnan2004polynomial,pukelsheim2009iterative, among others, and convergence of Sinkhorn's algorithm under these conditions is well-known. Connecting Sinkhorn's algorithm to dual coordinate descent for entropy optimization, luo1992convergence show that the dual optimality gap, defined in (ref), converges linearly globally with some unknown rate $\lambda$ when finite scalings $D^0,D^1$ exist. However, their result is implicit and there are no results that quantify the global linear rate $\lambda$, even for special classes of non-negative matrices. When convergence results for $A>0$ in previous works are applied to non-negative matrices, the bounds often blow up or become degenerate as soon as $\min_{ij}A_{ij} \downarrow 0$. For example, in (ref) the contraction factor $\lambda \rightarrow 1$ when $A$ contains zero entries. {When $\min_{ij}A_{ij} = 0$, complexity bounds on $\|r^{(t)}- p\|_1$ and $\|r^{(t)}- p\|_2$ have been established for Sinkhorn's algorithm, for example by the works of kalantari2008complexity and chakrabarty2021better, with polynomial dependence on $1/\varepsilon$, i.e., sub-linear convergence.} Under the minimal condition that guarantees the convergence of Sinkhorn's algorithm, leger2021gradient gives a {quantitative} global sub-linear bound on the KL divergence between $r^{(t)}$ and $p$ in the continuous setting for general probability distributions, which include non-negative matrices $A\geq 0$.

It therefore remains to reconcile the various results on Sinkhorn's algorithm for $A\geq 0$ and characterize the global and asymptotic linear convergence rates for non-negative $A$. Our results precisely fill these gaps left by previous works. The global linear convergence result in (ref) establishes a contraction like (ref) for the optimality gap whenever finite scalings $D^0,D^1$ exist, and characterize the convergence rate $\lambda$ in terms of the algebraic connectivity. Moreover, the asymptotic linear rate in (ref) directly extends the result of knight2008sinkhorn. See (ref) for a detailed summary and comparison of the convergence results in previous works and this paper. The dependence of Sinkhorn's convergence rate on spectral properties of graphs can be compared to convergence results in the literature on decentralized optimization and gossip algorithms, where a spectral gap quantifies the convergence rate boyd2006randomized,xiao2007distributed.

Lastly, we note that other algorithms with better complexities have been developed for the matrix balancing problem, utilizing e.g., the ellipsoid algorithm kalantari1996complexity and geometric programming nemirovski1999complexity, interior point algorithms cohen2017matrix,chen2022maximum, or customized first/second order techniques linial1998deterministic,allen2017much. However, despite having better theoretical complexities, most of these algorithms are yet to be implemented practically. Sinkhorn's algorithm, on the other hand, remains an attractive choice in practice due to its simplicity, robustness, and parallelization capabilities.

Preliminaries on Choice Modeling and Matrix Balancing

We start by providing brief but self-contained introductions to the two main subjects of this paper, choice modeling and matrix balancing, including their respective underlying mathematical problems and assumptions. Then, we formally establish their equivalence in (ref).

Maximum Likelihood Estimation of Luce Choice Models

In the Luce choice modeling framework, we have $n$ observations $\{(j_i,S_i)\}_{i=1}^n$, each consisting of a choice set $S_{i}\subseteq \{1,\dots,m\}=[m]$ that is a subset of the total $m$ alternatives/items/objects, and the alternative selected, denoted by $j_i \in S_i$. The choice probability is prescribed by Luce's axiom of choice given model parameter $s\in\mathbb{R}^m_{++}$ in the interior of the probability simplex $\Delta_m$:

and the likelihood of the observed data is thus given by

A popular method to estimate $s = \{s_{1},\dots,s_m\}$ is the maximum likelihood estimation approach, which maximizes the log-likelihood

over the interior of the probability simplex. Note that the choice sets $S_i$ can vary across $i$. In other words, in each observation the choice is made from a potentially distinct set of alternatives. This feature of the problem turns out to be important for both the algorithmic efficiency of computing the maximizers to (ref), as well as the statistical efficiency of the resulting maximum likelihood estimators (MLEs), which can be quantified by a measure of connectivity of the data structure. We will elaborate on these points shortly. For now, we focus on the existence and uniqueness of MLE.

If we reparameterize $\exp(u_{j})=s_{j}$, it is obvious that (ref) is concave in $u$. However, to ensure the log-likelihood (ref) has a unique maximizer in the interior of the simplex, additional assumptions on the comparison structure of the dataset $\{(j_i,S_i)\}_{i=1}^n$ are needed. The following classic condition is necessary and sufficient for the maximum likelihood problem to be well-posed.

(ref) is standard in the literature hunter2004mm,noothigattu2020axioms and appeared as early as the work of zermelo1929berechnung and ford1957solution for pairwise comparisons. hunter2004mm shows that (ref) is necessary and sufficient for the upper-compactness of (ref), which guarantees the existence of a maximizer in the interior of the probability simplex. In fact, when an interior maximizer exists, it is also unique, since (ref) implies the following weaker condition, which guarantees the strict concavity of (ref).

The intuitions provided by ford1957solution and hunter2004mm are helpful for understanding Assumptions (ref) and (ref). If items from some $S\subsetneq [m]$ are never compared with those in $S^C$, i.e., never appeared together in any choice set $S_i$, it is impossible to rank across the two subsets. In this case, we can rescale the relative weights of $S$ and $S^C$ of an interior maximizer and obtain another maximizer. On the other hand, if items in $S$ are always preferred to those in $S^C$, we can increase the likelihood by scaling $s_j$ for items $j\in S^C$ towards $0$, and no maximizer in the interior of the probability simplex exists. Nevertheless, a boundary solution can still exist. This case turns out to be important in the present work: in the equivalent matrix balancing problem, it corresponds to the slowdown regime of Sinkhorn's algorithm, where scalings diverge but the scaled matrix converges ((ref)).

(ref) also has a concise graph-theoretic interpretation. Define the weighted {undirected comparison graph} $G_c$ on $m$ vertices with adjacency matrix $A^c$ given by

In other words, there is an undirected edge between $j$ and $k$ if and only if they are both included in some choice set $S_i$, with the edge weight equal to the number of their co-occurrences, which could be zero. We can verify that (ref) precisely requires $G_c$ to be connected.

The Canonical Matrix Balancing Problem

Matrix balancing is a classic problem that shows up in a wide range of disciplines. See (ref) for a short survey on some applications. The underlying mathematical problem can be stated concisely in matrix form as:

We henceforth refer to the above as the “canonical” matrix balancing problem. Other variants of the problem replace the row and column sums (the 1-norm) with other norms bauer1963optimally,ruiz2001scaling. Note that for any $c>0$, $(D^0/c,cD^1)$ is also a solution whenever $(D^0,D^1)$ is. A finite positive solution $(D^{0},D^{1})$ to the canonical matrix balancing problem is often called a direct scaling.

The structure of the matrix balancing problem suggests a simple iterative scheme: starting from any initial positive diagonal $D^{0}$, invert $D^{1}AD^{0}\mathbf{1}_m = p$ using $p/(AD^{0}\mathbf{1}_m)$ to update $D^{1}$. Then invert $D^{0}A^{T}D^{1}\mathbf{1}_n=q$ using $q/(A^{T}D^{1}\mathbf{1}_n)$ to compute the new estimate of $D^{0}$, and repeat the procedure, leading to a solution if it converges. Here, divisions involving two vectors of the same length are entry-wise. This simple iterative scheme is precisely Sinkhorn's algorithm, described in (ref), where vectors $d^0,d^1$ are the diagonal elements of $D^0,D^1$.

An important dichotomy occurs depending on whether the entries of $A$ are strictly positive. If $A$ contains no zero entries, then direct scalings and a unique scaled matrix $D^1AD^0$ always exist sinkhorn1964relationship. Moreover, Sinkhorn's algorithm converges linearly franklin1989scaling. On the other hand, when $A$ contains zero entries, the problem becomes more complicated. Additional conditions are needed to guarantee meaningful solutions, and the convergence behavior of Sinkhorn's algorithm is less clearly understood. Well-posedness of the matrix balancing problem has been studied by brualdi1968convex,sinkhorn1974diagonal,pukelsheim2009iterative, among others, who characterize the following equivalent existence conditions.

It is well-known in the matrix balancing literature that the above two conditions are equivalent, and that a positive finite solution $(D^{0},D^{1})$ to the canonical problem exists if and only if they hold. See, for example, Theorem 6 in pukelsheim2009iterative. (ref) also guarantees the convergence of Sinkhorn's algorithm. However, it is not a necessary condition. In other words, Sinkhorn's algorithm could converge even if the matrix balancing problem does not admit a direct scaling. This phenomenon turns out to be important in characterizing the convergence rate, which we study in (ref).

Clearly, (ref)(a) is the minimal necessary condition when a solution to the matrix balancing problem exists, and trivially holds when $A>0$ (take for example $A'$ as the Kronecker product of $p,q$). (ref)(b) is closely connected to conditions for perfect matchings in bipartite graphs hall1935representatives,galichon2021matching. In flow networks gale1957theorem,ford1956maximal,ford1957simple, it is a capacity constraint that guarantees the maximum flow on a weighted bipartite graph is equal to $\sum_i p_i=\sum_j q_j$ and with positive flow on every edge idel2016review. The weighted bipartite graph, denoted by $G_b$, is important in this paper. Its adjacency matrix $A^b \in \mathbb{R}^{(n+m)\times(n+m)}$ can be represented concisely using $A$ as

and $A$ is sometimes called the biadjacency matrix of $G_b$. See (ref) for more information. Just like in the choice setting, where the connectivity of the undirected comparison graph $G_c$ plays an important role, the connectivity of $G_b$ turns out to be important for the linear convergence rate of Sinkhorn's algorithm (see (ref)).

Lastly, the necessary and sufficient condition for the uniqueness of finite scalings essentially requires that $A$ is not block-diagonal, and precisely guarantees that $G_b$ is connected.

With a proper introduction to both problems, we are now ready to establish the equivalence between Luce choice model estimation and matrix balancing. In (ref), we return to Sinkhorn's algorithm for the matrix balancing problem, and provide answers to open problems concerning its linear convergence for non-negative $A$, by leveraging the connections we establish next.

Connecting Choice Modeling and Matrix Balancing

In this section, we formally establish the connections between choice modeling and matrix balancing. We show that maximizing the log-likelihood (ref) is equivalent to solving a canonical matrix balancing problem. We also precisely describe the correspondence between the relevant conditions in the two problems. In view of this equivalence, we show that Sinkhorn's algorithm, when applied to estimate Luce choice models, is in fact a parallelized generalization of the classic iterative algorithm for choice models, dating back to zermelo1929berechnung,dykstra1956note,ford1957solution, and studied extensively also by hunter2004mm,vojnovic2020convergence,vojnovic2023accelerated.

Maximum Likelihood Estimation of Luce Choice Models as Matrix Balancing

The optimality conditions for maximizing the log-likelihood (ref) for each $s_j$ are given by

Multiplying by $s_{j}$ and dividing by $1/n$, we have

where $W_j:=|\{i\in[n]\mid (j,S_i)\}|$ is the number of observations where $j$ is selected.

Note that in the special case where $S_i\equiv [n]$, i.e., every choice set contains all items, the MLE simply reduces to the familiar empirical frequencies $\hat s_j = {W_j}/{n}$. However, when the choice sets $S_i$ vary, no closed form solution to (ref) exists, which is the primary motivation behind the long line of works on the algorithmic problem of solving (ref). With varying $S_i$, we can interpret the optimality condition as requiring the observed frequency of $j$ being chosen (left hand side) be equal to the average or expected probability of $j$ being selected (right hand side), which conditional on choice set $S_i$ is $\frac{s_{j}}{\sum_{k\in S_{i}}s_{k}}$ if $j\in S_i$ and 0 otherwise. In addition, note that since the optimality condition in (ref) only involves the frequency of selection, distinct datasets could yield the same optimality conditions and hence the same MLEs. For example, suppose that two alternatives $j$ and $k$ both appear in choice sets $S_{i}$ and $S_{i'}$, with $j$ selected in $S_{i}$ and $k$ selected in $S_{i'}$. Then switching the choices in $S_{i}$ and $S_{i'}$ does not alter the likelihood and optimality conditions. This feature holds more generally with longer cycles of items and choice sets, and can be viewed as a consequence of the context-independent nature of Luce's choice axiom, i.e., IIA. In some sense, it is also the underpinning of many works in economics that estimate choice models based on marginal sufficient statistics. A prominent example is berry1995automobile, which estimates consumer preferences using data on {aggregate} market shares of products.

We are now ready to reformulate (ref) as a canonical matrix balancing problem. Define $p\in\mathbb{R}^{n}$ as $p_i=R_i$, i.e., the number of times choice set $S_i$ appears in the data. Define $q\in\mathbb{R}^{m}$ as $q_j={W_j}$, i.e., the number of times item $j$ was selected in the data. By construction we have $\sum_i p_i=\sum_j q_j$, and $p_i,q_j>0$ whenever (ref) holds.

Now define the $n\times m$ binary matrix $A$ by $A_{ij}=1\{j\in S_{i}\}$, so the $i$-th row of $A$ is the indicator of which items appear in the (unique) choice set $S_i$, and the $j$-th column of $A$ is the indicator of which choice sets item $j$ appears in. We refer to this $A$ constructed from a choice dataset as the participation matrix. By construction, $A$ has distinct rows, but may still have identical columns. We may also remove repeated columns by “merging” items and their win counts. Their estimated scores can be computed from the score of the merged item proportional to their respective win counts. We do not require this reduction in our results. (ref) provides an illustration of the matrix balancing representation of the luce choice modeling problem with $(A,p,q)$ defined as above.

Let $D^{0}\in\mathbb{R}^{m\times m}$ be the diagonal matrix with $D_{j}^{0}=s_{j}$ and $D^{1}\in\mathbb{R}^{n\times n}$ be the diagonal matrix with $D_{i}^{1}={R_i}/{\sum_{k\in S_{i}}s_{k}}$, and define the scaled matrix

The matrices $D^{1}$ and $D^{0}$ are scalings of rows and columns of $A$, respectively, and

The key observation is that the optimality conditions (ref) can be rewritten as

Moreover, by construction $\hat{A}$ also satisfies

Therefore, if $s_j$'s satisfy the optimality conditions for maximizing (ref), then $D^0,D^1$ defined above solve the matrix balancing problem in (ref). Moreover, the converse is also true, and we thus establish the equivalence between choice maximum likelihood estimation and matrix balancing. All omitted proofs appear in (ref).

(ref) implies that (ref) has a unique maximizer $s$ in the interior of the probability simplex if and only if (ref) has a unique normalized solution $D^0$ as well. The next question, naturally, is then how (ref) and (ref) for choice modeling are connected to (ref) and (ref) for matrix balancing.

Thus when the choice maximum likelihood estimation problem is cast as a matrix balancing problem, (ref) exactly characterizes the gap between (ref) and (ref). We provide some intuition for (ref). When we construct a triplet $(A,p,q)$ from a choice dataset, with $p$ the numbers of appearances of unique choice sets and $q$ the winning counts of each item, (ref) precludes the possibility of partitioning the items into two subsets that never get compared with each other, i.e., (ref). Then (ref) requires that whenever a strict subset $M\subsetneq [m]$ of objects only appear in a strict subset $N\subsetneq [n]$ of the observations, their total winning counts are strictly smaller than the total number of these observations, i.e., there is some object $j\notin M$ that is selected in $S_i$ for some $i\in N$, which is required by (ref).

Interestingly, while (ref) requires the directed comparison graph, defined by the $m\times m$ matrix of counts of item $j$ being chosen over item $k$, to be strongly connected, the corresponding conditions for the equivalent matrix balancing problem concern the $n\times m$ participation matrix $A$ and positive vectors $p,q$, which do not explicitly encode the specific choice of each observation. This apparent discrepancy is due to the fact that $(A,p,q)$ form the sufficient statistics of the Luce choice model. In other words, there can be more than one choice dataset with the same optimality condition (ref) and $(A,p,q)$ defining the equivalent matrix balancing problem (see (ref)).

Having formulated a particular matrix balancing problem from the estimation problem given choice data, we may ask how one can go in the other direction. In other words, when/how can we construct a “choice dataset” whose sufficient statistics is a given triplet $(A,p,q)$? First off, for $(A,p,q)$ to be valid sufficient statistics of a Luce choice model, $p,q$ need to be positive integers. Moreover, $A$ must be a binary matrix with unique rows, each containing at least two non-zero elements (valid choice sets have at least two items). Given such a $(A,p,q)$ satisfying Assumptions (ref) and (ref), a choice dataset can be constructed efficiently. Such a procedure is described, for example, in kumar2015inverting, where $A$ is motivated by random walks on a graph instead of matrix balancing ((ref)). Their construction relies on finding the max flow on the bipartite graph $G_b$. For rational $p,q$, this max flow can be found efficiently in polynomial time balakrishnan2004polynomial,idel2016review. Moreover, the maximum flow implies a matrix $A'$ satisfying (ref)(a), thus providing a feasibility certificate for the matrix balancing problem as well.

We have thus closed the loop and fully established the equivalence of the maximum likelihood estimation of Luce choice models and the canonical matrix balancing problem.

Connections to discrete choice modeling have also been established for the related problem of regularized semi-discrete optimal transport tacskesen2023semi, although the problem setting and results are distinct from the ones studied in this paper. We next turn our attention to the algorithmic connections between choice modeling and matrix balancing.

Algorithmic Connections between Matrix Balancing and Choice Modeling

Given the equivalence between matrix balancing and choice modeling, we can naturally consider applying Sinkhorn's algorithm to maximize (ref). In this case, one can verify that the updates in each full iteration of (ref) reduce algebraically to

in the $t$-th iteration. Comparing (ref) to the optimality condition in (ref), which we recall is given by

we can therefore interpret the iterations as simply dividing the winning count $W_j$ by the coefficient of $s_j$ on the right repeatedly, in the hope of converging to a fixed point. A similar intuition was given by ford1957solution in the special case of pairwise comparisons. Indeed, the algorithm proposed in that paper is a cyclic variant of (ref) applied to pairwise comparisons. However, this connection is mainly algebraic, as the optimality condition in ford1957solution does not admit a reformulation as the matrix balancing problem in (ref).

In (ref), we provide further discussions on the connections of Sinkhorn's algorithm to existing frameworks and algorithms in the choice modeling literature and the optimization literature. We demonstrate that many existing algorithms for Luce choice model estimation are in fact special cases or analogs of Sinkhorn's algorithm. These connections also illustrate the many interpretations of Sinkhorn's algorithm, e.g., as a distributed optimization algorithm as well as a minorization-maximization or majorization-minimization (MM) algorithm lange2016mm. However, compared to most algorithms for choice modeling discussed in this work, Sinkhorn's algorithm is more general as it applies to non-binary $A$ and non-integer $p,q$, and has the additional advantage of being paralellized and distributed, hence more efficient in practice. To highlight the ubiquity of Sinkhorn's algorithm in the choice setting, we summarize these algorithmic connections below.

The mathematical and algorithmic connections between matrix balancing and choice modeling we establish in this paper allow the transmission of ideas in both directions. For example, inspired by regularized maximum likelihood estimation maystre2017choicerank, we propose a regularized version of Sinkhorn's algorithm in (ref), which is guaranteed to converge even when the original Sinkhorn's algorithm does not converge. Leveraging optimization connections between maximum likelihood estimation of choice models and matrix balancing, chang2024inferring propose a statistical framework for network traffic that justifies the popular use of Sinkhorn's algorithm to infer detailed dynamic networks from aggregate node-level activities. In the rest of the main text, we focus on resolving some interesting open problems on the convergence of Sinkhorn's algorithm, motivated by results in choice modeling on the importance of algebraic connectivity in quantifying statistical and computational efficiencies.

Linear Convergence of Sinkhorn's Algorithm for Non-negative Matrices

In this section, we turn our attention on matrix balancing and study the global and asymptotic linear convergence of Sinkhorn's algorithm for general non-negative matrices $A\geq 0$ and positive marginals $p,q>0$. We first present the relevant optimization principles behind matrix balancing and discuss some existing results in (ref), and then present the convergence results in Sections (ref)--(ref). Throughout, we use superscript $(t)$ to denote quantities after $t$ iterations of Sinkhorn's algorithm with normalized columns, described in (ref).

Preliminaries

We start with the optimization principles associated with matrix balancing and Sinkhorn's algorithm. Given a matrix balancing problem with $A\geq 0,\sum_{ij}A_{ij}=1$ and target marginals $p,q$ with $\sum_ip_i=\sum_jq_j=1$, consider the following KL divergence (relative entropy) minimization problem

It is well-known that when scalings $D^0,D^1$ solve the matrix balancing problem with $(A,p,q)$, the scaled matrix $\hat A=D^1AD^0$ is the unique minimizer of (ref) bregman1967proof,ireland1968contingency. Moreover, vector representations $d^0,d^1$ of the optimal scalings $D^0,D^1$ precisely minimize the following (negative) dual objective of (ref):

and Sinkhorn's algorithm is a block coordinate descent (BCD) type algorithm tseng2001convergence applied to minimize (ref). In luo1992convergence, the authors study the linear convergence of block coordinate descent algorithms for a general class of objectives that includes (ref). In particular, their result implies that the convergence of Sinkhorn's algorithm, measured in terms of the optimality gap of $g$, is linear with some unknown rate, as long as finite positive scalings $D^0,D^1$ exist that satisfy (ref). The function $g$, sometimes referred to as the potential function of Sinkhorn's algorithm, also turns out to be crucial in quantifying the global linear convergence rate in the present work.

Although convergence results on Sinkhorn's algorithm are abundant, they are often stated in different forms, developed under different assumptions, and apply to settings of varying degrees of generality. Next, we briefly discuss some existing works to clarify their connections and distinctions, which also helps motivate the technical results in this paper. They are summarized in (ref).

First and foremost, how to define and measure the convergence of Sinkhorn's algorithm is not entirely trivial. Because of the indeterminacy of scalings under the transformation $(D^0,D^1) \rightarrow (D^0/c,c\cdot D^1)$, most works define convergence using quantities that are invariant under this transformation. Let $D^{0(t)},D^{1(t)}$ be the scalings obtained after $t$ iterations of Sinkhorn's algorithm, and let ${A}^{(t)}:=D^{1(t)}AD^{0(t)}$ be the scaled matrix based on these scalings. Since ${A}^{(t)}$ is invariant under the transformation, some earlier works such as franklin1989scaling,soules1991rate measure convergence in terms of ${A}^{(t)}$ to the optimally scaled matrix $\hat A= D^{1}AD^{0}$ that has the target row and column sums $p,q$. Most later works focus instead on the convergence of the marginal quantities

to the target row and column sums $p,q$. Since after each iteration in (ref), the column constraint is always satisfied: ${A}^{(t)}\mathbf{1}_n=q$, it suffices to focus on the convergence of $r^{(t)}$ to $p$. For example, leger2021gradient uses the KL divergence $D_{\text{KL}}(r^{(t)}\| p)$, while altschuler2017near,chakrabarty2021better use the $\ell^1$ distance $\|r^{(t)}- p\|_1$, which is upper bounded by the KL divergence via Pinsker's inequality. Given the entropy optimization perspective of matrix balancing and Sinkhorn's algorithm in (ref) and (ref), it is also possible to measure convergence in terms of the dual optimality gap $g(d^{0(t)},d^{1(t)})-g(d^{0},d^{1})$, which in turn bounds KL divergences via

luo1992convergence show that this dual optimality gap converges linearly, but how the convergence rate depends on problem structure has been an open question. Our global linear convergence result in (ref) is the first to characterize this rate.

Next, we note the distinction between global and local (particularly asymptotic) convergence results. Global results hold for all iterations $t>0$. For example, leger2021gradient shows $D_{\text{KL}}(r^{(t)}\| p)\leq D_{\text{KL}}^\ast /t$ for all $t$, where $D_{\text{KL}}^\ast$ is the optimal value of the relative entropy minimization problem (ref). On the other hand, local convergence results pertain to the behavior of an algorithm in a neighborhood of the optimal solution, while asymptotic results only hold in the limit as $t\rightarrow \infty$. For example, knight2008sinkhorn characterizes the asymptotic rate $\lim_{t\rightarrow \infty}\|r^{(t+1)}-p\|_{\ast}/\|r^{(t)}-p\|_{\ast}$ of linear convergence for Sinkhorn's algorithm, where $\|\cdot\|_\ast$ is an implicitly defined norm. When $A\geq 0$, knight2008sinkhorn is the only work on the exact asymptotic convergence rate, for square $A$ and uniform $p,q$. Our (ref) is the first asymptotic result for general $A\geq0$ and non-uniform $p,q$, with an explicit norm ($\|\cdot\|_2$). While asymptotic results provide a more precise description of an algorithm's behavior near the optimal solution, global results are useful for obtaining complexity bounds on the number of iterations required to obtain $\varepsilon$ accuracy solutions. In fact, global results are often stated directly as complexity bounds. For example, the result in altschuler2017near is that for $A>0$, $t=4(\varepsilon)^{-2}\log(\sum_{ij}A_{ij}/\min_{ij}A_{ij})$ iterations of Sinkhorn's algorithm guarantee $\|r^{(t)}- p\|_1\leq \varepsilon$.

Lastly, we note the distinction between (global) linear and sub-linear convergence results. Linear convergence is often understood as successive improvements of the convergence metric by a constant factor. For example, franklin1989scaling show that for $A>0$ and the Hilbert metric $d$, $d(r^{(t+1)}, p) \leq \lambda \cdot d(r^{(t)}, p)$ for all $t>0$ for some $\lambda\in (0,1)$. As a result, $d(r^{(t)}, p) \leq \lambda^t \cdot d(r^{(0)}, p)$ decreases exponentially in $t$, so that $d(r^{(t)}, p)\leq \varepsilon$ in $\mathcal{O}(\log(1/\varepsilon))$ iterations. In contrast, in sub-linear results, such as leger2021gradient, the convergence metric $D_{\text{KL}}(r^{(t)}\| p)$ only decreases polynomially in $t$, requiring $\mathcal{O}(1/\varepsilon)$ iterations to guarantee $D_{\text{KL}}(r^{(t)}\| p) \leq \varepsilon$. While sub-linear complexity bounds have worse (polynomial) dependence on $1/\varepsilon$, they often focus on optimizing the dependence on problem size and dimension. Our main focus in this paper is on understanding the linear convergence behavior of Sinkhorn's algorithm when $A\geq0$, i.e., $\mathcal{O}(\log(1/\varepsilon))$ iteration complexity. Nevertheless, we also provide refined complexity bounds in (ref) that optimize dependence on problem constants.

As discussed before, when $A\geq 0$, only sub-linear convergence results with explicit rates are known kalantari2008complexity,chakrabarty2021better,leger2021gradient, while luo1992convergence implies global linear convergence with an unknown rate. We now characterize this global linear rate of convergence in terms of the algebraic connectivity of the bipartite graph defined in (ref).

Global Linear Convergence

Our analysis starts with the following change of variables to transform the potential function (ref):

resulting in the reparameterized potential function $g(u,v)$ of (ref):

Note first that $g(u,v)=g(u+a,v+a)$ for any constant $a\in \mathbb{R}$. We can verify that Sinkhorn's algorithm is equivalent to the alternating minimization algorithm bertsekas1997nonlinear for (ref), which alternates between minimizing with respect to $u$ and $v$, holding the other block fixed:

or written more explicitly element-wise,

A main reason to focus on (ref) instead of the log-barrier form (ref) is that (ref) has a Hessian with desirable properties for proving linear convergence. The Hessian of $g(u,v)$ is

where $\mathcal{D}$ converts a vector to a diagonal matrix, and $\hat{A}=\mathcal{D}(d^1)A\mathcal{D}(d^0)=\mathcal{D}(\exp(-v))A\mathcal{D}(\exp(u))$ is the matrix scaled by $u,v$. Note that the Hessian $\nabla^2g(u,v)$ always has $\mathbf{1}_{m+n}$ in its null space. On the surface, it may seem that standard linear convergence results for first-order methods, which require strong convexity (or the Polyak--{\L}ojasiewicz condition) of the objective function, do not apply to $g(u,v)$. However, we will show that whenever the matrix balancing problem has finite scaling solutions, $g(u,v)$ is in fact strongly convex when restricted to bounded subsets of the subspace

Moreover, the invariance of $g(u,v)$ and its gradient and Hessian under constant translations of $(u,v)$ by $\mathbf{1}_{m+n}$ guarantees that the strong convexity constant of $g(u,v)$ on $\mathbf{1}_{m+n}^\perp$ in fact quantifies the linear convergence of Sinkhorn's algorithm even if the iterates $u^{(t)},v^{(t)}$ are not in $\mathbf{1}_{m+n}$. Similar types of “restricted strong convexity” properties have been studied by for example agarwal2010fast. It also shares similarities with the exp-concavity property popular in online learning hazan2016introduction,orabona2019modern, which implies the strong convexity of a function in the direction of the gradient, evaluated at any point. In contrast, $g(u,v)$ is strongly convex along any direction orthogonal to $\mathbf{1}_{m+n}$, but its gradient evaluated at any $(u,v)$ is not necessarily orthogonal to $\mathbf{1}_{m+n}$. However, the key is that along the trajectory of iterates $(u^{(t)},v^{(t)})$ obtained by running Sinkhorn's algorithm, the gradients of $g$ evaluated at $(u^{(t)},v^{(t)})$ are indeed orthogonal to $\mathbf{1}_{m+n}$, which is sufficient to guarantee the linear convergence of Sinkhorn's algorithm.

We now introduce the key quantities and definitions used in our result. Let Sinkhorn's algorithm initialize with a $u^{(0)}$, with $v^{(0)}$ given by (ref). Define the constant $B$ as

In other words, $B$ is the diameter of the initial normalized sub-level set. We will show that $B$ is finite and that it bounds normalized Sinkhorn iterates $\|(u^{(t)},v^{(t)})\|_\infty$, since under (ref) the function $g(u,v)$ is coercive on the subspace $\mathbf{1}^\perp_{m+n}$. Coercivity is an important property, and we define it below, following bertsekas2016nonlinear.

Next, define the Laplcian matrix $\mathcal{L}:=\mathcal{L}(A)$ of the bipartite graph $G_b$ (see (ref)) as

and refer to the second smallest eigenvalue $\lambda_{-2}(\mathcal{L})$ as the Fiedler eigenvalue. $\lambda_{-2}(\mathcal{L})>0$ if and only if (ref) holds, and it quantifies “connectivity” of the data structure spielman2007spectral. Although an important quantity in the choice modeling literature, algebraic connectivity has not been previously used in the analysis of Sinkhorn's algorithm. For details on the graph Laplacian and the Fiedler eigenvalue, see (ref). Finally, define the smoothness parameters

which are used to quantify the smoothness of $g(u,v)$.

We can now state one of our main contributions to the study of Sinkhorn's algorithm.

The linear convergence rate of Sinkhorn's algorithm is therefore quantified by $\lambda_{-2}(\mathcal{L})/\min\{l_{0},l_{1}\}$, which is invariant under rescalings of $A\rightarrow c\cdot A$. Although the corresponding bipartite graph $G_b$ with biadjacency matrix $A$ is a natural object to consider in the study matrix balancing problems, to our knowledge, (ref) is the first to highlight the precise role of its spectral property, described by $\lambda_{-2}(\mathcal{L})$, in the linear convergence of Sinkhorn's algorithm. It fills the gap left by luo1992convergence who establish linear convergence with an implicit rate, and allows us to compute its dependence on problem parameters by applying established bounds on $\lambda_{-2}(\mathcal{L})$ from spectral graph theory spielman2007spectral.

With this auxiliary normalization procedure, the proof of (ref) then relies on the observation that the Hessian of $g(u,v)$ is precisely the graph {Laplacian} $\mathcal{L}(u,v)$ of the bipartite graph with biadjacency matrix $\hat A = \mathcal{D}(\exp(-v))A\mathcal{D}(\exp(u))$. As $(u,v)$ are {bounded} on normalized Sinkhorn iterates thanks to the coercivity of $g$, the Fiedler eigenvalue of $\mathcal{L}=\mathcal{L}(0,0)$ quantifies the strong convexity on $\mathbf{1}_{m+n}^\perp$. Linear convergence then follows from results on block coordinate descent and alternating minimization methods for strongly convex and smooth functions beck2013convergence. Typically, the leading eigenvalue of the Hessian quantifies the smoothness, which is bounded by $2\max\{l_0,l_1\}$ for $\mathcal{L}$. For alternating minimization methods, the better smoothness constant $\min\{l_0,l_1\}$ is available. {Thus the quantity $\min\{l_{0},l_{1}\}/\lambda_{-2}(\mathcal{L})$ in (ref) can be interpreted as a type of “condition number” of the graph Laplacian $\mathcal{L}$.} When $A$ is positive (not just non-negative), then the strong existence and uniqueness conditions are trivially satisfied, and our results continue to hold with the rate quantified by $\min\{l_{0},l_{1}\}/\lambda_{-2}(\mathcal{L})$. In this case, both $\min \{l_0,l_1\}$ and $\lambda_{-2}(\mathcal{L})$ are $\Theta(n)$ where $n$ is problem dimension, so $\min\{l_{0},l_{1}\}/\lambda_{-2}(\mathcal{L})$ does not increase with problem dimension.

Although (ref) implies an $\mathcal{O}(\log(1/\varepsilon))$ iteration complexity, the complexity bound's dependence on problem parameters can be further improved. In particular, the constant $B$ which bounds $\|(u^{(t)},v^{(t)})\|_\infty$ can be hard to compute for some problems. We next establish an iteration complexity bound that does not depend exponentially on the implicit constant $B$ defined in (ref).

(ref) relies on two technical innovations. First, we bound $e^B \leq C$ where $C$ is explicitly constructed from any optimal solution pair and is invariant under rescalings. Second, we improve the dependence from $e^{4B}$ to $e^{2B}$, using the target marginals $p,q$ to quantify the smoothness of $g(u,v)$ instead. Our message here is that the convergence behavior of Sinkhorn's algorithm has two phases. Initially, we can apply a sub-linear complexity bound with $\mathcal{O}(1)$ iterations to obtain Sinkhorn iterates sufficiently close to the optimal solution. Afterwards, the convergence can better be captured by a linear convergence with rate depending on the optimal solution and target marginals $p,q$. The dependence of $C$ on problem dimension is problem specific. In the worst case, it can be exponential kalantari1993rate. In Appendix (ref), we plot (ref) as a function of problem dimension on randomly generated data, and find the dependence is quadratic. In contrast, sub-linear bounds, such as altschuler2017near,chakrabarty2021better, have logarithmic dependence on problem dimension. It remains an interesting question to improve the dependence on problem dimension in (ref), and to study trade-offs with the dependence on $\varepsilon$.

Strong vs. Weak Convergence of Sinkhorn's Algorithm

We now discuss the two different convergence regimes of Sinkhorn's algorithm when $A\geq 0$. As mentioned in (ref), when $A\geq 0$, the canonical matrix balancing problem with target marginals $p,q$ has a finite positive solution pair $D^0,D^1$ if and only if (ref) holds (which trivially holds when $A>0$). In this case, Sinkhorn's algorithm converges to $D^1AD^0$, which also solves the KL minimization problem (ref). We call this case the strong convergence of Sinkhorn.

However, even if (ref) fails and no positive finite scalings $D^0,D^1$ exist that solve the matrix balancing problem, the sequence of scaled matrices ${A}^{(t)}=D^{1(t)}AD^{0(t)}$ based on Sinkhorn's algorithm can still converge entry-wise to the solution of (ref), whenever it has a finite solution. This apparent discrepancy is explained by the fact that the solution of problem (ref) requires a weaker condition than (ref) for the matrix balancing problem. It can be stated in the following equivalent forms.

The equivalence of the two conditions above follows from Theorem 4 in pukelsheim2009iterative, which also shows that they are the minimal requirements for the convergence of Sinkhorn's algorithm. (ref)(a) precisely guarantees that the constrained KL minimization problem (ref) is feasible and bounded. It relaxes (ref)(a) by allowing additional zeros in the matrix $A'$. Similarly, (ref)(b) relaxes (ref)(b) by allowing equality between $\sum_{i\in N}p_i$ and $\sum_{j\in M}q_j$ even when $M,N$ do not correspond to a block-diagonal structure.

The distinction between (ref) and (ref) is important for the matrix balancing problem and Sinkhorn's algorithm. (ref) guarantees the solutions of (ref) and (ref) coincide, and have exactly the same zero pattern as $A$. If (ref) holds but (ref) fails, then the solution $\hat A$ of (ref) has additional zeros relative to $A$, and no direct (finite and positive) scaling $(D^{0},D^{1})$ exists such that $\hat A=D^1AD^0$. However, the sequence of scaled matrices $\hat A^{(t)}$ still converges to $\hat A$. We call this case the weak convergence of Sinkhorn. In this case, the matrix balancing problem is said to have a limit scaling, where some entries of $D^{0},D^{1}$ in Sinkhorn iterations approach 0 or $\infty$, resulting in additional zeros in $\hat A$. Below we give an example adapted from pukelsheim2009iterative, where $p,q=(3,3)$ and the scaled matrices $\hat A^{(t)}$ converge but no direct scaling exists:

Under (ref), leger2021gradientproves the sub-linear convergence of Sinkhorn's algorithm, while it is known since at least soules1991rate,luo1992convergence that the convergence is linear under (ref). It is therefore important to clarify the convergence behaviors of Sinkhorn's algorithm in the two settings. We next show that if (ref) holds but (ref) fails, then there exists an entry of $A^{(t)}$ that converges at a lower bound rate $\Omega(1/t)$, i.e., sub-linear. Together with existing and new results in this paper, (ref) fully characterizes the following convergence behavior of Sinkhorn's algorithm: whenever a direct scaling exists for the matrix balancing problem, Sinkhorn's algorithm converges linearly. If only a limit scaling exists, then convergence deteriorates to sub-linear. This generalizes the observations made by sinkhorn1967concerning,achilles1993implications for square matrices and uniform marginals.

The regime of sub-linear convergence also has an interpretation in the choice modeling framework. The weak existence condition (ref), when applied to $(A,p,q)$ constructed from a choice dataset, allows the case where some subset $S$ of items is always preferred over $S^C$, which implies, as observed already by the early work of ford1957solution, that the log-likelihood function (ref) is only maximized at the boundary of the probability simplex, by shrinking $s_j$ for $j\in S^C$ towards 0, i.e., $D^0_j \rightarrow 0$. Incidentally, bacharach1965estimating also refers to the corresponding regime in matrix balancing as “boundary solutions”.

Sharp Asymptotic Rate

Having established the global convergence and iteration complexity of Sinkhorn's algorithm when finite scalings exist, we now turn to the problem of characterizing the sharp, i.e., best possible asymptotic linear convergence rate as $t\rightarrow \infty$, for general non-negative $A$ and non-uniform marginals $(p,q)$. knight2008sinkhorn computed this rate for uniform $(p,q)$ under an implicit metric. Our analysis is distinct from theirs and relies on an intrinsic orthogonality structure of Sinkhorn's algorithm, which is also different from the auxiliary normalization in our global linear convergence analysis. Note that unlike the global rate, which depends on the initial problem data $A$ and $(u^{(0)},v^{(0)})$, the asymptotic rate now depends on the optimal solution $\hat A=D^1AD^0$, as expected.

In the special case of $m=n$ and $p=q=\mathbf{1}$, the asymptotic rate in (ref) reduces to that in knight2008sinkhorn. Note, however, that the convergence metric is different: we use the $\ell^2$ norm $\|r^{(t)}/\sqrt{p}-\sqrt{p}\|_2$ while knight2008sinkhorn uses $\|r^{(t)}-p\|_\ast$ with an {implicit} norm $\|\cdot\|_\ast$ on $\mathbb{R}^n$. Moreover, one cannot directly extend results for square matrices such as theirs to non-square matrices by padding them with zeros, as doing so results in target marginals that are not strictly positive. See, however, knight2008sinkhorn for a symmetrization proposal.

The proof of (ref) relies on a sequence of novel data-dependent mappings associated with Sinkhorn's algorithm. Intuitively, the dependence of the asymptotic linear rate on the second largest eigenvalue of $\tilde{A}^T\tilde{A}$ (and $\tilde{A}\tilde{A}^T$) is due to the fact that near the fixed point $\sqrt{p}$ of the mapping associated with Sinkhorn iterations, $\tilde{A}\tilde{A}^T$ (which is the Jacobian at $\sqrt{p}$) approximates the first order change in $r^{(t)}/\sqrt{p}$. Normally, the leading eigenvalue quantifies this change. The unique leading eigenvalue of $\tilde{A}\tilde{A}^T$ is equal to 1 with eigenvector $\sqrt{p}$, which does not imply contraction. Fortunately, using the quantity $r^{(t)}/\sqrt{p}$ allows us to exploit the following orthogonality structure:

by virtue of Sinkhorn's algorithm preserving the quantities $r^{(t)T}\mathbf{1}_{n}$ for all $t$. Thus, the residual $r^{(t)}/\sqrt{p}-\sqrt{p}$ is always orthogonal to $\sqrt{p}$, which is both the leading eigenvector and the fixed point of the iteration. The convergence is then controlled by the second largest eigenvalue of $\tilde{A}\tilde{A}^T$. This proof approach echoes that of the global linear convergence result in (ref), where we also exploit an orthogonality condition to obtain a meaningful bound. In (ref) the bound depends on the second smallest eigenvalue of a Hessian matrix, while in (ref) the bound depends on the second largest eigenvalue of a Jacobian matrix.

{ Lastly, we note that the asymptotic rate $\lambda_\infty$ is itself a Fiedler eigenvalue, associated with the Laplacian that is the Schur complement of the scaled graph Laplacian

}

Conclusion

In this paper, we develop extensive connections between matrix balancing and choice modeling. We show that the maximum likelihood estimation of choice models based on the Luce axioms of choice is an instance of the canonical matrix balancing problem. Moreover, many algorithms in choice modeling can be viewed as special cases or analogs of Sinkhorn's algorithm for matrix balancing. These connections can potentially benefit multiple disciplines. For choice modeling, they open the door to tools and insights from well-studied topics in optimization and numerical linear algebra. For matrix balancing, the connections enable us to resolve some interesting open problems on the linear convergence of Sinkhorn's algorithm for non-negative matrices. We establish the first quantitative global linear convergence result for Sinkhorn's algorithm applied to general non-negative matrices. Our analysis reveals the importance of algebraic connectivity for matrix balancing. We also provide the first characterization of the exact asymptotic linear rate of convergence for general non-negative matrix and non-uniform target marginals. Lastly, we clarify the linear and sub-linear convergence behaviors of Sinkhorn's algorithm under the strong and weak existence assumptions for matrix balancing. Overall, we believe that the connections established in this paper are useful for researchers from different domains and can lead to further interesting results.

\ACKNOWLEDGMENT{This work is supported in part by NSF CAREER Award \#2143176 and European Research Council consolidator grant (ERC-CoG) No. 866274. We are very grateful for insightful comments and suggestions from Kurt Anstreicher, Dimitris Bertsimas, St\'ephane Bonhomme, Agostino Capponi, Serina Chang, Patrick Ding, Yanqin Fan, Robert Freund, Han Hong, Yuchen Hu, Guido Imbens, S{\"u}leyman Kerimov, Yongchan Kwon, Frederic Koehler, Flavien L\'eger, Gary Qian, Pavel Shibayev, Ruoxuan Xiong, and Yinyu Ye. We also thank the area editor, the associate editor, and the three anonymous referees for helping to improve the manuscript significantly.}