arXiv CSopen access

Block coordinate descent for joint delay-energy optimization in multi-hop D2D networks

Block coordinate descent for joint delay-energy optimization in multi-hop D2D networks Kai-Xiang Hua,b,∗, Jacek Gondziob , Caixia Koua a School of Mathematical Sciences, Beijing University of Posts and Telecommunications, Beijing, 100876, China b School of Mathematics and Maxwell Institute for Mathematical Sciences, The University of Edinburgh, Edinburgh, EH9 3FD, Scotland, United Kingdom

arXiv:2606.08544v1 [math.OC] 7 Jun 2026

Abstract In multi-hop device-to-device (D2D) networks, the optimization of network-level metrics is particularly difficult due to the tight coupling between network-layer routing and physical-layer resource allocation. Departing from traditional average-performance metrics, this paper addresses the joint optimization of routing paths, transmission power, and bandwidth allocation. We formulate a generalized cost function to minimize the maximum transmission time (i.e., the bottleneck delay) alongside the total energy consumption. To tackle the resulting highly non-convex formulation, we propose a novel block coordinate descent (BCD) framework. At the network layer, we develop two adaptive routing algorithms: a matrix-free Frank-Wolfe (MF-FW) algorithm for fast execution in dense topologies, and a low-rank primal-dual interior-point method (LR-PDIPM) that bypasses dense matrix inversions via the Sherman-Morrison formula for high-precision solutions. At the physical layer, we design a parallel dual ascent algorithm leveraging a time-domain perspective transformation to solve the resource allocation subproblem to global optimality. The proposed BCD framework is proven to converge to an ε-neighborhood of a stationary point. Through comprehensive experiments, the proposed BCD framework establishes its superiority in achieving the optimal delay-energy trade-off. Specifically, the LR-PDIPM variant achieves a maximum 9.14× reduction in total energy consumption and up to an order of magnitude improvement in energy efficiency, while maintaining a bounded maximum delay gap (up to 3.78×) relative to the best baseline. Meanwhile, the warm-start MF-FW variant identifies near-optimal solutions in mere seconds, serving as a highly practical engineering approach. Keywords: Joint routing and resource allocation, Block coordinate descent, Frank-Wolfe algorithm, Interior-point method, Parallel dual ascent, Device-to-device networks 1. Introduction Device-to-device (D2D) communications have emerged as a cornerstone for future B5G/6G networks, facilitating proximity-based data exchange without traversing the base station [1, 2, 3, 4]. While single-hop D2D offers spectral efficiency gains, it is limited by short transmission ranges and susceptibility to signal blockage. Consequently, multi-hop multi-path architectures are indispensable for extending network coverage and ensuring robust connectivity in dense urban or industrial environments [5, 6, 7]. However, realizing the potential of multi-hop multi-path networks requires navigating a trade-off between end-to-end delay and energy efficiency. Time-sensitive applications (e.g., industrial automation) demand low transmission time [8, 9]. Crucially, in a multi-path transmission setting, the overall service completion time is dictated by the slowest routing path (i.e., the bottleneck flow). Therefore, simply minimizing the average delay is insufficient; instead, minimizing the min-max transmission time across all concurrent routing paths is paramount to prevent any single path from becoming a performance bottleneck. Yet, minimizing this delay typically necessitates aggressive power allocation and bandwidth usage, which conflicts di∗ Corresponding author

Email addresses: [email protected] (Kai-Xiang Hu), [email protected] (Jacek Gondzio), [email protected] (Caixia Kou)

rectly with the limited battery capacity of relay nodes [10, 11]. Purely delay-centric strategies risk draining the power of critical nodes prematurely, while purely energy-centric strategies result in unacceptable transmission time. Beyond the conflicting nature of these objectives, optimizing such a network presents two mathematical challenges. First, because our goal is to minimize the maximum transmission time among all routing paths, the resulting objective function is non-smooth and intractable. This leads to a “jagged” optimization landscape where the identity of the slowest path can suddenly jump from one route to another as resources are adjusted. This non-smoothness prevents gradient-based optimization tools from finding an exact descent direction, as they often get stuck at cusps where differentiability fails. Second, the decision variables for routing and resource allocation are strongly coupled. The transmission time on any given path depends not only on the traffic load but also on the transmission rates of its constituent links. Crucially, these rates are not fixed; they are variables that change dynamically based on the allocation of limited power and bandwidth. This creates a complex interdependency where routing strategy and resource allocation are intertwined: altering the resource allocation changes the transmission rate, which in turn necessitates a re-evaluation of the optimal routing paths. Existing literature inadequately addresses these challenges, leaving a gap between theoretical requirements and practical

solutions.

to achieve rapid execution. For high-precision scenarios, we develop a low-rank primal-dual interior-point method (LR-PDIPM). The LR-PDIPM leverages the ShermanMorrison (SM) rank-1 formula to bypass the direct inversion of the dense reduced Newton system, reducing the per-iteration computational complexity from O(|K|3 |V|3 ) to O(|K||V|3 ). Regarding the physical-layer allocation, we exploit a time-domain perspective transformation and Lagrangian dual decomposition to split the problem into linklevel subproblems. These subproblems are then solved to global optimality in parallel via coordinate descent integrated with inner 1D bisection search.

• While some studies employ alternating strategies to simplify the strongly coupled problem, they typically optimize the original problem or its subproblems using heuristic ideas [12, 13]. These alternating updates lack rigorous theoretical guarantees, often leading to solution oscillation or failure to converge to a stationary point. • Some existing studies focus on optimizing average delay or sum throughput [14, 15, 16]. These aggregate metrics mask the bottleneck flows. The resulting algorithms might achieve excellent average performance while leaving the bottleneck path congested, thereby violating the strict service completion time required by time-sensitive applications.

• Through rigorous derivation and comprehensive experiments, we prove that the proposed BCD framework converges to an ε-neighborhood of a stationary point. Numerically, we validate its superiority in achieving the optimal delay-energy trade-off. The BCD variant integrated with LR-PDIPM achieves a maximum 9.14× reduction in total energy consumption and up to an order of magnitude improvement in energy efficiency relative to the best baseline. Meanwhile, the warm-started MF-FW variant identifies near-optimal solutions in mere seconds, serving as a highly practical engineering approach.

• Standard gradient-based methods fail at cusps of the minmax function, prompting most researchers to resort to heuristic algorithms [17, 18]. While flexible, these heuristics lack the theoretical convergence guarantees required for high-reliability communications. To overcome these challenges, we adopt an overlay D2D architecture with orthogonal resource slicing. Unlike underlay D2D, which suffers from severe mutual interference within cellular networks, the overlay architecture allocates dedicated spectrum bandwidth to D2D users. This architecture eliminates cross-tier interference. However, the total bandwidth is constrained. To prevent severe resource overutilization, the system must jointly optimize its network-layer routing and physicallayer resources. While the overlay mode removes external interference, the joint optimization of routing and Shannoncapacity-based resource allocation remains highly coupled and non-convex. To tackle this, we formulate a generalized cost function. This function unifies the conflicting objectives of worst-case transmission time and network energy consumption. Furthermore, rather than relying on heuristics, we apply a smooth approximation to navigate the non-differentiable minmax objective. This transformation facilitates the application of the gradient-based optimization techniques and allows us to provide theoretical guarantees of convergence to a stationary point. The main contributions of this paper are summarized as follows:

The remainder of this paper is organized as follows. Section 2 reviews the related work. Section 3 presents the system model and formulates the JRRA problem. Section 4 details the proposed BCD framework, including the adaptive network-layer routing algorithms and the physical-layer resource allocation algorithm. Section 5 provides a rigorous theoretical analysis regarding algorithmic convergence and computational complexity. Section 6 presents extensive numerical results and performance evaluations. Finally, Section 7 concludes the paper. 2. Related works 2.1. JRRA in multi-hop D2D networks Given the inherent complexity of joint multi-dimensional resource management in D2D communications, early literature predominantly relied on simplified models. To ensure tractability, early works often assumed predetermined link scheduling or single-channel scenarios [19]. Consequently, many studies separated the tightly coupled variables, focusing strictly on single-dimensional power control and ignoring the joint allocation of spectrum resources. For instance, the work in [20] addressed the non-convex sum-rate maximization problem by modeling power allocation as a potential game, achieving the convergence to local maxima. Similarly, the work in [21] derived closed-form optimal power control strategies by partitioning circuit power consumption into distinct operational regions. While effective for single-dimensional resource management, these decoupled approaches fall short in modeling the strict interdependency among routing, spectrum, and power in multihop D2D networks. To capture these physical-layer dependencies, subsequent research transitioned toward joint resource allocation. In [22], the

• We formulate the joint routing and resource allocation (JRRA) problem in overlay D2D networks as a highly coupled, non-convex optimization problem. To overcome the prohibitive computational overhead, we propose a block coordinate descent (BCD) framework that jointly optimizes the network-layer routing and the physical-layer resource allocation. This framework decouples the intractable formulation into a sequence of solvable convex subproblems. • To address the routing subproblem, we develop an adaptive algorithmic framework. For dense topologies, we propose a matrix-free Frank-Wolfe (MF-FW) algorithm 2

authors proposed an iterative combinatorial auction algorithm to allocate both power and radio resources efficiently. Furthermore, to address practical battery limitations, works such as [23] and [24] incorporated Peukerts law to model nonlinear battery lifetime, demonstrating the necessity of strict distance constraints for D2D pairs to maintain energy efficiency. Notably, the study in [25] investigated energy-efficient allocation specifically in D2D overlay networks. They highlighted that allocating spectrum resources orthogonally completely eliminates cross-tier interference, significantly simplifying interference management. However, these studies were confined strictly to physical-layer resource allocation for predefined routing paths, leaving the combinatorial complexities of network-layer routing completely unexplored. In multi-hop settings, the coupling between physical-layer resource allocation and network-layer routing becomes the dominant performance bottleneck, prompting a shift toward cross-layer joint optimization. For example, the work in [26] investigated joint spectrum and power optimization tailored for multi-hop multi-path D2D video delivery. In the context of D2D-assisted decentralized learning, the work in [27] formulated a joint problem encompassing computing power, wireless resource allocation, and link selection to minimize a weighted sum of learning latency and energy consumption. Despite these cross-layer advancements, a gap remains in fully capturing the intricate coupling between physical-layer resource allocation and network-layer routing, particularly when attempting to bound the worst-case delay for bottleneck flows. Most of these joint frameworks either rely on aggregate metrics (e.g., weighted sums [27]) that mask bottleneck flows, or they assume a static network topology where routing paths are predetermined before resource allocation [26]. They fall short in accurately modeling the dynamic dependency of routing decisions on the continuously updating transmission rate, which fluctuate based on power and bandwidth allocation. These inaccurate allocation and routing strategies result in suboptimal network utilization, particularly failing to guarantee strict service completion times in latency-sensitive applications.

by a greedy resource allocation based on quality-to-rate ratios. The work in [29] proposed an alternating decomposition framework: the discrete routing variables are determined via a sequential greedy strategy, while the continuous resource allocation is optimized through an iterative potential game. In the context of D2D communications within 5G heterogeneous networks, the authors in [30] proposed an exemplary handover scheme based on fuzzy logic to decentralize softwaredefined radio control. When addressing the mixed binarycontinuous nature of joint routing and resource allocation, genetic algorithms are frequently employed. For instance, the work in [26] utilized genetic algorithms to solve the spectrum and power optimization problem, demonstrating high empirical video throughput and energy efficiency. However, despite their flexibility, meta-heuristics and fuzzy logic approaches lack theoretical convergence guarantees. They are particularly vulnerable to the non-smooth, “jagged” landscape of min-max objective functions, where they easily stagnate at suboptimal points without the rigorous theory to verify optimality or stationarity. Learning-based and data-driven approaches. To handle the inherent uncertainties and high computational costs of D2D networks, learning-based approaches have recently gained traction for managing the intractable nonlinear coupling and combinatorial complexity of joint routing and resource allocation. For instance, the authors in [31] formulated the joint channel and power allocation problem as an interactive learning task, specifically for scenarios where channel state information is unknown a priori. By employing a recency-based Q-learning method, they successfully tracked the dynamic environment, achieving an order-optimal policy with strict mathematical convergence. Shifting the focus toward algorithmic efficiency, the work in [32] leveraged a hierarchical extreme learning machine to bypass the heavy computational overhead of traditional iterative power control solvers. Their approach achieved nearoptimal transmission rates while significantly reducing the algorithmic execution time. For fully decentralized generation of D2D networks, the study in [33] utilized extended beliefdesire-intention agents and the AI-based algorithm to enable autonomous transmission mode selection, effectively balancing data rate maximization with power minimization under local computational constraints. To explicitly handle multi-hop topologies, the work in [34] implemented federated learning with in-network aggregation to reduce mesh network bottlenecks. Furthermore, the work in [35] leveraged graph neural networks to capture spatial dependencies, predicting link usage for multi-commodity flows and drastically reducing computation time. While these data-driven paradigms offer rapid execution, they suffer from inherent black-box uninterpretability and require extensive retraining when network topologies dynamically change. More crucially, statistical learning methods cannot provide the strict, deterministic worst-case delay guarantees mandated by critical overlay D2D networks.

2.2. Algorithmic paradigms: heuristics, learning, and mathematical programming Traditionally, the joint optimization of resource allocation and multi-hop routing is formulated using explicit discrete link selection alongside continuous variables for power and bandwidth. This modeling approach results in a highly complicated mixed-integer nonlinear programming (MINLP) problem. To solve this problem, existing literature generally relies on three primary algorithmic paradigms: meta-heuristics, learning-based approaches, and deterministic mathematical programming. Heuristic and meta-heuristic approaches. To bypass the mathematical intractability of the original mixed-integer formulation, early implementations often resorted to heuristic rules and meta-heuristic algorithms. For instance, the authors in [28] resorted to a low-complexity heuristic that utilizes BreadthFirst Search (BFS) for multi-path routing discovery, followed

Deterministic and mathematical programming. Given the limitations of heuristic and learning-based models in highly dynamic environments, it is imperative to rely on deterministic 3

mathematical programming. For pure physical-layer resource allocation, researchers have extensively leveraged convex and generalized convex optimization tools. In [36], the authors developed an approximate interior point method (IPM) to maximize sum throughput by simplifying the Hessian matrix inversion. To address fractional objectives such as energy efficiency, the authors in [25] transformed the mixed-integer nonlinear fractional programming problem into a parametric subtractive form, solving it via Dinkelbach’s method and augmented Lagrangian techniques. Similarly, the work in [20] modeled power allocation as a potential game to ensure convergence to a local maximum. While mathematically rigorous, these approaches primarily focus on continuous variables in single-hop or fixed-routing scenarios, ignoring the topological complexities of multi-hop paths. To explicitly handle the integer constraints inherent in routing and link scheduling, classical operations research techniques have been applied. For instance, the pioneering research in [37] established a mixed-integer linear programming framework for service overlay networks and solved it via the solver CPLEX. In [38], the authors exploited the minimum-cost maximum-flow graph structure to design an exact branch-andbound algorithm, providing a strict global optimality benchmark. The study in [39] formulated the unsplittable multicommodity flow problem using mixed-integer programming and solved it via a branch-and-price algorithm. To accelerate these exact but computationally demanding methods, hybrid frameworks have emerged. For instance, the work in [40] proposed a column generation framework that provides not only an optimality benchmark but also structural support for a deep reinforcement learning agent in dynamic mmWave networks. However, despite their mathematical exactness, these integerfocused methods suffer from severe scalability issues. More crucially, to maintain computational tractability, they often require linearizing or oversimplifying the nonlinear physicallayer capacity models. Therefore, these methods have limited capability to capture the exact, tight coupling between routing decisions and continuous resource allocation. To bridge the gap between discrete routing and continuous resource allocation, recent works have adopted decomposition techniques for the resulting MINLP models. In [12], The authors utilized a two-layer approach, combining Dinkelbach’s algorithm in the outer loop with alternating optimization in the inner loop. The authors in [41] addressed a highly nonconvex MINLP by decomposing it into sequential subproblems, systematically integrating reinforcement learning for dynamic environment exploration with majorization-minimization and Dinkelbach’s method for tractable continuous resource allocation. Furthermore, dealing with long-term system stability, the work in [42] combined Lyapunov-based distributed schemes with successive convex approximation (SCA) to maximize the time-averaged network sum-rate.

frameworks for JRRA problem, they are challenged by the specific requirement of worst-case performance in multi-hop networks. As highlighted earlier, minimizing the maximum path delay introduces a non-differentiable objective function. Crucially, this min-max formulation intricately couples the resource variables across intersecting paths, rendering standard decoupled optimization techniques ineffective. Existing decomposition methods often bypass this mathematical bottleneck by reverting to tractable time-averaged metrics [42, 43] or sum-rate objectives [20, 36], masking the network’s most critical bottleneck flows. Historically, the majority of literature has favored aggregate metrics due to their smooth and differentiable nature. For instance, numerous studies focus on maximizing the sum rate or total throughput of D2D pairs [20, 36, 44]. To handle larger topologies, the work in [45] maximized total backhaul throughput in mmWave mesh networks by transforming the resource allocation into a scalable mixed-integer linear programming. For delay-centric applications, the work in [46] minimized the total weighted transmission time in air-ground mesh networks via a three-stage decomposition. While some works, such as [47], attempt to introduce proportional fairness by maximizing the logarithmic sum of data rates, these aggregate or weighted-sum objectives mask the performance of bottleneck flows. Because the objective function landscape is smooth, standard gradient descent or linear programming techniques apply seamlessly, but the resulting resource allocation often leaves bottleneck nodes or links severely congested. Recognizing the necessity of strict worst-case guarantees, recent literature has shifted toward min-max (or max-min) formulations. In [48], the authors investigated energy-efficient power control by maximizing the minimum individual energy efficiency among D2D pairs. In the realm of mobile edge computing, the study in [49] formulated a joint hybrid beamforming and resource allocation strategy explicitly designed to minimize the maximum system delay. However, explicitly targeting the worst-case path introduces a severe mathematical barrier: the max(·) operator renders the objective function nondifferentiable at points where multiple path delays intersect. To handle this non-smoothness, standard approaches typically rely on epigraph transformations, introducing auxiliary variables to shift the non-differentiability from the objective function into a set of inequality constraints. Although theoretically sound, this reformulation severely inflates the problem’s dimensionality. More critically, these joint constraints couple the feasible domains of the optimization variables, thereby destroying the Cartesian product structure essential for BCD. Consequently, when dealing with the highly nonlinear flow and resource variables inherent in multi-path routing, this forced constraint coupling renders standard alternating optimization ineffective, often trapping the algorithm in deadlocks or necessitating computationally prohibitive inner loops. To resolve this non-smooth intractability without incurring a dimensional inflation, optimization theory offers rigorous smooth approximation techniques [50, 51]. Recently, the LogSum-Exp (LSE) approximation has gained significant traction in signal processing and communications to handle complex

