倍增法

倍增法是一种用于求解凸优化问题的算法。假设我们有一个形如

\[\begin{split}\begin{array}{ll} \mbox{最小化} & f(x)\\ \mbox{约束} & Ax = b, \end{array}\end{split}\]

的问题，其中 \(f\) 是凸函数，\(x \in \mathbf{R}^n\) 是优化变量，\(A \in \mathbf{R}^{m \times n}\) 和 \(b \in \mathbf{R}^m\) 是问题数据。

为了应用倍增法，我们首先构造增广Lagrangian函数

\[L_{\rho}(x,y) = f(x) + y^T(Ax - b) + (\rho/2)\|Ax-b\|^2_2.\]

与增广Lagrangian函数相关的对偶函数是 \(g_{\rho}(y) = \inf_x L_{\rho}(x,y)\)。对偶函数 \(g_{\rho}(y)\) 是凹函数，其最大值与原始问题的最优值相同。

我们使用梯度上升法来最大化对偶函数。梯度上升的每一步都可以归结为 \(x\) 和 \(y\) 的更新

\[\begin{split}\begin{array}{lll} x^{k+1} & := & \mathop{\rm argmin}_{x}\left(f(x) + (y^k)^T(Ax - b) + (\rho/2)\left\|Ax-b\right\|^2_2 \right) \\ y^{k+1} & := & y^{k} + \rho(Ax^{k+1}-b) \end{array}\end{split}\]

下面的CVXPY脚本实现了乘子法并使用它来解决一个优化问题。

import cvxpy as cp
import numpy as np
np.random.seed(1)

# 初始化数据。
MAX_ITERS = 10
rho = 1.0
n = 20
m = 10
A = np.random.randn(m,n)
b = np.random.randn(m)

# 初始化问题。
x = cp.Variable(shape=n)
f = cp.norm(x, 1)

# 使用CVXPY求解。
cp.Problem(cp.Minimize(f), [A*x == b]).solve()
print("CVXPY得到的最优值: {}".format(f.value))

# 使用乘子法求解。
resid = A*x - b
y = cp.Parameter(shape=(m)); y.value = np.zeros(m)
aug_lagr = f + y.T*resid + (rho/2)*cp.sum_squares(resid)
for t in range(MAX_ITERS):
    cp.Problem(cp.Minimize(aug_lagr)).solve()
    y.value += rho*resid.value

print("乘子法得到的最优值: {}".format(f.value))

CVXPY得到的最优值: 5.5905035557463005
乘子法得到的最优值: 5.572761551213633