Interior Point Method (IPM)——内点法复杂度分析

chen_koen

已于 2022-08-10 18:39:13 修改

阅读量4.9k

点赞数 8

文章标签：算法机器学习人工智能

于 2022-08-08 22:15:20 首次发布

本文链接：https://blog.csdn.net/qq_40107455/article/details/126235096

版权

写在前面：当前不同文章对IPM的复杂度分析结果有所不同，本人在这里根据自己的理解进行了总结，读者如有不同的意见欢迎写在评论区。

首先需要说明的是IPM是用于求解凸问题的常用有效方法。当前在Boyd的凸优化书籍分了两类，即

1.障碍方法（Barrier method）

2.原对偶内点法（Primal-dual interior-point method）

由于原对偶内点法具有超线性收敛性质，它通常比障碍方法更有效，特别是在一些高精度求解需求时。现有一些凸问题的通用求解器也是基于原对偶内点法的，例如CVX的SDPT3。基于这个原因，现有文献在分析IPM算法的复杂度时，通常默认是对原对偶内点法进行复杂度分析。注意，下面没特别说明，也是基于原对偶内点法进行复杂度分析。

在进行IPM的复杂度分析时，自然会想到至少五点（后续根据读者意见进行修改或新增）：

1.IPM复杂度分析的思路应该是怎样的？

答：需要说明的是IPM的复杂度分析是基于worst-case。

首先，IPM是迭代算法，所以需要关心IPM的迭代次数。这里，有个小插曲，IPM是现在大家比较多的叫法，其实以前有个叫法——路径跟踪法（path-following），其原因是迭代过程中得到的点会紧密跟踪中心路径移动。至于什么是中心路径？可以参考Boyd的凸优化书籍。

其实现在IPM的迭代次数分析是源于路径跟踪法的迭代次数，根据Wright&1997&SIAM，路径跟踪法使得解的精度（duality gap）达到 $\epsilon$ ，一般需要执行的迭代次数（至多）为 $O(n^\omega \log(1/\epsilon ))$ 。其中， $n$ 是优化变量的维度，如果优化变量是方矩则 $n$ 是行或列数（其原因是每次迭代时只关注 $\omega$ 的取值去减少duality gap，并不关注 $n$ ）。

看来IPM的迭代次数跟 $\omega$ 有关，而且 $\omega$ 越小越好。为了让这个 $\omega$ 变小，人们做了很多努力。现在，已知最小 $\omega$ 的是 $\frac{1}{2}$ ，是通过short-step path-following method实现的。值得一提的是Long-step path-following method实现的 $\omega$ 的是 $1$ 。因此，short-step path-following method在很多地方被采用，IPM的迭代次数 $O(\sqrt{n} \log(1/\epsilon ))$ 也是基于此分析出来的。

然后，需要关心IPM每次迭代的计算复杂度。每次迭代的计算复杂度是通过浮点运算次数来分析的。需要注意的是，如果详细去分析每次迭代中每一次浮点运算是比较困难的，且在计算复杂度阶数（order）是没有必要的。为了简化分析，我们只关注主导阶数即可。在IPM的每次迭代中，求解Newton方程是经常耗费大部分计算量的，因此认为求解Newton方程是主导阶数的。注意IPM的Newton方程是一个方程组，具体可参考Boyd的凸优化书籍。所以在一些文献中（例如Pillo&2007&Springer, Th 3.12、Ben-Tal&2001&SIAM, Lecture 6），给出求解Newton方程的复杂度为 $O(mn^3+m^2n^2+m^3+n^3)$ ，其中 $m$ 为约束（含等式和不等式）的个数。根据m与n的相对大小，这个表达式三项中任何一项都有可能成为主导项。当前很多文献对此做了进一步缩放或处理，例如：

(i) Luo&2010&SPM放大为 $O(\max\{m,n\}^4)$ ，有很多文献不考虑 $m$ 或在 $m$ 很小时直接写为 $O(n^4)$ ；

(ii) 优化变量是维度为 $n$ 的向量时，有的文献直接写为 $O(n^3)$ ；

