凸优化学习笔记：Inner-point methods 内点法

2024年07月29日 2.9k 字大概 14 分钟

参考书目 convex optimization

$$\begin{align}minimize \ \ &f_0(x)\\ subject\ to\ \ &f_i(x)\le 0 \\ &Ax = b\end{align}$$ 其中$f_i:R^n\to R$是凸的且二次可微，$A\in R^{p\times n}, \textbf{rand}\ A=p<n$，假设问题可解，最优点为$x^*$，最优值$p^*=f_0(x^*)$。

内点法通过将牛顿法应用在一系列等式约束上来解最优化问题

首先应将不等式约束化为等式约束

$$\begin{align} minimize\ \ \ &f_0(x)+\Sigma_{i=1}^m I_-(f_i(x))\\ subject\ to\ &Ax=b \end{align}$$ $$I_-(u)=\begin{cases}0,u\le 0\\ \infty,u>0 \end{cases}$$ 当不等式不被满足时，目标函数会变的无穷大，因此当解该最优化问题时自然会满足不等式约束。此时虽然化为了等式约束，但是目标函数不可微。

牛顿法

牛顿步

二阶泰勒展开 $$f(x+\Delta x)\approx f(x)+\nabla f(x)^T\Delta x+\frac{1}{2}\Delta x^T\nabla^2 f(x)\Delta x$$ 是关于$\Delta x$的凸二次函数（默认$\nabla f(x)$是正定的）。为了找到极小值点，应使该函数对$\Delta x$的导数为0，可得当$\Delta x=-\nabla^2 f(x)^{-1}\nabla f(x)^T$时该式取极小值，是下降的方向。

我们将其称为牛顿步

$$\Delta x_{nt}=-\nabla^2f(x)^{-1}\nabla f(x)$$

此处详细解释一下，我们每次走一步牛顿步，都是将当前的二阶泰勒展开公式最小化，这是一个近似，但是每次都会使函数的数值变得越来越小。 ![Pasted_image_20240801144827.png](Pasted_image_Pasted image 20240801144827.png)

牛顿下降量

$$\lambda(x) = (\nabla f(x)^T\nabla^2f(x)^{-1}\nabla f(x))^{1/2}$$ 被称为牛顿下降量。假设$\hat f$为f在x处的二阶近似

$$f(x) -\inf_y\ \hat f(y)=f(x)-\hat f(x+\Delta x_{nt})=\frac{1}{2}\nabla f(x)^T\nabla^2f(x)^{-1}\nabla f(x)=\frac{1}{2}\lambda(x)^2$$ 也就是说，每次牛顿迭代后，二阶估计和目标函数的差值是$\frac{1}{2}\lambda(x)^2$ 将牛顿步代入

$$\lambda = (\nabla x_{nt}^T\nabla^2 f(x)\Delta x_{nt})^{1/2}$$ 换言之，$\lambda$是$\Delta x_{nt}$的Hessian范数

下降法搜索还需要确定步长，来实现$x = x+t\Delta x$的迭代，有以下方式实现步长的确定。

Exact line search

精确线搜索（Exact Line Search），在这种方法中，选择 tt 来最小化沿着射线${x+tΔx∣t≥0}$的函数 f 的值：

$$t=arg⁡min⁡s_{s≥0}f(x+sΔx)$$ 换言之，$\Delta x$决定了搜索的方向，而t则是寻找在这条射线上使函数值最小的自变量的数值。精确线搜索通常用于牛顿法（Newton’s Method）中，特别是在计算搜索方向 Δx的成本较高时。在这种情况下，精确线搜索可以提供比直接计算搜索方向更低的成本。

backtracking line search

大多数的搜索都是不精确的。已经提出了许多不精确线搜索方法。其中一个非常简单且相当有效的方法被称为回溯线搜索（Backtracking Line Search）。它依赖于两个常数 α 和 β，其中 0<α<0.5 和 0<β<1。其算法流程为

given a descent direction ∆x for f at x ∈ dom f, α ∈ (0,0.5), β ∈ (0,1).t := 1. while $f(x + t∆x) > f(x) + αt∇f(x)^T∆x$, t := βt.

也就是说，给定下降方向以及$\alpha,\beta$后，设置t的初始值为1。每次遍历时判断$f(x + t∆x) > f(x) + αt∇f(x)^T∆x$是否成立，成立就更新 t := βt

牛顿法

![Pasted_image_20240801150834.png](Pasted_image_Pasted image 20240801150834.png) 给出起点$x\in \textbf{dom}\ f,\epsilon >0$ 每次迭代牛顿步和牛顿下降量

$$\Delta x_{nt}=-\nabla^2f(x)^{-1}\nabla f(x);\lambda^2=\nabla f(x)^T\nabla^2f(x)^{-1}\nabla f(x)$$ 判断$\lambda^2/2\le \epsilon$，决定是否停止根据backtracking line search选择step size $t$ 更新$x = x+ t\Delta x_nt$。

Logarithmic barriers

