A Multi-Objective Interpretation of Optimal Transport

TL;DR

This paper links discrete optimal transport to multi-objective optimization, introducing balanced solutions that correspond to optimal transport solutions, and discusses their computation via regularization.

Contribution

It introduces the concept of balanced solutions in multi-objective optimization and establishes their equivalence with optimal transport solutions.

Findings

Balanced solutions correspond one-to-one with optimal transport solutions.

Examples from various fields illustrate the natural occurrence of balanced solutions.

Regularization methods can be used to compute these solutions.

Abstract

This paper connects discrete optimal transport to a certain class of multi-objective optimization problems. In both settings, the decision variables can be organized into a matrix. In the multi-objective problem, the notion of Pareto efficiency is defined in terms of the objectives together with non-negativity constraints and with equality constraints that are specified in terms of column sums. A second set of equality constraints, defined in terms of row sums, is used to single out particular points in the Pareto efficient set which are referred to as "balanced solutions". Examples from several fields are shown in which this solution concept appears naturally. Balanced solutions are shown to be in one-to-one correspondence with solutions of optimal transport problems. As an example of the use of alternative interpretations, the computation of solutions via regularization is discussed.

Equations

$${\Psi }_i(X_{i\cdot })=\sum _{j=1}^m{g}_{ij}(x_{ij})(i=1,\dots ,n)$$

$$\sum _{i=1}^n{x}_{ij}=c_j(j=1,\dots ,m),x_{ij}\ge 0(i=1,\dots ,n,j=1,\dots ,m)$$

$$\sum _{j=1}^m{x}_{ij}=r_i(i=1,\dots ,n)$$

$$\sum _{i=1}^n{r}_i=\sum _{j=1}^m{c}_j.$$

$$\sum _{i=1}^n{p}_i{x}_{ij}=c_j,\sum _{j=1}^m{q}_j{x}_{ij}=r_i$$

$${\tilde{x}}_{ij}=p_i{q}_j{x}_{ij},{\tilde{c}}_j=q_j{c}_j,{\tilde{r}}_i=p_i{r}_i.$$

$${\tilde{g}}_{ij}(x_{ij})=s_i{g}_{ij}(x_{ij})+t_{ij}$$

$$\sum _{j=1}^m{b}_{ij}{x}_{ij}$$

$$\sum _{i=1}^n{x}_{ij}=c_j$$

$$x_{ij}\ge 0$$

$$\sum _{j=1}^m{x}_{ij}=r_i$$

$$\sum _{i=1}^n{\alpha }_i\sum _{j=1}^m{b}_{ij}{x}_{ij}$$

$$\sum _{i=1}^n{x}_{ij}=c_j$$

$$x_{ij}\ge 0$$

$$\displaystyle\big{(}x_{ij}=0\text{ and }\alpha_{i}b_{ij}\leq\beta_{j}\big{)}\text{ or }\big{(}x_{ij}\geq 0\text{ and }\alpha_{i}b_{ij}=\beta_{j}\big{)}\quad(i=1,\dots,n;\;j=1,\dots,m)$$

$$\sum _{i=1}^n{x}_{ij}=c_j(j=1,\dots ,m)$$

$$\sum _{j=1}^m{x}_{ij}=r_i(i=1,\dots ,n).$$

$$\sum _{i=1}^n\sum _{j=1}^m{f}_{ij}$$

$$\text{subject to}\sum _{i=1}^n{x}_{ij}$$

$$\sum _{j=1}^m{x}_{ij}$$

$$x_{ij}$$

$$f_{ij}^{\prime }(x)=\log g_{ij}^{\prime }(x).$$

$${\hat{g}}_{ij}^{\prime }(x)=\frac{1}{g_{ij}^{\prime }(x)}.$$

$$\sum _{i=1}^n\sum _{j=1}^m{a}_{ij}$$

$$\text{subject to}\sum _{i=1}^n{x}_{ij}$$

$$\sum _{j=1}^m{x}_{ij}$$

$$x_{ij}$$

$$\sum _{i=1}^n$$

$${\lambda }_i+{\mu }_j\ge a_{ij}(i=1,\dots ,n;j=1,\dots ,m).$$

$$a_{ij}=\log b_{ij}.$$

$$\big{(}x_{ij}=0\text{ and }a_{ij}\leq\lambda_{i}+\mu_{j}\big{)}\text{ or }\big{(}x_{ij}\geq 0\text{ and }a_{ij}=\lambda_{i}+\mu_{j}\big{)}\qquad(i=1,\dots,n;\;j=1,\dots,m),\\ \sum_{i=1}^{n}x_{ij}=c_{j}\quad(j=1,\dots,m),\qquad\sum_{j=1}^{m}x_{ij}=r_{i}\quad(i=1,\dots,n).\hskip 56.9055pt$$

$$f_{ij}^{\prime }(x_{ij})={\lambda }_i+{\mu }_j(i=1,\dots ,n;j=1,\dots ,m)$$

$$\sum _{i=1}^n{F}_{ij}({\lambda }_i+{\mu }_j)$$

$$\sum _{j=1}^m{F}_{ij}({\lambda }_i+{\mu }_j)$$

$${\beta }_j=\underset{i}{\max }{\alpha }_i{b}_{ij}(j=1,\dots ,m),{\hat{\alpha }}_i=\underset{j}{\min }{\beta }_j/b_{ij}(i=1,\dots ,n).$$

$$\underset{j}{\min }\frac{{\max }_i{\alpha }_i^{(k)}{b}_{ij}}{{\alpha }_i^{(k+1)}{b}_{ij}}=1.$$

$$f_{ij}(x_{ij},\eta )=a_{ij}{x}_{ij}+\eta f_1(x_{ij})$$

$$g_{ij}(x_{ij},\eta )=b_{ij}{g}_{\eta }(x_{ij})$$

Full text

1 Introduction

The theory of optimal transport is an extremely rich research area, with applications ranging from probability theory to mathematical physics, from image processing to partial differential equations, and from city planning to differential geometry. The two volumes by Rachev and Rüschendorf [1, 2] and the celebrated books by Fields medalist Villani [3, 4] provide ample evidence of the spectacular growth of the field since the original work of Monge [5], Kantorovich [6], and Fréchet [7]. Specialized treatments for applied mathematicians and for economists, respectively, are provided in the recent books by Santambrogio [8] and Galichon [9].

The purpose of the present paper is to highlight a connection between discrete optimal transport and a certain class of multi-objective problems. This connection appears to have gone unnoticed so far, even though the proof is not difficult, and the corresponding result for strictly concave problems has been known for several decades [10, 11, 12]. By means of this connection, the techniques of optimal transport can be applied to certain multi-objective problems, while at the same time an alternative perspective on optimal transport is obtained. As an illustration of the latter point, it will be shown below how the popular entropic regularization of optimal transport problems [13, 14] has a natural interpretation as an isoelastic regularization within the multi-objective framework. The transfer from a single objective to multiple objectives also entails a change from an additive to a multiplicative environment, and in this way leads naturally to applications of nonlinear Perron-Frobenius theory. This theory has been used in [15] (cf. also [16]) to analyze a matrix scaling problem that can be related to single-objective concave optimization, and it comes up perhaps even more naturally in the context of multi-objective problems [17]. The application of nonlinear Perron-Frobenius theory will be illustrated below.

