通过GPU加速的近端消息传递实现大规模网络效用最大化

大家读完觉得有帮助记得关注和点赞！！！

摘要​

我们提出了一种GPU加速的近端消息传递算法（PMP），用于解决大规模网络效用最大化（NUM）问题。NUM是资源分配中的一个基本问题，其目标是在满足链路容量约束的条件下，将资源 across 网络中的各种流进行分配以最大化总效用。我们的方法是ADMM的一种变体，仅需要与链路-路由矩阵进行稀疏矩阵-向量乘法以及逐元素的近端算子求值，从而实现了跨流和链路的完全并行更新。它支持异构效用类型（包括NUM中常见的对数效用），并且不假设严格的凹性。我们在PyTorch中实现了该方法，并在包含数千万变量和约束的问题上展示了其性能，相比现有的CPU和GPU求解器实现了4倍到20倍的加速，并且能够解决耗尽基线方法内存的问题规模。此外，我们表明我们的算法对拥塞和链路容量下降具有鲁棒性。最后，通过一个时间扩展的交通座位分配案例研究，我们说明了该方法如何在现实网络中产生可解释的分配方案。

​I 引言​

网络效用最大化（NUM）是一个用于在网络中竞争用户之间分配资源的优化框架。Kelly等人研究并形式化了大型网络（如互联网）中拥塞控制和公平资源分配之间的联系。这种网络公平资源分配的范式是通用的，也可应用于传感器网络、缓存系统和交通系统[citation:2, 3, 4]。我们专注于互联网和交通网络，这些网络通常可扩展到 105−106个变量和约束。

解决NUM的方法范围从集中式内点法

到分布式算法，其中许多基于对偶分解[citation:7, 8, 9]。这些对偶分解方法通常要求网络中的效用函数是严格凹的。本工作中考虑的对数效用允许将NUM重新表述为带有指数锥约束的锥优化问题，以供现成的锥求解器使用。最先进的商业求解器，例如MOSEK，支持具有原生指数锥处理功能的内点算法，并且非常适合解决中小型NUM实例。最近，针对大规模凸优化问题出现了GPU加速算法。CuClarabel是一种用于解决锥优化问题的GPU加速内点法，能够处理指数锥约束。然而，CuClarabel依赖CUDSS来求解线性系统，这对于非常大尺度的问题可能慢得令人无法接受。文献[12]的作者开发了原始-对偶混合梯度（PDHG）算法的一种变体，以有效解决大规模线性规划问题。一些工作利用GPU加速重启PDHG来解决线性和二次规划问题[citation:13, 14]。虽然这些PDHG方法不涉及求解线性系统，只需要矩阵-向量乘积，但先前的求解器都不支持NUM中常用的对数目标。文献[15]的作者开发了一种基于PDHG的GPU兼容算法，用于解决大规模多商品网络流问题。Kraning等人开发了近端消息传递（PMP）算法，这是交替方向乘子法（ADMM）的一种变体，用于在电力网络上解决非常大尺度的直流最优潮流问题（DC-OPF）。Degleris等人扩展了的框架，通过在GPU上实现用于DC-OPF的PMP算法，仅使用稀疏关联矩阵乘法和向量化标量操作。

在这项工作中，我们采用的GPU加速PMP算法来有效解决大规模NUM。我们的方法完全是分布式的，处理极大规模的问题实例，并且可以支持多种效用函数。

​II 问题设置​

我们考虑一个涉及 n个流量流和 m条链路的NUM问题。每个流量流 j有一个由链路子集组成的固定路由和一个效用函数 Uj​:R→R，该函数是凹的、二次可微的，并且定义域包含在 R+​中。分配流速率 xj​所产生的效用由 Uj​(xj​)给出。在这项工作中，我们考虑形式为 Uj​(xj​)=wj​xj​和 Uj​(xj​)=wj​log(xj​)的效用函数。我们注意到后者对应于加权比例公平，是中研究的更一般的 α-公平效用函数族的一个特例。我们的框架也适用于其他效用函数的选择，只要它们具有易于计算的近端算子。一个这样的例子是 α=2公平性，它对应于通信网络中的最小潜在延迟公平性

