A Sparse Smoothing Newton Method for Solving Discrete Optimal Transport Problems

Recall that a primal-dual IPM solves the KKT system (3) indirectly by solving a sequence of perturbed KKT systems of the form:

via the Newton method, where $\mu>0$ is the barrier parameter that is driven to zero. At each IPM iteration, the following Newton equation is solved:

where $(x,y,z)\in\mathbb{R}^{mn}_{++}\times\mathbb{R}^{m+n}\times\mathbb{R}^{mn}_{++}$ is the current primal-dual iterate, $(\Delta x,\Delta y,\Delta z)\in\mathbb{R}^{mn}\times\mathbb{R}^{m+n}\times\mathbb{R}^{mn}$ is the Newton direction, $\mu$ is set to be $\langle x,z\rangle/(mn)$, and $\gamma\in[0,1]$ is a given parameter. To reduce the computational cost, the large-scale system of $2mn+m+n$ linear equations is typically reduced to the so-called normal system of $m+n$ equations in $\Delta y$ through the block elimination of $\Delta x$ and $\Delta z$. Specifically, we solve the normal equation:

where $\Theta:=\mathrm{Diag}(\theta)\in\mathbb{R}^{mn\times mn}$ is a diagonal matrix whose diagonal entries are $\theta_{i}=x_{i}/z_{i}$, for $i=1,\dots,mn$. After obtaining $\Delta y$ by solving (7), we can compute $\Delta x$ and $\Delta z$ accordingly. To solve (7) efficiently, the following structure of $A\Theta A^{T}$ is useful.

As we can see, the iterates $x$ and $z$ generated by an IPM should be positive because of the third equation in (5). As a consequence, the matrix $V=\mathrm{Mat}(\theta)$ is dense and $A\Theta A^{T}$ contains two fully dense blocks. Therefore, factorizing $A\Theta A^{T}$ in (7) can be very expensive when $m+n$ is large. As a result, solving the normal system by an iterative solver is the only viable option. Moreover, it is obvious that $A\Theta A^{T}$ is always singular, since the vector $(e_{m};-e_{n})\in\mathbb{R}^{m+n}$ is an eigenvector for $A\Theta A^{T}$ corresponding to the zero eigenvalue. Note that the singularity of the coefficient matrix may not be an issue since one can remove the last row of $A$ in a prepossessing phase, which results in a normal equation having a coefficient matrix of the form

where $\hat{V}\in\mathbb{R}^{m\times(n-1)}$ is formed from the first $n-1$ columns of $V$. Then, by [35, Lemma 2.3], we see that the resulting coefficient matrix is nonsingular. However, the issue of ill-conditioning could still be a critical difficulty that goes against an iterative solver for fast convergence. In view of the above discussions, IPMs may not perform as efficiently as one would expect when solving large-scale OT problems.

In this subsection, we shall discuss the SSN method. We first revisit the concept of semismoothness which was originally introduced in [36] for functionals and later extended for vector-valued functions in [34]. Let $\mathbb{X}$ and $\mathbb{Y}$ be two finite dimensional Euclidean spaces and $F:\mathbb{X}\rightarrow\mathbb{Y}$ be a locally Lipschitz continuous function. According to Rademacher’s Theorem [37], $F$ is F-differentiable almost everywhere, and the Clarke’s generalized Jacobian [38] of $F$ at $x\in\mathbb{X}$, denoted by $\partial F(x)$, is well-defined and it holds that

where “$\mathrm{conv}$” denotes the convex hull operation and $D_{F}\subset\mathbb{X}$ is the set of points at which $F$ is F-differentiable.

If $F$ is directionally differentiable at $x\in\mathbb{X}$, i.e., the limit

exists for all $d\in\mathbb{X}$, and $\left\lVert F(x+d)-F(x)-Vd\right\rVert=o(\left\lVert d\right\rVert),\;\forall V\in\partial F(x+d),\;\forall d\rightarrow 0,$ we say that $F$ is semismooth at $x$. If $F$ is semismooth at $x\in\mathbb{X}$ and it holds that

then $F$ is said to be strongly semismooth at $x$. Moreover, $F$ is said to be (strongly) semismooth everywhere on $\mathbb{X}$, if $F$ is (strongly) semismooth at any point $x\in\mathbb{X}$.

Convex functions, smooth functions, and piecewise linear functions are examples of semismooth functions. In particular, one can verify that the projection operator onto the nonnegative orthant is strongly semismooth everywhere. Hence, the KKT mapping $F(x,y,z)$ of the following equation is also strongly semismooth everywhere:

where $\sigma>0$ is a given constant. Therefore, a natural idea is to apply the semismooth Newton (SSN) method [34], i.e., a nonsmooth version of Newton method, for solving (8). In particular, at each iteration of the SSN method, the following linear system is solved

where $(x,y,z)$ denotes the current iterate and $V=\mathrm{Diag}(v)\in\mathbb{R}^{mn\times mn}$ is an element in the Clarke’s generalized Jacobian (see [38]) of $\Pi_{+}(x-\sigma z)$. Specifically, for given $(x,z)$, we may choose

