线性规划(一)——单纯形法

问题介绍

单纯形法(simplex method)是求解线性规划问题一种通用算法，在实际生产生活中有广泛的应用。有些教材在介绍单纯形法时使用了复杂的矩阵和下标运算，使得算法的脉络难以把握。本文试图从几何意义和简明的代数运算出发，解释单纯形法的运行逻辑。

标准形

• 标准形：线性规划问题的标准形形如：
  (a) 目标函数：\displaystyle\min~z=\sum_{j=1}^nc_jx_j
  (b) 线性约束： \displaystyle\sum_{j=1}^na_{ij}x_j=b_i~~(1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}m)
  (c) 非负约束： x_j\geqslant0~~(1\hspace{-1pt}\leqslant \hspace{-1pt}j\hspace{-1pt}\leqslant \hspace{-1pt}n)
  其中 b_i\geqslant0~~(1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}m).

标准形更加常用的是以下矩阵形式：

• 矩阵形式：矩阵形式的标准形形如：
  (a) 目标函数：\min z=c^\mathrm{T}x
  (b) 线性约束：Ax=b
  (c) 非负约束：x\geqslant0
  其中 A\in\mathbb{R}^{m\times n},~c\in\mathbb{R}^n,~b\in\mathbb{R^m}，且 b\geqslant0.

线性约束条件 Ax=b 由 m 个线性方程组成，为了理论上的方便，一般要求这 m 个方程线性无关，即 \mathrm{rank}~A=m . 以下均假设该条件成立。

---

下面来讨论标准形的几何意义。众所周知，下面的方程
a_1x_1+\cdots+a_nx_n=b 表示 n 维欧氏空间 \mathbb{R}^n 中的一个 n-1 维超平面，m 个线性无关方程构成的方程组 Ax=b 则表示 m 个超平面的交，也就是一个 n-m 维超平面。举例来说，若 (n,m)=(3,1), 则线性约束条件转化为 \mathbb{R}^3 中的一个平面(如下图所示)。

(n,m)=(3,1) 的线性规划问题

在 x\geqslant0 的非负约束条件下，n-m 维超平面 S 被坐标轴截为 n-m 维凸图形 S_0, 从而线性函数 z=c^\mathrm{T}x 的最小值一定在 S_0 的顶点上取。确切地说，这些顶点是 S 与某 m 条坐标轴张成的 m 维子空间的交点。在 (n,m)=(3,1) 的条件下，平面 S 被坐标轴截为三角形 S_0, 并且 S_0 的顶点恰好是 S 与三个一维坐标轴的交点。

S_0 的顶点由 S 和 m 维坐标轴子空间相交而成，而 m 条坐标轴又选自于 \mathbb{R}^n 的 n 条坐标轴，因此 S_0 的顶点数不超过 C_n^m. (实际可能会少的多)

---

回到矩阵形式。记原始的线性规划问题为 \mathcal{P}，则超平面 S 表示 Ax=b 的解，超平面凸图形 S_0 表示 \mathcal{P} 的可行域。由于线性函数 c^\mathrm{T}x 的最小值在 S_0 的顶点取到，因此 S_0 的顶点处的性质具有“重要的研究意义”。若 P 是 S_0 的一个顶点，则 P 的 n 个分量中只有坐标轴对应的 m 个分量可能取非零值，其余 n-m 个分量均必取 0。

对于“坐标轴”“顶点”的初步印象自然引出了“基向量组”“基本解”的概念：

• 基向量组： 设 A\in\mathbb{R}^{m\times n}，\mathrm{rank}~A=m，且 A=(P_1,\cdots,P_n)。 A 的 m 个线性无关的列向量 (P_{i_1},\cdots,P_{i_m}) 称为 A 的一个基向量组。
• 基变量： 线性方程组 Ax=b 的解 x 在 (i_1,\cdots,i_m) 这 m 个坐标上的分量称为基变量(或基分量)。对于给定的基向量组而言，基变量共有 m 个，剩下的 n-m 个分量则称为非基变量(或非基分量)。
• 基本解：设 B=(P_{i_1},\cdots,P_{i_m}) 是 A 的一个基向量组，则 Ax=b 有唯一的解 x, 使得它在非基分量上的取值均为 0，即 x 满足 x_j=0~~(j\neq i_1,\cdots,i_m)。 这一解 x 称为基向量组 B 对应的基本解。