2.3. Min-max objectives and the non-smooth coupling challenge While the aforementioned mathematical programming methods (e.g., alternating optimization and SCA) provide robust 4

Table 1: Summary of key notations

saddle point and min-max formulations [52]. For example, the work in [53] converted a non-convex, non-differentiable minmax fairness problem in an uplink edge computing system into a tractable form using LSE approximation, effectively managing the nonlinear coupling of design variables. Similarly, the study in [54] employed the LSE approximation combined with successive convex optimization for power control. By replacing the non-smooth discrete minimum objective with a concave LSE lower bound, they bypassed the need for iterative feasibility-check loops (e.g., bisection searches), thereby enabling gradient-based optimization. Despite these theoretical advancements, the application of smooth approximations in wireless networks remains largely confined to continuous physical-layer resource allocation or single-hop computation offloading. This min-max nondifferentiability is rarely addressed in multi-hop networks, where flow routing and resource allocation are tightly coupled. Motivated by this limitation, we introduce the LSE approximation into the tightly coupled domain of multi-path flow and resource allocation, thereby transforming the non-smooth worstcase objective into a mathematically tractable structure. Unlike conventional methods that struggle with the combinatorial explosion inherent to discrete network topologies, we adopt a node-arc continuous flow formulation. Within this purely continuous yet highly nonlinear formulation, our LSE approach fully unlocks the tractability of gradient-based alternating optimization for the JRRA problem in multi-hop D2D networks.

Symbol

Description

Sets and indices G(V, E) Directed network graph with node set V and link set E (i, j) ∈ E Directed link from node i to node j K Set of D2D communication commodities (traffic) k∈K Index of a specific commodity System parameters Mk Total traffic demand for commodity k (Mbits) B Total available bandwidth allocated to the overlay D2D network (MHz) Pmax Maximum transmission power budget at node i (dBm) i hi j Channel gain over link (i, j) N0 Noise power spectral density (dBm/MHz) α, (1 − α) Weighting coefficients for delay and energy in the generalized cost µ Smoothing parameter for the LSE worst-case delay approximation Decision variables xikj ∈ [0, 1] Flow proportion of commodity k routed over link (i, j) li j ≥ 0 Bandwidth allocated to link (i, j) pi j ≥ 0 Transmission power allocated to link (i, j) ti j ≥ 0 Transmission time over link (i, j) ri j ≥ 0 Transmission rate over link (i, j) mi j ≥ 0 Traffic load over link (i, j)

• Let li j denote the spectrum bandwidth allocated to link (i, j). The allocation is constrained by the available bandwidth budget B, which is exclusively partitioned among all active links.

3. System model and problem formulation The primary symbols and variables used in this paper are summarized in Table 1.

• Let pi j represent the transmission power assigned to link (i, j). Each node i is bounded by a maximum power budget Pmax , limiting the energy radiated by all its outgoing links. i

3.1. Network and traffic model Consider a multi-hop D2D network, modeled as a directed graph G = (V, E), where V denotes the set of nodes and E represents the set of directed wireless links (edges). The network serves a set of data commodities k. Each commodity k ∈ K is defined by a source-destination node pair, and a specific endto-end data demand M k (in Mbits). To capture the routing flows without enumerating all possible paths, we define continuous arc-based flow variables. Let xikj ∈ [0, 1] denote the flow proportion of commodity k routed through link (i, j). The topological structure of the network is governed by the node-arc incidence matrix A. These flow variables xikj must satisfy strict flow conservation constraints across all nodes, balancing the incoming and outgoing flows to match the source-destination requirements defined by vector bk .

In the considered D2D overlay network, the achievable transmission rate ri j over link (i, j) is governed by the allocated resources (l, p). Given the exclusive resource allocation inherent to the overlay architecture, multi-user interference is eliminated. Assuming an additive white Gaussian noise channel, the noise-limited link transmission rate is formulated as ! pi j hi j ri j = li j log2 1 + , (1) N0 li j where N0 denotes the noise power spectral density. The parameter hi j defines the channel gain of link (i, j), which follows the standard 3GPP propagation model hi j = 10−PL(di j )/10 . The distance-dependent path loss PL(di j ) is given by ! di j PL(di j ) = 128.1 + 37.6 log10 , (2) 1000

3.2. Physical layer and resource model The physical layer coordinates the allocation of spectrum bandwidth and transmission power. Under the overlay D2D network architecture, D2D communications are isolated from cellular traffic. Furthermore, the system employs exclusive subcarrier allocation to eliminate mutual co-channel interference among active D2D links. Within this interference-free environment, the physical resources are defined as follows:

where di j represents the Euclidean distance in meters between nodes i and j. 3.3. Problem formulation and arc-based approximation The primary objective is to strike an optimal trade-off between two competing performance metrics: the worst-case 5

transmission time and the total energy consumption. Following standard practices in multi-objective network optimization [27], we scalarize these metrics into a unified composite cost function. Specifically, we introduce a weighting parameter α ∈ [0, 1] to govern the trade-off between the two metrics.

the following continuous nonlinear programming model:

min

x,y,l,p

X (i, j)∈E

pi j

k k k∈K xi j M

ri j

Axk = bk ,

∀k ∈ K, ! pi j hi j ri j = li j log2 1 + , N0 li j X li j ≤ B,

k k k∈K xi j M

ri j

(4a) (4b)

∀(i, j) ∈ E,

(4c) (4d)

X

, pi j ≤ Pmax i

∀i ∈ V,

(4e)

∀k ∈ K, ∀(i, j) ∈ E,

(4f)

j∈Vout (i)

xikj ≥ 0,

(i, j)∈p

+ (1 − α)

s.t.

P

(i, j)∈E

   X M k  k   α max y p k∈K, p∈Pk ri j  P

min

X Mk α max xikj + (1 − α) pi j k∈K ri j (i, j)∈E (i, j)∈E

x,l,p

Ideally, the transmission time for commodity k must be evaluated over its explicit routed paths. Let ykp denote the flow proportion of commodity k routed through a specific path p belonging to the predefined path set Pk . The exact formulation of the objective function seeks to minimize the maximum path delay and the total energy simultaneously:

X

li j ≥ 0, pi j ≥ 0,

(3) .

∀(i, j) ∈ E.

(4g)

The objective (4a) minimizes a generalized cost function, balancing the maximum delay across all commodities and the total energy consumption. Constraints (4b) enforce the flow conservation, where A is the node-arc incidence matrix and bk denotes the normalized source-destination indicator vector for commodity k. The cross-layer coupling is captured by Shannon-capacity formula (4c), which maps the physical-layer resources (li j , pi j ) to the achievable network-layer rate ri j . Furthermore, (4d) and (4e) impose the resource budgets, bounding the aggregate bandwidth allocation by B and limiting the cumulative transmission power over the outgoing neighbour set Vout (i) of each node i to its budget Pmax . Finally, (4f) and (4g) i establish the non-negativity domain for all optimization variables. Although formulation (4) bypasses the exponential dimensionality of explicit path enumeration, direct optimization remains precluded by two mathematical challenges: nonsmoothness and non-convexity. First, the delay metric introduces the maxk∈K (·) operator, rendering the objective nondifferentiable at cusps where multiple commodities achieve identical transmission times. While subgradient methods can address non-smoothness, deploying them in highly non-convex domains typically induces severe iterative zig-zagging and fails to provide rigorous convergence guarantees to a stationary point. Second, the objective (4a) introduces a highly non-linear coupling between the routing variable xikj and physical-layer resources (li j , pi j ) through the transimission rate ri j . Although optimizing the physical-layer resources under a fixed routing strategy yields a convex subproblem (as proven in Section 4), the joint continuous domain remains highly non-convex. Consequently, achieving algorithmic tractability necessitates a smooth approximation of the max operator to ensure continuous differentiability, followed by a BCD framework to jointly optimize the cross-layer variables.

For clarity in presenting our numerical evaluations in Section 6, we introduce the shorthand notations T̂ = maxk∈K T̂ k to denote the network  bottleneck delay, where T̂ k =  P k max p∈Pk y p (i, j)∈p M k /ri j represents the maximum delay of commodity k. Accordingly, the total network energy consumpP P tion is denoted as Ê = (i, j)∈E pi j k∈K xikj M k /ri j . However, directly solving the formulation (3) is computationally intractable due to the exponential cardinality of the candidate path set Pk . Although modeled with continuous variables, this exponential dimensionality precludes direct optimization. Furthermore, to couple the path-based delay with the link-based energy consumption, the objective necessitates P a global mapping constraint, i.e., xikj = p∈Pk :(i, j)∈p ykp . Such explicit mapping breaks the block-separable structure inherent to node-arc formulations, thereby invalidating standard decomposition methods. From a numerical optimization perspective, applying an LSE smooth approximation to the max(·) operator over this exponentially large set Pk induces severe numerical overflow. To resolve these bottlenecks, we approximate the exact pathbased delay using arc-based variables. Specifically, we replace P the single-path delay ykp (i, j)∈p M k /ri j with the aggregate transP mission time (i, j)∈E xikj M k /ri j . Because the total time consumed across all utilized links is greater than or equal to the delay of any single end-to-end path, this summation establishes an upper bound for the original objective (3) and eliminates the dependency on the intractable path set Pk . Crucially, unlike traditional formulations that sum times across all commodities (which masks the bottleneck flow), our aggregate time is calculated independently for each individual commodity k. By preserving the outer maxk∈K operator, our objective identifies the bottleneck flow and minimizes the upper bound of its maximum transmission time.

4. BCD framework for joint optimization To tackle the intractable non-convexity and variable coupling in the joint optimization problem (4), we decouple this problem

Building upon this upper-bound approximation, we formulate the JRRA problem in multi-hop overlay D2D networks as 6

via a BCD framework, iteratively alternating between the network and physical layers. For the network-layer routing subproblem, we propose two customized algorithms to accommodate diverse scenario requirements: the MF-FW algorithm tailored for fast execution by leveraging the first-order gradient direction, and the LR-PDIPM that bypasses dense matrix inversions via SM formula for high-precision scenarios. For the physical-layer resource allocation, we employ a time-domain transformation to expose its hidden convexity. This subproblem is exactly solved via a parallel dual ascent (PDA) algorithm based on the Lagrangian decomposition. The rigorous theoretical analyses about these proposed algorithms are deferred to Section 5.

4.1.1. MF-FW algorithm for path routing For the requirement of fast execution, we propose a MF-FW algorithm to solve path routing subproblem (5), detailed in Algorithm 1 To prevent the zero-gradient stalling phenomenon associated with non-bottleneck commodities, we augment the LSE objective (4a) with a linear transmission time regularization term, scaled by ϵ > 0 (e.g., 10−6 ). The gradient-based link weights in (8) incorporate this regularization: wk,(n) = α(βk (x(n) ) + ϵ) ij

With the physical-layer resources (l, p) fixed, the overall joint optimization (4) reduces to a multi-commodity network flow (MCNF) model. As established previously, the maxk∈K (·) operator renders objective function (4a) non-differentiable. To restore gradient tractability without compromising the bottleneck-optimization objective, we employ LSE smoothing P
[50]. By approximating the max operator as µ1 ln k expµ(·) , we obtain the following differentiable and convex MCNF formulation:    X X α X min F(x) = ln  expµTk (x)  + cEi j xikj (5a) x µ k∈K

Axk = bk , xikj ≥ 0,

Gap(n) = −∇F(x(n) )T d(n) .

∀k ∈ K, ∀(i, j) ∈ E,

(9)

To suppress the zig-zagging phenomenon inherent to the standard FW algorithm near the optimum, we employ the directional curvature C (n) = (d(n) )T Hd(n) to determine the step size. Using the rank-one correction structure in (6), C (n) is calculated exactly in O(|K|) operations using the pre-aggregated P directional differences uk = (i, j)∈E dikj M k /ri j :

k∈K (i, j)∈E

∀k ∈ K,

(8)

Crucially, since this regularization introduces only a linear penalty, the second-order curvature of the objective remains unaffected. This preserves the exactness of the matrix-free Hessian computation in (10). Solving |K| independent shortest-path yields an auxiliary problems using these dynamic weights wk,(n) ij (n) (n) (n) flow y and a descent direction d = y − x(n) . The algorithm’s convergence is monitored via the standard Frank-Wolfe (FW) duality gap:

4.1. Path routing subproblem and smooth approximation

s.t.

Mk + cEi j . ri j

(5b) (5c)

  2  X X   C (n) = αµ  βk (uk )2 −  βk uk   .

where the parameter µ > 0 dictates the approximation tightness, P T k (x) = (i, j)∈E xikj M k /ri j is the aggregate transmission time for commodity k, and cEi j = (1 − α)M k pi j /ri j denotes the fixed link energy cost. To solve (5), we exploit the structural sparsity of its Hessian P matrix. Let βk (x) = expµTk / m expµTm define the LSE gradient component. Despite the exponential inter-flow coupling introduced by LSE, the exact Hessian H = ∇2 F(x) elegantly decomposes into a block-diagonal matrix minus a rank-one correction:   H = αµ W − vvT . (6)

k∈K

(10)

k∈K

