Efficient Global Optimal Resource Allocation in Non-Orthogonal   Interference Networks

Bho Matthiesen; Eduard A. Jorswieck

arXiv:1812.07253·cs.IT·October 17, 2019

Efficient Global Optimal Resource Allocation in Non-Orthogonal Interference Networks

Bho Matthiesen, Eduard A. Jorswieck

PDF

1 Repo

TL;DR

This paper introduces a novel, efficient algorithm for solving complex non-convex resource allocation problems in interference networks by differentiating variable types, leading to significant computational speed-ups and improved numerical stability.

Contribution

The authors develop the SIT algorithm that exploits variable convexity, introduces {}-essential feasibility, and handles fractional objectives without Dinkelbach's method, outperforming existing algorithms.

Findings

01

Achieves four orders of magnitude faster computation than state-of-the-art methods.

02

Provides a robust approach to non-convex resource allocation with better numerical stability.

03

Effectively handles fractional objectives without iterative Dinkelbach's algorithm.

Abstract

Many resource allocation tasks are challenging global (i.e., non-convex) optimization problems. The main issue is that the computational complexity of these problems grows exponentially in the number of variables instead of polynomially as for many convex optimization problems. However, often the non-convexity stems only from a subset of variables. Conventional global optimization frameworks like monotonic optimization or DC programming treat all variables as global variables and require complicated, problem specific decomposition approaches to exploit the convexity in some variables. To overcome this challenge, we develop an easy-to-use algorithm that inherently differentiates between convex and non-convex variables, preserving the low computational complexity in the number of convex variables. Another issue with these widely used frameworks is that they may suffer from severe…

Tables2

Table 1. TABLE I: Mean and median run times of throughput maximization for SND and “traditional” SND . For “traditional” SND , both the SIT and the Polyblock algorithm (PA) [ 4 ] are employed.

		SNR
		0 dB	15 dB	30 dB
SND	Mean	$385.78$ s	$351.92$ s	$5, 335.48$ s
SND	Median	$75.51$ s	$29.16$ s	$2, 494.54$ s
“trad.” SND	Mean	$19.77$ s	$6.23$ s	$3, 217.21$ s
“trad.” SND	Median	$9.5$ s	$1.84$ s	$831.71$ s
PA	Mean	$1.92 \cdot 10^{5}$ s	$77, 670.36$ s	$10, 602.06$ s
PA	Median	$46, 210.42$ s	$7, 091.09$ s	$35.67$ s

Table 2. TABLE II: Mean and median run times of EE computation for “traditional” SND and different solvers, all with precision η = 0.01 𝜂 0.01 \eta=0.01

		SNR
		0 dB	20 dB	40 dB
Gurobi	Mean	$5.1438$ s	$0.1771$ s	$0.155$ s
Gurobi	Median	$3.2781$ s	$0.0762$ s	$0.06$ s
Mosek	Mean	$15.0453$ s	$0.1368$ s	$0.0323$ s
Mosek	Median	$10.2756$ s	$0.0219$ s	$0.019$ s
Dinkelbach	Mean	$377.1501$ s	$145.4181$ s	$36.969$ s
Dinkelbach	Median	$162.811$ s	$23.027$ s	$16.9229$ s

Equations102

x \in C max f (x)

x \in C max f (x)

⎩ ⎨ ⎧ (x, ξ) \in C max s. t. \frac{f ^{+} ( x , ξ )}{f ^{-} ( x , ξ )} g_{i}^{+} (x, ξ) - g_{i}^{-} (x) \leq 0, i = 1, 2, \dots, m,

⎩ ⎨ ⎧ (x, ξ) \in C max s. t. \frac{f ^{+} ( x , ξ )}{f ^{-} ( x , ξ )} g_{i}^{+} (x, ξ) - g_{i}^{-} (x) \leq 0, i = 1, 2, \dots, m,

⎩ ⎨ ⎧ p, R max s. t. f (p, R) a_{i}^{T} R \leq lo g (1 + \frac{b _{i}^{T} p}{c _{i}^{T} p + σ _{i}}), i = 1, \dots, n R \geq 0, p \in [0, P]

⎩ ⎨ ⎧ p, R max s. t. f (p, R) a_{i}^{T} R \leq lo g (1 + \frac{b _{i}^{T} p}{c _{i}^{T} p + σ _{i}}), i = 1, \dots, n R \geq 0, p \in [0, P]

a_{i}^{T} R + lo g (c_{i}^{T} p + σ_{i}) - lo g ((b_{i}^{T} + c_{i}^{T}) p + σ_{i}) \leq 0.

a_{i}^{T} R + lo g (c_{i}^{T} p + σ_{i}) - lo g ((b_{i}^{T} + c_{i}^{T}) p + σ_{i}) \leq 0.

g_{i}^{+} (p, R)

g_{i}^{+} (p, R)

g_{i}^{-} (p)

⎩ ⎨ ⎧ p, R, t max s. t. t a_{i}^{T} R \leq lo g (1 + \frac{b _{i}^{T} p}{c _{i}^{T} p + σ _{i}}), i = 1, \dots, n R \geq t w, p \in [0, P]

⎩ ⎨ ⎧ p, R, t max s. t. t a_{i}^{T} R \leq lo g (1 + \frac{b _{i}^{T} p}{c _{i}^{T} p + σ _{i}}), i = 1, \dots, n R \geq t w, p \in [0, P]

GEE = \frac{\sum _{k} R _{k}}{ϕ ^{T} p + P _{c}},

GEE = \frac{\sum _{k} R _{k}}{ϕ ^{T} p + P _{c}},

g_{i}^{+} (p, R)

g_{i}^{+} (p, R)

g_{i}^{-} (p)

x \in [a, b] max f (x) s. t. g_{i} (x) \leq 0, i = 1, 2, \dots, m

x \in [a, b] max f (x) s. t. g_{i} (x) \leq 0, i = 1, 2, \dots, m

x \in [a, b] max f (x) s. t. g_{i} (x) \leq ε, i = 1, 2, \dots, m,

x \in [a, b] max f (x) s. t. g_{i} (x) \leq ε, i = 1, 2, \dots, m,

f (x^{*}) + η \geq sup {f (x) ∣ x \in [a, b], \forall i : g_{i} (x) \leq - ε},

f (x^{*}) + η \geq sup {f (x) ∣ x \in [a, b], \forall i : g_{i} (x) \leq - ε},

x \in [a, b] min i = 1, 2, \dots, m max g_{i} (x) s. t. f (x) \geq γ

x \in [a, b] min i = 1, 2, \dots, m max g_{i} (x) s. t. f (x) \geq γ

⎩ ⎨ ⎧ (x, ξ) \in C min s. t. i = 1, 2, \dots, m max (g_{i}^{+} (x, ξ) - g_{i}^{-} (x)) \frac{f ^{+} ( x , ξ )}{f ^{-} ( x , ξ )} \geq γ

⎩ ⎨ ⎧ (x, ξ) \in C min s. t. i = 1, 2, \dots, m max (g_{i}^{+} (x, ξ) - g_{i}^{-} (x)) \frac{f ^{+} ( x , ξ )}{f ^{-} ( x , ξ )} \geq γ

{(x, ξ) \in C min s. t. i = 1, 2, \dots, m max (g_{i}^{+} (x, ξ) - g_{i}^{-} (x)) γ f^{-} (x, ξ) - f^{+} (x, ξ) \leq 0

{(x, ξ) \in C min s. t. i = 1, 2, \dots, m max (g_{i}^{+} (x, ξ) - g_{i}^{-} (x)) γ f^{-} (x, ξ) - f^{+} (x, ξ) \leq 0

M_{-}^{k}

M_{-}^{k}

M_{+}^{k}

(x^{k}, ξ^{k}) \in D, ξ \in D_{y^{k}} min g (y^{k}, ξ) - β (M^{k}) \to 0 as x^{k} - y^{k} \to 0

(x^{k}, ξ^{k}) \in D, ξ \in D_{y^{k}} min g (y^{k}, ξ) - β (M^{k}) \to 0 as x^{k} - y^{k} \to 0

g (x, ξ) = max_{i = 1, 2, \dots, m} (g_{i}^{+} (x, ξ) - g_{i}^{-} (x)) .

g (x, ξ) = max_{i = 1, 2, \dots, m} (g_{i}^{+} (x, ξ) - g_{i}^{-} (x)) .

i = 1, 2, \dots, m max {g_{i}^{+} (x, ξ) - g_{i}^{-} (\overset{ˉ}{x}_{M}^{*})}

i = 1, 2, \dots, m max {g_{i}^{+} (x, ξ) - g_{i}^{-} (\overset{ˉ}{x}_{M}^{*})}

⎩ ⎨ ⎧ x, ξ min s. t. i = 1, 2, \dots, m max {g_{i}^{+} (x, ξ) - g_{i}^{-} (\overset{ˉ}{x}_{M_{i}}^{*})} γ f^{-} (x, ξ) - f^{+} (x, ξ) \leq 0 (x, ξ) \in C, x \in M_{i} .

⎩ ⎨ ⎧ x, ξ min s. t. i = 1, 2, \dots, m max {g_{i}^{+} (x, ξ) - g_{i}^{-} (\overset{ˉ}{x}_{M_{i}}^{*})} γ f^{-} (x, ξ) - f^{+} (x, ξ) \leq 0 (x, ξ) \in C, x \in M_{i} .

⎩ ⎨ ⎧ x, ξ min s. t. i = 1, 2, \dots, m max {g_{i, ξ}^{+} (ξ) + g_{i, x}^{+} (x) - g_{i}^{-} (\overset{ˉ}{x}_{M_{i}}^{*})} γ f_{ξ}^{-} (ξ) - f_{ξ}^{+} (ξ) + γ f_{x}^{-} (x) - f_{x}^{+} (x) \leq 0 ξ \in Ξ, x \in X \cap M_{i} .

⎩ ⎨ ⎧ x, ξ min s. t. i = 1, 2, \dots, m max {g_{i, ξ}^{+} (ξ) + g_{i, x}^{+} (x) - g_{i}^{-} (\overset{ˉ}{x}_{M_{i}}^{*})} γ f_{ξ}^{-} (ξ) - f_{ξ}^{+} (ξ) + γ f_{x}^{-} (x) - f_{x}^{+} (x) \leq 0 ξ \in Ξ, x \in X \cap M_{i} .

\left\{\begin{aligned} &\underset{\bm{\xi}}{\text{min}}&&\max_{i=1,2,\dots,m}\left\{g_{i,\xi}^{+}(\bm{\xi})+g_{i,x}^{+}(\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}})-g_{i}^{-}(\bar{\bm{x}}^{\ast}_{\mathcal{M}_{i}})\right\}\\ &\text{s.\,t.}&&\gamma f_{\xi}^{-}(\bm{\xi})-f_{\xi}^{+}(\bm{\xi})+\gamma f_{x}^{-}(\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}})-f_{x}^{+}(\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}})\leq 0\\ &&&\bm{\xi}\in\Xi.\end{aligned}\right.

\left\{\begin{aligned} &\underset{\bm{\xi}}{\text{min}}&&\max_{i=1,2,\dots,m}\left\{g_{i,\xi}^{+}(\bm{\xi})+g_{i,x}^{+}(\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}})-g_{i}^{-}(\bar{\bm{x}}^{\ast}_{\mathcal{M}_{i}})\right\}\\ &\text{s.\,t.}&&\gamma f_{\xi}^{-}(\bm{\xi})-f_{\xi}^{+}(\bm{\xi})+\gamma f_{x}^{-}(\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}})-f_{x}^{+}(\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}})\leq 0\\ &&&\bm{\xi}\in\Xi.\end{aligned}\right.

p_{i}^{0} = (x, ξ) \in C min x_{i} q_{i}^{0} = (x, ξ) \in C max x_{i} .

p_{i}^{0} = (x, ξ) \in C min x_{i} q_{i}^{0} = (x, ξ) \in C max x_{i} .