In this paper, the term “multi-objective problem” is used to cover multi-agent problems (several agents, each with their own objective), as well as multicriteria or vector optimization problems (one agent with several objectives). The solution concept that will be employed can be used in both contexts, as illustrated in Section 2.2. The examples in this section also show the relevance of the concept in different fields: operations research, statistics, and actuarial science.

From the multi-agent point of view, the multi-objective problems considered in this paper may be placed within the broad class of fair division problems [18, 19]. Typical features of problems considered in this active area of research include multiple agents, each equipped with their own utility functionals, and resources in limited supply. Solution concepts are based on notions of efficiency and fairness. While many different notions of fairness have been proposed, the one that appears to be used most often is envy-freeness: no agent would prefer anybody else’s share to her own [20]. In the present paper, a different concept of fairness is used, which is less standard in the literature on fair division, and which is based on bringing in a notion of size that is independent of any particular agent. This fairness concept can be viewed as a version of “exogenous rights”​, mentioned as one of four principles of distributive justice by Moulin [19, Ch. 2], or of “entitlements” as discussed by Brams [18]. In contrast to the use of entitlements in [18], the size constraint is taken in this paper as an equality constraint, rather than as an inequality constraint. It will be assumed below that allotments to agents can be represented as vectors, and that size is expressed as a single number determined by a linear functional.

Size constraints can be used as a way to single out unique solutions among the set of all efficient solutions. In that sense, their role is comparable to, for instance, the Nash collective utility function [21]. The use of size constraints presupposes the availability of a measurement functional by which the magnitudes of allotments can be determined. In this sense, size constraints call for a richer problem environment than the Nash bargaining solution does. In addition to the availability of a measurement functional, the problem data should provide, for each of the agents/objectives, the value of the right-hand side of the equation that expresses the corresponding size constraint. Below, these numbers appear as row sums that are prefixed.

In the literature on multicriteria optimization (see for instance [22]), several methods have been proposed to select particular solutions among the set of all Pareto efficient solution, such as scalarization and lexicographic ordering. Scalarization is a way to relate a multi-objective problem to a single-objective problem. While this paper also establishes a connection between single-objective and multi-objective problems, and weighted-sum scalarization will be used as part of the proof technique, it should be emphasized that the connection established below cannot be viewed as a scalarization. The objective function in the associated single-objective problem cannot be obtained through aggregation of the objective functions in the multicriteria problem by means of an achievement function, as in [22, Section 4.6]. Multi-objective optimal transport has been studied for instance by Isermann [23], but again this is different from the problem considered in this paper, which is to connect optimal solutions of single-objective optimal transport problems to balanced solutions (as introduced below) of a class of multi-objective problems.

In the context of equilibrium economics, efficient solutions that satisfy budgetary size constraints have been termed BCPE solutions (budget-constrained Pareto efficient solutions) by Balasko [24]. Gale and Sobel [11] do not introduce a separate term for the combination of Pareto efficiency and satisfaction of budget constraints; they speak of optimal and fair solutions (optimal and proportional in [25]). Bühlmann and Jewell [12] use the term fair Pareto optimal risk exchanges (FAIRPOREX). Neither of these terminologies appears to be entirely appropriate to cover the range of different situations exemplified below. For generality and simplicity, the term balanced solutions will be used in this paper.111It is difficult to find expressions that have not already been used elsewhere. The term “balanced outcome” has been used by [26] to refer to an extension of the Nash bargaining solution to networks, which is different from but somewhat related to the notion used in this paper. In the Shapley-Bondereva theorem (see for instance [27]), balancedness is used in a quite different sense. In connection with the transportation problem, the term “balanced” is sometimes used to refer to a property that is called the global feasibility condition in this paper; see (1).

Individual rationality and behavioral effects do not play a role in the definition of balanced solutions. In some (but not all) of the examples shown below, it would in fact be quite unnatural to introduce these notions. All in all, it may perhaps be said that the concept of a balanced solution is closer in spirit to the multicriteria setting than to the multi-agent setting. In the case where multiple agents are involved, the concept represents a social planner’s point of view.

The organization of the paper is as follows. Following this introduction, a formal definition of the notion of balanced solutions is given in Section 2, along with three examples of situations in which this notion may be used. The same section also provides a characterization of balanced solutions. This characterization is used in Section 3 to establish the connection to optimal transport. Regularization and the connection to Perron-Frobenius theory are discussed in Section 4. A discussion of numerical issues follows in Section 5, and Section 6 concludes.

2 Balanced solutions of linear multi-objective problems

2.1 Definitions

In the optimal transport problem, the decision variables are organized into a matrix, and there are constraints related to row sums, as well as constraints related to column sums. Rows and columns typically correspond to objects that are different in nature, for instance producers and consumers. The problem is formulated with a single objective. As will be illustrated below, situations arise in practice which similarly call for a table to be filled, taking into account row sum constraints, as well as column sum constraints, but which are different in the sense that there are multiple objectives, which relate (by convention) to the rows of the matrix that is to be formed. Additionally, as shown in the examples below, there is often a certain status difference between the row constraints and the column constraints. The column constraints express “feasibility”, whereas the row constraints express a notion of “fairness”​. While the column constraints can be thought of being imposed by Nature, the row constraints are man-made, and the right-hand sides used in these constraints are to some extent a matter of choice. Such a hierarchy leads to the notion of what we will call “balanced” solutions of the multi-objective problem. The idea is that the objective functions together with the feasibility constraints are used to define a Pareto efficient surface, and that the row constraints are then used to select specific solutions on the efficient surface. Under suitable circumstances, it may happen that there is a unique balanced solution. A precise definition is given below. In the interest of generality, the terms “primary” and “secondary” are used to refer to the two types of constraints, rather than “feasibility” and “fairness”​.

Definition 2.1**.**

A multi-objective matrix allocation (MOMA) problem is specified by

• (i)

the decision variables $x_{ij}$, which together form a matrix $X$ of size $n\times m$;

• (ii)

objective functions given by

[TABLE]

• (iii)

primary constraints given by

[TABLE]

where $c_{1},\dots,c_{m}$ are given positive constants;

• (iv)

secondary constraints given by

[TABLE]

where $r_{1},\dots,r_{n}$ are given positive constants.

A specification of the problem above can be given in terms of the triple $(g,c,r)$, where $g=(g_{ij}(\cdot))_{i=1,\dots,n;\,j=1,\dots,m}$ is the matrix of reward functions, $c\in\mathbb{R}^{m}_{++}$ is the vector of given column sums, and $r\in\mathbb{R}^{n}_{++}$ is the vector of given row sums.222The notation $\mathbb{R}_{++}$ is used in this paper to indicate the set of strictly positive real numbers. The set of non-negative real numbers is denoted by $\mathbb{R}_{+}$. The functions $g_{ij}:\mathbb{R}_{+}\rightarrow\mathbb{R}\cup\{-\infty\}$ will be assumed to be twice continuously differentiable (for convenience), strictly increasing, and concave. An obvious necessary condition for a MOMA problem to be feasible is that the sum of the row sums equals the sum of the column sums:

[TABLE]

This condition will be referred to as the global feasibility condition. We take it as a standing assumption in all MOMA problems to be discussed below.

Definition 2.2**.**

An $n\times m$ matrix $X$ is said to provide a balanced solution of a given MOMA problem if it satisfies the following two conditions.

• (i)

