大规模稀疏线性规划求解思路梳理

锅逗逗

发布于 2022-08-01 14:56:12

1.5K0

发布于 2022-08-01 14:56:12

文章被收录于专栏：锅逗逗的杂学笔记

背景

2019年初秋，我接到了如下需求：

已知现在有M个广告主和N个广告词，其中每个单位流量的（广告主，广告词）收益固定，且每个广告主/广告词均有流量分配限制，问如何给（广告主，广告词）分配流量，使得收益达到最大。

这个需求是一个大规模稀疏线性规划问题，接下来本文将就上述需求描述如何加速求解。

0. 方案调研：Mosek

线性规划问题的求解快慢，既与迭代收敛速度有关，又和每轮迭代更新的速度有关。

通过调研，首先将Primal-dual和Mosek作为候选的求解方法

锅逗逗：内点法初探——线性规划标准形式下的求解思路

对比求解相同线性规划问题两种方法的收敛情况

上图显示了在10^4求解变量规模上Mosek和Primal-dual方法的收敛情况，可以看到Mosek方法比Primal-dual方法更快收敛。

最终基于Mosek方法来求解线性规划问题。

1. 化解约束方程

问题

Mosek方法要求将输入的约束化为标准型：

在需求中只包含不等式约束，目标变量x的取值范围为x>=0，且存在x=0的情况。

原始线性规划问题格式问题举例：

分析scipy.optimize.linprog预处理过程，发现该过程充斥着大量冗余的循环判断操作，且化成的标准型也并非最简模式。随着目标变量的增加，该预处理过程占用的时间也明显增加。

上述例子经scipy.optimize.linprog预处理后得到的标准型如下：

优化方法

结合需求中x=0或x>=0的特殊性质，采用以下步骤将目标问题化简成标准型：

step1: 将x=0变量从约束方程中消除；

step2: 检查约束方程中是否存在单变量约束，若存在，则根据单变量约束条件重新确定待求解变量x的取值范围，并将该约束方程剔除；

step3: 根据剩下约束方程和变量取值范围化为标准型。

最终得到的标准型如下：

结果

[1] 化简形式对比

优化后的方案能够将原线性规划问题化简成最简形式的标准型，进而减少变量/约束个数

[2] 化简耗时对比

将原线性规划问题化简成最简形式的标准型，进而减少变量/约束个数。

变量个数（广告主+广告词）	scipy.optimize.linprog	优化方案
100+1,000	6.29s	0.01s
100+5,000	142.17s	0.03s
100+10,000	559.68s	0.03s
100+20,000	2228.7s	0.16s

随着求解变量个数的增加，scipy.optimize.linprog的预处理过程耗时明显增加，且耗时不可忽视。相比之下，改进方案预处理过程耗时非常低，在求解大规模线性规划问题时可忽略不计。

2. 单轮迭代过程优化：Eigen CG

问题

采用Mosek方法求解线性规划问题，需要经过若干轮迭代才能获得最优解，每轮迭代包括以下四个步骤：

step1: 计算残差res

step2: 采用牛顿步得到二阶导矩阵M，结合残差res计算更新方向△，即求解线性方程组M△=-res

step3: 确定步长alpha

step4: 更新当前迭代点

其中，求解更新方向的过程决定了单次迭代的速度快慢。

scipy.optimize.linprog中采用scipy.sparse.linalg.splu方法，首先对矩阵M进行LU分解，再求解M的逆矩阵M^{-1}，最后计算△=-M^{-1}res。该过程为单线程计算，计算速度低效、不能满足例行更新的耗时需求。

优化

分析发现在Mosek方法涉及到的二阶导矩阵M是一个对称、正定、稀疏的方阵，可以采用共轭梯度法（Conjugate Gradient），通过直接求解线性方程组M△=-res得到△的值，共轭梯度法相较直接求解法（LU分解）速度更快且可并行化求解。

经过调研，使用Eigen::ConjugateGradient类对象来完成求解线性方程组的工作。

结论：求解相同的线性方程组，使用Eigen::ConjugateGradient的比scipy.sparse.linalg.splu具有优先一个量级的求解精度。

2.+ 加速线性方程组的求解：DPCG+ICCG

通过分析计算时间发现，尽管使用了Eigen的共轭梯度法来求解线性方程组，这个过程依旧非常耗时，所以优化重点在于进一步加速线性方程组的求解。

Diagonal Preconditioned Conjugate Gradient（DPCG）

直接使用共轭梯度（Conjugate Gradient）方法求解线性方程组的收敛速度完全依赖于线性方程组稀疏矩阵的条件数，通常很难在理想的迭代次数（几到几十步）获得解向量，CG方法通常需要和Preconditioner一起使用。

之前的实现中引用了第三方库Eigen::ConjugateGradient实现方程组的求解，其中Eigen::ConjugateGradient默认采用Diagonal Preconditioner，该Preconditioner是一个对角矩阵，该矩阵的逆矩阵非常容易求解。

通过统计Mosek方法每轮迭代中求解线性方程组的难易程度发现，随着Mosek方法迭代轮数的增加，求解线性方程组越来越困难（获得解向量的迭代次数增加），后期甚至到了无法接受的上千次迭代次数。

Incomplete Cholesky Conjugate Gradient （ICCG）

在Mosek方法论文中采用Choleksy方法分解系数矩阵求解线性方程组。该方法为直接求解法，能够一次获得方程组的解向量，

结合Cholesky和Conjugate Gradient，在CG迭代过程中将Diagonal Preconditioner替换成Incomplete Cholesky Preconditioner。

构建Incomplete Cholesky的主要工作如下：

a. Incomplete Cholesky方法在分解过程中保留系数矩阵的稀疏性，忽略Cholesky分解过程中产生的填充元。为了使Cholesky和Incomplete Cholesky的分解结果尽可能接近，使用Approximate Minimum Degree Ordering Algorithm对系数矩阵进行重排；

b. 运用Multifrontal方法构建组装树，使用需求提供的数据，通过分析发现组装树的深度接近2，第一层（叶子节点）个数接近广告词数量M，第二层（根节点）个数接近广告主数量N（具体情况与系数矩阵重排结果有关）。

采用的策略是在每次求解中开辟一个N*N的连续空间，首先分解第一层节点，再在N*N的空间里分解第二层节点，最后再更新第二层节点对应的元素。

c. 采用icfm方法对系数矩阵进行缩放求解，不同之处在对每行/列进行分解时保留原始元素的位置而非不保留最大的p个元素，只在对角线的计算上考虑填充元的信息。

DPCG和ICCG对比和使用策略

上图表示了在PCG方法中，使用Diagonal Preconditioner和Incomplete Cholesky Preconditioner随着Mosek方法迭代轮数的增加，求解线性方程组的难易程度情况。

可以发现，通过设计合理的Incomplete Cholesky Preconditioner，可以有效缓解求解线性方程组迭代次数爆炸的现象。（注：模拟数据中，使用Incomplete Cholesky Preconditioner最多需要37次可以得到方程解向量，而Diagonal Preconditioner需要4045次）

同时考虑到Diagonal Preconditioner求解过程比Incomplete Cholesky分解过程更容易，最终策略：在Mosek迭代初期系数矩阵条件数较低的前提下，先采用DPCG求解，待求解过程中迭代次数超过一定阈值时，再切换成ICCG方法进行求解。