{find s. t. ξ \in C_{x^{k}} g_{i}^{+} (x^{k}, ξ) - g_{i}^{-} (x^{k}) \leq 0, i = 1, \dots, m .

{find s. t. ξ \in C_{x^{k}} g_{i}^{+} (x^{k}, ξ) - g_{i}^{-} (x^{k}) \leq 0, i = 1, \dots, m .

⎩ ⎨ ⎧ ξ \in C_{x^{k}} min s. t. \frac{f ^{+} ( x ^{k} , ξ )}{f ^{-} ( x ^{k} , ξ )} g_{i}^{+} (x^{k}, ξ) - g_{i}^{-} (x^{k}) \leq 0, i = 1, 2 \dots, m .

⎩ ⎨ ⎧ ξ \in C_{x^{k}} min s. t. \frac{f ^{+} ( x ^{k} , ξ )}{f ^{-} ( x ^{k} , ξ )} g_{i}^{+} (x^{k}, ξ) - g_{i}^{-} (x^{k}) \leq 0, i = 1, 2 \dots, m .

\left\{\begin{aligned} &\underset{(\bm{x},\bm{\xi})\in\mathcal{C}}{\text{max}}&&\min_{j}\bigg{\{}\frac{f_{j}^{+}(\bm{x},\bm{\xi})}{f_{j}^{-}(\bm{x},\bm{\xi})}\bigg{\}}\\ &\text{s.\,t.}&&g_{i}^{+}(\bm{x},\bm{\xi})-g_{i}^{-}(\bm{x})\leq 0,\quad i=1,2,\dots,m.\end{aligned}\right.

\left\{\begin{aligned} &\underset{(\bm{x},\bm{\xi})\in\mathcal{C}}{\text{max}}&&\min_{j}\bigg{\{}\frac{f_{j}^{+}(\bm{x},\bm{\xi})}{f_{j}^{-}(\bm{x},\bm{\xi})}\bigg{\}}\\ &\text{s.\,t.}&&g_{i}^{+}(\bm{x},\bm{\xi})-g_{i}^{-}(\bm{x})\leq 0,\quad i=1,2,\dots,m.\end{aligned}\right.

⎩ ⎨ ⎧ (x, ξ) \in C min s. t. i = 1, 2, \dots, m max (g_{i}^{+} (x, ξ) - g_{i}^{-} (x)) γ f_{j}^{-} (x, ξ) - f_{j}^{+} (x, ξ) \leq 0, j = 1, 2, \dots

⎩ ⎨ ⎧ (x, ξ) \in C min s. t. i = 1, 2, \dots, m max (g_{i}^{+} (x, ξ) - g_{i}^{-} (x)) γ f_{j}^{-} (x, ξ) - f_{j}^{+} (x, ξ) \leq 0, j = 1, 2, \dots

R_{k} \leq lo g (1 + \frac{∣ h _{k} ∣ ^{2} S _{k}}{γ _{k} ( S )})

R_{k} \leq lo g (1 + \frac{∣ h _{k} ∣ ^{2} S _{k}}{γ _{k} ( S )})

R_{k}

R_{k}

R_{k} + R_{l (k)}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Code & Models

Repositories

bmatthiesen/efficient-global-opt
noneOfficial

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Efficient Global Optimal Resource Allocation in Non-Orthogonal Interference Networks

Bho Matthiesen, and Eduard A. Jorswieck

Parts of this paper were presented at the 19th IEEE International Workshop on Signal Processing Advances in Wireless Communications, 2018 [1], and the 2019 IEEE International Conference on Acoustics, Speech and Signal Processing [2].

The authors are with the Chair for Communications Theory, Communications Laboratory, Technische Universität Dresden, Dresden, Germany (e-mail: [email protected], [email protected]).

This work is supported in part by the German Research Foundation (DFG) in the Collaborative Research Center 912 “Highly Adaptive Energy-Efficient Computing,” and under grant number JO 801/24-1.

Abstract

Many resource allocation tasks are challenging global (i.e., non-convex) optimization problems. The main issue is that the computational complexity of these problems grows exponentially in the number of variables instead of polynomially as for many convex optimization problems. However, often the non-convexity stems only from a subset of variables. Conventional global optimization frameworks like monotonic optimization or DC programming [3] treat all variables as global variables and require complicated, problem specific decomposition approaches to exploit the convexity in some variables [4]. To overcome this challenge, we develop an easy-to-use algorithm that inherently differentiates between convex and non-convex variables, preserving the low computational complexity in the number of convex variables. Another issue with these widely used frameworks is that they may suffer from severe numerical problems. We discuss this issue in detail and provide a clear motivating example. The solution to this problem is to replace the traditional approach of finding an $\varepsilon$ /̄approximate solution by the novel concept of $\varepsilon$ /̄essential feasibility. The underlying algorithmic approach is called successive incumbent transcending (SIT) algorithm and builds the foundation of our developed algorithm. A further highlight of this algorithm is that it inherently treats fractional objectives making the use of Dinkelbach’s iterative algorithm obsolete. Numerical experiments show a speed-up of four orders of magnitude over state-of-the-art algorithms and almost three orders of magnitude of additional speed-up over Dinkelbach’s algorithm for fractional programs.

Index Terms:

Resource allocation, global optimization, successive incumbent transcending, essential feasibility, multi-way relay channel, simultaneous non-unique decoding, interference networks

I Introduction

Resource allocation is essential in most communication systems [5]. Practical systems usually use algorithms with no or only weak optimality guarantees for performance reasons. Nevertheless, asserting the quality of these algorithms requires the knowledge of the optimal solution to these problems. The general optimization problem

[TABLE]

with $\mathcal{C}\in\mathds{R}^{n}$ and $f:\mathds{R}^{n}\mapsto\mathds{R}$ covers a large class of resource allocation problems. A point $\bm{x}^{\ast}\in\mathcal{C}$ satisfying $f(\bm{x}^{\ast})\geq f(\bm{x})$ for all $\bm{x}\in\mathcal{C}$ is called a global maximizer of $f$ . If $\bm{x}^{\ast}$ only satisfies this condition for all $\bm{x}$ in an open $\varepsilon$ -neighborhood of $\bm{x}^{\ast}$ for some $\varepsilon>0$ , i.e., for all $\bm{x}\in\{\bm{x}\in\mathds{R}^{n}:\left\lVert\bm{x}-\bm{x}^{\ast}\right\rVert<\varepsilon\}\cap\mathcal{C}$ , it is called a local minimizer. The difficulty in obtaining a global optimal solution to P1 is that all algorithms with polynomial computational complexity can at most obtain a local optimal solution. So, unless P1 belongs to the class of optimization problems with the property that every local maximum is a global maximum,111An important example are convex optimization problems where the objective of P1 is a concave function and $\mathcal{C}$ is a convex set. solving P1 has exponential computational complexity [6].

As an example, consider allocating the transmit power in an interference network. Albeit the capacity region of such an network is not known in general, the optimal decoder under the assumption of random codebooks222The random codebooks are restricted to superposition coding and time sharing. is known to be simultaneous non-unique decoding (SND) [7]. This leads to a global optimization problem where optimization is done jointly over the rates and powers. A close examination of this problem reveals that it is linear in the rate variables, i.e., for fixed power variables the problem can be solved in polynomial time [8, 9, 4]. Hence, the power variables are the only reason that the optimization problem is global and has exponential complexity. We call these variables global variables, while the remaining ones are named non-global.333A more precise definition is given in Section I-B. The most popular solution approaches for global resource allocation problems are monotonic optimization and DC programming. Both frameworks treat all variables as global variables which often results in unnecessary high numerical complexity. Moreover, transforming typical resource allocation problems to fit into these frameworks often requires auxiliary variables, which, of course, further increases computational complexity. Instead, in this paper we present a novel framework that preserves the computational complexity of the non-global variables and does not require any auxiliary variables.

Another often neglected issue with these algorithms is the assumption of a robust feasible set, i.e., a set with no isolated points. If this assumption does not hold, which might be the case for resource allocation problems, it leads to serious numerical problems. We avoid this problem entirely by using robust global optimization444There are, at least, two different meanings of “robust optimization:” the one discussed here that is robust against the effects of non-robust feasible sets and small changes in the tolerances, and the one that provides robustness against uncertainty in the input data [10], e.g., robust beamforming [11], or robust monotonic optimization [12]. [13, 14, 3]. The core idea is to shrink the feasible set by an infinitesimal amount and then solve a sequence of feasibility problems with a branch and bound (BB) procedure. This approach is called successive incumbent transcending (SIT) scheme and does not require any assumptions on the robustness of the feasible set $\mathcal{C}$ because it is designed to operate only on the accumulation points of $\mathcal{C}$ . The result is a numerically much stabler procedure than could be obtained using classical monotonic or DC programming algorithms. Moreover, the SIT approach always provides a good feasible solution even if stopped prematurely. Instead, conventional algorithms usually outer approximate the solution rendering intermediate solutions almost useless (because they are infeasbile).

Fractional objectives, which occur, e.g., in the optimization of the energy efficiency (EE) , can not be handled directly by monotonic optimization or DC programming. Instead, Dinkelbach’s algorithm is used where the original problem is transformed into an auxiliary problem which is then solved several times with one of these frameworks [15]. However, this approach has several drawbacks. First, convergence to the optimal solution of the original problem is only guaranteed if the auxiliary problem is solved exactly. In practice, this algorithm also works well for approximate solutions but the numerical accuracy should be sufficiently high. Second, the auxiliary problem needs to be solved several times, and, finally, the stopping criterion is unrelated to the distance of the obtained approximate optimal value to the true optimum. Especially the first two are critical for global optimization since they increase the computation time significantly. Instead, our framework is able to deal directly with fractional objectives avoiding these problems entirely.

Related work

The SIT approach was developed by Hoang Tuy in [13, 14, 3] and, to the best of our knowledge, has not been adopted for resource allocation problems yet. However, the importance of robust feasible sets has been noted in [16] where beamforming in a cognitive radio network is solved with DC programming. In [17], the basic principle of the SIT approach is used to solve a multi-objective optimization problem (MOP) .

Instead, decomposition approaches[18] are widely used. For example, in [19] the design of linear transceivers for multicarrier multiple-input multiple-output (MIMO) channels is considered. This challenging non-convex problem is solved by primal decomposition into a convex outer problem and inner problems with closed-form solutions. The authors of [20] combine successive convex approximation and primal decomposition to solve the sum rate maximization problem with Quality of Service (QoS) constraints for interfering broadcast channels with first order optimality. A distributed algorithm for coordinated beamforming in multicell multigroup multicast systems is developed in [21] based on primal decomposition and semidefinite relaxation. In [22] the partly convex-monotone structure of utility maximization problems in broadcast and interference channels is exploited via a BB procedure where branching is only performed over the global variables. This approach is similar to our framework but the problem setting is more specific and the optimization is over a convex set, i.e., the SIT approach is not needed. In [4], we solve a special case of the problems considered here by decomposing it into an inner linear and an outer monotonic program which is solved by the Polyblock algorithm [23].

Resource allocation for interference networks is mostly done under orthogonality assumptions to fall back to the noise-limited case. For example, in [24] monotonic optimization was first used to maximize the throughput in an interference network where interference is treated as noise. In [25] a monotonic optimization based framework for resource allocation in coordinated multi-cell systems is presented. Optimization of the EE in interference networks under the assumption that interference is treated as noise is considered in [15, 26]. In [27] energy-efficient coordinated beamforming in multi-cell, multi-user systems is considered under the assumption of realistic power consumption models. Energy-efficient resource allocation in OFDMA systems with and without wireless power transfer is studied in [28] and [29], respectively. This is, naturally, just a very incomplete list of papers dealing with resource allocation problems in interference networks under orthogonality constraints. Of course, the reasons for making orthogonality assumptions are manifold. One reason is surely that they are considerably easier to optimize than more involved systems relying on non-orthogonal communication schemes (e.g. SND ). The algorithm we develop here is well suited for exactly these more involved resource allocation problems. For further references on resource allocation problems please refer to Section II-A.

I-A Major Contributions

The key contributions of this paper are the following:

•

The concept of essential feasibility is introduced and its relevance for resource allocation problems is shown with a simple example. The numerical problems of conventional global optimization algorithms are discussed and the SIT scheme is introduced as a remedy for the these issues.

•

We extend the SIT algorithm as developed in [3] into an optimization framework able to deal with fractional objectives, non-convex variables, and multiple constraints. Specifically, we design an adaptive BB procedure that only branches over the global variables and leverages the power of commercially available state-of-the-art solvers to handle the non-global variables. To this end, we identify two different sets of sufficient conditions on the optimization problem. Fractional objectives are directly treated within the developed algorithm making the use of Dinkelbach’s iterative algorithm obsolete.

•

The framework is applied to optimize the throughput and energy efficiency of a multi-way relay channel (MWRC) . Extensive numerical results show a speed-up over state-of-the-art algorithms of four orders of magnitude for throughput maximization and almost three orders of magnitude over Dinkelbach’s Algorithm for EE maximization.

•

Code and data are made publicly available on GitHub [30]. This allows other researchers to easily verify our results and adapt the SIT algorithm for their own research.

I-B Notation & Preliminaries

A vector $\bm{x}\in\mathds{R}^{n}$ with components $(x_{1},\dots,x_{n})$ is said to dominate another vector $\bm{y}\in\mathds{R}^{n}$ , i.e., $\bm{y}\leq\bm{x}$ , if $y_{i}\leq x_{i}$ for all $i=1,\dots,n$ . For $\bm{a}\leq\bm{b}$ , the set $[\bm{a},\bm{b}]=\{\bm{x}\,|\,\bm{a}\leq\bm{x}\leq\bm{b}\}$ is called a box. A function $f:\mathds{R}^{n}_{\geq 0}\mapsto\mathds{R}$ is increasing if $f(\bm{x}^{\prime})\leq f(\bm{x})$ whenever $\bm{x}^{\prime}\leq\bm{x}$ , and decreasing if $-f$ is increasing. It is called mixed monotonic if it is increasing in the variables $(x_{i})_{i\in\mathcal{I}}$ and decreasing in $(x_{i})_{i\in\{1,2,\ldots,n\}\setminus\mathcal{I}}$ for some index set $\mathcal{I}$ . The functions $f_{1}(\bm{x}),\dots,f_{n}(\bm{x})$ are called jointly mixed monotonic if all functions are mixed monotonic with respect to the same index set $\mathcal{I}$ . A common minimizer (maximizer) of the functions $f_{1}(\bm{x}),\dots,f_{n}(\bm{x})$ over the set $\mathcal{X}$ is any $\bm{x}^{\ast}$ that satisfies $\bm{x}^{\ast}\in\bigcap_{i=1}^{n}\operatorname*{arg\,min}_{\bm{x}\in\mathcal{X}}f_{i}(\bm{x})$ ( $\bm{x}^{\ast}\in\bigcap_{i=1}^{n}\operatorname*{arg\,max}_{\bm{x}\in\mathcal{X}}f_{i}(\bm{x})$ ).

A set $\mathcal{G}\subseteq\mathds{R}^{n}_{\geq 0}$ is said to be normal if for $0\leq\bm{x}^{\prime}\leq\bm{x}$ , $\bm{x}\in\mathcal{G}$ $\Rightarrow$ $\bm{x}^{\prime}\in\mathcal{G}$ , and normal in a box $[\bm{a},\bm{b}]$ if the previous implication only holds for $\bm{a}\leq\bm{x}^{\prime}\leq\bm{x}\leq\bm{b}$ . A set $\mathcal{H}\subseteq\mathds{R}^{n}_{\geq 0}$ is called conormal if $\bm{x}+\mathds{R}_{\geq 0}^{n}\subseteq\mathcal{H}$ whenever $\bm{x}\in\mathcal{H}$ , and conormal in a box $[\bm{a},\bm{b}]$ if for $\bm{b}\geq\bm{x}^{\prime}\geq\bm{x}\geq\bm{a}$ , $\bm{x}\in\mathcal{H}\Rightarrow\bm{x}^{\prime}\in\mathcal{H}$ [3, Sec. 11.1.1]. Let $\mathcal{A}\subseteq\mathds{R}^{n}$ and $(\bm{x},\bm{y})\in\mathcal{A}$ . Then, $\operatorname{proj}_{\bm{x}}\mathcal{A}=\{\bm{x}\,|\,(\bm{x},\bm{y})\in\mathcal{C}\ \text{for some}\ \bm{y}\}$ , i.e., the projection of $\mathcal{C}$ onto the $\bm{x}$ coordinates; $\operatorname{diam}\mathcal{A}$ is the diameter of $\mathcal{A}$ , i.e., the maximum distance between two points in $\mathcal{A}$ ; and $\mathcal{A}_{\tilde{\bm{x}}}=\{\bm{y}|(\tilde{\bm{x}},\bm{y})\in\mathcal{D}\}$ is called the $\tilde{\bm{x}}$ -section of $\mathcal{A}$ .

Finally, consider P1 and let its optimal value be $v\eqref{opt}$ . To make the previously introduced notion of global and non-global variables more precise, split $\bm{x}$ into two vectors $\bm{y}=(x_{i})_{i\in\mathcal{I}}$ and $\bm{z}=(x_{i})_{i\in\{1,2,\ldots,n\}\setminus\mathcal{I}}$ and consider a modified version of P1 where we only optimize over $\bm{y}$ for some fixed $\bm{z}$ , i.e., $\max_{\bm{y}\in\mathcal{C}_{\bm{z}}}f(\bm{y},\bm{z})$ . If there exists an algorithm to solve this optimization problem with computational complexity significantly less than required for solving the global part of P1,555Ideally, the computational complexity for solving the non-global part of P1 grows polynomially in the number of non-global variables. This is, for example, the case if P1 is linear in $\bm{y}$ [8, 9]. It also holds in many cases were P1 is convex in $\bm{y}$ [31]. However, there exist convex optimization problems that are not solvable in polynomial time [32, Sect. 5.4] but still have computational complexity significantly less than general non-convex optimization problems [33]. Nevertheless, it might be helpful to think of the terms “global” and “non-global” variables as synonyms for “non-convex” and “convex” variables, respectively. the variables $\bm{y}$ and $\bm{z}$ are denoted as non-global and global variables, respectively.

I-C Outline

The remaining part of this paper is organized as follows. In the next section, we formally state the considered optimization problem and motivate it in the context of resource allocation problems in wireless interference networks. Section III introduces important mathematical preliminaries including the SIT scheme that forms the basis of the proposed algorithm. These developed principles are then applied in Section IV to construct Algorithm 2, one of the main contributions of this paper. In Section V, we apply the developed framework to a specific resource allocation problem that is used to benchmark our method against the state-of-the-art. Finally, we give our conclusions in Section VI.

II Problem statement

We consider the following global optimization problem

[TABLE]

with global variables $\bm{x}$ and non-global variables $\bm{\xi}$ . The functions $\{g_{i}^{-}(\bm{x})\}$ are required to have a common maximizer over every box $[\underaccent{\bar}{\vec x},\bar{\bm{x}}]\subseteq\mathcal{M}_{0}$ with $\mathcal{M}_{0}$ being a box enclosing the $\bm{x}$ dimensions of $\mathcal{C}$ , i.e., $\mathcal{M}^{0}\supset\operatorname{proj}_{\bm{x}}\mathcal{C}$ . A sufficient condition for the existence of this common maximizer is that the functions $g_{i}^{-}(\bm{x})$ are jointly mixed monotonic. This includes the case where all functions are either increasing or decreasing. Further, we assume the functions $f^{-}(\bm{x},\bm{\xi})$ , $g_{i}^{+}(\bm{x},\bm{\xi})$ , $i=1,\dots,m$ , to be lower semi-continuous (l.s.c.) , the functions $f^{+}(\bm{x},\bm{\xi})$ , $g_{i}^{-}(\bm{x})$ , $i=1,\dots,m$ , to be upper semi-continuous (u.s.c.) , and, without loss of generality (w.l.o.g.) , $f^{-}(\bm{x},\bm{\xi})>0$ .

The goal of this paper is to design a numerically stable BB procedure to solve P2 that preserves the computational complexity in the non-global variables $\bm{\xi}$ . This requires additional assumptions to those stated above. Specifically, we have identified two different sets of technical requirements that are stated below. Both of these cases contain conditions that depend on a constant $\gamma$ which will hold the current best known value in the developed algorithm. We will discuss the domain of $\gamma$ after the definition of both cases below.

Case A (DC problems)

If $\mathcal{C}$ is a closed convex set and $\gamma f^{-}(\bm{x},\bm{\xi})-f^{+}(\bm{x},\bm{\xi})$ , $g_{1}^{+}(\bm{x},\bm{\xi})$ , …, $g_{m}^{+}(\bm{x},\bm{\xi})$ are jointly convex in $(\bm{x},\bm{\xi})$ for all $\gamma$ , problem P2 resembles a DC optimization problem but with fractional objective and additional non-DC variables. $\lozenge$

Case B (Separable problems)

Let $\mathcal{C}=\mathcal{X}\times\Xi$ such that $\bm{x}\in\mathcal{X}$ and $\bm{\xi}\in\Xi$ with $\Xi$ being a closed convex set, and let each function of $(\bm{x},\bm{\xi})$ be separable in the sense that $h(\bm{x},\bm{\xi})=h_{x}(\bm{x})+h_{\xi}(\bm{\xi})$ . Further, let the functions $\gamma f^{-}_{\xi}(\bm{\xi})-f^{+}_{\xi}(\bm{\xi})$ , $g_{1,\xi}^{+}(\bm{\xi})$ , …, $g_{m,\xi}^{+}(\bm{\xi})$ be convex in $\bm{\xi}$ for all $\gamma$ , and let the functions $\gamma f^{-}_{x}(\bm{x})-f^{+}_{x}(\bm{x})$ , $g_{1,x}^{+}(\bm{x})$ , …, $g_{m,x}^{+}(\bm{x})$ have a common minimizer over $\mathcal{X}\cap\mathcal{M}$ for every box $\mathcal{M}\subseteq\mathcal{M}_{0}$ and all $\gamma$ . Finally, let the function $\gamma f^{-}_{x}(\bm{x})-f^{+}_{x}(\bm{x})$ be either increasing for all $\gamma$ with $\mathcal{X}$ being a closed normal set in some box, or decreasing for all $\gamma$ with $\mathcal{X}$ being a closed conormal set in some box. $\lozenge$

Remark 1

Separable problems often lead to linear auxiliary optimization problems with less variables instead of, typically, convex problems for Case A. They usually have lower computational complexity than DC problems and, thus, if the problem at hand falls into both cases, it is usually favorable to consider it as a separable problem.

With both cases defined we can continue our discussion of $\gamma$ . First, observe that $\gamma$ only appears as a factor to $f^{-}(\bm{x},\bm{\xi})$ . Thus, its value is only relevant if P2 is a fractional program, i.e., if $f^{-}(\bm{x},\bm{\xi})$ is not constant. In that case, the only relevant property of $\gamma$ is its sign and whether it may change during the algorithm. For example, in Case A the function $\gamma f^{-}_{\xi}(\bm{\xi})-f^{+}_{\xi}(\bm{\xi})$ is convex if $f^{+}_{\xi}(\bm{\xi})$ is concave and $\gamma f^{-}_{\xi}(\bm{\xi})$ is convex. The latter is the case if $\gamma\geq 0$ and $f^{-}_{\xi}(\bm{\xi})$ is convex, or if $\gamma\leq 0$ and $f^{-}_{\xi}(\bm{\xi})$ is concave. Thus, in most cases, we should ensure that the sign of $\gamma$ is constant. In general, $\gamma$ may take values between some $\gamma_{0}$ and $v(\ref{opt:genRA})+\eta$ for some small $\eta>0$ . The lower end of the range $\gamma_{0}$ is either the objective value of P2 for some preliminary known nonisolated feasible point $(\bm{x},\bm{\xi})$ or an arbitrary value satisfying $\gamma_{0}\leq\frac{f^{+}(\bm{x},\bm{\xi})}{f^{-}(\bm{x},\bm{\xi})}$ for all feasible $(\bm{x},\bm{\xi})$ . This implies, e.g., that $\gamma$ is non-negative if $f^{+}(\bm{x},\bm{\xi})$ is non-negative. Otherwise, it might be necessary to find a nonisolated feasible point such that $f^{+}(\bm{x},\bm{\xi})\geq 0$ or transform the problem.

II-A Application Example: Resource Allocation in Interference Channels

Determining an achievable rate region of a communication network usually involves two steps: first, characterizing the achievable rate region with information theoretical tools, and, second, finding Pareto-optimal resource allocations. The corresponding optimization problem in many Gaussian interference networks is

[TABLE]

for some performance function $f(\bm{p},\bm{R})$ and positive vectors $\bm{a}_{i},\bm{b}_{i},\bm{c}_{i}\geq 0$ , $i=1,\ldots,n$ . The optimization variables $\bm{p}$ are the allocated transmit powers and $\bm{R}$ are the achievable transmission rates for asymptotically error-free communication. Usually, the vectors $\bm{b}_{i}$ and $\bm{c}_{i}$ represent the effective channel gain, $\sigma_{i}$ is the variance of the Gaussian noise observed at receiver $i$ , $\bm{a}_{i}$ is a sparse vector where the non-zeros entries are small integers (mostly ones), and $\bm{P}$ are the maximum transmit powers. Applications of this model include multi-cell communication systems [25], heterogeneous dense small-cell networks [34], cognitive radio [35], and digital subscriber line (DSL) systems [36]. For the sake of simplicity, we assume the signal to interference plus noise ratio (SINR) in the right-hand side (RHS) expressions of the constraints to be linear fractions of the transmit powers. We note that this is not always the case and that our framework is not limited to this case.

The feasible set of LABEL:opt:exRA belongs to the class of considered problems. We identify the global variables as $\bm{p}$ and the non-globals as $\bm{R}$ . The rate constraints are equivalent to

[TABLE]

Since $\bm{a}_{i},\bm{b}_{i},\bm{c}_{i}\geq 0$ for all $i$ , the linear function $\bm{a}_{i}\bm{R}$ and the $\log$ -functions are increasing. Thus, we can identify $g^{+}_{i}(\bm{p},\bm{R})$ and $g^{-}_{i}(\bm{p})$ as:

[TABLE]

The functions $g_{i}^{-}(\bm{p})$ are decreasing, and, thus, are jointly maximized over the box $[\underaccent{\bar}{\vec p},\bar{\bm{p}}]$ by $\underaccent{\bar}{\vec p}$ . Further, $g^{+}_{i}(\bm{p},\bm{R})$ is separable in $\bm{p}$ and $\bm{R}$ with $g^{+}_{i,R}(\bm{R})$ being linear and $g^{+}_{i,p}(\bm{p})$ decreasing and convex in $\bm{p}$ . Thus, depending on $f(\bm{p},\bm{R})$ problem LABEL:opt:exRA qualifies for both, Cases A and B.

The complete Pareto boundary is characterized by all solutions to the MOP LABEL:opt:exRA with the vector objective $f(\bm{p},\bm{R})=[R_{1},R_{2},\ldots,R_{K}]$ [37]. Several approaches exist to transform this into a scalar optimization problem. Two very popular are the scalarization and the rate profile approach666This approach is also known as rate balancing [38]. [39].

Weighted Sum Rate

In the scalarization approach, the weighted sum of the objectives is optimized, i.e., $f(\bm{p},\bm{R})=f^{+}(\bm{p},\bm{R})\coloneqq\bm{w}^{T}\bm{R}$ . Obviously, with this objective, LABEL:opt:exRA belongs to the class of separable problems since Case B is satisfied. Varying the weights between 0 and 1 with $\sum_{k}w_{k}=1$ characterizes the convex hull of the Pareto boundary of the achievable rate region. This is by far the most widely used performance metric for power control in wireless communication systems due its clear operational meaning: with all weights $w_{k}=1$ , $f^{+}(\bm{p},\bm{R})$ is the total throughput in the network. Moreover, in several networks this condition leads to the characterization of the stability region. In this queueing theoretic setting, choosing the weights proportional to the queue lengths prioritizes longer queues and stabilizes the network [40, 41].

Many instances of the weighted sum rate maximization problem considered in the literature are either convex optimization problems [42, 43, 44], or have a rather good-natured rate region allowing to eliminate the rates $\bm{R}$ from the problem such that the resulting global optimization is only over the powers [24, 45, 25]. If neither is the case, variable reduction techniques are necessary to keep the computational complexity at a reasonable level. In [4], the problem at hand is transformed into a monotonic optimization problem [23] where the objective is a linear program with the transmit power as parameter. This approach is considerably slower, theoretically more involved, and less versatile than the proposed method.

Rate Profile Approach

The utility profile approach finds the intersection of a ray in the direction $\bm{w}$ and the Pareto boundary of the achievable rate region, i.e.,

[TABLE]

for some $\bm{w}\geq 0$ . W.l.o.g. one can assume $\left\lVert\bm{w}\right\rVert=1$ . By varying the direction of the ray, the complete Pareto boundary can be characterized. It is easily verified that LABEL:opt:exMaxmin satisfies Case B.

Problem LABEL:opt:exMaxmin is the epigraph form of LABEL:opt:exRA with $f(\bm{p},\bm{R})=\min_{k}\frac{R_{k}}{w_{k}}$ as long as $\bm{w}\neq\bm{0}$ . Since $f(\bm{p},\bm{R})$ is continuous and concave, it also satisfies Case B. Besides reducing the number of variables, this reformulation also makes it apparent that the rate profile approach is equivalent to the weighted max-min fairness performance function. For $\bm{w}=\bm{1}$ , it puts very strong emphasis on the weakest user and often results in low spectral efficiency. However, in some cases this approach might be computationally less challenging than the scalarization approach as we will discuss below.

Global Energy Efficiency

The global energy efficiency (GEE) is the most widely used metric to measure the network energy efficiency, a key performance metric in 5G and beyond networks [46, 26]. It is defined as the benefit-cost ratio of the total network throughput and the associated power consumption, i.e.,

[TABLE]

where $\bm{\phi}\geq 1$ are the inverses of the power amplifier efficiencies and $P_{c}$ is the total circuit power necessary to operate the network. Using the GEE as performance metric results in a fractional programming problem [47, 48, 49, 50]. Usually, it is solved iteratively with Dinkelbach’s Algorithm [48] where, in each iteration, an auxiliary optimization problem is solved. This auxiliary problem generally has very similar properties to the weighted sum rate maximization problem, i.e., for interference networks it is a global optimization problem with exponential complexity. This inner problem must be solved several times with high numerical accuracy for the convergence guarantees of Dinkelbach’s Algorithm to hold [15].

Instead, our proposed algorithm allows to solve these fractional programs directly resulting in significantly lower complexity. Consider the generic resource allocation problem LABEL:opt:exRA with objective 2. We identify $f^{+}(\bm{p},\bm{R})=\sum_{k}R_{k}$ and $f^{-}(\bm{p},\bm{R})=\bm{\phi}^{T}\bm{p}+P_{c}$ and observe that both are linear functions. Hence, $\gamma f^{-}(\bm{x},\bm{\xi})-f^{+}(\bm{x},\bm{\xi})$ is convex for all $\gamma$ and Case A is satisfied. Of course, $f^{+}(\bm{p},\bm{R})$ and $f^{-}(\bm{p},\bm{R})$ are also separable, and because the GEE is non-negative, $\gamma\geq 0$ and $\gamma f_{p}^{-}(\bm{p})$ is increasing. But since $g_{i,p}^{+}(\bm{p})$ is decreasing, the functions $\gamma f_{p}^{-}(\bm{p})$ , $g_{1,p}^{+}(\bm{p})$ , …, $g_{m,p}^{+}(\bm{p})$ do not have a common minimizer over $\mathcal{X}\cap\mathcal{M}$ for any box $\mathcal{M}$ and Case B does not apply. However, observe that instead of the choice in 1, we can also identify $g^{+}_{i}(\bm{p},\bm{R})$ and $g^{-}_{i}(\bm{p})$ as

[TABLE]

The functions $g^{+}_{i}(\bm{p},\bm{R})$ are separable in $\bm{p}$ and $\bm{R}$ , and $g^{+}_{i,p}(\bm{p})$ and $g^{-}_{i}(\bm{p})$ are increasing in $\bm{p}$ . Hence, $\gamma f_{p}^{-}(\bm{p})$ , $g_{1,p}^{+}(\bm{p})$ , …, $g_{m,p}^{+}(\bm{p})$ have a common minimizer over $\mathcal{X}\cap\mathcal{M}$ for every box $\mathcal{M}$ and Case B applies. In general, 1 will result in tighter bounds and, thus, faster convergence. On the other hand, Case B has significantly lower numerical complexity than Case A (cf. Remark 1) which should compensate for this drawback. However, a final assessment is only possible through numerical experimentation which will be carried out in Section V.

Energy Efficiency Region

Similar to the GEE , the individual energy efficiency of link $k$ is defined as the benefit-cost ratio of the link’s throughput divided by the link’s power consumption, i.e., $\mathrm{EE}_{k}=\frac{R_{k}}{\phi_{k}p_{k}+P_{c,k}}$ . Analogue to the achievable rate region, there is also an achievable EE region whose Pareto boundary is the solution to the MOP LABEL:opt:exRA with objective $f(\bm{p},\bm{R})=[\mathrm{EE}_{1},\ldots,\mathrm{EE}_{K}]$ [51]. The same approaches as discussed earlier can be used to transform this MOP into a scalar optimization problem. In this case, the rate profile (or rather utility profile) approach should be favored because the scalarization approach leads to a sum-of-ratios problem. The sum-of-ratios problem is one of the most difficult fractional programs and known to be essentially $\mathcal{NP}$ -hard [52, 53]. While for other EE problems (e.g., the GEE ) first-order optimal solutions are usually observed to be globally optimal [15] this does not hold for the sum-of-ratios case [34].

Instead, with the utility profile approach the objective in LABEL:opt:exRA is the weighted minimum EE $f(\bm{p},\bm{R})=\min_{k}\frac{\mathrm{EE}_{k}}{w_{k}}$ . This is, as is the sum-of-ratios problem, a multi-ratio optimization problem but numerically less challenging. Another benefit of using the utility profile approach over the scalarization approach is that the EE region is non-convex and the scalarization approach only obtains the convex hull of the Pareto region [15]. The general approach to solving such an optimization problem is to use the Generalized Dinkelbach’s Algorithm [54, 51, 15]. As in the GEE case, our proposed method solves this problem with much lower complexity. However, some minor modifications of the algorithm are necessary to deal with this kind of problem. Please refer to Section IV-C for more details.

Proportional Fairness

Another well known performance function is proportional fairness where the objective is the product of the rates, i.e., $\prod_{k}R_{k}$ [55]. The operating point achieved by this metric is usually almost as fair as the one obtained by max-min fairness but achieves significantly higher throughput. Observe that the objective is log-concave. Thus, we can determine the proportional fair operating point by solving LABEL:opt:exRA with objective $f(\bm{p},\bm{R})=\sum_{k}\log(R_{k})$ .

III Robust Approach to Global Optimization

In this section, we introduce some mathematical preliminaries essential for developing the proposed algorithm. This is done by means of a simpler problem than P2 to ease the exposition. Specifically, we consider the optimization problem

[TABLE]

where $f,g_{1},g_{2},\dots,g_{m}$ are non-convex continuous real-valued functions and $\bm{a},\bm{b}$ are real-valued vectors satisfying $\bm{a}\leq\bm{b}$ . This is a general non-convex optimization problem with possibly quite complicated feasible set.

III-A The Issue with $\varepsilon$ -Approximate Solutions

Most current solution methods for this problem are devised to compute a solution $\bar{\bm{x}}(\varepsilon)$ of the $\varepsilon$ -relaxed problem

[TABLE]

This solution $\bar{\bm{x}}(\varepsilon)$ is usually accepted as an $\varepsilon$ -approximate optimal solution to P5 since $\bar{\bm{x}}(\varepsilon)$ is almost feasible for small $\varepsilon>0$ and tends to a feasible solution as $\varepsilon\rightarrow 0$ . Since $f(\bar{\bm{x}}(\varepsilon))$ also tends to the optimal value $v(\ref{opt:generalNC})$ , $f(\bar{\bm{x}}(\varepsilon))$ should be close to $v(\ref{opt:generalNC})$ for a sufficiently small $\varepsilon=\varepsilon_{0}$ . The problem with this approach is that $\varepsilon_{0}$ is, in general, unknown and hard to determine. Thus, the obtained solution and $f(\bar{\bm{x}}(\varepsilon)$ ) can be quite far away from the true optimal value even for a small $\varepsilon$ [13]. This is apparent from the following example.

Example 1 (Issues with $\varepsilon$ -approximate solutions)

Consider a Gaussian 2-user single-input single-output (SISO) multiple-access channel (MAC) with channel gains $h_{1}$ and $h_{2}$ , transmit powers $p_{1}\leq P_{1}$ and $p_{2}\leq P_{2}$ , and a minimum total throughput of $Q$ , i.e., $\log_{2}(1+\left\lvert h_{1}\right\rvert^{2}p_{1}+\left\lvert h_{2}\right\rvert^{2}p_{2})\geq Q$ . Transmitter $i$ , $i=1,2$ , is eavesdropped by a single antenna adversary over a channel $g_{i}$ . The eavesdroppers are only able to overhear one of the transmitters and do not cooperate with each other [56]. The total information leakage to the eavesdroppers is limited by $L$ , i.e., $\log_{2}(1+\left\lvert g_{1}\right\rvert^{2}p_{1})+\log_{2}(1+\left\lvert g_{2}\right\rvert^{2}p_{2})\leq L$ . The transmit power of transmitter 1 should be minimized without violating these constraints. The resulting feasible set for $\left\lvert h_{1}\right\rvert^{2}=\left\lvert h_{2}\right\rvert^{2}=10$ , $\left\lvert g_{1}\right\rvert^{2}=\frac{1}{2}$ , $\left\lvert g_{2}\right\rvert^{2}=1$ , $Q=\log_{2}(61)$ , and $L=\log_{2}(8.99)$ is shown in Fig. 1. The true optimum solution is $\bm{p}^{\ast}=(4.00665,1.99335)$ with $f(\bm{p}^{\ast})=4.00665$ , while the $\varepsilon$ -approximate solution for $\varepsilon_{1}=10^{-3}$ is $\bar{\bm{p}}(\varepsilon_{1})=(0.995843,5)$ with $f(\bar{\bm{p}}(\varepsilon_{1}))=0.995843$ . This is, obviously, quite far away from both, the optimal solution and value. Instead, the $\varepsilon$ -approximate solution obtained for $\varepsilon_{2}=10^{-4}$ is $\bar{\bm{p}}(\varepsilon_{2})=(4.00541,1.99417)$ with $f(\bar{\bm{p}}(\varepsilon_{2}))=4.00541$ . $\lozenge$

Thus, even for well behaved problems the $\varepsilon$ -relaxation approach might fail. A slightly better approach to the approximate optimal solution of P5 is that of an $\eta$ -optimal solution. A feasible vector $\bar{\bm{x}}$ is called $\eta$ -optimal solution if it satisfies $f(\bar{\bm{x}})\geq f(\bm{x})-\eta$ for all feasible $\bm{x}$ . The downside of this approach is that it does not converge in finitely many steps if the optimal solution is an isolated feasible point, i.e., a point at the center of a ball containing no other feasible points [13].

Example 2 (Isolated optimal solution)

Consider 1 again. With $L=\log_{2}(9)$ there exists an isolated feasible point $\bm{p}=(1,5)$ (close to $\bar{\bm{p}}(\varepsilon_{1})$ in Fig. 1) that also happens to be the optimal solution. However, apart from the algorithmic difficulties in computing it, it might also be unstable under small pertubations of the data and is, thus, quite impractical from an engineering point of view. $\lozenge$

III-B $\varepsilon$ -Essential Feasibility

A common approach to this dilemma is to assume that the feasible set $\mathcal{F}$ is robust, i.e., it satisfies $\mathcal{F}^{\ast}=\operatorname{cl}(\operatorname{int}\mathcal{F})$ where $\operatorname{cl}$ and $\operatorname{int}$ denote the closure and interior, respectively. Unfortunately, this condition is generally very hard to check, so that, in practice, we have to deal with feasible sets where we do not know a priori whether they are robust or not.

This motivates the concept of $\varepsilon$ -essential optimality developed by Tuy [13, 14, 3]. A solution $\bm{x}^{\ast}\in\mathcal{F}^{\ast}$ is called essential optimal solution of P5 if $f(\bm{x}^{\ast})\geq f(\bm{x})$ for all $x\in\mathcal{F}^{\ast}$ . A point $\bm{x}\in[\bm{a},\bm{b}]$ satisfying $g_{i}(\bm{x})\leq-\varepsilon$ for all $i$ and some $\varepsilon>0$ is called $\varepsilon$ -essential feasible and a solution of P5 is said to be essential $(\varepsilon,\eta)$ -optimal if it satisfies

[TABLE]

for some $\eta>0$ . Clearly, for $\varepsilon,\eta\rightarrow 0$ an essential $(\varepsilon,\eta)$ -optimal solution is a nonisolated feasible point which is optimal.

III-C Successive Incumbent Transcending Scheme

The robust approach to global optimization employed here uses the SIT scheme in Algorithm 1 to generate a sequence of nonisolated feasible solutions converging to an essential optimal solution of P5. The core problem in the SIT scheme is, given a real number $\gamma$ , to check whether P5 has a nonisolated feasible solution $\bm{x}$ satisfying $f(\bm{x})\geq\gamma$ , or, else, establish that no such $\varepsilon$ -essential feasible $\bm{x}$ exists. Given that this subproblem is solved within finitely many steps, Algorithm 1 converges to the global optimal solution within finitely many iterations. Apart from the improved numerical stability and convergence, the SIT algorithm has another very desirable feature: it provides a good nonisolated feasible (but possibly suboptimal) point even if terminated prematurely. Instead, conventional algorithms usually outer approximate the solution rendering intermediate solutions almost useless.

Consider the optimization problem

[TABLE]

where we interchanged objective and constraints of P5. Very often $f(x)$ has, or could be modified easily to have, nice properties like being concave or increasing, such that the feasible set of P7 is nice, i.e., it is robust and a feasible point can be computed efficiently using an adaptive BB procedure [3, Prop. 6.2]. The following proposition, which is an adapted version of [3, Prop. 7.13], establishes a duality between P5 and P7 in the sense that the feasibility problem in Step 1 of Algorithm 1 is equivalent to solving P7.

Proposition 1

For every $\varepsilon>0$ , the $\varepsilon$ -essential optimal value of P5 is less than $\gamma$ if the optimal value of P7 is greater than $-\varepsilon$ .

Proof:

If the optimal value of P7 is greater than $-\varepsilon$ , then any $\bm{x}\in[\bm{a},\bm{b}]$ such that $g_{i}(\bm{x})\leq-\varepsilon$ , for all $i=1,2,\dots,m$ , must satisfy $f(\bm{x})<\gamma$ . Hence, by the compactness of the feasible set of P5, $\max\{f(\bm{x})|\bm{x}\in[\bm{a},\bm{b}],\ \forall i:g_{i}(\bm{x})\leq-\varepsilon\}<\gamma$ . ∎

Observe that every point $\bm{x}^{\prime}$ in the feasible set of P7 with objective value less or equal than zero is also a feasible point of P5 with objective value better than $\gamma$ . Thus, we can solve P5 sequentially by solving P7 with a BB method. Each time this BB algorithm finds a feasible point $\bm{x}^{\prime}$ with objective value less or equal than zero (in P7), the current best value $\gamma$ is updated with the objective value of $\bm{x}^{\prime}$ in the original problem P5 (plus the tolerance $\eta$ ). Then, the BB solver continues solving P7 with updated feasible set until it either finds a new point to update $\gamma$ or establishes that no solution to P7 with objective value less or equal than $-\varepsilon$ exists. By virtue of Proposition 1, the last feasible point $\bm{x}^{\prime}$ that was used to update $\gamma$ is an $(\varepsilon,\eta)$ -optimal solution of P5. This observation will be formalized in Proposition 4.

In the next section, we return our attention to P2 and use Algorithms 1 and 1 to solve it globally.

IV Robust Global Resource Allocation

We now apply the theory developed in the previous section to the solution of P2. Recall that a core idea is to exchange objective and constraints to obtain an optimization problem that is considerably easier to solve with a BB procedure than the original problem. This is mainly due to nice structural properties of the dual feasible set which facilitates an easy implementation of the feasibility checks required in BB methods and has no isolated points. These points are hard to compute and can lead to numerical instabilities.

Interchanging objective and constraints in P2 leads to

[TABLE]

or, equivalently,

[TABLE]

since $f^{-}(\bm{x},\bm{\xi})>0$ by assumption. The feasible set $\mathcal{D}$ of P8 is nice as long as the conditions in Case A or B is satisfied. This is formally established in the proposition below.

Proposition 2

The feasible set of P8 does not contain any isolated points if the conditions in Case A or B are satisfied.

Proof:

Please refer to the appendix. ∎

We now design a BB procedure to solve P8. Together with the SIT scheme in Algorithm 1 this will result in a method to solve P2 and is stated in Algorithm 2. The core idea of BB is to relax the feasible set and subsequently partition it such that lower bounds on the objective value can be determined easily. In our case, a rectangular subdivision procedure is a reasonable choice due to the required existence of a common maximizer of $\{g_{i}^{-}(\bm{x})\}$ over every box in the domain of P8.

Since the BB procedure is supposed to only operate on the global variables $\bm{x}$ , it partitions the $\bm{x}$ -dimensions of $\mathcal{D}$ successively into boxes $\{\mathcal{M}_{i}\}$ . Specifically, in iteration $k$ the algorithm selects a box $\mathcal{M}^{k}=[\bm{p}^{k},\bm{q}^{k}]$ with the lowest bound and bisects it via $(\bm{v}^{k},j_{k})$ , i.e., $\mathcal{M}^{k}$ is replaced by

[TABLE]

For each box $\mathcal{M}_{i}$ , a lower bound $\beta(\mathcal{M}_{i})$ for P8 with additional constraint $\bm{x}\in\mathcal{M}_{i}$ is computed. We will discuss the computation of $\beta(\mathcal{M}_{i})$ in the next subsection. For now it suffices to assume that it generates two points $\bm{x}^{k}\in\mathcal{M}^{k}$ , $\bm{y}^{k}\in\mathcal{M}^{k}$ satisfying

[TABLE]

for some $\bm{\xi}^{k}$ and with

[TABLE]

Thus, we can employ an adaptive bisection that exhibits much faster convergence than the common exhaustive subdivision, and choose the bisection parameters in 5 as $\bm{v}^{k}=\frac{1}{2}(\bm{x}^{k}+\bm{y}^{k})$ and $j_{k}=\operatorname*{arg\,max}_{j}\left\lvert y_{j}^{k}-x_{j}^{k}\right\rvert$ . For a formal proof of the convergence of such an adaptive BB procedure, we refer the interested reader to [3, Prop. 6.2].

IV-A Bounding

We now discuss the computation of lower bounds for P8 that satisfy 6. First, we establish a lower bound on the objective of P8.

Proposition 3

Let $\bar{\bm{x}}^{\ast}_{\mathcal{M}}$ be a common maximizer of $\{g_{i}^{-}(\bm{x})\}_{i=1,\dots,m}$ over the box $\mathcal{M}$ . Then, P8’s objective is lower bounded over $\mathcal{M}$ by

[TABLE]

This bound is tight at $\bar{\bm{x}}^{\ast}_{\mathcal{M}}$ .

Proof:

Please refer to the appendix. ∎

With Proposition 3 we can determine the lower bound $\beta(\mathcal{M}_{i})$ as the optimal value of

[TABLE]

This is a convex optimization problem if the conditions in Case A or B are satisfied and can be solved in polynomial time under very mild assumptions using standard tools [31]. For Case A, the objective of P9 is a convex function because $g_{i}^{-}(\bar{\bm{x}}^{\ast}_{\mathcal{M}_{i}})$ is constant and the feasible set is convex (cf. Proposition 2). Thus, P9 is convex given Case A. For Case B, P9 can be written as

[TABLE]

Let $\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}}$ be the common minimizer of $\gamma f^{-}_{x}(\bm{x})-f^{+}_{x}(\bm{x})$ , $g_{1,x}^{+}(\bm{x})$ , …, $g_{m,x}^{+}(\bm{x})$ over $\mathcal{X}\cap\mathcal{M}_{i}$ . Then P10 is equivalent to

[TABLE]

It is easy to see that $\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}}$ is the optimal solution of P10 since it jointly minimizes the objective and the first constraint. Problem P11 is convex due to the assumptions on $\Xi$ and the remaining functions of $\bm{\xi}$ made in Case B.