The decision matrix $X$ is Pareto efficient with respect to the objective functions $\Psi_{i}$ and the primary constraints of the MOMA problem. In other words, there does not exist a matrix $\tilde{X}$ such that $\tilde{X}$ satisfies the primary constraints, $\Psi_{i}(\tilde{X}_{i\cdot})\geq\Psi_{i}(X_{i\cdot})$ for all $i=1,\dots,n$, and $\Psi_{i}(\tilde{X}_{i\cdot})>\Psi_{i}(X_{i\cdot})$ for at least one $i\in\{1,\dots,n\}$.

• (ii)

The matrix $X$ satisfies the secondary constraints of the MOMA problem.

Remark 2.3**.**

The formulation as given above refers to the point of view of maximization. An analogous formulation can be given for minimization problems. In the case of minimization, the decision variables represent quantities that are preferred to be small, at least from the perspective of the $i$-th objective, such as “payment to be made by agent $i$ if outcome $j$ occurs” or “part of chore $j$ to be carried out by agent $i$”​. For such problems, natural assumptions on the functions $g_{ij}$ (cost functions in this case, rather than reward functions) are that they are increasing and convex.

Remark 2.4**.**

In the above, we work with unweighted sums. A similar problem with weighted row and column sums

[TABLE]

($p_{i}>0$, $q_{j}>0$ for all $i$ and $j$) can be transformed to a problem with unweighted sums by the substitutions

[TABLE]

In this sense, the use of unweighted sums entails no loss of generality. On the other hand, it might be argued that such usage may not offer the best preparation for the continuous case. For simplicity of notation, however, in this paper unweighted sums will be used.

Remark 2.5**.**

If the reward functions $g_{ij}(x_{ij})$ are modified to

[TABLE]

where $s_{i}>0$ and $t_{ij}$ are constants ($i=1,\dots,n$; $j=1,\dots,m$), then any balanced solution of the original problem is also a balanced solution of the modified problem, and vice versa. This is immediate from the definition of balanced solutions, and from the additive separability of the objective functions.

Remark 2.6**.**

If $X$ is a balanced solution of the MOMA problem specified by $(g,c,r)$, and $s$ is any positive number, then $\tilde{X}:=X/s$ is a balanced solution of the problem specified by $(\tilde{g},r/s,c/s)$ with $\tilde{g}_{ij}(x_{ij}):=g_{ij}(sx_{ij})$, and vice versa. Under the global feasibility constraint, it is therefore no restriction of generality to assume that $\sum_{i}r_{i}=\sum_{j}c_{j}=1$.

A linear MOMA problem is one in which the reward functions are linear, i.e. $g_{ij}(x_{ij})=b_{ij}x_{ij}$ with constant coefficients $b_{ij}>0$.

2.2 Examples

Three examples are given of situations in which MOMA problems, and in particular linear MOMA problems, arise naturally. These examples refer to task allocation, statistical classification, and risk sharing respectively.

Example 2.7**.**

Mathematical formalization of the problem of chore division appears to have been first undertaken by Martin Gardner [28]. The problem is naturally formulated as a minimization problem, and the exogenous rights that were mentioned in Section 1 are in this case exogenous obligations. Consider a situation in which $n$ agents must carry out $m$ different tasks. For an example that may be familiar to many readers, one can think of $m$ courses, which need to be taught by $n$ available lecturers. The size of the contribution to be made by each of the agents is supposed to be prefixed; in the example of course allocation, this corresponds to a teaching load which is given in advance for each lecturer. The contribution size is computed by means of a linear formula, in which each of the tasks has its own nominal weight. The decision variable $x_{ij}$ represents the fraction of task $j$ that is assigned to agent $i$. The objective functions (cost in this case) are the disutilities of the agents that result from their package of assigned tasks. The primary constraints express that the task fractions must be non-negative, and that each task must be fulfilled completely. The secondary constraints correspond to the predetermined magnitudes of the contributions of the agents. The column sums are all 1 (i.e. 100%), and the row sums are given by the teaching loads. In this way, a MOMA problem is defined. If the agents’ disutilities depend linearly on the task fractions, then we have a linear MOMA problem.

Example 2.8**.**

Statistical classification is often considered as part of the field of machine learning [29]. The central problem in this area is to place objects into mutually exclusive categories, on the basis of observed characteristics. Examples of such situations are numerous and include, for instance, various types of diagnosis. To stay within the discrete context, it will be assumed here that only a finite number $m$ of different characteristics can be observed (in other words, the “feature space” is finite), and that there is a finite number $n$ of categories. The probability that a given object has characteristic $j$ and belongs to category $i$ will be denoted by $p_{ij}$. It is assumed that these probabilities are known. The situation in which $n=2$ corresponds to the problem of testing a simple null hypothesis against a simple alternative. Already in this case, the multiple-objective nature of the problem is apparent; the probability of a type-I error must be weighed against the probability of a type-II error. For a general number of categories, we can define with each category an objective function which is given by the probability of correct classification for an object in that category. Decisions may be randomized; we denote by $x_{ij}$ the probability by which the classifier will place an object with characteristic $j$ into category $i$. A problem of MOMA type arises when the percentage of objects to be placed in each of the categories is determined in advance. Such a prescription may be reasonable in circumstances in which classification is carried out for the purpose of treatment (for instance, medical treatment, or participation in an educational program), and there are capacity constraints for each form of treatment. The decision variables are given by the numbers $x_{ij}$. The probability that an object will be classified into category $i$, given that it indeed belongs to that category, is given by $\sum_{j=1}^{m}p_{ij}x_{ij}/\sum_{j=1}^{m}p_{ij}$ ($i=1,\dots,n$). These probabilities may be taken as objective functions; in the case of two categories, they are known as specificity and sensitivity, and are in a one-to-one relationship with the probabilities of type-I errors and type-II errors. The primary constraints are given by the non-negativity constraints on the decision variables $x_{ij}$, and by the condition that all objects with a given feature $j$ should be classified into one of the categories ($\sum_{i=1}^{n}x_{ij}=1$ for all $j$). The secondary constraints correspond to the prescribed percentages of objects to be placed in each category. Since the objective functions are linear in this case, we have a linear MOMA problem.

Example 2.9**.**

Risk sharing is a classical topic in actuarial science (see for instance [30]) and remains a very active field of research. While part of the literature on risk sharing is focused on information asymmetries and behavioral effects which we do not consider here, the problem has been formulated also in the MOMA framework [12]. For a simple illustration, consider $n$ individuals who form a group that will carry out a project with an uncertain payoff (say, a fishing expedition). For the purpose of discreteness, assume that there are $m$ possible outcomes. The variable $x_{ij}$ denotes the share of revenue that agent $i$ will receive in case outcome $j$ materializes. The right-hand side of the column constraint $\sum_{i=1}^{n}x_{ij}=c_{j}$ refers to the total amount that is available in outcome $j$. It would also be reasonable to impose the non-negativity constraint $x_{ij}\geq 0$. The decision variables $x_{ij}$ are to be determined at the beginning of the project, before the “veil of uncertainty” has been lifted. To express a notion of fairness, the group might agree on state prices $q_{j}$ for each of the possible outcomes, and fix a value for the allotment to each agent, i.e. $\sum_{j=1}^{m}q_{j}x_{ij}=v_{i}$, where the numbers $v_{i}$ are chosen in relation to the contributions (in terms of effort and money) made by participants. If the participants are risk neutral, i.e. their objectives are given by their subjective expected payoffs (in other words, $\Psi_{i}(X_{i\cdot})=\sum_{j=1}^{m}p_{ij}x_{ij}$ where the coefficients $p_{ij}$ refer to subjective probabilities), then we obtain a linear MOMA problem. In a more actuarial context, one might consider a group of companies which agree to a mutual reinsurance agreement with respect to claims that will come in during a given period. Fairness may then for instance be expressed on the basis of the expected claim sizes of the portfolios that the companies bring into the pool. In this case, the problem is more naturally formulated as a minimization problem, since the agents prefer smaller realized claims to larger ones.

