1
% Run steepest descent (% Solver settings  % 求解器设置
2
max_iterations = 800; % 设置最大迭代次数为800
3
tol_g = 1e-5; % termination condition ||gradient|| <= tol  % 终止条件：梯度范数小于等于容差
4
alpha0 = 1; % initial step length  % 初始步长为1
5
tau = 0.5; % backtracking parameter  % 回溯参数为0.5
6
beta = 0.001; % for Armijo condition  % Armijo条件参数为0.001
7

8
% Useful data to see progress of solver  % 用于观察求解器进展的有用数据
9
n = numel(x); % 获取变量向量x的元素个数
10
xs = zeros(max_iterations+1, n); % iterate  % 初始化存储所有迭代点的矩阵
11
fs = zeros(max_iterations+1,1); % objective value at each iteration  % 初始化存储每次迭代函数值的向量
12
norm_gs = zeros(max_iterations+1,1); % ||gradient|| at each iteration  % 初始化存储每次迭代梯度范数的向量
13

14
% Set initial data  % 设置初始数据
15
xs(1,:) = x; % 将初始点存储到迭代点矩阵的第一行
16
[f, g] = objfun(x); % 计算初始点的函数值和梯度
17
fs(1,:) = f; % 将初始函数值存储到函数值向量的第一个元素
18
norm_gs(1,:) = norm(g); % 将初始梯度范数存储到梯度范数向量的第一个元素% 运行最速下降法（带线搜索）
19
% Lindon Roberts, 2019  % 作者：Lindon Roberts，2019年
20

21
clear, close all % 清除工作空间变量，关闭所有图形窗口
22

23
% Problem and x0  % 问题设置和初始点x0
24
if 0 % 如果条件为0（即false，不执行此分支）
25
    a = 10; % 设置参数a为10
26
    objfun = @(x)Scaled_quadratic_function(x,a); % 创建缩放二次函数句柄
27
    x = [1; a]; % 设置初始点为[1; a]
28
    xmin = [0;0]; % true minimiser  % 设置真实最小值点为[0;0]
29
else % 否则执行此分支
30
    objfun = @(x)Rosenbrock_function(x); % 创建Rosenbrock函数句柄
31
    x = [1.5; 1.5]; % 设置初始点为[1.5; 1.5]
32
    %x = [-1.2; 1]; % 另一个可选初始点[-1.2; 1]
33
    %x = [0.4; 0.2]; % start close to the solution  % 接近解的起始点[0.4; 0.2]
34
    xmin = [1;1]; % true minimiser  % 设置真实最小值点为[1;1]
35
end
36
nhistory = 10; % use last N iterates to check asymptotic rate  % 使用最后N次迭代来检查渐近收敛率
37

38
% True info  % 真实信息
39
[fmin, gmin, Hmin] = objfun(xmin); % true minimum  % 计算真实最小值点处的函数值、梯度和Hessian矩阵
40
kappa = cond(Hmin); % 计算Hessian矩阵的条件数
41

42
% Solver settings  % 求解器设置
43
max_iterations = 800;
44
tol_g = 1e-5; % termination condition ||gradient|| <= tol
45
alpha0 = 1; % initial step length
46
tau = 0.5; % backtracking parameter
47
beta = 0.001; % for Armijo condition
48

49
% Useful data to see progress of solver
50
n = numel(x);
51
xs = zeros(max_iterations+1, n); % iterate
52
fs = zeros(max_iterations+1,1); % objective value at each iteration
53
norm_gs = zeros(max_iterations+1,1); % ||gradient|| at each iteration
54

55
% Set initial data
56
xs(1,:) = x;
57
[f, g] = objfun(x);
58
fs(1,:) = f;
59
norm_gs(1,:) = norm(g);
60

61
k = 1; % 初始化迭代计数器为1
62
fprintf('  k  |  f(xk)       |  ||grad|| \n'); % 打印表头：迭代次数、函数值、梯度范数
63
fprintf('--------------------------------\n'); % 打印分隔线
64
fprintf('  %i  |  %.4e  |  %.4e  \n', k-1, f, norm(g)); % 打印初始状态信息
65
while k <= max_iterations && norm(g) >= tol_g % 当迭代次数未超过最大值且梯度范数大于容差时继续迭代
66
    s = -g; % steepest descent direction  % 设置最速下降方向为负梯度方向
67
    % Backtracking Armijo linesearch  % 回溯Armijo线搜索
68
    alpha = alpha0; % 设置初始步长
69
    xtest = x + alpha*s; % 计算测试点
70
    while objfun(xtest) > f + beta*alpha*(g'*s) % 当不满足Armijo条件时
71
        alpha = tau*alpha; % 缩小步长
72
        xtest = x + alpha*s; % 重新计算测试点
73
    end
74
    x = xtest; % 更新当前点为测试点
75
    [f, g] = objfun(x); % 计算新点的函数值和梯度
76
    if mod(k, 10) == 0 % 如果迭代次数是10的倍数
77
        fprintf('  %i  |  %.4e  |  %.4e  \n', k, f, norm(g)); % 打印当前状态信息
78
    end
79
    % Save info  % 保存信息
80
    xs(k+1,:) = x; % 将当前点存储到迭代点矩阵
81
    fs(k+1,:) = f; % 将当前函数值存储到函数值向量
82
    norm_gs(k+1,:) = norm(g); % 将当前梯度范数存储到梯度范数向量
83
    k = k + 1; % 迭代计数器加1
84
end
85
fprintf('  %i  |  %.4e  |  %.4e   <- finished\n', k-1, f, norm(g)); % 打印最终状态信息并标记结束
86
fprintf('Finished after %g iterations\n', k-1) % 打印总迭代次数
87
xs = xs(1:k, :); % 截取实际使用的迭代点数据
88
fs = fs(1:k, :); % 截取实际使用的函数值数据
89
norm_gs = norm_gs(1:k, :); % 截取实际使用的梯度范数数据
90
xdists = zeros(k,1); % 初始化存储到最优点距离的向量
91
for i=1:k % 对每个迭代点
92
    xdists(i) = norm(xs(i,:) - xmin); % 计算当前点到最优点的距离
93
end
94

95
% Check asymptotic order of convergence  % 检查渐近收敛阶
96
if numel(fs) < nhistory % 如果函数值个数少于历史记录数
97
    asym_fs = fs; % 使用所有函数值
98
else % 否则
99
    asym_fs = fs(end-nhistory:end); % 使用最后nhistory个函数值
100
end
101
fit_fs = polyfit(log(asym_fs(1:end-1)), log(asym_fs(2:end)), 1); % 对对数函数值进行一次多项式拟合
102
fprintf('Objective values converge with order %1.2f\n', fit_fs(1)); % 打印目标函数值的收敛阶
103

104
%=====================================================
105
% Plot iterates, objective decrease, gradient decrease  % 绘制迭代点、目标函数下降、梯度下降
106
%=====================================================
107

108
subplot(2,2,1); % 创建2×2子图的第1个子图
109
npts = 30; % 设置网格点数为30
110
xplt = linspace(min(min(xs)), max(max(xs)), npts); % 在迭代点范围内创建x轴网格
111
yplt = xplt; % y轴网格与x轴相同
112
[X,Y] = meshgrid(xplt,yplt); % 创建二维网格
113
Z = zeros(npts,npts); % 初始化函数值网格
114
for i=1:npts % 对每个网格点i
115
    for j=1:npts % 对每个网格点j
116
        Z(i,j) = log(objfun([X(i,j); Y(i,j)])); % 计算网格点处的对数函数值
117
    end
118
end
119
contour(X,Y,Z,20) % 绘制20条等高线
120
hold on % 保持当前图形
121
axis equal % 设置坐标轴等比例
122
plot(xs(:,1), xs(:,2), 'r.-', 'MarkerSize', 15); % 绘制迭代轨迹，红色点线，标记大小15
123
xlabel('x1'); % 设置x轴标签
124
ylabel('x2'); % 设置y轴标签
125
hold off % 释放图形保持
126

127
subplot(2,2,2); % 创建2×2子图的第2个子图
128
semilogy(fs-fmin, 'b-', 'Linewidth', 2); % 绘制目标函数值与最小值差的半对数图，蓝色实线，线宽2
129
hold on % 保持当前图形
130
xlabel('Iteration'); % 设置x轴标签为"迭代次数"
131
ylabel('Objective value - fmin'); % 设置y轴标签为"目标函数值 - 最小值"
132
rho = ((kappa-1)/(kappa+1))^2; % 计算最速下降法的理论收敛率
133
fprintf('rho_{SD} convergence rate <= %g (from kappa = %g)\n', rho, kappa); % 打印理论收敛率和条件数
134
semilogy(1:numel(fs), (fs(1)-fmin) * rho.^(0:1:numel(fs)-1), 'r--', 'Linewidth', 2); % 绘制理论收敛曲线，红色虚线，线宽2
135
grid on % 显示网格
136
hold off % 释放图形保持
137

138
subplot(2,2,3); % 创建2×2子图的第3个子图
139
semilogy(norm_gs, 'b-', 'Linewidth', 2); % 绘制梯度范数的半对数图，蓝色实线，线宽2
140
xlabel('Iteration'); % 设置x轴标签为"迭代次数"
141
ylabel('Norm of gradient'); % 设置y轴标签为"梯度范数"
142
grid on % 显示网格
143

144
subplot(2,2,4); % 创建2×2子图的第4个子图
145
semilogy(xdists, 'b-', 'Linewidth', 2); % 绘制到最优点距离的半对数图，蓝色实线，线宽2
146
hold on % 保持当前图形
147
xlabel('Iteration'); % 设置x轴标签为"迭代次数"
148
ylabel('||x-x*||'); % 设置y轴标签为"到最优点的距离"
149
if exist('a') % 如果变量a存在
150
    xs_rate = zeros(numel(xdists),1); % 初始化理论距离收敛率向量
151
    xs_rate(1) = xdists(1); % 设置初始距离
152
    for i=2:numel(fs) % 对每个后续迭代
153
        xs_rate(i) = (a-1)/(a+1) * xs_rate(i-1); % 计算理论距离收敛率
154
    end
155
    semilogy(xs_rate, 'r--', 'Linewidth', 2); % 绘制理论距离收敛曲线，红色虚线，线宽2
156
end
157
grid on % 显示网格
158
hold off % 释放图形保持

1
[f, g, H] = Rosenbrock_function(x)

1
function varargout = Rosenbrock_function(x) % 定义Rosenbrock函数，输入x为变量向量，varargout为可变输出参数
2
    % Rosenbrock test function  % Rosenbrock测试函数
3
    % An interesting start point is x = [-1; 0.8]  % 一个有趣的起始点是x = [-1; 0.8]
4
    % Unique local/global minimum is f=0 at [1;1]  % 唯一的局部/全局最小值点是[1;1]处的f=0
5
    % Lindon Roberts, 2019  % 作者：Lindon Roberts，2019年
6

7
    % Call test functions as:  % 调用测试函数的方式：
8
    %   f         = rosenbrock(x);  % function value only  % 仅获取函数值
9
    %   [f, g]    = rosenbrock(x);  % function value and gradient  % 获取函数值和梯度
10
    %   [f, g, H] = rosenbrock(x);  % function value, gradient and Hessian  % 获取函数值、梯度和Hessian矩阵
11

12
    if nargout > 3 % 如果输出参数个数大于3
13
        error('Invalid number of output arguments'); % 抛出错误：输出参数个数无效
14
    end
15

16
    % Function value  % 计算函数值
17
    a = 10; % 设置参数a为10
18
    f = a*(x(2)-x(1).^2).^2 + (x(1)-1).^2; % 计算Rosenbrock函数值：a*(x2-x1^2)^2 + (x1-1)^2
19
    varargout{1} = f; % 将函数值存储为第一个输出参数
20

21
    if nargout > 1 % 如果需要输出梯度
22
        % Gradient  % 计算梯度
23
        g = [2*(x(1)-1)-4*a*x(1).*(x(2)-x(1).^2); 2*a*(-x(1).^2 + x(2))]; % 计算梯度向量，包含对x1和x2的偏导数
24
        varargout{2} = g; % 将梯度存储为第二个输出参数
25
    end
26

27
    if nargout > 2 % 如果需要输出Hessian矩阵
28
        % Hessian  % 计算Hessian矩阵
29
        H = [2 + 12*a*x(1).^2 - 4*a*x(2), -4*a*x(1); -4*a*x(1), 2*a]; % 计算2×2的Hessian矩阵（二阶偏导数矩阵）
30
        varargout{3} = H; % 将Hessian矩阵存储为第三个输出参数
31
    end
32
end

1
[f, g, H] = Scaled_quadratic_function(x, a)

1
function varargout = Scaled_quadratic_function(x,a) % 定义缩放二次函数，输入x为变量向量，a为缩放参数
2
    % Scaled quadratic, parametrised by a > 0  % 缩放二次函数，由参数a > 0控制
3
    % Large a -> poorly scaled  % 较大的a值会导致函数缩放不良
4
    % An interesting start point is x = [1; a]  % 一个有趣的起始点是x = [1; a]
5
    % Unique local/global minimum is f=0 at [0;0]  % 唯一的局部/全局最小值点是[0;0]处的f=0
6
    % Lindon Roberts, 2019  % 作者：Lindon Roberts，2019年
7

8
    % Call test functions as:  % 调用测试函数的方式：
9
    %   a = 10;  % 设置参数a为10
10
    %   objective = @(x)scaled_quadratic(x,a);  % 创建目标函数句柄
11
    %   f         = objective(x);  % function value only  % 仅获取函数值
12
    %   [f, g]    = objective(x);  % function value and gradient  % 获取函数值和梯度
13
    %   [f, g, H] = objective(x);  % function value, gradient and Hessian  % 获取函数值、梯度和Hessian矩阵
14

15
    if a <= 0 % 如果参数a小于等于0
16
        error('Parameter a must be strictly positive'); % 抛出错误：参数a必须严格为正数
17
    end
18

19
    if nargout > 3 % 如果输出参数个数大于3
20
        error('Invalid number of output arguments'); % 抛出错误：输出参数个数无效
21
    end
22

23
    % Function value  % 计算函数值
24
    f = (a*x(1).^2 + x(2).^2) / 2; % 计算缩放二次函数值：(a*x1^2 + x2^2) / 2
25
    varargout{1} = f; % 将函数值存储为第一个输出参数
26

27
    if nargout > 1 % 如果需要输出梯度
28
        % Gradient  % 计算梯度
29
        g = [a*x(1); x(2)]; % 计算梯度向量：[a*x1; x2]
30
        varargout{2} = g; % 将梯度存储为第二个输出参数
31
    end
32

33
    if nargout > 2 % 如果需要输出Hessian矩阵
34
        % Hessian  % 计算Hessian矩阵
35
        H = [a, 0; 0, 1]; % 计算2×2的Hessian矩阵：对角矩阵[a, 0; 0, 1]
36
        varargout{3} = H; % 将Hessian矩阵存储为第三个输出参数
37
    end
38
end