(iii) 在优化变量是维度为 $n$ 的方阵时，有的文献Long&2020&IoT、Sidiropoulos&2006&TSP将方阵变量拉直为维度为 $n^2$ 的向量再处理，则将复杂度的最坏情况写为 $O(n^6)$ ；

(iv) 在优化变量是维度为 $n$ 的方阵时，有的文献Ma&2008&ICASSP为降低复杂度做了努力，则将复杂度变为了 $O(\max\{m,n\}^3)$ 。

当然，也有很多研究致力于根据不同的问题结构去设计减少这个复杂度的算法。

最后，对IPM的迭代次数和每次迭代的计算复杂度相乘便可得到IPM的复杂度order。

2.求解Newton方程的复杂度 $O(mn^3+m^2n^2+m^3+n^3)$ 是怎么算出来的？（基于不等式和等式约束共m个的线性规划）

答：前面的分析是基于路径跟踪法，这个复杂度也是基于该方法分析的。路径跟踪法里面的一个重点就是求解Newton方程得到Newton步径 $\bigtriangleup x_{\text{nt}}$ 以及相关的对偶变量。由于对偶变量的求解复杂度低于求解 $\bigtriangleup x_{\text{nt}}$ 复杂度，就只考虑后者了。

求解 $\bigtriangleup x_{\text{nt}}$ 复杂度依赖于优化变量的维度、约束的类型（等式、不等式or半定）和个数。

case1：对于只有 $m$ 个正半定约束的线性规划SDP，通过求解 $H \bigtriangleup x_{\text{nt}}=-g$ 来得到 $\bigtriangleup x_{\text{nt}}$ ，其中 $H$ 是Hessian矩阵， $g$ 是梯度。

对于半定约束 $S\succeq 0$ 且不存在等式约束的线性问题，路径跟踪法对不等式约束会采用罚函数 $\log\det(S)$ 。基于线性SDP， $H_{ij}=tr(S^{-1}F_iS^{-1}F_j)$ ， $g_i=tc_i+tr(S^{-1}F_i)$ 。计算 $H$ 比 $g$ 更复杂，所以我们接下来只关注 $H$ 的计算。

(i) 计算 $S^{-1}$ 的复杂度为 $O(m^3)$ ；

(ii) 计算 $S^{-1}F_i$ 的复杂度为 $O(m^3)$ ，为了得到 $H$ 需要算 $n$ 项，所以总计 $O(nm^3)$ ；

(iii) 已知 $S^{-1}F_i$ 计算 $tr(S^{-1}F_iS^{-1}F_j)$ ，即将矩阵 $S^{-1}F_i$ 拉直为维度为 $m^2$ 的向量后算内积，因此其复杂度为 $O(m^2)$ ，为了得到 $H$ 需要算 $\frac{n^2}{2}$ 项，忽略常数项所以总计 $O(n^2m^2)$ ；

(iv) 计算 $H^{-1}$ 和 $H^{-1}g$ 的复杂度为 $O(n^3)$ ；

综上，在考虑只有 $m$ 个正半定约束的线性规划情况，便得到了求解Newton方程的复杂度 $O(nm^3+m^2n^2+m^3+n^3)$ ，注意该复杂度是基于半定约束 $S\succeq 0$ 计算出来的。由于 $O(nm^3)$ 主导了 $O(m^3)$ ，有的文献（例如，Ben-Tal&2001&SIAM, Lecture 6、Wang&2014&TSP）直接给出求解Newton方程的复杂度为 $O(nm^3+m^2n^2+n^3)$ 。

case2：对于只有 $m$ 个等式约束的线性规划SDP，参考文献Pillo&2007&Springer以及Boyd的凸优化书籍，直接求解原问题复杂度较高，为了降低复杂度，可以先求解原问题的对偶问题再得到原问题的解。类似于case1，得到Hessian矩阵为 $H_{ij}=tr(A^{-1}A_iA^{-1}A_j)$ 和梯度为 $g_i=tr(A^{-1}A_i)-b_i$ ，其中 $A$ 是 $A_i$ 的线性和，从而分析出求解Newton方程的复杂度 $O(mn^3+m^2n^2+m^3+n^3)\approx O(mn^3+m^2n^2+m^3)$ 。这里值得说明的是，如果能够巧妙地利用原问题的结构直接求解原问题，也有可能得到上述的复杂度。

