Half-quadratic transportation problems

TL;DR

This paper introduces a memory-efficient primal-dual algorithm for solving relaxed, differentiable transportation problems by approximating the cost function with weighted quadratic problems, enabling solutions to non-convex cases.

Contribution

It proposes a novel approach that approximates the transportation problem with a differentiable cost, allowing for efficient solutions to non-convex problems.

Findings

Efficient memory usage in solving transportation problems.

Ability to handle non-convex, differentiable cost functions.

Approximation method via weighted quadratic transportation problems.

Abstract

We present a primal--dual memory efficient algorithm for solving a relaxed version of the general transportation problem. Our approach approximates the original cost function with a differentiable one that is solved as a sequence of weighted quadratic transportation problems. The new formulation allows us to solve differentiable, non-- convex transportation problems.

Equations

$${\min }_x$$

$${\sum }_j{x}_{ij}=p_i$$

$$f(t)=\underset{0<\omega \le M}{\inf }\left\{\omega t^2+\psi (\omega )\right\};$$

$${\omega }^{\ast }=f^{\prime }(t)/(2t).$$

$${\min }_x{\sum }_{ij}{f}_{ij}(x_{ij})={\min }_{x,w}{\sum }_{ij}\{{\omega }_{ij}{x}_{ij}^2+\psi ({\omega }_{ij})\}.$$

$${\omega }^{k+1}={\mathrm{argmin}}_{\omega }\mathcal{L}(x^k,\omega ,{\mathbf{y}}^k)$$

$$x^{k+1},{\mathbf{y}}^{k+1}=\underset{x}{\mathrm{argmin}}\underset{\mathbf{y}}{\mathrm{argmax}}\mathcal{L}(x,{\omega }^{k+1},\mathbf{y}).$$

$${\omega }_{ij}=f_{ij}^{\prime }(x_{ij})/(2x_{ij}).$$

$${\omega }_{ij}{x}_{ij}-{\lambda }_i-{\gamma }_j-s_{ij}$$

$${\sum }_j{x}_{ij}-p_i$$

$${\sum }_i{x}_{ij}-q_j$$

$$s_{ij}{x}_{ij}$$

$$s_{ij},x_{ij}$$

$$x_{ij}={\bar{\omega }}_{ij}\left({\lambda }_i+{\gamma }_j+s_{ij}\right),$$

$${\sum }_j{\bar{\omega }}_{ij}\left({\lambda }_i+{\gamma }_j+s_{ij}\right)=p_i.$$

$${\lambda }_i=\frac{p_i-{\sum }_j({\gamma }_j+s_{ij}){\bar{\omega }}_{ij}}{{\sum }_j{\bar{\omega }}_{ij}}.$$

$${\gamma }_j=\frac{q_j-{\sum }_i({\lambda }_i+s_{ij}){\bar{\omega }}_{ij}}{{\sum }_i{\bar{\omega }}_{ij}}.$$

$$s_{ij}=\max \{0,-({\lambda }_i+{\gamma }_j)\}.$$

$${\lim }_{\beta \to 0}{c}_{ij}{\left(x_{ij}^2+{\beta }^2\right)}^{\frac{1}{2}}=c_{ij}{x}_{ij};c_{ij},x_{ij}\ge 0.$$

$${\tilde{f}}_{ij}(x_{ij})=c_{ij}{x}_{ij}^2/({\beta }^2+x_{ij}^2)$$