基向量组、基变量、基本解

上述对于基本解的定义是合理的。基向量组 B=(P_{i_1},\cdots,P_{i_m}) 也可以理解为列向量的水平叠置，因而 B 也被称为基向量矩阵。在下面的叙述当中，基向量组和基向量矩阵的概念是等同的，不再加以区分。B 所对应的基本解仅可能在 i_1,\cdots,i_m 这 m 个分量上取非零值，将这 m 个分量构成的向量记为 x_B~(x_B\in\mathbb{R}^m)，则原 n 元线性约束条件可简化为 m 元的情形：
Ax=b\Rightarrow Bx_B=b由于 B 是 m 阶可逆矩阵，故 x_B=B^{-1}b。这说明基向量组 B 对应的基本解的确是唯一的。

更进一步，如果基向量组 B 对应的基本解 x 满足 x\geqslant0，则 x 位于可行域当中，因而 B 称为可行基向量组，x 称为基本可行解。从几何意义来看，基本可行解就是可行域 S_0 的顶点。如果 \mathcal{P} 的最优解(使目标函数取最小值的解)存在，则一定存在一个基本可行解是最优的，换言之，最优解(目标函数最小值)一定能在基本可行解(S_0 的顶点)上取到。

如上图所示，P_1,P_2,P_3 三个顶点都分别只有一个分量非零，因而它们都是基本解。由于它们都在各自坐标轴的正半轴上，因此都是基本可行解。如果移动 S 使得 S 与 x_1 轴的交点 P_1 在负半轴上，那么 P_1 就不是基本可行解了。

---

以上的讨论都是关于线性规划问题 \mathcal{P} 的静态性质。下面所要讨论的单纯形法，则是真正意义上给出了求出最优解的算法。尽管 \mathcal{P} 的几何意义对于单纯形法而言没有实质性的作用，但它提供了理解诸如“基本解”等概念的一种有效途径。

单纯形法的核心步骤是：(a) 最优性检验; (b) 基变换。简单来说，这两个步骤的主要作用是：

• 最优性检验： 给定可行基向量组 B，判断 B 对应的基本(可行)解是否是最优的。
• 基变换： 如果不是，构造另一个可行基向量组 B'，使得目标函数不增。

交替进行上述2个步骤，直至找到最优解或证明最优解不存在。

最优性检验

给定可行基向量组 B=(P_{\pi(1)},\cdots,P_{\pi(m)})，为了讨论基本解的最优性，需要计算目标函数 c^{\mathrm{T}}x 在基本解附近的值。

设 B\in\mathbb{R}^{m\times m} 为 A 的 \pi 分量部分构成的基向量组，而 N\in\mathbb{R}^{m\times(n-m)} 为 A 的非 \pi 分量部分构成的非基向量组。设 x\in S_0 是可行域内任一点，x 的 \pi 分量部分记为 x_B\in\mathbb{R}^m，非 \pi 分量部分记为 x_N\in\mathbb{R}^{n-m}，c_B,c_N 的意义类似。于是
Ax=b\Rightarrow (B~N) \begin{pmatrix} x_B\\ x_N \end{pmatrix}=b \Rightarrow x_B=B^{-1}(b-Nx_N) 上式表明可行解的基分量部分 x_B 可由非基分量部分 x_N 唯一线性表示。这一结论从几何意义上是容易理解的：S 是 n-m 维超平面，含有 n-m 个自由变量，因此 S 上的任一点的坐标可以仅由其中的 n-m 个线性表示。

于是，z=c^{\mathrm{T}}x 也应当仅由 n-m 个分量表示：
\begin{aligned} z&=c^{\mathrm{T}}x\\ &=(c_B^{\mathrm{T}}~c_N^\mathrm{T})\begin{pmatrix}x_B\\ x_N\end{pmatrix}\\ &=c_B^{\mathrm{T}}x_B+c_N^{\mathrm{T}}x_N\\ &=c_B^{\mathrm{T}}B^{-1}(b-Nx_N)+c_N^{\mathrm{T}}x_N\\ &=c_B^{\mathrm{T}}B^{-1}b+(c_N^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}N)x_N \end{aligned}上式表明 z 可仅由 n-m 个非基分量的值 x_N 线性表示。