case3：对于不等式约束和等式约束共 $m$ 个的线性规划，参考文献Ma&2008&ICASSP，可以通过引入松弛变量将不等式约束转化为等式约束，这样分析出求解Newton方程的复杂度和前面的结果相同，即 $O(mn^3+m^2n^2+m^3+n^3)$ 。由于这种case更具有普遍意义，求解Newton方程的复杂度就默认采用该式。

从上述几个case中可以发现:

在 $m< n$ 时，求解Newton方程的复杂度可以简写为 $O(n^3)$ ，一般会假设约束个数给定，因此这种写法最为常见;

在 $m>n$ 时，求解Newton方程的复杂度可以简写为 $O(m^3)$ ;

在 $m\approx n$ 时，求解Newton方程的复杂度可以简写为 $O(\max\{m,n\}^4)$ ;

3.优化变量分别为 $n$ 维向量和 $n$ 维矩阵时，求解Newton方程的复杂度都可以写为 $O(mn^3+m^2n^2+m^3+n^3)$ 吗？

答：显然，当考虑的变量是 $n$ 维向量，求解Newton方程的复杂度为 $O(mn^3+m^2n^2+m^3+n^3)$ ；

然而，文献Ma&2008&ICASSP、Pillo&2007&Springer、Yu&2020&JSAC考虑的是 $n$ 维矩阵，直接求解原问题复杂度较高，为了降低复杂度，可以先求解原问题的对偶问题再得到原问题的解，这样得到的复杂度仍为 $O(mn^3+m^2n^2+m^3+n^3)$ ；当然巧妙地利用原问题的结构直接求解原问题，也有可能得到该复杂度。值得说明的是，由于求解的是对偶问题，迭代降低的仍是对偶残差duality gap，则最大迭代次数仍可写为 $O(\sqrt{n} \log(1/\epsilon ))$ 。

文献Long&2020&IoT将 $n$ 维矩阵考虑为 $n^2$ 维向量，按 $n^2$ 的向量给出最坏情况的复杂度为 $O(n^6)$ 。

4.凸问题是线性的（linear）还是非线性的（nonlinear）会不会影响到前面分析的复杂度表达式？

答：由于IPM始终会求解Newton方程，因此无论凸问题是线性还是非线性，前面复杂度的分析思路仍然可以作为参考。值得说明的是，当凸问题是非线性（无论是目标函数还是约束），计算Hessian矩阵和梯度可能会需要更多的计算成本，因此在计算非线性问题的复杂度时需要考虑这部分。但有时对于复杂的非线性函数，难以去分析Hessian矩阵和梯度的复杂度，因此有的文献为了简便忽略掉这部分的分析，直接在已知 $H$ 和 $g$ 的前提下考虑 $\bigtriangleup x_{\text{nt}}=-H^{-1}g$ 的计算复杂度，其复杂度阶数为 $O(n^3)$ （对应于n维向量变量）和 $O(n^6)$ （对应于n维矩阵变量），然后将该式作为每次迭代复杂度的估算式。此外，根据Wright&1997&SIAM，“路径跟踪法使得解的精度（duality gap）达到 $\epsilon$ 一般需要执行的迭代次数是（至多）为 $O(n^\omega \log(1/\epsilon ))$ ”，这一结论与凸问题是否是线性还是非线性无关。 因此，无论凸问题是线性还是非线性，都不会改变IPM能够实现的最优迭代次数 $O(\sqrt{n} \log(1/\epsilon ))$ 。

5.优化变量在复数域会不会影响到前面分析的复杂度表达式？

答：前面复杂度的分析是浮点运算次数的，因此当优化变量在复数域时，参考Wang&2014&TSP，需要将将复值变量等效为实值变量。具体地，维度为 $n$ 的复向量或复矩阵等效为维度为 $2n$ 的实向量或实矩阵，这样将前面复杂度中的 $n$ 替换为 $2n$ 即可。然而，在表示复杂度order时可以将常数项去掉，因此在很多文献中在进行复数域优化变量分析时仍用前面实数域的复杂度。