Efficient Capacity Computation and Power Optimization for Relay Networks

The capacity or approximations to capacity of various single-source single-destination relay network models has been characterized in terms of the cut-set upper bound. In principle, a direct computation of this bound requires evaluating the cut capacity over exponentially many cuts. We show that the minimum cut capacity of a relay network under some special assumptions can be cast as a minimization of a submodular function, and as a result, can be computed efficiently. We use this result to show that the capacity, or an approximation to the capacity within a constant gap for the Gaussian, wireless erasure, and Avestimehr-Diggavi-Tse deterministic relay network models can be computed in polynomial time. We present some empirical results showing that computing constant-gap approximations to the capacity of Gaussian relay networks with around 300 nodes can be done in order of minutes. For Gaussian networks, cut-set capacities are also functions of the powers assigned to the nodes. We consider a family of power optimization problems and show that they can be solved in polynomial time. In particular, we show that the minimization of the sum of powers assigned to the nodes subject to a minimum rate constraint (measured in terms of cut-set bounds) can be computed in polynomial time. We propose an heuristic algorithm to solve this problem and measure its performance through simulations on random Gaussian networks. We observe that in the optimal allocations most of the power is assigned to a small subset of relays, which suggests that network simplification may be possible without excessive performance degradation.

is the Gaussian relay network, which models power limited transmitters and received signals corrupted by additive white Gaussian noise.
While the capacity of some network models (e.g. wireless erasure and ADT) is well characterized, the capacity of the Gaussian relay network, even in its simplest form with one transmitter, one relay, and one receiver, is in general unknown. The best known capacity upper bound is the so-called cut-set bound. A cut Ω of a network can be considered as a subset of nodes which includes the source node and excludes the destination node. For this cut, the capacity F (Ω) is defined as the maximum rate that information can be transferred form the nodes in Ω to the nodes that are not in Ω conditioned on the fact the information on Ω c (the nodes that are not in Ω) is known. The cut-set upper bound is the minimum cut capacity over all the possible cuts.
In the Gaussian setting, there are several capacity lower bounds based on different communication schemes, such as amplify-and-forward, decode-and-forward, compress-andforward, quantize-and-forward, etc. [5], [7], [18]. Recently, Avestimehr, et al. [2] made significant progress in the capacity characterization of Gaussian relay networks by showing that a quantization and coding communication scheme can achieve a communication rate within a constant gap of the cut-set upper bound, where the gap only depends on the number of nodes in the network (i.e. it is independent of the channel gains and power levels). However, the evaluation of the achievable communication rate, which is necessary to implement the scheme, requires the computation of the cut-set bound for the network. Assuming that for a given cut the cut capacity is easy to compute, finding the cut-set upper bound can be a challenging problem. For a network with n relays there are 2 n different cuts and a greedy algorithm needs exponential time in the number of relays to compute the cut-set bound.
In this work we show that the achievable rate of the scheme of [2] for the Gaussian relay network can be computed in polynomial time, and as a result, can be computed efficiently. This result is obtained by showing that the cut capacity of a fairly large class of networks under the assumption of independent encoding at the nodes in Ω is a submodular function. For the special case of layered relay networks, [27] showed the equivalent of our submodularity result simultaneously with our conference version of this paper [25]. Submodularity properties of conditional entropy (in terms of which cutcapacities are expressed) have also been used in [2], [3] to bound the cut-capacity of a network in terms of the cutcapacity of the corresponding unfolded graph 1 .
Existing results on minimization of submodular functions provide algorithms with polynomial running time O(n 5 α + n 6 ), where α is the time that it takes to compute F (Ω) and n is the number of nodes in the network [20]. In addition, there exist possibly faster algorithms without polynomial time performance guarantees based on Wolfe's minimization norm algorithm [12]. In Section VI, by simulations, we show that the cut-set bound for a Gaussian relay network with around 300 nodes can be computed on a laptop computer in about a minute using a Matlab package for submodular minimization provided in [19].
Our results, extend and generalize previous results for the ADT model. This model can be seen as a high signalto-noise-ratio (SNR) approximation of the Gaussian model, incorporating the effects of broadcasting and superposition of signals while de-emphasizing the effects of noise. Amaudruz et al. [1] showed that the cut-set bound for a layered 2 ADT model can be computed efficiently. They have extended graph flow algorithms such as Ford-Fulkerson's in a nontrivial way to find the maximum possible linearly independent (LI) paths in the network. They showed that the capacity of the network is equal to the maximum number of (LI) paths and can be computed in time O(M · |E| · C 5 ), where M is the maximum number of nodes per layer, |E| is the total number of edges and C is the capacity of the network. Moreover, they showed that the capacity can be achieved by using a very simple one-bit processing at the relay nodes. Later Goemans et al. [13] showed that the deterministic model is a special case of a flow model based on linking systems, a combinatorial structure with a tight connection to matroids. As a by-product, they obtained the submodularity of the cut capacity for layered ADT networks. Using this observation they provided various algorithms related to matroid theory to compute the cut capacity of the layered deterministic model based on finding intersection or partition of matroids. These results led to faster algorithms to compute the capacity of large layered ADT networks. In addition, there has been other extensions on improving the running time of the current algorithms for computing the capacity of ADT networks [9], [10], [22], [23].
In addition to showing that the capacity within a constant gap of the Gaussian relay network can be computed in polynomial time, our results allow us to compute in polynomial time the capacity of the wireless erasure network. Furthermore, we provide a simple proof for the computability in polynomial time of the capacity of the layered and non-layered ADT networks.
Building on the submodularity of the cut-capacity for independent encoding at the nodes, we show that, in the Gaussian setting, it is possible to efficiently optimize the power allocated to the source and relay nodes. We consider two power optimization problems: (i) minimize the total power satisfying a minimum source-destination data rate constraint and power constraints at each node; (ii) maximize the sourcedestination data rate satisfying total and individual power constraints at the nodes. Since the capacity of the Gaussian relay network is approximately given by the cut-set upper bound with independent encoding at the nodes, we use this cut-set bound to characterize data rate in the optimization problems. We show that these optimization problems can be solved in polynomial time and use simulations to get insights about some properties of the optimal power allocations for networks of various sizes. We observe that optimal power allocations assign most of the power to a small subset of nodes and that setting the power to zero in the remaining nodes (i.e. removing these nodes from the network) often results in a small rate loss. Nazaroglu, et al. showed in [26] that for the special case of the N -relay Gaussian diamond network a fraction k/(k + 1) of the total capacity can be approximately achieved by using only k of the total N relays. This suggests that the diamond network can be significantly simplified by tolerating a small performance loss. Our results provide a numerical counterpart to the fundamental performance bounds derived in [26] and suggest that network simplification may also be possible in more general Gaussian relay networks.
We obtain these results by considering a general framework to compute the cut-set bound. We assign transmit signal random variable X i to node i ∈ V and we assume the probability distribution over the signals X 1 , X 2 , . . . , X n to be independent, i.e p(X 1 , X 2 , . . . , X n ) = p 1 (X 1 )p 2 (X 2 ) · · · p n (X n ). We also assign received signal random variables Y i 's to each node. The network is defined by the transition probability function f (Y 1 , Y 2 , . . . , Y n |X 1 , X 2 , . . . , X n ). We further assume that the transition probability function is of the form f 1 (Y 1 |X 1 , . . . , X n ) · · · f n (Y n |X 1 , . . . , X n ), meaning that the received signals are independent conditioned on the transmitted signals in the network. For such networks we show that F (Ω) = I(Y Ω c ; X Ω |X Ω c ) 3 is submodular with respect to Ω. Later we show that for ADT networks, the Gaussian relay network and the wireless erasure network, we can find p 1 (X 1 ) · · · p n (X n ) such that min Ω F (Ω) becomes equal to the capacity or the capacity within a constant gap. In other words, the min-cut problem for these networks can be cast as a minimization of a submodular function.
The paper is organized as follows. In Section III we show that for specific type of networks the cut value, F (Ω), is a submodular function. We then show in Section IV that for many wireless network models such as the ADT deterministic network, Gaussian relay network and wireless erasure network the capacity or an approximation to the capacity can be cast as a minimization of F (Ω). In Section V we study two power optimization problems and show that they can be solved efficiently. Finally, in Section VI we describe results related to solving optimization problems involving submodular functions and perform power optimization in various randomly generated networks of different sizes. We start by introducing the notation used in the rest of the paper.

