Solving Splitted Multi-Commodity Flow Problem by Efficient Linear Programming Algorithm

TL;DR

This paper introduces two novel algorithms, locSolver and incSolver, that significantly improve the efficiency of solving linear equations in column generation for multi-commodity flow problems by exploiting sparsity and solution reuse.

Contribution

The paper presents new algorithms for solving sparse linear systems more efficiently within column generation, reducing computational time in multi-commodity flow problem solving.

Findings

incSolver is at least 37 times faster than LAPACK.

Algorithms effectively exploit sparsity and solution similarity.

Preliminary experiments demonstrate substantial speedups.

Abstract

Column generation is often used to solve multi-commodity flow problems. A program for column generation always includes a module that solves a linear equation. In this paper, we address three major issues in solving linear problem during column generation procedure which are (1) how to employ the sparse property of the coefficient matrix; (2) how to reduce the size of the coefficient matrix; and (3) how to reuse the solution to a similar equation. To this end, we first analyze the sparse property of coefficient matrix of linear equations and find that the matrices occurring in iteration are very sparse. Then, we present an algorithm locSolver (for localized system solver) for linear equations with sparse coefficient matrices and right-hand-sides. This algorithm can reduce the number of variables. After that, we present the algorithm incSolver (for incremental system solver) which utilizes similarity in the iterations of the program for a linear equation system. All three techniques can be used in column generation of multi-commodity problems. Preliminary numerical experiments show that the incSolver is significantly faster than the existing algorithms. For example, random test cases show that incSolver is at least 37 times and up to 341 times faster than popular solver LAPACK.

Tables

Table 1. Table 1 : Different parts time comparison of incCG ,  kluCG  and lapackCG .

Case	Total time(s)			Shortest path computing time(s)			Linear equation solving time(s)
	incCG	kluCG	lapackCG	incCG	kluCG	lapackCG	incCG	kluCG	lapackCG
$R\left(1000\right)$	1538.51	1831.76	29438.10	178.65	177.50	188.56	721.23	1431.20	27392.50
$R\left(1500\right)$	625.76	801.76	18970.80	171.21	176.44	184.13	225.17	503.32	17690.60
$R\left(2000\right)$	258.90	301.54	5720.87	146.49	154.18	145.39	50.48	105.68	5219.57
$R\left(3000\right)$	180.84	187.72	2157.97	150.84	153.10	153.97	12.17	22.26	1888.38
$R\left(5000\right)$	158.91	161.67	609.58	151.81	152.48	152.80	2.14	4.56	425.30
$R\left(7000\right)$	200.37	201.93	631.31	194.18	194.16	198.87	1.61	3.34	404.67
$R\left(9000\right)$	219.11	217.10	606.23	213.44	210.16	231.69	1.31	2.73	350.03
$R\left(11000\right)$	246.31	251.99	622.97	240.89	245.40	265.44	1.12	2.30	334.39
$R\left(13000\right)$	266.82	278.51	618.77	261.25	271.96	287.63	1.02	2.03	309.15
$R\left(15000\right)$	270.44	271.01	491.10	265.30	265.19	277.14	0.74	1.46	197.66
$R\left(17000\right)$	315.24	314.90	590.03	309.22	308.13	324.42	0.84	1.65	246.35
$R\left(19000\right)$	370.37	375.79	706.79	363.69	368.36	372.57	0.91	1.73	312.02
$R\left(21000\right)$	366.81	374.80	604.03	360.39	367.75	350.00	0.76	1.46	235.44
$R\left(23000\right)$	393.72	395.79	600.08	387.01	388.48	357.91	0.72	1.37	223.79
$R\left(25000\right)$	383.42	386.20	544.56	377.12	379.47	386.96	0.55	1.05	142.73
$R\left(27000\right)$	459.25	460.63	690.73	451.89	452.79	460.94	0.68	1.24	211.07
$R\left(29000\right)$	484.89	484.76	701.85	477.21	476.70	488.44	0.64	1.16	195.01
$R\left(31000\right)$	537.48	539.50	738.52	529.59	531.18	524.11	0.63	1.14	195.86
$R\left(33000\right)$	560.34	588.66	758.42	552.50	580.30	563.80	0.58	1.05	177.01
$R\left(35000\right)$	618.65	623.70	817.98	610.31	614.62	605.78	0.61	1.14	193.48
$R\left(37000\right)$	692.31	701.09	886.41	683.28	691.64	644.37	0.64	1.11	221.87
$R\left(39000\right)$	677.80	699.56	812.32	668.49	689.65	609.01	0.60	1.13	184.44

Equations

$$\min$$

$$\sum _{i=1}^l{f}_i(u,v)\le q(u,v)(u,v)\in E,(\text{mbox}Capacityconstaints)$$

$$\sum _{v\in V}{f}_i(u,v)=\sum _{v\in V}{f}_i(v,u),u\in V\setminus \{s_i,t_i\},i=1,\cdots ,l(\text{mbox}Flowconservation)$$

$$\sum _{v\in V}{f}_i(s_i,v)=\sum _{v\in V}{f}_i(v,t_i)=d_i,i=1,\cdots ,l(\text{mbox}Demandsatisfication)$$

$$f_i(u,v)\ge 0,i=1,\cdots ,l\text{mbox}and(u,v)\in E$$

$$\sum _{v\in V}{f}_i(s_i,v)=\sum _{v\in V}{f}_i(v,t_i)\le d_i$$

$$E=[(1,3),(3,5),(5,7),(7,9),(2,4),(4,6),(6,8),(8,10),(3,4),(5,4),(6,5),(6,7),(7,9)]$$

$$\min$$

$$\sum _{i=1}^l{f}_i(u,v)\le q(u,v),$$

$$\sum _{v\in V}{f}_i(u,v)=\sum _{v\in V}{f}_i(v,u),$$

$$\sum _{v\in V}{f}_i(s_i,v)=\sum _{v\in V}{f}_i(v,t_i)\le d_i,$$

$$f_i(u,v)\ge 0,$$