。总网络效用是 U(x)=∑j=1n​Uj​(xj​)，其中 x∈Rn是流速率向量。令 R∈Rm×n为链路-路由矩阵，定义如下：

每条链路 i∈{1,…,m}具有容量 ci​>0。链路 i上的流量不能超过 ci​，且流速率必须非负。

NUM问题是：

其中变量是 x∈Rn。我们可以将问题(2)重写为：

其中 s∈Rm是一个松弛变量，并且 x的非负约束隐含在目标函数的定义域中。当效用函数 Uj​(xj​)是凹函数时，问题(3)是一个凸最小化问题。

​III 近端消息传递​

为了推导出NUM问题的PMP算法，我们将问题(3)视为一个二分图，其中两个节点集分别是流和链路，连接它们的边是终端。令终端集合为 J={1,…,J}，流集合为 S={σ1​,…,σS​}，链路集合为 L={l1​,…,lL​}。每个终端 j∈J恰好连接一个流 σ∈S和一条链路 l∈L，因此流集合和链路集合各自划分了 J。我们将展示问题(3)的结构可以解释为这个具有流、链路和终端的二分图。二分图中的每个终端包含流和链路之间的流量。流是具有产生或消耗流量成本的节点，而链路是具有流量平衡约束的节点。

为了说明这一点，我们扩展等式约束的单行以准确理解在每个链路上求和为0的流量。对第 i行这样做，我们得到：

注意 x,s是我们问题中的优化变量。从(4)我们可以看到两种类型的流量项：Rij​xj​，代表连接到链路 i的流量流 j的流量；以及 si​−ci​，可以解释为链路 i处的“松弛”流流量。我们总是有两种形式的流：(i) 原始流量流和 (ii) 松弛流。松弛流允许我们人为地使链路的容量饱和，从而保持流量平衡。

我们定义变量 p∈R∣J∣来表示系统中的终端流量。我们写 j∈σ和 j∈l来表示终端 j连接到流 σ或链路 l。令 ∣σ∣和 ∣l∣表示连接到给定流或链路的终端数量。我们使用集合值索引，例如 pl​∈R∣l∣和 pσ​∈R∣σ∣，分别表示连接到链路 l和流 σ的所有终端的流量。

​III-A 一个小例子​

考虑一个具有三个流量流和三个链路的网络，其等式约束集如下：

图1显示了相应的二分图。

​图1: NUM的二分图​

​III-B 流​

每个流（无论是流量流还是松弛流）都有一个效用函数。我们为各种流定义成本函数。

​对数效用流​

这里，xσ​作为流量流的局部变量，控制终端流量。约束 pσ​=1∣σ∣​xσ​意味着对于一个具有多个终端的流量流，每个关联终端中的流量必须恰好等于 xσ​。标量 wσ​>0是流权重。

​线性效用流​

​松弛流​

该成本是强制执行 pσ​+cσ​≥0的指示函数。

​其他效用流​

我们的方法也适应更复杂的效用函数。我们只要求这些函数是凹的、二次可微的效用函数，并且具有简单的近端算子（参见第III-E节）。

​III-C 链路​

一条链路上所有流量的总和必须平衡到0。为了编码这一点，我们为每条链路 l∈L引入一个指示函数：

​III-D 来自ADMM的近端消息传递​

问题(3)可以重新表述为：

我们如下解释问题(9)：n个效用流在其路由链路上产生流量，而 m个松弛流注入所需的剩余流量以使每条链路的容量饱和。第一个目标项编码流成本，而第二个通过链路流量平衡约束强制执行可行性。