II. NOTATION
Let V denote the set of nodes in the network and |V| its cardinality. For any subset A of nodes we denote by V\A or A c the set of nodes in V that are not in A. We assume V\A ∪ B = V\(A ∪ B). A cut Ω is defined as a subset of nodes in V. A cut splits the nodes in the network into two groups, the nodes that are in Ω and the ones that belong to V\Ω. Random variables are shown in capital letters such as X i and Y i . We use boldface letter for vectors, e.g. x is a constant vector and X is a random vector. We use X Ω to denote (X v1 , X v2 , . . . , X v |Ω| ) with v i ∈ Ω. The function I(X; Y |Z) is the mutual information between random variables X and Y conditioned on random variable Z. With a slight abuse of notation we use H(X) to denote either the entropy or differential entropy of the discrete or continuous random variable X [6]. By F p we denote a finite field with p elements. Finally, all the log(·) functions are in base two.

III. SUBMODULARITY OF CUT-SET FUNCTION
Submodularity arises in many combinatorial optimization problems and large body of research has been developed on minimizing or maximizing submodular functions under various constraints.
A submodular function f : 2 V → R is defined as a function over subsets of V with diminishing marginal returns, i.e. if A, B ⊆ V with A ⊆ B and any v ∈ V\B, The theorem below establishes the submodularity of the cut capacity function of a general relay network under some special assumptions. This theorem will be used in Section IV to prove that the capacity or an approximation to the capacity of various specific relay network models can be computed by minimizing a submodular function.