This algorithm then determines a truncated Newton-based step size: ) ( Gap(n) (11) λ∗ = min 1, (n) . C Crucially, the strict convexity of the quadratic approximation for objective (5a) guarantees that this quadratic function of λ decreases monotonically on the interval [0, 1]. Thus, the truncation enforces λ∗ = 1 whenever Gap(n) /C (n) > 1. Due to the nature of convex combinations in standard FW, the flow x(n) may contain microscopic topological loops. In the output stage, we implement a standard BFS-based flow decomposition to extract loop-free flows.

Specifically, let the indices corresponding to the routing varik′ ables xikj and xmn be denoted by the tuples (k, i j) and (k′ , mn). The elements of the block-diagonal matrix W and the column vector v are given by:       βk M k /ri j M k /rmn , if k = k′ , W(k,i j),(k′ ,mn) =   0, otherwise (7)   k v(k,i j) = βk M /ri j .

4.1.2. LR-PDIPM for path routing For the requirement of high-precision solutions, we propose a tailored LR-PDIPM, detailed in Algorithm 2. The core computational bottleneck of conventional IPMs lies in solving the KKT Newton system [55, 56], which necessitates the dense matrix inversion of the modified Hessian Θ = H + X−1 S. Here, X = diag(x) and S = diag(s) denote the diagonal matrices constructed from the primal variable vector x and the corresponding dual slack vector s, respectively.

Here, the block-diagonal term W captures the independent intra-commodity routing curvature, while the rank-one correction vvT encapsulates the global cross-commodity coupling. Leveraging (6) and (7), we propose two customized algorithms. 7

√ where Λk = Xk −1 Sk is a diagonal matrix, and uk = αµβk dk with its elements defined by [dk ](i, j) = M k /ri j . Because Λk is diagonal, its inverse is trivial. Consequently, the block inverse D−1 k can be computed in linear O(|E|) operations via the SM formula:

Algorithm 1 MF-FW algorithm for path routing Input: Network graph G(V, E), demands {M k }, initial flow x(0) , physical resources (l, p), weight α, smoothing µ, tolerance ε, bounds Nmax . Initialize: Set n ← 0, Gap(0) ← ∞. 1: while n < Nmax do 2: Compute link weights wik,(n) via (8). j 3: For each k ∈ K, solve the shortest path problem using weights wk,(n) to construct auxiliary flow y(n) . ij 4: Determine the descent direction d(n) = y(n) − x(n) . 5: Compute the FW duality gap Gap(n) via (9). 6: if Gap(n) ≤ ε then break 7: Compute directional time variations uk and the directional curvature C (n) via (10). 8: Compute Newton-based step size λ∗ via (11). 9: x(n+1) ← x(n) + λ∗ d(n) , and n ← n + 1. 10: end while 11: Output: Extract optimal loop-free flows x∗ from x(n) via BFS.

−1 D−1 k = Λk −

∆x = D−1 vstep + ρ(v̂T vstep )v̂,

rd = −∇F(x) + AT y + s,

Consequently, the fraction-to-the-boundary rule confines the iterates within the wide neighborhood N−∞ (γ), formulated as o n N−∞ (γ) = (x, y, s) ∈ F 0 | xi si ≥ γµ(n) (20) gap , T where µ(n) gap = (x s)/|K||E| defines the current duality gap. Furthermore, the adaptive centering strategy is used to update the parameter σ(n) based on the step size λ. This design enforces strong centering when the search direction is blocked by the boundary (λ ≪ 1), while facilitating near-full Newton steps (λ ≈ 1) to achieve rapid convergence within N−∞ (γ).

(12)

rc = XS1 − τ1. (13)

Direct factorization of the system (12) is computationally prohibitive, demanding O(|K|3 |V|3 ) operations. Instead, we invert Θ via the SM formula: Θ−1 = D−1 + ρv̂v̂T , where v̂ = D−1 ṽ and ρ = (1 − ṽT v̂)−1 . This decouples the Schur matrix into M0 + ρzzT , where M0 = AD−1 AT and z = Av̂. To maintain this low-rank structure on the right-hand side, we construct the composite RHS vector g using the intermediate residual vrhs = rd − X −1 rc : g = −r p − AD−1 vrhs − ρ(v̂T vrhs )z.

4.2. Resource allocation subproblem and convex reformulation Given a fixed routing strategy x∗ obtained from the networklayer optimization (5), the joint formulation (4) reduces to the physical-layer resource allocation subproblem. Let mi j = P k k k∈K xi j M denote the aggregate traffic load assigned to link (i, j). The primal subproblem is formulated as: X pi j min αT + (1 − α) mi j (21a) l,p,T ri j (i, j)∈E

(14)

−1

s.t.

Crucially, because both D and A are block-diagonal, the matrix M0 decomposes into |K| independent sparse graph Laplacians. By solving M0 ĝ = g and M0 ẑ = z via parallel Cholesky factorizations, the exact dual step ∆y is recovered in merely O(|K||V|3 ) operations: ρ ∆y = ĝ − (zT ĝ)ẑ. (15) 1 + ρzT ẑ

X

xikj

(i, j)∈E

X

Mk ≤ T, ri j

∀k ∈ K,

li j ≤ B,

(21b) (21c)

(i, j)∈E

X

, pi j ≤ Pmax i

∀i ∈ V,

(21d)

pi j ≥ 0,

∀(i, j) ∈ E,

(21e)

j∈Vout (i)

li j ≥ 0,

Note that the block-diagonal matrix Dk embeds a local rankone correction structure: Dk = Λk + uk uTk ,

∆s = −X−1 rc − X−1 S∆x. (18)

Finally, to maintain strict positivity of (x, s) and prevent severe step size truncation, Algorithm 2 employs a wide neighborhood fraction-to-the-boundary rule. Let λ ∈ (0, 1] denote the primaldual step size and F 0 denote the strictly feasible primal-dual interior set, defined as n o F 0 = (x, y, s) | Ax = b, AT y + s = ∇ f (x), x > 0, s > 0 . (19)

where ∆y is the dual step and RHS vector g is defined in (14). Furthermore, the standard primal, dual, and central path KKT residuals are defined as: r p = Ax − b,

(17)

Since the IPM maintains X ≻ 0 and S ≻ 0, the scalar term uTk Λ−1 k uk is positive. Subsequently, the primal and dual slack steps are recovered via straightforward back-substitution. This computation decomposes into sparse matrix-vector multiplications and linear-time inner products, utilizing the auxiliary residual vstep = AT ∆y + vrhs :

Substituting (6), the modified Hessian inherits a rank-one correction: Θ = D − ṽṽT , where the base matrix D = αµW + √ X−1 S remains block-diagonal, and the scaled vector ṽ = αµv. Let A = I|K| ⊗ A denote the block-diagonal incidence matrix spanning all |K| commodities, distinguishing it from the singlecommodity incidence matrix A. Finding the Newton direction requires solving the Schur complement system: (AΘ−1 AT )∆y = g,

1 −1 T (Λ−1 k uk )(Λk uk ) . 1 + uk Λ−1 u k k T

where the transmission rate ri j = li j log2 (1 + pi j hi j /N0 li j ) couples the bandwidth and power variables. This equation renders both the transmission time 1/ri j in (21b) and the energy

(16) 8

Algorithm 2 LR-PDIPM for path routing Input: Network graph G(V, E), demands {M k }, initial flow x(0) > 0, physical resources (l, p), weight α, smoothing µ, tolerance ε, bound γ ∈ (0, 1), σmin , σmax ∈ (0, 1). Initialize: Set (y(0) , s(0) ) > 0, σ(0) ∈ (0, 1), n ← 0, µgap ← ∞. 1: while max(∥r p ∥∞ , ∥rd ∥∞ , µgap ) > ε do 2: Compute µgap = xT s/(|K||E|) and τ = σ(n) µgap . 3: Compute KKT residuals (r p , rd , rc ) via (13). 4: For each block k ∈ K, compute block inverse D−1 k via (17). , . . . , D−1 5: D−1 = blkdiag(D−1 K ). √1 6: Compute ṽ = αµv, v̂ = D−1 ṽ, scalar ρ = (1 − ṽT v̂)−1 and auxiliary vector z = Av̂ via (7). 7: Assemble the composite RHS g via (14). 8: For each block k ∈ K in parallel: T 9: Form M0,k = Ak D−1 k Ak . 10: Solve M0,k ĝk = gk and M0,k ẑk = zk via sparse Cholesky. 11: End For 12: Stitch vectors ĝk and ẑk to form global ĝ and ẑ. 13: Recover dual step ∆y via (15). 14: Recover primal and slack steps (∆x, ∆s) via (18). 15: Compute maximum step size λmax :
16: λmax = min 1, mini:∆xi <0 (−xi /∆xi ) , mini:∆si <0 (−si /∆si ) 17: Initialize λ ← 0.99λmax . 18: Compute µgap (λ) = (x + λ∆x)T (s + λ∆s)/(|K||E|) 19: while ∃i ∈ {1, . . . , |K||E|} s.t. (xi + λ∆xi )(si + λ∆si ) < γµgap (λ) do 20: Update λ ← c · λ, where c ∈ (0, 1) 21: end while    22: σ(n+1) = max σmin , min σmax , (1 − λ)2 . 23: (x, y, s) ← (x, y, s) + λ(∆x, ∆y, ∆s), and n ← n + 1. 24: end while 25: Output: Optimal primal-dual pair (x∗ , y∗ , s∗ ). consumption term mi j pi j /ri j highly non-convex with respect to (li j , pi j ), prohibiting direct global optimization. To overcome this intractability, we project the resource allocation space from the power domain into the time domain. Let ti j = mi j /ri j > 0 define the transmission time required to deliver the traffic load mi j over active link (i, j). Substituting p h ri j = mi j /ti j yields mi j /ti j = li j log2 (1 + Ni 0j liijj ). Inverting this relation isolates the transmission power as a function of allocated bandwidth and transmission time: pi j (li j , ti j ) =

 N0 li j  mi j /li j ti j 2 −1 . hi j

(22)

This time-domain transformation exposes a hidden convexity governed by the properties of perspective functions and convex composition. In the following content, we optimize the subproblem (21) over (li j , ti j ) and then recover the power variable pi j via (22).

vg(u/v) preserves joint convexity for all v > 0. Consequently, the bivariate function f (u, v) = uv(21/uv − 1) is strictly and jointly convex for all u, v > 0. Proof: See Appendix A.1. Leveraging Lemma 1, we linearize the delay constraints and convexify the power constraints, yielding the final convex formulation: X N0   min αT + (1 − α) li j ti j 2mi j /li j ti j − 1 (23a) l,t,T h (i, j)∈E i j X s.t. ρkij ti j ≤ T, ∀k ∈ K, (23b) (i, j)∈E

X

li j ≤ B,

(23c)

(i, j)∈E

X

pi j (li j , ti j ) ≤ Pmax , i

∀i ∈ V,

(23d)

j∈Vout (i)

li j ≥ 0,

ti j > 0,

∀(i, j) ∈ E,

(23e)

where the routing coefficient ρkij = xikj M k /mi j . To establish the theoretical foundation for global optimization, we formally prove the convex properties of (23). Theorem 1 (Convexity of the time-domain formulation). The transformed time-domain resource allocation formulation (23) is a continuous convex optimization problem. Proof sketch. The proof relies on establishing the joint convexity of the non-linear energy and power functions via affine domain mapping and the scalar composition theorem. Detailed derivations are provided in Appendix A.2. Despite the theoretical convexity established in Theorem 1, directly solving (23) using standard interior-point solvers (e.g., MOSEK) is practically infeasible. Modern solvers require models to be formulated exclusively using standard cone representations, primarily the exponential cone Kexp = {(x1 , x2 , x3 ) | x1 ≥ x2 exp(x3 /x2 ), x2 > 0} and the rotated quadratic cone Krquad = {(x1 , x2 , x3 ) | 2x1 x2 ≥ x32 }. To map the power and energy functions into these cones, the bandwidth-time product li j ti j must be decoupled. If we introduce an auxiliary variable Qi j to replace the product, the required constraint Qi j ≤ li j ti j defines a non-convex region. Alternatively, representing the product via a squared variable q2i j √ forms a valid rotated quadratic cone 2li j ti j ≥ ( 2qi j )2 . However, substituting this squared term q2i j into the energy function violates the requirement of the exponential cone Kexp , failing to obtain global optimum. Consequently, reformulating this model into solver-compatible cones inevitably introduces nonconvex constraints, making standard solvers inapplicable. To circumvent the solver’s incompatibility and exploit the inherent convexity and separability of the network resources, we develop a PDA algorithm via Lagrangian dual decomposition. 4.2.1. Parallel solution via dual decomposition By associating non-negative dual multipliers θk , ω, and φi with the delay, total bandwidth, and per-node power constraints

Lemma 1 (Convexity of perspective composite). If g(x) is a convex function, its perspective transformation P(u, v) = 9

P simplex k∈K θk = α. Specifically, the projection operator P∆ is evaluated by first computing the unprojected step θ̃k = P θk,(n) +δ(n) (i, j)∈Eact ρkij ti(n) j . The exact projection is then achieved by finding a unique root ν∗ such that: Xh i h i θ̃k − ν∗ = α =⇒ θk,(n+1) = θ̃k − ν∗ , ∀k ∈ K. (26)

respectively, the Lagrangian function associated with the convex formulation (23) is defined as:   X  X  k  L = α − θ  T − ωB − φi Pmax i i∈V k∈K     X  X   (1 − α)Ei j + φi pi j + ωli j +  θk ρkij  ti j , (24) + (i, j)∈Eact k∈K | {z }

k∈K

where Eact ⊆ E denotes the subset of active links with non-zero traffic load (mi j > 0), and Ei j = pi j ti j represents the energy consumption function of link (i, j). To prevent numerical singularities (e.g., division by zero), the resource allocation is exclusively executed over the active link set Eact , while inactive links are pruned. Applying the KKT stationarity condition with respect to the = 0, restricting the delay auxiliary variable T demands that ∂L P∂T multipliers to a scaled simplex: k∈K θk = α. Consequently, the global T term vanishes from the Lagrangian, decoupling the primal minimization into |E| independent, link-level 2D subproblems. Given fixed dual variables, each link minimizes its aggregate cost function Φi j (li j , ti j ) in parallel. Because Φi j is continuously differentiable and jointly convex in (li j , ti j ), its global minimum can be obtained via coordinate descent. Specifically, for a fixed transmission time ti j , the strict convexity of the link-level subproblem (Theorem 1) guarantees the strict monotonicity of its ∂Φ partial derivative ∂liijj . This property ensures a unique root for

 P  P    (i, j)∈Eact li(n)    j∈Vout (i) p(n) − B  − Pmax i j i j  , max    , ∆prim = max  max    i∈V  B Pi +

! ∥z(n+1) − z(n) ∥∞ , εnum > 0. ∆dual = max z={ω,φ,θ} ∥z(n) ∥∞ + ϵnum

The joint optimization problem (4) is non-convex due to the coupling between the routing variables x and the physical resources (l, p). To solve this problem while exploiting the biconvex structure of the decomposed two subproblems, we propose a BCD framework. This framework iteratively alternates between the network-layer routing and the physical-layer resource allocation. The strictly feasible initial point (l(0) , p(0) ) is generated via one iteration of PDA algorithm upon the K-shortest paths [57]. Let ñ denote the outer BCD iteration index. Each subsequent iteration n ≥ 0 proceeds through the following two phases: Phase 1: path routing optimization. Given the fixed resources (l(ñ) , p(ñ) ), the transmission rate ri(ñ) j and energy costs cE,(ñ) degenerate into constants. The joint problem collapses ij into a convex routing subproblem. We acquire the optimal flow x(ñ+1) by solving the routing subproblem (5). The proposed framework adaptively deploys either the MF-FW (Algorithm 1) for fast execution, or the LR-PDIPM (Algorithm 2) for high-precision solutions. Regardless of the deployed algorithm, the exact global optimum of this convex routing subproblem is guaranteed, thereby preserving the monotonic descent property of the overall BCD framework. Phase 2: resource allocation optimization. With the updated routing strategy x(ñ+1) fixed, the link traffic loads P m(ñ+1) = k∈K xik,(ñ+1) M k transition from decision variables to ij j fixed parameters. The optimization reduces to the convex timedomain formulation (23). By executing the PDA (Algorithm 3), we compute the optimal bandwidth and transmission time (l(ñ+1) , t(ñ+1) ) under current x(ñ+1) . Subsequently, the physical resource vector is recovered to (l(ñ+1) , p(ñ+1) ) via (22).

(25a) ∀i ∈ V,

+

θ



= P∆ θ

(n)



+ δ ∇θ L , (n)

(27b)

4.3. Overall BCD framework

∂Φ

(n+1)

+

(27a)

The iterative procedure terminates when max(∆prim , ∆dual ) ≤ ε. The detailed execution of PDA is presented in Algorithm 3.

the stationarity condition ∂liijj = 0, allowing the optimal bandwidth li j to be exactly located via 1D bisection search. A symmetric bisection procedure is subsequently applied to optimize ti j for a fixed li j . This coordinate descent method alternates iteratively until the link-level subproblem convergence is achieved. Once the primal variables (li j , ti j ) are optimized, the problem (23) progressively approaches the global optimum by updating the dual variables via projected gradient ascent. To mitigate the numerical ill-conditioning caused by heterogeneous parameter scales (e.g., MHz for bandwidth versus mW for power) and to suppress severe dual multiplier oscillations during early iterations, we employ normalized relative gradients and nonnegative projections. Let δ(n) denote the step size satisfying standard divergent sum conditions (e.g., δ(n) = δ0 /(n + 1)). The bandwidth, nodal power, and delay multipliers are updated as: P    (i, j)∈Eact li(n)  − B  j (n) (n)  , ω = ω + δ   B + P   max   j∈Vout (i) p(n)  − P  i ij (n) (n+1) (n)    , φi = φi + δ  max Pi

+

This scalar root-finding is a standard projection onto a simplex, which can be exactly solved in O(|K| log |K|) via classical sorting-based algorithms. To evaluate the convergence, we define the primal constraint violation ∆prim and the dual stability metric ∆dual based on the relative infinity norm:

Φi j (li j ,ti j )

(n+1)

+

(25b) (25c)

where [·]+ denotes the non-negative projection, ∇θk L = P k (n) (i, j)∈Eact ρi j ti j , and P∆ projects the vector θ onto the scaled 10

Algorithm 3 PDA algorithm for resource allocation Input: Initial (l(0) , t(0) ), active link set Eact , traffic loads mi j , QoS routing fractions ρkij , parameters B, Pmax , N0 , hi j , lmin , tmin , T max , i ε, εinner , Qmax , Nmax . Initialize: Initialize ω(0) = 0, φ(0) = 0, and θ(0) = α1/|K|. Set n ← 0, δ(0) ← δ0 , and ∆best ← ∞. 1: while max(∆dual , ∆prim ) > ε and n ≤ Nmax do 2: Parallel For (i, j) ∈ Eact : 3: Initialize q ← 0. (n) (0) (n) 4: Warm start: li(0) j ← li j and ti j ← ti j . 5: Repeat 6: q←q+1 ∂Φi j (q−1) via 1D 7: Update li(q) j by solving ∂li j = 0 given ti j bisection on [lmin , B]. ∂Φi j (q) 8: Update ti(q) j by solving ∂ti j = 0 given li j via 1D

Algorithm 4 BCD framework for joint optimization (4) Input: Network graph G(V, E), demands {M k }, physical parameters (h, N0 , B, Pmax ), objective weight α, smoothing µ, tolerances εout , max outer iterations Nout . Initialize: Initial solution (x(0) , l(0) , p(0) ). Set outer iteration ñ ← 0, global error ∆out ← ∞, pocket parameter J ∗ ← ∞. 1: while ∆out > εout and ñ < Nout do 2: Given (l(ñ) , p(ñ) ), compute transmission rate ri(ñ) j via (ñ) (1) and energy costs cE,(ñ) = (1 − α)M k p(ñ) ij i j /ri j . 3: Update routing variables x(ñ+1) via Algorithm 1 (MFFW) or Algorithm 2 (LR-PDIPM). P = k∈K xik,(ñ+1) Mk . 4: Compute link load m(ñ+1) ij j (ñ+1) (ñ+1) ,t ) via Al5: Update time-domain resources (l gorithm 3 (PDA). 6: Recover transmission power p(ñ+1) via (22). ij (ñ+1) via (4a). 7: Compute the objective J 8: if J (ñ+1) < J ∗ then J ∗ ← J (ñ+1) and (x∗ , l∗ , p∗ ) ← (x(ñ+1) , l(ñ+1) , p(ñ+1) ). 9: Update ∆out via (28). 10: Iteration Update: ñ ← ñ + 1. 11: end while 12: Output: Globally stationary routing and resource allocation (x∗ , l∗ , p∗ ).

QoS bisection on [tmin , T max ]. (q−1)
(q) | ≤ εinner or q ≥ 9: Until max |li j − li(q−1) |, |ti(q) j − ti j j Qmax . 10: Assign updated block solutions: li(n+1) ← li(q) j j and

ti(n+1) ← ti(q) j j . 11: End Parallel For 12: Update ω(n+1) , φ(n+1) , θ(n+1) via (25). 13: Compute ∆prim , ∆dual via (27). 14: If ∆prim ≤ ε and ∆dual < ∆best then 15: Update (l∗ , t∗ ) ← (l(n+1) , t(n+1) ) and ∆best ← ∆dual . 16: End If P 17: Compute T (n+1) = maxk∈K (i, j)∈Eact ρkij ti(n+1) . j √ (n+1) 18: Update δ = δ0 / n + 2. 19: n ← n + 1. 20: end while 21: Recover p∗ via (22). 22: Output: Optimal resources (l∗ , p∗ ).

T max (x) = maxk∈K T k (x) as the maximum aggregate transmission time across all commodities. Proposition 1 (LSE approximation bound). For any scaling parameter µ > 0, the LSE-smoothed bottleneck function  P T̃ max (x, µ) = µ1 ln k∈K expµTk (x) uniformly bounds the true maximum delay for any feasible routing flow x: T max (x) ≤ T̃ max (x, µ) ≤ T max (x) +

Termination. To prevent premature stopping of the heuristic, the BCD procedure evaluates the relative successive differences of both variable blocks. The algorithm terminates when: ! ∥z(ñ+1) − z(ñ) ∥∞ ∆out := max ≤ εout , (28) z={,l,p} ∥z(ñ) ∥∞ + εnum

ln |K| . µ

(29)

Proof. This bounding property is a fundamental result in convex analysis (see, e.g., Boyd and Vandenberghe [50, Sec. 3.1.5] and Nesterov’s smoothing techniques [51]). Thus, the detailed derivation is omitted here.

where εout is the outer tolerance, and εnum > 0. The execution of the BCD framework is summarized in Algorithm 4.

Remark 1 (Theoretical trade-off on smoothing parameter µ). Proposition 1 reveals a trade-off governed by the smoothing parameter µ. As µ → ∞, the approximation gap ln |K|/µ vanishes, recovering the exact non-smooth bottleneck objective. However, as established in Lemma 3, the Lipschitz constant of the Hessian scales quadratically with the smoothing parameter, i.e., LH = O(µ2 ). Consequently, an excessively large µ quadratically amplifies the local curvature variations, rendering the Newton system severely ill-conditioned. This forces the LR-PDIPM to take restricted step sizes (λ ≪ 1) to maintain interiority, thereby degrading the overall convergence rate. In practice, selecting a moderately large constant for µ strikes a necessary balance, guaranteeing numerical stability in the Newton directions while maintaining a sufficiently tight approximation of the true maximum delay.

5. Theoretical analysis and convergence guarantees In this section, we provide a rigorous theoretical analysis of the proposed BCD framework. We establish the convergence analysis for the customized subproblem algorithms. Specifically, MF-FW, LRPDIPM and PDA algorithms. Subsequently, we synthesize these components to prove the convergence of BCD algorithm. Finally, we conclude with an algorithmic complexity analysis. 5.1. Convergence of Network-Layer Routing Algorithms The network-layer optimization reduces to a convex routing subproblem (5) featuring the smoothed objective F(x). Define 11

5.1.1. Convergence of the MF-FW algorithm To establish the convergence rate of the MF-FW (Algorithm 1), we first prove that the LSE-smoothed objective exhibits Lipschitz continuous gradients (i.e., L-smoothness) with respect to the flow routing variables. Let X denote the feasible routing domain defined by the conn flow conservation and non-negativity o straints, i.e., X = x | Axk = bk , xk ≥ 0, ∀k ∈ K .

that the optimization error bound in Theorem 2 grows sublinearly with the network size, guaranteeing stable convergence even for dense D2D topologies. Ultimately, we bridge the gap between the solution obtained by MF-FW and the original non-smooth min-max formulation to quantify the total approximation error. Let J denote the nonsmooth objective (4a) and J ∗ be the minimum value of J.

Lemma 2 (Convexity and L-Smoothness). The smoothed objective function F(x) in (5) is convex and L-smooth over X. Since the individual commodity delay T k (x) is linear with respect to x, the Lipschitz constant L is bounded by the maximal gradient norm of the delay functions: 2 , L ≤ αµ max ∥∇T k (x)∥22 ≤ αµCmax k∈K

P

Theorem 3 (Total approximation error bound). The routing solution x(n) obtained at iteration n satisfies: J(x(n) ) − J ∗ ≤

2 D2 2αµCmax + +1 } | n{z Optimization Error

(30)

α ln |K| . µ | {z }

(32)

Smoothing Error

Proof. Let x∗F denote the optimum of formulation (5) with the smoothed objective F(x). By the approximation bound established in Proposition 1, we have J(x(n) ) ≤ F(x(n) ) and J ∗ ≥ F(x∗F ) − α ln |K|/µ. Subtracting the latter from the former yields:

where Cmax = maxk∈K (i, j)∈E M /ri j denotes the maximum cumulative path delay coefficient over the network. k

Proof. The detailed derivation is deferred to Appendix B.1.

α J(x(n) ) − J ∗ ≤ F(x(n) ) − F(x∗F ) + ln |K|. µ

Leveraging Lemma 2 and the standard Descent Lemma, we establish the sublinear convergence rate of the MF-FW algorithm. Let n denote the iteration index, and let D = maxx,y∈X ∥x − y∥2 denote the compact diameter of the flow conservation polytope.

Substituting the O(1/n) primal gap bound of F(x(n) ) − F(x∗F ) established in Theorem 2 directly completes the proof. This theorem characterizes the trade-off governed by the smoothing parameter µ: a larger µ suppresses the smoothing error but inflates the Lipschitz constant L, decelerating the MFFW convergence.

Theorem 2 (O(1/n) convergence rate of MF-FW). The primal optimality gap of the MF-FW algorithm at any iteration n satisfies: 2LD2 , ∀n ≥ 0, (31) F(x(n) ) − F(x∗ ) ≤ n+1 where x∗ = arg minx∈X F(x) is the exact global optimum of the smoothed subproblem (5), and L is the Lipschitz constant established in Lemma 2.

Remark 3 (Per-iteration computational complexity of MF-FW). In each inner MF-FW iteration, the gradient evaluation is computed in O(|K||E|) operations. The dominant computational overhead in each inner FW iteration is computing the descent direction d(n) , which requires executing Dijkstra’s algorithm for |K| commodities. Thus, assuming a Fibonacci heap implementation, the per-iteration complexity is bounded by O(|K|(|E| + |V| log |V|)), making the MF-FW algorithm exceptionally scalable for dense D2D networks.

Proof sketch. By invoking the standard descent lemma for Lsmooth functions, the objective reduction achieved by the truncated Newton step size λ∗ = min{1, Gap(n) /C (n) } is guaranteed to be no less than the reduction obtained via the classic diminishing step size γ(n) = 2/(n + 2). By applying the convexity of F(x), we establish a recursive contraction inequality for the primal optimality gap hn = F(x(n) ) − F(x∗ ). The global O(1/n) convergence rate is then obtained by resolving this recurrence via mathematical induction. Detailed derivations are deferred to Appendix B.2.

5.1.2. Convergence of the LR-PDIPM algorithm Unlike the first-order methods which typically exhibit sublinear convergence, the Newton-based LR-PDIPM achieves rapid, high-precision convergence for dense networks. Establishing the polynomial-time complexity of this algorithm requires bounding the nonlinearity introduced by the LSE smoothing and proving that the Newton steps decrease the duality gap. To guarantee the convergence and step-size stability of LRPDIPM, we must establish the second-order smoothness   of the smoothed objective F(x). Let Rmax = maxk,(i, j) M k /ri j denote the maximum path delay coefficient across the network.

Remark 2 (Polytope diameter bound). While the diameter D serves as a theoretical worst-case metric, it admits an upper bound governed by the network topology. Since the decision variables x represent routing proportions (xikj ∈ [0, 1]), any extreme point of the feasible polytope X corresponds to a set of single paths for the K commodities. Given that a loop-free path traverses at most |V|−1 edges, the squared L2 -norm of any path vector is strictly bounded by |V| − 1. Consequently, the maximum Euclidean distance between any two feasible flow assignp ments is bounded by D ≤ 2|K|(|V| − 1). This demonstrates

Lemma 3 (Lipschitz continuous hessian). The Hessian matrix ∇2 F(x) is LH -Lipschitz continuous over the feasible domain X. Specifically, there exists a constant LH > 0 such that: ∥∇2 F(x) − ∇2 F(y)∥2 ≤ LH ∥x − y∥2 , 12

∀x, y ∈ X,

(33)

where the Lipschitz constant scales quadratically with the smoothing parameter, bounded by LH = O(αµ2 |E|3/2 R3max ).

structure of the Schur complement without introducing numerical approximation errors, complexity is   √ the total arithmetic bounded by O |K||V|3 |K||E| ln(1/ε) , overcoming the prohibitive O(|K|3 |V|3 ) per-iteration cost of conventional dense Hessian inversion.

Proof sketch. Establishing the Lipschitz continuity of the Hessian requires bounding the spectral norm of the third-order directional derivative of the LSE function. The detailed derivation is provided in Appendix C.1.

5.2. Convergence of physical-layer resource allocation To establish a rigorous theoretical foundation, we analyze three different aspects of the proposed PDA algorithm: subproblem optimality, dual convergence, and computational complexity.

To preclude algorithmic stagnation, the backtracking line search yields a step size that is lower bounded by a positive constant. For any iterate within the wide neighborhood N−∞ (γ), the second-order cross terms of the Newton direction are quadratically bounded: ∥∆X∆s∥22 ≤ C(µgap )2 , where ∆X = diag(∆x). Exploiting this bound alongside the Lipschitz continuous Hessian (Lemma 3), we establish the following invariance and convergence guarantee.

5.2.1. Convergence of link-level subproblems The primal update involves a decoupled 2D subproblem Φi j (li j , ti j ) (defined in (24)) for each active link.

Lemma 4 (Wide-neighborhood invariance). The backtracking line search (Algorithm 2, Steps 17-21) is guaranteed to terminate with a step size λ ≥ λmin > 0, ensuring a monotonic gap reduction:
(n) µ(n+1) (34) gap ≤ 1 − λmin (1 − σmax ) µgap .

Lemma 5 (Optimality of link-level subproblem). The alternating coordinate descent procedure (Algorithm 3, Steps 1-9) converges to the unique global minimizer (li∗j , ti∗j ) of the Lagrangian subproblem Φi j . Proof sketch. According to (24), the subproblem objective Φi j is a positive linear combination of the energy consumption Ei j , transmission power pi j , and linear resource/time terms. By Lemma 1 and the proof of Theorem 1, Ei j is strictly convex and pi j is convex. Given weighting factor α ∈ (0, 1), the strict convexity is preserved in Φi j . Over the compact Cartesian product domain Fi j , the alternating minimization generates a descending sequence. The decoupled constraint structure ensures that every limit point is a stationary point, which is the unique global optimum. The detailed derivation is provided in Appendix D.1.

Proof sketch. The proof hinges on bounding the nonlinear Taylor residual of the Newton step. First, standard IPMs guarantee that the second-order cross terms are quadratically bounded by the duality gap. Second, because the Hessian is LH -Lipschitz continuous (Lemma 3), the deviation from the linear KKT prediction grows quadratically with the step size (i.e., O(λ2 )). Consequently, for sufficiently small λ, the linear descent of the duality gap (i.e., O(λ)) dominates the quadratic penalty, ensuring that the backtracking line search terminates at a positive step size λmin > 0 before violating the wide neighborhood conditions. The detailed derivation is deferred to Appendix C.2.

5.2.2. Convergence of PDA algorithm The convergence of the PDA algorithm hinges on the strict convexity of the Lagrangian function (24) and the boundedness of the resulting dual gradients. At the n-th iteration, let (n) T (n) T T g(n) ≜ [g(n) ω , (gφ ) , (gθ ) ] denote the gradient vector of (24). As formulated in the dual update rules (25), this vector evaluates the normalized residuals (i.e., relative violations) of the capacity and delay constraints at the current unique primal minimizer (l(n) , t(n) ). Note that the auxiliary delay variable T has P been eliminated via the dual simplex constraint k∈K θk = α.

Building upon Lemmas 3 and 4, we summarize the convergence and computational complexity of the proposed LRPDIPM. Theorem 4 (Computational complexity of LR-PDIPM). Under the strictly interior initialization (x(0) , s(0) ) > 0, the LR-PDIPM converges to an optimal solution of the convex smoothed routing subproblem (5). To achieve  √ an ε-optimal  solution, the algorithm requires at most O |K||E| ln(1/ε) iterations.

Lemma 6 (Bounded dual gradients). The dual gradients are uniformly bounded over all iterations. That is, there exists a finite constant G > 0 such that ∥g(n) ∥2 ≤ G, ∀n ≥ 0.

Proof sketch. The core logic relies on the interaction between the sequence feasibility and the monotonic gap reduction. First, Lemma 4 ensures that the line search generates a sequence bounded within the wide neighborhood, yielding a reduction rate of µ(n+1) ≤ (1 − δ)µ(n) . Second, guided by standard interiorpoint path-following theories, the single-step reduction factor δ is inversely proportional to the square root of the inequality constraint dimension |K||E|. Unrolling this reduction recursion to satisfy the √ ε-KKT tolerance yields the logarithmic iteration bound O( |K||E| ln(1/ε)). The detailed derivation is deferred to Appendix C.3.

Proof sketch. The gradient components correspond to the normalized constraint residuals of the bandwidth, power, and delay. Because the inner iterations confine the resource vari(n) ables (li(n) j , ti j ) within a compact Cartesian product domain QoS Fi j = [lmin , B] × [tmin , T max ], both the linear combinations and the continuous nonlinear power functions attain finite values. Dividing these finite values by their respective constraint limits ensures the relative violations remain finite. Consequently, the gradient vector is bounded. The detailed proof is provided in Appendix D.2.

Remark 4 (Computational complexity of LR-PDIPM). Since the nested SM rank-1 formula preserves the exact algebraic 13

To establish the overall convergence of the proposed PDA algorithm, we assume the original convex formulation (23) satisfies Slater’s condition, meaning there exists at least one strictly interior feasible point. This standard assumption guarantees strong duality.

solved via finite-step optimization rather than exact global oracles, the strict monotonic descent of objective J across outer BCD iterations cannot be guaranteed. Instead, the framework exhibits a bounded convergence behavior, necessitating the Pocket Mechanism to track the best feasible solution.

Theorem 5 (Global optimality of PDA). Suppose Slater’s condition holds. The dual sequence (ω(n) , φ(n) , θ(n) ) generated by the PDA algorithm converges to the optimal dual set. Consequently, the primal sequence (l(n) , t(n) ) converges to the optimum (l∗ , t∗ ) of the convex resource allocation subproblem (23).

Theorem 6 (Convergence of BCD framework). The sequence generated by the inexact BCD framework converges to an εneighborhood of a stationary point of the problem under the smoothed objective J̃. Proof sketch. The proof relies on the convergence properties of BCD framework under continuous differentiability and finite inner-loop truncation. First, the LSE approximation introduces a deterministic bounded gap sup |J̃−J| ≤ ln |K|/µ (Proposition 1). Second, within each BCD outer iteration, the routing and resource blocks are solved to their respective ε-optimalities rather than exact minima. According to standard inexact BCD theory for non-convex smooth functions [58], as long as the objective reduction achieved by each block update dominates the accumulated inner-loop truncation errors, the sequence guarantees monotonic descent. As the sequence approaches the stationary point and the gradients diminish, the truncation errors eventually lower-bound the descent, trapping the sequence within a bounded ε-neighborhood of the stationary point.

Proof sketch. In the dual domain, Lemma 6 guarantees that the dual gradients are uniformly bounded. According to standard gradient optimization theory, applying a square summable but not summable step size to bounded gradients ensures the asymptotic convergence of the dual objective to its optimum. In the primal domain, because the link-level Lagrangian subproblem Φi j is strictly convex (as established in Lemma 5), its primal minimizer is uniquely determined by the given dual variables. Therefore, the convergence of the dual sequence strictly forces the primal sequence to converge to the optimum without any duality gap. Remark 5 (Pocket mechanism and dual stability). Although the dual function is differentiable, the dual objective may still exhibit non-monotonic behavior due to the overshooting characteristic of predetermined step sizes. The incorporated Pocket Mechanism (Steps 14-16, Algorithm 3) tracks and preserves the feasible primal solution with the lowest objective value. This guarantees that the algorithm invariably terminates with the highest-quality feasible solution discovered across all iterations, neutralizing dual oscillations.

Remark 7 (Stabilization via Pocket Mechanism). While the sequence converges with respect to J̃, the oscillations of the nonsmooth objective J cannot be avoided due to the smoothing gap and finite inner-loop truncation. However, these are neutralized by the Pocket Mechanism (Step 8, Algorithm 4), which filters the non-monotonic solution ensuring that the output sequence of the objective is monotonically non-increasing.

Remark 6 (Computational complexity of PDA). The computational complexity of the PDA algorithm is evaluated based on its parallel execution time. In the primal phase, leveraging the decoupled structure, the 2D subproblems are solved in parallel across all active links. Thus, the time complexity for the primal update is O(log(1/εinner )), where εinner is the inner tolerance. In the dual phase, computing the gradient vector requires O(|Eact |) operations. Subsequently, updating the scalar multipliers ω and φ requires O(1) and O(|V|) operations respectively, while the simplex projection for θ takes O(|K| log |K|) operations. Given the network is connected (|Eact | ≥ |V| − 1), the per-iteration time complexity is bounded by O(log(1/εinner ) + |Eact | + |K| log |K|). Given the standard O(1/ε) iterations of the gradient method to achieve an ε-optimal   dual solution, the total time complexity  is bounded by O 1ε log(1/εinner ) + |Eact | + |K| log |K| . This linear-logarithmic scaling validates its practical efficiency for dense D2D deployments.

Building upon the established algorithmic analyses, the total computational complexity of the BCD framework (Algorithm 4) can be obtained. The complexity analysis begins with a one-time initialization phase (i.e., applying the PDA over the K-shortest paths). Assuming a Fibonacci heap implementation,  this procedure requires Cinit = O |K|Npath |V|(|E| + |V| log |V|) operations, where Npath denotes the predefined maximum number of initial routing paths for each commodity. Within the outer BCD loop, capped at a maximum of Nout outer iterations, each iteration alternates between the routing optimization and the resource allocation. Let Croute and Cpda denote their respective operational costs. As derived in Remark 3 and Remark 4, the routing complex
ity yields Croute = O Ifw |K|(|E| + |V| log |V|) when utilizing the MF-FW algorithm (assuming Ifw inner iterations), or   √ Croute = O |K||V|3 |K||E| ln(1/ε) via the LR-PDIPM. Meanwhile, as established in Remark 6, the operational cost for the gradient-based resource allocation is bounded by Cpda =   O 1ε log(1/εinner ) + |Eact | + |K| log |K| . Synthesizing these components, the overall worst-case time complexity of the proposed BCD framework is bounded by:    O Cinit + Nout Croute + Cpda . (35)

5.3. Convergence analysis of overall BCD framework Finally, we establish the theoretical convergence of the proposed inexact BCD framework. Let J̃ denote the LSE smoothed objective of the non-smooth objective function J defined in (4a). Because the inner block updates are iteratively 14

Here, T̂ , Ê denote the maximum delay and the total energy consumption defined in (3) respectively, and T 0 , E0 are the initial values of these two objectives obtained via the baseline KSP-PDA. Furthermore, we introduce the following multidimensional metrics to assess spatial utilization, energy efficiency, and commodity fairness:

6. Numerical results and performance evaluation In this section, we conduct comprehensive numerical experiments to validate the effectiveness of the proposed BCD framework. We first introduce the simulation parameters and performance metrics, followed by an evaluation of the Pareto tradeoff to determine the optimal weighting coefficient α. Under this optimal α, we conduct a sequence of experiments, focusing on the resilience, scalability, and robustness. For these experiments, 50 test instances are randomly generated. The code and all topology datasets are publicly available in our repository.1

• Active links & multi-path count: the total number of utilized links (li j > 0) and the number of commodities employing more than one routing path. • Energy efficiency (EE): defined as the ratio of total delivered data to the total energy consumed across the network, measured in Mbits/Joule: P Mk EE = k∈K . (37) Ê

Network topology and traffic generation. We consider a R × R m2 simulated area where |V| user nodes are randomly distributed within this area [27, 29]. Two directional links are established between any node pair if their distance satisfies di, j ≤ Dmax . To induce realistic multi-hop routing, the |K| commodities (source-destination pairs) are generated under a minimum distance constraint d src,dst ≥ 0.6R. The traffic demands M k are generated following a heterogeneous bimodal distribution: 80% of the streams are assigned light loads M k ∈ [0.1, 0.5] (e.g., periodic sensor telemetry), while the remaining 20% are assigned heavy loads M k ∈ [1.0, 2.0] (e.g., video streams or control commands). The key system and physical layer parameters are summarized in Table 2.

• Jain’s fairness index (JFI): evaluated based on the maximum transmission times T̂ k of all K commodities to quantify the severity of the bottleneck: P JFI =

|K|

k∈K T̂ k

2

P

2 k∈K T̂ k

.

(38)

Parameter

Value

System radius (R) Maximum D2D link distance (Dmax ) Number of nodes (|V|) Average node degree Number of D2D commodities (|K|) Commodity demand (M k ) Total available bandwidth (B) Max transmission power per node (Pmax ) i Noise power spectral density (N0 ) Channel path loss model (3GPP)

500 m 200 m 60 5∼7 20 0.1 ∼ 2.0 Mbits 100 MHz 23 dBm −174 dBm/Hz di j ) 128.1 + 37.6 log10 ( 1000

Total energy consumption (J)

Table 2: Simulation parameters

0.8

0.2

0.7 0.1

0.6 0.5

0.3 0.4 0.3 0.2 1

3

5

0.2

0.1

7

9

11

13

15

17

Maximum delay (s)

Figure 1: Pareto frontier illustrating the trade-off between maximum delay and total energy consumption.

We first investigate the Pareto optimal frontier of the proposed BCD-FW and BCD-IPM under various α values. As illustrated in Fig. 1, the knee point α = 0.4 achieves the optimal trade-off between maximum delay and total energy consumption. Any further reduction in either delay or energy consumption beyond this point results in a severe surge in the other metric. Consequently, we set α = 0.4 for all subsequent evaluations.

Performance metrics. We record the normalized values of original objective (3) alongside its individual components: the maximum transmission time and the total network energy consumption. The normalized value is computed as αT̂ /T 0 + (1 − α)Ê/E0 .

BCD-IPM

0.3

0.0

To evaluate the performance of our proposed BCD framework, we benchmark it against the state-of-the-art and conventional baselines. These evaluated algorithms are summarized in Table 3, categorized by their underlying path strategies, routing and resource allocation algorithms, optimization architectures, and initialization strategies. Notably, we introduce BCD-FWWARM and BCD-IPM-WARM, which leverage the heuristic solution of GG as a high-quality initial point to accelerate the solution progress.

BCD-FW

0.9

0.4

(36)

6.1. Case study To demonstrate the optimization mechanism from a microscopic perspective, we extract a case study (Topology 10 in

1 https://github.com/Gthu/JRRA-for-D2D-Networks

15

Table 3: Taxonomy of the evaluated algorithms

Algorithm

Path strategy

Routing solver

Resource allocation

Architecture

Initialization

BCD-FW BCD-FW-WARM BCD-IPM BCD-IPM-WARM GG [29] SP-SCA [46] SP-PDA KSP-PDA SP-SA

Multi-path Multi-path Multi-path Multi-path Single-path Single-path Single-path Multi-path Single-path

MF-FW MF-FW LR-PDIPM LR-PDIPM Greedy Dijkstra Dijkstra Dijkstra K-Shortest path Dijkstra

PDA PDA PDA PDA Potential Game SCA, Projected gradient PDA PDA Equal partition

Joint optimization Joint optimization Joint optimization Joint optimization Decoupled optimization Decoupled optimization Decoupled optimization Decoupled optimization Decoupled optimization

KSP-PDA solution GG solution KSP-PDA solution GG solution N/A N/A N/A N/A N/A

Coverage radius: 500 m

Coverage radius: 500 m

10 2

10 2

Power (W)

10 1

Power (W)

10 1

10 3

10 3

10 4

10 4

Shared with SP-SA BCD-IPM unique link Source nde Destination node

Shared with SP-SA GG unique link Source nde Destination node

(a) BCD-IPM: multi-path routing trajectories.

(b) GG: single-path routing trajectories.

Figure 2: Microscopic visualization of routing trajectories and power allocation across Topology 10. The thickness of each link is proportional to its allocated flow volume. Furthermore, the specific links constituting the multi-path routing generated by the BCD-IPM are displayed with a thermal colormap (transitioning from deep blue to dark red), to reflect the transmission power usage.

Section 6.4). In this instance comprising 60 nodes and 20 commodities, baselines fall short in achieving a good trade-off between the two objectives, whereas the BCD-IPM discovers a superior routing and resource allocation solution. To expose the underlying mechanisms driving this performance gap, Fig. 2 visualizes the precise routing trajectories generated by the BCDIPM and the GG heuristic. For clarity, the visualization isolates only the active nodes and links participating in data transmission. As depicted in Fig. 2, we establish a microscopic visualization of the routing trajectories mapped across Topology 10. To construct a routing strategy reference, physical links shared with the SP-SA shortest-path baseline are rendered in light grey, whereas activated topological links are delineated in dark charcoal. Furthermore, the thickness of each active link is scaled in proportion to its allocated data flow volume. To deconstruct the localized power allocation, we isolate three specific commodities ((10, 17), (45, 44), and (56, 23)) that trigger multi-path routing within the BCD-IPM algorithm. The routing paths servicing these targeted commodities are displayed with a thermal

colormap; the color continuously transitions from deep blue to dark red to reflect the escalating transmission power (in Watts) allocated in each link. As illustrated in Fig. 2a and 2b, both algorithms activate distinguished routing paths relative to the SP-SA baseline. The proposed BCD-IPM algorithm explores the spatial routing domain to achieve the optimal delay-energy trade-off. This algorithm dynamically activates multi-path offloading specifically when the energy-to-throughput translation rate on primary paths deteriorates. Conversely, the GG baseline aggressively selects single-path routes based on a greedy marginal-delay, evaluating the ratio between expected flows and transmission rates without considering the physical resource penalties. This algorithmic discrepancy is visually striking. In Fig. 2a, the thermal colormap reveals that the power allocation across the vast majority of BCD-IPM links is strictly bounded below 0.1 W, dominated by cool blue and light red links. In stark contrast, Fig. 2b is heavily saturated with dark red links, exposing the massive energy consumption of the GG heuristic. This severe power inflation is particularly severe along the backbone 16

Table 4: Microscopic resource allocation profile on the congested backbone trajectory (Commodity 17: Node 10 to 17)

Aggregated metrics (all commodities) Link

(10, 48) (48, 0) (0, 8) (8, 26) (26, 14)

Routed commodities

Bandwidth (MHz)

For Commodity 17 Power (W)

Delay (s)

GG

BCD-IPM

GG

BCD-IPM

GG

BCD-IPM

GG

BCD-IPM

5 5 7 9 12

2 2 4 6 9

1.909 1.909 3.687 5.887 10.113

0.546 0.466 1.386 1.430 3.005

0.050 0.067 0.046 0.033 0.033

0.002 0.001 0.004 0.002 0.005

0.031 0.023 0.017 0.010 0.007

0.155 0.132 0.059 0.050 0.026

route (10 → 48 → 0 → 8 → 26 → 14 → 59 → 42 → 54 → 24 → 21 → 45 → 17) for commodity 17. To deconstruct this phenomenon, Table 4 records the microscopic resource allocation profile along this specific backbone route. The empirical data exposes a severe resource overallocation within the GG baseline. For instance, across the link (48, 0), GG forces an aggregation of 5 distinct commodities. To push this massive load through the bottleneck and maintain a superficially low link delay of 0.023 s for Commodity 17, GG inflates this link’s transmission power and bandwidth to 0.067 W and 1.909 MHz, respectively. In stark contrast, BCD-IPM overcomes this resource over-allocation via load balancing. It traverses Commodity 17 through different paths in favor of underutilized topological regions (such as nodes 31 and 47), bypassing the heavily congested nodes (such as nodes 59, 42, and 54) . On this critical link (48, 0), BCD-IPM slashes the required transmission power to a mere 0.001 W by limiting the local aggregation to only 2 commodities. Scaling this microscopic observation to the network level, to service all 20 commodities, the GG baseline consumes all the 100 MHz of total network bandwidth and 1.679 W of accumulated transmission power. Conversely, the BCD-IPM algorithm utilizes merely 45.821 MHz of bandwidth and 0.136 W of total power to service these commodities. This constitutes an astonishing 12.3-fold reduction in total transmission power (an approximately 91.9% decrease) and a 54.2% saving in bandwidth consumption. While the BCD-IPM yields a higher maximum end-to-end delay (3.84 s versus the 1.78 s achieved by GG), it prevents the severe network-wide resource depletion.

with a high-quality initial solution derived from the GG heuristic. Conversely, while the GG heuristic exhibits relatively stable behavior under extreme stress (ξ ≥ 9), it operates on a suboptimal baseline. Under low-to-moderate conditions (ξ ≤ 7), GG consistently yields inferior objective values compared to the proposed BCD framework. As the stress factor increases (ξ ≥ 9), an intriguing phenomenon emerges. Although the warm-started algorithms stabilize the objective trajectory by leveraging the priors of the GG heuristic, they are occasionally susceptible to the local optimum trap. Specifically, at ξ = 11, both BCD-FW-WARM and BCD-IPM-WARM converge to an objective value of 0.573, identical to the GG baseline. This indicates that the GG initialization anchores these two algorithms within a local minimum. In contrast, the proposed BCD-IPM and the baseline SP-PDA discover superior solutions, yielding optimal objective values of 0.498 and 0.492, respectively. Despite these localized anomalies, the SP-PDA and KSP-PDA heuristics occupy the second performance tier in the normalized objective. This confirms that the optimal PDA resource allocation remains effective at mitigating congestion and compressing the energy consuption, even when constrained to single-path routing geometries. The SPSCA and SP-SA baselines consistently collapse into the lowest performance tier, exhibiting a massive optimality gap. Ultimately, this comprehensive resilience evaluation proves that the superiority of the proposed BCD framework stems from the synergy between the multi-path routing and optimal PDA resource allocation. Fig. 3b and 3c illustrate the individual objective performance for the evaluated algorithms. The proposed BCDFW and BCD-IPM dominate in total energy consumption. While these two algorithms exhibit the second performance tier regarding maximum delay, they consistently maintain a bounded gap between 1.55× and 3.78× relative to the GG baseline. To achieve its delay advantage, the GG heuristic consumes 3.41× to 9.14× more energy compared to the BCD-FW and BCD-IPM. By anchoring their initialization with the GG baseline, the warm-started BCD-FW-WARM and BCD-IPMWARM achieve highly competitive maximum delays. They consistently guarantee a worst-case trade-off (e.g. ξ = 7, 11) bounded by the GG heuristic. The SP-PDA and KSP-PDA consistently secure the minimum in total energy consumption. In contrast, the SP-SCA and SP-SA baselines exhibit massive

6.2. Algorithmic resilience against traffic stress We evaluate the algorithmic resilience by varying the stress factor ξ ∈ [1, 15] on the traffic demand M k . As illustrated in Fig. 3a, the proposed four BCD algorithms consistently dominates the optimal performance frontier. Under low-to-moderate traffic loads (ξ ≤ 7), the BCD-FW and BCD-IPM outperform the baseline heuristics, maintaining minimized normalized objective values as defined in (36). As the network stress factor intensifies (ξ ≥ 9), these algorithms experience slight inflation, caused by navigating a severely congested feasible region. Remarkably, the warm-started algorithms (BCD-FW-WARM and BCD-IPM-WARM) neutralize this degradation, securing superior minima in most cases. This stability stems from the synergy 17

GG

GG

GG

SP-SCA

SP-SCA

SP-SCA

SP-PDA

SP-PDA

SP-PDA

KSP-PDA

KSP-PDA

KSP-PDA

SP-SA

SP-SA

SP-SA

(a) Normalized objective value in (36)

(b) Maximum delay

(c) Total energy consumption

BCD-IPM

GG

GG SP-SCA SP-PDA

SP-SCA SP-PDA KSP-PDA SP-SA

KSP-PDA SP-SA

(d) Energy efficiency (EE)

(e) Jain’s fairness index (JFI)

(f) Convergence trajectory under ξ = 5

Figure 3: Comprehensive performance evaluation of the proposed four BCD algorithms and baselines under varying stress factors.

energy waste. While the failure of SP-SA is driven by its blind, static equal-partitioning of resources, the behavior of SP-SCA exposes the fragility of SCA strategy in heavily congested, nonconvex landscapes. The proximal gradient descent mechanism within SP-SCA struggles with compounding approximation errors and poor local optima, ultimately consuming immense energy (peaking at 3.183 J at ξ = 15). This proves that without the exact gradient-based resource optimization inherent to the BCD framework, the network loses its ability to efficiently translate consumed energy into data transmission.

An insightful phenomenon emerges among the baseline algorithms. The single-path SP-PDA algorithm exhibits remarkable energy efficiency (ranging from 322.80 to 998.04), significantly outperforming the multi-path KSP-PDA (which stagnates between 184.15 and 277.52). This discrepancy exposes the fatal flaw of heuristic multi-path routing: by partitioning traffic across all the K-shortest paths, KSP-PDA inevitably forces flows through secondary links with severe congestion. This flaw penalizes its energy efficiency. In contrast, SP-PDA restricts transmission to the shortest path with the highest channel gain (1/hi j ) and optimizes resource allocation via PDA, thereby minimizing energy waste. Finally, the conventional heuristics experience a complete EE collapse. The GG remains trapped at EE of approximately 112.50 regardless of load, further illustrating its strong reliance on aggressive power and bandwidth consumption. The SP-SCA and SP-SA yield the lowest EE values of approximately 35.26 and 30.82, respectively. This validates previous analyses: while SP-SA is restricted by its equalpartitioning of resources, SP-SCA struggles with compounding approximation errors and poor local optima. Both baselines exhibit immense energy waste.

As illustrated in Fig. 3d, an evaluation of EE across varying network stress unveils the resource-utilization characteristics of the evaluated algorithms. The proposed BCD-FW and BCD-IPM algorithms consistently dominate this metric, maintaining an exceptionally high EE (bounded between a minimum of 383.73 and a maximum of 1028.95) across the entire stress spectrum. Following closely are the warm-started BCD-FWWARM and BCD-IPM-WARM, which exhibit highly competitive EE profiles except at ξ = 7 and ξ = 11, where they degenerate into the initial GG solutions. This profound superiority demonstrates that the exact multi-path routing strategy does not blindly flood the network; instead, it maximizes the throughputto-energy ratio by allocating fractional flows only to the paths with superior transmission rates.

As illustrated in Fig. 3e, the JFI metric provides a microscopic view into the detailed delay distribution and service equity across all |K| commodities. The proposed BCD algorithms 18

consistently achieve near-perfect fairness (JFI > 0.894), with the BCD-FW and BCD-IPM peaking at 0.995. This demonstrates that the routing subproblem model (5) does not merely optimize the maximum delay; it guarantees the delay equalization for each commodity. An intriguing phenomenon emerges within the SP-PDA and KSP-PDA baselines, which initially appear to occupy the top performance tier. While the optimal PDA resource allocation enforces delay equalization across commodities, this fairness is achieved at the severe expense of the inflated delays previously observed in Fig. 3b. In the subsequent tiers, the greedy nature of the GG heuristic becomes evident. Stagnating at a JFI of approximately 0.894 under high stress (ξ ≥ 5), GG exhibits significant delay fluctuations. It prioritizes certain routes while starving others due to its greedy path selection based on the marginal-delay. Finally, the SP-SCA and SP-SA baselines suffer from severe unfairness, achieving similar JFI values bounded between 0.513 and 0.537. This massive disparity in end-to-end delay can be attributed respectively to compounding approximation errors and the equal partition strategy in single-path routing. As illustrated in Fig. 3f, the convergence trajectories of the proposed BCD-FW and BCD-IPM algorithms are evaluated under a medium congestion scenario (ξ = 5). Remarkably, both algorithms exhibit ultra-fast convergence, stabilizing at their respective optima within merely 4 outer iterations. Specifically, as denoted by the star markers, the optimal objective values are secured at the 3rd and 2nd iterations for BCD-FW and BCDIPM, respectively. A deeper inspection of the convergence trajectories reveals the theoretical distinction between the firstorder and second-order algorithms. Initially, BCD-IPM starts at a higher normalized objective value (0.4679) compared to BCD-FW (0.4276). This behavior is expected: unlike FW, which operates on the extreme vertices of the feasible polytope, IPM strictly enforces interior feasibility via barrier penalties. However, empowered by the exact second-order Hessian information, BCD-IPM captures the exact geometric curvature of the objective. BCD-IPM ultimately discovers a superior highquality solution (0.4155 versus 0.4206 of BCD-FW). This validates the advantage of exploiting Hessian-guided search directions in such highly coupled non-linear problems. It is crucial to evaluate the approximation gap between the original objective (3) and its upper bound surrogate (4a). Remarkably, across all 50 test instances, encompassing the Pareto frontier, resilience, scalability, and robustness scenarios, the approximation gap of the proposed BCD framework is predominantly zero. Specifically, the proposed BCD-FW and warmstarted BCD-FW-WARM and BCD-IPM-WARM algorithms achieve a zero gap across all 50 instances. The BCD-IPM yields a non-zero gap in only 4 out of 50 cases, with an average relative gap of 0.017. These minimal approximation gaps stem from the topological nature of the network bottleneck. In the vast majority of congestion scenarios, the maximum delay is dominated by a specific bottleneck commodity that is routed along a single path. Under such single-path bottleneck conditions, the delay relaxation on each commodity perfectly tightens, rendering the approximate objective (4a) exactly equivalent to the original objective (3). This confirms that the pro-

posed upper-bound relaxation is not merely a theoretical convenience, but a highly accurate representation of the underlying routing bottlenecks. In summary, the preceding resilience analysis establishes the decisive superiority of the proposed BCD framework in achieving the optimal delay-energy trade-off. This joint optimization framework maximizes the energy-to-throughput translation efficiency (yielding up to an order of magnitude improvement); it guarantees bounded transmission times (less than 3.78× relative to the baseline GG) while operating at the minimized energy expenditure (up to a 9.14× improvement versus GG). Having validated the algorithms’ extreme resilience under severe traffic congestion, we subsequently transition to scalability experiments. 6.3. Algorithmic scalability against network scale We expand the network load by first varying the number of commodities |K| ∈ [20, 50] in increments of 5, followed by scaling the device node cardinality |V| ∈ [60, 160] in increments of 20. It is imperative to note that SP-PDA and KSPPDA are omitted in subsequent analyses, as they exhibit maximum delays surging drastically above 89 s (e.g., 89.92 s for KSP-PDA at |K| = 30, and 402.41 s for SP-PDA at |K| = 50). As illustrated in Fig. 4a and 4b, BCD-IPM-WARM demonstrates profound scalability. It consistently dominates the EE metrics in most cases while occupying the second performance tier regarding maximum delay. Remarkably, even under the extreme network load of |K| = 50, its maximum delay experiences only a controlled inflation to 5.84 s, while preserving an exceptionally high EE of 765.08. Consequently, compared to the GG baseline at this peak stress, BCD-IPM-WARM confines its maximum delay gap to less than 2.89×, while achieving a massive 8.35× improvement in energy efficiency. BCD-FWWARM exhibits a similar performance, though it degenerates to the initial GG solution at |K| = 45 (yielding identical EE and delay values of 74.80 and 3.02 s, respectively). Furthermore, the BCD-FW and BCD-IPM algorithms exhibit the robust energy efficiency under high-density traffic, guaranteeing up to a 9.52× EE improvement over GG (|K| = 50). They secure reasonable maximum delays, bounding the gap relative to the GG baseline within a worst-case factor of 3.53× (|K| = 50). Conversely, while the GG heuristic manages to maintain a lower maximum delay, it suffers from severely degraded energy efficiency. This exposes the drawback of GG: to compensate for severe single-path congestion bottlenecks, it aggressively consumes energy to brute-force data transmission, resulting in massive and unsustainable resource waste. As illustrated in Fig. 4c, evaluating the number of commodities employing multi-path routing across varying |K| unveils the underlying topological behavior of the proposed algorithms. (Note: single-path baselines consistently yield zero and are therefore omitted). The proposed BCD framework exhibits a highly adaptive multi-path activation strategy. Rather than blindly splitting flows under light traffic, this framework exclusively activates secondary paths when localized bottlenecks experience severe congestion and resource starvation. For instance, in response to escalating network loads, the BCD-IPM 19

No. of multipath routing

BCD-IPM

GG

SP-SA

GG SP-SCA SP-SA

K

|

K

|

|

(a) Energy efficiency vs. |K|

(b) Maximum delay vs. |K|

GG

EE (Mbits/J)

BCD-IPM-WARM

SP-SCA SP-SA

GG SP-SCA SP-SA

|

(c) Multipath count vs. |K|

Computational time (ms)

EE (Mbits/J)

SP-SCA

GG SP-SCA SP-SA

V

|

V

|

|

(d) Energy efficiency vs. |V|

(e) Maximum delay vs. |V|

|

(f) Computational time vs. |V|

Figure 4: Scalability performance evaluation of the proposed BCD algorithms and baselines.

scales its multi-path routing from 1 (at |K| = 20) to a peak of 32 (at |K| = 40), before stabilizing at 21 under extreme stress (|K| = 50). Conversely, BCD-FW maintains an exceptionally low multi-path activation quantity (peaking marginally at 1 when |K| = 45). This algorithm prefers to explore singlepath routing to alleviate the energy consumption instead of activating multi-path routing. This stems from the nature of the first-order FW mechanism, which relies on sequential shortestpath evaluations over the gradient direction. The BCD-FW-WARM is constrained by initial single-path solution, yielding exactly 0 across all instances. However, BCD-IPM-WARM exhibits a superior capability to explore multi-path routing. Even when anchored by the GG initialization, it overcomes the local bottleneck, activating 2 and 20 multi-paths at |K| = 35 and |K| = 50, respectively. This stems from the exploitation of the exact second-order Hessian matrix coupled with the logarithmic barrier penalties within the LRPDIPM (Algorithm 2). These techniques capture the non-linear curvature of the congestion landscape, pushing the flow distribution away from the heavily congested solutions and deeper into the multi-path routing space. This unique characteristic guarantees the exceptional algorithmic stability previously validated in Fig. 4a and 4b. As illustrated in Fig. 4d and 4e, a comprehensive evaluation of EE and maximum delay is conducted under a scaling

|V| ∈ [60, 160]. As |V| expands, the average node degree surges from 6.1 at |V| = 60 to 15.8 at |V| = 160. This generates a highly dense topological structure. For the EE metric, the BCD-FW and BCD-IPM algorithms consistently define the optimal performance tier, securing peak EE values up to 616.48 and 622.46, respectively. The warm-started BCD-FW-WARM and BCD-IPM-WARM closely follow in the second tier; however, they degenerate into the initial GG solutions at specific topologies (|V| = 60, 120, and 140), yielding identical performance metrics. Conversely, the GG heuristic stagnates at an inferior EE tier, exhibiting a worst-case gap of 3.99× compared to the BCD-IPM optimal bound (specifically at |V| = 80). Furthermore, SP-SCA and SP-SA yield the lowest performance, demonstrating a significant inefficiency in translating consumed energy into effective data transmission. Concurrently, Fig. 4e illustrates the algorithmic performance in maximum delay. While the GG heuristic forces a lower bound in maximum delay by penalizing energy efficiency, the proposed BCD-FW and BCD-IPM achieve comparable transmission times. They maintain a controlled gap of no more than 3.00× relative to GG (peaking at |V| = 60 for BCD-FW). Finally, the BCD-FW-WARM and BCD-IPM-WARM deliver a highly competitive intermediate delay performance while preserving massive energy advantages over the baselines. As illustrated in Fig. 4f, an evaluation of computational 20

60 3.0

5

1.5 1.0

(b) Maximum delay distributions

-S

A

D A

P S

P K S

S

P

-P

-P

D A

C A

G G

-S P S

M P

A R M -W M

P

W

A R M -W W

B C D -I

B C D -I

B C D -I

B C D -F

B C D -F

P

-S

A

D A K S

S

P

-P P S

S

-P

D A

C A

G G

P

-S

M P

A R M -W M

P

B C D -I

A R M

A -S P

2.0

0.0

B C D -F

P K S

S

-P

D A

C A S

P

-P

-S P S

P

A R M -W

B C D -I

M P

B C D -I

A R M -W W

B C D -F

B C D -F

2.5

0.5

W

0 D A

0 G G

10

M

1

(a) Normalized objective value distributions

Total energy consumption (J)

20

-W

2

30

W

Maximum delay (s)

3

40

B C D -F

4

W

Normalized objective value

50

(c) Total energy consumption distributions

Figure 5: Algorithmic robustness evaluation across 20 randomly generated network topologies under the same scenario where |V| = 60, |K| = 20, ξ = 5 and α = 0.4.

scalability is conducted by recording the execution times (in milliseconds) across varying |V| ∈ [60, 160]. As anticipated, a clear gap in computational times emerges between the proposed BCD algorithms and baselines. The BCD-IPM incurs the most substantial computational overhead, peaking at an extreme 551, 196 ms under dense scenarios (|V| = 140). This is driven by the demanding complexity of computing and inverting the exact second-order Hessian matrix within highly dense topologies. However, the proposed BCD-FW-WARM emerges as a profound computational breakthrough. By anchoring the BCD-FW with a high-quality initial solution, and exclusively relying on computationally lightweight first-order gradient information, BCD-FW-WARM achieves its peak execution time of merely 4328 ms at (|V| = 160). Crucially, across the entire evaluated instances, this warm-started algorithm consistently executes faster than the baseline SP-SCA (5588 ms at |V| = 160). This algorithm even compresses the computational discrepancy with the GG heuristic to a single order of magnitude. This validates that BCD-FW-WARM achieves the ultimate synthesis of near-optimal solution and real-world deployment feasibility. Finally, while the heuristic baselines (GG and SP-SA) execute in strictly under 162 ms, this is achieved by abandoning exact iterative resource optimization and spatial multi-path routing. As demonstrated in the preceding resilience analyses, such compensations inevitably lead to severe resource waste, restricting their practical deployment. In summary, the comprehensive scalability analysis establishes the profound robustness and superiority of the proposed BCD framework. As the number of commodities or the overall topological scale expands, the BCD-FW and BCD-IPM algorithms guarantee near-optimal solutions. They bound the maximum delay within a controlled gap (less than 3.53×) relative to the baseline GG, while sustaining the peak in energy efficiency (up to a 9.52× improvement versus GG). Among the proposed algorithms, BCD-FW-WARM emerges as a computational breakthrough. By leveraging computationally lightweight firstorder gradient information under the high-quality initialization, it delivers a near-optimal solution within seconds. Conversely,

the GG baseline exposes the fatal flaw of heuristic designs. While GG offers computationally cheap execution and competitive maximum delay, it aggressively employs immense transmission power and bandwidth to force data through congested single-path bottlenecks. These results prove that the proposed BCD framework stands as a profoundly sound paradigm, capable of ultra-reliable, highly scalable, and energy-efficient communications in dense overlay D2D networks. 6.4. Algorithmic robustness against topological uncertainty To evaluate the algorithmic robustness against topological uncertainty, we conduct comprehensive empirical evaluations across a diverse set of 20 independently generated random topologies. The system settings are identical to the settings in Section 6.2 with ξ = 5. As illustrated in the boxplots of Fig. 5a, the distribution of the normalized objective values (36) establishes the superiority and stability of the proposed BCD framework. The interquartile ranges (IQRs) of the four BCD algorithms exhibit profound compression compared to the baselines. This minimized statistical dispersion demonstrates that the proposed algorithms maintain optimal delay-energy tradeoffs across random topologies. Furthermore, the upper whiskers and outlier data points reveal a critical performance boundary. While the BCD-FW algorithm experiences a localized objective inflation (peaking at 0.772 in Topology 12), the other three BCD algorithms completely neutralize this vulnerability. Specifically, the extreme values of BCD-IPM and both warmstarted algorithms are bounded below 0.668, outperforming the worst-case value of 0.722 yielded by the GG heuristic. In contrast, heuristic baselines such as KSP-PDA exhibit severe volatility, with extreme outliers surging to 5.09 (Topology 10) due to blind multi-path routing in random topologies. While the single-path GG and SP-PDA baselines maintain compact IQRs comparable to the proposed BCD framework, their median objective values and overall distributional boundaries are elevated. This demonstrates their inferior capabilities in achieving the optimal delay-energy trade-off. To further deconstruct the objective values, Fig. 5b and 5c depict the distributions of its two components: maximum delay 21

and total energy consumption, respectively. The proposed BCD framework consistently secures highly compact IQRs across both objectives, with merely a single outlier recorded for BCDFW-WARM in maximum delay (specifically 6.69 s at Topology 20). This profound consistency validates the superior robustness against topological uncertainty. While the GG and SP-SCA baselines exhibit compact IQRs regarding maximum delay, they yield relatively elongated boxes and higher median values for total energy consumption. This corroborates the conclusions drawn in Sections 6.2 and 6.3: these heuristic baselines achieve competitive transmission times through an aggressive over-consumption of transmission power and bandwidth. Furthermore, all algorithms integrating the optimal PDA algorithm (including the BCD framework, SP-PDA, and KSPPDA) consistently exhibit remarkably compressed IQRs in total energy consumption (e.g., SP-PDA consistently operating below 0.18 J). This validates that the PDA algorithm consistently bounds the energy waste, regardless of the underlying topological structures. However, as evidenced by the severe delay fluctuations of KSP-PDA in Fig. 5b (surging to 54.22 s at Topology 10), simply applying optimal resource allocation atop heuristic multi-path routing fails to stabilize overall network performance. This contrast validates the necessity of the proposed BCD framework, jointly optimizing the multi-path routing via MF-FW/LR-PDIPM and resource allocation via PDA algorithm.

BCD-IPM emerges as the most robust variant, guaranteeing the optimal or near-optimal trade-off even under extremely dense topologies. This algorithm achieves a maximum 9.14× reduction in total energy consumption and up to an order of magnitude improvement in energy efficiency, while maintaining a bounded maximum delay gap (up to 3.78×) relative to the best baseline GG. Meanwhile, the BCD-FW-WARM algorithm distinguishes itself as a highly practical engineering solution, synergizing heuristic initializations with lightweight first-order updates to deliver near-optimal solutions in mere seconds. Appendix A. Theoretical proofs for resource allocation Appendix A.1. Proof of Lemma 1 Proof. The proof proceeds in two distinct parts, establishing the joint convexity of the generalized perspective function and the strict joint convexity of the specific bivariate composite function, respectively. Joint convexity of the perspective function P(u, v). Given a convex function g(·) and the perspective mapping P(u, v) = vg(u/v) for v > 0, we introduce the auxiliary variable z = u/v. Computing the second-order partial derivatives yields the Hessian matrix of P(u, v): # " g′′ (z) 1 −z 2 . ∇ P(u, v) = v −z z2

7. Conclusion

To establish joint convexity, we verify that the Hessian is positive semi-definite (∇2 P ⪰ 0). For any arbitrary non-zero vector w = [w1 , w2 ]T , the associated quadratic form evaluates to:

In this paper, we investigated the joint routing and resource allocation problem in overlay D2D networks. Because the network-layer traffic flows and physical-layer resources are intricately intertwined via the non-convex Shannon capacity, the resulting formulation is computationally intractable. To resolve this, we proposed a novel BCD framework to iteratively optimize the routing and resource domains. For the routing subproblem, we first employed LSE smoothing to approximate the non-differentiable minimax delay objective. For the requirement of rapid execution, the MF-FW algorithm was developed to guide the descent direction and optimize step sizes based on the directional curvature. For scenarios demanding highprecision solutions, the LR-PDIPM was developed. By exploiting the Sherman-Morrison rank-1 formula to accelerate the Newton system inversion, this method reduces the per-iteration complexity from O(|K|3 |V|3 ) to O(|K||V|3 ). For the physicallayer resource allocation subproblem, we exploited the inherent convexity by utilizing a time-domain perspective transformation. To bypass the scalability bottlenecks of standard commercial solvers, we designed a PDA algorithm. This approach updates Lagrangian multipliers via gradients and solves the 2D link-level subproblems using an inner bisection search. Theoretical analysis guarantees that the proposed BCD framework converges to an ε-neighborhood of a stationary point in polynomial time. Through comprehensive empirical evaluations encompassing resilience, scalability and robustness, the proposed BCD framework establishes its superiority in achieving the optimal delay-energy trade-off. Specifically,

wT ∇2 P(u, v)w =

g′′ (z) 2 g′′ (z) (w1 −2zw1 w2 +z2 w22 ) = (w1 −zw2 )2 . v v

Since g(·) is convex, its second derivative satisfies g′′ (z) ≥ 0. Coupled with the domain definition v > 0 and the non-negative squared term (w1 −zw2 )2 ≥ 0, the quadratic form is non-negative for all w. Consequently, ∇2 P ⪰ 0, proving that P(u, v) is jointly convex over its domain. Strict and joint convexity of f (u, v). We now analyze the bivariate function f (u, v) = uv(21/uv − 1) for positive variables u, v > 0. Let x = uv and define the univariate core function q(x) = x(21/x − 1), such that f (u, v) = q(uv). Applying the multivariate chain rule, the Hessian matrix of f (u, v) is formulated as: " # v2 q′′ (x) xq′′ (x) + q′ (x) 2 ∇ f (u, v) = . xq′′ (x) + q′ (x) u2 q′′ (x) The first and second derivatives of q(x) are given by: q′ (x) = 21/x − 1 −

ln 2 1/x 2 , x

q′′ (x) =

(ln 2)2 1/x 2 . x3

(A.1)

To prove strict joint convexity, we must guarantee that ∇2 f (u, v) is positive definite (∇2 f ≻ 0), which holds if and only if both its trace and determinant are positive. 22

P P expµTk / j∈K expµT j satisfies βk > 0 and k∈K βk = 1. The Hessian matrix is derived as:  T    X    X X ∇2 F(x) = αµ  βk ∇T k ∇T kT −  βk ∇T k   βk ∇T k   .

1) Trace analysis: Evaluating the trace of the Hessian yields: Tr(∇2 f ) = (u2 + v2 )q′′ (x). Since u, v > 0 =⇒ x > 0, equation (A.1) guarantees q′′ (x) > 0. Thus, Tr(∇2 f ) > 0. 2) Determinant analysis: The determinant of the Hessian is computed as:

