第 8 章 凸优化

CHAPTER STATUS编写进度1%查看全书进度

8.1 凸集

8.1.1 基本凸集

8.1.1.1 直线

Definition 8.1. 设\(x_1,x_2\in\mathbb{R}^{n}\)且\(x_1\ne x_2\)。若\(\theta\in[0,1]\),称所有满足形式:

\[\begin{equation*} y=\theta x_1+(1-\theta)x_2 \end{equation*}\]

的点\(y\)构成的集合为线段(line segment)。若\(\theta\in\mathbb{R}\),称所有满足上述形式的点\(y\)构成的集合为直线(line)。

note 8.1. 上述形式等价于:

\[\begin{equation*} y=x_2+\theta(x_1-x_2) \end{equation*}\]

这给出了另一种直观:\(y\)是基\(x_2\)和被\(\theta\)缩放的方向\(x_1-x_2\)的和。

8.1.1.2 凸集

Definition 8.2. 设\(E\subseteq\mathbb{R}^{n}\)。若对于任意不同的\(x_1,x_2\in E\),由\(x_1,x_2\)确定的线段都在\(E\)中,则称\(E\)是凸集(convex set)。

Definition 8.3. 设\(n\in\mathbb{N}^+\),\(\theta_1, \theta_2, \dots, \theta_{n}\geqslant0\)且\(\sum\limits_{i=1}^{n}\theta_i=1\),\(x_1, x_2, \dots, x_{n}\in\mathbb{R}^{n}\)。称:

\[\begin{equation*} \theta_1x_1+\theta_2x_2+\cdots+\theta_nx_n \end{equation*}\]

为\(x_1, x_2, \dots, x_{n}\)构成的一个凸组合(convex combination)。

Definition 8.4. 设\(E\subseteq\mathbb{R}^{n}\)。定义:

\[\begin{equation*} \left\{\sum_{i=1}^{n}\theta_ix_i:x_i\in E,\;\theta_i\geqslant0,\;\sum_{i=1}^{n}\theta_i=1,\;n\in\mathbb{N}^+\right\} \end{equation*}\]

为\(E\)的凸包(convex hull),记作\(\operatorname{Conv}(E)\)。

TikZ 数学图示
图 8.1 凸集与非凸集的二维几何结构。

Property 8.1.1. 设\(E\subseteq\mathbb{R}^{n}\)。凸集具有如下性质:

  1. \(E\)是凸集当且仅当\(E\)中任意有限个元素的凸组合都在\(E\)中;

  2. \(E\)的凸包是包含\(E\)的最小凸集;

  3. 凸集的交集是凸集;

  4. 设\(A\in M_{m\times n}(\mathbb{R}),\;b\in\mathbb{R}^{n}\),\(E\)是凸集,则\(AE+b=\{Ax+b:x\in E\}\)也是凸集;

  5. 设\(E_1,E_2\)是凸集,则\(E_1\pm E_1=\{x\pm y:x\in E_1,\;y\in E_2\}\)也是凸集;

  6. 设\(E_1,E_2\)是凸集,则\(\{(x,y_1+y_2):(x,y_1)\in E_1,\;(x,y_2)\in E_2\}\)也是凸集;

  7. 若\(E\)是凸集,则\(\{x_1\in\mathbb{R}^{m}:(x_1,x_2)\in E,\;m\in\mathbb{N}^+,\;m\leqslant n\}\)也是凸集;

  8. 若\(E\)是凸集,则\(\overline{E}\)也是凸集;

证明. (1)必要性:当\(n=2\)时,由凸集的定义可知\(E\)中任意\(2\)个元素的凸组合都在\(E\)中。假设结论对\(n-1\)成立,下面证明对\(n\)个也成立。

任取\(\theta_1, \theta_2, \dots, \theta_{n}\geqslant0\)且\(\sum\limits_{i=1}^{n}\theta_i=1\),\(x_1, x_2, \dots, x_{n}\in E\),则:

\[\begin{equation*} \sum_{i=1}^{n}\theta_ix_i=\sum_{i=1}^{n-1}\theta_ix_i+\theta_nx_n=\sum_{i=1}^{n-1}\theta_ix_i+\left(1-\sum_{i=1}^{n-1}\theta_i\right)x_i \end{equation*}\]

若\(\sum\limits_{i=1}^{n-1}\theta_i=0\),则上式退化为\(x_n\),结论成立。若\(\sum\limits_{i=1}^{n-1}\theta_i\ne0\),则:

\[\begin{equation*} y=\sum_{i=1}^{n-1}\frac{\theta_i}{\sum\limits_{j=1}^{n-1}\theta_j}x_i \end{equation*}\]

是关于\(x_1, x_2, \dots, x_{n-1}\)的一个凸组合,由归纳假设可知\(y\in E\)。变形可得:

\[\begin{equation*} \sum_{i=1}^{n}\theta_ix_i=\sum_{i=1}^{n-1}\theta_ix_i+\left(1-\sum_{i=1}^{n-1}\theta_i\right)x_i=\sum\limits_{j=1}^{n-1}\theta_jy+\left(1-\sum_{i=1}^{n-1}\theta_i\right)x_i\in E \end{equation*}\]

必要性得证。

充分性:任意有限个包含\(2\)个。

(2)由(1)立即可得。

(3)任取指标集\(I\)和凸集\(C_i,i\in I\),令\(C=\underset{i\in I}{\overset{}{\cap}}C_i\)。任取\(x_1,x_2\in C\)和\(\theta\geqslant0\),则对任意的\(i\in I\)有\(x_1,x_2\in C_i,\;\theta x_1+(1-\theta)x_2\in C_i\),于是\(\theta x_1+(1-\theta)x_2\in C\),\(C\)是凸集。

(4)(5)(6)(7)利用凸集的定义可直接得到证明,略去。

(8)任取\(x,y\in\overline{E}\),于是存在\(E\)中的\(\{x_n\}\)和\(\{y_n\}\)满足\(\{x_n\}\to x\)和\(\{y_n\}\to y\)。任取\(\theta\in[0,1]\),因为\(E\)是凸集,所以\(\theta x_n+(1-\theta)y_n\in E\)。由性质 3.4.1(3)可知\(\theta x+(1-\theta)y\in\overline{E}\)。根据\(x,y\)的任意性可得\(\overline{E}\)是凸集。 ◻

8.1.1.3 仿射集

Definition 8.5. 设\(E\subseteq\mathbb{R}^{n}\)。若对于任意不同的\(x_1,x_2\in E\),由\(x_1,x_2\)确定的直线都在\(E\)中,则称\(E\)是仿射集(affine set)。

Definition 8.6. 设\(n\in\mathbb{N}^+\),\(\theta_1, \theta_2, \dots, \theta_{n}\in\mathbb{R}\)且\(\sum\limits_{i=1}^{n}\theta_i=1\),\(x_1, x_2, \dots, x_{n}\in\mathbb{R}^{n}\)。称:

\[\begin{equation*} \theta_1x_1+\theta_2x_2+\cdots+\theta_nx_n \end{equation*}\]

为\(x_1, x_2, \dots, x_{n}\)构成的一个仿射组合(affine combination)。

Definition 8.7. 设\(E\subseteq\mathbb{R}^{n}\)。定义:

\[\begin{equation*} \left\{\sum_{i=1}^{n}\theta_ix_i:x_i\in E,\;\theta_i\in\mathbb{R},\;\sum_{i=1}^{n}\theta_i=1,\;n\in\mathbb{N}^+\right\} \end{equation*}\]

为\(E\)的仿射包(affine hull),记作\(\operatorname{Aff}(E)\)。

Property 8.1.2. 设\(E\subseteq\mathbb{R}^{n}\)。仿射集具有如下性质:

  1. \(E\)是仿射集当且仅当\(E\)中任意有限个元素的仿射组合都在\(E\)中;

  2. \(E\)是仿射集当且仅当存在一个\(\mathbb{R}^{n}\)的子空间\(V\)使得\(E=\{x+x_0:x\in V,\;x_0\in E\}\),并且表达式与\(x_0\)的选择无关;

  3. \(E\)的仿射包是包含\(E\)的最小仿射集;

  4. 仿射集是凸集;

证明. (1)完全类似于性质 8.1.1(1)。

(2)必要性:设\(E\)是仿射集,任取\(x_0\in E\),令\(V=\{x-x_0:x\in E\}\),则\(E=\{x+x_0:x\in V,\;x_0\in E\}\)。

任取\(\alpha,\beta\in\mathbb{R}\)和\(x_1,x_2\in V\),则:

\[\begin{equation*} \alpha x_1+\beta x_2+x_0=\alpha(x_1+x_0)+\beta(x_2+x_0)+(1-\alpha-\beta)x_0\in E \end{equation*}\]

于是\(\alpha x_1+\beta x_2\in V\)。由定理 1.5可知\(V\)是子空间。