2.3 Characterization of balanced solutions

Consider a linear MOMA problem given by

[TABLE]

It will be assumed throughout that the coefficients $b_{ij}$ and the prescribed row and column sums $r_{i}$ and $c_{j}$ are positive, and that the global feasibility condition (1) is satisfied. Under these conditions, the feasible set in decision space corresponding to the primary constraints is nonempty. The expressions above can refer to maximization as well as to minimization. In both cases, the coefficients $b_{ij}$ are supposed to be positive. Of course, the two cases differ by a reversal of the direction of the inequality used in the definition of Pareto efficiency.

It follows from a theorem of Isermann [31] (see also [22, Thm. 6.6]) that a matrix $X$ is a Pareto efficient solution of the problem consisting of the objectives and the primary constraints, taken as a maximization problem, if and only if it is a solution of a weighted-sum optimization problem of the form

[TABLE]

where $\alpha_{1},\dots,\alpha_{n}$ are positive parameters. As a result, we obtain the following characterization of balanced solutions of linear MOMA problems.

Lemma 2.10**.**

A matrix $X\in\mathbb{R}_{+}^{n\times m}$ is a balanced solution of the linear MOMA problem specified by (2), taken as a maximization problem, if and only if there exist positive numbers $\alpha_{1},\dots,\alpha_{n}$ (weights) and positive numbers $\beta_{1},\dots,\beta_{m}$ (dual variables relating to equality constraints), such that the following conditions are satisfied:

[TABLE]

Proof.

The weighted-sum optimization problem (3) is a linear program with equality and inequality constraints. By linear programming duality (see for instance [32]), a matrix $X$ is optimal for (3) if and only if there exist $\beta_{j}\in\mathbb{R}$ such that the complementary slackness conditions (4a) and the column constraints (4b) are satisfied. The set of all Pareto efficient solutions of the multi-objective linear program specified by the objective functions and the primary constraints of (2) is therefore given by matrices $X\in\mathbb{R}_{+}^{n\times m}$ for which there exist $\alpha_{i}>0$ and $\beta_{j}\in\mathbb{R}$ such that (4a) and (4b) hold. By the standing assumption that $b_{ij}>0$, this can only happen for dual variables which are positive. ∎∎

Remark 2.11**.**

An analogous statement to Lemma 2.10 can be proven in the case of minimization; the only difference lies in the direction of the inequality in the first term of the disjunction in (4a). Because of the constraints $x_{ij}\geq 0$ and $\sum_{i=1}^{n}x_{ij}=c_{j}>0$, for every index $j$ there must be a least one index $i$ such that $x_{ij}>0$; for this index $i$, the equality $\alpha_{i}b_{ij}=\beta_{j}$ holds. It follows that, also in the MOMA minimization case, the dual variables $\beta_{j}$ are positive.

3 Connection to optimal transport

MOMA maximization problems with increasing and strictly concave utility functions were considered by [11]. It is shown in this paper that there is a connection between balanced solutions of such problems and optimal solutions of related single-objective matrix allocation (SOMA) problems of the form

[TABLE]

Specifically, under regularity conditions, it was shown by [11] that a matrix $X$ is a balanced solution of the MOMA problem of Def. 2.1 if and only if $X$ is an optimal solution of the SOMA problem (5), where the reward functions $f_{ij}$ of the SOMA problem relate to the utility functions $g_{ij}$ of the MOMA problem by the prescription

[TABLE]

In other words, if one starts with the utility functions of the SOMA problem, then the utility functions of the corresponding MOMA problem are found from these by executing the following operations in succession: differentiation, exponentiation, integration. If one applies this procedure to a linear function, then the result is again linear, since the exponential of a constant is a constant. This suggests that linear MOMA problems should correspond to SOMA problems of the form (2.1) with linear reward functions. Such linear SOMA problems are known as optimal transport problems. In this section, we confirm the conjecture that there is a one-to-one relationship between solutions of the discrete optimal transport problem on the one hand, and solutions of related linear MOMA problems on the other hand.

Remark 3.1**.**

For MOMA minimization problems with increasing and strictly convex disutility functions, a result that is analogous to the one in [11] can be proven. The relation between the cost functions $g_{ij}(x)$ in the MOMA minimization problem and the cost functions $f_{ij}(x)$ in the SOMA minimization problem is given by the same formula (6) as in the maximization case.

Remark 3.2**.**

In the single-objective case, maximization and minimization problems can be simply converted into each other by the usual device of changing the sign of the objective function. The same is not true for MOMA problems in which we are looking for balanced solutions, since here we require that objective functions are increasing, both in the case of maximization and in the case of minimization. Instead, one can make use of the connection via the equivalent SOMA problem as specified in (6). It is seen that a matrix $X$ is a balanced solution of the MOMA maximization problem with reward functions $g_{ij}$ if and only if $X$ is a balanced solution of the MOMA minimization problem with cost functions $\hat{g}_{ij}$ (and the same row and column constraints), when the functions $g_{ij}$ and $\hat{g}_{ij}$ are related via

[TABLE]

The conjugation defined by (7) relates for instance $\log x$ (increasing, concave) to $\frac{1}{2}x^{2}$ (increasing, convex). This means that a logarithmic MOMA maximization problem is equivalent to a quadratic MOMA minimization problem. For another example, the increasing and concave power functions $x^{1-\gamma}/(1-\gamma)$ with $\gamma>0$, popular in utility theory, are linked by (7) to the increasing and convex power functions $x^{1+\gamma}/(1+\gamma)$. A linear MOMA maximization problem with coefficients $b_{ij}$ is equivalent to a linear MOMA minimization problem with coefficients $b_{ij}^{-1}$.

The optimal transport problems that will be considered below, in the maximization case, are of the discrete form

[TABLE]

The given row sums $r_{i}$ and column sums $c_{j}$ are all supposed to be positive, and the global feasibility condition $\sum_{i}r_{i}=\sum_{j}c_{j}$ is supposed to be satisfied. Under these conditions, the feasible set is not empty; for instance, we can take $x_{ij}=r_{i}c_{j}/s$ where $s:=\sum_{i}r_{i}=\sum_{j}c_{j}$. The problem (8) is known as the Monge-Kantorovich problem. It is a linear programming problem; its dual is given by

[TABLE]

Theorem 3.3**.**

A matrix $X$ is a balanced solution of the linear MOMA problem (2), taken as a maximization problem, if and only if $X$ is an optimal solution of the optimal transport problem (8) with the same row and column constraints, and

[TABLE]

Proof.