k∈K

k∈K

k∈K

1) Positive semi-definiteness: For any arbitrary vector v, the quadratic form is given by:   2  X X    T 2 2 T T v ∇ F(x)v = αµ  βk (∇T k v) −  βk ∇T k v  . (B.1)

det(∇2 f ) = u2 v2 (q′′ (x))2 − (xq′′ (x) + q′ (x))2   = x2 (q′′ (x))2 − x2 (q′′ (x))2 + 2xq′′ (x)q′ (x) + (q′ (x))2   = −q′ (x) 2xq′′ (x) + q′ (x) .

k∈K

k∈K

P

Since k∈K βk = 1, we can apply Jensen’s inequality to the P P
2 ≤ k∈K βk z2k function f (z) = z2 , which yields k∈K βk zk by setting zk = ∇T kT v. Consequently, vT ∇2 F(x)v ≥ 0. Thus, ∇2 F(x) ⪰ 0, ensuring the convexity of F(x). 2) Lipschitz continuity: The L-smoothness constant corresponds to the spectral norm ∥∇2 F(x)∥2 . To bound this, we analyze the Rayleigh quotient for any unit vector v (i.e., ∥v∥2 = 1). From (B.1), we obtain an upper bound: X vT ∇2 F(x)v ≤ αµ βk (∇T kT v)2 .

