跳到正文
格致开物MATHWIKI

内点法

内点法(interior-point method)是一类在约束内部或通过原始与对偶变量的内部近似路径求解优化问题的方法。以下用凸问题的对数障碍法说明其基本机制:暂时避开不等式边界,求解一串较容易处理的光滑问题,再让障碍权重逐渐减小。其他内点实现并不都要求每个原始迭代点严格可行。

区间内的障碍与中心点

考虑 min⁡x,约束为 0≤x≤1。最优点是左端零,但为避免直接碰到边界,引入参数 μ>0,在 0<x<1 上最小化 ψμ(x)=x−μlog⁡x−μlog⁡(1−x). 靠近任一端点,障碍项趋向正无穷。其导数为 1−μ/x+μ/(1−x),二阶导数为 μ/x2+μ/(1−x)2>0,故导数零点是唯一最小点。乘以正数 x(1−x),得到 x2−(1+2μ)x+μ=0,xμ=1+2μ−1+4μ22. 另一根大于一,不在定义域。为避免浮点相减损失精度,还可将小根有理化为 2μ/(1+2μ+1+4μ2)。

不同障碍权重下区间内目标的谷底
障碍在两端发散;缩小权重使中心点靠近左端最优点。

取 μ=0.5,0.1,0.01,中心点分别约为 0.292893,0.090098,0.009900。图中障碍目标的谷底向左移动;每个有限参数的解仍在内部,极限才到达原问题的边界最优点。

凸模型与可计算的误差证书

一般考虑 min⁡f0(x),约束 fi(x)≤0(i=1,…,m)、Ax=b,其中各 fi 凸且在所用开域上二阶连续可微,等式为仿射。假定存在严格可行点,且当前障碍子问题达到最小值 xμ。子问题为 minAx=b f0(x)−μ∑i=1mlog⁡[−fi(x)],fi(x)<0. 假设中心点存在相应等式乘子,定义 λi=μ/[−fi(xμ)]>0。障碍的一阶条件便成为 ∇f0(xμ)+∑iλi∇fi(xμ)+A𝖳ν=0. 由凸性,xμ 全局最小化拉格朗日函数 L(x,λ,ν)=f0(x)+∑iλifi(x)+ν𝖳(Ax−b)。其下确界因此是 g(λ,ν)=L(xμ,λ,ν)=f0(xμ)−mμ. 对任何原问题可行点,非负乘子使 L≤f0,所以这个对偶值是原最优值 f∗ 的下界。又因中心点可行,得到 0≤f0(xμ)−f∗≤mμ. 这给出了完整的障碍误差论证。它要求精确中心条件;数值上尚未解完子问题时,不能不加残差修正就使用精确的 mμ 证书。教材也常用 t=1/μ 写成 m/t。[1]

中心点实际误差与可证明对偶差距
横轴依次对应 μ=0.5、0.1、0.01、0.001;实际误差小于证书 2μ。

在区间例中,两约束为 −x≤0、x−1≤0,相应乘子是 λ1=μ/xμ、λ2=μ/(1−xμ)。驻点条件为 1−λ1+λ2=0,对偶值为 −λ2;差距 xμ+λ2=2μ 可逐项复算。实际误差 xμ 通常小于这个上界。

路径跟踪与严格可行性的边界

障碍法从较大的 μ 开始,解当前子问题,再缩小参数并以前一中心点作为初值。牛顿更新还须控制步长,使 fi(x)<0 不被破坏。接近边界时障碍曲率变大,直接从远处求一个极小参数的中心点可能很困难。

变式将区间改为 [0,2]。一般区间 [0,b] 的小根为 [b+2μ−b2+4μ2]/2;取 b=2,μ=0.1 得 x=1.1−1.01≈0.095012。两约束仍给出 2μ 的对偶差距。

若约束同时为 x≤0 与 −x≤0,原问题有可行点零,但严格内部为空,对数障碍没有定义域。需要先消去固定变量、作可行性阶段或选择合适原始对偶方法。仅凭“原问题可行”不能启动这里的严格可行障碍法。

可运行的中心路径计算

下例只用标准库,以二分法解一维障碍的一阶条件,返回中心点及实际误差、对偶差距,不改写参数列表。每个正参数求解七十次区间折半,区间宽度至多 b2−70;计算量为 O(Klog⁡(b/ε)),返回 K 个参数的历史需 O(K)。一般内点法的主要成本是牛顿线性系统,不是这里的一维折半。

Python 3
def centers(parameters, b=1.0):
    if b <= 0 or any(mu <= 0 for mu in parameters):
        raise ValueError("positive interval and barrier weights required")
    rows = []
    for mu in parameters:
        lo, hi = 0.0, float(b)
        for _ in range(70):
            x = (lo+hi)/2
            derivative = 1-mu/x+mu/(b-x)
            if derivative < 0:
                lo = x
            else:
                hi = x
        x = (lo+hi)/2
        dual = -b*mu/(b-x)
        rows.append((mu,x,x-dual))
    return rows

for row in centers([0.5,0.1,0.01]):
    print(tuple(round(v,6) for v in row))
print("variant", centers([0.1],2)[0])
assert all(abs(gap-2*mu) < 1e-12 for mu,x,gap in centers([.5,.1,.01]))

参考文献

  1. ↑ Boyd, S. and Vandenberghe, L., Convex Optimization,§11.2.1—11.2.2及§11.3。