Finally, for each $\mathcal{M}^{k}$ , we identify the variables from 6 as $\bm{y}^{k}=\bar{\bm{x}}^{\ast}_{\mathcal{M}^{k}}$ and $(\bm{x}^{k},\bm{\xi}^{k})$ as the optimal solution of P9 or $\bm{x}^{k}=\underaccent{\bar}{\vec x}^{\ast}_{\mathcal{M}_{i}}$ if P11 is solved.

IV-B The SIT Algorithm

The original SIT algorithm from [3, Sect. 7.5.2] does not distinguish between global and non-global variables. We extent it such that branching is only performed over the global variables and state-of-the-art commercially available solvers can be used for the non-global variables. This preserves the computational complexity in the non-global variables, and increases computational performance and numerical accuracy compared to self-crafted algorithms due to the high maturity of these industry-grade solvers. We also extend it to fractional objectives which removes the necessity for Dinkelbach’s algorithm.

The BB procedure from the previous section solves P8, but that is not exactly what is required by the SIT scheme. Instead, Algorithm 1 requires the implementation of

Step 1

Check if P2 has a nonisolated feasible solution $\bm{x}$ satisfying $f(\bm{x})\geq\gamma$ ; otherwise, establish that no such $\varepsilon$ -essential feasible $\bm{x}$ exists and go to Step 3.