To determine the sign of this product, we define the positive auxiliary variable ω = ln 2/x > 0, mapping the exponential term to 21/x = expω . First, we analyze the sign of q′ (x). Let ϕ(ω) = expω (1−ω)−1. Its derivative is ϕ′ (ω) = −ω expω < 0 for all ω > 0. Given the initial condition ϕ(0) = 0, it follows that ϕ(ω) < 0. Therefore, q′ (x) < 0, which yields −q′ (x) > 0. Second, we analyze the sign of the bracketed term S (x) = 2xq′′ (x) + q′ (x). Substituting ω, we obtain the mapped function ψ(ω) = expω (2ω2 −ω+1)−1. Its derivative evaluates to ψ′ (ω) = expω (2ω2 + 3ω) > 0 for ω > 0. Given ψ(0) = 0, it follows that ψ(ω) > 0, implying S (x) > 0. Since both decomposed factors are positive, their product guarantees det(∇2 f ) > 0. Concurrently satisfying Tr(∇2 f ) > 0 and det(∇2 f ) > 0, the Hessian matrix is positive definite. This completes the proof that f (u, v) is strictly jointly convex.

k∈K

Applying the Cauchy-Schwarz inequality, (∇T kT v)2 ≤ ∥∇T k ∥22 ∥v∥22 = ∥∇T k ∥22 . This yields: vT ∇2 F(x)v ≤ αµ