为了解决不可微的问题，可以将分段函数转化为对数函数

$$\hat I_-(u)=-(1/t)\log(-u),\textbf{dom}\ {\hat{I}}_-=-R_{++},t > 0$$ 该函数为凸增函数，可微，调节t的值可以调整曲线形状。 ![Pasted_image_20240729113613.png](Pasted_image_Pasted image 20240729113613.png) 此时目标函数可以化为

$$f_0(x)+\Sigma_{i=1}^m -(1/t)\log(-f_i(x))$$ 其中$\phi(x) = - \Sigma_{i=1}^m \log(-f_i(x))$即为Logarithmic barrier。越大的t会使的函数模拟的质量更好，但是使用牛顿法求解$f_0+(1/t)\phi$的难度也会增大（Hessian阵在可行域边界处快速变化）。在目标函数处乘正数t，目标函数可化为

$$tf_0(x)+\phi(x)$$ logarithmic barrier函数的Hessian为：

$$\begin{align}\nabla \phi(x) =\Sigma_{i=1}^m \frac{1}{-f_i(x)}\nabla f_i(x)\\ \nabla^2\phi(x)=\Sigma_{i=1}^m\frac{1}{f_i(x)^2}\nabla f_i(x)\nabla f_i(x)^T+\Sigma_{i=1}^m\frac{1}{-f_i(x)}\nabla^2f_i(x) \end{align}$$

central points

对于不同的t取值，能找到不同的点$x^*(t)$，即为中心点，其集合为中心线。其满足条件 $Ax^*(t)=b,f_i(x^*(t))<0,i=1,…,m$ 且存在$\hat{v}\in R^p$

$$0 = t\nabla f_0(x^*(t))+\nabla \phi(x^*(t))+A^T\hat{v}$$ 其中$\hat{v}$为拉格朗日乘子向量，这一项 $A^Tv$ 的存在是为了满足等式约束 $Ax=b$ 在优化过程中的KKT（Karush-Kuhn-Tucker）条件。

例1:

将log barrier应用到LP问题中

$$\begin{align}minimize\ &c^Tx\\ subject\ to \ &Ax\preceq b\end{align}$$ $\phi(x)=-\Sigma_{i=1}^m\log(b_i-a_i^T x),dom\ \phi = \{x|Ax\prec b\}$ 其中$a_i^T$为A的行。可得对应梯度和的Hessian为

$$\nabla \phi(x)=\Sigma_{i=1}^m \frac{1}{b_i-a_i^Tx}a_i=A^Td,\nabla^2\phi(x)=\Sigma_{i=1}^m\frac{1}{(b_i-a_i^Tx)^2}a_ia_i^T=A^Tdiag(d)^2A$$ $d = 1/(b_i-a_i^Tx)$ 中心点条件为

$$tc+A^Td=0$$

![Pasted_image_20240729140838.png](Pasted_image_Pasted image 20240729140838.png)

Dual points from central path

每个中心点都有一个dual feasible point(相关知识见dual部分内容)，因此可一确定最优值$p^*$的下界。

定义

$$\lambda_i^*(t)=-\frac{1}{tf_i(x^*(t))},\nu^*(t)=\hat\nu /t$$ $λ_i^∗(t)$的值与目标函数 $f_i(x)$ 在当前迭代点 $x^∗(t)$处的负倒数成正比。随着迭代进行，$λi∗(t)$ 的值会逐渐减小，直到满足互补松弛性条件。$ν^*$ 是问题的界限，而 t 是迭代参数，通常初始值设为 1。随着迭代进行，t 的值会逐渐减小，以引导算法逐渐逼近问题的最优解。这些参数的定义是内点法能够有效解决线性规划问题的关键。通过迭代更新这些参数，内点法能够确保算法在每一步都满足KKT条件，从而确保算法能够找到问题的最优解。

此时，将其带入到之前所述的

$$0 = t\nabla f_0(x^*(t))+\nabla \phi(x^*(t))+A^T\hat{v}$$ ，可以得到KKT条件中的

$$\nabla f_0(x^*(t))+\Sigma_{i=1}^m\lambda_i^*(t)\nabla f_i(x^*(t))+A^Tv^*(t)=0$$ 对偶函数是有界的，有

$$g(\lambda^*(t),\nu^*(t))=f_0(x^*(t))+\Sigma_{i=1}^m\lambda_i^*(t)f_i(x^*(t))+\nu^T(Ax-b)=f_0(x^*(t))-m/t$$ 可以得到duality gap是m/t，有

$$f_0(x^*(t))-p^*\le m/t$$

例2：

求例1的对偶问题

$$L(x,\lambda,\nu)=c^Tx+\lambda^T(Ax-b)$$ 因此其下界为

$$g(\lambda,\nu)=\begin{cases}-b\lambda,A^T\lambda+c =0\\ -\infty,A^T\lambda+c \not=0\end{cases}$$ 因此该对偶问题为

