3 Optimal Flows

Up until now we have analyzed problems where network structure is either fixed or generated by some specified random process. In this chapter, we investigate networks where connections are determined endogenously via equilibrium or optimality conditions. In the process, we cover some of the most powerful methods available for solving optimization problems in networks and beyond, with applications ranging from traditional graph and network problems, such as trade, matching, and communication, through to machine learning, econometrics and finance.

3.1Shortest Paths¶

As a preliminary step, we study the shortest path problem—a topic that has applications in production, network design, artificial intelligence, transportation and many other fields. The solution method we adopt also happens to be one of the clearest illustrations of Bellman’s principle of optimality, which is one of the cornerstones of optimization theory and modern economic analysis.

3.1.1Definition and Examples¶

We start proceedings by introducing simple examples. (In the next section we will formalize the problem and consider solution methods.)

3.1.1.1Set Up¶

Consider a firm that wishes to ship a container from $A$ to $G$ at minimum cost, where $A$ and $G$ are vertices of the weighted digraph $\gG$ shown in Figure 3.1. Arrows (edges) indicate paths that can be used for freight, while weights indicate costs of traversing them. In this context, the weight function is also called the cost function, and we denote it by $c$. For example, $c(A, B)$ is the cost of traveling from vertex $A$ to vertex $B$.

Since this graph is small, we can find the minimum cost path visually. A quick scan shows that the minimum attainable cost is 8. Two paths realize this cost: $(A, C, F, G)$ and $(A, D, F, G)$, as shown in Figure 3.2 and Figure 3.3, respectively.

3.1.1.2A Recursive View¶

Let’s now consider a systematic solution that can be applied to larger graphs. Let $q^*(x)$ denote the minimum cost-to-go from vertex $x$. That is, $q^*(x)$ is the total cost of traveling from $x$ to $G$ if we take the best route. Its values are shown at each vertex in Figure 3.4. We can represent $q^*$ in vector form via

As is clear from studying Figure 3.4, once $q^*$ is known, the least cost path can be computed as follows: Start at $A$ and, from then on, at arbitrary vertex $x$, move to any $y$ that solves

Here $\oO(x) = \setntn{y \in V}{(x, y) \in E}$ is the set of direct successors of $x$, as defined in Section 1.4.1, while $c(x, y)$ is the cost of traveling from $x$ to $y$. In other words, to minimize the cost-to-go, we choose the next path to minimize current traveling cost plus cost-to-go from the resulting location.

Thus, if we know $q^*(x)$ at each $x$, then finding the best path reduces to the simple two stage optimization problem in (3.2).

But now another problem arises: how to find $q^*$ in more complicated cases, where the graph is large? One approach to this problem is to exploit the fact that

must hold for every vertex $x$ in the graph apart from $G$ (where $q^*(G)=0$).

Take the time to convince yourself that, for our example, the function $q^*$ satisfies (3.3). In particular, check that (3.3) holds at each nonterminal $x$ in Figure 3.4.

We can understand (3.3), which is usually called the Bellman equation, as a restriction on $q^*$ that helps us identify this vector. The main difficulty with our plan is that the Bellman equation is nonlinear in the unknown function $q^*$. Our strategy will be to convert this nonlinear equation into a fixed point problem, so that fixed point theory can be applied.

We do this in the context of a more general version of the shortest path problem.

3.1.2Bellman’s Method¶

We consider a generic shortest path problem on a flow network, which consists of a weighted digraph $\gG = (V, E, c)$ with a sink $d \in V$ called the destination and a weight function $c \colon E \to (0, \infty)$ that associates a positive cost to each edge $(x,y) \in E$. In what follows we take $|V|=n$.

We consider how to find the shortest (i.e., minimum cost) path from $x$ to $d$ for every $x \in V$. In doing so we impose the following assumption.

To simplify the arguments below, we consider a second digraph $\gG' = (V, E')$ where $V$ is unchanged and $E' = E \cup \{(d,d)\}$. In other words, $E'$ extends $E$ by adding a self-loop at $d$. We also extend $c$ to this self-loop by setting $c(d, d) = 0$. This just means that “terminating at $d$” is the same as “staying at $d$,” since no more cost accrues after arrival at the destination. Figure 3.5 illustrates in the context of the previous example.

We will use the intuition from the previous section, based around the Bellman equation, to construct a solution method.

3.1.2.1Policies¶

Instead of optimal paths, it turns out to be more convenient to aim for optimal policies. In general, a policy is a specification of how to act in every state. In the present setting, a policy is a map $\sigma \colon V \to V$ with $\sigma(x) = y$ understood as meaning “when at vertex $x$, go to $y$.” A policy is called feasible if $\sigma(x) \in \oO(x)$ for all $x \in V$.

For any feasible policy $\sigma$ and $x \in V$, the trajectory of $x$ under $\sigma$ is the path from $x$ to the destination indicated by the feasible policy. In other words, it is the sequence $(p_0, p_1, p_2, \ldots)$ defined by $p_0 = x$ and $p_i = \sigma(p_{i-1})$ for all $i \in \NN$. It can also be expressed as $(\sigma^i(x)) := (\sigma^i(x))_{i \geq 0}$, where $\sigma^i$ is the $i$-th composition of $\sigma$ with itself.

Let $\Sigma$ be the set of all feasible policies. In view of Exercise 3.1.1, policies in $\Sigma$ have the property that every trajectory they generate reaches $d$ in $n$ steps and stays there. In other words,

Given $q \in \RR^V_+$, we call $\sigma \in \Sigma$ $q$-greedy if

In essence, a greedy policy treats $q$ as the minimum cost-to-go function and picks out an optimal path under that assumption.