By linear programming duality, which in this case is Monge-Kantorovich duality (see for instance [9, Thm. 2.2]), a matrix $X\in\mathbb{R}^{n\times m}_{+}$ is a solution of the optimal transport problem (8) if and only if there are $\lambda_{i}\in\mathbb{R}$ ($i=1,\dots,n$) and $\mu_{j}\in\mathbb{R}$ ($j=1,\dots,m$) such that the following conditions hold:

[TABLE]

Under (10), these conditions are identical to the conditions (4), since we can take $\alpha_{i}=\exp(-\lambda_{i})$ and $\beta_{j}=\exp\mu_{j}$. The statement of the theorem follows. ∎∎

The value of the function $\sum_{i}\sum_{j}(\log b_{ij})x_{ij}$ for general $x_{ij}$ cannot be determined on the basis of information concerning only the values of the functions $\sum_{j}b_{ij}x_{ij}$. This shows that, even in the linear case, the relation between MOMA problems and SOMA problems as reflected in the theorem above cannot be viewed as a scalarization procedure. The term “social welfare function” used by Gale and Sobel [11] must be interpreted carefully; it should be noted that in the single-objective problem, the welfare is served of agents who are different from those that appear in the corresponding multi-objective problem.

Remark 3.4**.**

The relation (10) can be used to translate properties from linear SOMA problems to linear MOMA problems and vice versa. For instance, it is well known that the optimal transport problem can be solved very easily by a greedy algorithm when the weight matrix satisfies the Monge property, which requires that $a_{i_{1}j_{1}}+a_{i_{2}j_{2}}\geq a_{i_{1}j_{2}}+a_{i_{2}j_{1}}$ for all $1\leq i_{1}<i_{2}\leq n$ and $1\leq j_{1}<j_{2}\leq m$; see [33] for an extensive discussion.333The direction of the inequality is reversed here with respect to the convention in the cited paper, because we consider a maximization problem, rather than a minimization problem as in [33]. The corresponding condition for the related multi-objective problem is that the determinants of all 2-by-2 submatrices of the coefficient matrix $(b_{ij})$ should be non-negative.

4 Iterative solution method

4.1 Direct iteration

For strictly concave single-objective matrix allocation problems of the form (5), there is a natural iterative solution method, which may be viewed as an application of the alternating projections method due to Bregman [34]. The idea of solving matrix allocation problems by an iterative process that switches between row and column constraints already occurs in the early contributions of Kruithof [35] and Deming and Stephan [36]. The algorithm proposed in [35] and [36] is known as the Iterative Proportional Fitting Procedure (IPFP), or also as the RAS method. The iterative solution method for SOMA maximization problems that we describe below is a direct generalization of IPFP. For simplicity, it is assumed that the reward functions $f_{ij}(x)$ in (5) satisfy $\lim_{x\downarrow 0}f^{\prime}_{ij}(x)=\infty$ and $\lim_{x\rightarrow\infty}f_{ij}(x)=-\infty$; the corresponding conditions for the MOMA reward functions related via (6) are usually referred to as Inada conditions. These conditions imply, in particular, that the non-negativity constraints are redundant (optimal solutions must be strictly positive).

The standard Lagrange optimization method for the problem (5) gives rise to the equations

[TABLE]

where $\lambda_{i}$ ($i=1,\dots,n$) and $\mu_{j}$ ($j=1,\dots,m$) are Lagrange multipliers corresponding to row and column constraints respectively. Under the strict concavity assumptions on the functions $f_{ij}$, the derivatives $f_{ij}^{\prime}$ are strictly decreasing. Consequently, these functions have inverses, which will be denoted by $F_{ij}$. Under the Inada conditions, the inverse functions $F_{ij}$ are defined on all of the real line. We can then rewrite the row and column constraints in terms of the multipliers $\lambda_{i}$ and $\mu_{j}$:

[TABLE]

Starting with an initial guess for the values of $\lambda_{i}$, one can use the equations (13a) to compute corresponding multipliers $\mu_{j}$. Each multiplier $\mu_{j}$ appears in exactly one of the equations. Depending on the nature of the functions $F_{ij}$, it may be possible to write down the answer in analytic form, but even if this is not the case, the problem of finding $\mu_{j}$ only calls for determining the root of a scalar monotone function. When values for $\mu_{j}$ have been found in this way, one can use the equations (13b) to find updated values of the multipliers $\lambda_{i}$. Again the equation system is diagonal in the sense that each $\lambda_{i}$ occurs in exactly one of the equations. After solving for $\lambda_{i}$ ($i=1,\dots,n$), one can determine updates for $\mu_{j}$ ($j=1,\dots,m$), and so on.

Convergence of the iterative proportional fitting procedure and of its generalization to infinite dimensions has been studied extensively; see for instance [37, 38, 39, 40, 41, 42, 43]. Motivated by a MOMA problem in risk sharing, Pazdera et al. [17] proved convergence of the iterative procedure for a broad class of strictly concave utilities, on the basis of a nonlinear version of the Perron-Frobenius theorem due to [44]. One may write down an analogous iterative procedure for optimal transport problems, or equivalently for their multi-objective counterparts. While in this case the objective functions are linear rather than strictly concave, the iterative procedure can still be followed; however, as will be seen below, the procedure is not effective as a solution method. Nevertheless, the iteration in the linear case is worked out below, because it sheds light on the behavior of regularized versions.

We follow the MOMA viewpoint, which leads naturally to a multiplicative formulation. Iteration takes place between the utility weights $\alpha_{i}$ and the multipliers $\beta_{j}$. Starting with positive initial values for the weights $\alpha_{i}$ (for example, one might choose $\alpha_{i}=1$ for all $i$), the multipliers $\beta_{j}$ are to be found from the conditions for optimality of the weighted-sum problem (4a) and the column constraints. It is easily verified that, given a vector of weights $\alpha$, the vector of multipliers $\beta$ is such that (4a) as well as the column constraints are satisfied for some $x_{ij}$ if and only if $\beta_{j}=\max_{i}\alpha_{i}b_{ij}$ for $j=1,\dots,m$. Likewise, given a vector of multipliers $\beta$, the vector of weights $\alpha$ is such that (4a) as well as the row constraints are satisfied for some $x_{ij}$ if and only if $\alpha_{i}=\min_{j}\beta_{j}/b_{ij}$ for $i=1,\dots,n$. This leads to the iteration $\alpha^{(k+1)}=\varphi(\alpha^{(k)})$, where the mapping $\varphi$ from $\mathbb{R}^{n}_{+}$ into itself is defined by $\varphi(\alpha)=\hat{\alpha}$ with

[TABLE]

The behavior of this iterative procedure may be understood as follows. The update $\alpha_{i}^{(k+1)}$ is such that

[TABLE]

In other words, the new vector of weights $\alpha_{i}^{(k+1)}$ is determined such that $\alpha_{i}^{(k+1)}b_{ij}\leq\max_{i}\alpha_{i}^{(k)}b_{ij}$ for all $j$, with equality for at least one value of $j$. The quantity $\max_{i}\alpha_{i}^{(k)}b_{ij}$ may be called the column maximum for column $j$ in the row-scaled matrix with elements $\alpha_{i}^{(k)}b_{ij}$. If row $i$ in this matrix does not contain an element whose value is equal to the column maximum, then the row will receive a renewed scaling by multiplication by $\alpha_{i}^{(k+1)}$, so that the column maximum will be achieved. The renewed scaling is the minimal one which has this effect, and therefore it does not affect the value of any of the column maxima. After one iteration, each row in the row-scaled matrix contains a column maximal element, and the iteration arrives at a fixed point.