$$\begin{align}maximize\ &-b^T\lambda\\ subject\ to\ &A^T\lambda+c=0\end{align}$$ 定义

$$\begin{align} \lambda_i^*(t)&=-\frac{1}{tf_i(x^*(t))},\\ &=-\frac{1}{t(a_i^{T}x-b_i)} \\\nu^*(t)&=\hat\nu /t \end{align}$$ dual目标值为

$$-b^T\lambda^*=c^Tx^*(t)+1/t\Sigma_{i=1}^m 1+(Ax^*(t)-b)^T\lambda^*(t)=c^Tx^*(t)-m/t$$ 这是因为目标x*应该满足$(Ax^*(t)-b)=0$，而前一项因为是倒数相乘变成了1，再求和就是m/t 所以应有

$$f_0(x^*(t))-p^*\le m/t$$

再次说明的是，$p^*$是$f(x^*)$的下界，因此不等式左侧大于0，如果我们的算法迭代时能够不断缩小m/t直到小于允许误差$\epsilon$时，就可以近似认为两数值相等。

barrier method 障碍法

设精度为$\epsilon$，我们使$t = m/\epsilon$，然后使用牛顿法求解问题

$$\begin{align}minimize\ \ &tf_0(x)+\phi(x)\\ subject\ to\ \ &Ax=b \end{align}$$ 即为[Logarithmic barriers](## Logarithmic barriers)部分所转化出的问题。这种方法可以称为无约束最小化方法，因为它允许我们通过解决无约束或线性约束的问题来以保证的准确性解决不等式约束问题。虽然这种方法适用于小问题、良好的起点和中等精度（即ε不太小），但在其他情况下效果不佳。因此，它很少被使用。

实际算法不会直接使t满足精度要求，而是从一个值开始逐步增大。其算法流程为

given strictly feasible x, t := t(0) > 0, μ > 1, tolerance ε > 0. repeat

Centering step.
Compute $x^*(t)$ by minimizing tf(0) + φ, subject to Ax = b, starting at x.

Update. x := $x^⋆(t)$.

Stopping criterion. quit if m/t < ε. 4. Increase t. t := μt.

即在每次迭代时，我们计算新的$x^*$，然后以因子$\mu>0$放大t，直到满足精度要求。在这种求解方式中，我们有两重迭代，外部的centering step迭代和内部的牛顿求解迭代。

$\mu$的选择

当 μ 较小（即接近 1）时，每次外迭代 t 只增加一个小的因子，这意味着初始点（即前一个迭代点 x）对于牛顿过程来说是一个非常良好的起点，因此计算下一个迭代点所需的牛顿步数较少。因此，对于较小的 μ，我们期望每个外迭代只需要少量的牛顿步，但相应地需要大量的外迭代，因为每次外迭代只减少了少量的间隙。在这种情况下，迭代点（实际上，内迭代的迭代点也是）紧密跟随中心路径，这也解释了另一个名称“路径跟随法”。μ 值较大时，外层迭代次数少，内层迭代次数多，迭代点在中心路径上相隔较远；内迭代点会远离中心路径。大约从 3 到 100 左右，两种效果几乎相互抵消，因此牛顿步的总数保持大致恒定。这意味着 μ 的选择并不是特别关键。

$t^{(0)}$的选择

最好能使m/t和$f_0(x^{(0)})-p^*$的量级差不多。

![Pasted_image_20240801101150.png](Pasted_image_Pasted image 20240801101150.png)

示例

$$\begin{align} minimize\ \ &c^Tx\\ subject\ to\ &Ax\preceq b \end{align}$$ $A\in R^{100\times 50}$，可以用一个随机矩阵来进行实验。

#include <Eigen/Dense>

#include <Eigen/Sparse>

#include <Eigen/SparseCholesky>

#include <iostream>

  

using namespace Eigen;

using namespace std;

  

int main() {

// 假设我们有以下线性规划问题

// min c^Tx

// s.t. Ax <= b

// 其中 A 是 5x3 的稀疏矩阵，b 是 5x1 的向量，c 是 3x1 的向量

SparseMatrix<double> A(3, 3);

VectorXd b(3);

VectorXd c(3);

  

// 这里你需要填充 A, b, c 的具体值

// 例如，可以随机生成它们

A.resize(3, 3);

A.coeffRef(0, 0) = 1;

A.coeffRef(0, 1) = 2;

A.coeffRef(0, 2) = 3;

A.coeffRef(1, 0) = 4;

A.coeffRef(1, 1) = 5;

A.coeffRef(1, 2) = 6;

A.coeffRef(2, 0) = 7;

A.coeffRef(2, 1) = 8;

A.coeffRef(2, 2) = 9;

  
  
  

b << 10, 25, 30;

c << 1, 3, 7;

  

// 创建一个SimplicialLDLT对象

SimplicialLDLT<SparseMatrix<double>> solver(A);

  

// 求解问题

VectorXd x = solver.solve(c);

  

// 输出结果

cout << "Optimal solution: " << x.transpose() << endl;

  

return 0;

}