Category: 数学规划

优化算法在应用问题中的常见技巧

这是今天去Hulu交流时用的PDF,比较简短。总结了在做应用问题时的一点点感觉。
点击下面链接下载。
Techniques for Optimization Methods in Applications

0-1背包问题与子集合加总问题的近似算法

最近没有怎么更新博客,因为一直比较忙。最近发现所里在做的一个项目中,可以抽出一部分内容和0-1背包问题、子集合加总问题非常相似(虽然表面上不容易看出相似点),所以看了一些这方面的资料和论文,这里主要对问题特点和算法思想做一些整理。
这类问题其实很有意思,做数学和做计算机的人都会研究,而且我这里将要提到的论文都是做计算机的人所写的。

问题简述
0-1 Knapsack Problem (0-1背包问题,下面简称KP)和Subset Sum Problem (子集合加总问题,下面简称SSP)是经典的NP完全问题。
两个问题简要描述如下:

KP:有$n$个物品要放入背包,第$i$个物品的价值为$c_i$,占据体积为$v_i$,背包的总容积为$V$,要选择一部分物品放入背包,使得他们的总价值最大。对应的优化问题是
$\max_{x_i} \sum c_i * x_i $
$\mbox{s.t.} \sum v_i * x_i \le V, x_i \in \{ 0, 1 \} $
这里$x_i$代表是否选取第$i$个物品进背包,等于$1$就代表放入背包,等于$0$代表不放入背包。

SSP: 给一个集合$\{c_1, c_2, \ldots, c_n \}$,还有一个目标值$V$,问能否选出一个子集,使得子集中元素求和刚好等于$V$。我们一般考虑的是他的另一种表述方式:选出一个子集,使得子集中元素求和不超过$V$,且尽量大。对应的优化问题是
$\max_{x_i} \sum c_i * x_i $
$\mbox{s.t.} \sum c_i * x_i \le V, x_i \in \{ 0, 1 \} $
这里$x_i$代表是否选入子集,等于$1$就是选入子集,等于$0$就是不选入子集。

SSP是KP的特殊情况,也即当$c_i = v_i$的时候,KP退化为SSP,从问题式子上看,也完全一样了。
尽管如此,研究了KP不代表就不用研究SSP了,后面会说明这一点。

精确算法与近似算法
这两个问题都有很简单的动态规划算法可以精确求解,但可惜算法的时间复杂度是伪多项式的,也即和$V$相关,但$V$不是问题输入数据的规模,$n$才是。在ACM竞赛等算法比赛中,经常会遇到一些问题属于KP的变种,而伪多项式算法也就足够了。由于网上资料很多,而且难度不大,这里就不详细介绍了。如果你不知道,请你搜索“动态规划求解0-1背包问题”。
这里我们更关心多项式近似算法,也即PTAS(Polynomial Time Approximation Scheme),也即对任意给定的$\epsilon$,算法可以在关于$n$的多项式时间内求得一个解,且该解和真实最优解的最多相差$\epsilon$倍。
实际上这两个问题都有FPTAS(Fully PTAS)算法,也即时间复杂度不仅要是$n$的多项式,而且还要是$\frac{1}{\epsilon}$的多项式。

KP的近似算法
KP有非常简单的近似算法。
对任意$\epsilon$,令$P = \max_i c_i$,$K = \frac{\epsilon * P } { n}$,记$c_i’ = \lfloor \frac{c_i}{K} \rfloor$ (这里$\lfloor \cdot \rfloor$表示向下取整)
则使用$c_i’$代替$c_i$求0-1背包问题,求得的$x$即为$\epsilon$近似解。这是因为每一个物品由于取整而产生的误差都很小。
而这时候,由于$c_i’$都是整数,可以使用动态规划算法求解,容易算得动态规划的状态数至多$n \lfloor \frac{ P}{K} \rfloor$个,因此整个算法的时间复杂度为$\mbox{O}(n^2\lfloor \frac{ P}{K} \rfloor) = O(n^2 \lfloor \frac{n}{\epsilon} \rfloor ) $,
算法复杂度之所以能够达到多项式,是由于取整减小了动态规划的状态数,而$K$的出现是为了控制因取整而产生的误差,$K$足够小的时候就能保证解的相对误差不会超过$\epsilon$
此算法参考 Vijay V. Vazirani写的书 Approximation Algorithms(点击下载) , 第8.1,8.2章节(第69页-70页)
虽然这不一定是最快的多项式近似算法,但暂时没看到别的算法在任何情况下复杂度都能更低。
我自己也没有想到更好的。而有趣的是我最近在思考另一个等价的问题,但在发现等价性之前,我也想出了一个思路几乎相同,复杂度完全一致的算法。现在看来,其实本质是同一个算法。