Using our new terminology, we can rephrase the discussion in Section 3.1.1.2 as follows: the shortest path problem can be solved by finding the true minimum cost-to-go function $q^*$ and then following a $q^*$-greedy policy. In the remainder of this section, we prove this claim more carefully.

3.1.2.2Cost of Policies¶

We need to be able to assess the cost of any given policy. To this end, for each $x \in V$ and $\sigma \in \Sigma$, let $q_\sigma(x)$ denote the cost of following $\sigma$ from $x$. That is,

This second equality holds because $\sigma^i(x) = d$ for all $i \geq n$ and $c(d, d) =0$. The function $q_\sigma \in \RR^V_+$ is called the cost-to-go under $\sigma$.

It will be helpful in what follows to design an operator such that $q_\sigma$ is a fixed point. For this purpose we let $U$ be all $q \in \RR_+^V$ with $q(d)=0$ and define $T_\sigma \colon U \to U$ by

Here and below, with $k \in \NN$, the expression $T_\sigma^k$ indicates the $k$-th composition of $T_\sigma$ with itself (i.e., $T_\sigma$ is applied $k$ times).

3.1.2.3Optimality¶

The minimum cost-to-go function $q^*$ is defined by

The definition of $q^*$ matches our intuitive definition from Section 3.1.1.2, in the sense that $q^*(x)$ is, for each $x \in V$, the minimum cost of traveling from $x$ to the destination $d$.

A policy $\sigma^* \in \Sigma$ is called optimal if it attains the minimum in the definition of $q^*$, so that $q^* = q_{\sigma^*}$ on $V$. Our main aims are to

Regarding the first step, we claim that $q^*$ satisfies the Bellman equation (3.3). To prove this claim we introduce the Bellman operator $T$ via

By construction, $q^*$ satisfies the Bellman equation (3.3) if and only if $Tq^* = q^*$.

The next result is central. It confirms that the minimum cost-to-go function satisfies the Bellman equation and also provides us with a means to compute it: pick any $q$ in $U$ and then iterate with $T$.

We now have a means to compute the minimum cost-to-go function (by iterating with $T$, starting at any $q \in U$) and, in addition, a way to verify the following key result.

3.1.2.4An Implementation in Julia¶

Let’s use the ideas set out above to solve the original shortest path problem, concerning shipping, which we introduced in Section 3.1.1.1. We will implement in Julia.

Our first step is to set up the cost function, which we store as an array called c. We identify the vertices $A$, ..., $G$, with the integers $1, \ldots, 7$. We set c[i, j] = Inf when no edge exists from i to j, so that such a path is never chosen when evaluating the Bellman operator, as defined in (3.3). When an edge does exist, we enter the cost shown in Figure 3.1.

c = fill(Inf, (7, 7))

c[1, 2], c[1, 3], c[1, 4] = 1, 5, 3
c[2, 4], c[2, 5] = 9, 6
c[3, 6] = 2
c[4, 6] = 4
c[5, 7] = 4
c[6, 7] = 1
c[7, 7] = 0

Next we define the Bellman operator:

function T(q)
    Tq = similar(q)
    n = length(q)
    for x in 1:n
        Tq[x] = minimum(c[x, :] + q[:]) 
    end
    return Tq
end

Now we arbitrarily set $q \equiv 0$, generate the sequence of iterations $Tq, T^2 q, T^3 q$ and plot them:

using PyPlot
fig, ax = plt.subplots()

n = 7
q = zeros(n)
ax.plot(1:n, q)

for i in 1:3
    new_q = T(q)
    ax.plot(1:n, new_q, "-o", alpha=0.7)
    global q = new_q
end

After adding some labels, the output looks like the image in Figure 3.6. Notice that, by $T^3 q$, we have already converged on $q^*$. You can confirm this by checking that the values of $T^3 q$ line up with those we obtained manually in Figure 3.4.

3.1.3Betweenness Centrality¶

In Section 1.4.3 we discussed a range of centrality measures for networks, including degree, eigenvector and Katz centrality. Aside from these, there is another well-known centrality measure, called betweenness centrality, that builds on the notion of shortest paths and is particularly popular in analysis of social and peer networks.

Formally, for a given graph $\gG = (V, E)$, directed or undirected, the betweenness centrality of vertex $v \in V$ is

where $S(x,y)$ is the set of all shortest paths from $x$ to $y$ and $S(x, v, y)$ is the set of all shortest paths from $x$ to $y$ that pass through $v$. (As usual, $|A|$ is the cardinality of $A$.) Thus, $b(v)$ is proportional to the probability that, for a randomly selected pair of nodes $x, y$, a randomly selected shortest path from $x$ to $y$ contains $v$. Thus, the measure will be high for nodes that “lie between” a large number of node pairs.

For example, Networkx stores a graph called florentine_families_graph that details the marriage relationships between 15 prominent Florentine families during the 15th Century. The data can be accessed via

import networkx as nx
G = nx.florentine_families_graph()    

The network is shown in Figure 3.7, which was created using

nx.draw_networkx(G, [params])    

where [params] stands for parameters listing node size, node color and other features. For this graph, node size and node color are scaled by betweenness centrality, which is calculated via

nx.betweenness_centrality(G)    

Although this graph is very simple, the output helps to illustrate the prominent role of the Medici family, consistent with their great wealth and influence in Florence and beyond.

3.2Linear Programming and Duality¶

Our study of shortest paths in Section 3.1 used a relatively specialized optimization method. In this section we cover more general results in optimization and duality, which will then be applied to endogenous networks, optimal transport and optimal flows. Part of the strategy is to take challenging optimization problems and regularize them, often by some form of convexification.

Before diving into theory, we use Section 3.2.1 below to provide motivation via one very specific application. This application involves what is typically called matching in economics and linear assignment in mathematics. We illustrate the key ideas in the context of matching workers to jobs. Later, when we discuss how to solve the problem, we will see the power of linear programming, convexification and duality.