正如ADMM中的典型做法，我们引入一个复制变量 z。变量 p代表终端流量变量的流端副本，而 z是流量变量的链路端副本。我们在ADMM中的子问题跨流和链路是可分离的。遵循[citation:16, 17, 19]，我们可以将问题(9)的缩放增广拉格朗日函数写为：

其中 ρ>0是惩罚参数，u是缩放对偶变量，由 y/ρ给出，其中 y是未缩放的对偶变量。使用(10)中的增广拉格朗日函数，我们可以推导出ADMM迭代如下：

遵循，我们进一步将上述ADMM迭代简化为近端消息传递。

1. ​近端设备更新​

2. ​缩放价格更新​

   其中函数 f的近端算子由下式给出：

   并且 pˉ​σ​∈R∣σ∣和 pˉ​l​∈R∣l∣是向量，其中每个值分别是连接到流 σ或链路 l的所有终端流量的平均值。

​III-E 近端更新​

(11)中近端更新的有效性取决于我们计算系统中各种流的近端算子的效率。具有对数或线性效用的流量流以及松弛流都具有简单的近端算子，在点 z处求值：

1. ​对数效用流：​​

2. ​线性效用流：​​

3. ​松弛流：​​

​III-F 收敛性​

第 i次迭代的原始和对偶残差分别为。原始残差是所有链路上的净流量不平衡，这正好是问题(9)中原始不可行性的度量。对偶残差是流量与每个网络上平均流量之间连续迭代的差。我们根据以下标准终止PMP算法：

其中 εtol​=εabs​∣J∣ ​且 εabs​>0是绝对容差。因为PMP是ADMM的一种变体，关于ADMM收敛的更多细节我们请读者参考。

​III-G 加速收敛​

我们使用几种技术来加速PMP收敛。

​残差平衡​

我们自适应地调整惩罚参数 ρ，以平衡原始和对偶残差。我们利用中的简单更新方案：

其中 γ>1和 μ>1是更新规则的参数。

每次更新 ρ时，我们相应地重新缩放 u：

我们设置 μ=2和 γ=1.1，并且每50次迭代更新一次 ρ。

​过松弛​

过松弛通常用于提高ADMM的收敛速度。令 α∈[1,2]表示松弛参数。在 z,u更新中，我们将 plk+1​替换为：

使用过松弛时，我们的更新不再简化为上面推导的消息传递更新。相反，我们写：

在我们的实验中，我们发现 α=1.6效果很好。

​图2：​​ 不同问题规模上求解器的扩展性能

​图3：​​ 无拥塞与有拥塞的收敛情况

​IV GPU实现​

PMP中两个需要高效GPU实现的关键操作是：(i) 每个流的近端更新和 (ii) 计算跨链路的平均值。

​IV-A 流数据模型​

流按效用函数和终端数量分组为类型 c。这允许向量化的近端更新，因为相同类型的流共享相同的效用函数和流量维度。对于每种类型 c，我们将其变量存储在维度为 ∣c∣的 τc​张量中，其中 ∣c∣是该类型的流数量。每个这样的张量与类型 c设备的一个终端相关联。

​IV-B 向量化​

流的近端更新是闭式的，这意味着我们可以使用完全向量化、并行更新来计算它们。对于每种流类型 c，我们将其 ∣c∣个流进行批处理并同时更新它们。例如，对于一批 ∣c∣=k个流，在(15)中给出的线性效用更新是：

其中 Jτc​,k​∈Rτc​×k是全一矩阵，Z=[z1​,…,zk​]∈Rτc​×k是堆叠的评估点（每流一列），而 wτc​​∈Rk是流权重。方程(22)通过广播允许高效的GPU实现。

​IV-C 平均值​

PMP更新需要有效计算连接到每条链路的所有终端的平均值。这些平均值对于连接到同一链路的所有终端是相同的，因此我们将它们存储为大小为 m的张量。从概念上讲，平均是跨链路收集信息，可以表示为乘以稀疏关联矩阵，这可以使用散射核在GPU上实现。