注意到目标函数在 B 对应的基本可行解 x^{(0)} 处的取值是 z_0=c_B^{\mathrm{T}}B^{-1}b，有
\begin{aligned} z&=z_0+(c_N^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}N)x_N\\ &=z_0+(c_B^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}B)x_B+(c_N^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}N)x_N\\ &=z_0+(c^{\mathrm{T}}-c_B^{\mathrm{T}}B^{-1}A)x\\ :&=z_0+\lambda^{\mathrm{T}} x \end{aligned}其中 \lambda^{T}=c^{\mathrm{T}}-c_B^{T}B^{-1}A\in\mathbb{R}^n 为检验向量，其各个分量称为检验数。尽管我们最终将 z 重新写成了 x 的线性函数，但实际上 x_B 的系数(m 个检验数)均为 0，因此 z 一共只有 n-m 个自由变量(即 x_N)。基于此，目标函数 z=z_0+\lambda^{\mathrm{T}}x 的这种形式也被称为简化形式。

在可行域的限制之下，x\geqslant0。如果 \lambda 的每一个分量均非负(检验数均非负)，则 z\geqslant z_0 恒成立，且 x=x^{(0)} 时等号成立(注意 x_N^{(0)}=0)。\lambda 各个分量的正负性成为了判别最优解的基本依据。

在介绍下面的定理前，先引入如下记号：\alpha=(\alpha_{ij})_{m\times n}=B^{-1}A, P_j'=B^{-1}P_j=(\alpha_{1j},\cdots,\alpha_{mj})^{\mathrm{T}}, (1\hspace{-1pt}\leqslant \hspace{-1pt}j\hspace{-1pt}\leqslant \hspace{-1pt}n), \beta=B^{-1}b.

定理 (最优性检验)：
设可行基向量组 B 对应的基本解为 x^{(0)}，目标函数的简化形式为 z=z_0+\lambda^{\mathrm{T}}x。

• 若所有的检验数 \lambda_k\geqslant 0，则 x^{(0)} 是最优解，z_0 是目标函数的最小值。
• 若存在某检验数 \lambda_k<0，且 \alpha_{ik}\leqslant0~~(1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}m)，则无最优解。

证明：第一部分的证明已经在上面给出。考虑第二部分：
因为 \lambda_k\neq0，所以 x_k 是非基分量。取的 x 的 k 分量 x_k=m>0，而其它非基分量均为 0。此时 x_N=(0,\cdots,m,\cdots,0)^{\mathrm{T}}，代入 x_B 的表达式可得
x_B=B^{-1}(b-Nx_N)=\beta-B^{-1}P_km=\beta-P_k'm 由于 P_k'\leqslant0，故 x_B\geqslant0，从而 x 的各个分量非负，是 \mathcal{P} 的可行解。但此时 z=z_0+\lambda_km，m 可以取充分大的正数，因此目标函数无下界，不存在最优解。

---

总而言之，最优性检验利用目标函数的简化形式、检验数的正负来判断当前基本可行解的最优性。在最优性检验定理当中，还有一种情况未被涉及：存在某检验数 \lambda_k<0，且存在下标 i 使得 \alpha_{ik}>0。这正是下面的基变换所要解决的问题。

基变换

在最优性检验当中，我们利用可行基向量组 B 和对应的基本解 x^{(0)} 得到来目标函数的简化形式 z=z_0+\lambda^{\mathrm{T}}x 和检验向量 \lambda。在最优性检验定理的未涉及的情形中，存在检验数 \lambda_k<0，并且 \alpha 的第 k 列 P_k' 至少有一个分量为正。于是，我们希望替换基向量组中某一个基向量，得到一个新的可行基向量组，并且保证目标函数不增(然后重新评估新基对应的基本解的最优性)。这一步骤就是所谓的基变换。