3.2.1Linear Assignment¶

Behold the town of Springfield. A local employer is shutting down and 40 workers stand idled. Fortunately for these workers, Springfield lies within a political battleground state and, as a result, the mayor receives backing to attract a new employer. She succeeds by promising a certain firm that the 40 workers will be retrained for the new skills they require. For mathematical convenience, let us suppose that there are exactly 40 new positions, each with distinct skill requirements.

The team set up by the mayor records the individual skills of the 40 workers, along with the requirements of the new positions, and estimates the cost $c(i, j)$ of retraining individual $i$ for position $j$. The team’s challenge is to minimize the total cost of retraining. In particular, they wish to solve

where $\pP$ is the set of all permutations (i.e., bijective self-maps) on the integers $1, \ldots, 40$. Figure 3.8 illustrates one possible permutation.

The problem is discrete, so first order conditions are unavailable. Unsure of how to proceed but possessing a powerful computer, the team instructs its workstation to step through all possible permutations and record the one that leads to the lowest total retraining cost. The instruction set is given in Algorithm 2.

After five days of constant execution, the workstation is still running and the mayor grows impatient. The team starts to calculate how long execution will take. The main determinant is the size of the set $\pP$. Elementary combinatorics tells us that the number of permutations of a set of size $n$ is $n!$. Quick calculations show that $40!$ exceeds $8 \times 10^{47}$. A helpful team member points out that this is much less than the number of possible AES-256 password keys (approximately $10^{77}$). The mayor is not appeased and demands a runtime estimate.

Further calculations reveal the following: if, for each $\sigma$ in $\pP$, the workstation can evaluate the cost $t(\sigma) = \sum_{i=1}^{40} c(i, \sigma(i))$ in $10^{-10}$-th of a second (which is extremely optimistic), then the total run time would be

Another helpful team member provides perspective by noting that the sun will expand into a red giant and vaporize planet Earth in less than $10^{10}$ years.

The great computational cost of solving this problem by direct calculations is an example of what is often called the curse of dimensionality. This phrase, coined by Richard Bellman (1920-1984) during his fundamental research into dynamic optimization, refers to the exponential increase in processor cycles needed to solve computational problems to a given level of accuracy as the number of dimensions increases. The matching problem we have just described is high-dimensional because the choice variable $\sigma$, a permutation in $\pP$, is naturally associated with the vector $(\sigma(1), \ldots, \sigma(40))$. This, in turn, is a point in 40-dimensional vector space.^[1]

Fortunately, clever algorithms for this matching problem have been found and, for problems such as this one, useful approximations to the optimal allocation can be calculated relatively quickly. For example, Dantzig (1951) showed how such problems can be cast as a linear program, which we discuss below. Further progress has been made in recent years by adding regularization terms to the optimization problem that admit derivative-based methods.

In the sections that follow, we place matching problems in a more general setting and show how they can be solved efficiently. Our first step is to review the basics of linear programming.

3.2.2Linear Programming¶

A linear program is an optimization problem with a linear objective function and linear constraints. If your prior belief is that all interesting problems are nonlinear, then let us reassure you that linear programming is applicable to a vast array of interesting applications. One of these is, somewhat surprisingly, the assignment problem in Section 3.2.1, as famously demonstrated by the American mathematician George Bernard Dantzig (1914–2005) in Dantzig (1951). Other applications include optimal flows on networks and optimal transport problems, which in turn have diverse applications in economics, finance, engineering and machine learning.

We explain the key ideas below, beginning with an introduction to linear programming.

3.2.2.1A Firm Problem¶

One way to define linear programs is in terms of what open source and commercial solvers for linear programming problems handle. Typically, for these solvers, the problem takes the form

Here each $A_i$ is a matrix with $n$ columns and $c, b_1, b_2, d_\ell$ and $d_u$ are conformable column vectors.

To illustrate, let’s consider a simple example, which concerns a firm that fabricates products labeled $1, \ldots, n$. To make each product requires machine hours and labor. Product $i$ requires $m_i$ machine hours and $\ell_i$ labor hours per unit of output, as shown in the table below

A total of $M$ machine hours and $L$ labor hours are available. Letting $q_i$ denote output of product $i$, the firm’s problem is

subject to

Here $\pi$ is profits and $w$ and $r$ are the wage rate and rental rate of capital. We are taking these values as fixed, along with $L$ and $M$, so choosing the $q \in \RR^n_+$ that maximizes profits is the same as choosing the $q$ that maximizes revenue $\sum_{i=1}^n p_i q_i$. This will be our objective in what follows.

(Why are total labor supply and machine hours held fixed here? We might think of this problem as one of designing a daily or weekly production plan, which optimally allocates existing resources, given current prices.)

Figure 3.9 shows an illustration of a simple case with two products. Consistent with Exercise 3.2.3, the feasible set is a polyhedron, shaded in blue. The green lines are contour lines of the revenue function $(q_1, q_2) \mapsto p_1 q_1 + p_2 q_2$, with $p_1=3$ and $p_2=4$. By inspecting this problem visually, and recognizing that the contour lines are increasing as we move to the northeast, it is clear that the maximum is obtained at the extreme point indicated in the figure.

Maximizing revenue is equivalent to minimizing $\sum_i (- p_i) q_i$, so with $c = (-p_1, -p_2)$ we have a linear programming problem to which we can apply a solver.

3.2.2.2A Python Implementation¶

Let’s look at one option for solving this problem with Python, via SciPy’s open source solver linprog. For the simple two-product firm problem, we ignore the unnecessary equality constraint in (3.9) and set

The bound $0 \leq q$ is imposed by default so we do not need to specify $d_\ell$ and $d_u$.