X

βk ∥∇T k ∥22

k∈K

Appendix A.2. Proof of Theorem 1

≤ αµ max ∥∇T k ∥22

Proof. To establish (23) as a convex optimization problem, we verify the convexity of its objective and feasible region. 1) Objective function: The objective (23a) is a sum of a linear term αT and energy functions Ei j (li j , ti j ) = hNi0j li j ti j (2mi j /(li j ti j ) −1). Via the affine mapping u = li j and v = ti j /mi j , we obtain Ei j ∝ f (u, v). By Lemma 1, f (u, v) = uv(21/uv − 1) is strictly jointly convex for positive (u, v). Since convexity is preserved under affine mapping and positive scaling, the objective is convex. 2) Feasible region: Delay and bandwidth constraints (23b)(23c) define affine spaces. The power constraints (23d) depend on pi j (li j , ti j ) = hNi0j li j (2mi j /(li j ti j ) − 1). Let z(ti j ) = mi j /ti j , which is convex for ti j > 0. The perspective P(z, li j ) = li j (2z/li j − 1) is jointly convex for li j > 0 and monotonically increasing in z as ∂P ∂z > 0. By the scalar composition theorem, substituting the convex z(ti j ) into P(z, li j ) preserves joint convexity. Thus, the feasible region is an intersection of convex sets, completing the proof.