在进一步叙述基变换的具体操作之前，我们需要补充一些线性代数的知识。事实上，我更希望基变换的计算最后是由“纯粹的”向量和矩阵运算完成的。我喜欢赋予向量和矩阵鲜明的代数意义，通过简洁而又富有想象力的符号运算来获取结果，而讨厌与令人头疼的下标和臃肿的 n 阶矩阵展开式打交道。在我使用的《算法分析与设计》一书中，基变换的核心步骤就是以下标运算和矩阵展开写就的，矩阵的一些优美性质被掩盖了。因此，我决定使用更加富有技巧和表现力的矩阵运算来展现这一部分。

---

Gauss-Jordan 矩阵

设 n 为给定的正整数。记 e_l 是第 l 个分量为 1，其余分量为 0 的 n 维列向量，即
e_l=(0,\cdots,\underbrace{1}_{l\mathrm{-th}},\cdots,0)^{\mathrm{T}}考察 n 阶矩阵 H=I-ye_l^{\mathrm{T}}，其中 y=(y_1,\cdots,y_n)\in\mathbb{R}^n，则
H= \begin{pmatrix} 1& & &{}-y_1& & &\\[4pt] &\ddots&&\vdots&&&\\ &&1&{}-y_{l-1}\\[2pt] &&&1-y_l\\[2pt] &&&{}-y_{l+1}&1\\[4pt] &&&\vdots&&\ddots\\ &&&{}-y_n&&&1 \end{pmatrix}具有以上形式的矩阵称为Gauss-Jordan 矩阵。矩阵 H 拥有一些优良的性质：

• 左乘的效果：HA=A 的每一行减去第 l 行的若干倍
  更简洁地说，就是用 A 的第 l 行修正 A 的所有行。设 A 的第 l 行为 A_l，则
  
  \begin{aligned} HA&=(I-ye_l^{\mathrm{T}})A\\ &=A-ye_l^{\mathrm{T}}A\\ &=A-yA_l\\ &=A-\begin{pmatrix} y_1A_l\\[4pt] \vdots\\ y_nA_l \end{pmatrix} \end{aligned}
  由上式可知，在 n 阶矩阵 A 上左乘 H 的效果是：第 i 行减去 y_i 倍的 A_l, (1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}n)。

• 右乘的效果：AH=A 的第 l 列减去所有行的线性组合
  更简洁地说，就是用 A 的所有列修正 A 的第 l 列。设 A=(P_1,\cdots,P_n)，则
  
  \begin{aligned} AH&=A(I-ye_l^{\mathrm{T}})\\ &=A-Aye_l^{\mathrm{T}}\\ &=A-(y_1P_1+\cdots+y_nP_n)e_l^\mathrm{T} \end{aligned}
  由上式可知，在 n 阶矩阵 A 上右乘 H 的效果是：第 l 列减去 y_i 倍的 P_i, (1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}n)。

• 逆矩阵：
  H=I-ye_l^\mathrm{T} 为可逆矩阵当且仅当 y_l\neq1。并且当 y_l\neq1 时，
  H^{-1}=I+\frac{ye_l^{\mathrm{T}}}{1-y_l} 证明：线性代数有这样的结果：当 A\in\mathbb{R}^{m\times n}, B\in\mathbb{R}^{n\times m} 时，\det(I_m-AB)=\det(I_n-BA)。
  因此 \det(H)=\det(I-ye_l^\mathrm{T})=1-e_l^\mathrm Ty=1-y_l。于是 H 可逆 \Leftrightarrow y_l\neq1。
  
  在 y_l\neq1 的情况下，可直接验证 \displaystyle(I-ye_l^{\mathrm{T}})\bigg(I+\frac{ye_l^{\mathrm{T}}}{1-y_l}\bigg)=I，因此命题得证。

---