While convergence therefore takes place after at most one step, there are many fixed points. Any row scaling which is such that every row contains a column maximal element corresponds to a fixed point of the iteration. Moreover, the decision matrix cannot be written as a function of the weights $\alpha_{i}$ and the multipliers $\beta_{j}$, in contrast to the strictly concave case.

4.2 Regularization

A possible solution to remedy the unsatisfactory behavior of the direct iterative algorithm for the optimal transport problem as discussed above is to apply regularization. Here, we consider an additive regularization of the optimal transport problem, and we aim at finding the corresponding regularization in the multi-objective framework. In order to regularize, one can replace the linear reward functions $a_{ij}x_{ij}$ by concavified versions

[TABLE]

where the regularization parameter (or “temperature”) $\eta$ is a small positive number, and the regularization function $f_{1}(x)$ is strictly concave. The optimization problem defined by these regularized reward functions is strictly concave, and the relation (6) can be used to determine the corresponding multi-objective problem. Defining the function $g_{1}(x)$ (up to an irrelevant constant) from the regularization function $f_{1}(x)$ by $\log g_{1}^{\prime}(x)=f_{1}^{\prime}(x)$ as in (6), one finds that the MOMA regularization corresponding to (16) is given by

[TABLE]

where $b_{ij}=\exp(a_{ij})$ and the function $g_{\eta}(x)$ satisfies

[TABLE]

The relation (17) has a straightforward interpretation in terms of the Arrow-Pratt coefficient of risk aversion. The Arrow-Pratt risk aversion function $r_{\eta}(x)$ associated to the utility function $g_{\eta}(x)$ is given by $r_{\eta}(x)=-g^{\prime\prime}_{\eta}(x)/g^{\prime}_{\eta}(x)$. With this definition, the relation (17) is equivalent (up to an irrelevant constant) to

[TABLE]

This shows a rather natural relation between the risk aversion functions corresponding to $g_{\eta}$ and $g_{1}$.

A regularization function that is frequently chosen for the optimal transport problem is the entropic function given by

[TABLE]

[13, 14, 45, 9]. Up to a positive factor, which we can ignore on the basis of Remark 2.5, the corresponding regularization of the MOMA problem is $g_{1}(x)=1/x$. From (17), it follows that the MOMA regularization in this case can be written as

[TABLE]

The function that appears here is known as power utility, or also as isoelastic utility or CRRA utility (constant relative risk aversion). It is seen, therefore, that entropic regularization of the optimal transport problem corresponds to isoelastic regularization of the corresponding multi-objective matrix allocation problem.

In the case of entropic regularization of optimal transport problems, or, equivalently, isoelastic regularizations of linear MOMA problems, the iterative algorithm based on (13a–13b) specializes as follows. Working from the multiplicative form, the $k$-th iterate $x^{(k)}_{ij}$ is determined in terms of the corresponding iterates of the multiplicative parameters by444Usage of the notation $\alpha$ and $\beta$ here is in line with the notation in Section 2.3.

[TABLE]

In the special case (20), this produces

[TABLE]

Applying the column and row constraints successively leads to the iteration mapping $\varphi_{\eta}$ from $\mathbb{R}^{n}_{+}\setminus\{0\}$ into $\mathbb{R}^{n}_{++}$ defined by $\varphi_{\eta}(\alpha)=\hat{\alpha}$, where the vector $\hat{\alpha}$ is obtained from

[TABLE]

Here, $\|\cdot\|_{p}$ (with $p=1/\eta$) refers to the usual vector $p$-norm \|x\|_{p}=\big{(}\sum_{i=1}^{n}|x_{i}|^{p}\big{)}^{1/p}. The notations $\alpha_{\cdot}b_{\cdot j}$ and $b_{i\cdot}/\beta_{\cdot}$ are used to indicate the $n$-vector with entries $\alpha_{i}b_{ij}$ ($i=1,\dots,n$) and the $m$-vector with entries $b_{ij}/\beta_{j}$ ($j=1,\dots,m$) respectively. In the limit, as the regularization parameter $\eta$ tends to 0, the update mappings that are defined here converge to the ones in (14).

Remark 4.1**.**

The update equations in (21) can equivalently be written as follows:

[TABLE]

Instead of computing the sequences $\alpha_{i}^{(k)}$ and $\beta_{j}^{(k)}$, one can therefore also process the iteration with the sequences \tilde{\alpha}_{i}^{(k)}:=\big{(}\alpha_{i}^{(k)}\big{)}^{1/\eta} and \tilde{\beta}_{j}^{(k)}:=\big{(}\beta_{j}^{(k)}\big{)}^{1/\eta}. This alternative form of the iteration algorithm is the one presented by Benamou et al. [14]. It coincides with the IPFP algorithm, applied to the matrix with entries $b_{ij}^{1/\eta}=\exp(a_{ij}/\eta)$. The advantages and disadvantages of various possible formulations are discussed below in Section 5.

Remark 4.2**.**

As noted above, the IPFP algorithm corresponds to a SOMA minimization problem in which the cost function is given by relative entropy. It is readily verified, via the differentiation-exponentiation-integration route (6), that the related MOMA minimization problems refers to agents whose utility functions are quadratic. There is some historical interest to this observation. The IPFP method was presented in 1940 by Deming and Stephan [36] as an algorithm that achieves a best approximation in the quadratic sense to a given matrix, subject to the row and column constraints. The authors soon discovered, however, that their claim was mistaken [46]. Decades later, it was established by Ireland and Kullback [39] that the approximation computed via IPFP is optimal in the sense of Kullback-Leibler divergence (relative entropy), rather than in the quadratic sense. The error in [36] may be softened a bit by the observation that, in the MOMA interpretation, IPFP does solve a quadratic minimization problem.

4.3 Perron-Frobenius analysis

The mapping (21), when extended to take the value 0 at 0, is a homogeneous mapping of the nonnegative cone into itself. The convergence behavior of iterations defined by such mappings has been studied extensively, both in finite-dimensional and in infinite-dimensional spaces; see for instance [47, 48, 49, 16, 50]. A key tool in this literature is the Hilbert metric, which is defined for positive vectors of equal length by

[TABLE]

The Hilbert metric is in fact only a pseudo-metric on the positive cone, since points of the same ray have distance zero, but on the open unit simplex it is a true metric. The use of this metric for the study of eigenvectors of positive linear mappings was initiated in [51], while its application to nonlinear mappings starts with [52]. The purpose of this section is to show how a result of Oshime [44] in nonlinear Perron-Frobenius theory can be used to prove with little effort that the iteration (21) converges to a unique (up to a positive multiplicative factor) fixed point. This follows the development in [17], but the proof can be much shorter in the present case. The approach facilitates comparison with the ineffective iteration (14).

The notion of nonsectionality is needed, which may be thought of as a nonlinear version of the well known irreducibility condition of linear Perron-Frobenius theory. It is defined as follows. In the definition below, the notation $x_{R}$ where $x$ is an $n$-vector and $R$ is a subset of the index set $\{1,\dots,n\}$ denotes the subvector with entries $x_{i}$, $i\in R$.

Definition 4.3**.**

