[译]在R中使用quadprog包求解二次规划

概述

本文将探究一个被称为二次规划的优化问题，这是一种特殊形式的非线性约束优化问题。二次规划在许多领域都有运用，比如投资组合优化、求解支持向量机(SVM)分类问题等。在R中求解二次规划有许多包，这次，我们将讨论一下quadprog包。在我们开始讲解案例之前，我们将先简短地介绍一下二次规划的机理。

什么是二次规划

对于一个二次规划问题,首先要考虑的就是一个二次目标函数:$$ Q(x) = \frac{1}{2} x^T D x - d^T x + c. $$

这里 \(x\) 在 \(\mathbb{R}^n\) 中是一个向量, \(D\) 是一个\(n \times n\) 的对称正定矩阵，在 \(\mathbb{R}^n\) 中 \(d\) 是常数项约束，\(c\) 是一个标量常数。\(Q(x)\)函数通常以二次函数的形式出现，并且它高维的通项表达式是：$$ q(x) = ax^2 +bx + c$$\(Q(x)\)的关键特性在于这是一个凸函数。

我们也对向量\(x\)构造一个线性约束集合，即\(x\) \(\in \mathbb{R}^n\)。
我们把这些约束写成:$$Ax = f \qquad Bx \geq g$$

这里，\(A\) 是一个 \(m_1 \times n\) 的矩阵且约束为 \(m_1 \leq n \)，\(B\)B 是一个 \(m_2 \times n\) 的矩阵.向量 \(f\) 和向量 \(g\)的长度分别是\(m_1\)和\(m_2\).

这是一种让我们可以充分考虑实际条件的标准型。比如我们让 \(x\) 强制满足$$ \sum_{i=1}^n x_i = 1$$的求和条件，或者满足\(a_i \leq x_i \leq b_i\)的区间约束。接下来，我们将介绍如何将这些约束转化为矩阵表达。

用这个符号系统,我们可以简洁表示二次规划 \((QP)\):