形式上，对于流类型 c，令 Rci​∈Rm×∣c∣表示终端 i的关联矩阵，其元素为：

其中 i=1,…,τc​。矩阵 Rci​是每个流类型的每个终端的关联矩阵。平均量 pˉ​然后由下式给出：

其中 pci​是一个 ∣c∣维张量，pˉ​是一个 m维张量，∣l∣∈Rm是每条链路的终端计数向量，而 ⊙表示逐元素乘法。

​IV-D 软件实现​

我们为NUM问题开发了一个开源的PyTorch PMP算法实现。我们的实现是模块化的，支持线性和对数效用，并且可以通过指定所需的效用函数及其相应的近端算子扩展到其他流类型。在配备单个具有80GB内存的NVIDIA A100-SXM4 GPU的机器上，我们可以在1847秒内解决具有 m=1×107条链路和 n=5×106个流的问题。

​图4：(a) 无拥塞路由长度。(b) 无拥塞链路利用率。(c) 有拥塞链路利用率。​​

​图5：链路降级和故障下的热启动​

​V 性能结果​

我们进行了几个数值实验，以证明我们的方法在各种不同问题规模和场景下的有效性。所有实验均在单个具有80GB内存的NVIDIA A100-SXM4 GPU上进行，该GPU由32个虚拟CPU核心和64GB RAM支持。虽然我们的实验是合成构建的，但它们映射到互联网拥塞控制和交通座位分配的应用。我们首先将我们的方法与几种开源和商业求解器进行基准测试，以证明我们的方法在无拥塞网络上具有卓越的扩展性能。我们检查了PMP在更复杂的具有拥塞链路的网络中的性能，表明我们的方法是鲁棒的，并且仍然能够收敛到可接受的解。最后，我们观察到热启动我们的方法允许它在网络降级的情况下快速重新计算最优分配。

​V-A 在互联网和交通设置中的应用​

本节中的数值实验是合成的，但直接适用于两个实际领域。在互联网拥塞控制中，流是端到端的会话，链路是带宽受限的通信链路，链路降级模拟中断。目标是在带宽约束下最大化总效用。在交通中，流是行程-出发选项，链路是时间扩展的车辆-座位资源；当行程重叠时出现拥塞，降级反映了服务中断。

​V-B 扩展性​

我们按照[6, 9]的方法随机生成链路-路由矩阵 R。给定 m条总链路和 n=m/2个流量流，我们为流分配链路，使得平均流由10条链路组成。我们将我们的方法与商业求解器MOSEK以及开源锥求解器Clarabel和CuClarabel进行基准测试，后者为GPU加速而设计。所有求解器均在低精度设置下运行，除非另有说明。图3显示了挂钟时间与问题规模（即链路数 m）的关系。如果求解器未能在10小时内返回解，我们将其标记为未完成（DNF）；如果由于内存限制而失败，则标记为内存不足（OOM）。对于小问题规模，MOSEK和CuClarabel优于我们的方法，但在 m=2×105时，我们的方法比CuClarabel快约4倍，比MOSEK快约20倍。对于 m≥5×105，只有我们的方法能够解决实例而不会耗尽内存。我们还注意到，低精度下的MOSEK在三个实例（m=3×104,m=5×104,m=105）上失败，返回NaN值。

​V-C 拥塞网络​

我们通过使一小部分链路 heavily utilized 来模拟拥塞。我们随机选择0.1%的链路，并将每条链路连接到大约10%的流。在图3中，我们看到PMP在大约1000次迭代中达到中等精度（原始和对偶残差 <10−4）（无拥塞情况）。拥塞降低了 R的稀疏性，并略微减慢收敛速度。然而，只需额外几百次迭代即可达到类似的精度。