This is accomplished by a modified version of the adaptive BB algorithm from the previous section. Consider the following proposition which is adapted from [3, Prop. 7.14] and leverages the simple observation in Proposition 1.

Proposition 4

Let $\varepsilon>0$ be given. Either $g(\bm{x}^{k},\bm{\xi}^{\ast})<0$ for some $k$ and $\bm{\xi}^{\ast}$ or $\beta(\mathcal{M}^{k})>-\varepsilon$ for some $k$ . In the former case, $(\bm{x}^{k},\bm{\xi}^{\ast})$ is a nonisolated feasible solution of P5 satisfying $\frac{f^{+}(\bm{x}^{k},\bm{\xi}^{\ast})}{f^{-}(\bm{x}^{k},\bm{\xi}^{\ast})}\geq\gamma$ . In the latter case, no $\varepsilon$ -essential feasible solution $(\bm{x},\bm{\xi})$ of P5 exists such that $\frac{f^{+}(\bm{x},\bm{\xi})}{f^{-}(\bm{x},\bm{\xi})}\geq\gamma$ .

Proof:

Straightforward adaption of [3, Prop. 7.14]. ∎

Thus, an adaptive BB algorithm for solving P7 with deletion criterion $\beta(M)>-\varepsilon$ and stopping criterion $\operatorname{\text{min}}_{\bm{\xi}\in\mathcal{D}_{\bm{x}^{k}}}g(\bm{x}^{k},\bm{\xi})<0$ implements Step 1 in Algorithm 1: In the first case of Proposition 4 the incumbent feasible solution can be improved, in the latter case, if $\gamma=\frac{f^{+}(\bar{\bm{x}},\bar{\bm{\xi}})}{f^{-}(\bar{\bm{x}},\bar{\bm{\xi}})}+\eta$ for a given $\eta>0$ and a nonisolated feasible solution $(\bar{\bm{x}},\bar{\bm{\xi}})$ , the incumbent $(\bar{\bm{x}},\bar{\bm{\xi}})$ is an essential $(\varepsilon,\eta)$ -optimal solution of P5.

