Lecture 3： Resource allocation as utility maximization 笔记
本节课主要介绍了使用primal算法与dual算法最大化效用函数(utility function)的问题。
1 前置知识
1.1 符号约定
- V ˙ \dot{V} V˙ 表示函数 V V V的导数。
- ∇ a \nabla a ∇a 表示向量 a a a的梯度。
- f ( x ) f(x) f(x) 表示关于自变量 x x x的连续可导函数。
- f ( x ∗ ) f(x^*) f(x∗) 表示函数 f ( x ) f(x) f(x)的全局最大/小值，其中 x ∗ x^* x∗为使函数 f ( x ) f(x) f(x)取到全局最大/小值的向量。
- x ( t ) x(t) x(t)为 t t t时刻的数据传输速率，在不考虑时间的时候写作 x x x。
- x ˙ ( t ) \dot{x}(t) x˙(t)为 x ( t ) x(t) x(t)关于 t t t的导数，在不考虑时间的时候写作 x ˙ \dot{x} x˙。
1.2 凸优化知识回顾
1.2.1 凸函数与凹函数
凸函数：设函数
f
(
x
)
f(x)
f(x)的定义域为
R
n
R^n
Rn，则对于任意
x
1
,
x
2
∈
R
n
,
λ
∈
(
0
,
1
)
x_1,x_2\in R^n, \lambda \in(0,1)
x1,x2∈Rn,λ∈(0,1)，恒有不等式
f
(
λ
x
1
+
(
1
−
λ
)
x
2
)
≤
λ
f
(
x
1
)
+
(
1
−
λ
)
f
(
x
2
)
f(\lambda x_1 + (1-\lambda )x_2)\le\lambda f(x_1)+(1-\lambda)f(x_2)
f(λx1+(1−λ)x2)≤λf(x1)+(1−λ)f(x2)成立，那么称
f
(
x
)
f(x)
f(x)为定义域上的凸函数。
∇ 2 f ( x ) ≥ 0 \nabla^2 f(x)\ge0 ∇2f(x)≥0则函数为凸函数
凸函数在定义域上有全局最小值。
常见凸函数：
- f ( x ) = x a , a ∈ [ 1 , ∞ ) f(x) = x^a, a\in[1,\infty) f(x)=xa,a∈[1,∞)
- f ( x ) = e x f(x)= e^x f(x)=ex
- f ( x ) = x l o g x f(x) = xlogx f(x)=xlogx
- f ( x ) = ∥ x ∥ p f(x) = \|x\|_p f(x)=∥x∥p，任意范数均为凸函数
- f ( x ) = max { x 1 , x 2 , ⋅ ⋅ ⋅ , x n } f(x) = \max\{x_1,x_2,\cdot\cdot\cdot,x_n\} f(x)=max{x1,x2,⋅⋅⋅,xn}
凹函数：设函数
f
(
x
)
f(x)
f(x)的定义域为
R
n
R^n
Rn，则对于任意
x
1
,
x
2
∈
R
n
,
λ
∈
(
0
,
1
)
x_1,x_2\in R^n, \lambda \in(0,1)
x1,x2∈Rn,λ∈(0,1)，恒有不等式
f
(
λ
x
1
+
(
1
−
λ
)
x
2
)
≥
λ
f
(
x
1
)
+
(
1
−
λ
)
f
(
x
2
)
f(\lambda x_1 + (1-\lambda )x_2)\ge\lambda f(x_1)+(1-\lambda)f(x_2)
f(λx1+(1−λ)x2)≥λf(x1)+(1−λ)f(x2)成立，那么称
f
(
x
)
f(x)
f(x)为定义域上的凹函数。
∇
2
f
(
x
)
≤
0
\nabla^2 f(x)\le0
∇2f(x)≤0则函数为凹函数
对任意凸函数
f
(
x
)
f(x)
f(x)取相反数，
g
(
x
)
=
−
f
(
x
)
g(x)=-f(x)
g(x)=−f(x)一定为凹函数。
凹函数在定义域上有全局最大值。
1.2.2 最速下降算法
本文所提出的所有算法均假设目标函数为光滑凸函数(凹函数取负号)，不考虑非光滑或非凸优化。
考虑凸函数 f ( x ) f(x) f(x)，其Hessian矩阵半正定( ∇ 2 f ( x ) ≥ 0 \nabla^2 f(x)\ge0 ∇2f(x)≥0)，所以其当 ∇ f ( x ) ≥ 0 \nabla f(x)\ge0 ∇f(x)≥0， ∇ f ( x ) \nabla f(x) ∇f(x)递增， ∇ f ( x ) ≤ 0 \nabla f(x)\le0 ∇f(x)≤0时， ∇ f ( x ) \nabla f(x) ∇f(x)递减。考虑梯度的数学定义：某一函数在该点处的方向导数沿梯度方向取得最大值，即函数在该点处沿着梯度方向变化最快。在凸函数中要最小化凸函数，应沿梯度的反方向寻找最小点。同理对于凹函数，则应按梯度方向寻找函数的最大值点。
因此我们可以得到最速下降算法的基本过程：
- 初始化 x k x_k xk
- 计算 ∇ f ( x k ) \nabla f(x_k) ∇f(xk)
- 计算 a k = a r g m i n a f ( x k + a ∇ f ( x k ) ) a_k = argmin_a f(x_k+a\nabla f(x_k)) ak=argminaf(xk+a∇f(xk))
- x k + 1 = x k + a k ∇ f ( x k ) x_{k+1} = x_k+a_k\nabla f(x_k) xk+1=xk+ak∇f(xk)
- 如果 ∇ f ( x k + 1 ) = 0 \nabla f(x_{k+1})=0 ∇f(xk+1)=0， x ∗ = x k + 1 x^*=x_{k+1} x∗=xk+1，输出 x k + 1 x_{k+1} xk+1.否则 k = k + 1 k=k+1 k=k+1，返回到第二步。
1.2.3 梯度线搜索算法
最速下降法保证了每一步都以最大的速度降低函数值，但是无法保证收敛性。
步子迈太大了，咔！容易扯着蛋。——《让子弹飞》
而且如何获取满足步骤3的 a a a，本身也是一项计算量很大的工作，甚至在某些情况下，求解 a a a的计算量不亚于求解 min f ( x ) \min f(x) minf(x)。所以我们需要改进步骤3。
常用改进算法有两种，分别是Wolfe Conditions和Goldstein Conditions。这里只介绍Wolfe Conditions。其对 a k a_k ak的要求如下：
- f ( x k + a p k ) ≤ f ( x k ) + c 1 a ∇ f ( x k ) T p k f(x_k+ap_k)\le f(x_k)+c_1a\nabla f(x_k)^Tp_k f(xk+apk)≤f(xk)+c1a∇f(xk)Tpk
- f ( x k + a p k ) T p k ≥ f ( x k ) + c 2 ∇ f ( x k ) T p k f(x_k+ap_k)^Tp_k\ge f(x_k)+c_2\nabla f(x_k)^Tp_k f(xk+apk)Tpk≥f(xk)+c2∇f(xk)Tpk
- 0 < c 1 < c 2 < 1 0<c_1<c_2<1 0<c1<c2<1
其中 p k p_k pk为下降搜索方向，当采用一维方法时为梯度方向 − ∇ f ( x k ) -\nabla f(x_k) −∇f(xk)，采用牛顿算法时为 − H k ∇ f ( x k ) -H_k\nabla f(x_k) −Hk∇f(xk)，其中 H k ∇ 2 f ( x k ) = I H_k\nabla^2 f(x_k)=I Hk∇2f(xk)=I.
首先看 ϕ ( a ) f ( x k ) + c 1 a ∇ f ( x k ) T p k \phi(a) f(x_k)+c_1a\nabla f(x_k)^Tp_k ϕ(a)f(xk)+c1a∇f(xk)Tpk的含义。由于除了 a a a外所有的符号都是常数值，显然这是一条类似 y = a x + b y=ax+b y=ax+b定义方式的直线。其截距为 f ( x k ) f(x_k) f(xk)，斜率是 c 1 ∇ f ( x k ) T p k c_1\nabla f(x_k)^Tp_k c1∇f(xk)Tpk。这保证了 x k + 1 x_{k+1} xk+1起码要比直线下降的值要小，即取值点要在图中最长的那条射线的下面。
再来看第二个条件。很显然的一个问题是，在 a a a趋近于 0 0 0的时候，第一个条件显然是满足的。但是优化的步长也就太小了，那么我们希望步子稍微迈的大一点。那么什么时候步子迈的够长了呢，一个显而易见的标准就是：当增加 a a a不能明显降低函数值时，步子就够长了。也就是关于 a a a的一阶导足够小的时候就可以了。有多小呢，wolfe给出的条件是小于 c 2 ∇ f ( x k ) T c_2\nabla f(x_k)^T c2∇f(xk)T就可以了，也就是在 x k + a p k x_k+ap_k xk+apk处的导数，比 c 2 c_2 c2倍的 x k x_k xk处导数小就可以了。改进后的算法如下
- 初始化 x k x_k xk
- 计算 ∇ f ( x k ) \nabla f(x_k) ∇f(xk)
- 计算满足Wolfe Conditions的步长 a k a_k ak
- x k + 1 = x k + a k ∇ f ( x k ) x_{k+1} = x_k+a_k\nabla f(x_k) xk+1=xk+ak∇f(xk)
- 如果 ∇ f ( x k + 1 ) = 0 \nabla f(x_{k+1})=0 ∇f(xk+1)=0， x ∗ = x k + 1 x^*=x_{k+1} x∗=xk+1，输出 x k + 1 x_{k+1} xk+1.否则 k = k + 1 k=k+1 k=k+1，返回到第二步。
1.2.4 KKT条件
KKT条件是一个必要条件，即对带约束优化问题：
min
f
(
x
)
s
.
t
.
f
i
(
x
)
≤
0
,
i
∈
I
f
u
(
x
)
=
0
,
u
∈
U
\begin{aligned} &\min f(x) \\ &s.t. \,\,\, f_i(x)\le0,i\in I\\ &\quad\,\,\,\,\, f_u(x)=0,u\in U \end{aligned}
minf(x)s.t.fi(x)≤0,i∈Ifu(x)=0,u∈U
那么其一定要满足：
∇
x
L
(
x
)
=
∇
x
f
(
x
)
+
∑
i
∈
I
λ
i
∇
x
f
i
(
x
)
+
∑
u
∈
U
μ
u
∇
x
f
u
(
x
)
=
0
∀
μ
≠
0
∀
λ
≥
0
∀
λ
i
f
i
(
x
)
=
0
∀
f
i
(
x
)
≤
0
∀
f
u
(
x
)
=
0
\begin{aligned} \nabla_x L(x) = &\nabla_x f(x) + \sum_{i\in I}\lambda_i \nabla_xf_i(x)+\sum_{u\in U}\mu_u \nabla_xf_u(x) = 0\\ &\forall\mu \not=0\\ &\forall\lambda_ \ge 0\\ &\forall\lambda_i f_i(x)=0\\ &\forall f_i(x)\le0\\ &\forall f_u(x)=0\\ \end{aligned}
∇xL(x)=∇xf(x)+i∈I∑λi∇xfi(x)+u∈U∑μu∇xfu(x)=0∀μ=0∀λ≥0∀λifi(x)=0∀fi(x)≤0∀fu(x)=0
λ
\lambda
λ和
μ
\mu
μ都是拉格朗日算子。KKT只是必要条件，不是充分条件，但是一般满足KKT条件的点，都是最优点。
1.3 Lyapunov(李雅普诺)定理
1.3.1 Lyapunov边界定理
对于定义在 R n R^n Rn上的函数 V ( x ) V(x) V(x)，如果当 ∥ x ∥ → ∞ \|x\|\rightarrow\infin ∥x∥→∞时， V ( x ) → ∞ V(x)\rightarrow\infin V(x)→∞。那么 V ( x ) V(x) V(x)的导数 V ˙ = ∇ V ( x ) x ˙ = ∇ V ∇ f ( x ) ≤ 0 \dot{V}=\nabla V(x)\dot{x}=\nabla V\nabla f(x)\le0 V˙=∇V(x)x˙=∇V∇f(x)≤0，对所有的 x x x成立，那么存在一个常数 B B B，满足对所有 t t t恒有 ∥ x ( t ) ∥ ≤ B \|x(t)\|\le B ∥x(t)∥≤B。
1.3.2 Lyapunov全局渐进稳定定理
当函数
V
(
x
)
V(x)
V(x)连续可导时，若(1)
V
(
x
)
≥
0
V(x)\ge0
V(x)≥0对任何
x
x
x成立，且只在
x
=
0
x=0
x=0时，
V
(
x
)
=
0
V(x)=0
V(x)=0，(2)
V
˙
<
0
\dot{V}<0
V˙<0对任何
x
≠
0
x\not=0
x=0成立，且
V
˙
(
0
)
=
0
\dot{V}(0)=0
V˙(0)=0。
那么平衡点
x
e
=
0
x_e=0
xe=0是全局渐进稳定的。
2 Primal algorithm
2.1 Primal solution
首先需要明确使用场景。这里我们认为数据传输发生在源节点与目的节点之间。源节点和目的节点之间有许多路由器，用来转发数据。所有的路由器组成了集合
S
\mathcal{S}
S。假设每个路由器在传输数据的过程中可以在所有链路中，选择其需要一个或多个的链路传输数据，所有链路又构成了集合
L
\mathcal{L}
L。我们的目的是要提高路由器的利用率。
显而易见的一点是，当路由器单位时间转发的数据较少时，路由器更多的资源处于空闲状态，那么利用率便会降低。所以增加单位时间的数据量(传输速率
x
x
x)可以有效提高路由器的利用率。然而提高传输速率，又会提高传输链路的负载，进而提高使用成本。
所以我们希望能在提高路由器利用效率的同时，降低链路的使用成本，即：
max
x
W
(
x
)
=
∑
r
∈
S
U
r
(
x
)
−
∑
l
∈
L
B
i
(
∑
s
:
l
∈
s
x
s
)
\max_x W(x)=\sum_{r\in \mathcal{S}}U_r(x)-\sum_{l\in\mathcal{L}}B_i(\sum_{s:l\in s}x_s)
xmaxW(x)=r∈S∑Ur(x)−l∈L∑Bi(s:l∈s∑xs)
这里我们忽略
U
r
(
x
)
U_r(x)
Ur(x)和
B
i
(
x
)
B_i(x)
Bi(x)的显示表达式，只关心他们的凸性。这里
U
r
(
x
)
U_r(x)
Ur(x)是一个凹函数，其表示路由器的利用率关于传输速率的函数，
B
i
(
x
)
B_i(x)
Bi(x)是一个凸函数，其为传输速率的代价函数。这样是为了保证
W
(
x
)
W(x)
W(x)为一凹函数。
由于
W
(
x
)
W(x)
W(x)为凹函数，我们可以使用最简单有效的线搜索算法来获取函数的最优解。那么回顾1.2.3节中的算法。我们首先需要计算
W
(
x
)
W(x)
W(x)的导数，
∇
W
(
x
)
=
U
r
′
(
x
r
)
−
B
l
′
(
∑
s
:
l
∈
s
x
s
)
=
U
r
′
(
x
r
)
−
∑
l
:
l
∈
r
f
l
(
∑
s
:
l
∈
s
x
s
)
\nabla W(x)= U^{\prime}_r(x_r)-B_l^{\prime}(\sum_{s:l\in s}x_s)=U^{\prime}_r(x_r)-\sum_{l:l\in r}f_l(\sum_{s:l\in s}x_s)
∇W(x)=Ur′(xr)−Bl′(∑s:l∈sxs)=Ur′(xr)−∑l:l∈rfl(∑s:l∈sxs)。这里我们称连续递增函数
f
l
f_l
fl为链路
l
l
l的价格函数，反应单位传输速率链路
l
l
l的使用价格。即
B
l
(
∑
s
:
l
∈
s
x
s
)
=
∫
0
∑
s
:
l
∈
s
x
s
f
l
(
y
)
d
y
B_l(\sum_{s:l\in s}x_s)=\int_0^{\sum_{s:l\in s}x_s}f_l(y)dy
Bl(∑s:l∈sxs)=∫0∑s:l∈sxsfl(y)dy.
这样我们就可以得到传输速率的更新方程：
Primal Solution:
x
˙
=
k
r
(
x
r
)
(
U
r
′
(
x
r
)
−
∑
l
:
l
∈
r
f
l
(
∑
s
:
l
∈
s
x
s
)
)
\text{Primal Solution:}\quad \dot{x}=k_r(x_r)(U^{\prime}_r(x_r)-\sum_{l:l\in r}f_l(\sum_{s:l\in s}x_s))
Primal Solution:x˙=kr(xr)(Ur′(xr)−l:l∈r∑fl(s:l∈s∑xs))
其中
k
r
(
x
r
)
k_r(x_r)
kr(xr)是步长函数，其可以根据wolfe条件搜索得到。
这一个更新方程看似平常，但是笔者一开始也没想明白，为什么用的是 x x x关于时间的导数。其实回看梯度下降算法的第四步我们可以得到 Δ x = x k + 1 − x k = a ∇ f ( x ) \Delta x=x_{k+1}-x_k=a\nabla f(x) Δx=xk+1−xk=a∇f(x)。但是由于这里的 x x x是关于时间 t t t连续变化的，所以 Δ x = x k + 1 − x k = ∫ k k + 1 x ˙ \Delta x=x_{k+1}-x_k=\int_k^{k+1}\dot{x} Δx=xk+1−xk=∫kk+1x˙。所以更新方程的左侧是 x ˙ \dot{x} x˙。
2.2 Primal solution的渐进稳定性证明
设 V ( x ) = W ( x ∗ ) − W ( x ) V(x)=W(x^*)-W(x) V(x)=W(x∗)−W(x)。显然 W ( x ∗ ) W(x^*) W(x∗)是最高效用值，所以任何其他 W ( x ) W(x) W(x)都满足 W ( x ) ≤ W ( x ∗ ) W(x)\le W(x^*) W(x)≤W(x∗)。也就是说 V ( x ) ≥ 0 V(x)\ge0 V(x)≥0。并且可以计算得到仅在 x = x ∗ x=x^* x=x∗时 V ˙ = 0 \dot{V}=0 V˙=0。因此可以得到 x ∗ x^* x∗是全局渐进稳定点，也就是说不管传输速率初始值是多少，只要用primal算法对其进行优化，随着时间的流逝，传输速率都会到全局最优点。
构造 V ( x ) = W ( x ∗ ) − W ( x ) V(x)=W(x^*)-W(x) V(x)=W(x∗)−W(x)是个很有用的技巧。
3 Duel algorithm
Primal algorithm本质上是罚函数法，虽然很高效，但是构造合适的能使目标函数收敛到全局最优解的罚函数是很困难的。所以我们还是回到最基本的思路，求解带约束的优化方程。
max
x
∑
r
∈
S
U
r
(
x
r
)
s
.
t
.
∑
r
:
l
∈
r
x
r
≤
c
l
,
∀
l
∈
L
x
r
≥
0
,
∀
r
∈
S
\begin{aligned} &\max_x \sum_{r\in S}U_r(x_r)\\ &s.t. \sum_{r:l\in r}x_r\le c_l, \forall l\in \mathcal{L}\\ &\quad\,\,\,\, x_r\ge 0, \forall r\in \mathcal{S} \end{aligned}
xmaxr∈S∑Ur(xr)s.t.r:l∈r∑xr≤cl,∀l∈Lxr≥0,∀r∈S
这里
c
l
c_l
cl使信道容量。剩下的工作没有什么特别的，就是构造拉格朗在对偶函数
D
(
p
)
=
max
x
∑
x
r
≥
0
U
r
(
x
r
)
−
∑
l
p
l
(
∑
s
:
l
∈
s
x
s
−
c
l
)
D(p)=\max_x \sum_{x_r \ge 0}U_r(x_r)-\sum_lp_l(\sum_{s:l\in s}x_s-c_l)
D(p)=maxx∑xr≥0Ur(xr)−∑lpl(∑s:l∈sxs−cl)，其中
p
l
p_l
pl是拉格朗日算子，然后求解KKT。
需要说明的几点是
- 对拉格朗日对偶函数求导，并令导数为0，利用一阶最优条件可以得到 x r x_r xr需要满足 U r ′ = ∑ s : l ∈ s p l U_r^{\prime}=\sum_{s:l\in s}p_l Ur′=∑s:l∈spl
- 在更新过程中，对 x r x_r xr的更新由步骤1获取， x r = U r ′ , − 1 ( ∑ s : l ∈ s p l ) x_r=U_r^{\prime,-1}(\sum_{s:l\in s}p_l) xr=Ur′,−1(∑s:l∈spl)
- 对
D
(
p
)
D(p)
D(p)关于
p
p
p求导，并利用线搜索算法可以得到
p
l
p_l
pl的更新方程：
p
l
=
a
(
∑
s
:
l
∈
s
x
s
−
c
l
)
p
l
+
p_l=a(\sum_{s:l\in s}x_s-c_l)^+_{p_l}
pl=a(∑s:l∈sxs−cl)pl+。这里:
( g ( x ) ) y + = g ( x ) y > 0 = m a x ( g ( x ) , 0 ) y = 0 \begin{aligned} (g(x))^+_y&=g(x)\quad\quad\quad\quad\quad y>0\\ &=max(g(x),0)\,\,\,\,\,\,\,\, y=0 \end{aligned} (g(x))y+=g(x)y>0=max(g(x),0)y=0
为什么要用这么复杂的一个函数呢？主要还是因为KKT条件中要求拉格朗日算子非负，所以加这么一个约束。所以duel solution如下：
Duel Solution: x r = U r ′ , − 1 ( ∑ s : l ∈ s p l ) p l = a ( ∑ s : l ∈ s x s − c l ) p l + \begin{aligned} \text{Duel Solution: }\quad &x_r=U_r^{\prime,-1}(\sum_{s:l\in s}p_l)\\ &p_l=a(\sum_{s:l\in s}x_s-c_l)^+_{p_l} \end{aligned} Duel Solution: xr=Ur′,−1(s:l∈s∑pl)pl=a(s:l∈s∑xs−cl)pl+
4 计算Primal算法中罚函数的方法
4.1 16bit长header法
最直观的一个算法是这样的，既然我需要计算链路的使用代价，那么我就发一个的数据包试试。从源节点发到目的节点，再从目的节点发送ACK返回。ACK通过链路的时候就会把链路的代价信息收集到，并返回给发送节点。这样所有节点也就都知道链路价格了。
4.2 1bit长header法
长header法
最直观的一个算法是这样的，既然我需要计算链路的使用代价，那么我就发一个的数据包试试。从源节点发到目的节点，再从目的节点发送ACK返回。ACK通过链路的时候就会把链路的代价信息收集到，并返回给发送节点。这样所有节点也就都知道链路价格了。
4.2 1bit长header法
16bit的ACK确实很好，但是太长了。我们可以减小称1bit的ACK。初始时，每个ACK都是0，每条链路的代价分别为 p 1 , p 2 , ⋅ ⋅ ⋅ , p n p_1,p_2,\cdot\cdot\cdot,p_n p1,p2,⋅⋅⋅,pn。每经过一个链路，就以 p l p_l pl的概率把ACK置1。这样再多个数据发送后，就可以预估处 p l p_l pl。ACK置1的概率为 P ( 1 ) = 1 − P ( 0 ) = 1 − e − ∑ l p l = p ^ P(1)=1-P(0)=1-e^{-\sum_lp_l}=\hat{p} P(1)=1−P(0)=1−e−∑lpl=p^。 ∑ l p l = − l o g ( 1 − p ^ ) \sum_lp_l=-log(1-\hat{p}) ∑lpl=−log(1−p^).