$$\left\{ \begin{array}{l} \mathrm{minimize}_{x \in \mathbb{R}^n}: \qquad Q(x) = \frac{1}{2} x^T D x - d^T x + c \\ \mathrm{subject\; to:} \qquad Ax = f \qquad Bx \geq g \end{array} \right.$$

示例

目标函数

考虑目标函数：
$$
\begin{matrix} Q(x,y) &=& \frac{1}{2} \left[ \begin{matrix} x & y \end{matrix} \right] \left[ \begin{matrix} 2 & -1 \\ -1 & 2 \end{matrix} \right]\left[ \begin{matrix} x \\ y \end{matrix} \right] - \left[ \begin{matrix} -3 & 2 \end{matrix} \right] \left[ \begin{matrix} x \\ y \end{matrix} \right] + 4 \\ &=& x^2 + y^2 -xy +3x -2y + 4 . \end{matrix}
$$

约束条件

我们这个约束条件下的可行域内寻求最小化：

$$
\begin{matrix} y &\geq& 2 - x \\ y &\geq& -2 + x \\ y &\leq& 3. \end{matrix}
$$

我们可以找到这个可行域的顶点并在R画出整个可行域:

plot(0, 0, xlim = c(-2,5.5), ylim = c(-1,3.5), type = "n", 
     xlab = "x", ylab = "y", main="Feasible Region")
polygon(c(2,5,-1), c(0,3,3), border=TRUE, lwd=4, col="blue")

化为标准型

想要用quadprog包求解二次规划，我们需要同时转化我们的目标函数和约束条件为矩阵形式。这里是官方文档的说明：

This routine implements the dual method of Goldfarb and Idnani (1982, 1983) for solving quadratic programming problems of the form min(-d^T b + 1/2 b^T D b) with the constraints A^T b >= b_0.

可惜官方文档多可读性不高，我们很难得知如何准确地转化二次型\(Q(x,y)\)为一个矩阵形式。首先，我们观察到，对于任意常数 \(c\)， 都存在\(Min Q(x,y)+c\) 和 \(Q(x,y)\)的解相等。因此，我们可以忽略二次规划中的常数项:

$$D = \left[ \begin{matrix} 2 & -1 \\ -1 & 2 \end{matrix} \right] \qquad d = \left[ \begin{matrix} -3 \\ 2 \end{matrix} \right]. $$

我们可以写出约束方程的矩阵形式:

$$\left[ \begin{matrix} 1 & 1 \\ -1 & 1 \\ 0 & -1 \end{matrix} \right] \left[ \begin{matrix} x \\ y \end{matrix} \right] \geq \left[ \begin{matrix} 2 \\ -2 \\ -3 \end{matrix} \right] $$

因此：
$$A = \left[ \begin{matrix} 1 & 1 \\ -1 & 1 \\ 0 & -1 \end{matrix} \right]^T \qquad b_0 = \left[ \begin{matrix} 2 \\ -2 \\ -3 \end{matrix} \right]$$

具体实现

quadprog包默认是求解最小化问题，目标函数二次，约束一次。所以，我们的约束条件默认的形式也就是AX>=bvec。通常我们需要把一些原来是求极大值的问题或者<=约束通过乘以负号来转化。

这是R的完整实现:

require(quadprog)
Dmat <- 2*matrix(c(1,-1/2,-1/2,1), nrow = 2, byrow=TRUE)
dvec <- c(-3,2)
A    <- matrix(c(1,1,-1,1,0,-1), ncol = 2 , byrow=TRUE)
bvec <- c(2,-2,-3)
Amat <- t(A)
sol  <- solve.QP(Dmat, dvec, Amat, bvec, meq=0)

• 参数Dmat表示海赛矩阵

• 参数dvet表示一阶向量，和Dmat的维数要相对应。

• 参数Amat表示约束矩阵，默认的约束都是是>=。

• 参数bvet表示右边值，由向量，和Amat的维数要相对应。

• 参数 meq 表示从哪一行开始Amat矩阵中的约束是需要被当作等式约束的。

结论

让我们检查一下 quadprog 求解器求解的结果吧：

> sol
$solution
[1] 0.1666667 1.8333333
$value
[1] -0.08333333
$unconstrained.solution
[1] -1.3333333  0.3333333
$iterations
[1] 2 0
$Lagrangian
[1] 1.5 0.0 0.0
$iact
[1] 1

\((1/6,11/6)\) 点是唯一满足约束条件和 \(Q(x,y)\)的最小化目标，但 \((−4/3,1/3)\)点才是 \(Q(x,y)\) 的最小值点。iterations,Lagrangian 和 iact 都是用来描述quadprog算法性能的。对于这些值之后我们将进一步讨论。现在，让我们先可视化二次规划的解。为此，我们在\(Q(x,y)\)的可行域边界添加一个外侧的等高线图。

在图中，深绿色区域表示\(Q(x,y)\) 表面目标函数值较小的解,而亮色表示目标函数值较大的解。红点是\(Q(x,y)\)的全局最小值点，而黄点表示二次规划的解。

# 与可行域叠加等高线图
require(lattice)
qp_sol <- sol$solution
uc_sol <- sol$unconstrained.solution
x <- seq(-2, 5.5, length.out = 500)
y <- seq(-1, 3.5, length.out = 500)
grid <- expand.grid(x=x, y=y)
grid$z <- with(grid, x^2 + y^2 -x*y +3*x -2*y + 4)
levelplot(z~x*y, grid,  cuts=30,
    panel = function(...){
      panel.levelplot(...)
      panel.polygon(c(2,5,-1),c(0,3,3), border=TRUE, lwd=4, col="transparent")
      panel.points(c(uc_sol[1],qp_sol[1]),
             c(uc_sol[2],qp_sol[2]),
             lwd=5, col=c("red","yellow"), pch=19)},
    colorkey=FALSE,
    col.regions = terrain.colors(30))

• 英文原文地址

作为分享主义者(sharism)，本人所有互联网发布的图文均遵从CC版权，转载请保留作者信息并注明作者 Harry Zhu 的 FinanceR专栏:https://segmentfault.com/blog/harryprince，如果涉及源代码请注明GitHub地址：https://github.com/harryprince。微信号: harryzhustudio
商业使用请联系作者。