With the observation in Proposition 4, we can formally state the complete procedure for solving P2 with global optimality in Algorithm 2. It is initialized in Step ‣ 2 where an initial box $\mathcal{M}^{0}=[\bm{p}^{0},\bm{q}^{0}]$ is required that contains the $\bm{x}$ -dimensions of $\mathcal{C}$ , i.e.,

[TABLE]

The set $\mathscr{P}_{k}$ contains new boxes to be examined in Step 1, $\gamma$ holds the current best value adjusted by the tolerance $\eta$ , and $\mathscr{R}$ holds all boxes that are not yet eliminated. In Step 1 the bound is computed for each box in $\mathscr{P}_{k}$ . If it is less than $-\varepsilon$ , the box may contain a nonisolated feasible solution with objective value greater than $\gamma$ and is added to $\mathscr{R}$ . Then, in Step 3, the box with the smallest bound is taken out of $\mathscr{R}$ . If the point $\bm{x}^{k}$ attaining the bound is feasible in the original problem P2, it is a nonisolated feasible point and needs to be examined further: if the objective value for $\bm{x}^{k}$ is greater than the current best value $\gamma-\eta$ , $\bm{x}^{k}$ is the new current best solution and $\gamma$ is updated accordingly in Step 4. Irrespective of $\bm{x}^{k}$ ’s feasibility, the box selected in Step 3 is bisected via $(j_{k},\bm{v}^{k})$ adaptively. These new boxes are then passed to Step 1 and the algorithm repeats until $\mathscr{R}$ holds no more boxes which is checked in Step 2. Convergence of the algorithm is stated formally in the theorem below.