有了以上的准备工作后，就可以开始进行基变换了。设原来的可行基向量组是
B=(P_{\pi(1)},\cdots,P_{\pi(l)},\cdots,P_{\pi(m)}) 如果将第 l 个基向量 P_{\pi(l)} 替换为非基向量 P_k~(k 是使得 \lambda_k<0 的那个下标)，则得到的新的向量组为
B'=(P_{\pi(1)},\cdots,P_{k},\cdots,P_{\pi(m)}) 由于 B' 仅在第 l 列上与 B 不同，故可考虑 m 阶 Gauss-Jordan 矩阵 H=I-ye_l^\mathrm{T} 使得 B'=BH。于是
\begin{aligned} &B'=B(I-ye_l^\mathrm T)\\ \Rightarrow~&B'=B-Bye_l^\mathrm T\\ \Rightarrow~& Bye_l^\mathrm T=B-B'=(0,\cdots,P_{\pi(l)}-P_k,\cdots,0)\\ \Rightarrow~& By=P_{\pi(l)}-P_k\\ \Rightarrow~& y=e_l-P_k' \end{aligned} 根据最优性检验的假设，P_k' 至少有一分量为正。注意到 y_l=1-\alpha_{lk}，因此当 \alpha_{lk}>0 时，一定有 H 可逆，故 B'=BH 可逆，从而 B' 是 A 的一个基向量组。

---

关于基向量组 B'，需要证明以下结论：

• B' 为可行基向量组。
  只需证明 B' 对应的基本解的基变量部分 (B')^{-1}b\geqslant0。易见
  \begin{aligned} &(B')^{-1}b\geqslant0\\ \Leftrightarrow~&H^{-1}B^{-1}b\geqslant0\\ \Leftrightarrow~&\bigg(I+\frac{ye_l^\mathrm{T}}{1-y_l}\bigg)\beta\geqslant0\\ \Leftrightarrow~&\alpha_{lk}\beta+y\beta_l\geqslant0\\ \Leftrightarrow~&\alpha_{lk}\beta_i+y_i\beta_l\geqslant0~~(1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}m) \end{aligned} i=l 时，上式等价于 \beta_l\geqslant0，不等式成立。i\neq l 时，等价于
  \alpha_{lk}\beta_i-\alpha_{ik}\beta_l\geqslant0 当 \alpha_{ik}\leqslant0 时上式直接成立，而 \alpha_{ik}>0 时这又等价于
  \frac{\beta_{l}}{\alpha_{lk}}\leqslant\frac{\beta_i}{\alpha_{ik}} 因此，l 的选择应当使其满足： \frac{\beta_l}{\alpha_{lk}}=\min\bigg\{ \frac{\beta_i}{\alpha_{ik}}\bigg|~\alpha_{ik}>0,~1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}m \bigg\}
  这就是著名的比值最小原则。

• B' 对应的基本解处的目标函数值不大于 x^{(0)} 处的函数值。
  因为目标函数的简化形式为 z=z_0+\lambda^\mathrm{T}x，因此在 B' 处 z=z_0+\lambda_k((B')^{-1}b)_l=z_0+\lambda_k\dfrac{\beta_l}{\alpha_{lk}}\leqslant z_0。
  于是进行一轮基变换之后，目标函数值不会增加。

---

最后给出基变换之后的参数变化情况：
记 y=e_l-P_k'，H=I-ye_l^\mathrm T，则：

• \alpha=B^{-1}A：
  \alpha'=(B')^{-1}A=H^{-1}\alpha=\alpha+\dfrac{y\alpha_l}{\alpha_{lk}} (\alpha_l 表示 \alpha 的第 l 行)
• \beta=B^{-1}b：
  \beta'=(B')^{-1}b=\beta+\dfrac{y\beta_l}{\alpha_{lk}}
