内点法求解线性规划的中心路径与自协调障碍函数

1. 标准化线性规划问题与障碍函数

定义标准线性规划形式：最小化 $c^T x$，满足 $Ax = b$，$x \ge 0$，其中 $x \in \mathbb{R}^n$，$A \in \mathbb{R}^{m \times n}$，$b \in \mathbb{R}^m$，$c \in \mathbb{R}^n$。假设可行域非空且有界。
引入障碍函数：将非负约束 $x \ge 0$ 加入目标函数，构造对数障碍函数：
$$B(x, \mu) = c^T x - \mu \sum_{i=1}^n \ln x_i$$
其中 $\mu > 0$ 是障碍参数。当 $\mu \to 0$ 时，障碍函数的极小点趋近于原问题的最优解。
写出一阶最优性条件：对 $B(x, \mu)$ 在可行域 $Ax = b$ 上求极值，得到 KKT 条件：
$$c - \mu X^{-1} e + A^T y = 0, \quad Ax = b, \quad x > 0$$
其中 $X = \mathrm{diag}(x)$，$e = (1,\dots,1)^T$，$y$ 为拉格朗日乘子。

定义中心路径：对于每个 $\mu > 0$，上述 KKT 条件的解 $(x(\mu), y(\mu))$ 构成一条光滑曲线，称为中心路径。它严格位于可行域内部（$x>0$）。
重写 KKT 条件为等价形式：引入对偶变量 $s = \mu X^{-1} e$，则 $s > 0$ 且满足：
$$A^T y + s = c, \quad Ax = b, \quad X s = \mu e$$
其中 $X s = \mu e$ 表示分量乘积：$x_i s_i = \mu$，$i=1,\dots,n$。
参数化：将 $\mu$ 视为连续参数，当 $\mu$ 从 $\infty$ 递减到 $0$ 时，中心路径从初始内点趋向于原问题和对偶问题的最优解集。
核心性质：中心路径上的点 $(x, y, s)$ 满足互补松弛条件 $x_i s_i = \mu$（严格正），且原始可行 $Ax = b$，对偶可行 $A^T y + s = c$。因此中心路径是原始-对偶内点法的核心跟踪对象。

目的：为了使用牛顿法快速逼近中心路径上的点，需要障碍函数具有良好的凸性——即“自协调”性质。
定义：一个凸函数 $f(x)$ 称为自协调，如果对于任意 $x \in \mathrm{dom} f$ 和任意方向 $h$，其三阶导数被二阶导数控制：
$$|\nabla^3 f(x)[h, h, h]| \le 2 (\nabla^2 f(x)[h, h])^{3/2}$$
直观理解：自协调函数的 Hessian 矩阵变化平缓，保证牛顿法具有全局二次收敛性。
检查对数障碍函数 $-\ln x$ 的自协调性：对一维情形，$f(x) = -\ln x$，计算导数 $f'(x) = -1/x$，$f''(x) = 1/x^2$，$f'''(x) = -2/x^3$。则：
$$|f'''(x) h^3| = 2|x^{-3} h^3| = 2 (x^{-2} h^2)^{3/2} = 2 (f''(x) h^2)^{3/2}$$
满足自协调定义。多个变量的对数障碍函数也是自协调的（Hessian 对角且各分量独立）。
应用：在中心路径的追踪过程中，牛顿步长 $\Delta x$ 由以下方程组确定：
$$\begin{bmatrix} \nabla^2 B & A^T \\ A & 0 \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} = -\begin{bmatrix} \nabla B \\ 0 \end{bmatrix}$$
其中 $\nabla^2 B = \mu X^{-2}$，$\nabla B = c - \mu X^{-1} e$。由于自协调性，牛顿步长可直接取为 1 而不需要线搜索，保证迭代点始终落在可行域内部。

初始化：选取初始 $\mu_0 > 0$ 和初始内点 $x^{(0)} > 0$（满足 $Ax^{(0)} = b$）。通常取 $x^{(0)} = e$，$y^{(0)} = 0$，$s^{(0)} = c$，若不满足对偶可行则调整。
迭代：对于每个 $\mu_k$，执行如下步骤：
1. 计算牛顿方向：求解线性系统：
  $$\begin{bmatrix} \mu_k X^{-2} & A^T \\ A & 0 \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \end{bmatrix} = -\begin{bmatrix} c - \mu_k X^{-1} e \\ 0 \end{bmatrix}$$
  其中 $X = \mathrm{diag}(x^{(k)})$。实际计算中，通常先消去 $\Delta x$ 得到关于 $\Delta y$ 的线性方程组（对称正定）。
2. 步长缩放：由于自协调性，可直接取步长 $\alpha = 1$。但为了严格保持 $x>0$，需进行回溯：寻找最大 $\alpha \in (0,1]$ 使得 $x^{(k)} + \alpha \Delta x > 0$。实践中常用分数边界修正。
3. 更新：$x^{(k+1)} = x^{(k)} + \alpha \Delta x$，$y^{(k+1)} = y^{(k)} + \alpha \Delta y$。
4. 减小 $\mu$：设置 $\mu_{k+1} = \sigma \mu_k$，其中 $\sigma \in (0,1)$（如 $\sigma = 0.1$ 或自适应规则）。
终止：当互补间隙 $x^{(k)T} s^{(k)} < \epsilon$（如 $10^{-8}$）时停止，此时 $x^{(k)}$ 近似为最优解。

直接写出 牛顿步的显式公式：利用 Schur 补，先求解：
$$(A (\mu X^{-2})^{-1} A^T) \Delta y = -A (\mu X^{-2})^{-1} \nabla B$$
再计算 $\Delta x = -(\mu X^{-2})^{-1} (\nabla B + A^T \Delta y)$。
自协调优势：保证每次牛顿步后，点 $(x, s)$ 仍然满足 $x_i s_i \approx \mu$（误差控制在自协调域内），从而不需要额外的校正步骤。这意味着算法复杂度为 $O(\sqrt{n} \log(1/\epsilon))$ 次迭代。
等价形式（原始-对偶）：实际实现中常用原始-对偶对称形式：
$$\begin{bmatrix} 0 & A^T & I \\ A & 0 & 0 \\ S & 0 & X \end{bmatrix} \begin{bmatrix} \Delta x \\ \Delta y \\ \Delta s \end{bmatrix} = -\begin{bmatrix} c - A^T y - s \\ Ax - b \\ X s - \mu e \end{bmatrix}$$
其中 $S = \mathrm{diag}(s)$。该形式利用自协调性质可直接求解，无需显式计算 $X^{-2}$。
数值稳定处理：当 $\mu$ 很小时，$X^{-2}$ 元素极大，造成病态。使用缩放技巧：令 $D = X^{1/2} S^{-1/2}$，将系统化为对称正定形式，避免直接求逆。