Theorem 1

Algorithm 2* converges in finitely many steps to the $(\varepsilon,\eta)$ -optimal solution of P2 or establishes that no such solution exists.*

Proof:

Please refer to the appendix. ∎

Observe that the BB procedure is directly incorporated in Algorithm 1. Thus, except for the first cycle, the BB procedure is started from the boxes already in $\mathscr{R}$ instead of starting from scratch [13, Remark 2]. Also observe that $\mathscr{R}$ is a priority queue [57, Sect. 5.2.3], and that P12 is a convex optimization problem. Due to the assumption of $\bm{\xi}$ being non-global variables, problem P13 is efficiently solvable. For example, if $f^{+}(\bm{x},\bm{\xi})$ is non-negative and concave in $\bm{\xi}$ , and $f^{-}(\bm{x},\bm{\xi})$ is convex in $\bm{\xi}$ it is solvable with Dinkelbach’s algorithm. Furthermore, if $f^{-}(\bm{x},\bm{\xi})$ is affine, $f^{+}(\bm{x},\bm{\xi})$ may also be negative. Please refer to [58, 50, 51] for more detailed treatments of fractional programming.

Remark 2

The original algorithm in [3] contains an additional reduction step before computing the bound of a box. This step is optional and should only be included if it speeds up the algorithm. A direct application of the approach outlined in [3] requires several convex optimization problems to be solved for each reduction. In our numerical experiments this slowed down the algorithm and is thus omitted. We leave the design of an efficient reduction procedure open for future work.

IV-C Extension: Pointwise Minimum

Algorithm 2 is easily extended to the case were the objective is the pointwise minimum of several functions, i.e.,

[TABLE]

In that case, all fractions in the minimum need to be greater or equal than $\gamma$ and the dual problem P8 becomes

[TABLE]

Then, the branching and bounding procedures are easily adjusted to this extended dual problem LABEL:opt:dualRAmin.

V Numerical Evaluation

The application of any optimization framework to a specific optimization problem requires, in general, some transformation of the initial problem to bring it into a form suitable for the framework. In addition, modification of the problem often allows to reduce the computational complexity significantly.

In this section, we first present the system model of a specific interference network, namely the Gaussian MWRC , and two of its achievable rate region. We then formulate the resource allocation problems and discuss, based on Section II-A, the application of Algorithm 2. Subsequently, we employ Algorithm 2 to obtain throughput and GEE optimal resource allocations and compare the performance of Algorithm 2 to the state-of-the-art. We conclude this section by evaluating how the performance of Algorithm 2 scales with an increasing number of global and non-global variables in Section V-D.

V-A System Model

We consider a 3-user SISO Gaussian MWRC [59, 60] with amplify-and-forward (AF) relaying, multiple unicast transmissions and no direct user-to-user links. Users are indexed by $k$ , $k\in\mathcal{K}=\{1,2,3\}$ , and the relay is node 0. User $k$ transmits with power $P_{k}$ over the channel $h_{k}$ to the relay. The relay propagates the observed symbol back to the users with transmit power $P_{0}$ over the channels $g_{k}$ , $k\in\mathcal{K}$ . Each node observes the independent and identically distributed (i.i.d.) zero-mean circularly symmetric complex Gaussian noise with power $N_{k}$ , $k\in\mathcal{K}\cup\{0\}$ and is subject to an average power constraint $\bar{P}_{k}$ on $X_{k}$ , $k\in\mathcal{K}\cup\{0\}$ . The message exchange is defined by the functions $q:\mathcal{K}\mapsto\mathcal{K}$ and $l:\mathcal{K}\mapsto\mathcal{K}$ , where the receiver of node $k$ ’s message is $q(k)$ and the user not interested in it is $l(k)$ . Without loss of generality we assume $q(1)=l(3)=2$ , $q(2)=l(1)=3$ , and $q(3)=l(2)=1$ . Due to space constraints we refer the reader to [61, 4, 1] for a more detailed treatment of this system model.

The receiver uses SND which is the optimal decoder for interference networks when restricted to random codebooks with superposition coding and time sharing [7]. We consider two different codebook constructions: traditional single message encoding and Han-Kobayashi [62] inspired rate splitting.

V-A1 Single Message

The achievable rate region is given below where $S_{k}=\frac{P_{k}}{N_{0}}$ , $\bm{S}=(S_{k})_{k\in\mathcal{K}}$ , and $\bar{\bm{S}}=(\bar{S}_{k})_{k\in\mathcal{K}}$ . Observe that this region is strictly larger than previously published SND regions [61, 4] and includes treating interference as noise (IAN) as a special case. This is due to recent insights on SND decoders [7].

Lemma 1

A rate triple $(R_{1},R_{2},R_{3})$ is achievable for the Gaussian MWRC with AF and SND if, for each $k\in\mathcal{K}$ ,

[TABLE]

or

[TABLE]

where $S_{k}\leq\bar{S}_{k}$ , $\delta_{k}(\bm{S})=1+\widetilde{g}_{q(k)}^{-1}\left(1+\sum_{i\in\mathcal{K}}\left\lvert h_{i}\right\rvert^{2}S_{i}\right)$ with $\widetilde{g}_{k}=\left\lvert g_{k}\right\rvert^{2}\frac{\bar{P}_{0}}{N_{k}}$ , and $\gamma_{k}(\bm{S})=\delta_{k}(\bm{S})+\left\lvert h_{l(k)}\right\rvert^{2}S_{l(k)}$ .

Proof:

The proof follows along the lines of [7, Sect. II-A.] and is omitted due to space constraints. ∎

Let $\mathcal{R}_{k,\text{IAN}}$ and $\mathcal{R}_{k,\text{SND}}$ be the regions defined by 9 and 10, respectively. Then, the rate region in Lemma 1 is

[TABLE]

Since $\inf_{\bm{x}\in\bigcup_{i}\mathcal{D}_{i}}f(\bm{x})=\min_{i}\inf_{\bm{x}\in\mathcal{D}_{i}}f(\bm{x})$ , we can split the resource allocation problem for Lemma 1 into eight individual optimization problems. Each is easily identified as an instance of LABEL:opt:exRA and solvable with Algorithm 2 using the initial box $\mathcal{M}_{0}=[\bm{0},\bar{\bm{S}}]$ .

V-A2 Rate Splitting

Each message is divided into a common message to be decoded by all receivers and a private part that is treated as additional noise by unconcerned receivers. These messages are then encoded by individual Gaussian codebooks with powers $P_{k}^{c}$ and $P_{k}^{p}$ and linearly superposed to be transmitted in a single codeword with power $P_{k}=P_{k}^{c}+P_{k}^{p}$ . The achievable rate region is given below in terms of the signal-to-noise ratios (SNRs) $S_{k}^{c}=\frac{P_{k}^{c}}{N_{0}}$ and $S_{k}^{p}=\frac{P_{k}^{p}}{N_{0}}$ . Further, we define $S_{k}=S_{k}^{c}+S_{k}^{p}$ , $\bar{S}_{k}=\frac{\bar{P}_{k}}{N_{0}}$ , $\bm{S}^{c}=(S_{k}^{c})_{k\in\mathcal{K}}$ , $\bm{S}^{p}=(S_{k}^{p})_{k\in\mathcal{K}}$ , $\bm{S}=(\bm{S}^{c},\bm{S}^{p})$ , and $\bar{\bm{S}}=(\bar{S}_{k})_{k\in\mathcal{K}}$ .

Lemma 2

A rate triple $(R_{1},R_{2},R_{3})$ is achievable for the Gaussian MWRC with AF relaying if, for all $k\in\mathcal{K}$ ,

[TABLE]

where $S^{c}_{k}+S^{p}_{k}\leq\bar{S}_{k}$ and

[TABLE]

with $\widetilde{g}_{k}=\left\lvert g_{k}\right\rvert^{2}\frac{\bar{P}_{0}}{N_{k}}$ .

Proof:

Please refer to the appendix. ∎

Global optimal resource allocation for this scenario is a straightforward extension of LABEL:opt:exRA:

[TABLE]

where $\kappa(i,j)$ maps from $(i,j)$ to the correct $k\in\mathcal{K}$ and $\bm{a}_{i},\bm{b}_{i,j}\geq 0$ are easily identified from Lemma 2. From the discussion in Section II-A it is apparent that LABEL:opt:HK, depending on the identification of $g^{+}_{i}(\bm{S},\bm{R})$ and $g^{-}_{i}(\bm{S})$ , satisfies both, Cases A and B, and is solvable with Algorithm 2 where the global and non-global variables are $\bm{S}$ and $\bm{R}$ , respectively. However, we can reduce the number of global variables in LABEL:opt:HK from six to four: When using the identification $g^{-}_{i}(\bm{S})\coloneqq-\sum_{j}\log(\gamma_{\kappa(i,j)}(\bm{S}))$ , the non-convexity due to $\bm{S}^{c}$ stems only from the sum $\sum_{k\in\mathcal{K}}\left\lvert h_{k}\right\rvert^{2}S^{c}_{k}$ in $\gamma_{k}(\bm{S})$ . Thus, if we replace this sum by an auxiliary variable $y$ , the variables $\bm{S}^{c}$ can be treated as non-global. The resulting problem is

[TABLE]

