6.1.1 A basic instance of the NUM problem
In this section block-coordinate implementations of the mapping Gare tested on a benchmark instance of the NUM problem discussed in Section 3.4.5. We consider the formulation of the NUM problem studied in [ZRJO11] and already encountered in Example 2.1, where a network with node setN = {1, ..., n}and edge setEis represented by a directed graph G = (N;E). Each edge connects a pair of nodes i, j ∈ N and is denoted by the ordered pair (i, j)∈ E, where node i is arbitrarily chosen as the origin and j as the destination. A pair of nodes is connected by at most one edge. Direct transmissions of information are only possible between neighbours, i.e. pairs of nodes connected by an edge. The set of neighbours of a node i is denoted by Ni ⊂ N and we have j ∈ Ni i (i, j)∈ E or (j, i)∈ E.
The considered problem can be stated as minimise
x
∑n
i=1fi(xi)
subject to Ax= b . (6.1)
The objective is to minimise a convex functionf : Rp → Rof the information ows on the edges, gathered in the composite vector x = (x1, ..., xn), where, for i ∈ N, the vector xi ∈ Rpi is assigned to node i and contains the pi ow variables (arranged in arbitrary order) corresponding to all the edges with origin i (∑n
i=1pi = p with pi ≥ 0 for i = 1, ..., n). The objective function is assumed to be additively separable with respect to x1, ..., xn, and we write
6.1. Networks with constant properties 149
f = ∑n
i=1fi with fi : Rpi → R convex. The network supports a single information ow specied, for each node i ∈ N, by the incoming rate bi ≥ 0 if i is a source node, or by the outgoing rate bi < 0 if i is a sink node. Flow conservation at each node is guaranteed by the constraint Ax = b, where b = (b1, ..., bn) is a vector of Rp such that ∑n
i=1bi = 0, and A denotes the node-edge incidence block matrix. For i, j ∈ N, Aji is a 1×pi row vector such that
Aji =
⏐
⏐
⏐
⏐
⏐
⏐
(1, ...,1)′ if j =i,
(0, ...,0)′ if j ̸=i, (i, j)∈/ E, (0, ...,0
k−1
,−1,0, ...,0)′ if (i, j) is the kth edge with origin i.
The ow conservation constraint can be rewritten as hi(x) = 0 for i ∈ N, where we dene hi(x) =∑
j∈Ni∪{i}Aijxj −bi.
The problem (6.1) is an instance of Problem 3.3 in Section 3.4.5 with ui = 1, vi = 0, mi = 1 and thus Yi = R for i∈ N. From (B.41), we nd
∇ij2g(y) =−∑
k∈N˜ij
Aik[∇2fk(x∗k(y))]−1A′jk, i, j ∈ N, y ∈ Y. (6.2) Since ∇2f = diag(∇2f1, ...,∇2fn) is a block diagonal matrix, (6.2) can be rewritten as
∇2g(y) =−A[∇2f(x∗(y))]−1A′, y ∈ Y, (6.3) where [∇2f(x∗(y))]−1 is diagonal with positive elements. Hence −∇2g can be seen as a scaled version of the form AA′, which is is called the Laplacian matrix of the directed graph and such that
[AA′]ij =
⏐
⏐
⏐
⏐
⏐
⏐
deg(i) if i= j
−1 if i∈ Nj or j ∈ Ni 0 otherwise
, (6.4)
where deg(i) denotes the degree if node i, i.e. the number of nodes connected to i by an edge or, equivalently, the number of neighbours of i.
In our simulations, we use the objective function suggested in [ZRJO11]
and dened, for all i∈ N, by
fi : x ∈ Rpi → fi(x) =∑pi
j=1(eγxj +e−γxj), (6.5) where xj denotes the jth component of x and γ is a positive constant.
Ft(ρ)
ρ S
N(0) N(1)
N(2) N(3)
N(4) N(5) N(6) N(7)
100 10−1
10−2 10−3
10−4 0
1 2
1
Generation of t = 1000 networks with n = 50 and random links, and computation, for S (dotted line) andN(q) withq= 0, ...,7(plain lines), of Ft(ρ) = 1t∑t
s=1δρs≤ρ, where ρs denotes the spectral radius of the matrix convergence rate of the considered method observed for the sth sample.
Figure 6.1: Spectral radii of the Gauss-Seidl and accelerated methods Rates of convergence
Fig. 6.1 compares, for this particular formulation of the NUM problem, the spectral radii ρ of the asymptotic convergence rates of the Gauss-Seidel algo-rithm S and the family of accelerated synchronous algorithms {N(q)} given in (4.82).
We generate connected networks with 50 nodes where the existence of each communication edge is decided at random with probability 12. The empirical probability distribution of the spectral radii of the convergence rates derived from the simulations using (4.83) and (4.85) are depicted on the gure. One sees that ρ is smaller for S than for N(0) with similar amounts of informa-tion transmitted between the nodes. In other words, the Gauss-Seidel algo-rithm (S) used with local Newton scaling has not only the advantage over the Jacobi mode (G ≡N(0)) that it does not require global synchronism, but also converges asymptotically faster in terms of local gradient descents per node.
This result is in accordance with the Stein-Rosenberg theorem (Theorem A.1 in Appendix A). Indeed, it is easily seen that −∇2g meets the conditions of Theorem A.1 by using the decomposition proposed in (4.57), and writing
−∇2g = D −L−L′, where, for any y ∈ Y, we know from (6.2) that D(y) is diagonal and positive denite and L is strictly lower triangular with non-negative elements, and we note that the forms D−1(L+L′)and (D−L)−1L′ give the asymptotic matrix convergence rates of the Jacobi and Gauss-Seidel algorithms, respectively.
By increasing the value of the parameter q and the quantity of exchanged information, ρ can be further reduced for N(q) and becomes arbitrarily small for q large enough.
Constrained feasible sets and sequential implementations.
The sequential implementations of the mapping G are now tested in con-strained sets. Our problem is modied by considering, in each node i ∈ N,
6.1. Networks with constant properties 151
an additional inequality constraint of the type di(x) ≤ 0, where di(x) =
∑
j∈Ni∪{i}ζij(xj). We introduce an operator (·)+ which converts the scalar components of a vector v = (v1, ..., vl) into positive values, i.e, (v)+ .
= (|v1|, ...,|vl|). A limitation of the total activity of each node i ∈ N can then be obtained by setting
ζij(xj) =
⏐
⏐
⏐
⏐
Aii(xi)+−di if j = i,
−Aij(xj)+ if j ∈ Ni, (6.6) where di is a positive constant. The new problem can be stated as
minimise
x
∑n
i=1fi(xi) subject to ∑
j∈Ni∪{i}Aijxj −bi = 0, i ∈ N
∑
j∈Ni∪{i}ζij(xj)≤0, i ∈ N
(6.7) and is an instance of Problem 3.3 withui = vi = 1, mi = 2, andYi = R×R≥0
for i∈ N.
The parameters {di}i∈N are chosen large enough so that the problem is well-conditioned and has one nite solution y⋆ ∈ Y. We randomly generate the network with 50 nodes and 100 links depicted in Fig. 6.2. Then we run the Gauss-Seidel implementation of G from various starting points chosen at random in Y, as well as a sequential arbitrary implementation mode where it is assumed that the time t between two local projected gradient descents in a subset Yi is distributed identically and independently for all i ∈ N with exponential distribution F(t) = 1−exp(−t), and that the block-coordinate sequences (ρk) are generated accordingly. Both algorithms are successively tested with local Newton and diagonal scaling.
The assumption on {di}i∈N ensures that the nonnegativity constraints specifying Y are all active at y⋆. By Proposition 4.7, local Newton scaling and diagonal scaling are thus equivalent in the reduced space at y⋆, which is reached in nite time by all the algorithms in accordance with Proposition 4.6.
Table 6.1 displays the mean number k⋆ of gradient descents per node needed by the algorithms to uncover {Ai(y⋆)}i∈N together with an estimate of the standard deviation for random starting points. Thanks to second-order scal-ing, the Gauss-Seidel implementations require only 1 to 3 cycles to identify these constraints, while 4 to 5 updates per node are needed on average by the random algorithms. Hence the eort-saving diagonal scaling strategy seems more appropriate in this particular example, although local Newton scaling sometimes proves signicantly faster in more dicult problems.
The table also displays, for various precisions ϵ, the number of gradient descents per node observed before the norm of the projected gradient be-comes less than ϵ. This quantity is given for the Gauss-Seidel algorithms by
(a) Topology (b) Optimal routing
Random generation of a network with 50 nodes and 100 edges depicted in Figure 6.2(a), where the source nodes are depicted in white and sink nodes in grey, with diameters directly proportinal to the rate of incoming or outgoing information. Figure 6.2(b) displays the optimal information ows, where the thickness of an edge is directly proportional to the optimal rate of information transmitted over the edge.
Figure 6.2: Network with 50 nodes and 100 edges
Table 6.1: Convergence of the Gauss-Seidel and random algorithms
Statistics of k⋆ and τ(·) for the Gauss-Seidel mode with local Newton scaling (Sˇ) and diagonal scaling (S¨), and for the random mode with local Newton scaling (Aˇk) and diagonal scaling (A¨k) [estimated mean ±standard deviation].
k⋆ τ(10−1) τ(10−2) τ(10−3) τ(10−4) τ(10−5) τ(10−6) Sˇ 2.2±0.4 26.3±3.7 44.4±4.1 62.7±4.6 81.0±5.1 99.5±5.4 117.9±5.9 S¨ 3.0±0.5 24.2±3.4 42.4±3.6 61.0±3.6 79.4±3.8 98.1±3.9 117.2±4.3 Aˇk 4.3±1.2 48.7±5.2 86.6±5.4 122.9±5.9 161.0±5.6 199.0±6.9 238.2±8.1 A¨k 5.2±1.5 43.8±7.8 80.3±7.4 116.4±7.3 154.7±8.9 193.0±9.8 231.9±9.3
τ(η) = min¯k{¯k : ∥E(yk)′∇g(yk)∥ < η, ∀k ≥k¯}, whereE is dened in (4.56).
The magnitude of the projected gradient vanishes in a seemingly linear fashion for the algorithms. In this problem, we note that convergence is about twice faster for the Gauss-Seidel mode than the random mode.
6.1.2 Network lifetime-maximisation
In this section we address a more complex routing problem, similar to those treated in [ML06, BAR11, Bil12], where the aim is to optimise the lifetime, as previously dened in (2.14), of a network subject to constraints concerned with the capacity of the transmission channels and the maximum power con-sumption of the nodes.
We consider a connected network with n nodes distributed in a region.
The set of nodes is denoted by {1, ..., n} and the set of directed edges by E. Each node i has a limited battery supply bi, and either generates packets of information or collects the packets forwarded by other nodes. The
parame-6.1. Networks with constant properties 153
ter si is the rate of information generated (si ≥ 0) or collected (si < 0) by node i. The set of neighbours Ni of a node i is dened like in Section 6.1.1 as the set of the nodes located at short distances from i. Direct transmissions of information are only possible between neighbours.
In the present setting, a nonnegative variable x˘, already introduced in Ex-ample 2.4, stands for an upper bound on the `inverse lifetime' of the network.
Our objective is thus to minimise x˘. Since x˘ is a global variable, we intro-duce the n auxiliary variables x˘1, ...,x˘n, where x˘i is a copy of x˘ assigned to node i together with the consistency constraints of the type x˘i = ˘xj for any neighbouring node j ∈ Ni (cf. Example 3.9). We use nonnegative ow vari-ables as in Example 2.3, which facilitate the introduction of the constraints.
To each node i is assigned a vector x¯¯i containing the ow varibles ⃗xij of all the edges (i, j) ∈ E starting form i, as well as the variables ⃗xji for all the edges (j, i) leading to i. The problem variables form a composite vector x = ((¯x¯1,x˘1), ...,(¯x¯n,x˘n)), where (¯x¯i,x˘i) contains all the nonnegative vari-ables assigned to the sensor i ∈ N.
By denition of the network lifetime, the constraintpi(¯x¯i) ≤bix˘i should be introduced for every nodei∈ N, wherepi(¯x¯i)is an estimation of the expected power consumption at i under the local ow policy x¯¯i. Reception power con-sumption is not considered, and it is assumed that the power concon-sumptionpi is given by the linear model
pi(¯x¯i) =∑
j:(i,j)∈E¯eiij⃗xij +∑
j:(j,i)∈E¯eiij ji⃗x , ∀i ∈ N, (6.8) where, for every i ∈ Ni, e¯iij is the mean energy cost of the transmission of a packet from i to j.
Since our intention is to apply dual distributed methods, we need to make sure that the dual function g enjoys attractive properties. In the presence of continuous constraints, we know from Corollary 3.1 that the compact-ness of the feasible set and the strict convexity of the cost function are suf-cient conditions for the dual function to be dierentiable on its domain.
Compactness of the feasible set X is ensured by adding the additional con-straint 0 ≤ (¯x¯i,x˘i) ≤ x¯i at each node i ∈ N. Treating all the inequality constraints as side constraints, the feasible set reduces to the Cartesian prod-uct set X = ∏
i∈NXi, where we dene
Xi = {(¯y,¯ y)˘ ∈ R2|Ni|×R: 0 ≤(¯y,¯ y)˘ ≤x¯i, pi(¯x¯i)−bix˘i ≤ 0}, ∀i∈ N, (6.9) and |Ni| denotes the number of neighbours of i. The positive vector x¯i has 2|Ni| + 1 components (i ∈ N). It species both limitations on the ca-pacities of the communication channels of the sensor i and an arbitrary large upper bound on the inverse lifetime variable x˘i.
(a) Topology (b) Optimal routing
Figure 6.3: Random network with25 nodes
Strict convexity of the cost function in all the problem variables is ob-tained by regularisation. We proceed as in (3.59) by adding a separable term quadratic in all the variables with a small positive coecientϵ. All in all, the lifetime-maximising routing problem is formulated as the quadratic program
minimise
x∈X
∑
i∈Nx˘2i +ϵ∥x∥2
subject to Aix¯¯i = si ∀i ∈ N
˘
xi = ˘xj ∀i ∈ N, j ∈ Ni
(6.10) where Ai is the local incidence matrix at i, and X is dened in (6.9). Us-ing the results of Section 3.4.3, it is straightforward to show that the dual function g obtained after dualisation of the equality constraints in (6.10) is everywhere dierentiable and piecewise quadratic over the dual feasible setY (cf. Example 3.6), which for this problem is the whole dual space. The dual problem can be solved in a distributed manner by the sensors using the gra-dient methods explored in Chapter 4. In the rest of this section we report basic numerical experiments showing both the necessity of using line-search routines for the step-sizes and the importance of scaling.
A collection of sensor networks are randomly generated1, so that the net-works are connected and distributed uniformly in a circular region. Fig-ure 6.3(a) depicts the topology and connectivity of an instance of these net-works with 25 nodes and 104 edges. It is assumed that bi is equal for all the nodes and that each si is chosen randomly and uniformly in the inter-val [−1,1], except in one node j where sj = −∑
k̸=jsk. We make the as-sumption of low trac conditions by assigning large values to the channel capacities {x¯i}i∈N, and use, as a path loss model for the transmissions, the d−4 power law suggested Section 2.1.3. Figure 6.3(b) depicts the ow policy solution of (6.10), where the segment sections are proportional to the corre-sponding packet ows.
1Since the eects of scaling and line-search are clearly noticeable on small problems with only hundreds of variables, and increasing the complexity of the problem tends to prohibitively slow down the convergence steepest gradient method, the upcoming tests are restricted to small networks.
6.1. Networks with constant properties 155
Table 6.2: Convergence of the gradient projection and the Gauss-Seidel algorithms
Statistics of τ(·) andφ(·)for the gradient projection algorithm with decreasing step-size (GP), the Gauss-Seidel mode with unit scaling (S˙), and the Gauss-Seidel mode with local Newton scaling (Sˇ) [estimated mean ±standard deviation].
τ(10−1) τ(10−2) τ(10−3) τ(10−4)
GP (1.6±1.4)×103 (10.4±9.3)×103 (6.3±7.4)×104 (1.4±1.4)×105 S˙ 7.3±1.1 15.8±3.0 32.1±9.2 48.9±16.7 Sˇ 6.7±1.3 9.6±1.5 11.7±2.0 13.4±2.5
φ(10−1) φ(10−3) φ(10−5) φ(10−7)
GP (7.1±4.6)×104 (4.6±3.2)×105 (2.9±3.2)×106 (7.1±9.5)×106 S˙ (7.8±3.2)×102 (18.0±7.5)×102 (3.4±1.2)×103 (5.4±1.9)×103 Sˇ (6.6±3.0)×102 (9.0±4.6)×102 (10.2±4.8)×102 (11.2±4.9)×102
The dual of (6.10) is solved for all the network samples by successively applying the three following distributed optimisation methods: the gradient projection algorithm (3.80) with identity matrix scaling (steepest descent) and with decreasing step-size lawak = 0.1k−0.5 (GP), the Gauss-Seidel implemen-tation of G with identity scaling (S˙), and the Gauss-Seidel implementation of G with local Newton scaling (Sˇ).
Since the primal of (6.10) only explicitely displays equality constraints, the gradient of the dual function is expected to vanish to 0 at dual solutions and the function τ dened in the previous section may be used to analyse the speed of convergence of the algorithms. A second function φ is introduced to characterise the quantity of information exchanged by the sensors during the optimisation process. Given an optimisation algorithm and a positive scalarη > 0, φ(η)is dened as an estimate of the average number of messages each sensor exchanges by message passing up to timeτ(η) during a run of the algorithm. Table 6.2 displays, for various precisions η, the expectation and standard deviation of τ and φ for the three methods as estimated from the random experiments. We see that the number of iterations and information transfers needed by the steepest gradient method is highly prohibitive even for small problems. Convergence is, however, much faster when line-search routines are implemented. According to the table, the convergence times and numbers of data transfers for this problem are further reduced when line-search is combined with second-order scaling.
In view of these experiments, it appears that scaling and line-search ac-celerate the convergence of gradient methods by several orders of magnitude.
The implementation of such techniques is thus needed in practice. In the net-work depicted in Figure 6.3, for instance, Sˇ identied the active constraints
at the optimum after only 7 iterations and solved the problem in a few tens of node updates, while hundreds of thousands of iterations were needed by the steepest gradient algorithm, which is clearly unmanageable in practical problems.