本文介绍求解线性规划问题的内点法。它是一个多项式时间算法,在实际应用中效率也很高。尤其是对求解大规模线性规划,一些经验说,内点法比单纯形法更快。此外,内点法还可以被扩展,用来求解凸优化以及非线性规划问题。
考虑线性规划标准问题及其对偶问题:
原始问题(P)
min
c
T
x
s.t.
A
x
=
b
x
≥
0
\begin{aligned} \min~ & c^Tx\\ \text{s.t. } & Ax=b\\ & x\geq 0 \end{aligned}
min s.t. cTxAx=bx≥0
对偶问题(D)
max
b
T
y
s.t.
A
T
y
+
s
=
c
s
≥
0
\begin{aligned} \max~ & b^Ty\\ \text{s.t. } & A^Ty + s = c\\ & s \geq 0 \end{aligned}
max s.t. bTyATy+s=cs≥0
其中
A
∈
R
m
×
n
A\in\mathbb{R}^{m\times n}
A∈Rm×n,
c
,
x
,
s
∈
R
n
c, x, s\in\mathbb{R}^n
c,x,s∈Rn,
b
,
y
∈
R
m
b,y\in\mathbb{R}^m
b,y∈Rm,且矩阵
A
A
A 行满秩。
内点
先定义原始问题和对偶问题的可行域:
F
(
P
)
=
{
x
∈
R
n
:
A
x
=
b
,
x
≥
0
}
F
(
D
)
=
{
(
y
,
s
)
∈
R
m
×
R
n
:
A
T
y
+
s
=
c
,
s
≥
0
}
\begin{aligned} & \mathcal{F}(P) =\{x\in\mathbb{R}^n: Ax=b,x\geq0\}\\ & \mathcal{F}(D) = \{(y,s)\in\mathbb{R}^m\times\mathbb{R}^n: A^Ty+s=c, s\geq 0\} \end{aligned}
F(P)={x∈Rn:Ax=b,x≥0}F(D)={(y,s)∈Rm×Rn:ATy+s=c,s≥0}
接下来定义可行域的 内部(interior):
F
∘
(
P
)
=
{
x
∈
R
n
:
A
x
=
b
,
x
>
0
}
F
∘
(
D
)
=
{
(
y
,
s
)
∈
R
m
×
R
n
:
A
T
y
+
s
=
c
,
s
>
0
}
\begin{aligned} & \mathcal{F}^{\circ}(P) =\{x\in\mathbb{R}^n: Ax=b, x > 0\}\\ & \mathcal{F}^{\circ}(D) = \{(y,s)\in\mathbb{R}^m\times\mathbb{R}^n: A^Ty+s=c, s > 0\} \end{aligned}
F∘(P)={x∈Rn:Ax=b,x>0}F∘(D)={(y,s)∈Rm×Rn:ATy+s=c,s>0}
本文介绍原始对偶内点法(Primal Dual Interior Point Method),它的思路是在可行域的内部
F
∘
(
P
)
×
F
∘
(
D
)
\mathcal{F}^{\circ}(P) \times \mathcal{F}^{\circ}(D)
F∘(P)×F∘(D) 进行迭代,朝着目标优化的方向去逼近最优解。
在可行域内部迭代有什么好处?
给定内部的任意一点,存在它的一个 邻域。如果目标函数是光滑的,那么在内点可导,于是可以用牛顿法解方程。而这个待解的方程就是最优解的充分必要条件(见下文)。
罚函数
给定一个内点,我们想要迭代始终发生在可行域的内部。思路是引入 罚函数(也称障碍函数):
F
(
x
)
=
−
∑
j
=
1
n
ln
(
x
j
)
F(x) = -\sum_{j=1}^n \ln(x_j)
F(x)=−j=1∑nln(xj)
给定
μ
>
0
\mu > 0
μ>0,定义目标函数:
B
μ
(
x
)
=
c
T
x
+
μ
F
(
x
)
B_{\mu} (x) = c^Tx + \mu F(x)
Bμ(x)=cTx+μF(x)
当
x
→
0
x\rightarrow 0
x→0 时,
F
(
x
)
→
∞
F(x) \rightarrow \infty
F(x)→∞,换句话说,
B
μ
(
x
)
B_{\mu}(x)
Bμ(x) 的定义域保证了
x
>
0
x > 0
x>0。
接下来考虑在
F
∘
(
P
)
\mathcal{F}^{\circ}(P)
F∘(P) 中最小化
B
μ
(
x
)
B_{\mu}(x)
Bμ(x),即
min
c
T
x
+
μ
F
(
x
)
s.t.
A
x
=
b
x
>
0
(
PP
)
\begin{aligned} \min~ & c^Tx + \mu F(x)\\ \text{s.t. } & Ax = b\\ & x > 0 \end{aligned} \quad\quad (\text{PP})
min s.t. cTx+μF(x)Ax=bx>0(PP)
当
μ
→
0
\mu\rightarrow 0
μ→0 时,
B
(
μ
)
→
c
T
x
B(\mu)\rightarrow c^Tx
B(μ)→cTx,相当于最小化
c
T
x
c^Tx
cTx;当
μ
→
∞
\mu\rightarrow \infty
μ→∞ 时,相当于最小化
F
(
x
)
F(x)
F(x)。
下面介绍最优解的充要条件。
最优条件
设 F ∘ ( P ) \mathcal{F}^{\circ}(P) F∘(P) 和 F ∘ ( D ) \mathcal{F}^{\circ}(D) F∘(D) 非空。如果 x x x 是上述问题 ( PP ) (\text{PP}) (PP) 的最优解,那么当且仅当存在 ( y , s ) ∈ F ∘ ( D ) (y, s)\in\mathcal{F}^{\circ}(D) (y,s)∈F∘(D),满足如下条件:
-
原始可行: A x = b Ax=b Ax=b
-
对偶可行: A T y + s = c A^Ty + s = c ATy+s=c
-
互补松弛: X S e = μ e XSe = \mu e XSe=μe
其中 X = d i a g ( x ) X = diag(x) X=diag(x), S = d i a g ( s ) S=diag(s) S=diag(s), e = ( 1 , 1 , … , 1 ) T ∈ R n e=(1,1,\dots,1)^T\in\mathbb{R}^n e=(1,1,…,1)T∈Rn。
下面用等式描述上面的条件。
定义函数
F
(
x
,
y
,
s
)
=
[
A
x
−
b
A
T
y
+
s
−
c
X
S
e
−
μ
e
]
F(x,y,s) = \begin{bmatrix} Ax-b\\ A^Ty + s - c\\ XSe-\mu e \end{bmatrix}
F(x,y,s)=⎣⎡Ax−bATy+s−cXSe−μe⎦⎤
x
x
x 是
(PP)
\text{(PP)}
(PP) 的最优解等价于下面的方程有解:
F
(
x
,
y
,
s
)
=
0
F(x,y,s) = 0
F(x,y,s)=0
牛顿法
可以用牛顿法迭代法求解上面的方程 F ( x , y , s ) = 0 F(x,y,s)=0 F(x,y,s)=0。
这里只介绍牛顿法的基本步骤。记
x
k
,
y
k
,
s
k
x^k, y^k, s^k
xk,yk,sk 代表第
k
k
k 步的迭代结果,那么
(
x
k
+
1
,
y
k
+
1
,
s
k
+
1
)
=
(
x
k
,
y
k
,
s
k
)
+
α
k
⋅
(
Δ
x
k
,
Δ
y
k
,
Δ
s
k
)
(x^{k+1}, y^{k+1}, s^{k+1}) = (x^k,y^k, s^k) + \alpha^k \cdot (\Delta x^k, \Delta y^k, \Delta s^k)
(xk+1,yk+1,sk+1)=(xk,yk,sk)+αk⋅(Δxk,Δyk,Δsk)
其中
(
Δ
x
k
,
Δ
y
k
,
Δ
s
k
)
(\Delta x^k, \Delta y^k, \Delta s^k)
(Δxk,Δyk,Δsk) 代表迭代方向,
α
k
∈
(
0
,
1
]
\alpha^k\in (0,1]
αk∈(0,1] 代表步长。
令
J
(
x
,
y
,
s
)
J(x,y,s)
J(x,y,s) 代表
F
(
x
,
y
,
s
)
F(x,y,s)
F(x,y,s) 的雅可比矩阵,即
J
(
x
,
y
,
s
)
=
[
∂
F
1
∂
x
∂
F
1
∂
y
∂
F
1
∂
s
∂
F
2
∂
x
∂
F
2
∂
y
∂
F
2
∂
s
∂
F
3
∂
x
∂
F
3
∂
y
∂
F
3
∂
s
]
=
[
A
0
0
0
A
T
I
A
0
0
S
0
X
]
J(x,y,s) = \begin{bmatrix} \frac{\partial F_1}{\partial x} & \frac{\partial F_1}{\partial y} & \frac{\partial F_1}{\partial s}\\ \frac{\partial F_2}{\partial x} & \frac{\partial F_2}{\partial y} & \frac{\partial F_2}{\partial s}\\ \frac{\partial F_3}{\partial x} & \frac{\partial F_3}{\partial y} & \frac{\partial F_3}{\partial s} \end{bmatrix} = \begin{bmatrix} A & 0 & 0 \\ 0 & A^T & I \\ A & 0 & 0 \\ S & 0 & X \end{bmatrix}
J(x,y,s)=⎣⎢⎡∂x∂F1∂x∂F2∂x∂F3∂y∂F1∂y∂F2∂y∂F3∂s∂F1∂s∂F2∂s∂F3⎦⎥⎤=⎣⎢⎢⎡A0AS0AT000I0X⎦⎥⎥⎤
求解下面的线性方程组:
J
(
x
,
y
,
s
)
[
Δ
x
Δ
y
Δ
s
]
+
F
(
x
,
y
,
s
)
=
0
J(x,y,s)\begin{bmatrix} \Delta x\\ \Delta y\\ \Delta s \end{bmatrix} + F(x,y,s) = 0
J(x,y,s)⎣⎡ΔxΔyΔs⎦⎤+F(x,y,s)=0
即,
[
A
0
0
0
A
T
I
S
0
X
]
[
Δ
x
Δ
y
Δ
s
]
=
−
[
A
x
−
b
A
T
y
+
s
−
c
X
S
e
−
μ
e
]
\begin{bmatrix} A & 0 & 0 \\ 0 & A^T & I \\ S & 0 & X \end{bmatrix}\begin{bmatrix} \Delta x\\ \Delta y\\ \Delta s \end{bmatrix} = -\begin{bmatrix} Ax-b\\ A^Ty + s - c\\ XSe-\mu e \end{bmatrix}
⎣⎡A0S0AT00IX⎦⎤⎣⎡ΔxΔyΔs⎦⎤=−⎣⎡Ax−bATy+s−cXSe−μe⎦⎤
可以得到牛顿法的迭代方向。
由于
x
,
y
x, y
x,y 是可行解,上式可以写成
[
A
0
0
0
A
T
I
S
0
X
]
[
Δ
x
Δ
y
Δ
s
]
=
−
[
0
0
X
S
e
−
μ
e
]
\begin{bmatrix} A & 0 & 0 \\ 0 & A^T & I \\ S & 0 & X \end{bmatrix}\begin{bmatrix} \Delta x\\ \Delta y\\ \Delta s \end{bmatrix} = -\begin{bmatrix} 0\\ 0\\ XSe-\mu e \end{bmatrix}
⎣⎡A0S0AT00IX⎦⎤⎣⎡ΔxΔyΔs⎦⎤=−⎣⎡00XSe−μe⎦⎤
中心路径
回顾问题
(PP)
\text{(PP)}
(PP),它在
F
∘
(
P
)
\mathcal{F}^{\circ}(P)
F∘(P) 中最小化最小化
B
μ
(
x
)
B_{\mu}(x)
Bμ(x),即
min
c
T
x
+
μ
F
(
x
)
s.t.
A
x
=
b
x
>
0
(
PP
)
\begin{aligned} \min~ & c^Tx + \mu F(x)\\ \text{s.t. } & Ax = b\\ & x > 0 \end{aligned} \quad\quad (\text{PP})
min s.t. cTx+μF(x)Ax=bx>0(PP)
我们已知它最优解的充要条件,而且可以用牛顿法求它的最优解。
但是,给定 μ > 0 \mu > 0 μ>0, (PP) \text{(PP)} (PP) 的最优解并不是原问题 (P) \text{(P)} (P) 的最优解。因此,迭代的目标不是求解 (PP) \text{(PP)} (PP) ,而是去逼近原问题 (P) \text{(P)} (P) 的最优解。逼近的方式就是在迭代中不断减小 μ \mu μ 值,然后更新下降方向;当 μ \mu μ 足够小时,得到的解与原问题的最优解足够接近。
具体来说,给定初始值
μ
0
>
0
\mu^0 > 0
μ0>0 和初始点
x
0
(
μ
0
)
,
y
0
(
μ
0
)
,
s
0
(
μ
0
)
∈
F
∘
(
P
)
×
F
∘
(
D
)
x^0(\mu^0), y^0(\mu^0), s^0(\mu^0) \in \mathcal{F}^{\circ}(P)\times \mathcal{F}^{\circ}(D)
x0(μ0),y0(μ0),s0(μ0)∈F∘(P)×F∘(D) ,迭代过程得到如下的点列:
x
0
(
μ
0
)
,
y
0
(
μ
0
)
,
s
0
(
μ
0
)
,
x
1
(
μ
1
)
,
y
1
(
μ
1
)
,
s
1
(
μ
1
)
,
⋮
x
k
(
μ
k
)
,
y
k
(
μ
k
)
,
s
k
(
μ
k
)
x^0(\mu^0), y^0(\mu^0), s^0(\mu^0),\\ x^1(\mu^1), y^1(\mu^1), s^1(\mu^1),\\ \vdots\\ x^k(\mu^k), y^k(\mu^k), s^k(\mu^k)\\
x0(μ0),y0(μ0),s0(μ0),x1(μ1),y1(μ1),s1(μ1),⋮xk(μk),yk(μk),sk(μk)
其中
μ
0
>
μ
1
>
⋯
>
μ
k
\mu^0 > \mu^1 > \dots > \mu^k
μ0>μ1>⋯>μk。当
μ
k
<
ϵ
\mu^k < \epsilon
μk<ϵ 时,迭代停止,其中
ϵ
>
0
\epsilon > 0
ϵ>0 代表逼近的精度。
我们把 { x ( μ ) , y ( μ ) , s ( μ ) : μ > 0 } \{x(\mu),y(\mu),s(\mu): \mu > 0\} {x(μ),y(μ),s(μ):μ>0} 称为 中心路径(Central Path)。迭代的过程就是沿着中心路径,令 μ \mu μ 值不断减小的过程,这样的算法也称为 路径跟踪(Path Following)。
内点法(理论版)
接下来介绍具体的迭代步骤。
原始对偶内点法 (Primal-Dual Interior Point Method)
第一步,初始化 μ 0 , x 0 , y 0 , s 0 \mu^0,x^0,y^0,s^0 μ0,x0,y0,s0。注意: μ , x , s > 0 \mu, x, s >0 μ,x,s>0 且 x 0 , y 0 , s 0 x^0,y^0,s^0 x0,y0,s0 可行。
第二步,判断误差:如果 μ < ϵ \mu < \epsilon μ<ϵ 则停止。
第三步,计算牛顿方向:
[
A
0
0
0
A
T
I
S
0
X
]
[
Δ
x
Δ
y
Δ
s
]
=
−
[
0
0
X
S
e
−
μ
e
]
\begin{bmatrix} A & 0 & 0 \\ 0 & A^T & I \\ S & 0 & X \end{bmatrix}\begin{bmatrix} \Delta x\\ \Delta y\\ \Delta s \end{bmatrix} = -\begin{bmatrix} 0\\ 0\\ XSe-\mu e \end{bmatrix}
⎣⎡A0S0AT00IX⎦⎤⎣⎡ΔxΔyΔs⎦⎤=−⎣⎡00XSe−μe⎦⎤
第四步,更新
x
,
y
,
s
x, y, s
x,y,s:
(
x
,
y
,
s
)
←
(
x
,
y
,
s
)
+
(
Δ
x
,
Δ
y
,
Δ
s
)
(x,y,s) \leftarrow (x,y,s) + (\Delta x, \Delta y, \Delta s)
(x,y,s)←(x,y,s)+(Δx,Δy,Δs)
第五步,减小
μ
\mu
μ 值,令
μ
←
σ
⋅
μ
\mu \leftarrow \sigma \cdot \mu
μ←σ⋅μ,其中
σ
∈
[
0
,
1
]
\sigma \in [0,1]
σ∈[0,1], 然后执行第二步。
这里还有几个问题需要回答:
-
μ \mu μ 和 σ \sigma σ 值如何选取?
-
第四步更新 x , y , s x,y,s x,y,s 之后,它还是可行解吗?
-
如何计算初始可行解 x 0 , y 0 , s 0 x^0,y^0,s^0 x0,y0,s0 ?
先看第一个问题。
先计算原问题
(P)
\text{(P)}
(P) 和对偶问题
(D)
\text{(D)}
(D) 的对偶间隙 (Duality Gap):
c
T
x
−
b
T
y
=
(
A
T
y
+
s
)
T
x
−
b
T
y
=
s
T
x
c^Tx - b^Ty = (A^Ty + s)^Tx - b^Ty = s^Tx
cTx−bTy=(ATy+s)Tx−bTy=sTx
令
μ
=
s
T
x
n
\mu = \frac{s^Tx}{n}
μ=nsTx,当
μ
<
ϵ
\mu < \epsilon
μ<ϵ 时,我们有
c
T
x
−
b
T
y
≤
n
ϵ
c^Tx - b^Ty \leq n \epsilon
cTx−bTy≤nϵ,即目标函数值与最优值当误差不超过
n
ϵ
n\epsilon
nϵ。
此外,令 σ = 1 − 0.4 n \sigma = 1- \frac{0.4}{\sqrt{n}} σ=1−n0.4,可以证明最多 O ( n ln ( 1 / ϵ ) ) O(\sqrt{n}\ln(1/\epsilon)) O(nln(1/ϵ)) 次迭代,就可以使得 μ < ϵ \mu < \epsilon μ<ϵ。
再看第二个问题。
只需要验证
x
+
Δ
x
x+\Delta x
x+Δx 与
(
y
+
Δ
y
,
s
+
Δ
s
)
(y+\Delta y, s+\Delta s)
(y+Δy,s+Δs) 是否满足原始问题和对偶问题的约束:
A
(
x
+
Δ
x
)
=
b
+
A
Δ
x
=
b
A
T
(
y
+
Δ
y
)
+
(
s
+
Δ
s
)
=
c
+
A
T
Δ
y
+
Δ
s
=
c
\begin{aligned} & A(x+\Delta x) = b + A\Delta x = b\\ & A^T(y +\Delta y) + (s + \Delta s) =c + A^T\Delta y + \Delta s = c \end{aligned}
A(x+Δx)=b+AΔx=bAT(y+Δy)+(s+Δs)=c+ATΔy+Δs=c
因此答案是肯定的。
关于第三个问题,计算初始可行解并不容易,这里不做介绍。原因是在实际应用中,可以通过逼近的方式满足可行性,因而初始解只要保证 x , s > 0 x,s>0 x,s>0 即可。
内点法(实践版)
编程实现内点法时,一般不会严格按照上面的理论版,原因是效率不高。
下面介绍一个实践版本。
第一步,初始化 x 0 > 0 , s 0 > 0 , y 0 x^0 > 0,s^0 >0, y^0 x0>0,s0>0,y0,令 μ = s T x n \mu = \frac{s^Tx}{n} μ=nsTx。
第二步,判断误差:如果
μ
<
ϵ
\mu < \epsilon
μ<ϵ 则停止,其中
μ
=
max
{
∣
∣
A
x
−
b
∣
∣
,
∣
∣
A
T
y
+
s
−
c
∣
∣
,
∣
∣
s
T
x
∣
∣
}
\mu = \max \{||Ax-b||, ||A^Ty+s-c||, ||s^Tx||\}
μ=max{∣∣Ax−b∣∣,∣∣ATy+s−c∣∣,∣∣sTx∣∣}
第三步,计算牛顿方向:
[
A
0
0
0
A
T
I
S
0
X
]
[
Δ
x
Δ
y
Δ
s
]
=
−
[
A
x
−
b
A
T
y
+
s
−
c
X
S
e
−
μ
e
]
\begin{bmatrix} A & 0 & 0 \\ 0 & A^T & I \\ S & 0 & X \end{bmatrix}\begin{bmatrix} \Delta x\\ \Delta y\\ \Delta s \end{bmatrix} = -\begin{bmatrix} Ax-b\\ A^Ty + s - c\\ XSe-\mu e \end{bmatrix}
⎣⎡A0S0AT00IX⎦⎤⎣⎡ΔxΔyΔs⎦⎤=−⎣⎡Ax−bATy+s−cXSe−μe⎦⎤
第四步,更新
x
,
y
,
s
x, y, s
x,y,s:
x
←
x
+
α
P
⋅
Δ
x
(
y
,
s
)
←
(
y
,
s
)
+
α
D
⋅
(
Δ
y
,
Δ
s
)
\begin{aligned} & x \leftarrow x + \alpha_P \cdot \Delta x \\ & (y,s) \leftarrow (y,s) + \alpha_D \cdot (\Delta y, \Delta s) \end{aligned}
x←x+αP⋅Δx(y,s)←(y,s)+αD⋅(Δy,Δs)
其中
α
P
,
α
D
\alpha_P, \alpha_D
αP,αD 使得
x
,
s
>
0
x,s>0
x,s>0,即
α
P
=
min
{
1
,
min
Δ
x
j
<
0
{
x
j
−
Δ
x
j
}
}
α
D
=
min
{
1
,
min
Δ
s
j
<
0
{
s
j
−
Δ
s
j
}
}
\alpha_P = \min \left\{1, \min_{\Delta x_j < 0} \left\{\frac{x_j}{-\Delta x_j}\right\}\right\}\\ \alpha_D = \min \left\{1, \min_{\Delta s_j < 0} \left\{\frac{s_j}{-\Delta s_j}\right\}\right\}
αP=min{1,Δxj<0min{−Δxjxj}}αD=min{1,Δsj<0min{−Δsjsj}}
第五步,减小
μ
\mu
μ 值,令
μ
←
σ
⋅
μ
\mu \leftarrow \sigma \cdot \mu
μ←σ⋅μ,其中
σ
=
0.1
\sigma = 0.1
σ=0.1 (经验值), 然后执行第二步。
注意两点变化:
第一,初始解不要求可行,原因是可以通过迭代的方式保障可行性。它通过 μ \mu μ 的选择来实现,它是三种误差的最大值,即原问题的可行性误差 ∣ ∣ A x − b ∣ ∣ ||Ax-b|| ∣∣Ax−b∣∣, 对偶问题的可行性误差 ∣ ∣ A T y + s − c ∣ ∣ ||A^Ty + s-c|| ∣∣ATy+s−c∣∣,以及目标函数误差 s T x s^Tx sTx。 当 μ < ϵ \mu < \epsilon μ<ϵ 且 ϵ \epsilon ϵ 精度足够小时,原问题和对偶问题可行,并且目标函数值达到最优。
第二,由于迭代过程中,原始解和对偶解可能不满足可行性,因此更新 x , y , s x,y,s x,y,s 的步长 α \alpha α 需要保证 x , s > 0 x,s>0 x,s>0。其中 x x x 与 ( y , s ) (y,s) (y,s) 也可以采用相同的步长,例如 α = min ( α D , α P ) \alpha = \min (\alpha_D, \alpha_P) α=min(αD,αP)。
算法实现
下面用 Python 实现内点法。
先定义问题的输入和输出。
import numpy as np
class IPMlp(object):
""" A Primal Dual Interior Point Method for linear programming.
说明:
1. 最小化问题
2. 行满秩
"""
def __init__(self, c, A, b):
"""
:param c: n * 1 vector
:param A: m * n matrix
:param b: m * 1 vector
"""
# 输入
self._c = np.array(c)
self._A = np.array(A)
self._b = np.array(b)
# 输出
self._obj = None # objective function value
self._x = None # primal solution
self._y = None # dual solution
# 辅助变量
self._m = len(self._b)
self._n = len(self._c)
self._s = None # dual slack variable
self._alpha = None
self._sigma = None
self._mu = None
self._iter_num = 0
self._status = None
接下来是算法实现的思路。
class IPMlp(object):
# 其它函数省略
# ...
def solve(self):
self._init() # initialization
# 判断停止条件
while not self._is_stop_criterion_satisfied():
# 解线性方程组,得到牛顿下降方向
dx, dy, ds = self._solve_linear_system()
if dx is None:
self._status = 'UNBOUNDED or INFEASIBLE'
break
# 计算步长,保证 x, s > 0
self._alpha = self._get_alpha(dx, ds)
# 更新原始解和对偶解
# 这里采用了统一的步长
self._x += self._alpha * dx
self._y += self._alpha * dy
self._s += self._alpha * ds
# 减小 mu 值
self._mu = (self._x @ self._s / self._n) * self._sigma
# 更新目标函数
self._obj = self._c @ self._x
更多细节请参考源码:完整代码
参考文献
[1]. David P. Williamson. ORIE 6300 Mathematical Programming I (lecture notes). 2014.
[2]. R. M. Freund and J. Vera. Interior-Point Methods for Linear Optimization. 2014.