取不同于\(x_0\)的\(x_1\in E\)构成集合\(V'=\{x-x_1:x\in C\}\)。任取\(V'\)中的元素\(\alpha=x-x_1\),则\(\alpha=x+x_1-x_0\)是\(x,x_1,x_0\)构成的一个仿射组合,于是\(\alpha\in V\)。由\(\alpha\)的任意性可知\(V'\subseteq V\)。同理可得\(V\subseteq V'\),于是有\(V=V'\),即表达式与\(x_0\)的选择无关。

充分性:任取\(E\)中不同的两点\(x_1+x_0,x_2+x_0\)和\(\theta\in\mathbb{R}\),有:

\[\begin{equation*} \theta(x_1+x_0)+(1-\theta)(x_2+x_0)=\theta x_1+(1-\theta)x_2+x_0 \end{equation*}\]

因为\(V\)是子空间,由定理 1.5可知\(\theta x_1+(1-\theta)x_2\in V\),所以\(\theta x_1+(1-\theta)x_2+x_0\in E\)。由\(x_1+x_0,x_2+x_0\)和\(\theta\)的任意性可知\(E\)是仿射集。

(3)由(1)立即可得。

(4)由定义立即可得。 ◻

Definition 8.8. 设\(E\subseteq\mathbb{R}^{n}\)是仿射集。由性质 8.1.2(2)可知存在一个\(\mathbb{R}^{n}\)的子空间\(V\)使得\(E=\{x+x_0:x\in V,\;x_0\in E\}\),称\(V\)是\(E\)的方向子空间,记\(\dim E=\dim V\)。

Definition 8.9. 设\(x_0,x_1,\dots,x_m\in\mathbb{R}^n\)。若:

\[\begin{equation*} \sum_{i=0}^{m}\lambda_ix_i=0,\quad\sum_{i=0}^{m}\lambda_i=0 \end{equation*}\]

只能推出:

\[\begin{equation*} \lambda_0=\lambda_1=\cdots=\lambda_m=0 \end{equation*}\]

则称\(x_0,x_1,\dots,x_m\)仿射无关(affine independent)。

Property 8.1.3. 设\(x_0,x_1,\dots,x_m\in\mathbb{R}^n\)。仿射无关具有以下性质:

  1. \(x_0,x_1,\dots,x_m\)仿射无关当且仅当对任意的或某个\(j\in\{0,1,\dots,m\}\)有\(\{x_i-x_j:i\neq j\}\)线性无关;

  2. 对任意的\(a\in\mathbb{R}^n\),\(x_0,x_1,\dots,x_m\)仿射无关当且仅当\(x_0+a,x_1+a,\dots,x_m+a\)仿射无关;

  3. 若\(x_0,x_1,\dots,x_m\)仿射无关,则对于任意的\(j\in\{0,1,\dots,m\}\)有\(\operatorname{Aff}(\{x_0,x_1,\dots,x_m\})=x_j+\operatorname{span}\{x_i-x_j:i\ne j\}\)且\(\dim\operatorname{Aff}(\{x_0,x_1,\dots,x_m\})=m\);

  4. 若\(E\)是维数为\(m\)的仿射集,则存在\(m+1\)个仿射无关的点\(x_0,x_1,\dots,x_m\in E\)使得\(E=\operatorname{Aff}(\{x_0,x_1,\dots,x_m\})\);

证明. (1)必要性:任取\(j\in\{0,1,\dots,m\}\),设:

\[\begin{equation*} \sum_{i\ne j}\lambda_i(x_i-x_j)=0 \end{equation*}\]

于是:

\[\begin{equation*} \sum_{i\ne j}\lambda_i(x_i-x_j)=\sum_{i\ne j}\lambda_ix_i-\sum_{i\ne j}\lambda_ix_j \end{equation*}\]

记\(\lambda_j=-\sum\limits_{i\ne j}\lambda_i\),则;

\[\begin{equation*} \sum_{i\ne j}\lambda_i(x_i-x_j)=\sum_{i=0}^{m}\lambda_ix_i,\quad\sum_{i=0}^{m}\lambda_i=0 \end{equation*}\]

因为\(x_0,x_1,\dots,x_m\)仿射无关,所以对任意的\(i\ne j\)都有\(\lambda_i=0\),即\(\{x_i-x_j:i\neq j\}\)线性无关。由\(j\)的任意性可知必要性成立。

充分性:任取\(j\in\{0,1,\dots,m\}\),设:

\[\begin{equation*} \sum_{i=0}^{m}\lambda_ix_i,\quad\sum_{i=0}^{m}\lambda_i=0 \end{equation*}\]

则\(\lambda_j=-\sum\limits_{i\ne j}\lambda_i\),代入上式可得:

\[\begin{equation*} \sum_{i=0}^{m}\lambda_ix_i=\sum_{i\ne j}\lambda_ix_i-\sum_{i\ne j}\lambda_ix_j=\sum_{i\ne j}^{}\lambda_i(x_i-x_j) \end{equation*}\]

因为\(\{x_i-x_j:i\ne j\}\)线性无关,所以对任意的\(i\ne j\)都有\(\lambda_i=0\),于是\(\lambda_j=0\),即\(x_0,x_1,\dots,x_m\)仿射无关。

(2)必要性:若:

\[\begin{equation*} \sum_{i=0}^{m}\lambda_i(x_i+a)=0,\quad\sum_{i=0}^{m}\lambda_i=0 \end{equation*}\]

注意到:

\[\begin{equation*} \sum_{i=0}^{m}\lambda_i(x_i+a)=\sum_{i=0}^{m}\lambda_ix_i+\sum_{i=0}^{m}\lambda_ia=\sum_{i=0}^{m}\lambda_ix_i \end{equation*}\]

因为\(x_0,x_1,\dots,x_m\)仿射无关,所以\(\lambda_0=\lambda_1=\cdots=\lambda_m=0\),于是\(x_0+a,x_1+a,\dots,x_m+a\)仿射无关。

充分性:由定义立即可得。

(3)任取\(x\in\operatorname{Aff}(\{x_0,x_1,\dots,x_m\})\)和\(j\in\{0,1,\dots,m\}\),则存在\(\lambda_i\in\mathbb{R},\;i\ne j\)使得:

\[\begin{equation*} x=\sum_{i\ne j}\lambda_ix_i+\left(1-\sum_{i\ne j}\lambda_i\right)x_j=x_j+\sum_{i\ne j}\lambda_i(x_i-x_j) \end{equation*}\]

由\(x\)和\(j\)的任意性可知对于任意的\(j\in\{0,1,\dots,m\}\)有\(\operatorname{Aff}(\{x_0,x_1,\dots,x_m\})=x_j+\operatorname{span}\{x_i-x_j:i\ne j\}\)。

根据性质 1.1.5(1)、性质 8.1.2(2)、(1)和性质 1.1.5(2)可知\(\dim\operatorname{Aff}(\{x_0,x_1,\dots,x_m\})=m\)。

(4)任取\(x_0\in E\),由性质 8.1.2(2)可知存在一个\(\mathbb{R}^{n}\)的子空间\(V\)使得\(E=\{x+x_0:x\in V\}\)。取\(V\)的一组基\(\alpha_1, \alpha_2, \dots, \alpha_{m}\),记\(x_i=x_0+v_i\in E,\;i=1,2,\dots,m\)。根据(1)可知\(x_0,x_1,\dots,x_m\)仿射无关。由(3)和性质 8.1.2(2)可知结论成立。 ◻

8.1.1.4 相对内部

Definition 8.10. 设\(E\)是\(\mathbb{R}^{n}\)的非空子集,\(A=\operatorname{Aff}(E)\)。定义:

\[\begin{equation*} \operatorname{ri}(E)\coloneq\{x\in E:\exists\;r>0,\;U(x,r)\cap A\nsubseteq E\} \end{equation*}\]

为\(E\)的相对内部(relative interior),其中的点称为\(E\)的相对内点。约定\(\operatorname{ri}(\varnothing)=\varnothing\)。

Property 8.1.4. 设\(E\subseteq\mathbb{R}^{n}\)为非空凸集,则:

  1. \(\operatorname{ri}(E)\ne\varnothing\);

  2. 若\(x\in\operatorname{ri}(E),\;y\in\overline{E}\),则\((1-\theta)x+\theta y\in\operatorname{ri}(E)\),\(\forall\;\theta\in[0,1)\);特别地,\(\operatorname{ri}(E)\)是凸集;

  3. \(\operatorname{ri}(\overline{E})=\operatorname{ri}(E)\);

证明. (1)设\(A=\operatorname{Aff}(E)\),根据性质 8.1.2(2)可知存在一个\(\mathbb{R}^{n}\)的子空间\(V\)使得\(A=\{x+x_0:x\in V,x_0\in E\}\),记\(\dim E=\dim V=m\)。

若\(m=0\),则\(E\)为单点集,结论成立。

若\(m>0\), ◻

8.1.2 分离超平面定理

Definition 8.11. 设\(a\in\mathbb{R}^{n},\;b\in\mathbb{R}\)。称\(\{x\in\mathbb{R}^{n}:a^{\top}x=b\}\)为超平面(hyperplane),\(\{x\in\mathbb{R}^{n}:a^{\top}x\leqslant b\}\)为闭半空间(halfspace),\(\{x\in\mathbb{R}^{n}:a^{\top}x< b\}\)为开半空间。

引理 9.1 Lemma 8.1. 设\(X\)是实Hilbert空间,\(E\)是\(X\)中的凸闭集。对任意的\(x\in X\),\(y\)是\(x\)在\(E\)中的最佳逼近元当且仅当:

\[\begin{equation*} \forall\;z\in E,\;(x-y,z-y)\leqslant0 \end{equation*}\]

证明. (1)必要性:对任意的\(z\in E\)和\(\theta\in[0,1]\),因为\(E\)是\(X\)中的凸集,所以\(\theta z+(1-\theta)y\in E\)。由于\(y\)是\(x\)在\(E\)中的最佳逼近元,结合:

\[\begin{equation*} ||x-y||^2\leqslant||x-[\theta z+(1-\theta)y]||^2=||x-y-\theta(z-y)||^2=||x-y||^2-2\theta(x-y,z-y)+\theta^2||z-y||^2 \end{equation*}\]

可知:

\[\begin{equation*} -2\theta(x-y,z-y)+\theta^2||z-y||^2\geqslant0 \end{equation*}\]

取\(\theta\to0^+\),由性质 3.3.3(4)可知\((x-y,z-y)\leqslant0\)。

(2)充分性:此时对任意的\(z\in E\)有:

\[\begin{align*} &||x-z||=||x-y+y-z||^2=||x-y||^2+2(x-y,y-z)+||y-z||^2 \\ =&|x-y||^2-2(x-y,z-y)+||y-z||^2\geqslant||x-y|| \end{align*}\]

于是\(y\)是\(x\)在\(E\)中的最佳逼近元。 ◻

Theorem 8.1. (Hyperplane Separation Theorem)

设\(A\)和\(B\)是\(\mathbb{R}^{n}\)上不交的非空凸集,则存在非零的\(a\in\mathbb{R}^{n}\)和\(b\in\mathbb{R}\)使得:

\[\begin{equation*} \forall\;x\in A,\;(a,x)\geqslant b,\quad\forall\;y\in B,\;(a,y)\leqslant b \end{equation*}\]

记\(C=A-B=\{x-y:x\in A,\;y\in B\}\),若\(\mathbf{0}\notin\overline{C}\),则存在\(a\in\mathbb{R}^{n}\)和\(b\in\mathbb{R}\)使得上式取严格的不等号。

证明. 由性质 8.1.1(5)(8)可知\(C\)是凸集且\(\overline{C}\)是闭凸集。

若\(\mathbf{0}\notin\overline{C}\),由定理 3.27可知在\(\overline{C}\)中存在\(\mathbf{0}\)的唯一的最佳逼近元\(a\)。根据引理 8.1可得:

\[\begin{equation*} \forall\;z\in\overline{C},\;(-a,z-a)\leqslant0 \end{equation*}\]

也即\((a,z)\geqslant||a||^2\)。对任意的\(x\in A\)和\(y\in B\)有\(x-y\in\overline{C}\),所以:

\[\begin{equation*} (a,x)\geqslant(a,y)+||a||^2 \end{equation*}\]

根据上式与确界原理,可记:

\[\begin{equation*} \alpha=\inf_{x\in A}(a,x),\quad\beta=\sup_{y\in B}(a,y) \end{equation*}\]

由上下确界的不等式性可知:

\[\begin{equation*} \alpha\geqslant\beta+||a||^2>\beta \end{equation*}\]

任取\(b\in(\beta,\alpha)\)可得:

\[\begin{equation*} \forall\;x\in A,\;(a,x)\geqslant\alpha>b,\quad\forall\;y\in B,\;(a,y)\leqslant\beta<b \end{equation*}\]

若\(\mathbf{0}\in\overline{C}\),因为\(\mathbf{0}\not\in C\)且\(\operatorname{ri}(C)\subseteq C\),所以\(\mathbf{0}\notin\operatorname{ri}(C)\)。由性质 8.1.4(3)可知\(\mathbf{0}\notin\operatorname{ri}(\overline{C})\),即:

\[\begin{equation*} \forall\;r>0,\;U(\mathbf{0},r)\cap\operatorname{Aff}(\overline{C})\nsubseteq\overline{C} \end{equation*}\]

于是对任意的\(n\in\mathbb{N}^+\),存在\(x_n\in\operatorname{Aff}(\overline{C})\backslash\overline{C}\)满足\(||x_n||\leqslant\dfrac{1}{n}\)。根据定理 3.27可知在\(\overline{C}\)中存在\(x_n\)的唯一的最佳逼近元\(y_n\)。由引理 8.1可得:

\[\begin{equation*} \forall\;z\in\overline{C},\;(x_n-y_n,z-y_n)\leqslant0 \end{equation*}\]

即:

\[\begin{equation*} \forall\;z\in\overline{C},\;\left(\frac{y_n-x_n}{||y_n-x_n||},z\right)\geqslant\left(\frac{y_n-x_n}{||y_n-x_n||},y_n\right) \end{equation*}\]

因为单位球面是有界闭集(因为\(\{1\}\)是闭集,根据性质 3.4.1(2)和定理 3.18可知\(\{x;||x||_2=1\}\)是闭集。),由定理 3.9可知\(\left\{\dfrac{y_n-x_n}{||y_n-x_n||}\right\}\to a\)且\(||a||=1\)。根据最佳逼近元的定义可知:

\[\begin{equation*} ||y_n||=||y_n-x_n+x_n-\mathbf{0}||\leqslant||x_n-y_n||+||x_n||\leqslant||x_n-\mathbf{0}||+||x_n||\to0 \end{equation*}\]

所以\(\{y_n\}\to\mathbf{0}\)。由性质 3.4.5(2)(1)可知:

\[\begin{equation*} \forall\;z\in\overline{C},\;(a,z)\geqslant0 \end{equation*}\]

对任意的\(x\in A\)和\(y\in B\)有\(x-y\in\overline{C}\),所以:

\[\begin{equation*} (a,x)\geqslant(a,y) \end{equation*}\]

根据上式与确界原理,可记:

\[\begin{equation*} \alpha=\inf_{x\in A}(a,x),\quad\beta=\sup_{y\in B}(a,y) \end{equation*}\]

由上下确界的不等式性可知\(\alpha\geqslant\beta\)。任取\(b\in[\beta,\alpha]\)可得:

\[\begin{equation*} \forall\;x\in A,\;(a,x)\geqslant b,\quad\forall\;y\in B,\;(a,y)\leqslant b \end{equation*}\]

◻

Definition 8.12. 设\(E\subseteq\mathbb{R}^{n}\),\(x_0\in\partial E\)。若非零的\(\alpha\in\mathbb{R}^{n}\)和\(b\in\mathbb{R}\)满足对任意的\(x\in E\)有\((\alpha,x)\leqslant(\alpha,x_0)\),则称超平面\(\{x\in\mathbb{R}^{n}:(\alpha,x)=(\alpha,x_0)\}\)为\(E\)在\(x_0\)处的支撑超平面(supporting hyperplane)。

Theorem 8.2. (Supporting Hyperplane Theorem)

若\(E\subseteq\mathbb{R}^{n}\)是凸集,则在\(E\)的任意边界点处都存在支撑超平面。

8.1.2.1 锥

Definition 8.13. 设\(E\subseteq\mathbb{R}^{n}\)。若对于任意的\(x\in E\)和\(\theta\geqslant0\),有\(\theta x\in E\),则称\(E\)是锥(cone)。若锥\(E\)还满足对任意的\(x_1,x_2\in E\)和\(\theta_1,\theta_2\geqslant0\)有\(\theta_1x_1+\theta_2x_2\in E\),则称\(E\)为凸锥(convex cone)。

Definition 8.14. 设\(n\in\mathbb{N}^+\),\(\theta_1, \theta_2, \dots, \theta_{n}\geqslant0\),\(x_1, x_2, \dots, x_{n}\in\mathbb{R}^{n}\)。称:

\[\begin{equation*} \theta_1x_1+\theta_2x_2+\cdots+\theta_nx_n \end{equation*}\]

为\(x_1, x_2, \dots, x_{n}\)构成的一个锥组合(Conic combination)。

Definition 8.15. 设\(E\subseteq\mathbb{R}^{n}\)。定义:

\[\begin{equation*} \left\{\sum_{i=1}^{n}\theta_ix_i:x_i\in E,\;\theta_i\geqslant0,\;n\in\mathbb{N}^+\right\} \end{equation*}\]

为\(E\)的锥包(conic hull),记作\(\operatorname{Cone}(E)\)。

TikZ 数学图示
图 8.2 锥只要求沿每条射线对非负缩放封闭;凸锥还要求对加法封闭,因此包含两条射线之间的全部非负线性组合。

Property 8.1.5. 设\(E\subseteq\mathbb{R}^{n}\)。锥具有如下性质:

  1. \(E\)是锥包当且仅当\(E\)中任意有限个元素的锥组合都在\(E\)中;

  2. \(E\)的锥包是包含\(E\)的最小锥;

证明. (1)完全类似于性质 8.1.1(1)。

(2)由(1)立即可得。 ◻

8.1.2.2 超平面与半空间

Property 8.1.6. 超平面与半空间具有如下性质:

  1. 超平面是仿射集;

  2. 半空间是凸集;

证明. (1)任取超平面\(\{x\in\mathbb{R}^{n}:a^{\top}x=b\}\)和其中任意两点\(x_1,x_2\),再任取\(\theta\in\mathbb{R}\),有:

\[\begin{equation*} a^{\top}[\theta x_1+(1-\theta)x_2]=\theta b+(1-\theta)b=b \end{equation*}\]

所以\(\theta x_1+(1-\theta)x_2\in\{x\in\mathbb{R}^{n}:a^{\top}x=b\}\),超平面是仿射集。

(2)任取闭半空间\(\{x\in\mathbb{R}^{n}:a^{\top}x\leqslant b\}\)和其中任意两点\(x_1,x_2\),再任取\(\theta\geqslant0\),有:

\[\begin{equation*} a^{\top}[\theta x_1+(1-\theta)x_2]\leqslant\theta b+(1-\theta)b=b \end{equation*}\]

所以\(\theta x_1+(1-\theta)x_2\in\{x\in\mathbb{R}^{n}:a^{\top}x\leqslant b\}\),闭半空间是凸集。开半空间同理可得。 ◻

8.1.2.3 多面体

Definition 8.16. 设\(A\in M_{m\times n}(\mathbb{R}),\;b\in\mathbb{R}^{n},\;C\in M_{p\times q}(\mathbb{R}),\;d\in\mathbb{R}^{q}\)。称\(\{x\in\mathbb{R}^{n}:Ax\leqslant b,\;Cx=d\}\)为多面体(polyhedron)。

Property 8.1.7. 多面体是凸集。

证明. 由性质 8.1.1(3)、性质 8.1.6(1)、性质 8.1.2(4)和性质 8.1.6(2)立即可得。 ◻

8.1.2.4 球与椭球

Definition 8.17. 设\(||\cdot||\)是\(\mathbb{R}^{n}\)上的一个范数,\(x_c\in\mathbb{R}^{n},\;r\in\mathbb{R}\)。称\(\{x\in\mathbb{R}^{n}:||x-x_c||\leqslant r\}\)为球(ball)。

Definition 8.18. 设\(P\in M_{n}(\mathbb{R})\)是正定矩阵,\(x_c\in\mathbb{R}^{n}\)。由定理 2.17(6)可得\(P\)的可逆性,称\(\{x\in\mathbb{R}^{n}:(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1\}\)为椭球(ellipsoid)。

Property 8.1.8. 球与椭球具有如下性质:

  1. 球是凸集;

  2. 椭球\(\{x\in\mathbb{R}^{n}:(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1\}\)等价于\(\{x_c+P^{\frac{1}{2}}u:u\in\mathbb{R}^{n},\;||u||_2\leqslant1\}\),由\(P\)的正定性可得\(P^{\frac{1}{2}}\)的存在性;

  3. 椭球是凸集;

证明. (1)任取球\(\{x\in\mathbb{R}^{n}:||x-x_c||\leqslant r\}\)和其中任意两点\(x_1,x_2\),再任取\(\theta\geqslant0\),有:

\[\begin{equation*} ||\theta x_1+(1-\theta)x_2-x_c||=||\theta(x_1-x_c)+(1-\theta)(x_2-x_c)||\leqslant\theta||x_1-x_c||+(1-\theta)||x_2-x_c||\leqslant\theta r+(1-\theta)r=r \end{equation*}\]

所以\(\theta x_1+(1-\theta)x_2\in\{x\in\mathbb{R}^{n}:||x-x_c||\leqslant r\}\),球是凸集。

(2)任取\(x\in\{x_c+P^{\frac{1}{2}}u:u\in\mathbb{R}^{n},\;||u||_2\leqslant1\}\),由性质 2.1.2(4)可得:

\[\begin{equation*} (x-x_c)^{\top}P^{-1}(x-x_c)=(P^{\frac{1}{2}}u)^{\top}P^{-1}P^{\frac{1}{2}}u=||u||_2\leqslant1 \end{equation*}\]

所以\(\{x_c+P^{\frac{1}{2}}u:u\in\mathbb{R}^{n},\;||u||_2\leqslant1\}\subseteq\{x\in\mathbb{R}^{n}:(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1\}\)。

任取\(x\in\{x\in\mathbb{R}^{n}:(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1\}\)。由\(P\)的正定性可得\(P^{-\frac{1}{2}}\)的存在性,取\(u=P^{-\frac{1}{2}}(x-x_c)\),根据性质 2.1.2(4)可得:

\[\begin{equation*} ||u||_2=[P^{-\frac{1}{2}}(x-x_c)]^{\top}P^{-\frac{1}{2}}(x-x_c)=(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1 \end{equation*}\]

所以\(\{x\in\mathbb{R}^{n}:(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1\}\subseteq\{x_c+P^{\frac{1}{2}}u:u\in\mathbb{R}^{n},\;||u||_2\leqslant1\}\)。

综上,\(\{x\in\mathbb{R}^{n}:(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1\}=\{x_c+P^{\frac{1}{2}}u:u\in\mathbb{R}^{n},\;||u||_2\leqslant1\}\)。

(3)任取椭球\(\{x\in\mathbb{R}^{n}:(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1\}\)和其中任意两点\(x_1=x_c+P^{\frac{1}{2}}u_1,x_2=x_c+P^{\frac{1}{2}}u_2\),再任取\(\theta\geqslant0\),有:

\[\begin{equation*} \theta x_1+(1-\theta)x_2=x_c+P^{\frac{1}{2}}[\theta u_1+(1-\theta)u_2] \end{equation*}\]

由(1)可知\(||\theta u_1+(1-\theta)u_2||_2\leqslant1\),所以\(\theta x_1+(1-\theta)x_2\in\{x\in\mathbb{R}^{n}:(x-x_c)^{\top}P^{-1}(x-x_c)\leqslant1\}\),椭球是凸集。 ◻

Definition 8.19. 设\(f:\mathbb{R}^{n}\to\mathbb{R}^{m}\)是函数。若\(f\)具有形式\(f(x)=Ax+b\),其中\(A\in M_{m\times n}(\mathbb{R}),\;b\in\mathbb{R}^{m}\),则称\(f\)是仿射函数(affine function)。

Property 8.1.9. 设\(f:\mathbb{R}^{n}\to\mathbb{R}^{m}\)是仿射函数,\(f=Ax+b\),其中\(A\in M_{m\times n}(\mathbb{R}),\;b\in\mathbb{R}^{m}\)。仿射函数具有如下性质:

  1. 若\(E\subseteq\mathbb{R}^{n}\)是凸集,则\(f(E)\)是\(\mathbb{R}^{m}\)中的凸集;

  2. 若\(E\subseteq\mathbb{R}^{m}\)是凸集,则\(f^{-1}(E)=\{x\in\mathbb{R}^{n}:f(x)\in E\}\)是\(\mathbb{R}^{n}\)中的凸集。

证明. (1)任取\(y_1,y_2\in f(E)\)和\(\theta\in[0,1]\),则存在\(x_1,x_2\in E\)使得\(f(x_1)=y_1,f(x_2)=y_2\),于是有:

\[\begin{align*} \theta y_1+(1-\theta)y_2&=\theta f(x_1)+(1-\theta)f(x_2)=\theta(Ax_1+b)+(1-\theta)(Ax_2+b) \\ &=A[\theta x_1+(1-\theta)x_2]+b=f[\theta x_1+(1-\theta)x_2] \end{align*}\]

因为\(E\)是凸集,所以\(\theta x_1+(1-\theta)x_2\in E\),于是\(\theta y_1+(1-\theta)y_2\in f(E)\)。由\(y_1,y_2,\theta\)的任意性可知\(f(E)\)是\(\mathbb{R}^{m}\)中的凸集。

(2)任取\(x_1,x_2\in f^{-1}(E)\)和\(\theta\in[0,1]\),则存在\(y_1,y_2\in E\)使得\(f(x_1)=y_1,f(x_2)=y_2\),于是有:

\[\begin{equation*} f[\theta x_1+(1-\theta)x_2]=A[\theta x_1+(1-\theta)x_2]+b=\theta(Ax_1+b)+(1-\theta)(Ax_2+b)=\theta y_1+(1-\theta)y_2 \end{equation*}\]

因为\(E\)是凸集,所以\(\theta y_1+(1-\theta)y_2\in E\),于是\(\theta x_1+(1-\theta)x_2\in f^{-1}(E)\)。由\(x_1,x_2,\theta\)的任意性可知\(f^{-1}(E)\)是\(\mathbb{R}^{n}\)中的凸集。 ◻

8.1.3 凸函数

Definition 8.20. 设\(X\)是一个非空集合,\(f:X\to\overline{R}^n\)。若存在\(x\in X\)使得\(f(x)<+\infty\)且对任意的\(x\in X\)有\(f(x)>-\infty\),则称\(f\)关于\(X\)是适当函数(proper function)。

Definition 8.21. 设\(f\)是关于凸集\(E\subseteq\mathbb{R}^{n}\)的适当函数。若对任意的\(x,y\in E\)和任意的\(\alpha\in[0,1]\)都有:

\[\begin{equation*} f[\alpha x+(1-\alpha)y]\leqslant\alpha f(x)+(1-\alpha)f(y) \end{equation*}\]

则称\(f\)是凸函数(convex function),若上式取严格的小于号则称\(f\)严格凸(strictly convex)。若\(-f\)是凸函数,则称\(f\)为凹函数(concave function)。若\(-f\)严格凸,则称\(f\)严格凹(strictly concave)。

Property 8.1.10. 设\(f\)是\(E\subseteq\mathbb{R}^{n}\)上的凸函数。凸函数具有如下性质:

  1. 若\(A\in M_{n}(\mathbb{R}),\;b\in\mathbb{R}^{n}\),对任意的\(x\in E\)有\(Ax+b\in E\),则\(f(Ax+b)\)也是\(E\)上的凸函数;

  2. 若\(f\)在\(E\)上Fréchet可微,根据性质 4.1.2(1),则\(f\)是\(E\)上的凸函数的充分必要条件为对任意的\(x,y\in E\)有23严格的证明:

    \[\begin{equation*} f(y)\geqslant f(x)+(\nabla fx,y-x) \end{equation*}\]

    \(f\)是\(E\)上的严格凸函数的充分必要条件为对任意不等的\(x,y\in E\)有:

    \[\begin{equation*} f(y)>f(x)+(\nabla fx,y-x) \end{equation*}\]

  3. 若\(f\)在\(E\)上二次可微,则\(f\)是凸函数的充分必要条件为对任意的\(x,h\in E\)有\(\operatorname{D}^2fx(h,h)\geqslant0\),对非零的\(h\)有\(\operatorname{D}^2fx(h,h)>0\)是\(f\)严格凸的充分条件;

  4. 若\(f\)是凸函数,\(t\in(a,b)\),则存在经过点\((t,f(t))\)的直线\(g(x)=c(x-t)+f(t)\)满足对任意的\(x\in(a,b)\)有\(g(x)\leqslant f(x)\)。当\(f\)为严格凸函数时,\(g(x)=f(x)\)当且仅当\(x=t\)。

证明. (1)任取\(x_1,x_2\in E\),则对任意的\(\alpha\in[0,1]\)有:

\[\begin{align*} &f\{A[\alpha x_1+(1-\alpha)x_2]+b\}=f[A\alpha x_1+\alpha b+A(1-\alpha)x_2+(1-\alpha)b] \\ =&f[\alpha(Ax_1+b)+(1-\alpha)(Ax_2+b)]\leqslant\alpha f(Ax_1+b)+(1-\alpha)f(x_2+b) \end{align*}\]

(2)必要性:对任意的\(x,y\in E\),因为\(f\)是凸函数,由(1)可知:

\[\begin{equation*} \varphi(t)=f[x+t(y-x)] \end{equation*}\]

在\([0,1]\)上是凸函数,所以对任意的\(t\in(0,1]\)有:

\[\begin{equation*} \varphi(t)\leqslant(1-t)\varphi(0)+t\varphi(1) \end{equation*}\]

即:

\[\begin{equation*} \frac{\varphi(t)-\varphi(0)}{t}\leqslant\varphi(1)-\varphi(0) \end{equation*}\]

因为\(f\)在\(x\)处Fréchet可微,由性质 4.1.2(1)可知\(\varphi\)在\(t=0\)处可微,于是取\(t\downarrow0\),根据性质 3.3.3(4)可知:

\[\begin{equation*} \varphi'(0)=\operatorname{D}_Gfx(y-x)=(\nabla fx,y-x)\leqslant f(y)-f(x) \end{equation*}\]

充分性:任取\(x,y\in E\)和\(\alpha\in[0,1]\),记\(z=\alpha x+(1-\alpha)y\),由条件可得:

\[\begin{equation*} f(x)\geqslant f(z)+(\nabla fz,x-z),\quad f(y)\geqslant f(z)+(\nabla fz,y-z) \end{equation*}\]

即:

\[\begin{equation*} \alpha f(x)\geqslant\alpha f(z)+\alpha(\nabla fz,x-z),\quad(1-\alpha)f(y)\geqslant(1-\alpha)f(z)+(1-\alpha)(\nabla fz,y-z) \end{equation*}\]

所以:

\[\begin{equation*} \alpha f(x)+(1-\alpha)f(y)\geqslant f(z)+\alpha(\nabla fz,x-z)+(1-\alpha)(\nabla fz,y-z) \end{equation*}\]

注意到:

\[\begin{align*} &\alpha(\nabla fz,x-z)+(1-\alpha)(\nabla fz,y-z)=\Big(\nabla fz,\alpha(x-z)+(1-\alpha)(y-z)\Big) \\ =&\Big(\nabla fz,\alpha x+(1-\alpha)y-\alpha z-(1-\alpha)z\Big)=\Big(\nabla fz,z-\alpha z-(1-\alpha)z\Big) \\ =&\Big(\nabla fz,\mathbf{0}\Big)=0 \end{align*}\]

所以:

\[\begin{equation*} \alpha f(x)+(1-\alpha)f(y)\geqslant f[\alpha x+(1-\alpha)y] \end{equation*}\]

由\(x,y,\alpha\)的任意性可知\(f\)是\(E\)上的凸函数。

(3)必要性:对任意的\(x\in E\),若存在\(h\in E\)使得\(\operatorname{D}^2f(h,h)<0\),由定理 4.6(1)可得:

\[\begin{equation*} f(x_0+th)=f(x_0)+t\operatorname{D}fx(h)+\frac{t^2}{2}\operatorname{D}^2fx(h,h)+\operatorname{o}(t^2),\;(t\to0) \end{equation*}\]

即:

\[\begin{equation*} \frac{f(x_0+th)-f(x_0)-t\operatorname{D}fx(h)}{t^2}=\frac{1}{2}\operatorname{D}^2fx(h,h)+\frac{\operatorname{o}(t^2)}{t^2},\;(t\to0) \end{equation*}\]

由性质 3.3.3(5.b)(3)可知当\(t\)足够小时有:

\[\begin{equation*} \frac{f(x_0+th)-f(x_0)-t(\operatorname{D}fx)h}{t^2}<0 \end{equation*}\]

根据性质 4.1.2(1)可知这与(2)矛盾。

充分性:由定理 4.6(2)和(2)立即可得。 ◻

8.2 最优性理论

8.2.0.1 可微函数的最优性

由定理 4.7给出。

8.2.0.2 不可微函数的最优性

Definition 8.22. 设\(f\)关于\(\mathbb{R}^{n}\)是适当函数,\(x\in\mathbb{R}^{n}\)。若\(g\in\mathbb{R}^{n}\)对任意的\(y\in\mathbb{R}^{n}\)满足:

\[\begin{equation*} f(y)\geqslant f(x)+g^T(y-x) \end{equation*}\]

则称\(g\)为\(f\)在\(x\)处的一个次梯度(subgradient),称\(f\)在\(x\)处所有次梯度构成的集合为次微分(subdifferential)。

8.2.1 映射的上下极限与半连续性

Definition 8.23. \((X,\rho)\)是度量空间,\(E\)是\(X\)的子空间,\(f\)是\(E\)到\(\mathbb{R}\)上的映射。对于\(E\)中的任一聚点\(a\),定义:

\[\begin{gather*} \liminf_{x\to a}f(x)=\lim_{\varepsilon\to0}\Bigl(\inf\{f(x):x\in U(a,\varepsilon)\}\setminus\{a\}\Bigr) \\ \limsup_{x\to a}f(x)=\lim_{\varepsilon\to0}\Bigl(\sup\{f(x):x\in U(a,\varepsilon)\}\setminus\{a\}\Bigr) \end{gather*}\]

Definition 8.24. \((X,\rho)\)是度量空间,\(f\)是\(X\)到\(\overline{\mathbb{R}}\)上的映射。若:

\[\begin{equation*} \limsup_{x\to a}f(x)\leqslant f(x_0) \end{equation*}\]

则称\(f(x)\)在\(a\)点上半连续(upper semicontinuous)。 若:

\[\begin{equation*} \liminf_{x\to a}f(x)\leqslant f(x_0) \end{equation*}\]

则称\(f(x)\)在\(a\)点下半连续(lower semicontinuous)。

8.2.2 线搜索法

在迭代优化算法中,一次迭代通常包含两个基本步骤:首先确定一个能够使损失函数下降的搜索方向,其次确定沿该方向移动的步长。搜索方向决定参数更新的方向,而步长决定参数沿该方向移动的距离。若步长过小,则每次迭代的推进幅度有限,算法的收敛速度可能较慢;若步长过大,则可能越过适宜的更新区域,使损失函数值无法下降,甚至反而增大。因此,步长的选择需要在保证损失函数充分下降与维持足够迭代进展之间取得平衡。

线搜索法正是为了解决这一问题而引入的。所谓线搜索法,是指在当前点\(x_n\)和搜索方向\(d_n\)确定之后,沿着射线\(\{x_n+\alpha d_n:\alpha>0\}\)寻找一个合适的步长\(\alpha_n\),并用它更新迭代点。

TikZ 数学图示
图 8.3 沿固定下降方向的一维函数。步长过小会使每次下降有限;步长过大可能越过低点,甚至使函数值高于当前点。

Definition 8.25. 给定初始步长\(\alpha>0\)和收缩因子\(\rho\in(0,1)\),定义候选步长序列:

\[\begin{equation*} \alpha_n=\rho^n\alpha,\quad n\in\mathbb{N} \end{equation*}\]

从\(n=0\)开始依次检验\(\alpha_n\)是否满足预先给定的步长接受准则。若取第一个满足该准则的步长\(\alpha_n\)作为当前迭代步长,则称这种步长选择方法为回溯法(backtracking method)。

Definition 8.26. 设\(X\)为实赋范线性空间,\(E\subseteq X\)为开集,函数\(f:E\to\mathbb{R}\)在\(x\in E\)处Fréchet可微,\(c\in(0,1)\),\(d\in X\)是一个下降方向。若\(\alpha>0\)满足\(x+\alpha d\in E\),并且

\[\begin{equation*} f(x+\alpha d)\leqslant f(x)+c\alpha\operatorname{D}fx(d) \end{equation*}\]

则称\(\alpha\)满足沿方向\(d\)的Armijo准则。

note 8.2. 因为\(d\in X\)为下降方向,所以\(\operatorname{D}fx(d)<0\)。由定理 4.6(1)可得:

\[\begin{equation*} f(x+\alpha d)=f(x)+\alpha\operatorname{D}fx(d)+\operatorname{o}(\alpha),\;(\alpha\to0) \end{equation*}\]

因此一阶线性近似预示着,沿方向\(d\)移动一个小步长\(\alpha\)后,函数值应当下降,其线性预测下降量为\(-\alpha\operatorname{D}fx(d)>0\)。Armijo准则要求步长\(\alpha\)满足:

\[\begin{equation*} f(x)-f(x+\alpha d)\geqslant-c\alpha\operatorname{D}fx(d) \end{equation*}\]

因此,Armijo准则的含义是实际下降量至少要达到一阶Taylor线性预测下降量的\(c\)倍。参数\(c\in(0,1)\)控制这个要求的严格程度,\(c\)越大,要求实际下降越接近线性预测下降量。在应用中,\(c\)常取较小的数值,如\(c=10^{-4}\)。

Armijo准则的缺陷在于它只控制了步长\(\alpha\)不能太大,因为太大的步长通常会违反下降条件,但它没有控制”步长不能太小”。如果步长非常非常小,函数值虽然下降得很少,但也可能满足Armijo准则。

Definition 8.27. 设\(X\)为实赋范线性空间,\(E\subseteq X\)为开集,函数\(f:E\to\mathbb{R}\)在\(x\in E\)处Fréchet可微,\(c\in(0,0.5)\),\(d\in X\)是一个下降方向。若\(\alpha>0\)满足\(x+\alpha d\in E\),并且

\[\begin{equation*} f(x)+(1-c)\alpha\operatorname{D}fx(d) \leqslant f(x+\alpha d) \leqslant f(x)+c\alpha\operatorname{D}fx(d), \end{equation*}\]

则称\(\alpha\)满足沿方向\(d\)的Goldstein准则。

note 8.3. Goldstein准则可以看作是在Armijo准则的基础上进一步加入了对步长下界的控制。

Goldstein准则中的右侧不等式正是Armijo准则,它要求函数值具有充分下降。Goldstein准则中的左侧不等式为:

\[\begin{equation*} f(x)-f(x+\alpha d) \leqslant -(1-c)\alpha\operatorname{D}fx(d) \end{equation*}\]

由于\(-\alpha\operatorname{D}fx(d)>0\)是一阶Taylor线性近似预示的下降量,因此左侧不等式要求实际下降量不能超过线性预测下降量的\((1-c)\)倍。由定理 4.6(1)可得:

\[\begin{equation*} f(x+\alpha d)=f(x)+\alpha\operatorname{D}fx(d)+\operatorname{o}(\alpha),\;(\alpha\to0) \end{equation*}\]

因此当\(\alpha\)过小时,实际下降量会过于接近完整的一阶线性预测下降量,从而违反Goldstein左侧条件。直观地说,这表示沿方向\(d\)还没有走得足够远,函数沿该方向的下降趋势尚未明显减弱。

综上,Goldstein准则通过两个不等式同时控制步长:右侧的Armijo条件保证步长不能太大,左侧条件保证步长不能太小。

Goldstein的缺陷出于其优势,它用Taylor展开的原理保证了所取的步长\(\alpha\)不会太小,但是对应的新增不等式在逻辑上却对应着要求下降量不能太大,导致它可能会错杀好步长。

Definition 8.28. 设\(X\)为实赋范线性空间,\(E\subseteq X\)为开集,函数\(f:E\to\mathbb{R}\)在\(x\)与\(x+\alpha d\)处Fréchet可微,\(0<c_1<c_2<1\),\(d\in X\)是一个下降方向。若\(\alpha>0\)满足\(x+\alpha d\in E\),并且:

\[\begin{gather*} f(x+\alpha d)\leqslant f(x)+c_1\alpha\operatorname{D}fx(d),\\ \operatorname{D}f(x+\alpha d)(d)\geqslant c_2\operatorname{D}fx(d) \end{gather*}\]

则称\(\alpha\)满足沿方向\(d\)的Wolfe准则。

note 8.4. Wolfe准则可以看作是对Goldstein准则的改进。它保留了Armijo准则中的充分下降条件,但将Goldstein准则中用于排除过小步长的函数值下界替换为了方向导数条件。

Wolfe准则中的第二个不等式称为曲率条件。由于\(\operatorname{D}fx(d)<0\),该条件要求新点\(x+\alpha d\)处沿方向\(d\)的方向导数不能过小。换言之,沿方向\(d\)的下降趋势必须已经有所减弱。因此,曲率条件的作用是排除过小的步长。

Wolfe准则的缺陷在于所取的\(\alpha\)可能导致\(x+\alpha d\)的下降趋势反向增长得过强,即已经从另一个方向远离了损失函数的极小点。

Definition 8.29. 设\(X\)为实赋范线性空间,\(E\subseteq X\)为开集,函数\(f:E\to\mathbb{R}\)在\(x\)与\(x+\alpha d\)处Fréchet可微,\(0<c_1<c_2<1\),\(d\in X\)是一个下降方向。若\(\alpha>0\)满足\(x+\alpha d\in E\),并且:

\[\begin{gather*} f(x+\alpha d)\leqslant f(x)+c_1\alpha\operatorname{D}fx(d) \\ \left|\operatorname{D}f(x+\alpha d)(d)\right|\leqslant c_2\left|\operatorname{D}fx(d)\right| \end{gather*}\]

则称\(\alpha\)满足沿方向\(d\)的强Wolfe准则。

note 8.5. 强Wolfe准则是Wolfe准则的加强形式,其第二个不等式要求沿方向\(d\)的下降趋势不仅不能过强,而且也不能反向增长得过强。

Theorem 8.3. 设\(X\)为实Hilbert空间,\(E\subseteq X\)为开集,函数\(f\in C^1:E\to\mathbb{R}\),\(f\)在\(E\)上有下界,\(\nabla f\)在\(E\)上Lipschitz连续,即存在常数\(L>0\),使得:

\[\begin{equation*} ||\nabla fx-\nabla fy||\leqslant L||x-y|| \end{equation*}\]

给定初始点\(x_0\in E\),考虑由线搜索方法生成的迭代序列\(\{x_{n+1}=x_n+\alpha_nd_n\}\),其中\(d_n\)是\(f\)在\(x_n\)处的下降方向。若步长\(\alpha_k>0\)满足Wolfe准则,则有:

\[\begin{equation*} \sum_{n=0}^{+\infty}\cos^2\theta_n||\nabla fx_n||^2<+\infty \end{equation*}\]

其中\(\theta_n\)为\(-\nabla fx_n\)与\(d_n\)之间的夹角,即:

\[\begin{equation*} \cos\theta_n=\frac{-(\nabla fx_n,d_n)}{||\nabla fx_n||\;||d_n||} \end{equation*}\]

证明. 由Wolfe准则可知:

\[\begin{equation*} (\nabla fx_{n+1},d_n)\geqslant c_2(\nabla fx_n,d_n) \end{equation*}\]

所以:

\[\begin{equation*} (\nabla fx_{n+1}-\nabla fx_n,d_n)\geqslant(c_2-1)(\nabla fx_n,d_n) \end{equation*}\]

根据不等式 1可知:

\[\begin{equation*} (\nabla fx_{n+1}-\nabla fx_n,d_n)\leqslant||\nabla fx_{n+1}-\nabla fx_n||\;||d_n||\leqslant L||x_{n+1}-x_n||\;||d_n|=\alpha_nL||d_n||^2 \end{equation*}\]

联合上两式可得:

\[\begin{equation*} \alpha_nL||d_n||^2\geqslant(c_2-1)(\nabla fx_n,d_n) \\ \alpha_n\geqslant\frac{(c_2-1)(\nabla fx_n,d_n)}{L||d_n||^2} \end{equation*}\]

由Wolfe准则可知:

\[\begin{equation*} f(x_{n+1})\leqslant f(x_n)+c_1\alpha_n(\nabla fx_n,d_n)\leqslant f(x_n)+c_1\frac{(c_2-1)(\nabla fx_n,d_n)}{L||d_n||^2}(\nabla fx_n,d_n) \end{equation*}\]

即:

\[\begin{equation*} f(x_{n+1})\leqslant f(x_n)+c_1\frac{(c_2-1)(\nabla fx_n,d_n)^2}{L||d_n||^2}=f(x_n)+c_1\frac{(c_2-1)\cos^2\theta_n}{L}||\nabla fx_n||^2 \end{equation*}\]

对上式求和即可得到:

\[\begin{equation*} f(x_{n+1})\leqslant f(x_0)-\frac{c_1(1-c_2)}{L}\sum_{i=0}^{n}\cos^2\theta_i||\nabla fx_i||^2 \end{equation*}\]

因为\(0<c_1<c_2<1\),所以\(c_1(1-c_2)>0\),即对任意的\(n\in\mathbb{N}\)有\(f(x_{n+1})\leqslant f(x_0)\)。因为\(f\)在\(E\)上有下界,所以结论成立。 ◻

Theorem 8.4. 设\(X\)为实Hilbert空间,\(E\subseteq X\)为开集,函数\(f\in C^1(E):E\to\mathbb{R}\),\(f\)在\(E\)上有下界,\(\nabla f\)在\(E\)上Lipschitz连续,即存在常数\(L>0\),使得:

\[\begin{equation*} ||\nabla fx-\nabla fy||\leqslant L||x-y|| \end{equation*}\]

给定初始点\(x_0\in E\),考虑由线搜索方法生成的迭代序列\(\{x_{n+1}=x_n+\alpha_nd_n\}\),其中\(d_n\)是\(f\)在\(x_n\)处的下降方向。若步长\(\alpha_k>0\)满足Wolfe准则,且对任意的\(n\in\mathbb{N}\),存在\(\beta>0\)使得:

\[\begin{equation*} \theta_n<\frac{\pi}{2}-\beta \end{equation*}\]

则:

\[\begin{equation*} \lim_{n\to\infty}||\nabla fx_n||=0 \end{equation*}\]

证明. 若结论不成立,由性质 3.1.1(4)可知存在\(\{\nabla fx_n\}\)的子列\(\{\nabla fx_{n_k}\}\)和\(\varepsilon>0\)使得:

\[\begin{equation*} \forall\;k\in\mathbb{N}^+,\;||\nabla fx_{n_k}||>\varepsilon \end{equation*}\]

注意到:

\[\begin{equation*} \cos(\theta_n)>\cos\left(\frac{\pi}{2}-\beta\right)=\sin\beta \end{equation*}\]

根据定理 8.3可得:

\[\begin{equation*} +\infty>\sum_{n=0}^{+\infty}\cos^2\theta_n||\nabla fx_n||^2>\sum_{k=0}^{+\infty}\cos^2\theta_{n_k}||\nabla fx_{n_k}||^2>\sum_{k=0}^{+\infty}\sin^2\beta\varepsilon^2 \end{equation*}\]

矛盾,所以结论成立。 ◻

note 8.6. 上述定理中的条件是很直观的,如果\(\theta_n\)很接近\(\dfrac{\pi}{2}\),就意味着梯度方向和下降方向垂直,根据定理 4.6(1)可知此此时损失函数值几乎不发生变化。

8.2.3 梯度下降法

Method 8.1. (Gradient Descent)

设\(X\)为实Hilbert空间,开集\(E\subseteq X\),函数\(f:E\to\mathbb{R}\)在\(E\)上Gâteaux可微,任取初始点\(x_0\in E\)。在第\(n\)次迭代中计算梯度\(\nabla fx_n\);若\(\nabla fx_n=\mathbf{0}\),则停止迭代,否则选择步长\(\alpha_n>0\)使得\(x_n-\alpha_n\nabla fx_n\in E\),并令:

\[\begin{equation*} x_{n+1}=x_n-\alpha_n\nabla fx_n \end{equation*}\]

如此反复更新得到迭代点列\(\{x_n\}\),这种方法称为梯度下降(gradient descent),步长\(\alpha_n\)被称为学习率(learning rate, LR)。

note 8.7. 由性质 4.1.4可知,在所有单位方向中负梯度方向:

\[\begin{equation*} -\frac{\nabla fx_n}{||\nabla fx_n||} \end{equation*}\]

使\(f\)在\(x_n\)处的方向导数最小,由方向导数的定义可知它是当前点处使得函数值下降最快的单位方向。

8.2.4 Newton法

Method 8.2. (Newton Method)

设\(E\subseteq \mathbb{R}^m\)为开集,函数\(f\in C^2(E):E\to\mathbb{R}\)。给定当前迭代点\(x_n\in E\),\(\nabla^2fx_n\)可逆,根据定理 4.6(1),考虑在\(x_n\)处沿方向\(d_n\)的二阶Taylor近似:

\[\begin{align*} f(x_n+d_n)&\approx f(x_n)+\operatorname{D}fx_n(d_n)+\frac{1}{2}\operatorname{D}^2fx_n(d_n,d_n) \\ &=f(x_n)+(\nabla fx_n)^Td_n+\frac{1}{2}d_n^T\nabla^2fx_nd_n \end{align*}\]

Newton法的基本思想是:不直接最小化原函数\(f(x_n+d_n)\),而是在当前点\(x_n\)处求解方向\(d_n\)使得其二次近似最小化。因此,根据定理 4.7(1),令上右式关于\(d_n\)的一阶Fréchet微分为零可得:

\[\begin{equation*} \nabla^2fx_nd_n=-\nabla fx_n \end{equation*}\]

从而得到:

\[\begin{equation*} d_n=-\left(\nabla^2fx_n\right)^{-1}\nabla fx_n \end{equation*}\]

于是Newton法的迭代格式为:

\[\begin{equation*} x_{n+1}=x_n+d_n=x_n-\left(\nabla^2fx_n\right)^{-1}\nabla fx_n \end{equation*}\]

其中\(d_n\)称为\(f\)在\(x_n\)处的Newton方向。

Theorem 8.5. 设\(E\subseteq \mathbb{R}^m\)为开集,函数\(f\in C^2(E):E\to\mathbb{R}\),\(f\)在\(E\)上的最小值点\(x\)满足:

\[\begin{equation*} \nabla fx=\mathbf{0},\quad\nabla^2fx>0 \end{equation*}\]

\(\nabla^2fx\)在\(x\)的\(\delta\)邻域(\(\delta>0\))\(U(x,\delta)\)上满足Lipschitz条件:

\[\begin{equation*} ||\nabla^2fx-\nabla^2fy||\leqslant L||x-y||,\quad\forall x,y\in U(x,\delta) \end{equation*}\]

考虑Newton迭代:

\[\begin{equation*} x_{n+1}=x_n+d_n=x_n-\left(\nabla^2fx_n\right)^{-1}\nabla fx_n \end{equation*}\]

则:

  1. 如果初始点\(x_0\)充分接近\(x\),则Newton迭代良定,并且\(\{x_n\}\)Q-二次收敛到\(x\);

  2. \(\{||\nabla fx_n||\}\)Q-二次收敛到\(0\)。

证明. (1)注意到:

\[\begin{align*} x_{n+1}-x&=x_n-\left(\nabla^2fx_n\right)^{-1}\nabla fx_n-x=\left(\nabla^2fx_n\right)^{-1}\left[\nabla^2fx_n(x_n-x)-\nabla fx_n\right] \\ &=\left(\nabla^2fx_n\right)^{-1}\left[\nabla^2fx_n(x_n-x)-\left(\nabla fx_n-\nabla fx\right)\right] \end{align*}\]

由Newton-Leibniz公式、性质 4.1.1(6)和性质 5.4.3(6)可得:

\[\begin{equation*} \nabla fx_n-\nabla fx=\int_{0}^{1}\nabla^2f[x_n+t(x-x_n)](x_n-x)\mathop{}\!\mathrm{d}t \end{equation*}\]

根据性质 5.4.3(6)(7)积分的范数小于等于范数的积分,算子范数可得:

\[\begin{align*} &\left\|\left[\nabla^2fx_n(x_n-x)-\left(\nabla fx_n-\nabla fx\right)\right]\right\| \\ =&\left\|\left[\nabla^2fx_n(x_n-x)-\int_{0}^{1}\nabla^2f[x_n+t(x-x_n)](x_n-x)\mathop{}\!\mathrm{d}t\right]\right\| \\ =&\left\|\int_{0}^{1}\left\{\nabla^2fx_n-\nabla^2f[x_n+t(x-x_n)]\right\}(x_n-x)\mathop{}\!\mathrm{d}t\right\| \\ \leqslant&\int_{0}^{1}\left\|\left\{\nabla^2fx_n-\nabla^2f[x_n+t(x-x_n)]\right\}(x_n-x)\right\|\mathop{}\!\mathrm{d}t \\ \leqslant&\int_{0}^{1}\left\|\left\{\nabla^2fx_n-\nabla^2f[x_n+t(x-x_n)]\right\}\right\|||x_n-x||\mathop{}\!\mathrm{d}t \\ \leqslant&||x_n-x||\int_{0}^{1}Lt||x-x_n||\mathop{}\!\mathrm{d}t=\frac{||x_n-x||^2L}{2} \end{align*}\]

于是由算子范数可得:

\[\begin{align*} &||x_{n+1}-x||=\left\|\left(\nabla^2fx_n\right)^{-1}\left[\nabla^2fx_n(x_n-x)-\left(\nabla fx_n-\nabla fx\right)\right]\right\| \\ \leqslant&\left\|\left(\nabla^2fx_n\right)^{-1}\right\|\left\|\left[\nabla^2fx_n(x_n-x)-\left(\nabla fx_n-\nabla fx\right)\right]\right\| \\ \leqslant&\left\|\left(\nabla^2fx_n\right)^{-1}\right\|\frac{||x_n-x||^2L}{2} \end{align*}\]

可以证明存在\(r>0\),当\(||x_n-x||\leqslant r\)时有\(\left\|\left(\nabla^2fx_n\right)^{-1}\right\|\leqslant2\left\|\left(\nabla^2fx\right)^{-1}\right\|\)1,所以:

\[\begin{equation*} \frac{||x_{n+1}-x||}{||x_n-x||^2}\leqslant L\left\|\left(\nabla^2fx\right)^{-1}\right\| \end{equation*}\]

当:

\[\begin{equation*} ||x_0-x||\leqslant\min\left\{\delta,r,\frac{1}{2L\left\|\left(\nabla^2fx\right)^{-1}\right\|}\right\} \end{equation*}\]

时结论成立,第三式是因为:

\[\begin{equation*} ||x_1-x||\leqslant L\left\|\left(\nabla^2fx\right)^{-1}\right\|||x_0-x||^2 \end{equation*}\]

右式未必满足小于等于\(||x_0-x||\),这将导致序列\(\{x_n:n\geqslant1\}\)不在\(U(x,\delta)\)中。一旦控制了该不等式,则\(||x_1-x||\leqslant\dfrac{1}{2}||x_0-x||\)。设\(n\)时有:

\[\begin{equation*} ||x_n-x||\leqslant\frac{1}{2L\left\|\left(\nabla^2fx\right)^{-1}\right\|} \end{equation*}\]

则\(n+1\)时有:

\[\begin{equation*} ||x_{n+1}-x||\leqslant L\left\|\left(\nabla^2fx\right)^{-1}\right\|||x_n-x||^2\leqslant\frac{1}{2}||x_n-x|| \end{equation*}\]

这保证了\(\{x_n\}\subseteq U(x,\delta)\),并且每进行一步范数误差至少减少一半。

(2)由Newton方程、Newton-Leibniz公式,积分的范数小于等于范数的积分,算子范数、性质 4.1.1(6)可得:

\[\begin{align*} &||\nabla fx_{n+1}||=||\nabla fx_{n+1}-\nabla fx_n-\nabla^2fx_nd_n\||=\left\|\int_{0}^{1}\nabla^2f(x_n+td_n)d_n\mathop{}\!\mathrm{d}t-\nabla^2fx_nd_n\right\| \\ =&\left\|\int_{0}^{1}\left[\nabla^2f(x_n+td_n)-\nabla^2fx_n\right]d_n\mathop{}\!\mathrm{d}t\right\|\leqslant\int_{0}^{1}\left\|\left[\nabla^2f(x_n+td_n)-\nabla^2fx_n\right]d_n\right\|\mathop{}\!\mathrm{d}t \\ \leqslant&\int_{0}^{1}\left\|\nabla^2f(x_n+td_n)-\nabla^2fx_n\right\|||d_n||\mathop{}\!\mathrm{d}t\leqslant\int_{0}^{1}L||td_n||\;||d_n||\mathop{}\!\mathrm{d}t=\frac{L||d_n||^2}{2} \\ =&\frac{L}{2}\left\|\left(\nabla^2fx_n\right)^{-1}\nabla fx_n\right\|^2\leqslant\frac{L}{2}\left\|\left(\nabla^2fx_n\right)^{-1}\right\|^2||\nabla fx_n||^2\leqslant2L\left\|\left(\nabla^2fx\right)^{-1}\right\|^2||\nabla fx_n||^2 \end{align*}\]

所以\(\{||\nabla fx_n||\}\)Q-二次收敛到\(0\)。 ◻

8.3 启发式算法

前面介绍的梯度下降法、牛顿法等方法都属于基于局部信息的确定性优化方法。它们利用目标函数在当前点附近的一阶或二阶导数信息构造搜索方向,并通过不断迭代寻找目标函数的极小点。当优化问题具有较好的数学性质时(如目标函数可微、是凸函数等),这类方法具有清晰的理论基础,并且往往能够获得较好的收敛性质。

然而,在实际优化问题中,我们经常会遇到更加复杂的情形。例如,目标函数可能是非凸的,从而存在大量局部极小点;目标函数可能不可微;决策变量还可能是离散的,使梯度本身失去意义。此时,仅依赖局部导数信息的优化方法会受到明显限制。

这促使我们考虑一种不同于经典数值优化的方法,也即启发式算法(heuristic algorithm)。启发式算法则不再要求每一步搜索都由梯度、Hessian矩阵等解析信息所决定,而是通过人为设计的搜索机制,在解空间中不断产生、评价和更新候选解。其核心不再是如何精确地确定当前点的下降方向,而是如何在有限的计算资源下有效地搜索整个解空间。

因此,启发式算法所面临的一个基本问题是如何协调探索与利用:一方面,需要充分利用当前已经发现的优质解,在其附近进行更细致的搜索;另一方面,又必须保留探索其他区域的能力,以避免搜索过程过早集中在某个局部最优解附近。不同的启发式算法,本质上就是通过不同的搜索机制来实现这两者之间的平衡。

note 8.8. 启发式算法中有相当一部分方法通过模仿自然界中的某种现象来构造搜索过程,例如生物进化、群体协作、觅食行为、物理退火过程等。这种思想本身并无问题:一个自然现象如果能够抽象出具有明确作用的搜索机制,例如随机扰动、信息共享、选择、记忆或探索与利用之间的动态平衡,那么它完全可以成为设计优化算法的有效启发。

但是,本人对于近年来大量所谓的”动物园算法”持较为抵触的态度。其典型模式是选择某种动物、植物或自然现象作为新的叙事背景,再将若干已有的随机搜索、局部扰动、个体更新和最优解保留机制重新组合,随后赋予其新的生物学名称,并据此宣称提出了一种新的优化算法。随着这种研究范式不断扩张,文献中逐渐出现了大量以不同动物行为命名、但在数学结构和搜索机制上高度相似的算法。

从优化的角度看,一个启发式算法是否有价值,并不取决于它所模拟的自然现象是否新颖,而应当取决于它是否引入了新的搜索机制,是否能够解释这种机制为何有效,以及它相对于已有方法究竟解决了什么问题。因此,本节无意对层出不穷的自然启发式算法进行系统罗列,而只选择少数已经被长期研究、广泛应用,并在大量实际问题与工业实践中反复验证有效的经典方法。介绍这些算法的目的也并非记忆某一种具体的自然界类比,而是借助它们理解启发式优化中若干具有普遍意义的思想。

8.3.1 模拟退火算法

模拟退火(simulated annealing, SA)是最具代表性的启发式优化方法之一。它所体现的核心思想并不复杂:如果搜索过程始终只接受能够立即改善目标函数的移动,那么算法很容易被局部最优解困住;因此,在搜索的早期阶段,应当允许算法以一定概率接受较差的解,以换取离开当前局部区域、探索其他潜在优质区域的机会,而随着搜索的进行,再逐渐降低这种冒险的程度,使算法最终转向更加稳定的局部改进。

这一思想来自统计物理中的退火过程。高温状态下,粒子具有较强的热运动能力,系统可以在不同能量状态之间频繁跃迁;随着温度缓慢降低,系统逐渐趋向低能量状态。模拟退火将目标函数值类比为系统能量,将温度引入优化过程,并通过一个随温度变化的概率机制决定是否接受使目标函数上升的候选解。

模拟退火真正值得关注的是它给出了一个非常典型的搜索设计原则:通过受控地接受劣化解,在全局探索能力与局部优化能力之间建立动态平衡。高温阶段强调探索,允许搜索越过局部极小值之间的障碍;低温阶段强调利用,使搜索逐渐集中于已经发现的优质区域。温度的下降过程实际上控制着算法从”广泛探索”向”局部收敛”的转变。

从这个角度看,模拟退火可以被理解为对局部搜索的一种随机化扩展。其基本框架由三个核心部分构成:候选解的生成机制、劣化解的接受机制以及温度的下降策略。下面将依次介绍这三个部分,并在此基础上给出模拟退火算法的完整形式。

8.3.1.1 候选解的生成机制

考虑优化问题:

\[\begin{equation*} \min_{x\in D}\;f(x) \end{equation*}\]

其中\(D\)为非空解空间,\(f:D\to\mathbb{R}^{}\)为目标函数。模拟退火每个\(x\in D\)指定一个邻域\(U(x)\subseteq D\),并在当前解\(x_n\)处按照概率分布\(q_n(\cdot|x_n)\)生成候选解\(y_n\in U(x_n)\)。

这里的概率分布\(q_n(\cdot|x)\)并不是由模拟退火算法本身固定给出的,而需要使用者根据具体问题的解空间结构以及所采用的邻域定义自行设计。其作用是规定在给定当前解\(x\)后,如何在邻域\(U(x)\)中随机选择下一步的候选解。不同的候选解生成机制对应不同的\(q_n(\cdot|x)\),并会直接影响算法对解空间的搜索方式。通常希望\(q_n(\cdot|x)\)能够使邻域中的候选解具有充分的被访问机会,从而避免搜索被人为限制在解空间的某个局部区域。

8.3.1.2 接受机制

在第\(n\)次迭代中,设当前解为\(x_n\),当前温度为\(T_n>0\),并由候选解生成机制得到候选解\(y_n\)。接受机制决定候选解是否成为下一次迭代的当前解:若\(y_n\)被接受,则令\(x_{n+1}=y_n\);否则保持\(x_{n+1}=x_n\)。记此次移动所引起的目标函数变化为:

\[\begin{equation*} \Delta_n=f(y_n)-f(x_n) \end{equation*}\]

当\(\Delta_n\leqslant0\)时,候选解不劣于当前解,因此令\(x_{n+1}=y_n\);当\(\Delta_n>0\)时,候选解虽然使目标函数增大,但仍以概率:

\[\begin{equation*} \exp\left(-\frac{\Delta_n}{T_n}\right) \end{equation*}\]

被接受。因此,接受概率可以统一写为:

\[\begin{equation*} a_{T_n}(x_n,y_n)=\min\left\{1,\exp\left[-\frac{f(y_n)-f(x_n)}{T_n}\right]\right\} \end{equation*}\]

具体实现时,可以独立生成\(u_n\sim\operatorname{U}(0,1)\),并按照:

\[\begin{equation*} x_{n+1}= \begin{cases} y_n,&u_n\leqslant a_{T_n}(x_n,y_n) \\ x_n,&u_n>a_{T_n}(x_n,y_n) \end{cases} \end{equation*}\]

更新当前解。

当\(\Delta_n>0\)固定时,\(\exp(-\Delta_n/T)\)关于\(T\)单调递增,并且:

\[\begin{equation*} \lim_{T\to+\infty}\exp\left(-\frac{\Delta_n}{T}\right)=1,\quad\lim_{T\downarrow0}\exp\left(-\frac{\Delta_n}{T}\right)=0 \end{equation*}\]

因此,在高温阶段,算法对劣化移动具有较高的容忍度,可以较自由地离开当前区域并探索其他可能的优质区域;随着温度降低,劣化移动被接受的概率逐渐减小,搜索越来越集中于已经发现的优质区域。在相同温度下,目标函数增加得越多,候选解被接受的概率越小,因此这种随机接受也并非无差别地允许劣化,而是更加倾向于接受程度较小的劣化移动。由此,温度实际上控制着算法在探索与利用之间的权衡:较高的温度强调探索,较低的温度强调利用,而逐渐降温的过程则使搜索从广泛探索平稳过渡到局部精细搜索。

8.3.1.3 温度的下降策略

温度下降策略只需要保证温度随迭代不断降低,并逐渐趋近于\(0\),即:

\[\begin{equation*} T_{k+1}<T_k,\qquad \lim_{k\to+\infty}T_k=0 \end{equation*}\]

常用的温度变化形式包括:

\[\begin{gather*} T_k=T_0-\alpha k,\quad\alpha>0 \\ T_k=T_0\rho^k,\quad0<\rho<1 \\ T_k=\frac{T_0}{1+\alpha k},\quad\alpha>0 \\ T_k=\frac{T_0}{(1+\alpha k)^\gamma},\quad\alpha,\gamma>0 \\ T_k=\frac{c}{\log(k+k_0)},\quad c>0,\;k_0>1 \end{gather*}\]

8.3.1.4 完整算法

Algorithm 8.1 Simulated Annealing

Input: Feasible set \(D\); objective function \(f:D\to\mathbb{R}^{}\); initial solution \(x_0\in D\); initial temperature \(T_0>0\); proposal distributions \(q_k(\cdot|x)\); number of temperature levels \(K\); inner-loop lengths \(M_k\); cooling rule \(T_{k+1}=C(T_k)\)

Output: Best solution \(x_{\mathrm{best}}\)

\(x\gets x_0\), \(x_{\mathrm{best}}\gets x_0\), \(T\gets T_0\)

for \(k=0,1,\dots,K-1\) do

for \(j=1,2,\dots,M_k\) do

Generate \(y\sim q_k(\cdot|x)\) such that \(y\in D\)

\(\Delta\gets f(y)-f(x)\)

Generate \(u\sim\operatorname{U}(0,1)\)

if \(\Delta\leqslant0\) or \(u\leqslant\exp(-\Delta/T)\) then

\(x\gets y\)

end if

if \(f(x)<f(x_{\mathrm{best}})\) then

\(x_{\mathrm{best}}\gets x\)

end if

end for

\(T\gets C(T)\)

end for

return \(x_{\mathrm{best}}\)

8.3.2 遗传算法

遗传算法(genetic algorithm, GA)是另一类经典的启发式优化方法。相比于模拟退火沿着单个候选解构成的搜索轨迹逐步迭代,遗传算法采用基于种群的搜索机制,在每一代同时保留多个候选解,并通过选择、交叉与变异等遗传操作产生下一代种群。种群中的不同个体可以分布在解空间的不同区域,因此算法能够同时探索多个可能的优化方向,而不必将全部搜索过程集中在某一个当前解附近。

遗传算法的思想来源于生物进化过程:适应度较高的个体具有更大的机会参与繁殖,两个亲本的部分信息通过交叉被组合到子代中,而随机变异则不断向种群引入新的信息。经过多代更新后,优质个体所包含的结构逐渐在种群中积累,从而推动候选解向目标函数较优的区域移动。

遗传算法真正重要的是它所形成的一套种群搜索机制:选择操作利用当前已经发现的优质信息,交叉操作重组不同个体中的已有信息,变异操作则维持种群的多样性并探索新的区域。这三个操作共同决定了算法在探索与利用之间的平衡。

8.3.2.1 个体编码与种群

仍考虑优化问题:

\[\begin{equation*} \min_{x\in D}\;f(x) \end{equation*}\]

其中\(D\)为非空解空间,\(f:D\to\mathbb{R}\)为目标函数。遗传算法将每个候选解\(x\in D\)称为一个个体,并用适合进行交叉与变异操作的数据结构表示该个体,这一表示过程称为编码。

由同一代中的若干个体组成的整体称为一个种群。设种群规模为\(N\in\mathbb{N}^+\),第\(t\)代种群记为:

\[\begin{equation*} \mathcal{D}_t=\left(x_1^{(t)},x_2^{(t)},\dots,x_N^{(t)}\right)\in D^N \end{equation*}\]

初始种群\(\mathcal{D}_0\)通常从\(D\)中随机生成;若已经具有若干质量较高的可行解,也可以将其加入初始种群,但仍应保留足够的随机个体以维持初始多样性。

8.3.2.2 适应度

遗传算法会为每个个体赋予一个用于衡量其优劣程度的数值,称为适应度。设适应度函数为:

\[\begin{equation*} F:D\to\mathbb{R} \end{equation*}\]

其中\(F(x)\)表示个体\(x\in D\)的适应度。适应度越高,通常表示该个体在当前优化问题中越优,并在后续选择过程中具有更大的保留或繁殖机会。

适应度函数通常由目标函数\(f\)构造。例如,对于最小化问题,应使得目标函数值较小的个体具有较高的适应度,即对于\(x_1,x_2\in D\),若:

\[\begin{equation*} f(x_1)<f(x_2) \end{equation*}\]

则通常要求:

\[\begin{equation*} F(x_1)>F(x_2) \end{equation*}\]

适应度的具体构造需要根据具体的问题背景去设计。合理的适应度设计既需要反映不同个体之间的质量差异,也应避免少数个体因适应度过高而过早支配整个种群,从而导致种群多样性迅速下降。

8.3.2.3 选择操作

选择操作从当前种群\(\mathcal{D}_t\)中选取个体作为亲本,用于产生下一代个体。其基本原则是使适应度较高的个体具有更大的被选择概率,从而使较优个体所携带的信息更容易保留下来并传递到后续种群。

设第\(t\)代种群中第\(i\)个个体的适应度为\(F_i^{(t)}\),其被选择的概率为\(p_i^{(t)}\)。选择机制通常应满足:若:

\[\begin{equation*} F_i^{(t)}>F_j^{(t)} \end{equation*}\]

则相应地有:

\[\begin{equation*} p_i^{(t)}>p_j^{(t)} \end{equation*}\]

具体的选择概率如何由适应度确定并无唯一形式,只要能够体现这种”适应度越高,被选择机会越大”的倾向即可。

8.3.2.4 交叉操作

交叉操作将两个亲本的部分编码组合起来生成子代。设两个亲本分别为:

\[\begin{equation*} x=(x_1,x_2,\dots,x_m),\qquad z=(z_1,z_2,\dots,z_m) \end{equation*}\]

对于二进制编码或一般的序列编码,可以随机选择交叉位置\(r\in\{1,2,\dots,m-1\}\),并生成:

\[\begin{align*} y_1&=(x_1,x_2,\dots,x_r,z_{r+1},z_{r+2},\dots,z_m) \\ y_2&=(z_1,z_2,\dots,z_r,x_{r+1},x_{r+2},\dots,x_m) \end{align*}\]

这种操作称为单点交叉。对于实数向量编码,则可以生成算术组合:

\[\begin{align*} y_1&=\lambda x+(1-\lambda)z \\ y_2&=(1-\lambda)x+\lambda z \end{align*}\]

其中\(\lambda\in[0,1]\)可以固定,也可以在每次交叉时随机生成。交叉通常以概率\(p_c\in[0,1]\)执行;若没有执行交叉,则直接复制两个亲本作为子代。

交叉操作只负责重组亲本中已经存在的信息,并不保证生成的子代一定属于可行域\(D\)。例如,单点交叉可能破坏排列编码的合法性,算术交叉在\(D\)非凸时也可能产生不可行解。因此,实际算法必须针对编码方式设计保持可行性的交叉算子,或者在交叉后通过修复、投影或重新生成等方式将子代变为可行解。

8.3.2.5 变异操作

变异操作以较小概率随机改变个体编码中的部分分量。对于二进制编码\(x=(x_1,x_2,\dots,x_m)\in\{0,1\}^m\),可以对每个分量独立生成\(u_j\sim\operatorname{U}(0,1)\),并按照:

\[\begin{equation*} y_j= \begin{cases} 1-x_j,&u_j\leqslant p_m \\ x_j,&u_j>p_m \end{cases} \quad j=1,2,\dots,m \end{equation*}\]

进行变异,其中\(p_m\in[0,1]\)为变异概率。对于实数向量,可以在随机选中的分量上加入随机扰动;对于排列编码,则可以交换两个位置、反转一个区间或将某个元素插入新的位置。

变异能够产生亲本中原本不存在的信息,是防止种群多样性完全消失的重要机制。若\(p_m\)过小,种群可能在较早阶段集中到某个局部最优区域;若\(p_m\)过大,子代与亲本之间的关联又会被大量破坏,使遗传算法逐渐退化为随机搜索。因此,变异概率通常明显小于交叉概率。

8.3.2.6 种群更新与精英保留

在完成选择、交叉和变异后,算法需要用生成的子代构造下一代种群。最简单的更新方式是直接以全部子代替换当前种群,但由于选择与遗传操作都具有随机性,当前种群中的最优个体可能在更新过程中丢失。为避免已经获得的最优目标函数值变差,可以采用精英保留策略:从\(\mathcal{D}_t\)中选出目标函数值最小的\(e\)个个体,直接复制到\(\mathcal{D}_{t+1}\)中,再用新生成的子代补足其余\(N-e\)个位置,其中\(e\in\{0,1,\dots,N-1\}\)。

精英保留保证了历代种群中的最优目标函数值单调不增,但\(e\)也不宜过大。若大量个体被直接复制,种群中可用于产生新组合的位置会减少,并且相同或相近的精英个体可能迅速占据整个种群,从而削弱探索能力。

8.3.2.7 完整算法

Algorithm 8.2 Genetic Algorithm

Input: Feasible set \(D\); objective function \(f:D\to\mathbb{R}\); fitness function \(F\); population size \(N\); number of generations \(G\); selection operator \(\mathcal{S}\); crossover operator \(\mathcal{C}\); mutation operator \(\mathcal{M}\); crossover probability \(p_c\); mutation probability \(p_m\); number of elites \(e\)

Output: Best solution \(x_{\mathrm{best}}\)

Generate an initial population \(\mathcal{D}_0\in D^N\)

\(x_{\mathrm{best}}\gets\underset{x\in\mathcal{D}_0}{\arg\min}\;f(x)\)

for \(t=0,1,\dots,G-1\) do

Evaluate the fitness of all individuals in \(\mathcal{D}_t\)

Let \(\mathcal{E}_t\) contain the \(e\) best individuals in \(\mathcal{D}_t\)

\(\mathcal{Q}_t\gets\varnothing\)

while \(|\mathcal{Q}_t|<N-e\) do

Select two parents \((x_1,x_2)\) from \(\mathcal{D}_t\) according to \(\mathcal{S}\)

if \(\operatorname{U}(0,1)\leqslant p_c\) then

\((y_1,y_2)\gets\mathcal{C}(x_1,x_2)\)

else

\((y_1,y_2)\gets(x_1,x_2)\)

end if

Apply mutation \(\mathcal{M}\) to \(y_1\) and \(y_2\) with probability \(p_m\)

Ensure that the resulting offspring belong to \(D\)

Add offspring to \(\mathcal{Q}_t\)

end while

\(\mathcal{D}_{t+1}\gets(\mathcal{E}_t,\mathcal{Q}_t)\)

\(x^*\gets\underset{x\in\mathcal{D}_{t+1}}{\arg\min}\;f(x)\)

if \(f(x^*)<f(x_{\mathrm{best}})\) then

\(x_{\mathrm{best}}\gets x^*\)

end if

end for

return \(x_{\mathrm{best}}\)


  1. 证明过于复杂,略去。↩︎

COMPUTATIONAL NOTES

计算实验

实验01

Newton 法

Python · Jupyter

优化 Newton 法的局部二次收敛

考虑严格凸函数

f(x, y) = ex − x + e2y − 2y,

其唯一最小值点为 x* = (0, 0)。实验记录 Newton 迭代轨迹、梯度范数和误差平方比。

import numpy as np
import matplotlib.pyplot as plt
from figure_settings import configure_matplotlib
configure_matplotlib()
from scipy.linalg import solve
from scipy.optimize import minimize

def objective(point):
    x, y = point
    return np.exp(x) - x + np.exp(2.0 * y) - 2.0 * y

def gradient(point):
    x, y = point
    return np.array([np.exp(x) - 1.0, 2.0 * np.exp(2.0 * y) - 2.0])

def hessian(point):
    x, y = point
    return np.diag([np.exp(x), 4.0 * np.exp(2.0 * y)])

start = np.array([1.5, 0.75])
verification = minimize(objective, start, jac=gradient, method="BFGS", tol=1e-13)
print("SciPy verification point:", verification.x)
print("Theoretical minimizer: [0. 0.]")
SciPy verification point: [ 3.51528465e-14 -7.89814117e-15]
Theoretical minimizer: [0. 0.]
def newton_history(x0, maximum_iterations=12, tolerance=1e-14):
    point = np.asarray(x0, dtype=float)
    history = [point.copy()]
    for _ in range(maximum_iterations):
        step = solve(hessian(point), gradient(point), assume_a="pos")
        point = point - step
        history.append(point.copy())
        if np.linalg.norm(gradient(point)) <= tolerance:
            break
    return np.asarray(history)

history = newton_history(start)
errors = np.linalg.norm(history, axis=1)
gradient_norms = np.asarray([np.linalg.norm(gradient(point)) for point in history])

print("k        x_k                 y_k             ||x_k-x*||       ||grad f(x_k)||")
for k, (point, error, grad_norm) in enumerate(zip(history, errors, gradient_norms)):
    print(f"{k:>2d}  {point[0]:>18.10e}  {point[1]:>18.10e}  {error:>14.6e}  {grad_norm:>14.6e}")

print("\nSquared-error ratio ||e_(k+1)|| / ||e_k||^2:")
for k in range(len(errors) - 1):
    if errors[k] > 0:
        print(f"k={k}: {errors[k + 1] / errors[k]**2:.6e}")
k        x_k                 y_k             ||x_k-x*||       ||grad f(x_k)||
 0    1.5000000000e+00    7.5000000000e-01    1.677051e+00    7.785293e+00
 1    7.2313016015e-01    3.6156508007e-01    8.084841e-01    2.372186e+00
 2    2.0836118568e-01    1.0418059284e-01    2.329549e-01    5.180029e-01
 3    2.0274917814e-02    1.0137458907e-02    2.266805e-02    4.579881e-02
 4    2.0415408242e-04    1.0207704121e-04    2.282512e-04    4.565490e-04
 5    2.0838026512e-08    1.0419013256e-08    2.329762e-08    4.659524e-08
 6    3.0837457679e-16    1.5418728840e-16    3.447733e-16    4.965068e-16

Squared-error ratio ||e_(k+1)|| / ||e_k||^2:
k=0: 2.874610e-01
k=1: 3.563928e-01
k=2: 4.177058e-01
k=3: 4.442064e-01
k=4: 4.471832e-01
k=5: 6.351999e-01
x_grid = np.linspace(-0.35, 1.65, 320)
y_grid = np.linspace(-0.25, 0.85, 320)
X, Y = np.meshgrid(x_grid, y_grid)
Z = np.exp(X) - X + np.exp(2.0 * Y) - 2.0 * Y

fig, axes = plt.subplots(1, 2, figsize=(12, 4.9), constrained_layout=True)

levels = 2.0 + np.geomspace(1e-4, 8.0, 16)
axes[0].contour(X, Y, Z, levels=levels, colors="#9aa2b1", linewidths=0.8)
axes[0].plot(history[:, 0], history[:, 1], "o-", color="#397d78", linewidth=1.8)
axes[0].scatter([0.0], [0.0], marker="*", s=120, color="#c26b4a", zorder=4)
for k, point in enumerate(history):
    axes[0].annotate(str(k), point, xytext=(5, 5), textcoords="offset points", fontsize=9)
axes[0].set(
    title="Newton Iteration Trajectory",
    xlabel="x",
    ylabel="y",
)

tiny = np.finfo(float).tiny
iterations = np.arange(len(history))
axes[1].semilogy(iterations, np.maximum(errors, tiny), "o-", color="#52649a", label=r"$||x_k-x^*||$")
axes[1].semilogy(iterations, np.maximum(gradient_norms, tiny), "s-", color="#397d78", label=r"$||\nabla f(x_k)||$")
axes[1].set(
    title="Error and Gradient Norm",
    xlabel="Iteration $k$",
    ylabel="Log scale",
)
axes[1].legend()
plt.show()

实验 01「Newton 法」结果 1(图 1)

观察

当迭代点进入最小值邻域后,误差平方比趋于有限常数,误差曲线迅速下降。这对应 Newton 法的局部 Q-二次收敛。

实验02

模拟退火

Python · Jupyter

模拟退火

考虑优化问题:

min(x, y) ∈ Df(x, y),   D = [−5.12, 5.12]2

其中目标函数为二维 Rastrigin 函数:

f(x, y) = 20 + x2 + y2 − 10[cos (2πx) + cos (2πy)]

该函数在 (0, 0) 处取得全局最小值 f(0, 0) = 0,但定义域内还密集分布着许多局部极小值点。先画出目标函数的曲面,观察其非凸结构。

import math
import random

import matplotlib.pyplot as plt
import numpy as np
from simanneal import Annealer

from figure_settings import configure_matplotlib

configure_matplotlib()

LOWER_BOUND, UPPER_BOUND = -5.12, 5.12


def rastrigin(state):
    x, y = state
    return 20.0 + x**2 + y**2 - 10.0 * (
        math.cos(2.0 * math.pi * x) + math.cos(2.0 * math.pi * y)
    )


grid = np.linspace(LOWER_BOUND, UPPER_BOUND, 180)
grid_x, grid_y = np.meshgrid(grid, grid)
grid_z = 20.0 + grid_x**2 + grid_y**2 - 10.0 * (
    np.cos(2.0 * np.pi * grid_x) + np.cos(2.0 * np.pi * grid_y)
)

fig = plt.figure(figsize=(9.0, 6.0))
axis = fig.add_subplot(111, projection="3d")
surface = axis.plot_surface(
    grid_x, grid_y, grid_z, cmap="viridis", rstride=2, cstride=2,
    linewidth=0.0, antialiased=True, rasterized=True
)
axis.set(xlabel="x", ylabel="y", zlabel="f(x, y)",
         title="Rastrigin function: many local minima")
axis.view_init(elev=37, azim=-55)
fig.colorbar(surface, ax=axis, shrink=0.62, pad=0.08, label="f(x, y)")
fig.tight_layout()
plt.show()

实验 02「模拟退火」结果 1(图 1)

使用 simanneal 求解

取初始解 x0 = (4.5, 4.5),并使用以下机制:

  • 候选解生成机制: 每次等概率选择一个坐标,在该坐标上加入标准差为 0.18 的 Gaussian 扰动;若扰动后的坐标越过 D 的边界,则将其反射回可行域。

  • 退火机制: 使用 simanneal 的指数降温。设总迭代次数为 K,则第 k 次迭代的温度为:

    $$T_k=T_{\max}\exp\left[-\frac{k}{K}\log\left(\frac{T_{\max}}{T_{\min}}\right)\right]$$

    其中 Tmax = 30,Tmin = 10−4,K = 60000。对候选解产生的目标函数增量 Δ,当 Δ ≤ 0 时直接接受;当 Δ > 0 时,以概率 exp (−Δ/Tk) 接受。

下面的类只定义目标函数与候选解的生成方式;Metropolis 接受规则和指数降温由 simanneal.Annealer 实现。

INITIAL_STATE = [4.5, 4.5]
STEP_SCALE = 0.18
ANNEALING_STEPS = 60_000


def reflect_to_domain(value):
    # 边界处理:把越过边界的坐标反射回可行域,而不是直接截断。
    while value < LOWER_BOUND or value > UPPER_BOUND:
        if value < LOWER_BOUND:
            value = 2.0 * LOWER_BOUND - value
        if value > UPPER_BOUND:
            value = 2.0 * UPPER_BOUND - value
    return value


class RastriginAnnealer(Annealer):
    copy_strategy = "slice"

    def __init__(self, initial_state):
        self.history = []
        super().__init__(initial_state)

    def move(self):
        # 候选解生成机制:随机选择一个坐标,并加入 N(0, 0.18^2) 扰动。
        coordinate = random.randrange(2)
        proposal = self.state[coordinate] + random.gauss(0.0, STEP_SCALE)
        self.state[coordinate] = reflect_to_domain(proposal)

    def energy(self):
        return rastrigin(self.state)

    def update(self, step, temperature, energy, acceptance, improvement):
        # 每隔固定迭代数记录当前解;不调用默认 update,以免打印进度表。
        self.history.append(
            {
                "step": step,
                "temperature": temperature,
                "energy": energy,
                "best_energy": self.best_energy,
                "x": self.state[0],
                "y": self.state[1],
            }
        )


# 退火机制:simanneal 在 Tmax 与 Tmin 之间按指数规律降温,并使用
# Metropolis 概率 exp(-Delta/T) 决定是否接受目标函数上升的候选解。
random.seed(19)
annealer = RastriginAnnealer(INITIAL_STATE)
annealer.Tmax = 30.0
annealer.Tmin = 1.0e-4
annealer.steps = ANNEALING_STEPS
annealer.updates = 3_000
best_state, best_energy = annealer.anneal()

best_state = np.asarray(best_state)
history = {
    key: np.asarray([record[key] for record in annealer.history])
    for key in annealer.history[0]
}
print(f"Best state: {best_state}")
print(f"Best objective value: {best_energy:.8g}")
Best state: [-1.51453786e-05 -2.79294234e-05]
Best objective value: 2.0026395e-07

解的移动与目标函数值的下降

左图在目标函数等高线上画出当前解的移动。为避免轨迹点过密,只显示每 300 次迭代的一个位置;连线仍由每 20 次迭代记录的状态构成。右图同时显示当前解的目标函数值和迄今为止找到的最优值:高温阶段允许当前值暂时上升,而最优值保持单调不增。

contour_levels = [1, 5, 10, 15, 20, 30, 40, 50, 60, 80]
fig, axes = plt.subplots(1, 2, figsize=(12.0, 4.8))

axes[0].contour(
    grid_x, grid_y, grid_z, levels=contour_levels, colors="0.72", linewidths=0.6
)
axes[0].plot(history["x"], history["y"], color="#2878B5",
             linewidth=0.65, alpha=0.55)
sample = slice(None, None, 15)
path_points = axes[0].scatter(
    history["x"][sample], history["y"][sample],
    c=history["step"][sample], cmap="plasma", s=12, alpha=0.8, zorder=3
)
axes[0].scatter(*INITIAL_STATE, marker="o", s=48, color="#D95319",
                label="start", zorder=4)
axes[0].scatter(*best_state, marker="*", s=125, color="#EDB120",
                edgecolor="black", linewidth=0.5, label="best", zorder=5)
axes[0].set(
    xlim=(LOWER_BOUND, UPPER_BOUND), ylim=(LOWER_BOUND, UPPER_BOUND),
    xlabel="x", ylabel="y", title="Movement of the current solution"
)
axes[0].legend(loc="lower right")
fig.colorbar(path_points, ax=axes[0], label="iteration")

positive_floor = 1.0e-7
axes[1].semilogy(
    history["step"], np.maximum(history["energy"], positive_floor),
    color="0.65", linewidth=0.8, label="current objective value"
)
axes[1].semilogy(
    history["step"], np.maximum(history["best_energy"], positive_floor),
    color="#2878B5", linewidth=1.6, label="best objective value"
)
axes[1].set(
    xlabel="iteration", ylabel="objective value (log scale)",
    title="Decrease of the objective value"
)
axes[1].legend()

fig.tight_layout()
plt.show()

实验 02「模拟退火」结果 1(图 2)

结果说明

搜索初期的温度较高,当前解能够接受一部分目标函数值上升的移动,因此轨迹可以离开初始点附近的局部极小区域,并在多个盆地之间探索。随着温度按指数规律降低,劣化移动越来越难被接受,搜索逐渐集中到已经发现的低目标函数值区域。本次固定随机种子的运行最终到达原点附近,所得最优目标函数值接近理论最小值 0。

实验03

遗传算法

Python · Jupyter

遗传算法求解 Rastrigin 函数

仍考虑二维 Rastrigin 函数的约束优化问题

min(x, y) ∈ Df(x, y),   D = [−5.12, 5.12]2,

f(x, y) = 20 + x2 + y2 − 10{cos (2πx) + cos (2πy)}.

其全局最小值为 f(0, 0) = 0。上一实验已经展示了目标函数的多峰结构,这里直接使用 DEAP 实现遗传算法,并观察种群如何逐代集中到全局最优点附近。

编码与遗传机制

实验采用与正文一致的实数向量编码。每个个体为 z = (x, y) ∈ D,种群规模为 N = 120,进化代数为 G = 200,初始个体在 D上均匀生成。下面依次给出本实验对选择、交叉、变异和种群更新的具体实现。

选择操作。由于 Rastrigin 函数满足 f(z) ≥ 0,定义正适应度

$$F(z)=\frac{1}{1+f(z)}$$

目标函数值越小,适应度越大。设第 t 代种群为 𝒟t = (z1(t), z2(t), …, zN(t)),令

$$p_i^{(t)}=\frac{F(z_i^{(t)})}{\sum_{j=1}^{N}F(z_j^{(t)})},\qquad i=1,2,\dots,N$$

每次按照概率 p1(t), p2(t), …, pN(t) 从当前种群中有放回地抽取一个亲本,重复抽取直到获得产生子代所需的全部亲本。

交叉操作。将依次抽取的两个亲本记为 x, z ∈ D。以概率 pc = 0.9 生成 λ ∼ U (0, 1),并按照正文中的算术交叉生成

y1 = λx + (1 − λ)z,   y2 = (1 − λ)x + λz.

若没有执行交叉,则直接复制两个亲本。由于 D 是凸的盒约束,算术交叉所得子代仍属于 D。

变异操作。每个子代以概率 pm = 0.2 发生变异。发生变异时,等概率选择一个坐标并加入 η ∼ N (0, 0.352);若新坐标越过边界,则将其截断到 [−5.12, 5.12],从而保证变异后的个体仍属于 D。

种群更新。每一代直接复制目标函数值最小的 e = 2 个个体,再用经过选择、交叉与变异产生的 N − e 个子代补足下一代种群。

import math
import random

import matplotlib.pyplot as plt
import numpy as np
from deap import base, creator, tools

from figure_settings import configure_matplotlib

configure_matplotlib()

LOWER_BOUND, UPPER_BOUND = -5.12, 5.12
POPULATION_SIZE = 120
GENERATIONS = 200
CROSSOVER_PROBABILITY = 0.9
MUTATION_PROBABILITY = 0.2
MUTATION_SCALE = 0.35
ELITE_COUNT = 2


def rastrigin(individual):
    x, y = individual
    value = 20.0 + x**2 + y**2 - 10.0 * (
        math.cos(2.0 * math.pi * x) + math.cos(2.0 * math.pi * y)
    )
    return (value,)


def arithmetic_crossover(first, second):
    # 算术交叉:两个子代分别取 lambda 与 1-lambda 的对称组合。
    # 可行域 D 是凸的盒约束,因此交叉后的子代仍然属于 D。
    mixing_weight = random.random()
    first_parent = first[:]
    second_parent = second[:]
    for coordinate in range(2):
        first[coordinate] = (
            mixing_weight * first_parent[coordinate]
            + (1.0 - mixing_weight) * second_parent[coordinate]
        )
        second[coordinate] = (
            (1.0 - mixing_weight) * first_parent[coordinate]
            + mixing_weight * second_parent[coordinate]
        )
    return first, second


def gaussian_mutation(individual):
    # 实数变异:随机选择一个坐标并加入 N(0, 0.35^2) 扰动。
    # 扰动后使用投影(坐标截断)修复越过盒约束边界的个体。
    coordinate = random.randrange(2)
    individual[coordinate] += random.gauss(0.0, MUTATION_SCALE)
    individual[coordinate] = min(
        UPPER_BOUND, max(LOWER_BOUND, individual[coordinate])
    )
    return (individual,)


def fitness_proportionate_selection(individuals, count):
    # 选择操作:F=1/(1+f),再按 p_i=F_i/sum(F_j) 有放回抽取亲本。
    weights = [1.0 / (1.0 + item.fitness.values[0]) for item in individuals]
    return random.choices(individuals, weights=weights, k=count)


# DEAP 的负权重只用于精英和历史最优个体的大小比较。
# 亲本选择使用上面明确定义的正适应度 F=1/(1+f)。
if not hasattr(creator, "RastriginFitnessMin"):
    creator.create("RastriginFitnessMin", base.Fitness, weights=(-1.0,))
if not hasattr(creator, "RastriginIndividual"):
    creator.create(
        "RastriginIndividual", list, fitness=creator.RastriginFitnessMin
    )

toolbox = base.Toolbox()
toolbox.register("coordinate", random.uniform, LOWER_BOUND, UPPER_BOUND)
toolbox.register(
    "individual", tools.initRepeat, creator.RastriginIndividual,
    toolbox.coordinate, n=2
)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)
toolbox.register("evaluate", rastrigin)
toolbox.register("select", fitness_proportionate_selection)
toolbox.register("mate", arithmetic_crossover)
toolbox.register("mutate", gaussian_mutation)
random.seed(2026)
population = toolbox.population(n=POPULATION_SIZE)

# 初始种群需要先计算适应度。后续只重新计算经过交叉或变异的个体。
for individual in population:
    individual.fitness.values = toolbox.evaluate(individual)

hall_of_fame = tools.HallOfFame(1)
hall_of_fame.update(population)
history = []


def record_generation(generation, current_population):
    objective_values = np.asarray(
        [individual.fitness.values[0] for individual in current_population]
    )
    coordinates = np.asarray(current_population, dtype=float)
    history.append(
        {
            "generation": generation,
            "best": objective_values.min(),
            "lower_quartile": np.quantile(objective_values, 0.25),
            "median": np.median(objective_values),
            "upper_quartile": np.quantile(objective_values, 0.75),
            "diversity": np.mean(np.std(coordinates, axis=0)),
        }
    )


record_generation(0, population)
for generation in range(1, GENERATIONS + 1):
    # 精英保留:先复制当前种群中目标函数值最小的 2 个个体。
    elites = [toolbox.clone(individual) for individual in tools.selBest(
        population, ELITE_COUNT
    )]

    # 按适应度比例有放回地选择亲本,再复制得到可以修改的子代。
    offspring = toolbox.select(population, POPULATION_SIZE - ELITE_COUNT)
    offspring = [toolbox.clone(individual) for individual in offspring]

    # 以概率 p_c 对相邻的两个子代执行算术交叉。
    for first, second in zip(offspring[::2], offspring[1::2]):
        if random.random() < CROSSOVER_PROBABILITY:
            toolbox.mate(first, second)
            del first.fitness.values
            del second.fitness.values

    # 以概率 p_m 对每个子代执行 Gaussian 变异与边界修复。
    for individual in offspring:
        if random.random() < MUTATION_PROBABILITY:
            toolbox.mutate(individual)
            del individual.fitness.values

    for individual in offspring:
        if not individual.fitness.valid:
            individual.fitness.values = toolbox.evaluate(individual)

    population = elites + offspring
    hall_of_fame.update(population)
    record_generation(generation, population)

best_individual = np.asarray(hall_of_fame[0], dtype=float)
best_objective = rastrigin(best_individual)[0]
print(f"Best individual: {best_individual}")
print(f"Best objective value: {best_objective:.8g}")
print("Theoretical minimizer: [0. 0.]")
Best individual: [-3.33386638e-05  2.66288495e-05]
Best objective value: 3.6118512e-07
Theoretical minimizer: [0. 0.]

优化过程

左图显示每一代种群中目标函数值的最优值、中位数与四分位区间。纵轴采用对数刻度,但刻度仍表示原始目标函数值,而不是其对数。右图以两个坐标标准差的平均值衡量种群多样性;选择与精英保留推动目标函数值下降,而交叉和变异在搜索过程中维持一定的个体差异。

generation = np.asarray([record["generation"] for record in history])
best = np.asarray([record["best"] for record in history])
lower_quartile = np.asarray([record["lower_quartile"] for record in history])
median = np.asarray([record["median"] for record in history])
upper_quartile = np.asarray([record["upper_quartile"] for record in history])
diversity = np.asarray([record["diversity"] for record in history])
positive_floor = 1.0e-14

fig, axes = plt.subplots(1, 2, figsize=(12.0, 4.6))
axes[0].fill_between(
    generation, np.maximum(lower_quartile, positive_floor),
    np.maximum(upper_quartile, positive_floor), color="#9ECAE1",
    alpha=0.45, label="interquartile range"
)
axes[0].semilogy(
    generation, np.maximum(median, positive_floor), color="#D95319",
    linewidth=1.2, label="population median"
)
axes[0].semilogy(
    generation, np.maximum(best, positive_floor), color="#2878B5",
    linewidth=1.7, label="best individual"
)
axes[0].set(
    xlabel="generation", ylabel="original objective value f(x, y)",
    title="Objective values (logarithmic y-axis)"
)
axes[0].legend()

axes[1].plot(generation, diversity, color="#7E2F8E", linewidth=1.5)
axes[1].set(
    xlabel="generation", ylabel="mean coordinate standard deviation",
    title="Population diversity"
)
fig.tight_layout()
plt.show()

实验 03「遗传算法」结果 1(图 1)

结果说明

初始种群分布在整个可行域中。适应度比例选择使目标函数值较小的个体具有更大的繁殖机会,算术交叉重组两个亲本的坐标信息,精英保留则避免已经找到的最优个体在随机更新中丢失;Gaussian 变异继续向种群引入新的坐标值,使算法仍有机会探索其他局部区域。最终最优个体位于原点附近,其目标函数值接近理论最小值 0。