Theorem 1. Consider a network consisting of nodes in V.
Each node sends a message X i , i ∈ V and receives Y i , i ∈ V. If the messages are independent p(X 1 , X 2 , . . . , X |V| ) = p 1 (X 1 )p 2 (X 2 ) · · · p |V| (X |V| ) and conditioned on the sent messages the received messages are independent, then the function where (a) is the definition of mutual information and (b) is from the chain rule for the entropy function. Therefore, In the following example we show that if the signals at the nodes are correlated then F (A) is not necessarily a submodular function. Example. Consider a symmetric Gaussian diamond network with two relays such that the channel gains from source to relays are equal to one and from relays to destination are equal to three. Letting X s , X r1 , and X r2 be the signals transmitted at the source and relay nodes, then the received signals at relays and destination are given by 1). For this example, we set the probability distribution of X s , X r1 , X r2 to be jointly Gaussian with zero mean and covariance matrix Finally, consider the sets A = {s, r 1 } and B = {s, r 2 }. Figure 1 shows how the function depending on the value of ρ. It follows that in general F (·) is not a submodular or a supermodular function when the there is correlation among the signals at the nodes.

IV. WIRELESS NETWORK MODELS
In this section, by applying the result of Theorem 1, we show that the capacity or an approximation to the capacity for the ADT deterministic network, Gaussian relay network, and wireless erasure network can be cast as a minimization of a submodular function.

A. Deterministic model (ADT)
We start by briefly describing the network model of [4] 4 . In this model, each link from node i to node j has an associated non-negative integer gain n ij . Each node i ∈ V transmits a signal X i and receives a signal Y i , both in F q p where q = max i,j n ij . At any given time, the received signal at node j is given by where d is the destination node, the shifting matrix S is given by and the sums and products are in F p . For a given cut Ω of the network, where Ω includes the source node and excludes the destination node, we can stack together the input vectors X i , i ∈ Ω and output vectors Y i , i ∈ Ω c , and define a transition matrix Λ Ω that gives the input-output relationship of these vectors according to (1). It is shown in [2] that the capacity of the deterministic network is equal to min Ω rank(Λ Ω ). We show next in Theorem 2 that rank(Λ Ω ) is submodular, and hence the capacity can be computed by minimizing a submodular function.
, and let G be the set of cosets of N in F n p . Pickx i to be an element in the ith coset of N for i = 1, 2, . . . , |G |, and set y i = Ax i . Notice that y i = y j if i = j. Now, if we choose x uniformly at random from elements of F n p with probability 1/|F n p |, then the mapping Ax maps x to {y 1 , y 2 , . . . , y |G | } uniformly at random with probability 1/|G |. In addition, the cosets of N form a partition of F n p into p n /|N | sets. Also rank(A) + log p (|N |) = n. Thus, log p |G | = rank(A).