• \lambda^\mathrm T=c^\mathrm T-c_B^\mathrm T\alpha
  由于 c_B^\mathrm T=(c_{\pi(1)},\cdots,c_{\pi(l)},\cdots,c_{\pi(m)})^\mathrm T，c_{B'}^\mathrm T=(c_{\pi(1)},\cdots,c_{k},\cdots,c_{\pi(m)})^\mathrm T，
  故 c_{B'}^\mathrm T-c_B^\mathrm T=\delta e_l^\mathrm T，其中 \delta=c_k-c_{\pi(l)}。于是
  \begin{aligned} (\lambda')^\mathrm T&=c^\mathrm T-c_{B'}^\mathrm T\alpha'\\ &=c^\mathrm T-(c_B^\mathrm T+\delta e_l^\mathrm T)\bigg(\alpha+\frac{y\alpha_l}{\alpha_{lk}}\bigg)\\ &=c^\mathrm T-c_B^\mathrm T\alpha-\delta e_l^\mathrm T\bigg(\alpha+\frac{y\alpha_l}{\alpha_{lk}}\bigg)-c_B^\mathrm Ty\cdot\frac{\alpha_l}{\alpha_{lk}}\\ &=\lambda^\mathrm T-\delta\bigg(\alpha_l+\frac{y_l\alpha_l}{\alpha_{lk}}\bigg)-c_B^\mathrm{T}y\cdot\frac{\alpha_l}{\alpha_{lk}}\\ &=\lambda^\mathrm T-(c_B^\mathrm Ty+\delta)\frac{\alpha_l}{\alpha_{lk}} \end{aligned} 而 \begin{aligned} c_B^\mathrm Ty+\delta&=c_B^\mathrm T(e_l-P_k')+c_k-c_{\pi(l)}\\ &=c_{\pi(l)}-c_B^\mathrm TP_k'+c_k-c_{\pi(l)}\\ &=c_k-c_B^\mathrm TP_k'=\lambda_k \end{aligned} 故
  (\lambda')^\mathrm T=\lambda^\mathrm T-\lambda_k\frac{\alpha_l}{\alpha_{lk}}
• z_0=c_B^\mathrm T\beta：
  z_0'=z_0+\lambda_k\dfrac{\beta_l}{\alpha_{lk}} (上文已经证明的结果)

---

总而言之，基变换从可行基 B 出发构造了一组新的可行基 B'，并且保证目标函数在基本解处的值不增。每一轮基变换之后，原先的参数 (\alpha,\beta,\lambda,z_0) 就会更新为 (\alpha',\beta',\lambda',z_0') 。

Bland 规则

单纯形法交替执行最优性检测和基变换，Bland 规则提供了一种简单的办法使算法不会陷入循环：

• 有多个 \lambda_k<0 时，取其中 k 最小的。
• 有多个 \dfrac{\beta_l}{\alpha_{kl}} 取最小值时，取其中 l 最小的。

---

综合以上的讨论，单纯形法的伪代码如下：

算法(单纯形法)： 在 Ax=b 和 x\geqslant0 的条件下，求 c^\mathrm{T}x 的最小值。

• 输入：A\in\mathbb{R}^{m\times n},~b\in\mathbb{R}^{m\times1},~c\in\mathbb{R}^{n\times1},~\pi\in\mathbb{R}^m
  其中 \pi 指定了 m 个基分量的下标。
• 输出：\pi\in\mathbb{R}^m,~\beta\in\mathbb{R}^m,~z_0\in\mathbb{R}
  其中 \beta 是基本可行解的 m 个基分量，z_0 是目标函数最小值。

---

1. 初始化参数
B=A(:,\pi) % 基向量组，m*m
\alpha=B^{-1}A % m*n
\beta=B^{-1}b % 约化的基本可行解，m*1
c_B=c(\pi) % m*1
\lambda=c-\alpha^\mathrm Tc_B % 检验向量，n*1
z_0=c_B^\mathrm T\beta % 目标函数最小值

---

2. 最优性检验/基变换
找到最小的 k 使得 \lambda(k)<0

k 不存在：返回 \pi,~\beta,~z_0

k 存在：检查 \alpha(:,k) 中是否有正数

不存在：返回错误信息

存在：选择 l 使得 \displaystyle\frac{\beta(l)}{\alpha(l,k)}=\min\bigg\{\frac{\beta(i)}{\alpha(i,k)}\bigg|~\alpha(i,k)>0,~1\hspace{-1pt}\leqslant \hspace{-1pt}i\hspace{-1pt}\leqslant \hspace{-1pt}m \bigg\}
\pi(l)=k % 更换基分量
y=e_l-\alpha(:,k) % 辅助计算

z_0=z_0+\lambda(k)\cdot\dfrac{\beta(l)}{\alpha(l,k)} % 更新z0

\lambda=\lambda-\lambda(k)\dfrac{\alpha(l,:)^\mathrm T}{\alpha(l,k)} % 更新lambda

\beta=\beta+\dfrac{y\beta(l)}{\alpha(l,k)} % 更新beta

\alpha=\alpha+\dfrac{y\alpha(l,:)}{\alpha(l,k)} % 更新alpha

goto 2