[44] A mapping $\varphi$ of $\mathbb{R}^{n}_{+}$ into itself is nonsectional if, for every decomposition of the index set $\{1,\dots,n\}$ into two complementary nonempty subsets $R$ and $S$, there exists $s\in S$ such that

• (i)

$(\varphi(x))_{s}>(\varphi(y))_{s}$ for all $x,\,y\in\mathbb{R}^{n}_{+}$ such that $x_{R}>y_{R}$ and $x_{S}=y_{S}>0$;

• (ii)

$(\varphi(x^{k}))_{s}\rightarrow\infty$ for all sequences $(x^{k})_{k=1,2,\dots}$ in $\mathbb{R}^{n}_{+}$ such that $x^{k}_{R}\rightarrow\infty$, while $x^{k}_{S}$ is fixed and positive.

The following theorem, as cited in [17], can be used to prove uniqueness of the fixed points of the iteration defined by (21).

Theorem 4.4**.**

[44, Thm. 8, Remark 2]* If a mapping $\varphi$ from $\mathbb{R}^{n}_{+}$ into itself is continuous, monotone, homogeneous, and nonsectional, then the mapping $\varphi$ has a positive eigenvector, which is unique up to scalar multiplication, with a positive associated eigenvalue. In other words, there exist a constant $\theta^{*}>0$ and a vector $x^{*}\in\mathbb{R}^{n}_{++}$ such that $\varphi(x^{*})=\theta^{*}x^{*}$, and if $\theta>0$ and $x\in\mathbb{R}^{n}_{++}$ are such that $\varphi(x)=\theta x$, then $x$ is a scalar multiple of $x^{*}$.*

The existence of a unique positive eigenvector to the mapping defined by (21) does not immediately imply the existence of a fixed point; it needs to be shown that the corresponding eigenvalue is 1. This follows from the global feasibility constraint (cf. similar arguments in [53, p. 246] and [17]).

Lemma 4.5**.**

If, in the iteration mapping (21), we have $\hat{\alpha}=\theta\alpha$ for $\alpha\in\mathbb{R}^{n}\setminus\{0\}$, then $\theta=1$.

Proof.

The equations that define $\hat{\alpha}$ from $\beta$ and $\beta$ from $\alpha$, written in implicit form, are

[TABLE]

The global feasibility constraint $\sum_{i}r_{i}=\sum_{j}c_{j}$ therefore implies that

[TABLE]

Consequently, if $\hat{\alpha}=\theta\alpha$, then $\theta$ must be equal to 1. ∎∎

We can verify the conditions in Thm. 4.4 for the mapping defined in (21), and compare them to the properties of the update mapping defined by (14). A mapping $\varphi$ from (a subset of) $\mathbb{R}^{n_{1}}$ to $\mathbb{R}^{n_{2}}$ is said to be strongly monotone if the conditions $x_{i}\geq y_{i}$ for $i=1,\dots,n_{1}$ and $x_{i}>y_{i}$ for some $i\in\{1,\dots,n_{1}\}$ imply $(\varphi(x))_{i}>(\varphi(y))_{i}$ for all $i=1,\dots,n_{2}$.

Proposition 4.6**.**

The mapping from the non-negative cone $\mathbb{R}_{+}$ to itself defined by (14) is continuous, monotone, and homogeneous. The mapping (21) from $\mathbb{R}_{+}$ into itself is continuous, strongly monotone, homogeneous, and nonsectional.

Proof.

Continuity, monotonicity, and homogeneity are immediate from the definitions, both for the non-regularized mapping and for the regularized version. To show strong monotonicity in the case of (21), note that, due to the assumed positiveness of all coefficients $b_{ij}$, an increase of any of the components of the weight vector $\alpha$ implies an increase of all of the components of the vector $\beta$, which in its turn implies an increase of all components of the image vector $\hat{\alpha}$. Strong monotonicity implies that condition (i) in the definition of nonsectionality is satisfied. Moreover, when one of the components of $\alpha$ tends to infinity, the same is true for all of the components of $\beta$, and consequently for all of the components of $\hat{\alpha}$. This shows that condition (ii) is satisfied by the mapping (21). ∎∎

Corollary 4.7**.**

The mapping $\varphi_{\eta}$ defined by (21) has a unique (up to multiplication by a positive constant) fixed point in $\mathbb{R}^{n}_{++}$. For any initial point $\alpha^{(0)}$ in $\mathbb{R}^{n}_{+}\setminus\{0\}$, the sequence of normalized iterates defined by $\alpha^{(k+1)}=\varphi_{\eta}(\alpha^{(k)})/\|\varphi_{\eta}(\alpha^{(k)})\|_{1}$ converges to the unique fixed point on the open unit simplex.

Proof.

The first claim follows from Prop. 4.6 and Thm. 4.4. By a basic result in Perron-Frobenius theory (see for instance [16, Ch. 2] and [17, Lemma 4.4]), a homogeneous and strongly monotone mapping from the positive cone into itself is contractive in the sense of the Hilbert metric. Given the uniqueness of the fixed point, convergence of the iteration then follows from a result of Nadler [54] (see also [17, Cor. 5.4]). ∎∎

The non-regularized mapping (14) does not satisfy the property of nonsectionality. Indeed, otherwise uniqueness of a fixed point would follow as in the corollary above, whereas it has already been argued that the mapping (14) has many fixed points, which are not related to each other by multiplication by a positive constant. To see the violation of nonsectionality more directly, consider a point $\alpha\in\mathbb{R}_{++}^{n}$ in which the first component $\alpha_{1}$ is large relative to the other components, so that, for all $j$, $\alpha_{1}b_{1j}>\alpha_{i}b_{ij}$ for $i=2,\dots,n$. A sufficiently small increase of the components $\alpha_{2},\dots,\alpha_{n}$ will then have no effect on any of the parameters $\beta_{j}$ defined in (14), and consequently there is no effect on the update $\hat{\alpha}$ either. This means that condition (i) in the definition of nonsectionality is not satisfied with $S=\{1\}$ and $R=\{2,\dots,n\}$. It also follows that the mapping defined by (14) is not strongly monotone.

5 Numerics

Alternative implementations of the regularization method for optimal transport problems can be derived from (21) and from (22). These formulas are mathematically equivalent, but the corresponding implementations may differ in their numerical properties, especially when the regularization parameter is close to 0. As noted in [14], when the parameter $\eta$ is very small, the computation of a number of the form $b^{1/\eta}$ may lead to numerical overflow or underflow. Such problems are less likely in the implementation (21), which uses the operations of raising a number to the power $\eta$ and computing the $p$-norm of a vector with $p=1/\eta$; both of these operations can be implemented without the creation of intermediate results that are very large or very small.555The size of the machine epsilon still imposes a bound on how small the parameter $\eta$ can be made. As soon as the $1/\eta$-norms that are computed become numerically equal to $\infty$-norms, and numbers raised to the power $\eta$ become numerically equal to 1, the algorithm (21) degenerates into the ineffective non-regularized version (14).

There are different forms of the IPFP algorithm. It can be written in a vector iteration form, starting from a given positive initial matrix $(x_{ij}^{(0)})$ and a positive initial vector $u^{(0)}$:

[TABLE]

Alternatively, a matrix iteration can be used:

[TABLE]