Theorem 2. For a deterministic model, given a cut Ω assume
Λ Ω is the transition matrix form nodes in Ω to nodes in Ω c . Set D(Ω) = rank(Λ Ω ), then D(Ω) is submodular.

Remark 1.
A special case of Theorem 2 for layered ADT networks was proved in earlier works [13], [23].
From the definition of transition matrix, Λ Ω , if we assume for the cut Ω, s = (s 1 , s 2 , . . . , s k ) t symbols are being sent from nodes in Ω and r = (r 1 , r 2 , . . . , r ℓ ) t symbols are being received by nodes in Ω c then r = Λ Ω s. Then we can write is the result of Proposition 1 and the fact the s has uniform probability distribution.
Notice that for the independent probability distribution on the sources the received signals are independent conditioned on transmitted signals so, based on Theorem 1,

B. Gaussian relay network
The Gaussian network model captures the effects of broadcasting, superposition and noise of power constrained wireless networks. In this model, at any time index (which we omit) the received signal at node j ∈ V\{s} is given by where X i ∈ C is the transmitted signal at node i, subject to an average power constraint E(|X i | 2 ) 1, h ij ∈ C is the channel gain from node i to node j, and N j ∈ CN (0, 1) is additive white circularly symmetric complex Gaussian noise, independent for different j.
It has been show in [21, Theorem 2.1] that using lattice codes for transmission and quantization at the relays, all rates R between source {s} and destination {d} satisfying can be achieved, where Ω is a source-destination cut of the network and X Ω = {X i , i ∈ Ω} are i.i.d . CN (0, 1). In addition, the restriction to i.i.d. Gaussian input distributions is within |V| bits/s/Hz of the cut-set upper bound [2]. Therefore the rate achieved using lattice codes in the above result is within 2|V| bits/s/Hz of the capacity of the network. The following corollary is an immediate consequence of Theorem 1.
Due to Corollary 1 the minimization in (3) is the minimization of a submodular function and the resulting optimal value is within 2|V| of the capacity of the network. 5 .

C. Wireless erasure network
In [8] the authors introduce a special class of wireless networks, called wireless erasure networks. In these networks, a directed graph G = (V, E) defines the interconnections between nodes. To model the broadcast effect of wireless networks, the signals on all outgoing arcs of any given node are equal to each other. There is no interference among multiple arcs arriving at a given node in this model, and the signals on the various arcs are erased independently of each other. We assume binary transmitted signals at each node, i.e. X i ∈ {0, 1}, i ∈ V\{d}, but all the results can be extended to models with larger input alphabets. It has been shown in [8] that the capacity of the network is where ǫ ij is the probability of erasure when node i is sending information to node j. We show in the following theorem that F (Ω) is submodular.
5 Notice that I(X Ω ; Y Ω c |X Ω c ) = log det(I + HH † ) where H is the matrix of channel gains from nodes in Ω to nodes in Ω c and H † is the conjugate transpose of H. Therefore, it is easy to compute the capacity of each cut.
We used in (a) the independence among X i and the channel erasures, in (b) the fact that for X i ∼ Bernoulli(1/2), H(X i ) = 1, and in (c) the fact that for X i ∼ Bernoulli(1/2), H(X i |Y j = e, j ∈ Ω c ) = 1 and for independent erasures we have p(Y j = e, j ∈ Ω c ) = j∈Ω c ǫ ij . Theorem 1 can be applied to conclude that F (Ω) is submodular.