where $(\bm{S}^{p},y)$ are the global and $(\bm{S}^{c},\bm{R})$ are the non-global variables. Since the computational complexity grows exponentially in the number of global variables, solving LABEL:opt:HK2 instead of LABEL:opt:HK reduces the complexity significantly. The drawback of this approach is that Case B is no longer satisfied and the bounding problem becomes convex instead of linear. However, the performance gain outweighs these, comparatively small, performances losses significantly. The proposition below formally states the equivalence of LABEL:opt:HK and LABEL:opt:HK2.

Proposition 5

Problems LABEL:opt:HK and LABEL:opt:HK2 are equivalent in the sense that $v\eqref{opt:HK}=v\eqref{opt:HK2}$ if $f(\bm{S},\bm{R})$ is increasing in $\bm{R}$ for fixed $\bm{S}$ .

Proof:

Please refer to the appendix. ∎

Finally, note that all constraints in LABEL:opt:HK2 except the first are linear and, thus, $\mathcal{C}$ is a convex set. The initial box $\mathcal{M}_{0}$ required by Algorithm 2 is easily identified as $[\bm{0},\bar{\bm{S}}]\times[0,\sum_{k\in\mathcal{K}}\left\lvert h_{k}\right\rvert^{2}\bar{S}_{k}]$ .

V-B Throughput

We employ Algorithm 2 to compute the throughput optimal resource allocation for the rate regions in Lemmas 1 and 2, i.e., the objective is $f(\bm{S},\bm{R})=\bm{w}^{T}\bm{R}$ with $w_{i}=1$ for all $i$ and obviously fulfills Proposition 5. We assume equal maximum power constraints and noise power at the users and the relay, and no minimum rate constraint, i.e. $\underaccent{\bar}{\vec R}=\bm{0}$ . Channels are assumed reciprocal and chosen i.i.d. with circular symmetric complex Gaussian distribution, i.e., $h_{k}\sim\mathcal{CN}(0,1)$ and $g_{k}=h_{k}^{*}$ . Results are averaged over 1,000 channel realizations. The precision of the objective value is $\eta=0.01$ , $\varepsilon$ is chosen as $10^{-5}$ , and the algorithm is started with $\gamma=0$ .

Figure 2 displays the maximum throughput for SND as in Lemma 1 and rate splitting (RS) as in Lemma 2. Before [7], common wisdom was that neither IAN nor SND dominates the other rate-wise, with IAN generally better in noise limited scenarios and SND superior when interference is the limiting factor. This misconception is due to an longstanding oversight in the SND proof that was clarified in [7]. Figure 2 also includes results for “traditional” SND and IAN as obtained rate regions defined by 10 and 9, respectively. First, observe that, in accordance with conventional wisdom, neither “traditional” SND nor IAN dominates the other. Instead, “extended” SND clearly dominates the other two where the gain is solely due to allowing each receiver to either use IAN or “traditional” SND . The average gain of SND over the other two is approximately 11 % and 22 % at 10 dB, respectively, or 0.29 bpcu and 0.54 bpcu. Note that this gain is only achieved by allowing each receiver to choose between IAN and joint decoding which does not result in higher decoding complexity than “traditional” SND . The average gain observed for RS over single message SND is rather small, e.g., at 25 dB it is only 2 % (or 0.2 bpcu). However, depending on the channel realization we observed gains up to half a bit (or 4.6 %) at 20 dB. With spectrum being an increasingly scarce and expensive resource this occasional gain might very well justify the slightly higher coding complexity.

A popular approach to solving non-convex resource allocation problems with global optimality is to use monotonic optimization [25, 15]. As already pointed out in the introduction, the main challenge is that for SND the optimization is over the rates and powers (instead of just the powers for IAN ). Thus, every approach that treats the rates as global variables will suffer from very long computation time. Instead, in [4] the problem is solved by decomposing it into an inner linear and an outer monotonic program. We use this approach as the state-of-the-art reference for our performance comparison. The inner linear problem is solved with Gurobi 8 [63] and the monotonic program with the Polyblock algorithm777The Polyblock algorithm was implemented in C++ with similar techniques as in the implementation of Algorithm 2. Thus, performance differences should be mostly due to algorithmic differences. [23] with a tolerance of 0.01. Results match those obtained with our SIT algorithm for “traditional” SND and, thus, are not displayed in Fig. 2. Run times for both algorithms are reported in Table I. About 5.4 % of the computations with the Polyblock algorithm did not complete within one week. For these, a run time of one week was assumed for the computation of the mean in Table I. Thus, the reported run times for the Polyblock algorithm are an underestimate. Nevertheless, we observe, on average, roughly 10,000 $\times$ faster convergence for the SIT algorithm for low to medium SNRs . Interestingly, the median run time for high SNRs is lower for the Polyblock algorithm while the SIT is clearly faster on average.

V-C Energy Efficiency

The EE as defined in 2 is computed with the same parameters as before except for the precision which is chosen as $\eta=10^{-3}$ for Fig. 3 and $\eta=10^{-2}$ for Table II. Additionally, we assume the static circuit power $P_{c}=1\,\mathrm{W}$ , the power amplifier inefficiencies $\phi_{i}=4$ , and that the relay always transmits at maximum power. Results for SND , “traditional” SND , and IAN are displayed in Fig. 3. First, observe that the curves saturate starting from 30 dB as is common for EE maximization. In this saturation region, all three approaches achieve the same EE . However, for lower SNRs , IAN and, thus also SND , outperform “traditional” SND by 24 % on average at 10 dB. Of course, the EE performance depends quite a lot on the choice of real-world simulation parameters [64, 65], so further work is necessary to draw final conclusions in this regard.

The main point of this subsection, however, is the performance comparison of different computational approaches to the EE computation. Recall from Section II-A that there are two possible identifications of $g_{i}^{+}$ and $g_{i}^{-}$ that result in the problem either belonging to Case A or Case B. As discussed before, identification 1 (resp. Case A) leads to tighter bounds than 3 but requires the solution of a convex optimization problem in the bounding step. Instead, 3 (resp. Case B) leads to a linear bounding problem which is considerably easier to solve. Table II summarizes mean and median computations times for both approaches. Each bounding problem is solved with the fastest solver available: the convex problem that stems from Case A with Mosek [66] and the linear from Case B with Gurobi [63]. Observe, that Gurobi is, on average, thrice as fast as Mosek at 0 dB and Mosek is almost five times faster than Gurobi at 40 dB. However, these are only relative numbers. The average total computation time per channel realization is 8 s for Gurobi and 19 s for Mosek. So, despite the faster convergence speed of Case A (i.e., Mosek), Case B is much faster due to the lower computational complexity of the linear bounding problem.

The state-of-the-art approach to compute the EE is Dinkelbach’s Algorithm [48, 51]. This requires the global solution of a sequence of auxiliary problems with very similar properties to the throughput maximization problem in the previous section. Here, we solve this inner problem with the fastest method available: the SIT algorithm developed in this paper (cf. Table I). Hence, the differences in the run time are solely due to the use of Dinkelbach’s Algorithm. It can be observed from Table II that the inherent treatment of fractional objectives in our algorithm is always significantly faster (up to $800\times$ on average at 20 dB) than Dinkelbach’s Algorithm. Moreover, the obtained result is guaranteed to lie within an $\eta$ -region around the true essential optimal value.

V-D Benchmark of Algorithm 2

Above, we evaluated the performance of Algorithm 2 against the state-of-the-art on a real-world resource allocation problem. The goal of this subsection is to assess how Algorithm 2 scales with an increasing number of global and non-global variables. This is done by means of an idealized example. Namely, consider a $K$ -user Gaussian interference channel (GIC) with input-output relation

[TABLE]

where $h_{kj}\in\mathds{C}$ is the (effective) channel gain from transmitter $j$ to receiver $k$ , $x_{j}$ is the complex-valued channel input of transmitter $j$ that is subject to an average power constraint $\bar{P}_{j}$ , $y_{k}$ is the received symbol at receiver $k$ , and $z_{k}$ is circularly-symmetric complex Gaussian noise with zero mean and power $N_{k}$ . Consider just the $k$ -th receiver which is interested in decoding $x_{k}$ and assume it jointly decodes a subset $\mathcal{S}_{k}\subseteq\{1,2,\ldots,K\}$ of messages and treats the remaining $x_{j}$ , $j\in\mathcal{S}_{k}^{c}=\{1,2,\ldots,K\}\setminus\mathcal{S}_{k}$ , as noise. This is a MAC with capacity region

[TABLE]

where $P_{k}(\mathcal{S})=\sum_{j\in\mathcal{S}}\left\lvert h_{kj}\right\rvert^{2}P_{j}$ . The achievable rate region for the complete system then is

[TABLE]

For the benchmark, assume that receiver $i$ jointly decodes its own message and the interfering message of transmitter $i+1\mod K$ . All other messages are treated as noise. Moreover, assume that $h_{ij}=0$ for $i\neq j$ , $j\neq i+1$ , $j>\kappa$ , and some positive $\kappa$ . These assumptions allow to precisely control the complexities in the global and non-global variables. Thus, $\mathcal{S}_{k}=\{k,k+1\mod K\}$ for all $k$ and $\mathcal{S}_{k}^{c}=\{1,\dots,\kappa\}\setminus\mathcal{S}_{k}$ . Maximizing the throughput in this system requires the solution of

[TABLE]

First, observe that this is a non-convex optimization problem due to $P_{k}(\mathcal{S}_{k}^{c})$ in the denominators of the $\log$ -terms in 14. Further, observe that $\mathcal{S}_{k}^{c}\subseteq\{1,\dots,\kappa\}$ for all $k=1,\dots,K$ . Hence, we can identify the global variables as $(p_{k})_{k=1}^{\kappa}$ and the non-global variables as $(p_{k})_{\kappa+1}^{K}$ and $\bm{R}$ . Moreover, the powers $(p_{k})_{\kappa+1}^{K}$ only occur in the numerators of the $\log$ -terms in 14 and these $\log$ -terms are increasing functions in $(p_{k})_{\kappa+1}^{K}$ . Thus, $p_{k}=\bar{P}_{k}$ is the optimal solution for all $k=\kappa+1,\dots,K$ . Plugging-in this partial solutions in LABEL:opt:gic, we obtain a simplified version of LABEL:opt:gic with $\kappa$ global variables $(p_{k})_{k=1}^{\kappa}$ , $K$ linear variables $\bm{R}$ , and $3K$ inequality constraints. Thus, we can indepently control the number of global and non-global variables. Since the number of constraints grows linearly in the number of global variables and the bounding problem is a linear program, the computational complexity of solving LABEL:opt:gic should grow polynomially in the number of non-global variables.

This rationale is verified numerically in Fig. 4 where LABEL:opt:gic is solved for a fixed number of $\kappa$ global variables and an increasing number of $K$ total users in the system. We employ Algorithm 2 with Gurobi as inner solver and choose $\bar{P}_{k}=1$ and $N_{k}=0.01$ for all $k=1,\dots,K$ , i.e., a transmit SNR of 20 dB, and i.i.d. $h_{kj}\sim\mathcal{CN}(0,1)$ . Each data point is obtained by averaging over 1,000 independent channel realizations. Algorithm 2 was started with parameters $\eta=0.01$ , $\varepsilon=10^{-5}$ , and $\gamma=0$ . Since the abcissa is directly proportional to the number of non-global variables ( $K-\kappa$ ) and the number of inequality constraints ( $3K$ ), we expect to observe a polynomial growth of the run time. Indeed, Fig. 4 displays a roughly linear growth expect for small $K$ . This is because, for big $K$ , the optimal power allocation converges to setting all global variables to zero. Instead, for small $K$ these powers are generally non-zero and the global part of the optimization problem requires more time to converge to the optimal solution.

Finally, Fig. 5 verifies the exponential complexity in the number of global variables. With the same paramters as in Fig. 4, Fig. 5 displays the run time for a fixed number of $K=7$ users over the number of global variables $\kappa$ . Due to the logarithmical y-axis, exponential growth corresponds to (at least) a linear function in Fig. 5.

V-E Implementation Details

Algorithm 2 was implemented in C++ and compiled with GCC 7.3. Linear optimization problems are solved with Gurobi 8 [63] and non-linear problems with Mosek 8.1 [66] using the kindly provided academic licenses. The complete source code is available on GitHub [30].

All computations were done on TU Dresden’s Bull HPC-Cluster Taurus. Reported performance results were obtained on Intel Haswell nodes with Xeon E5-2680 v3 CPUs running at 2.50 GHz. We thank the Center for Information Services and High Performance Computing (ZIH) at TU Dresden for generous allocations of computer time.

