基础的派生指南

本文档将介绍计算优化问题解的导数的基础知识。导数可用于进行**敏感性分析**，以查看在参数发生小变化时解会如何变化，以及计算解的标量函数的**梯度**。

在本文档中，我们将考虑一个简单的有规范几何规划。考虑的几何规划问题为

\[\begin{split}\begin{equation} \begin{array}{ll} \mbox{最小化} & 1/(xyz) \\ \mbox{满足} & a(xy + xz + yz) \leq b\\ & x \geq y^c, \end{array} \end{equation}\end{split}\]

其中 \(x \in \mathbf{R}_{++}\), \(y \in \mathbf{R}_{++}\), 和 \(z \in \mathbf{R}_{++}\) 是变量，而 \(a \in \mathbf{R}_{++}\), \(b \in \mathbf{R}_{++}\) 和 \(c \in \mathbf{R}\) 是参数。向量

\[\begin{split}\alpha = \begin{bmatrix} a \\ b \\ c \end{bmatrix}\end{split}\]

是参数向量。

import cvxpy as cp


x = cp.Variable(pos=True)
y = cp.Variable(pos=True)
z = cp.Variable(pos=True)

a = cp.Parameter(pos=True)
b = cp.Parameter(pos=True)
c = cp.Parameter()

objective_fn = 1/(x*y*z)
objective = cp.Minimize(objective_fn)
constraints = [a*(x*y + x*z + y*z) <= b, x >= y**c]
problem = cp.Problem(objective, constraints)

problem.is_dgp(dpp=True)

True

请注意关键字参数``dpp=True``。参数必须按照特定规则输入DGP问题中，我们将其称为``dpp``。DPP规则在一个`在线教程<https://www.cvxpy.org/tutorial/advanced/index.html#disciplined-parametrized-programming>`__中有描述。

接下来，我们解决这个问题，并将参数:a、:b和:c设为2、1和0.5。

a.value = 2.0
b.value = 1.0
c.value = 0.5
problem.solve(gp=True, requires_grad=True)

print(x.value)
print(y.value)
print(z.value)

5612147353889386
31496200373359456
36892055859991446

请注意关键字参数 requires_grad=True；这是后续计算导数所必需的。

解决方案映射

上述问题的**解决方案映射**是一个函数

\[\mathcal{S} : \mathbf{R}^2_{++} \times \mathbf{R} \to \mathbf{R}^3_{++}\]

它将参数向量映射到最优解向量

\[\begin{split}\mathcal S(\alpha) = \begin{bmatrix} x(\alpha) \\ y(\alpha) \\ z(\alpha)\end{bmatrix}.\end{split}\]

这里，\(x(\alpha)\)、\(y(\alpha)\) 和 \(z(\alpha)\) 是与参数向量相对应的变量的最优值。

例如，我们刚刚看到

\[\begin{split}\mathcal S((2.0, 1.0, 0.5)) = \begin{bmatrix} 0.5612 \\ 0.3150 \\ 0.3690 \end{bmatrix}.\end{split}\]

敏感性分析

当解映射可微分时，我们可以使用其导数

\[\mathsf{D}\mathcal{S}(\alpha) \in \mathbf{R}^{3 \times 3}\]

进行**敏感性分析**，即研究在参数微小变化时解会如何变化。

假设我们将参数扰动一个小幅度的向量 \(\mathsf{d}\alpha \in \mathbf{R}^3\)。我们可以通过导数来近似计算由于扰动而导致的解的变化 \(\Delta\)，即

\[\Delta = \mathcal{S}(\alpha + \mathsf{d}\alpha) - \mathcal{S}(\alpha) \approx \mathsf{D}\mathcal{S}(\alpha) \mathsf{d}\alpha.\]

我们可以使用CVXPY来计算这个问题，如下所示。

将扰动分解为：

\[\begin{split}\mathsf{d}\alpha = \begin{bmatrix} \mathsf{d}a \\ \mathsf{d}b \\ \mathsf{d}c\end{bmatrix}.\end{split}\]

我们将参数的“delta”属性设置为它们的扰动值，然后调用“derivative”方法。

da, db, dc = 1e-2, 1e-2, 1e-2

a.delta = da
b.delta = db
c.delta = dc
problem.derivative()

“derivative”方法通过副作用将变量的“delta”属性填充为变量的预测变化。我们可以将预测值与扰动问题的实际解进行比较。

x_hat = x.value + x.delta
y_hat = y.value + y.delta
z_hat = z.value + z.delta

a.value += da
b.value += db
c.value += dc

problem.solve(gp=True)
print('x: 预测值 {0:.5f} 实际值 {1:.5f}'.format(x_hat, x.value))
print('y: 预测值 {0:.5f} 实际值 {1:.5f}'.format(y_hat, y.value))
print('z: 预测值 {0:.5f} 实际值 {1:.5f}'.format(z_hat, z.value))

a.value -= da
b.value -= db
c.value -= dc

x: 预测值 0.55729 实际值 0.55734
y: 预测值 0.31783 实际值 0.31783
z: 预测值 0.37179 实际值 0.37175

在这种情况下，预测值和实际解非常接近。

梯度

我们可以计算解的标量值函数相对于参数的梯度。设 \(f : \mathbf{R}^{3} \to \mathbf{R}\)，假设我们想计算组合函数 \(f \circ \mathcal S\) 的梯度。根据链式法则，

\[\begin{split}\nabla f(S(\alpha)) = \mathsf{D}^T\mathcal{S}(\alpha) \begin{bmatrix}\mathsf{d}x \\ \mathsf{d}y \\ \mathsf{d}z\end{bmatrix},\end{split}\]

其中 \(\mathsf{D}^T\mathcal{S}\) 是导数算子的伴随（或转置），而 \(\mathsf{d}x\)、\(\mathsf{d}y\) 和 \(\mathsf{d}z\) 分别是 \(f\) 相对于其参数的偏导数。

我们可以使用 CVXPY 的 backward 方法计算梯度。例如，假设

\[f(x, y, z) = \frac{1}{2}(x^2 + y^2 + z^2),\]

使得 \(\mathsf{d}x = x\)， \(\mathsf{d}y = y\)，和 \(\mathsf{d}z = z\)。假设 \(\mathsf{d}\alpha = \nabla f(S(\alpha))\)，并且假设我们从参数中减去 \(\eta \mathsf{d}\alpha\)，其中 \(\eta\) 是一个正常数。使用下面的代码，我们可以比较 :math:`f(mathcal S(alpha - eta mathsf{d}alpha))`与梯度预测的值，

\[f(\mathcal S(\alpha - \eta \mathsf{d}\alpha)) \approx f(\mathcal S(\alpha)) - \eta \mathsf{d}\alpha^T\mathsf{d}\alpha.\]

problem.solve(gp=True, requires_grad=True)

def f(x, y, z):
    return 1/2*(x**2 + y**2 + z**2)

original = f(x, y, z).value

x.gradient = x.value
y.gradient = y.value
z.gradient = z.value
problem.backward()

eta = 0.5
dalpha = cp.vstack([a.gradient, b.gradient, c.gradient])
predicted = float((original - eta*dalpha.T @ dalpha).value)

a.value -= eta*a.gradient
b.value -= eta*b.gradient
c.value -= eta*c.gradient
problem.solve(gp=True)
actual = f(x, y, z).value

print('original {0:.5f} predicted {1:.5f} actual {2:.5f}'.format(
       original, predicted, actual))

原始值 0.27513 预测值 0.22709 实际值 0.22942