V. POWER OPTIMIZATION
In the previous section, for the Gaussian relay network model, we considered fixed power assignments to the different nodes in the network, and have shown that a constant gap approximation to the capacity can be efficiently computed. In many applications it is of interest to allocate the nodes' transmission powers to optimize a given objective. For example, in a network where the nodes are battery powered, it may be of interest to maximize the network lifetime while satisfying a baseline quality of service. Alternatively, it may be desirable to maximize the network throughput for a given total power budget. This total power budget may arise due to, e.g., a maximum system weight constraint which is dominated by the battery weight, or a total system cost, which may be heavily influenced by the cost of the batteries. Power allocation optimization may also naturally arise in situations where the channel gains, while known to all the nodes, slowly vary over time. In this case, it may be desirable to optimally allocate power for the current channel condition.
As before, we characterize communication rates in terms of cut-set capacities. We consider a model where the cut-set capacities are functions of the cuts and powers assigned to the nodes in the network: F (Ω, p) : V × R |V| → R, and we focus on the Gaussian model where this function depends explicitly on the power assignment, F (Ω, p) = I(X Ω ; Y Ω c |X Ω c ) = log det(I + H Ω P Ω H † Ω ) (5) where H Ω is the matrix of channel gains from nodes in Ω to nodes in Ω c , H † is the conjugate transpose of H, and P Ω is a diagonal matrix where the diagonal elements are the powers of the nodes in Ω.
We will show in Lemma 1 below that in the Gaussian case F (Ω, p) is a concave function of p. While the results of this section are stated and proved for the Gaussian model, we conjecture that similar results should hold for other models in which F (Ω, p) is a concave function of p.
For Gaussian relay networks, we show that the following optimization problem can be solved in polynomial time, minimize R,P,p for fixed constants µ 1 , µ 2 , R 0 , p max and P tot . In the rest of the section, we denote the feasible set of the optimization (6) by K.
We use the Ellipsoid method [14], [15] to show that the optimization (6) can be solved efficiently. We will use the following definitions and result. The reader is referred to [15] for more details.

Definition 1 (Polynomial computability). A family of optimization programs is polynomially computable if:
(i) for any instance of the program and any point x in the domain, the objective and its subgradient can be computed in polynomial time in the size of the instance. (ii) for a given measure of infeasibility Infeas(·), it should be possible to determine if Infeas(x) ε in polynomial time, and when this inequality is not satisfied, it should be possible to find in polynomial time a vector c such that c T x > c T y, ∀y : Infeas(y) ε.

Definition 2 (Polynomial growth). A family of optimization programs has polynomial growth if the objectives and the infeasibility measures as functions of points
x in the domain grow polynomially with x 1 .

Proposition 2 ([15, Theorem 5.3.1]). Let P be a family of convex optimization programs equipped with infeasibility measure Infeas(·). Assume that the family is polynomially computable with polynomial growth and with polynomially bounded feasible sets. Then P is polynomially solvable.
In order to use Proposition 2 we need to check that the optimization (6) is a convex program. Since the objective function is linear, we only need to check that the feasible set K is convex.
It is easy to check that the sets P 1 = {(R, P, p) : 0 p p max }, P 2 = {(R, P, p) : p i P }, P 3 = {(R, P, p) : R 0 R}, and P 4 = {(R, P, p) : P P tot } are also convex sets, and, as a result, Having proved that (6) is a convex program, in order to use Proposition 2 we need to check that the conditions of Definitions 1, 2, and 3 are satisfied. Part (i) of Definition 1 follows from the linearity of the objective in (6). For part (ii) of Definition 1 we specify an infeasibility measure Infeas(·) : R |V|+2 → R as follows 6 : The conditions of part (ii) of Definition 1 are verified in the following theorem.
Proof: Part (a) requires checking that each of the arguments of the max of (7) is smaller than or equal to ε in polynomial time. The first six terms are linear functions and can be easily computed. The last term can be compared to ε by performing a minimization of a submodular function, which as was shown in Section IV-B, can also be computed in polynomial time.
We focus on condition (b). In this case Infeas((R, P, p)) > ε, meaning that at least one of the arguments of the max of (7) is larger than ε. We consider each case separately.
If −p i > ε we set c = −e i+2 where e i has a one in the i th position and zeros everywhere else, which can be easily checked to satisfy the condition of part (b).
Having proved these preliminary results, we are ready to prove the main result of this section.
Theorem 5. The optimization in (6) can be solved in polynomial time on the size of the problem.
Proof: The proof uses Proposition 2, which requires verifying the convexity of the problem together with the conditions of polynomial computability, polynomial growth, and polynomial boundedness of the feasible set. Convexity was proved in Lemma 1, while polynomial computability was shown in Theorem 4. Polynomial growth follows from the fact that F (Ω, p) is the log of a polynomial on p with degree at most |V|, while the objective and remaining terms that define the infeasibility measure are linear on (R, P, p). Finally, to check that feasible set if polynomially bounded, we note that the feasible set is a subset of the hypercube . It follows that the feasible set is contained in the Euclidean ball centered at the origin with radius (R max , P tot , p max ) 2 , which can be easily checked to grow polynomially on the size of the problem.
The general optimization problem (6) can be specialized to a power minimization with a minimum rate constraint, and to a rate maximization with a total power constraint. Both problems can be solved in polynomial time, as stated in the following corollaries to Theorem 5.