VI Conclusions

We established $\varepsilon$ -essential feasibility as an important concept towards numerical stable global optimization algorithms and introduced the accompanying successive incumbent transcending (SIT)approach. Based on these concepts we built a novel global optimization framework tailored to resource allocation problems that preserves the computational complexity in the number of non-global variables, inherently supports fractional objectives, avoids numerical problems with non-robust feasible sets, and is four orders of magnitude faster than state-of-the-art algorithms. We applied the proposed scheme to a resource allocation problem for the Gaussian MWRC with AF relaying. Reproducible research and easy adoption of this algorithm is enabled by releasing the complete source code on GitHub [30].

-A Proof of Proposition 2

First consider Case A. The level set $\mathcal{F}=\{\gamma f^{-}(\bm{x},\bm{\xi})-f^{+}(\bm{x},\bm{\xi})\leq 0\}$ is closed and convex [67, Thms. 4.6 & 7.1]. Hence, $\mathcal{D}=\mathcal{C}\cap\mathcal{F}$ is also closed and convex. Every nonempty closed convex set is robust [67, Thm. 6.3]. Thus, $\mathcal{D}$ is robust or P8 is infeasible.

Now for Case B. The feasible set of P8 is $\mathcal{D}=\{\bm{x}\in\mathcal{X},\bm{\xi}\in\Xi:f_{\xi}(\bm{\xi})+f_{x}(\bm{x})\leq 0\}$ with $f_{\xi}(\bm{\xi})=\gamma f^{-}_{\xi}(\bm{\xi})-f^{+}_{\xi}(\bm{\xi})$ and $f_{x}(\bm{x})=\gamma f^{-}_{x}(\bm{x})-f^{+}_{x}(\bm{x})$ . By assumption, $f_{\xi}(\bm{\xi})$ is an l.s.c. convex and $f_{x}(\bm{x})$ an l.s.c. increasing (decreasing) function. Further, $\mathcal{X}$ is normal (conormal) within a box and $\Xi$ is convex. Observe that $\mathcal{D}$ is a convex set in $\bm{\xi}$ since for fixed $\bm{x}$ , $f_{x}(\bm{x})$ is a constant, $\tilde{f}_{\xi}(\bm{\xi})=f_{\xi}(\bm{\xi})+\mathrm{const.}$ is a convex function, and $\{\bm{\xi}:\tilde{f}_{\xi}(\bm{\xi})\leq 0\}$ is a closed convex set [67, Thms. 4.6 & 7.1]. By the same argument, $\tilde{f}_{x}(\bm{x})$ is an increasing (decreasing) function and $\{\bm{x}:\tilde{f}_{x}(\bm{x})\leq 0\}$ is a closed normal (conormal) set [3, Prop. 11.2]. Thus, $\mathcal{D}$ is normal (conormal) in a box in $\bm{x}$ . Neither closed convex nor closed (co-)normal sets have any isolated feasible points [14]. Since $\mathcal{D}$ is either convex or (co-)normal in each coordinate the proposition is proven. ∎

-B Proof of Proposition 3

For all $\bm{x}\in\mathcal{M}$ , $g_{i}^{-}(\bar{\bm{x}}^{\ast})\geq g_{i}^{-}(\bm{x})$ and hence also $g_{i}^{+}(\bm{x},\bm{\xi})-g_{i}^{-}(\bar{\bm{x}}^{\ast})\leq g_{i}^{+}(\bm{x},\bm{\xi})-g_{i}^{-}(\bm{x})$ with equality at $\bm{x}=\bar{\bm{x}}^{\ast}$ . It remains to show that for all real-valued functions $h_{1},h_{2},\dots,\underaccent{\bar}{h}_{1},\underaccent{\bar}{h}_{2},\dots$ satisfying $h_{i}\geq\underaccent{\bar}{h}_{i}$ and $h_{i}(\bm{y})=\underaccent{\bar}{h}_{i}(\bm{y})$ for all $i$ and some point $\bm{y}$ , $\max_{i}\{h_{i}\}\geq\max_{i}\{\underaccent{\bar}{h}_{i}\}$ and $\max_{i}\{h_{i}(\bm{y})\}=\max_{i}\{\underaccent{\bar}{h}_{i}(\bm{y})\}$ holds.

Consider the case with two functions and assume that $\max\{h_{1},h_{2}\}<\max\{\underaccent{\bar}{h}_{1},\underaccent{\bar}{h}_{2}\}$ . Since $h_{i}\geq\underaccent{\bar}{h}_{i}$ , this can only hold if $\max\{h_{1},h_{2}\}=h_{1}$ and $\max\{\underaccent{\bar}{h}_{1},\underaccent{\bar}{h}_{2}\}=\underaccent{\bar}{h}_{2}$ or vice versa. This implies $h_{1}\geq h_{2}\geq\underaccent{\bar}{h}_{2}$ which contradicts the assumption. The generalization to arbitrarily many functions follows by induction. Finally, if $\underaccent{\bar}{h}_{i}(\bm{y})=h_{i}(\bm{y})$ for all $i$ , then $\min_{i}\{\underaccent{\bar}{h}_{i}(\bm{y})\}=\min_{i}\{h_{i}(\bm{y})\}$ . ∎

-C Proof of Theorem 1

Convergence of the BB procedure is mostly established by Proposition 4 and [3, Prop. 6.2]. It remains to show that choosing $\beta(\mathcal{M}_{i})$ as the optimal value of P9 or P11 satisfies 6 and that continuing the BB procedure after updating $\gamma$ preserves convergence.

First, observe that neither P9 nor P11 relax the feasible set $\mathcal{D}$ of P8. Hence, the first condition in 6 is always satisfied. Further,

[TABLE]

Thus, $\beta(\mathcal{M}^{k})\rightarrow\operatorname{\text{min}}_{\bm{\xi}\in\mathcal{D}_{\bm{v}^{k}}}g(\bm{v}^{k},\bm{\xi})$ as $\left\lVert\bm{x}^{k}-\bm{v}^{k}\right\rVert=\left\lVert\bm{x}^{k}-\bar{\bm{x}}^{\ast}_{\mathcal{M}^{k}}\right\rVert\rightarrow 0$ and 6 is satisfied.

Next, let $\{\gamma_{k}\}$ be the sequence of updated gammas, and observe that this sequence is increasing, i.e., $\gamma_{k+1}\geq\gamma_{k}$ . Thus, the feasible sets of P9 form a decreasing sequence of sets

[TABLE]

with $\tilde{\mathcal{C}}=\{(\bm{x},\bm{\xi})\,|\,(\bm{x},\bm{\xi})\in\mathcal{C},\ \bm{x}\in\mathcal{M}\}$ . Therefore, the optimal value of the bound P9 is increasing with $k$

[TABLE]

where $\tilde{g}(\bm{x},\bm{\xi})=\max_{i=1,2,\dots,m}\left\{g_{i}^{+}(\bm{x},\bm{\xi})-g_{i}^{-}(\bar{\bm{x}}^{\ast}_{\mathcal{M}_{i}})\right\}$ . Hence, every box eliminated due to the deletion criterion $\beta(\mathcal{M})>-\varepsilon$ in a BB procedure with $\gamma_{k}$ , would also be eliminated in a procedure with $\gamma_{k+1}$ . It follows that the set holding the boxes $\mathscr{R}_{k}$ remains valid after updating $\gamma$ . In particular, no box is eliminated prematurely. Thus, restarting the BB procedure after updating $\gamma$ is not necessary. Convergence of Algorithm 2 is finite since the underlying BB procedure is finite and $\{\gamma_{k}\}$ is bounded. The same argument can be made for P11. ∎

-D Proof sketch for Lemma 2

Extend [61, Thm. 1] to Gaussian channels using the standard procedure in [68, Sect. 3.4.1]. Evaluate it with Gaussian inputs $U_{k}\sim\mathcal{CN}(0,P_{k}^{c})$ and $X_{k}=U_{k}+V_{k}$ with $V_{k}\sim\mathcal{CN}(0,P_{k}^{p})$ , and $\mathds{E}[X_{0}^{2}]=P_{0}$ to obtain the rate expressions above with $\widetilde{g}_{k}=\left\lvert g_{k}\right\rvert^{2}\frac{P_{0}}{N_{k}}$ . The achievable rates are increasing in $P_{0}$ . Thus, $P_{0}=\bar{P}_{0}$ is rate-optimal. ∎

-E Proof of Proposition 5

Let $\mathcal{F}$ be the feasible set of LABEL:opt:HK without the first constraint. Then, we have

[TABLE]

since relaxing a constraint does not decrease the optimal value. Conversely, $\gamma_{k}(\bm{S}^{p},y)$ is increasing in $y$ and, thus, the RHSs of the constraints are decreasing in $y$ . Since $f(\bm{S},\bm{R})$ and $\bm{a}_{i}^{T}\bm{R}$ are increasing in $\bm{R}$ , the RHS should be as large as possible. Thus, the optimal $y$ is as small as possible and

[TABLE]

This establishes Proposition 5. ∎

Bibliography68

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] B. Matthiesen and E. A. Jorswieck, “Optimal resource allocation for non-regenerative multiway relaying with rate splitting,” in Proc. IEEE 19th Int. Workshop Signal Process. Adv. Wireless Commun. (SPAWC) , Kalamata, Greece, Jun. 2018.
2[2] ——, “Global energy efficiency maximization in non-orthogonal interference networks,” in Proc. IEEE Int. Conf. Acoustics, Speech and Signal Processing (ICASSP) , Brighton, United Kingdom, May 2019.
3[3] H. Tuy, Convex Analysis and Global Optimization , ser. Springer Optimization and Its Applications. Springer, 2016.
4[4] B. Matthiesen and E. A. Jorswieck, “Weighted sum rate maximization for non-regenerative multi-way relay channels with multi-user decoding,” in Proc. IEEE 7th Int. Workshop Comput. Adv. Multi-Sensor Adaptive Process. (CAMSAP) , Curaçao, Dutch Antilles, Dec. 2017.
5[5] Z. Han and K. J. R. Liu, Resource Allocation for Wireless Networks: Basics, Techniques, and Applications . Cambridge University Press, 2008.
6[6] R. Horst and H. Tuy, Global Optimization: Deterministic Approaches , 3rd ed. Springer, 1996.
7[7] B. Bandemer, A. El Gamal, and Y.-H. Kim, “Optimal achievable rates for interference networks with random codes,” IEEE Trans. Inf. Theory , vol. 61, no. 12, pp. 6536–6549, Oct. 2015.
8[8] L. G. Khachiyan, “A polynomial algorithm in linear programming,” Doklady Academii Nauk SSSR , vol. 244, pp. 1093–1096, 1979.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Code & Models

Videos

Efficient Global Optimal Resource Allocation in Non-Orthogonal Interference Networks

Abstract

Index Terms:

I Introduction

Related work

I-A Major Contributions

I-B Notation & Preliminaries

I-C Outline

II Problem statement

Case A** (DC problems)**

Case B** (Separable problems)**

Remark 1

II-A Application Example: Resource Allocation in Interference Channels

Weighted Sum Rate

Rate Profile Approach

Global Energy Efficiency

Energy Efficiency Region

Proportional Fairness

III Robust Approach to Global Optimization

III-A The Issue with ε\varepsilonε-Approximate Solutions

Example 1** (Issues with ε\varepsilonε-approximate solutions)**

Example 2** (Isolated optimal solution)**

III-B ε\varepsilonε-Essential Feasibility

III-C Successive Incumbent Transcending Scheme

Proposition 1

Proof:

IV Robust Global Resource Allocation

Proposition 2

Proof:

IV-A Bounding

Proposition 3

Proof:

IV-B The SIT Algorithm

Proposition 4

Proof:

Theorem 1

Proof:

Remark 2

IV-C Extension: Pointwise Minimum

V Numerical Evaluation

V-A System Model

V-A1 Single Message

Lemma 1

Proof:

V-A2 Rate Splitting

Lemma 2

Proof:

Proposition 5

Proof:

V-B Throughput

V-C Energy Efficiency

V-D Benchmark of Algorithm 2

V-E Implementation Details

VI Conclusions

-A Proof of Proposition 2

-B Proof of Proposition 3

-C Proof of Theorem 1

-D Proof sketch for Lemma 2

-E Proof of Proposition 5

Case A (DC problems)

Case B (Separable problems)

III-A The Issue with $\varepsilon$ -Approximate Solutions

Example 1 (Issues with $\varepsilon$ -approximate solutions)

Example 2 (Isolated optimal solution)

III-B $\varepsilon$ -Essential Feasibility