We apply the numbers in Figure 3.9, so the maximization version of the problem is

Now we set up primitives and call linprog as follows:

A = ((2, 5),
     (4, 2))
b = (30, 20)
c = (-3, -4)  # minus in order to minimize

from scipy.optimize import linprog
result = linprog(c, A_ub=A, b_ub=b)
print(result.x)    

The output is [2.5, 5.0], indicating that we should set $q_1=2.5$ and $q_2=5.0$. The result is obtained via the simplex algorithm, discussed in Section 3.2.2.8.

3.2.2.3A Julia Implementation¶

When solving linear programs, one option is to use a domain specific modeling language to set out the objective and constraints in the optimization problem. In Python this can be accomplished using the open source libraries such as Pyomo and Google’s OR-Tools. In Julia we can use JuMP.

The following code illustrates the Julia case, using JuMP, applied to the firm problem (3.11).

using JuMP
using GLPK

m = Model()
set_optimizer(m, GLPK.Optimizer)

@variable(m, q1 >= 0)
@variable(m, q2 >= 0)
@constraint(m, 2q1 + 5q2 <= 30)
@constraint(m, 4q1 + 2q2 <= 20)
@objective(m, Max, 3q1 + 4q2)

optimize!(m)    

Notice how the JuMP modeling language allows us to write objectives and constraints as expressions, such as 2q1 + 5q2 <= 30. This brings the code closer to the mathematics and makes it highly readable.

If we now run value.(q1) and value.(q2), we get 2.5 and 5.0, respectively, which are the same as our previous solution.

3.2.2.4Standard Linear Programs¶

The programming problem in (3.8)–(3.9) is convenient for applications but somewhat cumbersome for theory. A more canonical version of the linear programming problem is

Here $x$ is interpreted as a column vector in $\RR^n$, $c$ is also a column vector in $\RR^n$, $A$ is $m \times n$, and $b$ is $m \times 1$. A linear program in the form of (3.12) is said to be in standard equality form. In preparation for our discussion of duality below, we also call (3.12) the primal problem.

Standard equality form is more general than it first appears. In fact the original formulation (3.8)–(3.9) can be manipulated into standard equality form via a sequence of transformations. Hence, when treating theory below, we can specialize to standard equality form without losing generality.

Although we omit full details on the set of necessary transformations (which can be found in Bertsimas & Tsitsiklis (1997) and many other sources), let’s gain some understanding by converting the firm optimization problem into standard equality form. To simplify notation, we address this task when $n=3$, although the general case is almost identical.

As above, we switch to minimization of a linear constraint by using the fact that maximizing revenue is equivalent to minimizing $\sum_i (- p_i) q_i$. Next, we need to convert the two inequality constraints (3.10) into equality constraints. We do this by introducing slack variables $s_m$ and $s_\ell$ and rewriting the constraints as

Indeed, we can see that requiring $\sum_{i=1}^3 m_i q_i + s_m = M$ and $s_m \geq 0$ is the same as imposing $\sum_{i=1}^3 m_i q_i \leq M$, and the same logic extends to the labor constraint.

Setting $x := (q_1, q_2, q_3, s_m, s_\ell)$, we can now express the problem as

subject to

This is a linear program in standard equality form.

3.2.2.5Duality for Linear Programs¶

One of the most important facts concerning linear programming is that strong duality always holds. Let us state the key results. The dual problem corresponding to the standard equality form linear program (3.12) is

Readers who have covered Section 6.1.7 in the appendix will be able to see the origins of this expression. In particular, by formula (6.19) in the appendix, the dual problem corresponding to the standard equality form linear program (3.12) can be expressed as

and $E = \RR^n_+$. (We can also treat the inequality $x \geq 0$ via a multiplier but this turns out to be unnecessary.) Now observe that

Since the dual problem is to maximize this expression over $\theta \in \RR^n$, we see immediately that a $\theta$ violating $A^\top \theta \leq c$ will never be chosen. Hence the dual to the primal problem (3.12) is (3.13).

3.2.2.6Strong Duality¶

As shown in Section 6.1.7.1 of the appendix, the inequality $D \leq P$ always holds, where $P$ is the primal value in (3.12). This is called weak duality. If $P = D$, then strong duality is said to hold. Unlike weak duality, strong duality requires conditions on the primitives.

The next theorem states that, for linear programs, strong duality holds whenever a solution exists. A proof can be obtained either through Dantzig’s simplex method or via Farkas’ Lemma. See, for example, Ch. 4 of Bertsimas & Tsitsiklis (1997) or Ch. 6 of Matousek & Gärtner (2007).

Strong duality of linear programs has many important roles. Some of these are algorithmic: duality can be used to devise efficient solution methods for linear programming problems. Another way that duality matters for economists is that it provides deep links between optimality and competitive equilibria, as we show below.

3.2.2.7Complementary Slackness¶

We say that $x^* \geq 0$ and $\theta^* \in \RR^m$ satisfy the complementary slackness conditions for the linear program (3.12) when

While it is possible to derive these expressions from the complementary slackness in the KKT conditions in Section 6.1.7.3, a better approach is to connect them directly to the saddle point condition.

To see how this works, suppose that $x^* \geq 0$ is feasible for the primal problem and $\theta^* \in \RR^m$ is feasible for the dual problem. If $(x^*, \theta^*)$ is a saddle point of $L$, then the complementary slackness conditions (3.15)–(3.16) must hold. Indeed, (3.15) is trivial when $x^*$ is feasible, since $Ax^* = b$. At the same time, (3.16) must be true because violation implies that

due to dual feasibility (i.e., $A^\top \theta^* \leq c$) and $x^* \geq 0$. But then $(x^*, \theta^*)$ is not a saddle point of $L$, since changing the $j$-th element of $x^*$ to zero strictly reduces the Lagrangian