Corollary 2.
The following power minimization problem can be solved in polynomial time.

VI. ALGORITHMS AND SIMULATIONS
In this section we study different algorithms and and provide simulation results regarding submodular function minimization and power allocation problems. In the first subsection we look at minimum norm algorithm for submodular function minimization and we show that this algorithm can find the approximation to capacity of layered Gaussian relay networks with more that 300 nodes in a couple of minutes. In the second subsection we propose a heuristic algorithm to find the optimum power allocation for Gaussian relay networks.

A. submodular function minimization
One approach to solve the submodular minimization problem due to Lovász is based on extension of the set function f : 2 V → R to a convex function g : [0, 1] |V| → R that agrees with f on the vertices of the hypercube [0, 1] |V| , with a guarantee that min A⊆V f (A) is equal to min x g(x) for x ∈ [0, 1] |V| . In this section we assume the normalization f (∅) = 0.
The Lovász extension g of any set function f can be defined as follows. For a given x ∈ [0, 1] |V| order the elements of A key result is that if f is submodular its Lovász extension g is a convex function [14], [11]. In addition, finding the minimum of the submodular function f over subsets of V is equivalent to finding the minimum of the convex function g in the hypercube [0, 1] |V| . The optimization can be done in polynomial time using Ellipsoid algorithm [14].
There are other algorithms with faster running time to solve the submodular minimization problem [16], [17], [20]. To the best of our knowledge, the running time of the fastest algorithm is in the order of O(n 5 α + n 6 ), where α is the time that the algorithms takes to compute f (A) for any subset A ⊆ V [20]. For ADT networks, Gaussian relay networks, and erasure networks, α is the time to compute: the rank of n × n matrices, the determinant of n × n matrices, and equation (4), respectively.
However, for networks of large size, a complexity of O(n 5 α + n 6 ) may still be computationally cumbersome. As a result, in these cases it is desirable to have faster algorithms. Recently, Fujishing [11], [12] showed that the minimization of any submodular function can be cast as a and the corresponding minimum norm optimization is Letting x * be the solution of this minimization, the set A * = {v i : x * (v i ) < 0} is the solution to min A f (A). Whether the above optimization problem can be solved in polynomial time is an open problem. However empirical studies [12] have shown that this algorithm has comparable or even faster running times than the other algorithms with polynomial time performance guarantees. In our specific setting, for layered Gaussian relay networks of size up to around 300 nodes with 4 nodes per layer, we were able to find the approximate capacity (cf. (3)) in order of minutes (see Figure 2). In order to solve the minimization (10) we used the Matlab package provided in [19].