$${\delta }_{p,e}=\{\begin{array}{ll}1& \text{ if link }e\text{ belongs to path }p\\ 0& \text{ otherwise.}\end{array}$$

$$\min$$

$$\sum _{p\in P_i}{\mathbf{x}}_p+y_i=d_i,$$

$$\sum _{i=1}^l\sum _{p\in P_i}{\mathbf{x}}_p{\delta }_{p,e}\le q(e),$$

$${\mathbf{x}}_p\ge 0,y_i\ge 0,$$

$$\begin{array}[]{lcl}{\mathbf{\beta}}_{e}[j]&=&\begin{cases}1&\mbox{ if }j=l+e^{\mbox{th}}\\ 0&\mbox{ otherwise}\end{cases}\end{array}$$

$$\begin{array}[]{lcl}{\mathbf{\beta}}_{p}[j]&=&\begin{cases}1&\mbox{ if }j=i\\ 1&\mbox{ if }j=l+e^{\mbox{th}}\mbox{ and }\delta_{p_{i},e}=1\\ 0&\mbox{ otherwise}\end{cases}\end{array}$$

$$\begin{array}[]{lcl}{\mathbf{\beta}}_{dummy_{i}}[j]&=&\begin{cases}1&\mbox{ if }j=i\\ 0&\mbox{ otherwise}\end{cases}\end{array}$$

$${\mathbf{A}}_k=[{\beta }_{p_{k,1}},\underset{p\in Q_{k,1}}{\underbrace{{\beta }_p}},\cdots ,{\beta }_{p_{k,l}},\underset{p\in Q_{k,l}}{\underbrace{{\beta }_p}},\underset{e\in N_k}{\underbrace{{\beta }_e}}]$$

$$\min$$

$${\mathbf{x}}_{p_i}+\sum _{p\in Q_i}{\mathbf{x}}_p=d_i,$$

$$\sum _{i=1}^l\left({\mathbf{x}}_{p_i}{\delta }_{p_i,e}+\sum _{p\in Q_i}{\mathbf{x}}_p{\delta }_{p,e}\right)=q(e),$$

$$\sum _{i=1}^l\left({\mathbf{x}}_{p_i}{\delta }_{p_i,e}+\sum _{p\in Q_i}{\mathbf{x}}_p{\delta }_{p,e}\right)+z_e=q(e),$$

$${\mathbf{x}}_p\ge 0,$$

$$z_e\ge 0,$$

$$\min$$

$$s.t.$$

$${\mathbf{A}}_k\mathbf{x}=\mathbf{b}$$

$$\mathbf{x}=[{\mathbf{x}}_{p_{k,1}},\underset{p\in Q_{k,1}}{\underbrace{{\mathbf{x}}_p}},\cdots ,{\mathbf{x}}_{p_{k,l}},\underset{p\in Q_{k,l}}{\underbrace{{\mathbf{x}}_p}},\underset{e\in N_k}{\underbrace{z_e}}]$$

$${\mathbf{A}}_k\lambda =\beta$$

$$j=\underset{i=1,{\lambda }_i>0}{\overset{l+|E|}{\mathrm{argmin}}}\frac{{\mathbf{x}}_i}{{\lambda }_i}.$$

$${\left({\mathbf{A}}_k\right)}^{⊺}\mu =-{\mathbf{c}}_k$$

$$O\left(h\left(l\left(|E|+|V|\log (|V|)\right)+{\left(l+|E|\right)}^{2.376}\right)\right).$$

$$\min$$

$$s.t.$$

$${\mathbf{A}}_k$$

$${\mathbf{c}}_k$$

$${\mathbf{b}}_k$$

Full text

1 Introduction

The multi-commodity flow problem (MCF) is a network flow problem with multiple commodities (flow demands) between different source and target nodes. Solving this problem is to find an assignment to all the flow variables such that certain given constraints are satisfied ford2004a . Many application problems can be reduced to MCF. Examples of these applications include the vehicle routing problem 　(VRP) letchford2015stronger ; cattaruzza2014an ; moshrefjavadi2016the , the traveling salesman problem (TSP) hernandezperez2009the , and problems of routing and wavelength assignment (RWA) leesutthipornchai2010solving ; patel2012routing . While it is well known that offline network resource optimization and planing in traditional network is a typical MCF Ahuja:1993 ; ford2004a ; Jajszczyk2005 , online network resource optimization and planing, which are now widely regarded as software defined network (SDN), are also treated as MCF hong2013achieving ; guo2014traffic ; kandula2015calendaring .

Because of its importance, there have been a sizable body of work on MCF, e.g. mahey2001capacity ; DBLP:GargK07 ; holmberg2003a ; Huisman2005 ; ZHU2012164 ; barnhart1998branch-and-price: ; degraeve2007a ; holmberg2003a ; huisman2007a ; salmasi2010total ; Briant2008 , in which colummn generation is widely used. A survey on column generation is given in lubbecke2005selected . There, the algorithms are divided in two classes, which are called exact algorithms dantzig1967generalized ; mccallum1977a ; barnhart1994a ; mamer2000a ; holmberg2003a ; dinh2013combining and approximation algorithms goldberg1992a ; grigoriadis1994fast ; awerbuch1994improved ; Grigoriadis96 ; even1999fast ; fleischer1999approximating ; fleischer2002fast ; bienstock2006approximating ; karakostas2008faster , respectively. In this paper, we will focus on the　 exact algorithm for the splitted multi-commodity flow problem in which the flow demands can be splitted among multiple paths for one commodity.

Organization: After this introduction, we introduce MCF and three different models for it in Section 2. In Section 3, we give a summary on column generation for MCF. We show in Section 4 how we apply the result in dantzig1967generalized to MCF, and present a concrete block structure of the basic matrix of column generation. In Section 5, we present the properties of the coefficient matrix. The test results show that the number of nonzero elements in each row of the coefficient matrix is less than 5 even when the length of the row is greater than $1000$. Thus, the matrix is very sparse. We devote Section 6 to present the two algorithms that are our main contribution in this paper. The first algorithm, called locSolver, is a localized system solver. This algorithm can reduce the number of variables in solving a linear equation when both its coefficient matrix and right-hand-sides are sparse. The second algorithm, called incSolver, is an incremental system solver which utilizes similarity during the iterations in solving linear equations. We present our experiment test results in Section 7, and conclusions in Section 8.

2 Model for Multi-Commodity Flow Problem

In this section, we define the basic formulation of multi-commodity flow problem (Model 1 below). We then present two more models, which are called Node-Link Formulation and Link-Path Formulation of MCF  respectively. Both are linear programming models with a large numbers of variables and constraints.

2.1 The Basic Model of MCF

Graphs are the most fundamental mathematical models for networks, and their edges and/or nodes are associated with numerical functions for quantity based network control and management. The basic graph model used to represent a MCF is a direct graph with capacities and weights assigned to its edges, which are used to represent factors and elements of “effectiveness” and “cost elements” of network resources, respectively.

A capacitated and weighted network is a triple $\mathcal{N}=\left(G(V,E),q,w\right)$, where

• •

$G(V,E)$ is a directed graph with the set $V$ of nodes (or vertices) and the set $E$ of links (or edges). A link $e\in E$ from node $u$ to node $v$ is denoted by $(u,v)$, where $u,v\in V$.

• •

$q$ and $w$ are mappings from $E$ to non-negative real numbers. For each edge $e\in E$, function $q$ assigns $e$ with a capacity $q(e)$, and function $w$ assigns $e$ with a weight $w(e)$, respectively.

A commodity is a measure of the demand in a network. Formally, for capacitated and weighted network $\mathcal{N}$, a commodity (or demand) is a triple $D=(s,t,d)$, where $s$ and $t$ are nodes of $\mathcal{N}$, and $d$ the bandwidth of non-negative value. The nodes $s$ and $t$ are called source and * target* of commodity $D$, respectively. We are now ready to formulate the basic model of MCF below.

Model 1** (MCF)**

Given a capacitated and weighed network $\mathcal{N}$, let $K=\{D_{1},D_{2},\cdots,D_{l}\}$ be a set of $l$ commodities, where $D_{i}=(s_{i},t_{i},d_{i})$ on $\mathcal{N}$, and $f_{i}(u,v)$ be a variable for each link $(u,v)$ of $\mathcal{N}$ that takes values in the interval $[0,d_{i}]$, for $i=1,\cdots,l$. The basic multi-commodity flow problem is to solve the following linear equation for the flow variables $f_{i}(u,v)$ with four constraints:

[TABLE]

Notice that constraint (1) is an objective function. The basic MCF formulation, the flow variables $f_{i}(u,v)$ of the commodities of $K$ represents the fraction of flow for commodity $D_{i}$ along edge $(u,v)$. Thus, $f_{i}(u,v)\in[0,d_{i}]$ in the general case when the commodity $d_{i}$ can be split among the flows of multiple paths, and $f_{i}(u,v)$ can only take one of the two possible valued $\{0,d_{i}\}$ otherwise (i.e. “single path routing”). In this paper, we focus on $f_{i}(u,v)\in[0,d_{i}]$. Taking the capacities and weights $q(u,v)$ and $w(u,v)$ of the edges $(u,v)\in E$ as the cost element, finding an assignment $f=(f_{1},\cdots,f_{l})$ in the above linear equation problem is called the minimum cost multi-commodity flow problem (min-MCF), indicated by constraint (1).

2.2 Node-Link Formulation

In Model 1, constraint (4) requires that the demand $d_{i}$ of each commodity is fully delivered through the flows along the paths from the source to the target. However, in general, only a part of the demand of a commodity can be “successfully” delivered, which means that constraints (4) become

[TABLE]

where $i=1,\cdots,l$.

Then it is desirable to seek the maximum portion of the command of each commodity to be successfully delivered with minimum cost. This case of MCF is called the maximal multi-commodity flow problem (MMCF). The primary requirement of MMC is to try to deliver all the demand, and the secondary requirement is to minimize the total cost.

We use $|S|$ to denote the cardinality of set $S$, and $|{\mathbf{A}}|$ to denote the dimension of a square matrix ${\mathbf{A}}$.

Model 2** (Node-Link Formulation Jajszczyk2005 )**

The formal description of MMCF  is defined as follows:

[TABLE]

where $W$ is a nonnegative real number that satisfies $W>\max\{\omega^{w}_{p}\mid p\in P_{i},\mbox{ for }i=1\cdots,l\}$ and $\omega^{w}_{p}=\sum_{(u,v)\in p}{w(u,v)}$.

$W\sum_{i=1}^{l}\left(d_{i}-\sum_{v\in V}f_{i}(s_{i},v)\right)$is the penalty term in the objective function. Node-Link Formulation is a linear programming model with $l|E|$ variables and $|E|+l(|V|-1)$ constraints. It is easy to see that MCF is a special case of MMCF  when $\sum_{i=1}^{l}\left(d_{i}-\sum_{v\in V}f_{i}(s_{i},v)\right)=0$, which means that all commodities are successfully delivered.

Example 1

In Fig. 1, we can choose $W$ as sum of all links’ weight, which is $34$.

2.3 Link-Path Formulation

In the previous two models of linear equations, the variables are the accounts of flows of links. We now present a formulation based on the accounts of flows of paths. For a path $p$, we denote the account of flow along path $p$ as a variable $\mathbf{x}_{p}$. For an arbitrary path $p$ and an edge $e$, we define the following (characteristic) function

[TABLE]

For a precise formulation of MMCF, we introduce the following notations below for a given set $K$ of commodity.

• •

Let $P_{i}$ denote an enumeration of the set of paths from $s_{i}$ to $t_{i}$ without loops (called simple paths), for $D_{i}=(s_{i},t_{i},d_{i})$ and $i=1,\cdots,l$.

• •

Given a path $p$, let $(u,v)\in p$ denote that edge $(u,v)$ is in path $p$ and path is along edge $(u,v)$.

Model 3** (Link-Path Formulation Jajszczyk2005 )**

MMCF* can be described as a problem of finding an assignment to the variables $\mathbf{x}_{p}$ for $p\in P_{i},\ i=1\cdots,l$, satisfying the following constraints.*

[TABLE]

In this model, (5) is the objective function, $y_{i}$ are slack variables which represent the portion of demand for commodity $D_{i}$ that fails to be delivered, and $W\sum_{i=1}^{l}y_{i}$ is the penalty term to objective function. Link-Path Formulation is a linear programming model with $l+|E|$ constraints and $\sum_{i=1}^{l}|P_{i}|+l$ variables. It is easy to see that $\sum_{i=1}^{l}|P_{i}|+l$ might become very large even for a small network.

Example 2

Fig. 2 is the topology (a bidirectional graph) of one backbone network of USA. This topology has 18 nodes and 52 links. Even in this small topology there are 97 different simple paths that connect Hawaii and Hartford.

In summary, we can see

both Node-Link Formulation and Link-Path Formulation are linear programming model. 2. 2.

Node-Link Formulation has fewer variables than Link-Path Formulation, while Link-Path Formulation has fewer constraints than Node-Link. 3. 3.

In general, both models either too many variables or too many constraints in practice.

3 The Column Generation Algorithm for Multi-Commodity Flow Problem

In this section, we first review the classical column generation. We then introduce a transition system model for understanding and analysis of this algorithm and the improved algorithm that we propose later. Finally in this section, we present the matrix formulation of classical column generation.

3.1 The Algorithm of Column Generation

The variables in Model 3 are often too many to be dealt with explicitly. Luckily, column generation ford2004a treats non-basic variables implicitly. It replaces the traditional method for determining a vector to entering basic by finding a shortest path which connects commodity source and target. It has better performance than the simplex method for Link-Path Formulation because both the number of variables and constraints are reduced to $|E|+l$ during every iteration. The basic idea of the algorithm is as follows.

In order to design an algorithm for full deliver of each demand $d_{i}$, we introduce a dummy path for for each commodity $D_{i}$, denoted by $dummy_{i}$. Let the capacity of $dummy_{i}$ be $d_{i}$ and $\omega^{w}_{dummy_{i}}=W$, where $W$ is value defined in Model 2. We call the original network extended with the dummy paths $dummy_{i},i=1,\cdots,l,$, the augmenting network, and define $\mathcal{P}=\bigcup_{i=1}^{l}\{P_{i}\cup\{dummy_{i}\}\}$ to denote the set of all paths of the commodities of the augmenting network.

The algorithm iteratively updates the load flow $\mathbf{x}_{p}$ for every path $p\in\mathcal{P}$, where $y_{i}=\mathbf{x}_{dummy_{i}}$ for the variables $y_{i}$ in Model 3, $i=1,\cdots,l$. . When it terminates, the values of path load flows $\mathbf{x}_{p}$ for all $p\in\mathcal{P}$ give an optimal solution for linear programming problem in Model 3.

Definition 1

Let $e^{\mbox{th}}$ is the index of edge $e$ in $E$. We introduce edge $e$’s basic vector ${\mathbf{\beta}}_{e}$ for $e\in E$, as follows:

[TABLE]

$e^{\mbox{th}}$* is the index of edge $e$ in $E$. In addition, we introduce path $p$’s basic vector ${\mathbf{\beta}}_{p}$ for $p\in P_{i},i=1,\cdots,l$, as follows:*

[TABLE]

We define ${\mathbf{\beta}}_{dummy_{i}}$ as follows:

[TABLE]

Example 3

In Fig. 1, let $e=(3,5)$, then ${\mathbf{\beta}}_{e}=[0,0,0,1,0,0,0,0,0,0,0,0,0,0,0]$. Let $p={1\rightarrow 3\rightarrow 5\rightarrow 7\rightarrow 9}\in P_{1}$, then ${\mathbf{\beta}}_{p}=[1,0,1,1,1,1,0,0,0,0,0,0,0,0,0]$.

For an assignment $\mathbf{x}_{p}$ of $p\in P_{i}$ and $i=1,\cdots,l$, the value $\min\{q(e)-\sum_{i=1}^{l}\sum_{p_{1}\in P_{i}}\delta_{p_{1},e}\mathbf{x}_{p_{1}}\mid\delta_{p,e}=1\}$ is called the remaining capacity of $p$, denoted by $\textit{RemainCapacity}(p)$. We say that $p$’s remaining capacity carries commodity $D_{i}$ if and only if $p\in P_{i}$ and $d_{i}$ is less or equal to $\textit{RemainCapacity}(p)$.

3.2 Transition System Model

To help the understanding and analysis of Algorithm SCG, we introduce a state transition system that models the state change by each iteration of the main loop of the algorithm, i.e. lines 1 - 1 of Algorithm SCG. To define the abstract states of the transition system, we need the invariant property of the algorithm in the following lemma.

Lemma 1

Constraint (6) in Model 3 is an invariant of the main loop in Algorithm SCG (lines 1 - 1).

Proof

The lemma holds because of the fact that the values of the variables $\{\mathbf{x}_{p}\leavevmode\nobreak\ |\leavevmode\nobreak\ p\in\mathcal{P}\}$ are alway kept in their feasible area is an invariant of the simplex method.∎

Since $d_{i}>0$, Lemma 1 implies that for each iteration, say the $k$th iteration, there is at least one $p\in(P_{i}\cup\{dummy_{i}\})$ for $i=1,\cdots,l$ such that $\mathbf{x}_{p}>0$. For the $k$th iteration and commodity $D_{i}$, a path $p_{k,i}\in P_{i}\cup\{dummy_{i}\}$ which has positive flow can be selected as the primary path and the subset $Q_{k,i}\subseteq P_{i}\setminus\{p_{k,i}\}$ of paths which have non-zero flow as the secondary paths of $D_{i}$, where $k,i=1,\ldots,l$.

We now describe the main loop of Algorithm SCG as the transition system such that the $k$th iteration changes from a state of the form $([(p_{k,1},Q_{k,1}),\cdots,(p_{k,l},Q_{k,l})],N_{k})$ to a state $([(p_{k+1,1},Q_{k+1,1}),\cdots,(p_{k+1,l},Q_{k+1,l})],N_{k+1})$ where $N_{k},N_{k+1}\subseteq E$.

After initial solution steps in Algorithm SCG  (line 1 to line 1), the system state is $p_{1,i}=p_{i},\ Q_{1,i}=\emptyset$ for $i=1,\cdots,l$ and $N_{1}=E$. The transition rules are defined in the following way.

When the entering variable is a link $e^{*}$ :

1. (a)

When the leaving variable is a path $p_{k,i}$:

By Lemma 1, there is a $p\in Q_{k,i}$. Let $p_{k+1,i}=p,$ $Q_{k+1,i}=Q_{k,i}\setminus\{p\}$ and $N_{k+1}=N_{k}\cup\{e^{*}\}$, the other $(k+1)$th’s state are the same as $k$th’s state. In the following description, without loss of generality we do not mention the unchanged state part. 2. (b)

When the leaving variable is a path $p\in Q_{k,i}$:

Let $Q_{k+1,i}=Q_{k,i}\setminus\{p\},$ $N_{k+1}=N_{k}\cup\{e^{*}\}$. 3. (c)

When the leaving variable is a link $e$:

Let $N_{k+1}=(N_{k}\cup\{e^{*}\})\setminus\{e\}$. 2. 2.

When the entering variable is a path $p^{\prime}_{j}$:

1. (a)

When the leaving variable is a path $p_{k,i}$:

1. i.

When $i=j$:

Let $p_{k+1,i}=p^{\prime}_{j}$. 2. ii.

When $i\neq j$:

By Lemma 1, there is a $p\in Q_{k,i}$.

Let $p_{k+1,i}=p,$ $Q_{,k+1,i}=Q_{k,i}\setminus\{p\},$ $Q_{k+1,j}=Q_{k,j}\cup\{p_{j}^{\prime}\}$. 2. (b)

When the leaving variable is a path $p\in Q_{k,i}$:

1. i.

When $i=j$:

Let $Q_{k+1,i}=(Q_{k,i}\cup\{p_{j}^{\prime}\})\setminus\{p\}$. 2. ii.

When $i\neq j$:

Let $Q_{k+1,j}=Q_{k,j}\cup\{p_{j}^{\prime}\},$ $Q_{k+1,i}=Q_{k,i}\setminus\{p\}$. 3. (c)

When the leaving variable is a link $e$:

Let $Q_{k+1,j}=Q_{k,j}\cup\{p_{j}^{\prime}\},$ $N_{k+1}=N_{k}\setminus\{e\}$.

It is easy to see that state $([(p_{k,1},Q_{k,1}),\cdots,(p_{k,l},Q_{k,l})],N_{k})$ represents the basic matrix ${\mathbf{A}}_{k}$ in the $k$th iteration, where

[TABLE]

In other words, ${\mathbf{A}}_{k}$ is the incidence matrix of paths $p_{k,1},Q_{k,1},\cdots,p_{k,l},Q_{k,l}$ and edges in $N_{k}$

Definition 2

In the above transition system, if the entering variable is a path $p$ and the leaving variable is a link $e$, then we call $p$ a basic variable which corresponds to $e$, denoted as $p_{e}$.

The variable $p_{e}$ has some update rules. When the entering variable is a path $p$ and the leaving variable is a path $p_{e}$, then we update $p_{e}=p$. If the entering variable is a link $e_{1}$ and leaving variable is a path $p_{e}$, then update $p_{e}=p_{e_{1}}$.

Note 1

Let $\SS_{k}=E\setminus N_{k},\ Q_{k}=\bigcup_{i=1}^{l}Q_{k,i}.$ We call $\SS_{k}$ is the set of saturated link.

The intuitive meaning of a saturated link is that its bandwidth has been fully taken up and its bandwidth restricts the objective function to further decrease under current basis.

Lemma 2

$|Q_{k}|=|\SS_{k}|$* is an invariant of the main loop, in other words, there is a path $p_{e}\in Q_{k}$ for each $e\in\SS_{k}$.*

Proof

We prove it by induction. When $k=0$, it obviously holds. Assume that the conclusion holds when $k\leq K_{1}$.

When $k=K_{1}+1$, if rules 1-(a) and 1-(b) are used in Fig. 3, then both cardinal of $Q_{k}$ and $\SS_{k}$ decrease by $1$ compared with last iteration. Hence conclusion holds. If rules 1-(c), 2-(a)-i and 2-(b) are used, then both cardinal of $Q_{k}$ and $\SS_{k}$ are unchanged. Hence conclusion holds. If rule 2-(a) is used, then both cardinal of $Q_{k}$ and $\SS_{k}$ increase by $1$. Thus, conclusion holds. In summary, no matter what rule is used in $k$-th, $|Q_{k}|=|\SS_{k}|$ holds. ∎

3.3 Matrix Formulation

We fix working paths on $p_{1},\cdots,p_{l},Q_{1},\cdots,Q_{l}$ and add slack variables $z_{e}$ for constraint (7) where $e\in N_{k}$, then MMCF can be described as follows:

Model 4** (Link-Path Formulation for augmenting network)**

[TABLE]

Note 2

$\mathbf{c}_{k}=[w_{p_{k,1}},\underbrace{w_{p}}_{p\in Q_{k,1}},\cdots,w_{p_{k,l}},\underbrace{w_{p}}_{p\in Q_{k,l}},\underbrace{0,\cdots,0}_{|N_{k}|}].$ ${\mathbf{b}}=[d_{1},\cdots,d_{l},\underbrace{q(e)}_{e\in E}]$.

Example 4

In Fig. 1, ${\mathbf{b}}=[10,11,10,10,10,10,10,15,8,10,10,7,10,5,10]$

In other words, Model 4 can be written as

Model 5** (Matrix formulation )**

[TABLE]

where ${\mathbf{A}}_{k}$ is defined as (8).

In $k$th iteration the leaving variable selection procedure in Algorithm SCG 　 can been described as follows:

[TABLE]

Firstly, solve equation (15), and obtain

[TABLE]

which satisfy constraints (10)-(14).

Secondly, solve equation (16), and obtain ${\mathbf{\lambda}}$.

[TABLE]

where vector ${\mathbf{\beta}}$ is the basic vector corresponding to entering variable.

Finally, choose a leaving variable by a pivot rule.

Note 3 (pivot rule)

There are many different mays to choose leaving variable. In this paper, we apply classical pivot rule to pick $j$th as leaving variable where

[TABLE]

The solution ${\mathbf{\mu}}$ of equation (17) are dual values of constraints (10) and (12). Then let $w^{\prime}=w+{\mathbf{\mu}}$ be newly updated link weights.

[TABLE]

3.4 Classical Column Generation Complexity Analysis

Suppose Algorithm SCG does $h$ main iterations before termination, then Algorithm SCG computes $(h+1)l$ shortest path and solves $3h$ linear equation systems in the form of (15), (16) and (17) where ${\mathbf{A}}_{k}$’s size is $(l+|E|)\times(l+|E|)$. As the authors know that the best shortest path algorithm complexity is $O\left(|E|+|V|\log(|V|)\right)$ which is given by Dijkstra’s algorithm based on Fibonacci heap and the best linear system solving algorithm complexity is $O\left(\left(l+|E|\right)^{2.376}\right)$. Hence, the Algorithm SCG’s complexity is

[TABLE]

4 Speedup Through Employing ${\mathbf{A}}_{k}$’s Structure

The complexity of (18) can not be accepted in reasonable time when the size of $G(V,E)$ is large. This hinders SCG’s use in some applications e.g. online load balance in SDN and large scale problem offline. Hence, how to improve the efficiency of column generation is a problem considered in many works dantzig1967generalized ; mccallum1977a ; barnhart1994a ; mamer2000a . The complexity (18) only has two parts, i.e. computing shortest path and solving linear equation systems (15), (16) and (17). 　Hence, reducing coefficient matrix size is a feasible approach. Luckily, the primal partitioning procedure, a specialization of the generalized upper bounding procedure developed by Dantzig and Van Slyke dantzig1967generalized , involves the determination at each iteration of the inverse of a basis containing only one row for each saturated link. In other words, we can reduce matrix size to the number of saturated link. In the following, we will concretely show how to apply conclusion of dantzig1967generalized on MCF. Through reordering column of basis matrix to obtain a special structure in resulted basis matrix ${\mathbf{A}}_{k}$, we give bellow a method called structured matrix method (SMCG). By this way we can reduce the size of linear equation to be solved in general.

4.1 Structured Matrix Method for Column Generation

After we reorder basic variable in $k$th iteration by $p_{k,1},\cdots,p_{k,l},\underbrace{p}_{p\in Q_{k,1}},\cdots,\underbrace{p}_{p\in Q_{k,l}},\underbrace{e}_{e\in N_{k}}$, Model 4 can be rewritten as:

Model 6** (Structure matrix model )**

[TABLE]

where

[TABLE]

Mathematically, we can rewrite equation (15) into:

[TABLE]

Then we have

[TABLE]

Hence, we can firstly solve equation (23). Secondly, substituting $\mathbf{x}_{\SS_{k}}$ in (20) to obtain $\mathbf{x}_{K}$. Finally, substituting $\mathbf{x}_{K},\mathbf{x}_{\SS_{k}}$ in (22) to obtain $\mathbf{x}_{N_{k}}$. In this way can solve equation system (15).

Note 4

Let

[TABLE]

Now equation (23) can be written as:

[TABLE]

Lemma 3

${\mathbf{M}}_{k}$* is a non-singular sparse matrix.*

Proof

By simplex method theory, ${\mathbf{A}}_{k}$ is a non-singular matrix. By structure of ${\mathbf{A}}_{k}$ in (19), $\det\left({\mathbf{A}}_{k}\right)=\det\left({\mathbf{M}}_{k}\right)\neq 0$. Therefore, ${\mathbf{M}}_{k}$ is a non-singular matrix.∎

Since equations (15) and (16) have the same coefficient matrix. Hence, employing the way to solve equation (15), we can solve equation (16). Through the same method we obtain

[TABLE]

${\mathbf{\lambda}}_{K}={\mathbf{\beta}}_{K}-{\mathbf{B}}_{k}{\mathbf{\lambda}}_{\SS_{k}}\mbox{ and }{\mathbf{\lambda}}_{N_{k}}={\mathbf{\beta}}_{N_{k}}-{\mathbf{H}}_{k}{\mathbf{\lambda}}_{K}-{\mathbf{F}}_{k}{\mathbf{\lambda}}_{\SS_{k}}.$

In equation (17), mathematically, $\left({\mathbf{A}}_{k}\right)^{T}{\mathbf{\mu}}=\mathbf{c}_{k}$ can be rewritten as

[TABLE]

We substitute ${\mathbf{\mu}}_{N_{k}}={\mathbf{0}}_{N_{k}}$ in (27) and (28) to obtain

[TABLE]

We simplify equation systems (29) and (30) by

[TABLE]

[TABLE]

Through the above simplification, we can solve equation (17) too.

4.2 Structure Matrix Method’s Complexity Analysis

Suppose SCG does $h$ main iterations before termination, then it computes $(h+1)l$ shortest path and in $k$th iteration we need to solve $3$ linear equation systems in the form of (25), (26) and (31) where ${\mathbf{M}}_{k}$’s size is $|\SS_{k}|\times|\SS_{k}|$. By the same discussion in Section 3.4, we can obtain that the Structure Matrix Method complexity is

[TABLE]

As given by the analysis in Section 3.4, the standard column generation method is

[TABLE]

by (18). By the Note 1, it is easy to see that $|\SS_{k}|<l+|E|$. And in general, $|\SS_{k}|<|E|$, hence SMCG  is better than the classical one (SCG).

5 Speedup Through Employing ${\mathbf{M}}_{k}$’s Sparse Structure

The sparse property is very useful when solving linear equation, in the following, to show the sparse property of matrix ${\mathbf{M}}_{k}$ we will discuss the element of matrix ${\mathbf{M}}_{k}$ in detail. In ${\mathbf{M}}_{k}={\mathbf{C}}_{k}{\mathbf{B}}_{k}-{\mathbf{D}}_{k}$, ${\mathbf{C}}_{k}[i]$ is a vector denoting whether path $p_{k,i}$ crosses each edge in $\SS_{k}$, i.e.

[TABLE]

${\mathbf{B}}_{k}[i]$ is a vector associated with $i$th path of $[Q_{k,1},\cdots,Q_{k,l}]$ and its value indicates which commodity the $i$th path belong to. Let $i$th path of $[Q_{k,1},\cdots,Q_{k,l}]$ be path for commodity $D_{h}$. Then

[TABLE]

Hence, ${\mathbf{C}}_{k}{\mathbf{B}}_{k}[i]={\mathbf{C}}_{k}[h]$, where ${\mathbf{C}}_{k}[h]$ is associated with path $p_{k,h}$ and

[TABLE]

${\mathbf{D}}_{k}[i]$　 is also a vector associated with $i$th path of $[Q_{k,1},\cdots,Q_{k,l}]$ and its value indicates which edge it crosses. Let path $p$ be the $i$th path of $[Q_{k,1},\cdots,Q_{k,l}]$. Then,

[TABLE]

By conclusion of fronczak2004average the ratio between path length and $|E|$ is very small when $|E|$ is large in general. When we see the graph as only consisted of saturated links, then number of nonzero elements in vector ${\mathbf{C}}_{k}{\mathbf{B}}_{k}[i]$ and ${\mathbf{D}}_{k}[i]$ are identical with the length of associated paths $p_{k,h}$ and $p$. Therefore, ${\mathbf{C}}_{k}{\mathbf{B}}_{k}[i]$ and ${\mathbf{D}}_{k}[i]$ are two sparse vectors. As discussed above, ${\mathbf{M}}_{k}$ is statistically a very sparse matrix.

In the following, we list some experiment results of matrix ${\mathbf{M}}_{k}$. We record the dimension and number of ${\mathbf{M}}_{k}$’s nonzero elements in every iteration for some random cases. Let $N({\mathbf{M}}_{k})$ be the number of nonzero element in matrix ${\mathbf{M}}_{k}$. The detail of cases’ configuration can be found in section 7. The dimension of ${\mathbf{M}}_{k}$ is equal to number of saturated link in $k$th iteration. In Fig. 4 we can see that the dimension starts from [math] to a large number (more than $1000$), which indicates that the number of saturated links is more and more larger as iteration proceeding, and the resource competition of different commodities is more and more intense. But the growth of the ratio between nonzero coefficients of matrix ${\mathbf{M}}_{k}$ and its dimension is very slow. In Fig. 4(a) when $k>500000$, the value $\frac{N({\mathbf{M}}_{k})}{|{\mathbf{M}}_{k}|}$ is still less than $5$ while $|{\mathbf{M}}_{k}|$ is larger than $1000$. Hence ${\mathbf{M}}_{k}$ is a very sparse matrix.

According to vanderbei1998linear ’s suggestions, we use LU decomposition to solve equations (25), (26) and (31). But because of the high sparsity of matrix ${\mathbf{M}}_{k}$, LAPACK  anderson1990lapack kernels are not applicable. Hence, we can use the linear solver KLU davis2010algorithm , which has high performance for sparse matrix, to solve equations (25), (26) and (31) instead of LAPACK.

In Table 1 we find an interesting phenomenon that when the value $\frac{N({\mathbf{M}}_{k})}{|{\mathbf{M}}_{k}|}$ is greater than or equal to $3$ (case $R(1000)$), solving linear equation will become dominating part of total time consumption. And when $\frac{N({\mathbf{M}}_{k})}{|{\mathbf{M}}_{k}|}$ is less than $3$, sparse linear solvers will greatly reduce the time of linear equation solving.

6 Speedup Through Sparse and Similar Properties

After we employ KLU to solve linear equations occurring in iteration of SMCG, when the number of saturated link is small, then linear equation solving step is not a dominating part of time. But when there are many saturated links the complexity of SMCG  is almost the same as classical one. When structure of matrix ${\mathbf{M}}$ is complex (i.e. nonzero elements of matrix ${\mathbf{M}}_{k}$ is more than $3|{\mathbf{M}}_{k}|$), then linear equation solving step dominates the entire algorithm time, even employing KLU. Therefore, in the following, we do not invoke KLU to solve equations 　but directly use results in previous iteration to incrementally solve equations (15), (16) and (17).

For keeping speedup, in this section we firstly give a fast method locSolver which reduces Problem 1 to a small one in Section 6.1. Secondly, we provide an incremental method incSolver to solve equations (15), (16) and (17) in Section 6.2. Thirdly, in the final, we discuss why locSolver and incSolver are proper solvers for equation during iteration in Section 6.4.

6.1 A Fast Method to Solve Sparse Linear Equation System

Problem 1

Solve linear equation system

[TABLE]

where ${\mathbf{A}}$ is a $n\times n$ matrix and ${\mathbf{b}}$ is a vector.

We will provide a fast algorithm to solve Problem 1. This method can reduce linear equation system to a small one. Especially, when ${\mathbf{A}}$ is a very sparse matrix and ${\mathbf{b}}$ is also a sparse vector, this method is very powerful. For presenting this fast method we first give following definitions and lemmas.

Definition 3

For $n\times n$ matrix ${\mathbf{A}}$, let $G({\mathbf{A}})$ be the undirected graph of matrix ${\mathbf{A}}$ with $2n$ nodes. $G({\mathbf{A}})$ has a link $(i,n+j)$ iff ${\mathbf{A}}[i,j]\neq 0$. Let ${\tt reach}_{G({\mathbf{A}})}(B)$ be the set of nodes reachable from element of $B$ through $G({\mathbf{A}})$.

Lemma 4

In Problem 1, let $B=\{i\mid{\mathbf{b}}[i]\neq 0\},I=\{i\mid i\in{\tt reach}_{G({\mathbf{A}})}(B),i<n\}$. If $B\not\subseteq I$, then Problem 1 has no solution.

Proof

Set $h\in B,h\not\in I$. Thus, all the elements in row $h$ are zeros. Therefore, we have ${\mathbf{A}}[h,1]{\mathbf{\alpha}}[1]+\cdots+{\mathbf{A}}[h,n]{\mathbf{\alpha}}[n]=0{\mathbf{\alpha}}[1]+\cdots+0{\mathbf{\alpha}}[n]=0$ for every ${\mathbf{\alpha}}\in\mathbb{R}^{n}$. But ${\mathbf{b}}[h]\neq 0$, so there is no ${\mathbf{\alpha}}\in\mathbb{R}^{n}$ satisfying Problem 1. ∎

Definition 4

$I,J$* are two sub-sequences of $1,2,\cdots,n$. ${{\mathbf{A}}_{I,J}}$ is called $(I,J)$-projection of ${\mathbf{A}}$ (briefly projection of ${\mathbf{A}}$) if ${{\mathbf{A}}_{I,J}}[i,j]={\mathbf{A}}[I[i],J[j]],$ for $i=1,\cdots,|I|,j=1,\cdots,|J|$*

Definition 5

In Problem 1, let $B=\{i\mid{\mathbf{b}}[i]\neq 0\}$. ${{\mathbf{A}}_{I,J}}$ is called a computable projection of Problem 1 if

1. (i)

$B\subseteq I$; 2. (ii)

$\{i\mid{\mathbf{A}}[i,j]\neq 0,j\in J\}\subseteq I$. 3. (iii)

$\{j\mid{\mathbf{A}}[i,j]\neq 0,i\in I\}\subseteq J$.

Definition 6

$I$* is a sub-sequence of $1,2,\cdots,n$. ${\mathbf{b}}_{I}$ is called an $I$-projection of ${\mathbf{b}}$ if ${\mathbf{b}}_{I}[i]={\mathbf{b}}[I[i]]$ for $i=1,\cdots,|I|$*

Definition 7

$I$* is a sub-sequence of $1,\cdots,n$ and ${\mathbf{\alpha}}$ is a vector such that $|{\mathbf{\alpha}}|=|I|$. ${\tt lift}({\mathbf{\alpha}},I,n)$ denotes a lifting vector where*

[TABLE]

Lemma 5

Let ${{\mathbf{A}}_{I,J}}$ be a computable projection of Problem 1. If there exists a vector ${\mathbf{\alpha}}$ satisfying that ${{\mathbf{A}}_{I,J}}{\mathbf{\alpha}}={\mathbf{b}}_{I}$, then ${\tt lift}({\mathbf{\alpha}},I,n)$ is a solution of Problem 1.

Proof

Let ${\mathbf{\xi}}={\tt lift}({\mathbf{\alpha}},I,n)$. In the following we want to prove that ${\mathbf{A}}{\mathbf{\xi}}={\mathbf{b}}$. W.l.o.g. set $I=\{1,2,\cdots,n_{1}\},\ J=\{1,2,\cdots,m_{1}\}$.

First, we will prove that ${\mathbf{A}}[i]{\mathbf{\xi}}={\mathbf{b}}[i]$ for $i=1,\cdots,n_{1}$. By definition of ${\mathbf{\xi}}$, ${\mathbf{\xi}}[i]=0$ for $i>m_{1}$. Thus,

[TABLE]

for $i=1,\cdots,n_{1}$.

Second, we will prove that ${\mathbf{A}}[i]{\mathbf{\xi}}={\mathbf{b}}[i]=0$ for $i=n_{1}+1,\cdots,n$. Since ${{\mathbf{A}}_{I,J}}$ is a computable projection, ${\mathbf{A}}[i,j]=0$ for $i>n_{1},j\leq m_{1}$. Thus,

[TABLE]

for $i=n_{1}+1,\cdots,n$. By the definition of ${\tt lift}$, ${\mathbf{\xi}}[i]=0$ for $i>m_{1}$. Therefore,

[TABLE]

for $i=n_{1}+1,\cdots,n$.

In summary, ${\mathbf{A}}[i]{\mathbf{\xi}}={\mathbf{b}}[i]$ for $i=1,\cdots,n$. So ${\tt lift}({\mathbf{\alpha}},I,n)$ is a solution to Problem 1. ∎

Theorem 6.1

If ${\mathbf{A}}_{I,J}$ is a computable projection of Problem 1, then the system has solution ${\mathbf{\xi}}$ iff there exists a vector ${\mathbf{\alpha}}$ such that ${{\mathbf{A}}_{I,J}}{\mathbf{\alpha}}={\mathbf{b}}_{I}$. Furthermore, ${\tt lift}({\mathbf{\alpha}},I,n)$ is a solution of Problem 1.

Proof

W.l.o.g. we set $I=\{1,2,\cdots,n_{1}\},\ J=\{1,2,\cdots,m_{1}\}$. Let us assume that ${\mathbf{\xi}}$ is a solution of Problem 1. By condition (iii) of Definition 5, we have $A[i,j]=0$ when $i\leq n_{1}$ and $j>m_{1}$. Thus,

[TABLE]

for $1,\cdots,n_{1}$. In other words, ${{\mathbf{A}}_{I,J}}{\mathbf{\xi}}_{J}={\mathbf{b}}_{I}$.

Since ${{\mathbf{A}}_{I,J}}$ is one of computable projections of Problem 1, then applying Lemma 5, ${\tt lift}({\mathbf{\alpha}},J,n)$ is one solution of Problem 1.∎

Lemma 6

In Problem 1, let $B=\{i\mid{\mathbf{b}}[i]\neq 0\},\ I=\{i\mid i\in{\tt reach}_{G({\mathbf{A}})}(B),\ i<n\},\ J=\{j-n\mid j\in{\tt reach}_{G({\mathbf{A}})}(B),\ j\geq n\}$. Then ${{\mathbf{A}}_{I,J}}$ is a computable projection of Problem 1, if it has a solution.

Proof

Since Problem 1 has a solution, we have $B\subseteq I$ because otherwise it will conflict with Lemma 4. By definition of $I$ and $J$, $I=\{i\mid{\mathbf{A}}[i,j]\neq 0,j\in J\}$ and $J=\{j\mid{\mathbf{A}}[i,j]\neq 0,i\in I\}$. In summary, $I,J$ satisfy condition (i) (ii) (iii) of Definition 5. So, ${{\mathbf{A}}_{I,J}}$ is a computable projection of Problem 1.∎

Corollary 1

In Problem 1, let $B=\{i\mid{\mathbf{b}}[i]\neq 0\},\ I=\{i\mid i\in{\tt reach}_{G({\mathbf{A}})}(B),\ i<n\},\ J=\{j-n\mid j\in{\tt reach}_{G({\mathbf{A}})}(B),\ j\geq n\}$. Then ${{\mathbf{A}}_{I,J}}\mathbf{y}={\mathbf{b}}_{I}$ has solution and ${\tt lift}(\mathbf{y},I,n)$ is a solution to Problem 1, if Problem 1 is feasible.

Proof

When there is an $\mathbf{x}$ satisfying Problem 1, employing Lemma 6, ${{\mathbf{A}}_{I,J}}$ is a computable projection of ${\mathbf{A}}$. By Theorem 6.1, the conclusion holds.∎

Theorem 6.2

In Problem 1, let $B=\{i\mid{\mathbf{b}}[i]\neq 0\},\ I=\{i\mid i\in{\tt reach}_{G({\mathbf{A}})}(B),i<n\},\ J=\{j-n\mid j\in{\tt reach}_{G({\mathbf{A}})}(B),\ j\geq n\}$. If ${\mathbf{A}}$ is a non-singular matrix then there is a unique vector $\mathbf{y}$ such that ${{\mathbf{A}}_{I,J}}\mathbf{y}={\mathbf{b}}_{I}$.

Proof

Since ${\mathbf{A}}$ is a non-singular matrix, Problem 1 has a unique solution $\mathbf{x}$. Employing Corollary 1, there is a vector $\mathbf{y}$ such that ${{\mathbf{A}}_{I,J}}\mathbf{y}={\mathbf{b}}_{I}$. Suppose that there is another $\mathbf{y}^{\prime}\neq\mathbf{y}$ such that ${{\mathbf{A}}_{I,J}}\mathbf{y}^{\prime}={\mathbf{b}}_{I}$. By Lemma 5, $\mathbf{x}={\tt lift}(\mathbf{y},J,n),\mathbf{x}^{\prime}={\tt lift}(\mathbf{y}^{\prime},J,n)$ are two solutions of Problem 1. In other words, ${\mathbf{A}}(\mathbf{x}-\mathbf{x}^{\prime})=0$. It is easy to check that $\mathbf{x}\neq\mathbf{x}^{\prime}$. Hence, ${\mathbf{A}}$ is a singular matrix, which conflicts with the fact that ${\mathbf{A}}$ is a non-singular matrix. So, ${{\mathbf{A}}_{I,J}}\mathbf{y}={\mathbf{b}}_{I}$ has a unique solution. ∎

Corollary 2

In Problem 1, let $B=\{i\mid{\mathbf{b}}[i]\neq 0\},\ I=\{i\mid i\in{\tt reach}_{G({\mathbf{A}})}(B),\ i<n\},\ J=\{j-n\mid j\in{\tt reach}_{G({\mathbf{A}})}(B),\ j\geq n\}$. If ${\mathbf{A}}$ is a non-singular matrix then $|I|=|J|$ and ${{\mathbf{A}}_{I,J}}$ is a non-singular matrix.

Proof

Let ${\mathbf{S}}=\left({{\mathbf{A}}_{I,J}}\right)^{T}\left({{\mathbf{A}}_{I,J}}\right)$. It is easy to know that ${\mathbf{S}}$ is a symmetric matrix. Proving ${\mathbf{S}}$ is a non-singular matrix is equivalent to check that $\mathbf{x}={\mathbf{0}}$ is a unique solution of $\mathbf{x}^{T}{\mathbf{S}}\mathbf{x}=0$. By Theorem 6.2, ${\mathbf{0}}$ is a unique solution of equation system ${{\mathbf{A}}_{I,J}}\mathbf{x}={\mathbf{0}}$. Thus, $\mathbf{x}={\mathbf{0}}$ is a unique solution of $\mathbf{x}^{T}{\mathbf{S}}\mathbf{x}=0$. In other words, ${\mathbf{S}}$ is a non-singular matrix and $|I|\geq|J|$ and ${\tt rank}({{\mathbf{A}}_{I,J}})=|J|$.

Let $b^{\prime}$ be a vector which satisfies that $\{i\mid b^{\prime}[i]\neq 0\}=J$. When solving equation ${\mathbf{A}}^{T}\mathbf{x}=b^{\prime}$, by the same way of above discussion we can obtain that $|J|\geq|I|$. Therefore, $|I|=|J|$ and ${{\mathbf{A}}_{I,J}}$ is a non-singular matrix. ∎

6.1.1 Impoving algorithm locSolver during iteration

Computing reachable edges for a given node set $B$ is a key step of locSolver. Although computing reachable set of a given graph is of linear complexity, but in this case we need to construct a new graph per iteration. Luckily, by Theorem 6.1, we can use any computable projection to replace ${{\mathbf{A}}_{I,J}}$. Thus we can present a fast method instead of explicitly computing ${\tt reach}_{G({\mathbf{A}})}(B)$. As discussed in Section 6.2, $G({\mathbf{M}}_{k+1})$ is very similar to $G({\mathbf{M}}_{k})$, so this approach is feasible. Thus, we utilize the information of $G({\mathbf{M}}_{k})$ to construct computable projection of $G({\mathbf{M}}_{k+1})$.

Note 5

For a graph $G$, $V(G)$ denotes set of nodes in $G$.

Definition 8

$G$* is a graph with nodes $1,\cdots,2n$. $\{G_{1},\cdots,G_{s}\}$ are graphs such that $V(G_{i})\subseteq\{1,\cdots,2n\}$. We call $\{G_{1},\cdots,G_{s}\}$ an over disjoint cover of $G$ if*

1. (i)

$V(G_{i})\cap V(G_{j})=\emptyset$* for $i\neq j$;* 2. (ii)

there is $G_{i}$ such that $e\in G_{i}$ for edge $e\in G$.

For a given graph $G$, $G$’s different connected components $C_{1},\cdots,C_{s}$ is one of its over disjoint cover.

Theorem 6.3

In Problem 1, let $B=\{i\mid{\mathbf{b}}[i]\neq 0\}$, $\{G_{1},\cdots,G_{s}\}$ be an over disjoint cover of $G({\mathbf{A}})$. Let $E=\{(i,j)\mid(i,j)\in G_{i},V(G_{i})\cap B\neq\emptyset\}$. Let $I=\{i\mid(i,j)\in E\},J=\{j-n\mid(i,j)\in E\}$. If Problem 1 has a solution then ${\mathbf{A}}_{I,J}$ is a computable projection.

Proof

Since Problem 1 has a solution, by Lemma 4, there is $(i,j)\in G({\mathbf{A}})$ for $i\in B$. Thus, $B\subseteq I$.

Assume that $\{i\mid A[i,j]\neq 0,j\in J\}\not\subseteq I$, in other words, there is a $h$ such that $h\in\left(\{i\mid A[i,j]\neq 0,j\in J\}\setminus I\right)$. In other words, $h\in\{i\mid A[i,j]\neq 0,j\in J\}$ and $h\not\in I$. Let $t\in J$ such that ${\mathbf{A}}[h,t]\neq 0$ and $(h,t+n)\in G_{v}$. By definition of $J$ there is an edge $(u,t+n)\in E$ because $t\in J$.

By Definition 8, all $G({\mathbf{A}})$’s edges whose nodes contain $t$ must be completely contained in $G_{v}$. Therefore, the assumption can not hold. Hence, $(u,t+n)\in G_{v}$.

We want to prove that $(u,t+n)\not\in G_{v}$.　We prove it by contradiction. If $(u,t+n)\in G_{v}$, then by definition of $E$ all the edges of $G_{v}$ will belong to $E$, in particular, $(h,t+n)\in E$. Thus, $h\in I$. This conflicts with $h\not\in I$, so $(u,t+n)\not\in G_{v}$.

Thus, $\{i\mid A[i,j]\neq 0,j\in J\}\subseteq I$, by the same way we can prove (iii) of Definition 5. Hence ${\mathbf{A}}_{I,J}$ is a computable projection. ∎

Theorem 6.3 provides a new approach to constructing computable projection. This approach can be used to replace the computation of ${\tt reach}_{G({\mathbf{A}})}(B)$ in line 2 in locSolver.

Lemma 7

${\mathbf{M}},{\mathbf{M}}^{\prime}$* are two $n\times n$ matrices. $\{G_{1},\cdots,G_{s}\}$ is an over disjoint cover of $G({\mathbf{M}})$. Denote the set of nonzero elements in matrix $({\mathbf{M}}-{\mathbf{M}}^{\prime})$ by $\{({\mathbf{M}}-{\mathbf{M}}^{\prime})[i_{k},j_{k}]\mid k=1,\cdots,m\}$. We iteratively update $\{G_{1},\cdots,G_{s}\}$ by the following operation: merging $G_{i},G_{j}$ to $G^{\prime}$ where $G^{\prime}=G_{i}\cup G_{i}\cup\{(i_{k},j_{k}+n)\}$ if there is a link $(i_{k},j_{k}+n)$ connecting $G_{i}$ and $G_{j}$. Let $G^{\prime}_{1},\cdots,G^{\prime}_{s^{\prime}}$ be finally resulted graphs of the above iteration. Then $\{G^{\prime}_{1},\cdots,G^{\prime}_{s^{\prime}}\}$ is an over disjoint cover of $G({\mathbf{M}}^{\prime})$.*

Proof

We prove it by induction.

When $m=1$. If both ${\mathbf{M}}[i_{1},j_{1}]$ and ${\mathbf{M}}^{\prime}[i_{1},j_{1}]$ are nonzeros, then $\{G^{\prime}_{i}\}=\{G_{i}\}$ and $G({\mathbf{M}})=G({\mathbf{M}}^{\prime})$. Thus, conclusion holds.

If only ${\mathbf{M}}[i_{1},j_{1}]\neq 0$, then $\{G^{\prime}_{i}\}=\{G_{i}\}$ and $V(G({\mathbf{M}}^{\prime}))\subseteq V(G({\mathbf{M}}))$. Thus conclusion holds.

If only ${\mathbf{M}}^{\prime}[i_{1},j_{1}]\neq 0$, then $G({\mathbf{M}}^{\prime})$ only has one more link $(i_{1},j_{1}+n)$ compared with $G({\mathbf{M}})$. If there exist $G_{i},G_{j}$ such that $i_{1}\in V(G_{i}),j_{1}\in V(G_{j})$, then merge $G_{i},G_{j}$ into a graph $G^{\prime}=G_{i}\cup G_{i}\cup\{(i_{1},j_{1}+n)\}$. It is easy to check that $\left(\{G_{1},\cdots,G_{s}\}\setminus\{G_{i},G_{j}\}\right)\cup\{G^{\prime}\}$ satisfies (i)-(ii) of Definition 8. Hence the conclusion holds when $m=1$.

Assuming that the conclusion holds when $m\leq K_{1}$. When $m=K_{1}+1$, let ${\mathbf{M}}^{\prime\prime}$ be a matrix such that $({\mathbf{M}}^{\prime\prime}-{\mathbf{M}}^{\prime})$ has only one nonzero element $({\mathbf{M}}^{\prime\prime}-{\mathbf{M}}^{\prime})[i_{1},j_{1}]$, and ${\mathbf{M}}[i_{1},j_{1}]={\mathbf{M}}^{\prime\prime}[i_{1},j_{1}]$. It is easy to know that matrix $({\mathbf{M}}^{\prime\prime}-{\mathbf{M}})$ has only $K_{1}$ nonzero elements. By assumption, ${\mathbf{M}}^{\prime}$’s over disjoint cover can be constructed from ${\mathbf{M}}^{\prime\prime}$’s, which can be constructed from ${\mathbf{M}}$’s .

In summary, conclusion holds for any $m\geq 0$.∎

Through Lemma 7, we can construct over disjoint cover of $G({\mathbf{M}}_{k+1})$ from $G({\mathbf{M}}_{k})$’s. This can be used to fast compute computable projection in $(k+1)$th iteration from $k$th’s.

6.2 Incremental Change Property of ${\mathbf{M}}_{k}$’s Nonzero Pattern

Fast solving equations (25), (26) and (31) is a feasible way of improving efficiency of SMCG. In the following section we will first give an Algorithm incSolver  which utilizes the sparse and incremental change properties of matrices and vectors occurring in two consecutive equation systems to fast solve target equation. Second we will describe the three interesting phenomenons during SMCG’s iteration. And these phenomenons can let us directly employ incSolver instead of other solvers. By this way, we can fast solve equations (25), (26) and (31).

6.2.1 Fast method of solving similar linear equations

Problem 2

${\mathbf{A}}$ is a non-singular sparse matrix. ${\mathbf{A}},{\mathbf{A}}^{\prime}$ are two very similar matrices, in other words, $({\mathbf{A}}-{\mathbf{A}}^{\prime})$ has few nonzero elements. ${\mathbf{b}},{\mathbf{b}}^{\prime}$ are two very similar vectors. When there is a vector ${\mathbf{\xi}}^{\prime}$ such that

[TABLE]

we want to give an efficient algorithm to solve equation

[TABLE]

Assume that ${\mathbf{\xi}}$ is a solution of (34). Because the coefficient matrices and right-hand-sides of (33) and (34) are very similar. It is reasonable to believe that ${\mathbf{\xi}},{\mathbf{\xi}}^{\prime}$ are very similar. Thus, we only need to compute the different part for these two solutions when solving (34). The concrete algorithm of this idea is listed in Algorithm incSolver.

The outline of Algorithm incSolver is as follows: When we want to solve equation (34) when there is a vector ${\mathbf{\xi}}^{\prime}$ satisfying (33). In this case, we can believe that $({\mathbf{b}}-{\mathbf{A}}{\mathbf{\xi}}^{\prime})$ is a sparse vector. So, we firstly use Algorithm locSolver to find a vector $\Delta{\mathbf{\xi}}$ which is a solution of equation ${\mathbf{A}}\mathbf{x}=({\mathbf{b}}-{\mathbf{A}}{\mathbf{\xi}}^{\prime})$. It is easy to check that ${\mathbf{A}}({\mathbf{\xi}}^{\prime}+\Delta{\mathbf{\xi}})={\mathbf{b}}$.

6.3 Incremental Change Property of ${\mathbf{M}}_{k}$’s Nonzero Pattern

In this section we will list the three interesting phenomenons. Firstly, matrices ${\mathbf{M}}_{k},{\mathbf{M}}_{k+1}$ have little difference. Secondly, right-hand-sides in (25) and (31) also have little difference between $k$th and $(k+1)$th iteration. Thirdly, right-hand-side of (26) is very sparse. All of these phenomenons indicate that locSolver and incSolver are the proper solvers for (25), (26) and (31). So, we can use Algorithm incSolver to quickly construct solutions of (25), (26) and (31).

As described in Model 6, ${\mathbf{M}}_{k}$ is entirely defined by $K,\SS_{k}$ and their elements’ order. Therefore, changing order of $K,\SS_{k}$’s elements can give a better incremental property of matrices and vectors occurring in iteration. In the following description, without special statement the same elements of $K$ have the same matrix index between $k$th and $(k+1)$th iteration, and the same elements of $\SS_{k}$ and $\SS_{k+1}$ have the same matrix index. For obtaining incremental property, we firstly redefine transition system rules as follows:

When the entering variable is a link $e^{*}$ :

1. (a)

The same as Fig. 3 1-(a) 2. (b)

The same as Fig. 3 1-(b) 3. (c)

When the leaving variable is a link $e$:

Let $N_{k+1}=\{N_{k}\cup\{e^{*}\}\}\setminus\{e\}$. And let index of $e$ in $\SS_{k+1}$ be the same as $e^{*}$ in $\SS_{k}$, and the other links’ indices are kept the same between $\SS_{k}$ and $\SS_{k+1}$. 2. 2.

When the entering variable is a $p^{\prime}_{j}$:

1. (a)

The same as Fig. 3 2-(a) 2. (b)

The same as Fig. 3 2-(b) 3. (c)

When the leaving variable is a link $e$:

Let $Q_{k+1,j}=Q_{k,j}\cup\{p_{j}^{\prime}\},N_{k+1}=N_{k}\setminus\{e\}$. $p^{\prime}_{j}$** is a path corresponding to saturate link $e$. Append $e$ to $\SS_{k}$ to obtain $\SS_{k+1}$. And other links’ indices are kept the same between $\SS_{k}$ and $\SS_{k+1}$.**

First phenomenon.

From transition system rules in Fig. 5, we can find that size of matrix ${\mathbf{M}}_{k+1}$ only has three cases, i.e. $|{\mathbf{M}}_{k}|,|{\mathbf{M}}_{k}|+1and|{\mathbf{M}}_{k}|-1$. Below we will discuss relation between ${\mathbf{M}}_{k+1}$ and ${\mathbf{M}}_{k}$ under these three cases.

First, when $|{\mathbf{M}}_{k+1}|=|{\mathbf{M}}_{k}|$, we will give a useful fact that the number of $\left({\mathbf{M}}_{k+1}-{\mathbf{M}}_{k}\right)$’s columns which has nonzero elements is very few.

Lemma 8

In Note 4, if matrix ${\mathbf{C}}^{\prime}$ only has one column corresponding to commodity $i$ different from ${\mathbf{C}}_{k}$, then there are at most $|Q_{k,i}|$ columns of ${\mathbf{C}}^{\prime}{\mathbf{B}}_{k}-{\mathbf{D}}_{k}$ different from ${\mathbf{M}}_{k}$.

Proof

By the definition of ${\mathbf{B}}_{k}$ in (19), every column of matrix ${\mathbf{B}}_{k}$ corresponding to a commodity $i$. Let $c$ be $i$th column of ${\mathbf{C}}_{k}$. By the order of $K$, $c$ is corresponding to primary path $p_{k,i}$. Let $\beta$ be $j$th column of matrix ${\mathbf{B}}_{k}$ which is corresponding to a commodity $i$. By equation (19), the form of $\beta$ is as follows

[TABLE]

Hence, in product ${\mathbf{C}}_{k}{\mathbf{B}}_{k}$ only column $c$ of matrix ${\mathbf{C}}_{k}$ affects $\beta$. And ${\mathbf{C}}_{k}\beta=c$ is the $j$th column of product ${\mathbf{C}}_{k}{\mathbf{B}}_{k}$. Thus, in product ${\mathbf{C}}_{k}{\mathbf{B}}_{k}$, $i$th column of ${\mathbf{C}}_{k}$ only affect columns corresponding to $Q_{k,i}$. Therefore, changing column $c$’s value only changes $|Q_{k,i}|$ columns’ values of product ${\mathbf{C}}_{k}{\mathbf{B}}_{k}$, since the number of ${\mathbf{B}}_{k}$’s columns are the same as $\beta$ is $|Q_{k,i}|$. ∎

Employing Lemma 8, we can give relation between ${\mathbf{M}}_{k+1}$ and ${\mathbf{M}}_{k}$ as follows:

If we use rule 1-(c) of Fig. 5 during iteration and let $e$ have the same index of $\SS_{k+1}$ as $e^{*}$ in $\SS_{k}$, then ${\mathbf{M}}_{k+1}$ has at most

[TABLE]

columns different from ${\mathbf{M}}_{k}$. 2. 2.

If we use rule 2-(a) of Fig. 5 during iteration, then ${\mathbf{M}}_{k+1}$ at most has $|Q_{k,i}|$ columns different from ${\mathbf{M}}_{k}$. 3. 3.

If we use rule 2-(b) of Fig. 5 during iteration, then ${\mathbf{M}}_{k+1}$ at most has $1$ column different from ${\mathbf{M}}_{k}$.

In summary, when size of ${\mathbf{M}}_{k+1}$ equals size of ${\mathbf{M}}_{k}$, matrix $({\mathbf{M}}_{k+1}-{\mathbf{M}}_{k})$ has few nonzero columns.

Second, when $|{\mathbf{M}}_{k+1}|=|{\mathbf{M}}_{k}|+1$. Only after transition system rule 2-(c), this case can occur. And matrix ${\mathbf{M}}_{k+1}$ can be written as

[TABLE]

where $\theta,\rho$ are vectors and $a$ is scalar. So ${\mathbf{M}}_{k+1}$ and ${\mathbf{M}}_{k}$ are very similar.

Third, when $|{\mathbf{M}}_{k+1}|=|{\mathbf{M}}_{k}|-1$. Only after transition system rule 1-(a), 1-(b), this case can occur. And relation between ${\mathbf{M}}_{k+1}$ and ${\mathbf{M}}_{k}$ can be written as

[TABLE]

where ${\mathbf{M}}^{(i)}_{k}$ are matrices, $\theta_{i},\rho_{i}$ are vectors and $a$ is a scalar. So ${\mathbf{M}}_{k+1}$ and ${\mathbf{M}}_{k}$ are very similar.

Second phenomenon.

By the same analysis procedure as above, we can easily check that

${\mathbf{A}}_{k}$ and ${\mathbf{A}}_{k+1}$ are similar; 2. 2.

$\mathbf{c}_{k}$ and $\mathbf{c}_{k+1}$ are similar; 3. 3.

${\mathbf{\beta}}_{k}$ and ${\mathbf{\beta}}_{k+1}$ are also similar.

Thus, right-hand-sides in (25) and (31) in $k$th and $(k+1)$th iterations are similar too.

Third phenomenon.

In equation (26), if associate entering variable of ${\mathbf{\beta}}$ is an edge $e^{*}$, the form of $\beta$ is as follows

[TABLE]

By the discussion in the beginning of Section 5, every column of ${\mathbf{C}}_{k}$ is a sparse column in general. So, ${\mathbf{C}}_{k}{\mathbf{\beta}}_{K}-{\mathbf{\beta}}_{\SS_{k}}$ is also a sparse vector in general at Model 6.

Otherwise, ${\mathbf{\beta}}$ is a basic vector associated with a path $p$. Then, in Model 6 ${\mathbf{\beta}}_{K}$ only has one nonzero element which is $1$. And ${\mathbf{\beta}}_{\SS_{k}}$ is also sparse in general by discussion in the beginning of Section 5. Thus, ${\mathbf{C}}_{k}{\mathbf{\beta}}_{K}-{\mathbf{\beta}}_{\SS_{k}}$ is the subtraction of two sparse vectors, so is also a sparse vector in general.

6.4 Fast Solving Equations During Iteration

By discussion in Section 6.3, coefficient matrix and right-hand-side of (26) are both sparse. Hence, we can employ algorithm locSolver to solve (26).

In addition, coefficient matrices and right-hand-sides of (25) and (31) are very similar between $k$th and $(k+1)$th iteration. So, in the following we will employ algorithm incSolver to solve equations in $(k+1)$th iteration from solution of $k$th ones.

1. (1)

When $|{\mathbf{M}}_{k+1}|=|{\mathbf{M}}_{k}|$, we directly employ incSolver to solve (25) and (31). 2. (2)

When $|{\mathbf{M}}_{k+1}|=|{\mathbf{M}}_{k}|+1$, we firstly extend solution of corresponding equation of $k$th by setting the element of new index as [math]. After this extension, we employ incSolver to solve (25) and (31). 3. (3)

When $|{\mathbf{M}}_{k+1}|=|{\mathbf{M}}_{k}|-1$, we firstly narrow solution of corresponding equation of $k$th by deleting element of lacking index. After this narrowing, we employ incSolver to solve (25) and (31).

7 Experiments

7.1 Environment

The algorithm of SMCG is implemented as a C++ program. Compilation was done using g++ version 5.4.0 with optimization flags -O2. We use latest LAPACK (version 3) and latest KLU  which is contained in tool SuiteSparse 4.5.6111 http://faculty.cse.tamu.edu/davis/suitesparse.html. All tests are done on a 64-bit Intel(R) Core(TM) i5 CPU 7400 @ 3.00GHz with 8GB RAM memory and Ubuntu 16.04 GNU/Linux.

We use incCG, kluCG and lapackCG to denote implementations of SMCG with incSolver, KLU and LAPACK  as linear equation solver, respectively. In other words, except for linear equation solver, the other parts of incCG, kluCG and lapackCG  are the same.

Random test cases are created by generator $R(n)$ where $n$ is the number of nodes. The average node degree (sum of in degree and out degree) is $10$. Each edge is generated by two random integers between $1$ and $n$ as its source and target node indices. The edge capacity is a random integer between $1$ and $300$ and edge weight is a random integer between $1$ and $10$. The source and target indices of commodity are two random integers between $1$ and $n$. Commodity demand is a random integer between $1$ and $100$. Every case has $1000$ commodities.

In Table 1, you can see that shortest path computing and linear equation solving are two major time consuming parts of implementations. And except for $R(1000)$, $R(1500)$ and $R(2000)$, the total time of incCG, kluCG and lapackCG  are almost equal to sum of shortest path computing time and linear equation solving time. And the shortest path computing time of different implementations are almost the same. In addition, considering total time, when linear equation solving is dominating part, incSolver will achieve high speedup. For example, we can see in case $R(1000)$ using incSolver instead of LAPACK will achieve $19\times$ improvement. On the other hand, when linear equation solving costs less time, incSolver can reduce linear equation solving to a negligible fraction. For example, in case $R(7000),\cdots,R(39000)$, using incSolver instead of LAPACK will reduce linear equation solving time to less than $1\%$ of the total.

In Fig. 6, we can see that incSolver outperforms KLU and LAPACK on all the test cases. Among these cases, KLU’s speedup is between $19$ and $199$ compared with LAPACK, while incSolver achieves a speedup from $37$ to $341$. When comparing incSolver with KLU, incSolver’s speedup is between $1.7$ and $2.1$. As incSolver is a prototype implementation, we believe that incSolver has great potential for improvement.

8 Conclusion

In this paper, for speeding up linear equation solving part in column generation for multi-commodity flow problem, firstly, we use transition system view to describe the procedure of column generation. This view can help us better understand the procedure of column generation and it also helps us conveniently present following improvement. Secondly, we discuss the sparse property of coefficient matrix. In the SMCG the average number of nonzero coefficient in each row of coefficient matrix is very few. In our test it is less than $5$, even when the dimension of matrix is more than $1000$. Finally, we present two algorithms. The first is a fast algorithm locSolver  (for localized system solver) which can reduce the number of variables in solving a linear equation system when both the coefficient matrix and right-hand-side are sparse. The other is an algorithm incSolver  (for incremental system solver) which utilizes similarity during the iteration in solving a linear equation system. All algorithms can be used in column generation of multi-commodity problem. Preliminary numerical experiments show that the algorithms are significantly faster than existing algorithms. For example, under random test cases incSolver delivers up to $341\times$ (from $37\times$) improvement in the linear equation solving part compared with LAPACK. In addition, considering total time, when linear equation solving is dominating part, incSolver will achieve high speedup. For example in some tests using incSolver instead of LAPACK will achieve $19\times$ improvement. On the other hand, when linear equation solving costs less time, incSolver can reduce linear equation solving to a negligible fraction. For example in some cases using incSolver instead of LAPACK will reduce linear equation solving time to less than $1\%$ of the total.

9 Acknowledgements

The authors would like to thank Prof. Zongyan Qiu who helps to improve the presentation of the article.