k∈K

X

βk = αµ max ∥∇T k ∥22 .

k∈K

k∈K

Since T k (x) is linear, ∇T k is a constant vector of path delay coefficients. Therefore, the spectral norm is bounded by:   X M k 2   = αµC 2 . ∥∇2 F(x)∥2 ≤ αµ max max k∈K ri j  (i, j)∈E

2 Setting L = αµCmax completes the proof.

Appendix B.2. Proof of Theorem 2 Proof. Let hn = F(x(n) ) − F(x∗ ) denote the primal optimality gap at iteration n, and let C = LD2 . Because the objective function F(x) is L-smooth (as established in Lemma 2), the standard descent lemma bounds the objective value at the next iteration for any γ ∈ [0, 1]: γ2 L∥y(n) − x(n) ∥22 . 2 (B.2) Due to D = maxx,y∈X ∥x − y∥2 , we have ∥y(n) − x(n) ∥22 ≤ D2 . Furthermore, due to the convexity of F(x) and the optimality of y(n) under current gradient ∇F(x(n) ), the gradient term bounds the negative primal gap: ∇F(x(n) )T (y(n) − x(n) ) ≤ F(x∗ ) − F(x(n) ) = −hn . F(x(n+1) ) ≤ F(x(n) ) + γ∇F(x(n) )T (y(n) − x(n) ) +

Appendix B. Theoretical proofs for MF-FW Appendix B.1. Proof of Lemma 2 Proof. The curvature of the smoothed objective F(x) is governed by its LSE component. By the multivariate chain P rule, the gradient is ∇F(x) = α k∈K βk ∇T k (x), where βk = 23

Substituting these bounds into (B.2) yields: F(x(n) + γd(n) ) ≤ F(x(n) ) − γhn +

Recall from the proof of Lemma 2 that the Hessian is written as:    X T T 2  ¯ ¯ ∇ F(x) = αµ  βk ∇T k ∇T k − ∇T ∇T  ,

γ2 C. 2

k∈K

¯ = Pk∈K βk ∇T k . where ∇T To obtain the third-order derivative, we compute the directional derivative of ∇2 F(x) along an arbitrary unit vector d ∈ R|E||K| (with ∥d∥2 = 1). By applying the chain rule, the ¯ ). Consederivative of βk is given by ∇βk = µβk (∇T k − ∇T quently, the third-order directional derivative is calculated as X   ¯ )(∇T k −∇T ¯ )T (∇T k − ∇T ¯ )T d . βk (∇T k −∇T ∇3 F(x)[d] = αµ2

The truncated Newton step size λ∗ is the exact minimizer of the local quadratic surrogate defined by the directional curvature C (n) . Because F(x) is L-smooth over the compact domain X, the local curvature is upper-bounded by the global curvature constant, i.e., C (n) ≤ C. Thus, for ∀γ ∈ [0, 1], F(x(n) + λ∗ d(n) ) ≤ F(x(n) ) − γhn +

γ2 C. 2

k∈K

Then, γ C. (B.3) 2 This inequality guarantees a primal gap hn+1 that is no greater than the gap produced by the standard diminishing step size γ(n) = 2/(n + 2). Therefore, it is sufficient to prove the theoretical upper bound hn ≤ 2C/(n + 1) via mathematical induction using the recurrence (B.3) generated by γ(n) . Base case (n = 0): For the initial iteration, applying γ(0) = 1 to (B.3) gives h1 ≤ (1 − 1)h0 + C2 (1)2 = C2 . The theorem 2C necessitates h1 ≤ 0+2 = C. Since C2 ≤ C strictly holds, the base case is verified. Inductive step: Assume the bound holds for iteration n, i.e., hn ≤ 2C/(n + 1). Substituting this inductive hypothesis and γ(n) = 2/(n + 2) into (B.3) for hn+1 yields: hn+1 ≤ (1 − γ)hn +

2

Then, ∇3 F(x)[d, d, d] = αµ2

X

  ¯ )T d 3 . βk (∇T k − ∇T

k∈K

P

Recall βk ∈ (0, 1) and βk = 1. Applying the triangle inequality and the standard operator norm bound, we obtain: sup ∇3 F(x)[d, d, d] ≤ 2αµ2 max ∥∇T k ∥32 . k∈K

∥d∥2 =1

(C.1)

Given that T k (x) is linear, its gradient vector ∇T k ∈ R|K||E| is highly sparse and the non-zero entries are the constant M k /ri j . The L2 -norm of the gradient is bounded by the network scale, √ specifically ∥∇T k ∥2 ≤ |E|Rmax and Rmax = maxk,(i, j) (M k /ri j ). Substituting this bound into (C.1) yields:

! !2 2 2C C 2 + n+2 n+1 2 n+2 2C 2C n · + = n + 2 n + 1 (n + 2)2 " # 2C n 1 = + . n+2 n+1 n+2

hn+1 ≤ 1 −

∥∇3 F(x)∥2 ≤ 2αµ2

p

3 |E|Rmax .

Setting the Lipschitz constant to LH = 2αµ2 |E|3/2 R3max , we obtain LH = O(αµ2 |E|3/2 R3max ). This guarantees that the Hessian is Lipschitz continuous, revealing that the Newton-step curvature variation is governed jointly by the smoothing parameter and network scale.

Analyzing the bracket term, we derive:

Appendix C.2. Proof of Lemma 4 (n) Proof. Let µ̃ = µgap denote the current duality gap. For any iterate within the wide neighborhood N−∞ (γ), the variables satisfy the centrality condition xi si ≥ γµ̃. 1) Bounding the second-order terms. From the Newton system of the PDIPM, the search directions satisfy: S∆x + X∆s = −Xs + σµ̃1. (C.2)

1 n(n + 2) + (n + 1) n2 + 3n + 1 n + = = 2 < 1. n+1 n+2 (n + 1)(n + 2) n + 3n + 2 Therefore, 2C . n+2 This completes the mathematical induction, establishing the global O(1/n) sublinear convergence rate of MF-FW algorithm. hn+1 ≤

Because the linear equality constraints dictate A∆x = 0, the primal direction lies in the null space of A. Correspondingly, the dual direction satisfies ∆s = ∇2 F(x)∆x − AT ∆y. Therefore, the inner product yields ∆xT ∆s = ∆xT ∇2 F(x)∆x. Since the Hessian is positive semi-definite (as established in Lemma 2), we have ∆xT ∆s ≥ 0. By dividing both sides of (C.2) by (XS)1/2 and applying the triangle inequality along with the wide neighborhood lower bound xi si ≥ γµ̃, standard interior-point algebraic manipulations [59] ensure that the scaled Newton directions are bounded.

Appendix C. Theoretical proofs for LR-PDIPM Appendix C.1. Proof of Lemma 3 Proof. By the Mean Value Theorem for matrix-valued functions, establishing the Lipschitz continuity of the Hessian matrix ∇2 F(x) is equivalent to bounding the spectral norm of the third-order derivative ∇3 F(x) over the feasible domain X. 24

Consequently, the Euclidean norm of the product ∆X∆s is bounded by the duality gap: X (∆xi ∆si )2 ≤ Cµ̃2 , ∥∆X∆s∥22 =

where the single-iteration reduction factor is bounded by δ = λmin (1 − σmax ). According to the standard path-following interior-point theory [59], to maintain the iterates within the neighborhood and ensure the quadratic boundedness of the cross terms, the allowable step size λmin scales inversely with the square root of the inequality dimension, yielding δ = √ Ω(1/ |K||E|). 2) Iteration complexity bound. Applying the recursive relation (C.4) over n iterations yields:

i

where C > 0 is a constant dependent on γ and the network scale, independent of the iteration count. 2) Guaranteeing a positive step size. To maintain the iterates within N−∞ (γ) at the next iteration, the condition x(λ) ◦ s(λ) = (x + λ∆x) ◦ (s + λ∆s) ≥ γµ̃(λ)1 must hold. Because the Hessian is LH -Lipschitz continuous (Lemma 3), the deviation of the nonlinear KKT residuals (13) from the linear Newton prediction is bounded by O(λ2 LH ∥∆x∥22 ). The bounded cross terms and the finite Lipschitz constant LH guarantee that the neighborhood violation grows quadratically with λ. Thus, the backtracking line search will accept a step size λ before it shrinks to zero, proving the existence of a lower bound λ ≥ λmin > 0. 3) Monotonic gap reduction. The duality gap at the candidate step λ is calculated as:

µ(n) ≤ (1 − δ)n µ(0) . To achieve an ε-optimal solution, the algorithm must satisfy µ(n) ≤ ε. Taking the natural logarithm of both sides gives: ! ε n ln(1 − δ) ≤ ln (0) . µ Using the standard logarithmic inequality ln(1 − δ) ≤ −δ for δ ∈ (0, 1), we obtain the required number of iterations: ! 1 µ(0) n ≥ ln . δ ε

(x + λ∆x)T (s + λ∆s) xT s + λ(xT ∆s + sT ∆x) + λ2 ∆xT ∆s = . |K||E| |K||E| √ (C.3) Since δ−1 = O( |K||E|), and the initial gap µ(0) is a given conMultiplying (C.2) by 1T , we obtain the sum of the cross terms: stant, the maximum number of iterations is bounded by: xT ∆s + sT ∆x = −xT s + σµ̃|K||E| = −|K||E|µ(1 − σ). Substituting !! p 1 this into (C.3) yields: . n = O |K||E| ln ε 2
λ ∆xT ∆s. µ̃(λ) = µ̃ 1 − λ(1 − σ) + Because the smoothed objective F(x) is convex over the fea|K||E| sible domain X, the satisfaction of the ε-KKT residual condiThe backtracking line search controls the second-order term tions guarantees that the limit of the iterate sequence is a global λ2 ∆xT ∆s. Enforcing the search criteria yields the upper bound ε-optimal solution, concluding the proof. for the updated gap: µ̃(λ) =

(n) µ(n+1) gap ≤ 1 − λmin (1 − σmax ) µgap .

Appendix D. Theoretical proofs for PDA Appendix D.1. Proof of Lemma 5

This confirms the monotonic reduction of the optimality gap, concluding the proof.

Proof. The time complexity of the LR-PDIPM is determined by the dimension of the inequality constraints and the reduction rate of the duality gap. We establish the proof in two sequential steps. 1) Gap reduction rate. For the smoothed routing subproblem (5), the inequality bounds originate from the non-negativity constraints x ≥ 0. Since x encompasses the flow of each commodity k ∈ K over each edge (i, j) ∈ E, the total number of inequality constraints is |K||E|. By initializing the IPM at (x(0) , s(0) ) > 0 within the wide neighborhood N−∞ (γ), Lemma 4 guarantees that the algorithm generates a sequence of iterates that are strictly feasible with respect to the inequality bounds. Furthermore, it establishes a monotonic reduction of the duality gap:

Proof. To establish the convergence of the alternating coordinate descent for the link-level subproblem, we analyze the structural properties of Φi j defined in (24). 1) Strict convexity and differentiability. The subproblem objective is a positive linear combination: P Φi j (li j , ti j ) = (1 − α)Ei j + φi pi j + ωli j + ( k∈K θk ρkij )ti j . By Lemma 1 and Theorem 1, the composite energy function Ei j is strictly jointly convex, and the transmission power function pi j is jointly convex with respect to (li j , ti j ). Given that α ∈ (0, 1) and the dual multipliers are non-negative, the strict convexity is preserved in their sum. Thus, Φi j is strictly convex and continuously differentiable. The feasible domain Fi j = QoS [lmin , B] × [tmin , T max ] is a compact Cartesian product domain, ensuring the constraints on li j and ti j are completely decoupled. 2) Monotonic descent. (q) Let z(q) = (li(q) j , ti j ) denote the iterate at the q-th inner iteration. The exact 1D bisection updates yield:

µ(n+1) ≤ (1 − δ) µ(n) ,

(q) (q) (q−1) Φi j (li(q) ) ≤ Φi j (li(q−1) , ti(q−1) ). j , ti j ) ≤ Φi j (li j , ti j j j

Appendix C.3. Proof of Theorem 4

(C.4) 25

Due to the strict convexity of Φi j , the sequence of objective values {Φi j (z(q) )} is monotonically decreasing until the minimum is reached. Bounded below over the compact set Fi j , the sequence is guaranteed to converge to a limit point. 3) Uniqueness of the limit point. For continuously differentiable functions optimized over a Cartesian product domain, any limit point of the coordinate descent sequence is guaranteed to be a stationary point. Because Φi j is strictly convex, this point is the global minimum. Therefore, the sequence converges to the unique minimum, concluding the proof.