By sharpening this argument, it is possible to show that, for linear programs, the complementary slackness conditions are exact necessary and sufficient conditions for a saddle point of the Lagrangian. This leads to the next theorem.

Chapter 4 of Bertsimas & Tsitsiklis (1997) provides more discussion and a full proof of Theorem 3.2.2. Below, we illustrate some of the elegant connections between complementary slackness and equilibria in competitive economic environments.

3.2.2.8The Simplex Algorithm¶

Linear programming is challenging in high-dimensional settings, partly because the linear objective function implies that solutions are not interior. The first efficient solution methods for linear programs appeared in the 1930s and 1940s, starting with the simplex method of Kantorovich and Dantzig. As the simplex method began to prove its worth, linear programming grew into a technique of enormous practical importance. Operations research, communication and production systems began to make heavy use of linear programs.

The simplex algorithm makes use of the following result, which is proved in Theorem 2.7 of Bertsimas & Tsitsiklis (1997).

(An extreme point of a polyhedron was defined in Section 6.1.5.1. The feasible set was shown to be a polyhedron in Exercise 3.2.3.)

The simplex algorithm makes use of Theorem 3.2.3 by walking along edges of the polyhedron that forms the feasible set, from one extreme point to another, seeking at each step a new extreme point (which coincides with a vertex of the polyhedron) that strictly lowers the value of the objective function. Details can be found in Bertsimas & Tsitsiklis (1997) and Matousek & Gärtner (2007).

3.3Optimal Transport¶

Next we turn to optimal flows across networks. One simple—but computationally nontrivial—example of a network flow problem is the linear assignment problem we analyzed in Section 3.2.1. There, the vertices are either workers or jobs and the edges are assignments, chosen optimally to minimize aggregate cost. More general network flow problems extend these ideas, allowing endogenous formation of networks in more sophisticated settings. The general structure is that vertices are given, while edges and weights are chosen to optimize some criterion. There are natural applications of these ideas in trade, transportation and communication, as well as less obvious applications within economics, finance, statistics and machine learning.

In our study of network flows, we begin with the optimal transport problem, which is the most important special case of the general network flow problem. (In fact, as we show in Section 3.3.4.3, there exists a technique by which the general network flow problem can always be reduced to an optimal transport problem.)

3.3.1The Monge-Kantorovich Problem¶

Optimal transport is a classic problem dating back to the work of Gaspard Monge (1746–1818), who studied, among many other things, the transport of earth for construction of forts. This simple-sounding problem—how to transport a pile of given shape into a new pile of given shape at minimum cost—can, after normalization, be identified with the problem of least cost transformation of one distribution into another distribution. Figure 3.10 gives a visualization of transforming one distribution into another in one dimension (although the cost function is not specified).

It turns out that, by varying the notion of cost, this transportation problem provides a highly flexible method for comparing the distance between two distributions, with the essential idea being that distributions are regarded as “close” if one can be transformed into the other at low cost. The resulting distance metric finds wide-ranging and important applications in statistics, machine learning and various branches of applied mathematics.^[2]

In economics, optimal transport has important applications in transportation and trade networks, as well as in matching problems, econometrics, finance and so on (see, e.g., Galichon (2018)). As such, it is not surprising that economists have contributed a great deal to the optimal transport problem, with deep and fundamental work being accomplished by the great Russian mathematical economist Leonid Kantorovich (1912–1986), as well as the Dutch economist Tjalling Koopmans (1910–1985), who shared the Nobel Prize with Kantorovich in 1975 for their work on linear programming and optimal transport.^[3]

We start our discussion with a straightforward presentation of the mathematics. Then, in Section 3.3.3, we will connect the mathematics to economic problems, and show the deep connections between transport, duality, complementary slackness and competitive equilibria.

3.3.1.1Monge’s Formulation¶

We start with the classical problem of Monge, which is simple to explain. We take as given two finite sets $\Xsf$ and $\Ysf$, paired with distributions $\phi \in \dD(\Xsf)$ and $\psi \in \dD(\Ysf)$. Elements of $\Xsf$ and $\Ysf$ are called locations. To avoid tedious side cases, we assume throughout that $\phi$ and $\psi$ are strictly positive on their domains. In addition, we are supplied with a cost function $c \colon \Xsf \times \Ysf \to \RR_+$. Our goal is to “transport” $\phi$ into $\psi$ at minimum cost. That is, we seek to solve

over the set of all maps $T$ from $\Xsf$ onto $\Ysf$ satisfying

The constraint says $T$ must be such that, for each target location $y$, the sum of all probability mass sent to $y$ is equal to the specified quantity $\psi(y)$. The symbol $\sum_x$ is short for $\sum_{x \in \Xsf}$. In this context, $T$ is often called a Monge map.

The problem is easily illustrated in the current discrete setting. Figure 3.11 gives a visualization when locations are enumerated as $\Xsf = \{x_1, \ldots, x_7\}$ and $\Ysf = \{y_1, y_2, y_3, y_4\}$, with both $\Xsf$ and $\Ysf$ being subsets of $\RR^2$. For simplicity, $\phi(x_i)$ is written as $\phi_i$ and similarly for $\psi(y_j)$. Vertex size is proportional to probability mass assigned to the vertex. The edges represent one feasible Monge map.

Discreteness and lack of convexity in the constraint (3.18) imply that the Monge problem is, in general, hard to solve. Truly fundamental progress had to wait until Kantorovich showed how convexification can greatly simplify the problem. The convexification process requires shifting the problem to a higher dimensional space, but the cost of higher dimensions is outweighed by the regularization provided by convexity and lack of discreteness. We study the Kantorovich formulation in Section 3.3.1.3.