It can be verified immediately that, if the initial vector $u^{(0)}$ is chosen such that $u_{i}^{(0)}=1$ for all $i$, the two recursions produce the same sequence of matrices $x_{ij}^{(k)}$. One might say that (24) presents a cumulative form of the computation, whereas (25) is the corresponding incremental form. Corresponding to the vector recursion (21), a matrix iteration can analogously be defined as follows. Starting from a given matrix $(z_{ij}^{(0)})$, define

[TABLE]

where the incremental multipliers $s^{(k)}\in\mathbb{R}^{m}$ and $t^{(k)}\in\mathbb{R}^{n}$ are given by

[TABLE]

The matrices $(z_{ij}^{(k)})$ relate to the matrices $(x_{ij}^{(k)})$ as defined in (21) via $x_{ij}=z_{ij}^{1/\eta}$. Starting from the initialization $z_{ij}=b_{ij}$, one can therefore iterate the matrices $(z_{ij})$, and extract the approximate solution of the optimal transport problem $(x_{ij})$ at the end.

Extracting the approximate solution only at the end raises the question how to monitor the quality of the solution in intermediate stages, in order to obtain a stopping criterion for the iteration. By construction of the algorithm, after each full step, the row sums of the matrix $(x_{ij}^{(k)})$ are equal to the prescribed row sums $r_{i}$, up to machine precision. Therefore, it is natural to measure the accuracy of the solution by the deviation of the column sums from the prescribed values $c_{j}$. From the equality $\sum_{j=1}^{m}x_{ij}=r_{i}$ for all $i$, it follows that $\sum_{j=1}^{m}\sum_{i=1}^{n}x_{ij}=\sum_{j=1}^{m}c_{j}$. Consequently, the Hilbert metric can be used as a valid distance measure between the vector of column sums of the approximate solution and the vector $c$. We have

[TABLE]

where $s_{j}^{(k)}$ is defined in (27), and $\mathbbm{1}$ is the all-ones vector. Computation of the Hilbert metric is therefore an easy byproduct of the updating rules (26–27), and one can monitor the progress of the convergence without actually computing the intermediate approximate solutions $(x_{ij}^{(k)})$.

As noted above, the non-regularized vector iteration (14) can have at most one nontrivial step. In this step, it produces a row balancing ensuring that every row of the weight matrix $(b_{ij})$ contains at least one element that is maximal among all elements in the column that it belongs to. From a numerical perspective, this type of row equilibration is possibly useful.

An iterative algorithm to solve the regularized version of the optimal transport problem may now be written as follows. Inputs to the computation are: the weight matrix $(a_{ij})$; the column sums $c_{j}$; the row sums $r_{i}$; the chosen regularization parameter $\eta$; and, finally, the tolerance level that is chosen as a stopping criterion for the iteration.

Define $b_{ij}=\exp(a_{ij})$. 2. 2.

Define $\hat{b}_{ij}=\alpha_{i}b_{ij}$, with \alpha_{i}=\min_{j}\!\big{(}(\max_{i}b_{ij})/b_{ij}\big{)} (preliminary row equilibration). 3. 3.

Initialize \big{(}z_{ij}^{(0)}\big{)}=(\hat{b}_{ij}). 4. 4.

Carry out the iteration (26–27), until the quantity \frac{1}{\eta}\log\!\big{(}(\max_{j}s_{j}^{(k)})/(\min_{j}s_{j}^{(k)})\big{)} is less than the chosen tolerance level. 5. 5.

Compute the approximate solution by x_{ij}=\big{(}z^{(k)}_{ij}\big{)}^{1/\eta}, where $k$ is the index of the final iteration step.

Numerical results for a test problem are shown in the figures below. The weight matrix is defined by a discretization, on a grid of size $256\times 256$, of the function a(x,y)=\sin\!\big{(}4\pi\big{(}(x-0.5)^{2}+(y-0.5)^{2}\big{)}\big{)} defined on the unit square. Marginals are given by scaled discretizations of the functions $c(x)=|x-0.5|$ and $r(y)=|y-0.5|$. Fig. 1 gives a visual description of the problem data (darker cells correspond to higher values of $a_{ij}$). The results of computations are shown in Fig. 2 for $\eta=10^{-2}$, $\eta=10^{-3}$, and $\eta=10^{-4}$. The tolerance level is set at $0.01$ in each case. CPU time per iteration step is mainly determined by the grid size, and is not influenced much by the choice of the regularization parameter $\eta$. However, the number of iteration steps required to reach a given level of convergence appears to be approximately inversely proportional to $\eta$. The computation time is therefore also inversely proportional to $\eta$.

Remark 5.1**.**

The computation time can be reduced by carrying out the iteration in stages, starting with a relatively high value of the regularization parameter, and reducing it stepwise. Such stagewise calculation is standard in the interior point method of linear programming (see for instance [55]). A simple reduction scheme was applied to the problem data given above, consisting of twelve stages in which the regularization parameter was reduced by the factor 1.5 at each step to reach the final level $10^{-4}$, while the tolerance level was kept at $0.01$ in all stages. The stagewise computation required less than 1000 iterations to converge. Compared to the single-stage calculation that starts immediately with the final value of the regularization parameter, this implies a reduction of computation time by approximately a factor 25.

It is seen from the graphs of the criterion value (28) that the convergence is somewhat irregular. This behavior becomes even more pronounced in the case of small examples. Of course, in such examples there is not really a reason to apply regularization, since small instances can be easily solved by any linear programming code. However, one may be able to get more insight into the behavior of the iterative algorithm from simple examples. Consider a problem with the following data:

[TABLE]

The convergence behavior of the iterative algorithm is shown in Fig. 3. As is evident from the figure, the convergence towards the true solution is not smooth at all. There are long periods in which little progress is made, and the algorithm appears almost trapped. On the other hand, once the algorithm is close to the true solution, progress is quick; reducing the tolerance level from $10^{-2}$, as used here, to levels close to machine accuracy would increase the total number of iterations just by a few percent. Upon inspection of the intermediate results, it appears that the phases of stagnation are linked to approximate solutions that, when used as a starting point for the IPFP algorithm, would lead to cycling, due to the presence of too many zeros.666Precise conditions for convergence of the IPFP algorithm for non-negative weight matrices are given in [56]. The algorithm passes closely by the following row-balanced matrices in succession:777More precisely, the matrices given are of the form $(x_{ij}^{(k)})$ where $x^{(k)}_{ij}$ = $(z^{(k)}_{ij})^{1/\eta}$, and the matrices \big{(}z_{ij}^{(k)}\big{)} are the ones that are actually computed in the proposed algorithm.

[TABLE]

Each of the first three matrices above has the property that it would cause the IPFP algorithm to cycle; the final matrix corresponds to the true solution. When the regularization parameter is reduced, the algorithm follows essentially the same path, but takes smaller steps, and hence requires more iterations.

6 Conclusions

The purpose of this paper has been to establish a relationship between the optimal transport problem on the one hand, and a certain class of multi-objective equilibrium problems with linear reward functions on the other hand. This connection forms a natural complement to a similar relationship which exists in the strictly concave / strictly convex case. By the availability of this connection, it is possible to apply different perspectives to a given problem. This variety of viewpoints has been illustrated in the context of regularization, where it was shown that entropic regularization in the single-objective problem corresponds to isoelastic regularization in the multi-objective framework. The ability to switch between different interpretations may be useful in the design of algorithms. The curious behavior of the iterative algorithm for small values of the regularization parameter seems worthy of further investigation.