dictated by the network topology (|Vout (i)|, |Eact |) and the finite domain bounds (Pbound /Pmax ). Therefore, their L2 -norm is i uniformly bounded, verifying the existence of a finite constant G > 0 such that ∥g(n) ∥2 ≤ G. References [1] A. Asadi, Q. Wang, V. Mancuso, A survey on device-todevice communication in cellular networks, IEEE Communications Surveys & Tutorials 16 (4) (2014) 1801– 1819. doi:10.1109/comst.2014.2319555. [2] G. P. de Freitas Cardoso, P. H. P. De Carvalho, P. R. de Lira Gondim, Joint spectrum allocation and power control for D2D communication and sensing in 6G networks using DRL-based hyper-heuristics, Computer Networks 276 (2026) 111969. doi:10.1016/j.comnet.2025.111969.

Appendix D.2. Proof of Lemma 6 Proof. To prove the uniform boundedness of the dual gradient ∥g(n) ∥2 , we formulate its normalized components based on the dual update rules (25) and bound them individually. 1) Bounding normalized linear residuals. According to the bandwidth update (25a) and the routing delay update (25c), the corresponding gradient components are formulated as: P (n) (i, j)∈Eact li j − B (n) , (D.1) gω = X B g(n) = ρkij ti(n) ∀k ∈ K. (D.2) j , θk

[3] R. I. Ansari, C. Chrysostomou, S. A. Hassan, M. Guizani, S. Mumtaz, J. Rodriguez, J. J. Rodrigues, 5G D2D networks: Techniques, challenges, and future prospects, IEEE Systems Journal 12 (4) (2017) 3970–3984. doi:10.1109/jsyst.2017.2773633. [4] Y.-C. Wang, W.-T. Chen, DACS: Efficient resource allocation and power control for D2D communication considering RAN sharing, Computer Networks 274 (2025) 111799. doi:10.1016/j.comnet.2025.111799.

(i, j)∈Eact

During the primal update, the allocated bandwidth variables li(n) j are restricted within [lmin , B]. Thus, the maximal theoretical sum is bounded by |Eact |B, yielding a strict bound on the relative bandwidth violation: |g(n) ω | ≤ |Eact |. Similarly, because the time QoS variables are bounded by the QoS delay tolerance ti(n) j ≤ T max , and the routing fractions ρkij ∈ [0, 1], the delay gradient is

[5] J. Xie, Y. Jia, W. Wen, Z. Chen, L. Liang, Dynamic D2D multihop offloading in multi-access edge computing from the perspective of learning theory in games, IEEE Transactions on Network and Service Management 20 (1) (2022) 305–318. doi:10.1109/tnsm.2022.3201470.

QoS bounded by a finite constant: |g(n) | ≤ |Eact |T max . θk 2) Bounding the normalized power residual. According to the node power update (25b), the normalized power gradient for node i is formulated as:

P gφ(n)i =

(n) (n) max j∈Vout (i) pi j (li j , ti j ) − Pi . Pmax i

[6] Y. Deng, H. Zhang, X. Chen, Y. Fang, Multi-hop task routing in vehicle-assisted collaborative edge computing, IEEE Transactions on Vehicular Technology 73 (2) (2023) 2444–2455. doi:10.1109/tvt.2023.3312142. [7] W. Huang, Z. Zhao, G. Min, J. Chen, Distributed multihop task offloading in massive heterogeneous IoT systems, IEEE Transactions on Computers 73 (4) (2024) 1126– 1137. doi:10.1109/tc.2024.3355767.

(D.3)

(n) During the primal update, the resource variables (li(n) j , ti j ) are

[8] I. Behnke, H. Austad, Real-time performance of industrial IoT communication technologies: A review, IEEE Internet of Things Journal 11 (5) (2023) 7399–7410. doi:10.1109/jiot.2023.3332507.

QoS restricted within the compact set Fi j = [lmin , B] × [tmin , T max ]. Because the transmission power pi j (l, t) is a continuous function over these bounded variables, its output cannot diverge to infinity. Consequently, the maximal power consumed by any single link is bounded by a finite constant, defined as Pbound = max(l,t)∈Fi j pi j (l, t). Thus, the maximal theoretical power sum for node i is bounded by |Vout (i)|Pbound , yielding a bound on the relative power violation:

|g(n) φi | ≤

[9] N. Hu, Z. Tian, X. Du, N. Guizani, Z. Zhu, Deepgreen: A dispersed energy-efficiency computing paradigm for green industrial IoT, IEEE Transactions on Green Communications and Networking 5 (2) (2021) 750–764. doi:10.1109/tgcn.2021.3064683.

|Vout (i)|Pbound + 1. Pmax i

[10] J. Zhao, F. Shen, J. Joung, Throughput maximization with rate-dependent power consumption in battery-limited multiuser networks, IEEE Transactions on Vehicular Technology 69 (1) (2019) 1141–1146. doi:10.1109/tvt.2019.2953711.

3) Global gradient bound. Equations (D.1)-(D.3) demonstrate that every scalar element of the normalized gradient vector g(n) is bounded by constants 26

[11] H. T. Nguyen, H. D. Tuan, T. Q. Duong, H. V. Poor, W.J. Hwang, Collaborative multicast beamforming for content delivery by cache-enabled ultra dense networks, IEEE Transactions on Communications 67 (5) (2019) 3396– 3406. doi:10.1109/tcomm.2019.2894797.

[21] Y. Wu, J. Wang, L. Qian, R. Schober, Optimal power control for energy efficient D2D communication and its distributed implementation, IEEE Communications Letters 19 (5) (2015) 815–818. doi:10.1109/lcomm.2015.2407871.

[12] P. Chen, X. Zhou, J. Zhao, F. Shen, S. Sun, Energyefficient resource allocation for secure D2D communications underlaying UAV-enabled networks, IEEE Transactions on Vehicular Technology 71 (7) (2022) 7519–7531. doi:10.1109/tvt.2022.3168277.

[22] F. Wang, C. Xu, L. Song, Z. Han, Energy-efficient resource allocation for device-to-device underlay communication, IEEE Transactions on Wireless Communications 14 (4) (2014) 2082–2092. doi:10.1109/twc.2014.2379653.

[13] A. Paul, S. P. Maity, Machine learning for spectrum information and routing in multihop green cognitive radio networks, IEEE Transactions on Green Communications and Networking 6 (2) (2021) 825–835. doi:10.1109/tgcn.2021.3127308.

[23] C. Liu, X. Wang, X. Wu, J. Guo, Economic scheduling model of microgrid considering the lifetime of batteries, IET Generation, Transmission & Distribution 11 (3) (2017) 759–767. doi:10.1049/iet-gtd.2016.0772. [24] L. Tao, J. Ma, Y. Cheng, A. Noktehdan, J. Chong, C. Lu, A review of stochastic battery models and health management, Renewable and Sustainable Energy Reviews 80 (2017) 716–732. doi:10.1016/j.rser.2017.05.127.

[14] D. Feng, L. Lu, Y. Yuan-Wu, G. Y. Li, G. Feng, S. Li, Device-to-device communications underlaying cellular networks, IEEE Transactions on communications 61 (8) (2013) 3541–3551. doi:10.1109/tcomm.2013.071013.120787.

[25] K. Yang, S. Martin, L. Boukhatem, J. Wu, X. Bu, Energyefficient resource allocation for device-to-device communications overlaying LTE networks, in: 2015 IEEE 82nd Vehicular Technology Conference (VTC2015-Fall), IEEE, 2015, pp. 1–6. doi:10.1109/vtcfall.2015.7390922.

[15] J. Zhao, Y. Liu, K. K. Chai, Y. Chen, M. Elkashlan, Joint subchannel and power allocation for NOMA enhanced D2D communications, IEEE Transactions on Communications 65 (11) (2017) 5081–5094. doi:10.1109/tcomm.2017.2741941.

[26] Q.-N. Tran, N.-S. Vo, M.-P. Bui, T.-M. Phan, Q.-A. Nguyen, T. Q. Duong, Spectrum sharing and power allocation optimised multihop multipath D2D video delivery in beyond 5G networks, IEEE Transactions on Cognitive Communications and Networking 8 (2) (2021) 919–930. doi:10.1109/tccn.2021.3133838.

[16] M. S. Al-Abiad, M. Z. Hassan, M. J. Hossain, A joint reinforcement-learning enabled caching and cross-layer network code in F-RAN with D2D communications, IEEE Transactions on Communications 70 (7) (2022) 4400– 4416. doi:10.1109/tcomm.2022.3168058.

[27] S. Liu, G. Yu, D. Wen, X. Chen, M. Bennis, H. Chen, Communication and energy efficient decentralized learning over D2D networks, IEEE Transactions on Wireless Communications 22 (12) (2023) 9549–9563. doi:10.1109/twc.2023.3271854.

[17] A. A. Al-habob, J. Lin, O. A. Dobre, Y. Jing, Minmax latency minimization for energy-constrained multiUAV mobile edge computing, IEEE Transactions on Network Science and Engineering 11 (5) (2024) 4577–4590. doi:10.1109/tnse.2024.3409207.

[28] Y. Ergiz, A. M. Demirtas, T. Girici, Joint multipath flow and layer allocation for scalable video streaming, Computer Networks 191 (2021) 107995. doi:10.1016/j.comnet.2021.107995.

[18] N. Li, W. Hao, F. Zhou, Z. Chu, S. Yang, O. Muta, H. Gacanin, Min–max latency optimization for IRSaided cell-free mobile edge computing systems, IEEE Internet of Things Journal 11 (5) (2023) 8757–8770. doi:10.1109/jiot.2023.3322751.

[29] E. Gures, P. Mach, Z. Becvar, Joint route selection and radio resources allocation for caching in multi-hop networks, IEEE Transactions on Communications (2026). doi:10.1109/TCOMM.2026.3664448.

[19] F. Jameel, Z. Hamid, F. Jabeen, S. Zeadally, M. A. Javed, A survey of device-to-device communications: Research issues and challenges, IEEE Communications Surveys & Tutorials 20 (3) (2018) 2133–2168. doi:10.1109/comst.2018.2828120.

[30] H. B. Valiveti, P. T. Rao, EHSD: An exemplary handover scheme during D2D communication based on decentralization of SDN, Wireless Personal Communications 94 (4) (2017) 2393–2416. doi:10.1007/s11277-0163490-7.

[20] A. Abrardo, M. Moretti, Distributed power allocation for D2D communications underlaying/overlaying OFDMA cellular networks, IEEE Transactions on Wireless Communications 16 (3) (2016) 1466–1479. doi:10.1109/twc.2016.2646360.

[31] S. Kazemi Rashed, R. Shahbazian, S. A. Ghorashi, Learning-based resource allocation in D2D communications with QoS and fairness considerations, Transactions 27

on Emerging Telecommunications Technologies 29 (1) (2018) e3249. doi:10.1002/ett.3249.

[42] I. Budhiraja, N. Kumar, D. Garg, M. Guizani, G. Kaddoum, Joint traffic admission, resource allocation and mode selection protocol for NOMAbased D2D users underlaying cellular network for 5G and beyond networks, IEEE Transactions on Network Science and Engineering 11 (5) (2024) 4371–4383. doi:10.1109/tnse.2024.3418135.

[32] J. Xu, X. Gu, Z. Fan, D2D power control based on hierarchical extreme learning machine, in: 2018 IEEE 29th Annual International Symposium on Personal, Indoor and Mobile Radio Communications (PIMRC), IEEE, 2018, pp. 1–7. doi:10.1109/pimrc.2018.8580872. [33] I. Ioannou, V. Vassiliou, C. Christophorou, A. Pitsillides, Distributed artificial intelligence solution for D2D communication in 5G networks, IEEE Systems Journal 14 (3) (2020) 4232–4241. doi:10.1109/jsyst.2020.2979044.

[43] L. Pu, X. Chen, J. Xu, X. Fu, D2D fogging: An energyefficient and incentive-aware task offloading framework via network-assisted D2D collaboration, IEEE Journal on Selected Areas in Communications 34 (12) (2016) 3887– 3901. doi:10.1109/jsac.2016.2624118.

[34] X. Chen, G. Zhu, Y. Deng, Y. Fang, Federated learning over multihop wireless networks with innetwork aggregation, IEEE Transactions on Wireless Communications 21 (6) (2022) 4622–4634. doi:10.1109/twc.2022.3168538.

[44] Z. Zhou, C. Gao, C. Xu, Joint peer discovery and resource allocation for social-aware D2D communications: A matching approach, in: 2016 IEEE international conference on communication systems (ICCS), IEEE, 2016, pp. 1–6. doi:10.1109/iccs.2016.7833601. [45] M. E. Rasekh, D. Guo, U. Madhow, Joint routing and resource allocation for millimeter wave picocellular backhaul, IEEE Transactions on Wireless Communications 19 (2) (2019) 783–794. doi:10.1109/twc.2019.2948624.

[35] S. Zhang, B. Yin, W. Zhang, Y. Cheng, Topology aware deep learning for wireless network optimization, IEEE Transactions on Wireless Communications 21 (11) (2022) 9791–9805. doi:10.1109/twc.2022.3179352.

[46] Y. Liu, H. Mao, L. Zhu, Z. Xiao, Z. Han, X.G. Xia, Routing and resource scheduling for airground integrated mesh networks, IEEE Transactions on Wireless Communications 22 (6) (2022) 4090–4105. doi:10.1109/twc.2022.3223152.

[36] H. Wang, G. Ding, J. Wang, S. Wang, L. Wang, Power control for multiple interfering D2D communications underlaying cellular networks: An approximate interior point approach, in: 2017 IEEE International Conference on Communications Workshops (ICC Workshops), IEEE, 2017, pp. 1346–1351. doi:10.1109/iccw.2017.7962846.

[47] J. Gu, S. J. Bae, S. F. Hasan, M. Y. Chung, Heuristic algorithm for proportional fair scheduling in D2D-cellular systems, IEEE Transactions on Wireless Communications 15 (1) (2015) 769–780. doi:10.1109/twc.2015.2477998.

[37] A. Capone, J. Elias, F. Martignon, Routing and resource optimization in service overlay networks, Computer Networks 53 (2) (2009) 180–190. doi:10.1016/j.comnet.2008.09.011.

[48] K. Yang, J. Wu, X. Gao, X. Bu, S. Guo, Energy-efficient power control for device-to-device communications with max-min fairness, in: 2016 IEEE 84th vehicular technology conference (VTC-Fall), IEEE, 2016, pp. 1–5. doi:10.1109/vtcfall.2016.7880995.

[38] D. Zhai, M. Sheng, X. Wang, Z. Sun, C. Xu, J. Li, Energy-saving resource management for D2D and cellular coexisting networks enhanced by hybrid multiple access technologies, IEEE Transactions on Wireless Communications 16 (4) (2017) 2678–2692. doi:10.1109/twc.2017.2671863.

[49] C. Zhao, Y. Cai, A. Liu, M. Zhao, L. Hanzo, Mobile edge computing meets mmWave communications: Joint beamforming and resource allocation for system delay minimization, IEEE Transactions on Wireless Communications 19 (4) (2020) 2382–2396. doi:10.1109/twc.2020.2964543.

[39] L. Xu, S. Haddad Vanier, Branch-and-price for energy optimization in multi-hop wireless sensor networks, Networks 80 (1) (2022) 123–148. doi:10.1002/net.22083.

[50] S. Boyd, L. Vandenberghe, mization, Cambridge university doi:10.1017/cbo9780511804441.

[40] B. Zhang, F. Devoti, I. Filippini, D. De Donno, Resource allocation in mmWave 5G IAB networks: A reinforcement learning approach based on column generation, Computer Networks 196 (2021) 108248. doi:10.1016/j.comnet.2021.108248.

Convex press,

opti2004.

[51] Y. Nesterov, Smooth minimization of non-smooth functions, Mathematical programming 103 (1) (2005) 127– 152. doi:10.1007/s10107-004-0552-5.

[41] A. Omidkar, A. Khalili, H. H. Nguyen, H. Shafiei, Reinforcement-learning-based resource allocation for energy-harvesting-aided D2D communications in IoT networks, IEEE Internet of Things Journal 9 (17) (2022) 16521–16531. doi:10.1109/jiot.2022.3151001.

[52] S. Lu, I. Tsaknakis, M. Hong, Y. Chen, Hybrid block successive approximation for one-sided non-convex minmax problems: Algorithms and applications, IEEE Transactions on Signal Processing 68 (2020) 3676–3691. doi:10.1109/tsp.2020.2986363. 28

[53] X. Chen, Y. Cai, M. Zhao, M.-M. Zhao, Joint computation offloading and resource allocation for min-max fairness in MEC systems, in: 2019 IEEE Wireless Communications and Networking Conference (WCNC), IEEE, 2019, pp. 1–6. doi:10.1109/wcnc.2019.8885984. [54] S. Liesegang, S. Buzzi, EMF-compliant power control in cell-free massive MIMO: Model-based and data-driven approaches, IEEE Transactions on Wireless Communications 25 (2026) 12246–12262. doi:10.1109/twc.2026.3664677. [55] J. Gondzio, Interior point methods 25 years later, European Journal of Operational Research 218 (3) (2012) 587– 601. doi:10.1016/j.ejor.2011.09.017. [56] J. Gondzio, Interior point methods in the year 2025, EURO Journal on Computational Optimization 13 (2025) 100105. doi:10.1016/j.ejco.2025.100105. [57] J. Y. Yen, Finding the K shortest loopless paths in a network, management Science 17 (11) (1971) 712–716. doi:10.1287/mnsc.17.11.712. [58] S. J. Wright, Coordinate descent algorithms, Mathematical programming 151 (1) (2015) 3–34. doi:10.1007/s10107-015-0892-3. [59] S. J. Wright, Primal-dual interior-point methods, SIAM, 1997. doi:10.1137/1.9781611971453.

29