B. Power allocation
In Section V we have shown that the Ellipsoid method can be used to solve the optimization in (8) in polynomial time.
While in theory this result shows that the optimization in (8) is tractable, in practice the Ellipsoid method has a number of shortcomings that limit its usability. On the one hand the running time of the Ellipsoid method can be large compared to alternative algorithms, and on the other hand, for high dimensional problems it has shown numerical instability. In this section, we propose a heuristic algorithm to solve the optimization in (8) and show that this algorithm converges to the right solution. While the running time can be exponential on the network size, we show through simulations that the algorithm often converges within a time proportional to the network size.
Our proposed algorithm, Algorithm 1 (see pseudo-code below), is based on the cutting plane methods [24] in convex optimization. The optimization in (8) contains exponentially many constraints of the form In Algorithm 1 we first find the min-cut corresponding to assigning maximum power to all nodes in the network: Then, we modify the optimization in (8) by replacing the constraint (11) with The resulting convex program can be easily solved since it contains few constraints. After optimization, we let p * be the optimum power allocation for the current set of constraints and set We iteratively add the constraint to our set of constraints, and solve the optimization again. We stop if the new constraint is already in the set of constraints.

Algorithm 1 Power minimization
Input: Channel gain matrix H, desired rate R, vector of nodes' power constraints p max . Output: Min-cut, Power assignment p * that achieves approximate to rate R with minimum sum of powers.
Since in each iteration the algorithm adds a new constraint to the constraint set and the number of constraints in (11) is finite, the algorithm is guaranteed to find the optimum power in a finite number of iterations, which can be exponential.
We used simulations to test the performance of the algorithm for networks of varying size n ranging from 10 to 40. For each n, we generated 300 random networks with channels gains drawn i.i.d. using a N (0, 1) distribution. We set the desired transmission rate R 0 = 4 in (8) and set a maximum power constraint in each node to p max = 100. The results are  shown in Figures 3-6, where the the vertical bars represent ±1 standard deviation around the mean computed over the 300 random networks. Figure 3 shows the number of iterations of Algorithm 1 as a function of the number of nodes in the network. We see that the number of iterations grows O(n 3/2 ) with the number of nodes.
In Figure 4, we present the simulation time as a function of the network size. We observe a running time that grows slower than O(n 6 ). The figure also shows that the power optimization of networks of 40 nodes completes in less than two hours on an Intel Xeon quad core CPU running at 2.33 GHz.
The minimum sum of powers to approximately achieve R 0 = 4 for networks of different size is presented in Figure 5 with a blue solid line. Interestingly, the plot shows that the minimum power concentrates around the mean. In addition, the figure shows diminishing returns in total power savings resulting from increasing the network size.   For the special case of the diamond relay network with N relays, [26] shows that a fraction k/(k + 1) of the network capacity can be approximately achieved by using only k relays. In our setting of a general Gaussian relay network and sumpower minimization, we are interested in investigating whether it is possible to remove a large fraction of the nodes from the network without significantly affecting its performance.
In order to determine which nodes to remove from the network, we solve the sum-power minimization and compare the optimal power allocation p * i of each node i to a threshold P th . All relays with p * i < P th are removed from the network. Let N be the set of nodes with p * i P th . We optimize the power allocation for the network with node set {s} ∪ N ∪ {d} and determine whether the problem is feasible for R 0 = 4. If the problem is infeasible, we enlarge N by adding more relay nodes in decreasing order of p * i , until the problem becomes feasible. Figure 5 shows in the red dashed curve the resulting minimum sum of powers obtained by setting P th = 1. Figure 6 shows the corresponding size of the set N . We observe in Figure 6 that the number of nodes with allocated power p * i exceeding P th = 1 (possibly including more relays to make the problem feasible) remains fairly constant as the size of the network n increases. This means that most of the power is allocated to a small subset of the nodes.
Removing the remaining nodes from the network and optimizing the power allocation again over the resulting simplified network results in the minimum total power plotted in Figure 5. This figure shows that even though the number of nodes in the simplified network remains approximately constant as n increases, the total power required to approximately achieve R 0 = 4 decreases with n. This is due to the fact that larger n allows to choose the best relays for the simplified network. There is some performance loss in terms of total power due to network simplification but this loss may be compensated by power savings arising from turning off some of the relays. While we have not modeled the power consumption of the clocks, CPU and other subsystems required to keep a relay active, in practice they may become comparable to the power consumed by radio transmissions. This makes network simplification very useful in practice.