3.3.1.2Assignment as Optimal Transport¶

The linear assignment problem studied in Section 3.2.1, with cost $c(i, j)$ of training worker $i$ for job $j$, is a special case of optimal transport. All we have to do is set $\Xsf = \Ysf = \natset{n}$ and take $\phi$ and $\psi$ to be discrete uniform distributions on $\natset{n}$.

Since $T$ must be a bijection on $\natset{n}$, which is also a permutation of $\natset{n}$, the objective of the optimal transport problem under the current configuration is

where $\pP$ is the set of all permutations of $\natset{n}$. This is the same optimization problem as the linear assignment problem in (3.7).

3.3.1.3Kantorovich’s Relaxation of the Monge Problem¶

The basic idea in Kantorovich’s relaxation of the Monge problem is to allow mass located at arbitrary $x \in \Xsf$ to be mapped to multiple locations in $\Ysf$, rather than just one. This means that we are no longer seeking a function $T$, since, by definition, a function can map a given point to only one image. Instead, we seek a “transport plan” that sends some fraction $\pi(x, y)$ of the mass at $x$ to $y$. The plan is constrained by the requirement that, for all $x$ and $y$,

These constraints on this transport plan mean that it takes the form of a “coupling,” which we define below.

Figure 3.12 illustrates a feasible transport plan in the discrete setting. As with Figure 3.11, $\phi(x_i)$ is written as $\phi_i$ and similarly for $\psi(y_j)$, while vertex size is proportional to probability mass. Unlike the Monge setting of Figure 3.11, the mass at each vertex $\phi_i$ can be shared across multiple $\psi_j$, as long as the constraints are respected.

Let’s write the constraints more carefully. We recall that, in probability theory, a coupling is a joint distribution with specific marginals. More precisely, given $\phi$ in $\dD(\Xsf)$ and $\psi$ in $\dD(\Ysf)$, a coupling of $(\phi, \psi)$ is an element $\pi$ of $\dD(\Xsf \times \Ysf)$ with marginals $\phi$ and $\psi$. This restriction on marginals means that

The constraints in (3.20)–(3.21) require that

In the present setting, a coupling is also called a transport plan.

Let $\Pi(\phi, \psi)$ be the set of all couplings of $\psi$ and $\phi$. Taking $\phi, \psi$ and the cost function $c$ as given, the general Monge–Kantorovich problem, also called the optimal transport problem, is to solve

where

is the Frobenius inner product of $c$ and $\pi$, treated as $|\Xsf| \times |\Ysf|$ matrices. The sum in $\inner{c, \pi}_F$ measures the total cost of transporting $\phi$ into $\psi$ under the plan $\pi$. There is linearity embedded in this cost formulation, since doubling the amount sent from $x$ to $y$ scales the associated cost at rate $c(x, y)$.

We call any $\pi$ solving (3.22) an optimal plan. Since we are maximizing over a finite set, at least one such plan exists.

3.3.1.4Optimal Transport as a Linear Program¶

With some relatively simple manipulations, the general optimal transport problem can be mapped into a standard equality form linear program. This provides two significant benefits. First, we can apply duality theory, which yields important insights. Second, on the computational side, we can use linear program solvers to calculate optimal plans.

To map the optimal transport problem into a linear program, we need to convert matrices into vectors. Will use the $\vecop$ operator, which takes an arbitrary $A \in \matset{n}{m}$ and maps it to a vector in $\RR^{nm}$ by stacking its columns vertically. For example,

In this section we adopt the notational conventions in Exercise 3.3.5. The objective function $\inner{c, \pi}_F$ for the optimal transport problem can now be expressed as $\vecop(c)^\top \vecop(\pi)$.

To rewrite the constraints in terms of $\vecop(\pi)$, we use the Kronecker product, which is denoted by $\otimes$ and defined as follows. Suppose $A$ is an $m \times s$ matrix with entries $(a_{ij})$ and that $B$ is an $n \times t$ matrix. The Kronecker product $A \otimes B$ of $A$ and $B$ is the $mn \times st$ matrix defined, in block matrix form, by

It can be shown that Kronecker products and the $\vecop$ operator are connected by the following relationship: for conformable matrices $A$, $B$ and $M$, we have

Using (3.24) and the symbol $I_k$ for the $k \times k$ identity matrix, we can rewrite the first constraint in (3.23) as

Now, using block matrix notation and setting

the optimal transport problem can be expressed as the standard equality form linear program

Finally, for a given solution $x$, the transport plan is recovered by inverting $x = vec(\pi)$.

3.3.1.5Implementation¶

Listing 1 is a function that implements the above steps, given flat (one-dimensional) arrays phi and psi representing the distributions over the source and target locations, plus a two-dimensional array c representing transport costs. (The method argument highs-ipm tells linprog to use a particular interior point method, details of which can be found in the linprog documentation. Simplex and other methods give similar results.)

import numpy as np
from scipy.optimize import linprog

def ot_solver(phi, psi, c, method='highs-ipm'):
    """
    Solve the OT problem associated with distributions phi, psi 
    and cost matrix c.

    Parameters
    ----------
    phi : 1-D array
        Distribution over the source locations.
    psi : 1-D array
        Distribution over the target locations.
    c : 2-D array
        Cost matrix.
    """
    n, m = len(phi), len(psi)

    # vectorize c
    c_vec = c.reshape((m * n, 1), order='F') 

    # Construct A and b
    A1 = np.kron(np.ones((1, m)), np.identity(n))
    A2 = np.kron(np.identity(m), np.ones((1, n)))
    A  = np.vstack((A1, A2))
    b  = np.hstack((phi, psi))

    # Call solver
    res = linprog(c_vec, A_eq=A, b_eq=b, method=method)

    # Invert the vec operation to get the solution as a matrix
    pi = res.x.reshape((n, m), order='F')
    return pi