$$\lim_{\beta\rightarrow 0}\frac{t^{2}}{\beta^{2}+t^{2}}~{}=~{}1-\hat{\delta}(t)~{}=~{}\left\{\begin{array}[]{rcl}0&\mathrm{for}&t=0,\\ 1&\mathrm{for}&t\in\mathbb{R}\backslash\{0\}.\end{array}\right.$$

$${\omega }_{ij}{x}_{ij}-{\lambda }_i-{\gamma }_j-s_{ij}=b_{ij}$$

Full text

Keywords: General transportation problem, Half–quadratic potentials, Sequential quadratic programming.

1 Introduction

The general transportation problem (GTP) deals with the distribution of goods from $m$ suppliers with production capacities $p=\{p_{i}\}_{i=1,\ldots,m}$ to $n$ destinations with demands $q=\{q_{i}\}_{i=1,\ldots,n}$. Without loss of generality, we assume balanced production and demand: $\sum_{i}p_{i}=\sum_{j}q_{j}$. A classical approach to this problem assumes that the cost of transport remains constant, independently of the the quantity to be transported. In real problems, this is not the case. The cost may increase or decrease according the volume of the transported good. We can write the general transportation problem as follows:

[TABLE]

where $x_{ij}$ denotes the quantity of goods to be transported from the $i$th supplier to the $j$th destination, $f:[m\times n]\times\mathbb{R^{+}}\rightarrow\mathbb{R^{+}}$ is a continuous cost function that depends on the supplier $i$, the destination $j$ and the volume $x_{ij}$.

The first case of a transportation problem was formulated by Hitchcock assuming that the cost functions are linear: $f_{ij}(x_{ij})=c_{ij}x_{ij}$ [1]. Another popular model is the quadratic one: $f_{ij}(x_{ij})=a_{ij}x_{ij}^{2}+b_{ij}x_{ij}$. According to Ref. [2], the quadratic model is popular because it can approximate other cost functions. Despite such flexibility, the limitation of quadratic models has been well documented in the context of robust statistics, and its applications to image processing and computer vision [3, 4]. In our opinion, the main limitations of quadratic models are the following:

The impossibility of limiting the cost for large values; i.e., one has $\lim_{x\rightarrow+\infty}|q_{ij}(x_{ij})|=\infty$. 2. 2.

The limitation of promoting sparse solutions; i.e., solutions that use a reduced number of routes.

In this work, we present an approximation scheme that allows us to define new cost functions that overcome the aforementioned limitations. We also present a primal–dual algorithm with limited memory requirements.

The GTP is relevant in modern computer science applications such as computer vision [5], machine learning [6] and data analysis [7, 8]. The Earth Mover Distance (EMD) is an interesting application of the transportation problem where the optimum cost is used as a metric between the histograms, vectors, $p$ and $q$. In Ref. [5], EMD is used as a metric for image retrieval in the context of computer vision. Recently, the EMD was proposed as a measure of reconstruction error for non–negative matrix factorisation [9]. The Word Moving Distance is the metric version for comparing documents based on the transportation problem [8]. Recent EMD applications include the quantification of biological differences in flow cytometry samples [10]. In addition, there is current interest in the learning of metrics for particular problems [11]; since our proposal is parametrised, the parameters involved can be learned.

2 Preliminaries

Before presenting our transportation formulation, we review an important result reported in the context of robust statistics and continuous optimisation applied to image processing [4]. The purpose of such work was to transform some non–linear cost functions to half–quadratic functions [3, 4]. A half–quadratic function is quadratic in the original variable and convex in a new auxiliary variable, where the minima of the auxiliar variable can be computed with a closed formula. The next proposition resumes the conditions imposed on the cost function $f$ and the transformed half–quadratic function.

Proposition 1**.**

Let $f:\mathbb{R^{+}}\rightarrow\mathbb{R^{+}}$ be a function that fulfils the following conditions:

$f(t)\geq m$* with $f(0)=m$, for $t\geq 0$ and $m>-\infty$.* 2. 2.

$f$* is continuously differentiable.* 3. 3.

$f^{\prime}(t)\geq 0$. 4. 4.

$\lim_{t\rightarrow+\infty}f^{\prime}(t)/(2t)=0$. 5. 5.

$\lim_{t\rightarrow+0^{+}}f^{\prime}(t)/(2t)=M$, $0<M<+\infty$.

Then,

there exists a strictly convex and decreasing function $\psi:(0,M]\rightarrow[0,\beta)$, where $\beta=\textstyle\lim_{t\rightarrow+\infty}\left\{f(t)-t^{2}f^{\prime}(t)/(2t)\right\}$ such that

[TABLE] 2. 2.

the solution to $\inf_{0<\omega\leq M}\{\omega t^{2}+\psi(\omega)\}$ is unique and given by

[TABLE]

Proof.

The proof is presented in Ref. [4]. ∎

Observe that our version of the half–quadratic Proposition assumes a non-negativity constraint on the primal variables.

3 Half–quadratic transportation problem

In this section we present a memory efficient primal–dual algorithm for solving GTP which cost functions satisfy Proposition 1.

Proposition 2**.**

Let $f_{ij}$ be a cost function in (1), such that each $f$ satisfies Proposition 1; then, a solution to the transportation problem can be computed with Algorithm 1.

Proof.

On the half quadratic transportation problem. By (5), the cost (1) can be rewritten as

[TABLE] 2. 2.

On the algorithm convergence. Let $\mathcal{L}$ denotes the Lagrangian of the half–quadratic transportation problem, then one can interchange the order of the minimisations; i.e., $\min_{x,\omega}\max_{\bf y}\mathcal{L}(x,\omega,{\bf y})=\min_{x}\min_{\omega}\max_{\bf y}\mathcal{L}(\cdot)=\min_{\omega}\min_{x}\max_{\bf y}\mathcal{L}(\cdot);$ where we denote with $\bf y$ the Lagrange’s multiplies vectors. This suggests an alternating minimisation scheme w.r.t. $\omega$ and $(x,{\bf y}$). Let $x^{k}$, $\omega^{k}$ and ${\bf y}^{k}$ be the current feasible values, then we define

[TABLE]

to be the updated $\omega$ value. Thus, $x$ and ${\bf y}$ are updated by solving the quadratic transportation problem:

[TABLE]

We define $F(x)=\sum_{ij}f_{ij}(x_{ij})$ and $\hat{F}(x,\omega)=\sum_{ij}\{\omega_{ij}x^{2}_{ij}+\psi(\omega_{ij})\}$ and observe that $F(x^{k})=\hat{F}(x^{k},\omega^{k})\geq\hat{F}(x^{k},\omega^{k+1})\geq\hat{F}(x^{k+1},\omega^{k+1})=F(x^{k+1})$. Then, the alternated minimisations w.r.t. $\omega$ and $x$ produce a feasible convergent sequence $\{x^{k},x^{k+1},x^{k+2},\ldots\}$ that reduces the cost of the GTP: $F(x^{k})\geq F(x^{k+1})\geq F(x^{k+2})\geq\ldots$. 3. 3.

On the alternated minimisations. From (6), the optimum $\omega$ in (8) is computed as

[TABLE]

for a given $x$. We define ${\bf y}$ equal to $(\lambda,\gamma,s)$ where $\lambda$ and $\gamma$ are the Lagrange’s multipliers for the equality constraints (12) and (13), respectively; and $s$ are the Lagrange’s multipliers for the non–negativity constraint. Then, the minimisation (9) corresponds to finding the vectors $(x,\lambda,\gamma,s)$ that solve the Karush-Kuhn-Tucker conditions (KKTs) with $\omega$ fixed:

[TABLE]

A strategy for solving the KKTs is to use an iterative Projected Gauss–Seidel scheme [12]. Thus, from (11):

[TABLE]

where we define $\bar{\omega}_{ij}=1/\omega_{ij}$. Substituting $x_{ij}$ in (12), we have

[TABLE]

We solve for $\lambda_{i}$ and obtain

[TABLE]

Similarly, we substitute $x_{ij}$ in (13) and solve for $\gamma_{j}$:

[TABLE]

From (11), (14) and (15), we see two cases: $x_{ij}=0$ and $s_{ij}=-(\lambda_{i}+\gamma_{j})\geq 0$; or $x_{ij}=\bar{\omega}_{ij}(\lambda_{i}+\gamma_{j})\geq 0$ and $s_{ij}=0$. Thus

[TABLE]

The complete procedure is shown in Algorithm 1.

∎

In practice, we observed that the internal loop in Algorithm 1 requires only a few iterarations to approximate the dual variables and is not necessary to achieve convergence; in our experiments we used five iterations. The solution is the global minima or a local one by depending if the cost is or is not convex. Following, we present two particular and interesting cases of half-quadratic transportation problems.

Proposition 3** (Approximation for Linear Transportation).**

The linear cost can be approximated with the differential function [4]: $f_{ij}(ij)=c_{ij}({x_{ij}^{2}+\beta^{2}})^{\frac{1}{2}},$ with $\beta\approx 0$. Thus, $\omega$ is computed with $\textstyle\omega_{ij}={\frac{1}{2}}{c_{ij}}/({x_{ij}^{2}+\beta^{2}})^{\frac{1}{2}}$ in Algorithm 1.

Proof.

It follows from

[TABLE]

Hence, the formula for $\omega$ follows directly from (10). ∎

Proposition 4** (Approximation for $L_{0}$ Transportation).**

The cost function of the form $f_{ij}(x_{ij})=c_{ij}[1-\hat{\delta}(x_{ij})]$, for $c_{ij},x_{ij}>0$ (where $\hat{\delta}(t)$ is the Kroneker’s delta) can be approximated with the half-quadratic function:

[TABLE]

with $\beta\approx 0$. Thus, $\omega_{ij}={c_{ij}\beta^{2}}/{(\beta^{2}+x_{ij}^{2})^{2}}$.

Proof.

It follows from

[TABLE]

The formula for $\omega$ follows directly from (10). ∎

Figure 1 plots the half–quadratic functions that approximate the $L_{1}$ and $L_{0}$ norms.

4 Relationship with the quadratic transportation problem

Proposition 5** (Simple Quadratic Transport (SQT)).**

The SQT problem is defined by the cost function $f_{ij}(ij)=c_{ij}x_{ij}^{2}$; thus, $\omega_{ij}=c_{ij}$.

Proof.

It follows directly from (10). ∎

Proposition 6** (Simple Quadratic Transportation Algorithm).**

In the case of the QT problem, the Algorithm 1 is reduced to the dual Algorithm 2.

Proof.

From Proposition 5, we note that $\omega$ is constant, and that the computations of $\lambda$, $\gamma$ and $s$ are independent of $x$. ∎

Dorigo and Tobler discussed the relationship between the QTP and the push–pull migration laws implemented in Algorithm 2 [13].

Proposition 7** (Quadratic Transportation (QT)).**

The QT is defined by a cost function of the form $f_{ij}(ij)=a_{ij}x_{ij}^{2}+b_{ij}x_{ij}$. Thus, the dual algorithm is derived with $\omega=c_{ij}$ and using the condition

[TABLE]

instead of (11).

Proof.

It follows directly from the KKTs. ∎

Remark. An alternative to Proposition 7 is given by the half–quadratic approximation $a_{ij}x_{ij}^{2}+b_{ij}x_{ij}\approx a_{ij}x_{ij}^{2}+2b_{ij}({{x_{ij}^{2}+\beta^{2}}})^{\frac{1}{2}}$, with $\beta\approx 0$; thus, $\omega_{ij}=a_{ij}+\frac{1}{2}{b_{ij}}/({{x_{ij}^{2}+\beta^{2}}})^{\frac{1}{2}}$. This approximation is presented with the sole aim of illustrating the potential of our approach. It is clear that the dual algorithm derived according to Proposition 7 is more accurate, faster and requires less memory to be implemented.

5 Discussion and Conclusions

The transportation problem is the base of the Earth Mover Distance which has become a relevant metric to for compare distributions in applications to data analysis and computer vision. The presented technique can motivate the design of new algorithms in those areas.

In order to demonstrate the versatility of our proposal, we generate two random vectors $p$ and $q$ (depicted in Figure 2); we compute the optimum transported volumes $x$ with three cost function models: the quadratic ($c_{ij}x_{ij}^{2}$), the approximation $L_{1}$ ($c_{ij}\sqrt{x_{ij}^{2}+\beta^{2}}$, with $\beta^{2}=1\times 10^{-3}$) and the approximation $L_{0}$ ($c_{ij}x_{ij}^{2}/(\beta^{2}+x_{ij}^{2})$, with $\beta^{2}=1\times 10^{-1}$). In all the cases, we use $c_{ij}=|i-j|+1$. Figure 3 depicts the computed $x$ values. One can observe that the quadratic cost function promotes dense solutions; i.e., there are many $x$’s with small values. On the other hand, one can observe the sparseness of the solution is induced with the use of the approximated $L_{1}$–norm. Such sparsity is emphasised with the approximated $L_{0}$–norm.

We have presented a model to approximate solutions to the general transportation problems by approximating the transportation cost functions with half–quadratic functions. The approach guarantees convergence using an alternated minimisation scheme. In the case of a non–convex cost function $f$ the convergence is guaranteed to a local minimum. Although we present a minimisation algorithm with reduced memory requirements, our scheme accepts other efficient solvers for the quadratic transportation subproblem; such as those reported in Refs. [2, 14, 15].