本文将探究一个被称为二次规划的优化问题,这是一种特殊形式的非线性约束优化问题。二次规划在许多领域都有运用,比如投资组合优化、求解支持向量机(SVM)分类问题等。在R中求解二次规划有许多包,这次,我们将讨论一下quadprog包。在我们开始讲解案例之前,我们将先简短地介绍一下二次规划的机理。
对于一个二次规划问题,首先要考虑的就是一个二次目标函数:
Q(x)=12xTDx−dTx+c.
这里 x 在 ℝn 中是一个向量, D 是一个n×n 的对称正定矩阵,在 ℝn 中 d 是常数项约束,c 是一个标量常数。Q(x)函数通常以二次函数的形式出现,并且它高维的通项表达式是:
q(x)=ax2+bx+c
Q(x)的关键特性在于这是一个凸函数。
我们也对向量x构造一个线性约束集合,即x ∈ℝn。
我们把这些约束写成:
Ax=fBx≥g
这里,A 是一个 m1×n 的矩阵且约束为 m1≤n,BB 是一个 m2×n 的矩阵.向量 f 和向量 g的长度分别是m1和m2.
这是一种让我们可以充分考虑实际条件的标准型。比如我们让 x 强制满足
∑i=1nxi=1
的求和条件,或者满足ai≤xi≤bi的区间约束。接下来,我们将介绍如何将这些约束转化为矩阵表达。
用这个符号系统,我们可以简洁表示二次规划 (QP):
{minimizex∈ℝn:Q(x)=12xTDx−dTx+csubjectto:Ax=fBx≥g
考虑目标函数:
Q(x,y)==12[xy][2−1−12][xy]−[−32][xy]+4x2+y2−xy+3x−2y+4.
我们这个约束条件下的可行域内寻求最小化:
yyy≥≥≤2−x−2+x3.
我们可以找到这个可行域的顶点并在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")
SHAPE \* MERGEFORMAT
想要用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, 都存在MinQ(x,y)+c 和 Q(x,y)的解相等。因此,我们可以忽略二次规划中的常数项:
D=[2−1−12]d=[−32].
我们可以写出约束方程的矩阵形式:
⎡⎣⎢⎢1−1011−1⎤⎦⎥⎥[xy]≥⎡⎣⎢⎢2−2−3⎤⎦⎥⎥
因此:
A=⎡⎣⎢⎢1−1011−1⎤⎦⎥⎥Tb0=⎡⎣⎢⎢2−2−3⎤⎦⎥⎥
quadprog包默认是求解最小化问题,目标函数二次,约束一次。所以,我们的约束条件默认的形式也就是AX>=bvec。通常我们需要把一些原来是求极大值的问题或者<=约束通过乘以负号来转化。
这是R的完整实现:
· 参数Dmat表示海赛矩阵
· 参数dvet表示一阶向量,和Dmat的维数要相对应。
· 参数Amat表示约束矩阵,默认的约束都是是>=。
· 参数bvet表示右边值,由向量,和Amat的维数要相对应。
· 参数 meq 表示从哪一行开始Amat矩阵中的约束是需要被当作等式约束的。
(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)的全局最小值点,而黄点表示二次规划的解。
数据分析咨询请扫描二维码