Notice that in Listing 1, the reshape order is specified to F. This tells NumPy to reshape with Fortran column-major order, which coincides with the definition of the $\vecop$ operator described in Section 3.3.1.4. (The Python vectorize operation defaults to row-major order, which concatenates rows rather than stacking columns. In contrast, Julia uses column-major by default.)

Let’s call this function for the very simple problem

With these primitives, all mass from $\phi_1$ and $\phi_2$ should be sent to $\psi_1$. To implement this problem we set

phi = np.array((0.5, 0.5))
psi = np.array((1, 0))
c = np.ones((2, 2))

and then call ot_solver via

ot_solver(phi, psi, c)

The output is as expected:

array([[0.5, 0. ],
       [0.5, 0. ]])    

3.3.1.6Python Optimal Transport¶

In the case of Python, the steps above have been automated by the Python Optimal Transport package, due to Flamary et al. (2021). For the simple problem from Section 3.3.1.5 we run

import ot
ot.emd(phi, psi, c)   # Use simplex method via the emd solver

The output is again equal to

array([[0.5, 0. ],
       [0.5, 0. ]])    

Figure 3.13 shows an example of an optimal transport problem solved using the Python Optimal Transport package. The interpretation is similar to Figure 3.12, although the number of vertices is larger. In addition, the edges show the optimal transport configuration, in the sense that $\pi^*$, the optimal transport plan, is treated as the adjacency matrix for the graph. (The figure shows the unweighted graph, with an arrow drawn from $\phi_i$ to $\psi_j$ whenever $\pi^*_{ij} > 0$.) The optimal transport plan is obtained by converting the transport problem into a linear program, as just described, and applying the simplex method. Although there are 32 nodes of each type, the problem is solved by the simplex routine in less than one millisecond.

3.3.1.7Kantorovich Relaxation and Linear Assignment¶

We showed in Section 3.3.1.2 that the linear assignment problem studied in Section 3.2.1 is a special case of optimal transport, where $\phi$ and $\psi$ become discrete uniform distributions on $\natset{n}$. Moreover, as we have just seen, Kantorovich’s relaxation method allows us to apply linear programming. This leads to fast solutions.

In the discussion in Section 3.2.1 we used $n=40$. Let’s start here with $n=4$, to illustrate the method, and then try with $n=40$. The matrix $c(i,j)$ of costs, which is the only other primitive, will be generated randomly as an independent array of uniform random variables: Then we apply the Python Optimal Transport (POT) library, as in Section 3.3.1.6.

Here is our set up:

import numpy as np
import ot

n = 4
phi = np.ones(n) 
psi = np.ones(n)

We have broken the rule that $\phi$ and $\psi$ should sum to one. This could be fixed easily by using np.ones(n)/n instead of np.ones(n), but the POT library does not care (as long as np.sum(phi) equals np.sum(psi)) and, moreover, the idea of putting unit mass everywhere is natural, since each element of $\phi$ represents one worker, and each element of $\psi$ represents one job.

Now we build the cost matrix:

c = np.random.uniform(size=(n, n))

The output is

array([[0.03841, 0.32896, 0.55989, 0.41682],
       [0.91527, 0.24566, 0.26022, 0.64510],
       [0.96275, 0.44089, 0.79274, 0.93065],
       [0.40454, 0.87307, 0.43555, 0.54903]])

(For example, the cost of retraining worker 1 for job 2 is 0.32896.) Finally, we call the solver:

ot.emd(phi, psi, c)    

The output is

array([[1., 0., 0., 0.],
       [0., 0., 1., 0.],
       [0., 1., 0., 0.],
       [0., 0., 0., 1.]])    

This is a permutation matrix, which provides another way to express a permutation of $\natset{n}$. The first row tells us that worker 1 is assigned to job 1, the second tells us that worker 2 is assigned to job 3, and so on.

If we now set $n=40$ and rerun the code, the line ot.emd(phi, psi, c), which calls the simplex-based solver, runs in less than 1 millisecond on a mid-range laptop. This is a remarkable improvement on the $2.5 \times 10^{30}$ year estimate for the brute force solver we obtained in Section 3.2.1.

3.3.1.8Tight Relaxation¶

Notice that the solution we obtained for the linear assignment problem using the simplex method does not split mass, as permitted by the Kantorovich relaxation. For example, we do not send half of a worker to one job and the other half to another. This is convenient but why does it hold?

While we omit the details, the basic idea is that a general Kantorovich transport plan is a bistochastic matrix, and all such matrices can be formed as convex combinations of permutation matrices. (This is called Birkhoff’s Theorem.) In other words, the permutation matrices are extreme points of the set of bistochastic matrices. Moreover, Theorem 3.2.3 tells us that any optimizer of a linear program will be an extreme point—in this case, a permutation matrix.

3.3.2Kantorovich Duality¶

One of the greatest achievements of Kantorovich was to show that the optimal transport problem can be connected to a dual problem, and how that dual problem can be used to characterize solutions. This work anticipated much of the later development of duality theory for arbitrary linear programs.

Throughout this section, in stating the main results, we use the notation

This is just the usual inner product, when we think of $f$ and $\phi$ as vectors in $\RR^{|\Xsf|}$. Also, given a cost function $c$ on $\Xsf \times \Ysf$, let $\fF_c$ be all pairs $(w, p)$ in $\RR^\Xsf \times \RR^\Ysf$ such that

One part of Kantorovich’s duality results runs as follows.

Theorem 3.3.1 can now be understood as a special case of the more general result that strong duality holds for linear programs, which we stated in Theorem 3.2.1. In that spirit, let us verify Theorem 3.3.1 using Theorem 3.2.1, by working with the linear programming formulation of the optimal transport problem provided in Section 3.3.1.4.