SSP的近似算法
根据前面所述,SSP是KP的特殊情况。所以KP的近似算法就可以认为是SSP的近似算法了。但是,之所以单独考虑SSP,是因为他可以求得更快。
SSP也有思路非常类似的算法,可以在一些国外的教学资料上看到。如这一篇(点击下载)。但是时间复杂度都比较高,事实上,SSP有$\mbox{O}(\min \{ \frac{n}{\epsilon} , n + \frac{1}{\epsilon^2} \log(\frac{1}{\epsilon}) \} )$的算法。我认为这个结论是比较漂亮的,复杂度严格比上面提到的KP的算法要低,而且低很多。
为了说明算法思想,我们先再深入理解一下前面的算法。
前面的算法是用取整来控制复杂度,取整的本质事实上是将状态空间划分成若干段,每一段只取一个端点作为“代表”,这样状态数就会是多项式的,总的算法复杂度也就是多项式的。在SSP中,我们可以不事先取整,但仍然将状态空间分段,让每个段记录两个数值(也即选两个“代表”),一个是落在该段上的最小子集和(最小代表),一个是落在该段上的最大子集和(最大代表)。这样我们的状态数相当于比前面KP的算法多了一倍,而且是需要动态维护的,但是由于误差控制的更好(极端的例子是如果集合只有一个数,这么做就没有误差),没有取整那么大的误差,所以不需要分那么多段就能保证解的相对误差不超过$\epsilon$,因此最终算法复杂度会低。
此算法参考 Hans Kellerer等人的论文 An efficient fully polynomial approximation scheme for the Subset-Sum Problem (点击查看) Journal of Computer and System Sciences 66 (2003) 349–370
可以发现,这个算法是不能简单推广到0-1背包问题上的。

SSP的拓展问题——ISSP
Anshul Kothari等人在2005年研究Uniform-Price Auction Clearing(其实我也不知道这是啥问题)的论文中提出了Interval Subset Sum Problem (下面简称ISSP),他是SSP的一个推广。这个问题没有前面的问题这么有名(事实上我没看到第二篇论文讨论这个问题),但刚好我也发现了另外一个实际问题可以抽象成这个模型,所以我觉得还是很有意义的。问题简要描述如下:
给定$n$个区间$[a_i,b_i]$,对每个区间可以选择上面的一个数,也可以什么都不选,使这些数的和不超过目标值$V$,且尽量接近。对应的优化问题是
$\max_{x_i,y_i} \sum y_i * x_i $
$\mbox{s.t.} \sum y_i * x_i \le V, x_i \in \{ 0, 1 \}, y_i \in [a_i, b_i] $
至少从形式上看,他不再是线性的了,而且当$a_i=b_i$时退化为SSP。所以似乎更难了。
实际上这个问题也能改写为线性模型。
$\max_{x_i,y_i} \sum y_i $
$\mbox{s.t.} \sum y_i \le V, x_i \in \{ 0, 1 \}, y_i \in [x_i * a_i, x_i * b_i] $
从直观上看,如果$a_i$很小,$b_i – a_i$很大,这个问题似乎很容易解。我已经可以给出看上去不是太强的条件,在该条件下,ISSP是多项式可解的,具体这里就不给出了。但无论如何,In general,ISSP也是NP难问题。
Anshul Kothari等人的论文中也给出了一个FPTAS算法,比较有趣的是,这个算法和上面的SSP的算法思想是一致的,时间复杂度也是一致的。因此这一算法应当是受到上面所述的SSP的算法启发而设计的,而且是一个比较简单的推广。
此ISSP算法参考论文为 Interval Subset Sum and Uniform-Price Auction Clearing (点击查看) Computing and Combinatorics, Lecture Notes in Computer Science Volume 3595, 2005, pp 608-620

一些启示
1,SSP是KP的一个特殊情况,但是他们的近似算法复杂度却可以完全不同,这说明特殊问题可以有特殊算法。使用General的算法是可以的,但需要仔细想想问题是否有特殊的地方没发现,是否有更具针对性的算法没找到。
2,ISSP是SSP的一个推广,也即SSP是ISSP的特殊情况。但事实上ISSP并没有更难,甚至某些时候很可能有多项式算法,这说明推广了,也有可能变简单,特殊情况也有可能更难。但是特殊情况的算法思想可以借鉴用来解决推广了的问题。

如何快速计算交叉项求和——从libFM联想到的一类数学问题