By the expression of $V$, we see that when $(x,z)$ are nearly optimal, $V$ could have a large number of zeros diagonal elements. Thus, the sparsity structure of the optimal solution can be utilized in the SSN method, which is an attractive feature when compared with the IPM. However, there are two issues that prevent us from applying the SSN method directly:

1. 1.

   The mapping $F$ is nonsmooth due to the nonsmoothness of $\Pi_{+}$. Hence, it is difficult to find a valid merit function for which a line search scheme for the SSN is well-defined. As a consequence, the globalization of the SSN method is a challenging difficulty that has yet to be resolved.

2. 2.

   It is not clear how to reduce the large-scale linear system (9), which is typically singular, to a certain normal linear system to save computational cost as in the case of the interior point method.

The key idea of a smoothing method [43] is to construct a smoothing approximation function $\mathcal{E}:\mathbb{R}_{++}\times\mathbb{R}^{mn}\times\mathbb{R}^{m+n}\rightarrow\mathbb{R}^{mn}\times\mathbb{R}^{m+n}$ for the mapping $E$, and then apply the Newton method for solving the smoothed system $\mathcal{E}(\epsilon,x,y)=0$, for $(\epsilon,x,y)\in\mathbb{R}\times\mathbb{R}^{mn}\times\mathbb{R}^{m+n}$. To this end, it is essential to construct a smoothing approximation of the projection operator $\Pi_{+}(\cdot)$, namely, $\Phi:\mathbb{R}_{++}\times\mathbb{R}^{mn}\rightarrow\mathbb{R}^{mn}$, such that for any $\epsilon>0$, $\Phi(\epsilon,\cdot)$ is continuously differentiable on $\mathbb{R}^{mn}$ and for any $x\in\mathbb{R}^{mn}$, it holds that $\left\lVert\Phi(\epsilon,x)-\Pi_{+}(x)\right\rVert\rightarrow 0$ as $\epsilon\downarrow 0$.

Note that a smoothing approximation function of $\Pi_{+}(\cdot)$ strongly depends on the smoothing function for the plus function

To the best of our knowledge, the Chen-Harker-Kanzow-Smale (CHKS) smoothing function [44, 45, 46], defined as

is the most commonly used smoothing function for $\pi(\cdot)$ in the literature. Readers are referred to [41] for other smoothing functions, whose properties are also well studied. It is easy to see that $\xi(\epsilon,t)$ maps any negative number $t$ to a positive value when $\epsilon\neq 0$, while the function $\pi(\cdot)$ maps any negative number $t$ to zero. Hence, the CHKS smoothing function does not preserve the sparsity structure of $\pi(\cdot)$. As a result, a smoothing Newton method developed based on the CHKS smoothing function is not able to exploit the sparsity in the solutions of the OT problem to reduce the cost of solving the linear system at each iteration of the algorithm.

To resolve the issue just mentioned, we propose to use the Huber smoothing function (a translation of the function considered in [47]) which also maps any negative number to zero. In particular, the Huber smoothing function we shall use in this paper is defined as:

Clearly, $h(\epsilon,t)$ is continuously differentiable for $\epsilon\neq 0$ and $t\in\mathbb{R}$. In fact, we see that for $\epsilon\neq 0$,

Moreover, it is easy to check that for any $t\in\mathbb{R}$,

and consequently,

where

Recall that the generalized Jacobian of the plus function $\pi$ has the following form

Therefore, we can see that for any $v\in\partial h(0,t)$, there exists a $u\in\partial\pi(t)$ such that $v(0,h)\equiv u(h)$ for any $h\in\mathbb{R}^{mn}$. One can also verify that the converse statement is true.

Using the Huber function $h(\cdot,\cdot)$, we can construct the smoothing approximation for the projection operator $\Pi_{+}(\cdot)$ as

Given the above smoothing function $\Phi(\cdot)$ for $\Pi_{+}(\cdot)$, a smoothing approximation of $E(\cdot,\cdot)$ can be constructed accordingly as follows:

where $\kappa_{p}>0$ and $\kappa_{c}>0$ are two given parameters. Note that adding the perturbation terms $\kappa_{p}\epsilon y$ and $\kappa_{c}\epsilon x$ is essential in our algorithmic design, since these perturbations will ensure that the Jacobian of $\mathcal{E}$ at $(\epsilon,x,y)$ is nonsingular when $\epsilon>0$.

In the current literature, smoothing Newton methods can be roughly divided into two groups, the Jacobian smoothing Newton methods and the squared smoothing Newton methods (see [48] for a comprehensive survey). The global convergence of a Jacobian smoothing Newton method strongly relies on the so-called Jacobian consistency property and the existence of a positive constant $\kappa$ such that

Unfortunately, the property (15) does not hold for (14). Therefore, in our setting, the Jacobian smoothing Newton method is not applicable. On the other hand, the squared smoothing Newton method tries to find a solution of $\mathcal{E}(\epsilon,x,y)=0$ by solving the following system