在图4a中，我们展示了一个无拥塞网络（m=1×106且 n=5×105）的每条流的链路数分布。在图4b–4c中，我们 then 绘制了相同维度下无拥塞和拥塞网络的每条链路的流数分布。相对于无拥塞情况，拥塞在每条链路的流数分布中引入了重尾。这些 heavily used 的链路对应于 R中的密集行，并且是收敛速度适度减慢的原因。

​V-D 链路降级下的重新分配​

接下来，我们考虑在链路容量降级下的最优分配。从链路-路由矩阵 R和容量 c的解 x0⋆​开始，我们在降级容量 cd​下重新求解，其中每条链路独立地有25%的概率发生50%的容量减少。在图5中，我们看到从 x0⋆​热启动产生大约5倍的加速。

在极端情况下，链路完全故障。鉴于NUM的结构，故障链路会移除其 R中的行，并消除所有经过它的流。然后，在以25%的概率故障每条链路后，我们用修剪后的 Rpruned​,cpruned​解决修剪后的实例。在图5中，我们看到结合修剪和热启动在求解时间上产生大约10倍的加速。

​VI 案例研究：交通座位分配​

​VI-A 问题设置​

我们考虑一个具有多个行程-出发选项的行程-预订设置[4]。我们将交通系统建模为一个时间扩展网络，其中包含站点 V、有向空间边 E⊆V×V和离散时间仓 t=0,…,T−1，长度 Δt。链路是时间戳边，L=E×{0,…,T−1}，其中 ℓ=(e,t)。所有链路的旅行时间为 τℓ​=Δt.

一个行程-出发流是 j=(k,r,t0​)，其中 (ok​,dk​)是一个OD对，r∈Rk​是一条空间路由（Rk​是OD对 k的所有空间路由的集合），而 t0​是一个出发仓。每个流接收分配 xj​>0和效用 Uj​(xj​)=wj​logxj​，其中 wj​编码行程偏好（例如，非高峰旅行）。链路代表预定的车辆-座位容量，流代表乘客选项。我们解决连续松弛问题并根据需要进行舍入。

​​(a)​​

​​(b)​​

​图6：​​ 在 t0​=0出发的两条路由的流量和链路价格。

​VI-B 实例生成与求解​

我们构建了一个合成的时间扩展网络，其中 S=100个站点，时间步长 Δt=5分钟，T=192。这产生 m=∣E∣T=182,784条链路和 n=19,800个流，具有加权对数效用（wj​=1）。为了解释解的一个切片，我们比较同一OD和出发时间的两条行程。图6(a)显示了空间路由。线宽与分配的流量成正比。这里，蓝色路由接收了近5倍的流量。

这种差异源于链路影子价格 λ(ℓ,t)≥0。我们将流 j的路径价格定义为：

对于加权对数效用，KKT平稳性给出 wj​/xj​=πj​，因此价格较低的路径接收更多流量。在图6(b)中，每条路由上的归一化价格显示橙色路径上有一个额外的瓶颈，抬高了 πj​，而蓝色路径只有相同的初始高价格链路。NUM将价格集中在拥塞链路上，并将分配转移到更便宜的路由。

​VII 结论​

我们提出了一种GPU加速的近端消息传递求解器，用于解决大规模网络效用最大化（NUM）问题。我们的方法仅使用与链路-路由矩阵的稀疏矩阵-向量乘法以及闭式近端更新。我们在纯PyTorch中实现了该方法，并表明它可以扩展到数百万个变量，比强大的CPU/GPU基线提供4倍到20倍的加速，并且可以解决那些方法无法解决的实例。然后我们证明我们的方法在拥塞和链路容量降级下是鲁棒的，并且在扰动后从热启动中显著受益。最后，我们使用一个交通座位分配案例研究来演示如何将NUM应用于现实的操作问题，在重叠的行程选项中产生可解释的分配，并强调我们方法更广泛的适用性。