libFM里面有一个很好的idea是遍历特征的交互作用。也即$\sum_{i \ne j} x_i * x_j$。
但是遍历交互作用需要计算$O(n^2)$次乘法,于是作者做了一个变换,成为$ 0.5 * ((\sum_i x_i)^2 – \sum_i x_i^2) $。变换后只需要线性次的乘法和加法即可。
听严强说,实际使用的时候常常不需要遍历所有交叉项,因为很多特征之间是没有关系的。如果那样,问题就变成了$v = \sum_{(i,j)\in S} x_i * x_j $,其中$S$是给定的指标集。可能在实际中$S$里面元素是比较少的,但无论如何,这引发一个数学问题:如何给出一个等价变换,使得用尽量少的乘法和加法就能求得这个式子?如果设计一个智能算法,自动选择有用的特征对,这时候可能就真的需要一个算法找到相应的等价变换式。
确定这个定价变换可以不用特别快速,因为只需计算一次,但是得到的结果应该是可以尽量快速计算的。

这和我以前考虑过的一个问题非常类似,但是不完全一样。
首先我们考虑“相对精确”地去找计算量最少的等价变换。说“相对精确”是因为我还没想到以一个便捷的方式处理完美的精确。我们考虑近似化为一个优化问题。
假设我们要找到一个需要至多$m$次乘法的等价表达式,则有$v = \sum_{k=1}^m \left( \sum_{i=1}^n (a_{ki} * x_i) * \sum_{j=1}^n (b_{kj} * x_j) \right)$,这里$a_{ki}$和$b_{kj}$都只从$1$,$-1$,$0$之间选择,所以他们与$x$的乘法实际是不需计算的。这样的话,根据对应项系数相等(也可以是不等号,只要覆盖到交叉项也可以认为考虑了这种特征组合,不一定非要系数为1),我们会得到$n*(n+1)/2$个约束,未知变量数目是$m*2*n$个,从自由度的角度考虑,$m$至少需要和$n$一个量级才有可行解,实际上只要$m$取$n$就一定有可行解了,可以举出例子来,只是那样可能需要计算非常多次加法,所以不是好办法。我们可以用运算量作为目标,假设一次加法的代价为$C_{add}$,一次乘法的代价为$C_{mul}$,则目标函数为$\min. C_{add}*\sum_{k,i} |a_{ki}| + C_{add}*\sum_{k,j} |b_{kj}| + C_{mul}*\sum_{k=1}^n \max_{i,j} ( |a_{ki}|, |b_{kj}| )$。这个目标其实并不精确,不过不妨碍理解问题。
可以看到,目标是线性的,约束是二次的,且为整数规划,所以非常难。如果放缩到$a$, $b$都为实数,则问题简化了一些。但是仍然是非常困难的问题,非线性,非凸,不过可以尝试计算。虽然不能保证得到最优解,但是可能会得到较好的解。一个小问题是如果$a$, $b$为实数则他们和$x$乘积的运算量稍微不好量化一些。这里面很多地方描述不精确,比如$a$, $b$确实也不需要是整数,只是取特殊值相当于少了一次计算。
一句话描述这个优化问题,就是:给定一个0-1矩阵,把他化成尽量少的rank=1的小矩阵相加,且小矩阵等于两个尽量稀疏的向量相乘。

也可以有启发式的办法,例如交叉项相似的分到一组,可以用一次乘法计算出来。但效果很难说。

和这个问题很相似的问题还有:使用最少的乘法,求矩阵乘。比如Strassen algorithm就使用了7次乘法计算2*2的矩阵乘,而传统方法需要8次乘法。如何使用更少的乘法计算矩阵乘,和上述问题很类似。这也是一个非常困难的问题,目前还没有人能精确指出任意阶矩阵乘法所需的最少乘积次数。

VPS、全局优化、Python、并行

说说最近的事,找不到一个合适的标题,就以若干关键字为题好了。

最近购买了一个VPS,新年特价,30刀一年,2G内存,续费仍然是30刀一年。对,你没看错,没有少写一个0
据说这家超售很严重,但是技术水平不错,不大看得出来。
有了VPS就可以继续肆无忌惮的挂网站、爬网页什么的了。实测CPU很给力,跑程序很快。但是毕竟是便宜货,刚入手两星期左右,今天差不多挂了10个小时才修复,看在价格上也就不说什么了。
科院选课助手 ishangke.net 目前已经迁移过去了。由于配置的提升,一些耗费内存、CPU的操作有明显变快。
预计一个月内可能会上线一个和音乐有关的小站,请期待。兴奋的是,以前从未发布过类似的东西。不过悲观的认为,不会有几个用户吧。
最近在想一个相关的问题,我如何通过网上能获得到的数据,进行音乐推荐呢?目前还没好主意。

这段时间在外访问,生活上比国内要无聊很多。这边的老师在做一些全局优化的问题,针对的是实际工程上的问题。我完全不懂背景,单从数学上来看,就是说有一个非常复杂的函数,复杂到求一次函数值需要解一大堆东西。函数本身也写不出表达式,只是给自变量能求出函数值而已。目的是用最短的时间尽可能找到函数最小值。由于求一次函数值花费时间很多,比如要将近一分钟。所以花费的时间就约等于函数值计算的次数。也就是要尽量少的找一些点求值,却希望找到最小值点。基本思想就是用一个函数去插值,但是为了求解到全局最优,又不能一直去找一堆距离很近的点,也要去尝试一些“未探索的区域”。所以对这两个目标进行一种权衡。

