conjugate gradient method
Brief description
\(\qquad\)The previous article introduced the conjugate direction method. It uses parallel search to construct the conjugate direction: starting from two points and searching once in the same direction, connecting the two one-dimensional minimum points, a new direction conjugate to the search direction is obtained. This method only requires function values, but it requires two one-dimensional searches for each new direction constructed; it is necessary to use it to systematically obtain \(n\) A conjugate direction must be replaced with the corresponding direction.
\(\qquad\)The foreshadowing has already been laid at the end of the previous article: the gradient difference of the quadratic function satisfies \(\nabla f(X^g_1\) – \nabla f(X^^2) = H(X^grad f – X^g_1\)\)so the matrix in the conjugate condition \(H\) It can be replaced by the gradient difference between two adjacent points. As long as the gradient at the iteration point is used, the conjugate direction can be derived directly without calculating \(H\)and there is no need to perform two parallel searches in order to construct a direction. this isconjugate gradient method(Conjugate Gradient Method)。
\(\qquad\)The conjugate gradient method is one of the conjugate direction methods, because each of its conjugate vectors is constructed relying on the negative gradient at the iteration point. Its first search direction takes the negative gradient direction of the initial point, which is the steepest descent method; the search direction of each subsequent step is to deflect the negative gradient by an angle, which is equivalent to correcting the negative gradient. Therefore, the conjugate gradient method is essentially an improvement on the steepest descent method, also known asspin gradient method。
The relationship between conjugate direction and gradient
Consider a positive definite quadratic function
\[f(X) = \frac\X^T H X + B^T X + C
\]
at point \(X^^2\)、\(X^\\) The gradients at are respectively
\[\nabla f(X^) = HX^g_1\ + B, \qquad \nabla f(X^^2) = HX^ < ” << eps << “,停止迭代” << endl;
break;
+ B
\]
Subtracting the two equations, we get
\[\nabla f(X^{(k+1)}) – \nabla f(X^{(k)}) = H\left(X^{(k+1)} – X^{(k)}\right)
\]
like \(X^{(k+1)}\) is from \(X^{(k)}\) along \(S^{(k)}\) Obtained by one-dimensional search, then \(X^{(k+1)} – X^{(k)} = \alpha^{(k)}S^{(k)}\)the above formula is written as
\[\nabla f(X^{(k+1)}) – \nabla f(X^{(k)}) = \alpha^{(k)} H S^{(k)}
\]
Use a certain conjugate direction obtained previously \(S^{(j)}\ (j < k)\) Multiply the formula on the left:
\[[S^{(j)}]^T\left[\nabla f(X^{(k+1)}) – \nabla f(X^{(k)})\right] = \alpha^{(k)}[S^{(j)}]^T H S^{(k)}
\]
and \(S^{(j)}\) and \(S^{(k)}\) Conjugate, the right end is zero, so
\[[S^{(j)}]^T\left[\nabla f(X^{(k+1)}) – \nabla f(X^{(k)})\right] = 0
\]
\(\qquad\)This is the relationship between conjugate directions and gradients: the gradient difference between two adjacent points is orthogonal to all used conjugate directions. Conversely, as long as the gradient difference is equal to \(S^{(j)}\) Orthogonal, due to \(\alpha^{(k)} \ne 0\),\(S^{(j)}\) and \(S^{(k)}\) necessarily about \(H\) Conjugate, the two are equivalent. on the other hand,\(X^{(k+1)}\) along \(S^{(k)}\) The one-dimensional minimum point of , so \([S^{(k)}]^T\nabla f(X^{(k+1)}) = 0\). Putting the two together, the conjugation condition can be completely expressed in terms of gradient,\(H\) It doesn’t have to be – that’s why the conjugate gradient method is able to construct directions from gradients alone.
\(\qquad\)Figure 4-10 illustrates this relationship:\(g_k\)、\(g_{k+1}\) They are \(X^k\)、\(X^{k+1}\) The gradient at\(g_{k+1} – g_k\) Perpendicular to the existing conjugate direction \(d^j\),and \(g_{k+1}\) And perpendicular to the direction of this search \(d^k\)。
Symbol description: in the picture \(x^k\)、\(x^{k+1}\) i.e. textual \(X^{(k)}\)、\(X^{(k+1)}\);\(g_k\)、\(g_{k+1}\) Right now \(\nabla f(X^{(k)})\)、\(\nabla f(X^{(k+1)})\);\(d^j\) is the obtained conjugate direction \(S^{(j)}\),\(d^k\) This is the search direction \(S^{(k)}\)。
recursion formula
Since the conjugation condition can be expressed in terms of gradient differences, a correction can be made in the negative gradient direction to construct a new direction that is conjugate to the previous direction. Pick
\[S^{(k+1)} = -\nabla f(X^{(k+1)}) + \beta^{(k)} S^{(k)}
\]
That is, on the negative gradient of the new point, deflect it along the search direction of the previous step by \(\beta^{(k)}\) share of control. coefficient \(\beta^{(k)}\) By the conjugation condition \([S^{(k+1)}]^T H S^{(k)} = 0\) Sure.
Let’s look at the first step. remember \(g_0 = \equation f(X^{(0)})\)、\(g_1 = \equation f(X^{(1)})\),at this time \(S^{(0)} = -g_0\). Depend on
\[g_1 – g_0 = \alpha^{(0)} H S^{(0)}
\]
Bundle \(S^{(1)} = -g_1 + \beta^{(0)}S^{(0)}\) Substituting the conjugation condition, we get
\[[S^{(1)}]^T H S^{(0)} = \frac{1}{\alpha^{(0)}}\left[-g_1 + \beta^{(0)}S^{(0)}\right]^T\left(g_1 – g_0\right) = 0
\]
Expand and pay attention \(S^{(0)} = -g_0\),and \(g_1\) and \(S^{(0)}\) Orthogonal (one-dimensional search), so \(g_1\) and \(g_0\) It is also orthogonal, so
\[-\|g_1\|^2 + \beta^{(0)}\|g_0\|^2 = 0
\]
Right now
\[\beta^{(0)} = \frac{\|g_1\|^2}{\|g_0\|^2} = \frac{\|\nabla f(X^{(1)})\|^2}{\|\nabla f(X^{(0)})\|^2}
\]
Right \(k\) Repeat the same process step by step to get the general recursive formula
\[\beta^{(k)} = \frac{\|\nabla f(X^{(k+1)})\|^2}{\|\nabla f(X^{(k)})\|^2}
\]
\[S^{(k+1)} = -\nabla f(X^{(k+1)}) + \beta^{(k)}S^{(k)}
\]
\(\qquad\)This coefficient only uses the gradient module of two adjacent points and does not involve \(H\)nor does it involve any matrix operations. This is the difference between the conjugate gradient method and the general conjugate direction method: the direction is still about \(H\) Conjugate, but only gradients are used in the construction process.
\(\qquad\)Notice: Some materials write the second item as \(-\beta^{(k)}\nabla f(X^{(k)})\). Because the first step \(S^{(0)} = -\nabla f(X^{(0)})\)the two writing methods are completely consistent in the first step; but what is inherited during the recursion is the search direction of the previous step, so it is generally written as \(S^{(k)}\)。
specific steps
For a given initial point \(X^{(0)}\) and convergence accuracy \(\varepsilon\):
① Order \(k = 0\);
② Take the negative gradient direction of the initial point as the first search direction
\[S^{(0)} = -\nabla f(X^{(0)})
\]
③ Along \(S^{(k)}\) Perform a one-dimensional search to find the optimal step size \(\alpha^{(k)}\),make
\[f(X^{(k)} + \alpha^{(k)}S^{(k)}) = \min_{\alpha} f(X^{(k)} + \alpha S^{(k)})
\]
And order
\[X^{(k+1)} = X^{(k)} + \alpha^{(k)}S^{(k)}
\]
④ If \(\|\nabla f(X^{(k+1)})\| < \varepsilon\)then stop and output \(X^* = X^{(k+1)}\); Otherwise, go to ⑤;
⑤ If \(k = n\)then order \(X^{(0)} = X^{(k+1)}\)go to ② to re-take the negative gradient direction (reset); otherwise go to ⑥;
⑥ Calculation
\[\beta^{(k)} = \frac{\|\nabla f(X^{(k+1)})\|^2}{\|\nabla f(X^{(k)})\|^2}, \qquad
S^{(k+1)} = -\nabla f(X^{(k+1)}) + \beta^{(k)}S^{(k)}
\]
make \(k = k + 1\)turn to ③.
Symbol explanation:
- \(X^{(k)}\) It’s the first \(k\) The iteration point is also the \(k-1\) The end point of the search;
- \(S^{(k)}\) It’s the first \(k\) search direction,\(n\) direction about \(H\) Pairwise conjugation;
- \(\alpha^{(k)}\) along \(S^{(k)}\) The optimal step length of\(\beta^{(k)}\) is the conjugate coefficient to construct the next direction.
\(\qquad\)Step 5 is usually not reached for quadratic functions:\(n\) The conjugate directions are exhausted,\(n\) The quadratic function has reached the minimum point. Its effect is on non-quadratic functions – at this time \(n\) After this search, the conjugation has been destroyed by the non-quadratic terms of the function, and there is no point in continuing to use the old direction. It is better to start from the negative gradient again and accumulate a set of conjugate directions again.
Block diagram
Symbol description: in the picture \(x\)、\(d^k\)、\(g_k\)、\(\alpha_k\)、\(\beta_k\) corresponding to the main text respectively \(X\)、\(S^{(k)}\)、\(\nabla f(X^{(k)})\)、\(\alpha^{(k)}\)、\(\beta^{(k)}\);The criterion is \(\|\nabla f(X^{(k+1)})\| < \varepsilon\);“\(k = n\)” box is the reset in step 5. After reset, return to the \(S^{(0)} = -\nabla f(X^{(0)})\) box.
Example calculation
Still using the standard functions from the previous articles, complete the calculation in two steps.
topic
\[f(X) = x_1^2 + 2x_2^2 – 4x_1 – 2x_1x_2,\quad X^{(0)} = [1,1]^T,\quad \varepsilon = 0.01
\]
Its gradient and Hessian matrix are
\[\nabla f(X) = \begin{bmatrix} 2x_1 – 2x_2 – 4 \\ 4x_2 – 2x_1 \end{bmatrix},\qquad
H = \begin{bmatrix} 2 & -2 \\ -2 & 4 \end{bmatrix},\qquad
B = \begin{bmatrix} -4 \\ 0 \end{bmatrix}
\]
The theoretical optimal solution is \(X^* = [4,2]^T\),\(f(X^*) = -8\)。
answer
1. First search
\[\nabla f(X^{(0)}) = \begin{bmatrix} 2 – 2 – 4 \\ 4 – 2 \end{bmatrix} = \begin{bmatrix} -4 \\ 2 \end{bmatrix},\qquad
S^{(0)} = -\nabla f(X^{(0)}) = \begin{bmatrix} 4 \\ -2 \end{bmatrix}
\]
along \(S^{(0)}\) Search, at this time \(X = [1+4\alpha,\ 1-2\alpha]^T\),
\[\begin{aligned}
\varphi(\alpha) &= (1+4\alpha)^2 + 2(1-2\alpha)^2 – 4(1+4\alpha) – 2(1+4\alpha)(1-2\alpha) \\
&= 40\alpha^2 – 20\alpha – 3
\end{aligned}
\]
make \(\varphi'(\alpha) = 80\alpha – 20 = 0\),have to \(\alpha^{(0)} = 0.25\),so
\[X^{(1)} = \begin{bmatrix} 1 \\ 1 \end{bmatrix} + 0.25\begin{bmatrix} 4 \\ -2 \end{bmatrix} = \begin{bmatrix} 2 \\ 0.5 \end{bmatrix},\qquad f(X^{(1)}) = -5.5
\]
The gradient at this point is
\[\nabla f(X^{(1)}) = \begin{bmatrix} 4 – 1 – 4 \\ 2 – 4 \end{bmatrix} = \begin{bmatrix} -1 \\ -2 \end{bmatrix},\qquad
\|\nabla f(X^{(1)})\| = \sqrt{(-1)^2 + (-2)^2} = \sqrt{5} \approx 2.236 > 0.01
\]
Precision not reached, continue.
2. Construct the conjugate direction
\[\beta^{(0)} = \frac{\|\nabla f(X^{(1)})\|^2}{\|\nabla f(X^{(0)})\|^2} = \frac{(-1)^2 + (-2)^2}{(-4)^2 + 2^2} = \frac{5}{20} = 0.25
\]
\[S^{(1)} = -\nabla f(X^{(1)}) + \beta^{(0)}S^{(0)} = \begin{bmatrix} 1 \\ 2 \end{bmatrix} + 0.25\begin{bmatrix} 4 \\ -2 \end{bmatrix} = \begin{bmatrix} 2 \\ 1.5 \end{bmatrix}
\]
Check conjugation:
\[H S^{(1)} = \begin{bmatrix} 2 & -2 \\ -2 & 4 \end{bmatrix}\begin{bmatrix} 2 \\ 1.5 \end{bmatrix} = \begin{bmatrix} 1 \\ 2 \end{bmatrix},\qquad
[S^{(0)}]^T H S^{(1)} = [4,-2]\begin{bmatrix} 1 \\ 2 \end{bmatrix} = 0
\]
Conjugate indeed.
3. Second search
along \(S^{(1)}\) Search, at this time \(X = [2+2\alpha,\ 0.5+1.5\alpha]^T\),
\[\begin{aligned}
\varphi(\alpha) &= (2+2\alpha)^2 + 2(0.5+1.5\alpha)^2 – 4(2+2\alpha) – 2(2+2\alpha)(0.5+1.5\alpha) \\
&= 2.5\alpha^2 – 5\alpha – 5.5
\end{aligned}
\]
make \(\varphi'(\alpha) = 5\alpha – 5 = 0\),have to \(\alpha^{(1)} = 1\),so
\[X^{(2)} = \begin{bmatrix} 2 \\ 0.5 \end{bmatrix} + 1\cdot\begin{bmatrix} 2 \\ 1.5 \end{bmatrix} = \begin{bmatrix} 4 \\ 2 \end{bmatrix}
\]
at this time
\[\nabla f(X^{(2)}) = \begin{bmatrix} 8 – 4 – 4 \\ 8 – 8 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix},\qquad f(X^{(2)}) = -8
\]
satisfy \(\|\nabla f(X^{(2)})\| = 0 < \varepsilon\),stop,\(X^* = [4,2]^T\)。
Tabulate the two searches:
| order | starting point | \(\nabla f(X^{(k)})\) | \(\beta^{(k-1)}\) | \(S^{(k)}\) | \(\alpha^{(k)}\) | end | \(f\) |
|---|---|---|---|---|---|---|---|
| 1 | \((1,1)\) | \([-4,2]^T\) | — | \([4,-2]^T\) | \(0.25\) | \((2,0.5)\) | \(-5.5\) |
| 2 | \((2,0.5)\) | \([-1,-2]^T\) | \(0.25\) | \([2,1.5]^T\) | \(1\) | \((4,2)\) | \(-8\) |
\(\qquad\)After two searches, the exact minimum point is reached, which is exactly \(n=2\) step. The first step is to take the steepest descent direction. If you continue to follow the negative gradient, sawtooth will appear like the steepest descent method; the conjugate gradient method takes the negative gradient in the second step. \([1,2]^T\) deflected into \([2,1.5]^T\)in one step.
Comparison with the methods of the previous two articles:
| method | Whether to use gradient | Number of one-dimensional searches | result | \(f\) |
|---|---|---|---|---|
| Coordinate rotation method (\(\varepsilon = 0.01\)) | no | \(8 \times 2 = 16\) | \((3.992188, 1.996094)\) | \(-7.999969\) |
| conjugate direction method | no | \(3\) | \((4,2)\) | \(-8\) |
| conjugate gradient method | yes | \(2\) | \((4,2)\) | \(-8\) |
code example
\(\qquad\)The program follows the forward and backward method and the golden section method for one-dimensional search. Compared with the previous article, there is no need for a second initial point and no parallel search; each iteration only does three things: one-dimensional search along the current direction, calculating the gradient of the new point, and using the ratio of the gradient modules to construct the next direction. Direction is full \(n\) After that, press step 5 to reset to the negative gradient direction.
The core recursion part is as follows:
C++ code example (conjugate direction recursion):
// 3. 构造新的共轭方向
if ((k + 1) % n == 0) {
// 已用满 n 个方向,重置为负梯度方向
for (int i = 0; i < n; ++i) s[i] = -g_new[i];
cout << " k = n,重置 S = -grad f" << endl;
} else {
double g2_old = 0.0;
for (double v : g) g2_old += v * v;
double beta = g2 / g2_old;
for (int i = 0; i < n; ++i) s[i] = -g_new[i] + beta * s[i];
cout << " beta = " << beta << ", ";
print_vec("S", s);
cout << endl;
}
The complete program is as follows, including the objective function, gradient, one-dimensional search and output.
C++ code implementation (complete):
#include
#include
#include
#include
#include
#include
#include
using namespace std;
// 进退法求初始区间 (复用)
pair bracket_minimum(function func, double a0, double h) {
const int MAX_ITER = 1e4;
if (h <= 0) throw invalid_argument("Step size h must be positive.");
if (func(a0) <= func(a0 + h)) {
double x0 = a0 + h, x1 = a0;
h *= 2;
double x2 = x1 - h;
int iter = 0;
while (func(x2) < func(x1)) {
if (iter++ > MAX_ITER) throw runtime_error("Max iterations exceeded");
x0 = x1; x1 = x2; h *= 2; x2 = x1 - h;
}
return make_pair(x2, x0);
} else {
double x0 = a0;
double x1 = a0 + h;
h *= 2;
double x2 = x1 + h;
int iter = 0;
while (func(x2) < func(x1)) {
if (iter++ > MAX_ITER) throw runtime_error("Max iterations exceeded");
x0 = x1; x1 = x2; h *= 2; x2 = x1 + h;
}
return make_pair(x0, x2);
}
}
// 黄金分割法 (复用)
double golden_section(function func, double l, double r, double eps) {
double a = l, b = r, aa, bb;
while (b - a > eps) {
aa = b - 0.618 * (b - a);
bb = a + 0.618 * (b - a);
if (func(aa) < func(bb)) b = bb;
else a = aa;
}
return (a + b) / 2;
}
// 从点 x 出发,沿方向 s 作一维搜索
vector line_search(
function&)> func,
const vector& x,
const vector& s,
double eps
) {
auto phi = [&](double alpha) {
vector y(x.size());
for (size_t i = 0; i < x.size(); ++i) y[i] = x[i] + alpha * s[i];
return func(y);
};
pair bracket;
try { bracket = bracket_minimum(phi, 0.0, 0.1); }
catch (...) { bracket = make_pair(0.0, 10.0); }
double alpha_star = golden_section(phi, bracket.first, bracket.second, eps);
vector y(x.size());
for (size_t i = 0; i < x.size(); ++i) y[i] = x[i] + alpha_star * s[i];
return y;
}
// 输出一个向量
void print_vec(const char* name, const vector& v) {
cout << name << " = [";
for (size_t i = 0; i < v.size(); ++i) cout << v[i] << (i == v.size() - 1 ? "" : ", ");
cout << "]^T";
}
// 共轭梯度法:f 为目标函数,grad 为梯度函数,x0 为初始点
vector conjugate_gradient(
function&)> f,
function(const vector&)> grad,
const vector& x0,
double eps
) {
int n = x0.size();
vector x = x0;
vector g = grad(x);
// S^(0) = -grad f(X^(0))
vector s(n);
for (int i = 0; i < n; ++i) s[i] = -g[i];
for (int k = 0; k < 1000; ++k) {
// 1. 沿 S^(k) 作一维搜索,得到 X^(k+1)
x = line_search(f, x, s, eps);
vector g_new = grad(x);
double g2 = 0.0;
for (double v : g_new) g2 += v * v;
double norm_g = sqrt(g2);
cout << "第 " << k + 1 << " 次搜索: ";
print_vec("X", x);
cout << ", ||grad f|| = " << norm_g << endl;
// 2. 收敛判断
if (norm_g < eps) {
cout << " 满足 ||grad f|| < " << eps << ",停止迭代" << endl;
break;
}
// 3. 构造新的共轭方向
if ((k + 1) % n == 0) {
// 已用满 n 个方向,重置为负梯度方向
for (int i = 0; i < n; ++i) s[i] = -g_new[i];
cout << " k = n,重置 S = -grad f" << endl;
} else {
double g2_old = 0.0;
for (double v : g) g2_old += v * v;
double beta = g2 / g2_old;
for (int i = 0; i < n; ++i) s[i] = -g_new[i] + beta * s[i];
cout << " beta = " << beta << ", ";
print_vec("S", s);
cout << endl;
}
g = g_new;
}
return x;
}
int main() {
cout << fixed << setprecision(6);
cout << "共轭梯度法求解二元二次函数极小值" << endl;
cout << "目标函数: f(X) = x1^2 + 2*x2^2 - 4*x1 - 2*x1*x2" << endl;
cout << "理论最优解: X* = [4.0, 2.0]^T, f(X*) = -8.0" << endl << endl;
auto f = [](const vector& x) {
return x[0] * x[0] + 2 * x[1] * x[1] - 4 * x[0] - 2 * x[0] * x[1];
};
auto grad = [](const vector& x) {
return vector{2 * x[0] - 2 * x[1] - 4, 4 * x[1] - 2 * x[0]};
};
vector x0 = {1.0, 1.0};
double eps = 0.01;
vector result = conjugate_gradient(f, grad, x0, eps);
cout << "\n最终结果:" << endl;
print_vec("最优解 X*", result);
cout << endl;
cout << "目标函数值 f(X*) = " << f(result) << endl;
return 0;
}
\(\qquad\)The actual measurement process of the program is as follows: the first search results in \((2.001808, 0.499096)\),\(\|\nabla f\| = 2.240126\)calculated \(\beta^{(0)} = 0.250908\)the direction is \([1.998208, 1.505417]^T\); The second search results in \((3.990933, 1.997670)\),at this time \(\|\nabla f\| = 0.016100\)has not yet reached \(0.01\) The accuracy of the two directions has been used up, so follow step 5 to reset to the negative gradient direction, search again, and get \(X^* = (3.993963, 1.995688)\),\(f = -7.999978\)。
\(\qquad\)Just two steps to get there \((4,2)\)because the one-dimensional search takes the precise minimum point; the program uses the golden section method, and the accuracy \(0.01\)so the second step stops at \((3.990933, 1.997670)\)which differs from the theoretical value by about \(0.009\)the final error comes from here.
Method analysis
-
Just need gradient: Each step performs a one-dimensional search and calculates the gradient, without forming or inverting the Hessian matrix. The amount of storage is just a few \(n\) dimensional vector, same as steepest descent, much less expensive than Newton’s method.
-
\(n\) step convergence:right \(n\) A positive definite quadratic function, as long as the one-dimensional search is exact,\(n\) When all conjugate directions are used up, the minimum point can be reached. This is the embodiment of the nature of the conjugate direction under the gradient structure.
-
It is an improvement on the steepest descent method.: The first step is the same as the steepest descent method. After each step, the negative gradient is deflected by an angle to avoid sawtooth caused by orthogonal adjacent directions, so it is also called the rotating gradient method. The direction is “effective” each time, unlike the coordinate rotation method, which repeatedly searches along the fixed coordinate axis direction.
-
The difference from the conjugate direction method lies in the construction method: Parallel search requires two one-dimensional searches for each direction constructed, and another initial point must be selected; the conjugate gradient method uses gradient recursion, and advances in one direction in one search, at the cost of being able to calculate the gradient.
-
Conjugation depends on the accuracy of the one-dimensional search: One-dimensional search is not precise, and the conjugation of directions is only approximate.\(n\) The step convergence is no longer strict. This is especially true for non-quadratic functions, so the algorithm retains the reset in step 5: use full \(n\) This direction starts from the negative gradient again and re-accumulates the conjugate direction.
Summarize
\(\qquad\)The conjugate gradient method changes the method of “constructing the conjugate direction” in the previous article: instead of using two parallel searches, it uses the gradient difference between two adjacent points. \(H\) Eliminate and make corrections directly on the negative gradient. Orientation is still about \(H\) Conjugate, quadratic functions are still \(n\) It converges step by step, but each step requires only one one-dimensional search and one gradient calculation.
\(\qquad\)Its cost is also clear: it must be able to calculate gradients. Next article introductionPowell Actit goes back to the way of using only function values, using parallel search plus direction substitution to organize conjugate directions, without gradients or \(H\)is a very representative one of the direct search methods.
Reference materials for this study note:
[1] Bai Qingshun, Sun Jingmin, Liang Yingchun. Mechanical Optimization Design 7th Edition[M]. Beijing: Machinery Industry Press, 2024. ISBN: 978-7-111-75103-8 Online link
[2] Wuhan University of Technology “Optimization Technology B” courseware, taught by teacher: Mr. Yan Bin