To do this, we take that formulation, which is stated in (3.27), and apply the dual formula in (3.13), which yields

If we write the argument $\theta \in \RR^{n+m}$ as $(-w, p)$, so that we now maximize over the two components $-w \in \RR^n$ and $p \in \RR^m$, as well as transposing the constraint, we get

where $w \in \RR^n$ and $p \in \RR^m$. By using the definition of the Kronecker product and carefully writing out the individual terms, it can be shown that the constraint in this expression is equivalent to requiring that $p_j - w_i \leq c_{ij}$ for all $(i,j) \in \natset{n} \times \natset{m}$. Recalling that $\Xsf$ has been mapped to $\natset{n}$ and $\Ysf$ has been mapped to $\natset{m}$, this is exactly the same restriction as (3.28).

At this point it is clear that (3.29) is nothing but the dual of the linear program formed from the optimal transport problem. The claims in Theorem 3.3.1 now follow directly from the strong duality of linear programs (Theorem 3.2.1).

3.3.3Optimal Transport and Competitive Equilibria¶

The other major achievement of Kantorovich in the context of duality theory for optimal transport was to connect optimality of transport plans with the existence of functions $w, p$ from the dual problem such that a version of the complementary slackness conditions holds. Here we present this result, not in the original direct formulation, but rather through the lens of a competitive equilibrium problem. In doing so, we illustrate some of the deep connections between prices, decentralized equilibria and efficient allocations.

3.3.3.1The Advisor’s Problem¶

We imagine the following scenario. Iron is mined at a finite collection of sites, which we denote by $\Xsf$. We identify an element $x \in \Xsf$ with a point $(a, b) \in \RR^2$, which can be understood as the location of the mine in question on a map. At the wish of the queen, who seeks to defend the empire from greedy rivals, this iron is converted to swords by blacksmiths. There are a number of talented blacksmiths in this country, located at sites given by $\Ysf$. As for $\Xsf$, each $y \in \Ysf$ indicates a point in $\RR^2$. Henceforth, we refer to “mine $x$” rather than “the mine at $x$” and so on.

Each month, mine $x$ produces $\phi(x)$ ounces of iron ore, while blacksmith $y$ consumes $\psi(y)$ ounces. We take these quantities as fixed. We assume that total supply equals total demand, so that $\sum_x \phi(x) = \sum_y \psi(y)$. For convenience, we normalize this sum to unity. As a result, $\phi$ and $\psi$ are elements of $\dD(\Xsf)$ and $\dD(\Ysf)$ respectively.

The cost of transporting from $x$ to $y$ is known and given by $c(x,y)$ per ounce. The king’s chief advisor is tasked with allocating and transporting iron from the mines to the blacksmiths, such that each blacksmith $y$ receives their desired quantity $\psi(y)$, at minimum cost. A small amount of thought will convince you that the advisor’s problem is a version of the optimal transport problem (3.22). We call this the primal problem in what follows.

Operating in the days before Kantorovich, Dantzig and the electronic computer, the advisor employs a large team of bean counters, instructing them to find the allocation with least cost by trying different combinations. However, after a few days, she realizes the futility of the task. (With infinite divisibility, which corresponds to our mathematical model, the number of allocations is infinite. If we replace infinite divisibility with a finite approximation, the scale can easily be as large as that of the matching problem discussed in Section 3.2.1, with only a moderate number of mines and blacksmiths.)

3.3.3.2The Guild’s Problem¶

At this point she has another idea. There is a guild of traveling salesmen, who buy goods in one town and sell them in another. She seeks out the guild master and asks him to bid for the project along the following lines. The guild will pay the queen’s treasury $w(x)$ per ounce for iron ore at mine $x$. It will then sell the iron at price $p(y)$ per ounce to the queen’s representative at blacksmith $y$. The difference can be pocketed by the guild, as long as all blacksmiths are provided with their desired quantities. The guild master is asked to propose price functions $w$ and $p$.

The guild master sees at once that $p$ and $w$ must satisfy $p(y) - w(x) \leq c(x, y)$ at each $x, y$, for otherwise the advisor, who is no ones fool, will see immediately that money could be saved by organizing the transportation herself. Given this constraint, the guild master seeks to maximize aggregate profits, which is $\sum_y p(y) \psi(y) - \sum_x w(x) \phi(x)$. At this point it will be clear to you that the problem of the guild master is exactly that of Kantorovich’s dual problem, as given in Theorem 3.3.1.

Since the advisor has given up on her team of bean counters, the guild master employs them, and asks them to produce the optimal pair of prices. The bean counters set to work, trying different combinations of prices that satisfy the constraints. However, without a systematic methodology to follow or fast computers to turn to, their progress is slow. The advisor begins to fear that the coming war will be over before the guild master replies.

3.3.3.3Decentralized Equilibrium¶

At this point, it occurs to the advisor that yet another approach exists: privatize the mines, abolish the guild, and let the traveling salesmen, mine owners and blacksmiths make individual choices in order to maximize their profits. Purchase and sales prices, as well as quantities transported from each mine to each blacksmith, will be determined by the free market.

Although the advisor predates Kantorovich, she reasons that competition will prevent each salesman from profiteering, while the desire for profits will encourage high levels of transportation and minimal waste. It turns out that this idea works amazingly well, in the sense that we now describe.

For the record, we define a competitive equilibrium for this market as pair of price vectors $(w, p)$ in $\RR^\Xsf \times \RR^\Ysf$ and a set of quantities $\pi \colon \Xsf \times \Ysf \to \RR_+$ such that the following three conditions hold: For all $(x, y)$ in $\Xsf \times \Ysf$,