这边的程序都有Matlab和Python两个版本。老师说是Matlab要付费所以后来的程序都选择编写Python版本,于是果然看到了国外开始用Python进行科学计算的实例了。但是我用的太少,稍微补了一下Python和NumPy的基本知识,略有收获。第一步的任务只是用现有的程序去测试一个实际问题,感觉难度不大,具体意义大不大,由于不了解问题背景,我也没有发言权。我觉得有必要做点什么在国内宣传一下使用Python做科学计算,但是还没找到力所能及且自己觉得有意思的事情。而且我自己也才入门水平。

我感觉他这里的并行应该是希望能到集群上跑,但是他这里没有集群环境,用的是一些配置很不错的多核服务器。然后程序采用的是单机多进程,不是多线程。所以逻辑上来说移植到集群上并行跑问题不是太大,只要找到相应的工具应该就能办到。

随机样本选择——快速求解机器学习中的优化问题

前阵子去参加了数学规划会议,报告很多,人也很多。或者说报告和人过多了……
有少数感兴趣的报告,这里谈一下全场最后一个报告。报告人是Jorge Nocedal,就是著名的LBFGS的作者。
他关注的问题是一类机器学习中非常常见的优化模型:
$J(w) = \frac{1}{N} \sum_{i=1}^N l(f(w; x_i); y_i)$
他谈到了机器学习里面常用的随机梯度下降法(stochastic gradient descent, SGD)。在实际大数据的情况下,他比经典优化方法效果要好。因为经典优化方法算一次下降方向就已经需要很大计算量了。
SGD需要的一次迭代计算量非常小,但是代价是他非常可能不是下降方向,而且需要一个有经验的人选择很好的步长(学习速率)。并且,在理论上,SGD的性质少,所以在优化这个方向里没什么人提到。即使是在计算上,Nocedal也提到SGD有个很大的缺点是不容易并行。
以上这些现象我在学习推荐系统SVD模型的时候也有一些感触。
所以他采用了另外一个思路,思想简单,效果却很好。而且能很好的利用上经典优化算法。经典算法都需要计算梯度。他的思想就是:只选取一部分样本来计算梯度。
这样有一个近似的梯度计算方法,用
$J_S(w)=\frac{1}{|S|} \sum_{i \in S} l(f(w; x_i); y_i)$
代替原来目标函数,求它的梯度。
其中S是样本集,样本集的大小相对于N来说很小。但是他也不能太小。只要有了梯度,很多算法都能用上了,比如说他文章中提到梯度方法,Newton-CG方法。S中的元素可以在每次迭代的时候随机选取。
但是需要面临的一个问题就是如何确定S的大小。他给出了一个公式用来确定S的大小是否合适。
$\| \nabla J_S(w) – \nabla J(w) \|_2 \le \theta \| \nabla J(w) \|_2 $
这个式子的含义就是样本集上的近似梯度和真实梯度足够接近。当然这个又回到了老问题,就是真实梯度计算量太大。于是他又做了一些近似变换,使得只需要计算在样本集S上的方差即可近似估计上式是否满足。
$ \frac{ \| Var^2_S( l(w; i) ) \|_1 }{ \| S \| } \le \theta ^2 \| \nabla J(w) \|_2 ^2 $
这样就可以每次迭代验证一下S是否足够大,如果不够大就扩大一下。
这个技巧部分解决了我对大规模优化计算的一些困惑,他可以利用上很多经典优化算法,并且又能让每次迭代的运算量不会太大。这个技巧也使得并行计算容易用上去了。
后面考虑去使用一下,算算某个实际例子测试一下效果。
Nocedal这篇论文的链接:http://users.eecs.northwestern.edu/~nocedal/PDFfiles/dss.pdf

为什么很多网络流问题总有整数解——全单模矩阵的性质

以前在Online Judge上遇到过很多问题需要转化为网络流模型,原问题一般是要求有整数解的。但是转化为网络流问题后虽然边上的最大流上限都是整数,其实并没有要求最大流也是整数,可求出来的最大流就是整数。我以前一直迷茫这个现象,并且我是很懒的人,一直没有严谨的想过这个问题或者查过这个问题。我们课题组平时很少讨论整数规划的问题,所以这个其实比较基本的问题我一直都不知道。其实你可以从直观感觉如果有非整数的流,那总应该能再怎么样多流一点变成整数,不过这显然是直观感觉,什么也说明不了。昨天看书看到有讲这些,终于理解了,在这总结一下。

Continue reading »