via the Newton method. As we shall see shortly, the global convergence of the squared smoothing Newton (SqSN) method does not rely on strong conditions such as (15). This explains why we choose the SqSN method in the present paper.

Associated with the system (16), we can define the merit function $\phi:\mathbb{R}\times\mathbb{R}^{mn}\times\mathbb{R}^{m+n}\rightarrow\mathbb{R}$:

to ensure that the line search procedure within the SqSN method is well-defined. We also define an auxiliary function $\zeta:\mathbb{R}\times\mathbb{R}^{mn}\times\mathbb{R}^{m+n}\rightarrow\mathbb{R}$:

which controls how the smoothing parameter $\epsilon$ is driven to zero, where $r\in(0,1)$ and $\tau\in(0,1]$ are two given constants. Using these notations, we can now describe our squared smoothing Newton method via the Huber function in Algorithm 1.

Before establishing the convergence properties of Algorithm 1, let us first see how to solve the linear system (17) and also explain the potential advantage of the SqSN method compared to the IPM and the SSN method. Suppose $\epsilon^{k}>0$, then, we see that

where $V_{1}^{k}:=\Phi_{1}^{\prime}(\epsilon^{k},w^{k})$ and $V_{2}^{k}:=\Phi_{2}^{\prime}(\epsilon^{k},w^{k})$ denote the partial derivatives of $\Phi$ with respect to the first and second arguments at the referenced point with $w^{k}:=x^{k}-\sigma(c-A^{T}y^{k})$, respectively. Note also that $\left\lvert(V_{1}^{k})_{ii}\right\rvert\in[0,1/2]$ and $\left\lvert(V_{2}^{k})_{ii}\right\rvert\in[0,1]$ for any $i=1,\dots,mn$. From (17), we get $\Delta\epsilon^{k}=-\epsilon^{k}+\zeta_{k}\epsilon^{0}$ and

where

It then follows that $\Delta x^{k}=\left[(1+\kappa_{c}\epsilon^{k})I_{mn}-V_{2}^{k}\right]^{-1}\left(r_{c}^{k}+\sigma V_{2}^{k}A^{T}\Delta y^{k}\right).$ Furthermore, we have that

Denote $v^{k}:=\mathrm{diag}\left[\left((1+\kappa_{c}\epsilon^{k})I_{mn}-V_{2}^{k}\right)^{-1}V_{2}^{k}\right]\in\mathbb{R}^{mn}$ and recall that $w^{k}:=x^{k}-\sigma(c-A^{T}y^{k})$. It is clear that

Keeping in mind that $z^{k}=c-A^{T}y^{k}$ for $k\geq 0$, we see that when $(x^{k},y^{k})$ is close to the optimal solution pair, $w^{k}$ is expected to have a large number of nonpositive elements and hence $v^{k}$ to be highly sparse. Moreover, let $V^{k}:=\mathrm{mat}(v^{k})\in\mathbb{R}^{m\times n}$. By Fact 1, we see that

When $V^{k}$ is sparse, it is clear that the coefficient matrix for computing $\Delta y^{k}$ is also sparse. This explains the appealing feature of Algorithm 1 as compared to the IPM, since the attractive sparsity structure of the optimal solution is exploited in Algorithm 1. Moreover, we are always able to reduce (17) to the smaller-scale normal system(18) in Algorithm 1, which is not the case in the SSN method. We can also see that the matrix $A$ is not explicitly required in the implementation of Algorithm 1, which is another advantage of the proposed method compared to an IPM. Lastly, let $J:=\mathrm{diag}(I_{m},-I_{n})$, then one can see that

which corresponds to the Laplacian matrix of a sparse bipartite graph defined by $V^{k}$ ²²2The bipartite graph has $m$ sources and $n$ sinks. There is an edge from source $i$ to sink $j$ if and only if $V^{k}_{ij}>0$ for $1\leq i\leq m$ and $1\leq j\leq n$.. First, one can identify all the connected components of the bipartite graph that are disjoint from each other. The matrix entries linking different components are always zero. Next, by permuting the coefficient matrix of the Newton equation, it can be converted into a block diagonal matrix. Finally, the Newton equation can be decomposed into several smaller equations. Each coefficient matrix of the smaller equation has a similar pattern as (20). In this way, the computational effort for computing the search direction may be further reduced.

For the rest of this section, we shall establish the convergence properties of Algorithm 1 which allow us to resolve the two difficulties of the SSN method as mentioned at the end of the last section. The tools for our analysis are adapted from those in the literature (see, e.g., [42]).

We first show that the coefficient matrix in (17) is nonsingular whenever $\epsilon^{k}>0$.

Lemma 1 shows that the linear system in (17) is always solvable, and the Newton direction $(\Delta^{k},\Delta x^{k},\Delta y^{k})$ is well-defined. Next, we shall show that the line search procedure in Algorithm 1 is well-defined, i.e., $\ell_{k}$ is always finite as long as $\epsilon^{k}$ is positive.