上一节笔记:数值优化(A)——线性规划中的单纯形法与内点法
————————————————————————————————————
大家好!
这一节我们开始介绍二次规划的相关内容。二次规划也是一类具体的非线性规划的问题,也有对应的方法。
二次规划也是我们这一个系列正课的最后一部分,但这并不意味着最优化这一门学科的结束。在之后还会有一节比较贴近现代优化和机器学习领域的一些算法的思想简介,对它们有一定的了解,才算是对优化这一门学科的一个正式的涉猎。
那么我们开始吧。
二次规划其实本质上就是为了解决这么一个问题
其中
,
是标量,且
是
的一个向量,并且我们设
(当然了这就意味着
,
是一个对称矩阵。
要理解二次规划就必然需要理解二次函数和凸的概念。比方说在这里,如果我们设
是一个对称半正定矩阵,那么我们认为它就是一个凸的二次规划问题,如果
是一个对称正定矩阵,那么它就是一个严格凸的二次规划问题。当然了如果它是一个对称不定的矩阵,这个问题就是一个非凸的二次规划问题。总结起来可以说,这个问题是什么样的问题,完全取决于
的性质,这也不难理解,因为这个目标函数的海塞矩阵就是
。在这一节中,我们主要关心的是凸二次规划问题,也就是
对称半正定。
使用和上一节完全一样的策略,我们可以把问题再做一步标准化,就会得到下面的形式
其中
,
是行满秩的。
要求解这个问题,依然需要依赖KKT条件,根据KKT条件,可以得到
如果矩阵是可逆的话,那么这个方程组有解,但是遗憾的是并不是所有的情况下这个矩阵都可逆。那么究竟什么时候可逆呢?这就是下面这个引理了。
Lemma 1: 设
满足
,也即
的列组成了
的零空间,那么只要
是一个正定阵,
行满秩,则矩阵可逆,也即KKT条件的解
唯一。
我们证明一下这个结论。事实上要证明说矩阵可逆,其实就是要说明,如果我们有
那么一定会有
。这里用了一个线性代数的小技巧,就是左边再乘一个式子
,就可以得到
最后一个等式得到是因为
,而这个条件其实就对应了一开始问题中的
。而事实上,这个条件意味着
本身是落在
张成的线性空间内的,也就是存在
,使得
,所以我们可以根据上面的等式可以得到
,所以
(矩阵的正定性),所以
。
有了
之后,根据
,就可以得到
,就可以得到
(这是因为
是列满秩的)。
这个性质就意味着其实最关键的部分在于
的零空间的性质,这个地方的思路可以对比上一节的单纯形法,在单纯形法中我们要求松弛变量的分量为0,这二者也算是有一点可以联想的地方(虽然不是很严格就是了……)
那么有了这个之后,我们下面说的就是大家比较熟悉的二次规划问题中,解唯一的这么一个性质了。
Proposition 1: 设
行满秩,
正定,那么
是二次规划问题的唯一全局解。
这里
就是KKT条件得到的解。
我们证明一下这个结论。对于任意的
,考虑设
,那么我们的目的就是证明目标函数
一定比
要大,而引入
可以很方便的解决这个问题。
注意到我们有
这里要注意到我们的最后一个等号其实有一个跳步,要注意到
这个式子,而这可以通过
(想想为什么?)和KKT条件得到。
那么注意到
,说明
。这样的话,根据
是正定的,就可以得到
的结论,换句话说
是恒成立的,只要
,这就是我们需要证明的结论。
千万不要以为到上面就算结束了,对于方程组
的求解可没有那么容易。硬解自然是没问题的毕竟理论上都有保障,但问题在于数值优化是数值方法,如果没有一个好的计算方法支撑,就谈不上是计算数学的主题了。事实上关于这个方程组的求解,还是有很多值得说的地方,同时在这一部分我们也会详细的分析数值方法每一步的计算复杂度,毕竟这个不说的话,值得说的地方就又少了一大半……
Schur补方法也可以形象的被称为打洞法,基本的思想就是对当前的线性方程组的分块系数矩阵,我通过一个高斯消元的方式,使得我的分块系数矩阵的左下角为0(当然这一块还是要小心表述,毕竟我没有碰到过其它的情况,也不知道是否有在其它地方打0的操作……),所以面对现在的方程组,假如说我们修改它的形式为
(只是做了一个换元而已)那么我们的思路就是消去左下角的矩阵
,换句话说就是两边乘一个变换矩阵
也就是
把它化简一下可以得到
解这个方程组有两种思路,一种是分解法(Factorization Method),一种是迭代法(Iterative Method)。我们下面分别来解释这两种方法。
分解法的关键是通过已有的数值方法尽可能的减少时间复杂度。在这里注意到
是对称半正定阵,所以可以考虑Cholesky分解,也就是说可以设
它的复杂度是
,这是因为
。
所以有了
这个设定之后,我们就可以设
,所以下一步就是解
,这个解
可以通过
来做到。这里要注意一个计算顺序,我们对于这个方程组要按照这样的顺序来解
当然了如果你学过数值分析那么不会对这个拆解感到陌生。这里每一个方程组求解需要
的计算次数,所以合在一起就是
。
同样的代入
的表达式,我们可以得到
这里会牵涉到一个矩阵
,这个东西不是一个方阵,所以如果希望做数值分解,一个常用的方法就是QR分解,也就是说设
其中
是一个正交阵。在这一步的复杂度是
。当然了因为你要计算
这两个矩阵相乘,所以还会多一个
的复杂度。
这样的话有一个好处,就是我们可以发现最终解
的时候,
是用不上的(想想为什么?),最后就会变成解下面这个方程组
这一步的复杂度为
,这是因为
。
解
也很简单,直接考虑解
即可,复杂度为
。
我们一定要强调的是,这里的每一步的算法都要遵守数值算法的计算方式和计算顺序,比方说QR分解和Cholesky分解都是有自己的一套算法可以保证较低的计算复杂度的。所以如果不按照我们刚才说的每一步的顺序去求解,那么按照数值算法这样的拆解就没有任何的实际意义。
这里还有一个细节是,对于矩阵
,它本身是一个对称半正定阵,所以是可以使用Cholesky分解的,这样的话就不会用上QR分解,复杂度会稍微低一些。但是这里有一个麻烦是数值误差。简单来说,在矩阵条件数很高的时候,也就是矩阵病态的时候,Cholesky分解是容易崩溃的,这会给算法带来不稳定性。而QR分解就不会有这个麻烦,如果对这方面的细节感兴趣,可以查看一些矩阵论或数值线性代数的书籍,这里我们不再赘述,大家当作结论记住就好。
迭代法的思路就是不停的考虑使用共轭梯度(CG)方法来解线性方程组(见第3节的线性共轭梯度法,链接如下:
)。可以这么做的原因是优化问题和解方程组问题很多时候可以看成是同一种问题,比方说
等价于解方程组
。
因为我们的目的也是要尽可能的避免一些繁杂的计算(比方说求逆),所以在这里我们设
,那么我们只需要解
即可。同样的,对于同样包含
的矩阵
,我们代入
,就可以得到下面一个方程组
最后就解一下
即可。所以一共有三个方程组,每一个方程组都可以使用一个CG方法来进行求解。
当然了其实严格来说,第二个方程组也需要一步CG方法,这是因为我们按照顺序计算,应该是
中间的那个方程组依然是会依赖到CG算法的。
当然了有的人可能会问,为什么我们要考虑使用CG算法来进行计算求解,这个地方的原因是,CG算法的复杂度为
,其中
为迭代次数,取决于我们方程组的系数矩阵不相同特征值的个数。而对于Cholesky分解来说,这个复杂度为
(之前已经计算过),所以如果
,那么CG其实反而会慢一些。不过往往问题中的矩阵
的性质都比较好,所以迭代的时候往往CG这样的算法就会快一些。
最后提一下,Schur补方法在多元统计和回归分析中也有所应用,这里给大家放几篇文章供参考和引申阅读。
https://zhuanlan.zhihu.com/p/47122783
https://zhuanlan.zhihu.com/p/78884647
零空间法的本质是考虑到我们之前的理论需要依赖到
这样的表达式,因此必然会有一些方法会利用上这个性质。对于方程组
我们设
,并且设
,那么这样的话,
就是一个非奇异的矩阵,因为
相当于占据了
的行空间,也就是像空间。像空间和核空间拼在一起就是完整的一组基(当然这一部分工科的线性代数说的不是很多,可以参考线性映射中,像空间和核空间的相关理论)。
既然我们可以得到这样的一种分解,那么自然也可以做分解为
那么这样的话,就会得到我们的式子
第一个式子可以解出
(因为
是满秩的,对比一下矩阵的大小就可以看出来了),第二个式子左右两边都乘上一个
就可以得到
注意到
根据假设是正定的,这就意味着可以通过这个方程来解
。
都得到之后,就可以来解
了,而解
就可以根据方程
来得到。
所以这个方法的关键就是在于这个
应该怎么取。这里的一个方案就是考虑取
为对
做QR分解得到的结果。当然了这里的QR分解中
是列满秩的,不是满秩的,那么这个时候根据高等代数的知识,我们可以得到
与
所张成的空间相同。
当然了,这个方法也有计算比较复杂的地方,这就是方程
。这个方程的求解的简单方案自然是依赖
的,也就是对它做Cholesky分解。但是问题在于我们需要
的信息,而
的信息不是特别好找。这里有一个思路是考虑使用两边乘上
,利用
然后使用投影共轭梯度法(projected CG)来解决这个问题就好。可以参见Numerical Optimization的第461页。
到此为止,你也许明白了求解二次规划问题的两种方案,但是它们俩究竟哪个好呢?这就要看我们不同方法究竟会在哪里会花很多时间了。比方说对于Schur补方法,我们的计算依赖的是对
的处理,同时需要
是对称正定阵(否则的话有的方程就无法用CG求解)。所以如果这个处理不耗时,
又满足这个要求,那么Schur补方法就更好一些。不然的话零空间方法会更好,这是因为
一定要正定的,不然根据我们之前证明的性质,KKT条件的解就不好说唯不唯一了233。
激活集(Active Set)方法的思考角度有点类似上一节的单纯形法。我们考虑把我们的约束按照等式约束和不等式约束区分,也即是考虑
我们在上一节有介绍过说,可以通过一些转换把问题转成标准形式。所以这个形式和我们之前讨论的
是不矛盾的。
同样的我们这里假设我们的激活集为
(别忘了,这些集合都是指标集)所以解KKT条件可以得到
所以其实解这个KKT条件,如果我们的激活集已知,那么问题其实就是一个线性规划问题,用之前的方法求解即可。但是很多时候,激活集是不知道的,所以我们的思路是考虑取等式约束和一部分不等式约束,作为我们的工作集(working set)。这样取的思路就是,考虑以每一个工作集的约束下的极小值,然后一直更新工作集,一直到KKT条件满足为止(也就是对应的系数
均非负)。也就是说每一次都会考虑一个等式约束问题
这个方法的操作步骤如下:
为一开始的工作集,其中
是线性无关的,并且对
,有
,对
,有
。
,使得它满足
而把这个式子化简可以得到
,其中
是可以满足所有约束条件的最大的那个。如果存在
但是
,那么把
指标加入工作集,也就是
,否则工作集不变。
,如果存在某一个
,那么将
从工作集中移除,重新解
,否则工作集保持不变。
这里的一个关键是对于第3点的理解。自然的,如果说
,这就相当于说我们在当前约束下取得的极小值是可以无障碍达到的。但是有些情况是做不到这样的,可以画一张图来解释这个
在什么情况下是不能够取到1的。
所以可以看出,能够取到
自然是最好的,但是这里因为有一个约束(边界)挡住了,所以我们的
就是这个方向走到边界上的这个点。所以很多时候,这个要新添加进工作集的约束也会有叫阻挡约束(blocking constraint)的说法。
这么一说还是挺抽象的,不如我们来看一个具体的例子。
Example 1: 考虑问题
,约束条件为
我们可以大致画出这些不等式所划定的范围
其中的标号就代表上面按顺序写的约束的顺序。这几个边界组成的区域内部就是我们的约束集。
假如说我们一开始的迭代点为
(这里因为上下标会同时出现,所以我们用下标来表示分量,上标表示迭代的次数,和之前的含义不同,注意区分),那么对应的工作集就是
,所以这样的话我们就可以找到
,这是因为对应的约束下,只有
这一个点符合条件,没有办法迭代到任何下一个不同的点。
既然这个确定完之后,我们开始计算KKT条件的系数,以观察是否它是极小值点。解
(其中
为对应约束条件的系数,
,大家可以想想为什么是这个值)可以得到
,不满足非负的要求,说明没有找到最小值,要移除工作集的元素,这里我们移除
,因为它对应的拉格朗日系数的绝对值更大一点,所以这一步我们就完成了。
下一步,我们的工作集就只有一个了,也就是说
,
(因为前一步其实什么都没干),那么这个时候我们还是一样要找一个最优的下降方向,这个下降方向对应的解的问题就是
这里目标函数其实有个化简的技巧,就是我们要解的问题的格式为
,但是如果我们去掉了
这一项,剩下的部分就是原来的目标函数把
改成
而已,所以其实只需要把
算出来就好了。这一项就是
那么自然的,我们的解就是
,并且我们发现这个解可以一直走下去,也就是说
,工作集保持不变。
下一步,我们会发现,
,因为之前已经求到了最小值,工作集又没有发生变化。因此这个时候要开始验证KKT条件,以判断是否找到了极小值。计算
我们可以得到
,不是非负数,说明没有找到极小值,并且我们要把这最后的指标从工作集中剔除出去。所以
下一步,我们要在工作集为空,也就是无约束的情况下找搜索方向。注意到这里我们要解的方程就是
所以这会得到
,但是注意到这个时候这个方向是走不完的,因为当我们走到
的时候会受到第一个约束条件的限制。换句话说,这个时候应该把第一个约束添加到工作集中,所以这个时候
。
下一步,注意到这个时候要解的搜索方向满足
这就意味着
,并且这个方向可以走完,也就是说
,工作集不变。
下一步,因为工作集没变,刚才又找到了极小值,所以对应的
,验证KKT条件
可以得到
,是大于0的,因此解就算是找到了。
解释这个例子确实也花了挺长时间的,但是相信大家如果跟完了这个例子,对于激活集方法的具体步骤也不太可能陌生了。
关于二次规划的内容,其实还没有完全结束,不过碍于篇幅我们就这一节就先到这里,剩下的放到下一节再说。
这一节我们主要介绍了二次规划的几种依赖线性方程组求解的方法和与单纯形法的出发点类似的激活集方法。在这一节,我们也会发现,对于依赖线性方程组求解的方法,它们各有千秋,要理解它们计算的思路需要对数值分析的一些算法和它们的复杂度有比较深刻的了解。对于激活集方法,则更需要从实操出发,来理清算法求解一步步的思路。
在下一节我们会结束二次规划的内容,然后我们会开始做一些现代优化方法的概述和理解。了解了现代优化方法,才算是真正的在优化这个领域迈下坚实的一步~