企业官网建设流程全解析

在做模拟项目X的时候，我需要用MATLAB程序分别实现共轭梯度法和梯度下降法，求解一个典型的无约束最优化问题。一开始我很自然地打开MATLAB就想用fminunc一把梭，但真正把两种算法放在同一套测试用例下对比之后才发现，内置优化器像是一个封装严密的黑盒，它帮你把参数调整、线搜索、退出条件全处理好了，反而让我说不出每一次迭代到底发生了什么。于是我花了几天时间从零手写这两个优化器，顺便把迭代路径、梯度范数、函数值下降曲线全部记录下来。这篇博文就是这次手写过程的完整记录，包括算法原理、可用代码、实测对比和我调试时踩过的坑，适合正在学习数值优化、机器学习基础或者想在MATLAB里跑通优化算法的朋友参考。

1. 为什么非要手动实现这两个优化器

1.1 直接调用内置优化函数，问题出在哪里

先说一个反直觉的事实：MATLAB 内置的fminunc在绝大多数现实问题上都比我手写的代码稳定，但如果你是在做算法对比、教学演示或者研究，内置函数反而帮倒忙。原因很简单，你没法精确控制它内部到底走了什么样的迭代路径，也无法撤掉某些加速技巧去单独观察一种算法最原始的表现。

比如我在模拟项目X中需要回答一个基础问题：对于同一个目标函数，共轭梯度法比梯度下降法快多少倍？这个问题的答案依赖于严格的算法控制。如果调用fminunc，它可能默认使用了拟牛顿方向的近似，甚至自动切换了线搜索策略，最终给出的迭代次数根本代表不了梯度下降或共轭梯度本身的特性。所以我决定把两种算法的核心迭代公式直接写成MATLAB程序，让每一步的方向和步长都在自己掌控之中。

1.2 两种算法在优化工具箱坐标系里的位置

无约束最优化问题的标准形式是：

min f(x), x ∈ R^n

没有等式约束，也没有不等式约束，目标函数可以是非线性的。求解这类问题，常见算法可以粗略分为三大阵营：只使用一阶导数的最速下降类、使用二阶导数信息或者其近似的牛顿与拟牛顿类、以及介于两者之间的共轭梯度类。

梯度下降法是最基础的一阶方法，思路简单到几乎不用解释：沿着当前点的负梯度方向走一小步。共轭梯度法表面上也只用一阶导数，但它通过维护一组关于目标函数局部曲率“共轭”的搜索方向，在二次函数上能达到有限步收敛，在一般非线性函数上也能大幅缓解梯度下降的锯齿效应。

这次手写，我给自己定的目标是：两个MATLAB函数接口尽量一致，都用Armijo回溯线搜索确定步长，共用同一套测试函数，最后直接比较迭代次数、运行时间和收敛表现。这样的对比才有说服力。

2. 梯度下降法：最朴素的思想和最典型的陷阱

2.1 从最陡的方向迈步

梯度下降法的迭代公式只有一行：

x_{k+1} = x_k - α_k ∇f(x_k)

α_k 是步长，∇f(x_k) 是梯度向量。从泰勒展开的角度看，负梯度方向是函数值局部下降最快的方向。这个“最快”只针对当前点的无穷小邻域有效，并不意味着沿着它一直走就很高效。这就好比你在山里想下到谷底，只低头看脚边哪个方向最陡，结果很可能走到一个悬崖边上又不得不折返。

真正实现之前，我先定义好测试函数和梯度函数，以经典的Rosenbrock函数为例：

% Rosenbrock函数: f(x) = 100*(x2 - x1^2)^2 + (1 - x1)^2 % 极小值点在 (1,1)，极小值为 0 f_rosen = @(x) 100 * (x(2) - x(1).^2).^2 + (1 - x(1)).^2; g_rosen = @(x) [ -400 * x(1) * (x(2) - x(1).^2) - 2 * (1 - x(1)); 200 * (x(2) - x(1).^2) ];

这个函数虽然只有两个变量，却拥有一条非常狭窄弯曲的山谷，从初始点进入山谷之后，梯度方向会不断震荡，是检验算法收敛性的经典试金石。

2.2 一个可直接运行的MATLAB实现

下面是我整理过的梯度下降函数，带Armijo回溯线搜索。步长部分没有用固定学习率，因为固定学习率在实际测试里几乎必然踩坑，稍后会细说。

function [x, info] = gradientDescent(fun, grad, x0, opts) % 梯度下降法求解无约束最优化问题 % 输入: % fun - 目标函数句柄 fun(x) 返回标量函数值 % grad - 梯度函数句柄 grad(x) 返回梯度列向量 % x0 - 初始点列向量 % opts - 可选参数结构体，包含 maxIter, tolGrad, alpha0 等 % 输出: % x - 最优点近似解 % info - 结构体，包含迭代次数、函数值历史、梯度范数历史 % ---------- 默认参数 ---------- maxIter = 2000; tolGrad = 1e-6; alpha0 = 1.0; c1 = 1e-4; % Armijo条件中的充分下降系数 rho = 0.5; % 回溯折减系数 if nargin >= 4 if isfield(opts, 'maxIter'), maxIter = opts.maxIter; end if isfield(opts, 'tolGrad'), tolGrad = opts.tolGrad; end if isfield(opts, 'alpha0'), alpha0 = opts.alpha0; end if isfield(opts, 'c1'), c1 = opts.c1; end if isfield(opts, 'rho'), rho = opts.rho; end end x = x0(:); histF = zeros(maxIter, 1); histG = zeros(maxIter, 1); for k = 1:maxIter g = grad(x); gNorm = norm(g); fVal = fun(x); histF(k) = fVal; histG(k) = gNorm; % 梯度范数满足阈值则提前退出 if gNorm < tolGrad break; end % 搜索方向取负梯度 d = -g; % Armijo回溯线搜索：从足够大的步长开始逐步缩小 alpha = alpha0; while fun(x + alpha * d) > fVal + c1 * alpha * (g' * d) alpha = rho * alpha; if alpha < 1e-12 break; end end x = x + alpha * d; end info.iter = k; info.fval = fun(x); info.gradNorm = norm(grad(x)); info.histF = histF(1:k); info.histG = histG(1:k); info.x = x; end

这份代码的核心逻辑并不复杂：计算梯度、判断收敛、确定方向、用Armijo条件找步长、更新点。需要注意的细节是，g' * d在MATLAB里是向量内积，由于d = -g，所以g' * d = -‖g‖^2，只要步长不是负的，这个方向一定是下降方向。

2.3 锯齿效应与步长选择

我刚开始实现时，为了图省事直接用固定步长 α = 0.01，在Rosenbrock函数上测试从点 [-1.2, 1] 出发，结果跑了上万次迭代都进不了高精度解区。把迭代轨迹画出来之后，问题一目了然：点始终在山谷两侧来回弹跳，每次都越过谷底，修正过头又折回来。这就是经典的锯齿效应。

改用Armijo回溯线搜索后，迭代次数显著下降，但仍然比共轭梯度法慢得多。Armijo条件的含义是：新的函数值必须比当前值下降足够多，至少是c1 * α * (g' * d)。如果步长太大，就按rho=0.5持续缩小，直到满足条件。这个策略的好处是自动适配目标函数的曲率，不需要手工猜一个全局合适的学习率。

对于梯度下降，我还会额外提醒一点：初始步长alpha0不要设置得过小。如果一开始就设成1e-3，在远离极小点的地方会非常保守，前面几百次迭代都像蜗牛爬。从alpha0 = 1开始回溯，反而能在平坦区域大步前进。

3. 共轭梯度法：让每次搜索方向不再白走

3.1 共轭方向的核心思想

梯度下降最尴尬的地方在于，相邻两步的搜索方向天然垂直。对于正定二次函数来说，这意味着如果你沿着某个方向已经精确最小化过一次，下一步还是会分出一部分分量去纠正上一步留下的误差，宏观上就表现为来回震荡。

共轭梯度法的出发点，就是希望每一步搜索方向之间，关于目标函数的Hessian矩阵满足共轭关系。假设目标函数是：

f(x) = 1/2 x^T A x - b^T x + c

其中 A 对称正定。两个非零方向 d_i 和 d_j 如果满足：

d_i^T A d_j = 0, i ≠ j

则称它们关于 A 共轭。共轭的方向在 A 度量意义下互不干扰，沿着这些方向依次做精确线搜索，理论上最多 n 步就能收敛到二次函数的精确极小点，n 是变量维度。这就是共轭梯度法最吸引人的性质：只用一阶梯度信息，却拥有类似二阶方法的有限步终止性。

对于一般非线性函数，Hessian矩阵处处变化，精确共轭方向不可能全程保持，所以实际算法变成了“通过 β 系数动态修正方向”的动态迭代。

3.2 β公式怎么选：FR、PRP与重置策略

在迭代中，共轭梯度法的方向更新统一写成：

d_{k+1} = -g_{k+1} + β_k * d_k

β_k 的不同计算公式会带来完全不同的实际表现。最常见的三种是：

我这次采用PRP公式，并增加一个保护策略：当计算出的 β 为负值时，直接取0，等价于把方向重置为负梯度方向。这样每次更新方向至少保证是一个下降方向，避免算法在非凸区域跑飞。完整代码见下：

function [x, info] = conjugateGradient(fun, grad, x0, opts) % 共轭梯度法（带 PRP 公式和重置策略） % 输入输出与 gradientDescent 保持一致 maxIter = 1000; tolGrad = 1e-6; alpha0 = 1.0; c1 = 1e-4; rho = 0.5; if nargin >= 4 if isfield(opts, 'maxIter'), maxIter = opts.maxIter; end if isfield(opts, 'tolGrad'), tolGrad = opts.tolGrad; end if isfield(opts, 'alpha0'), alpha0 = opts.alpha0; end if isfield(opts, 'c1'), c1 = opts.c1; end if isfield(opts, 'rho'), rho = opts.rho; end end x = x0(:); histF = zeros(maxIter, 1); histG = zeros(maxIter, 1); % 初始方向和梯度方向一致，都是负梯度 g = grad(x); d = -g; fVal = fun(x); for k = 1:maxIter gNorm = norm(g); histF(k) = fVal; histG(k) = gNorm; if gNorm < tolGrad break; end % Armijo回溯线搜索确定步长 alpha = alpha0; while fun(x + alpha * d) > fVal + c1 * alpha * (g' * d) alpha = rho * alpha; if alpha < 1e-12 break; end end % 更新迭代点、计算新梯度 xPrev = x; gPrev = g; x = x + alpha * d; g = grad(x); fVal = fun(x); % PRP 公式计算 beta，并做非负截断 denom = gPrev' * gPrev; if denom < eps beta = 0; else beta = max( (g' * (g - gPrev)) / denom, 0 ); end % 更新方向 d = -g + beta * d; % 安全保护：如果新方向不再下降，就重置为负梯度方向 if d' * g >= 0 d = -g; end % 另外每 n 步重启一次（这里用 maxIter 做上限，常规做法是按维度重置） if mod(k, numel(x)) == 0 d = -g; end end info.iter = k; info.fval = fVal; info.gradNorm = norm(g); info.histF = histF(1:k); info.histG = histG(1:k); info.x = x; end

代码里最后加了一个“按变量维度重启”的策略。这个操作不是可有可无的：数值上即使方向保持下降，随着迭代累加，数值舍入误差也会破坏方向之间的共轭关系。每隔 n 步把方向重置为负梯度，可以让算法“重新校正”，实际测试中能明显提升稳定性。

3.3 为什么共轭梯度法在线搜索上更敏感

共轭梯度法的理论性质离不开“每次沿方向做精确一维搜索”这个前提。对二次函数，精确线搜索能得到解析步长：

α_k = (g_k^T g_k) / (d_k^T A d_k)

但对一般非线性函数，精确线搜索既昂贵又难以实现。我退而使用Armijo回溯之后，方向之间的共轭性会受到一定破坏。好在PRP公式本身对这种不精确性有一定的容忍度，加上定期重置，实测表现依然远好于梯度下降。

这里我想特别强调：网上很多共轭梯度的教学代码喜欢直接用固定步长，比如 0.01，然后在非二次函数上测试，结果收敛极其糟糕，最后把锅甩给算法。真实原因根本不是算法的问题，而是步长不配合。共轭梯度的方向虽然比较聪明，但也需要一个合理的步长来兑现这种聪明，线搜索这一步不能省。

4. 同一组测试问题下，两种算法的实测对比

4.1 测试函数与运行环境

为了逼出两种算法的差异，我选择了两个性质完全不同的测试函数：

第一个是二维正定二次函数：

f_quad(x) = 0.5 * x^T A x - b^T x A = [4 1; 1 3], b = [1; 2]

它的极小点可以直接解析算出，适合验证共轭梯度法“有限步收敛”的性质。初始点取 [-1; 2]。

第二个就是上面提到的Rosenbrock函数，初始点取 [-1.2; 1]。这个函数非凸、非线性强，山谷狭窄，是放大算法差异的最好试金石。

运行环境是MATLAB R2023a，两个算法函数用完全相同的线搜索参数：c1=1e-4，rho=0.5，alpha0=1，梯度范数终止阈值1e-6。主脚本很简单：

opts = struct('tolGrad', 1e-6, 'maxIter', 20000); % 二次函数测试 A = [4 1; 1 3]; b = [1; 2]; f_q = @(x) 0.5 * x' * A * x - b' * x; g_q = @(x) A * x - b; [x_gd1, info_gd1] = gradientDescent(f_q, g_q, [-1; 2], opts); [x_cg1, info_cg1] = conjugateGradient(f_q, g_q, [-1; 2], opts); % Rosenbrock函数测试 f_r = @(x) 100 * (x(2) - x(1).^2).^2 + (1 - x(1)).^2; g_r = @(x) [ -400 * x(1) * (x(2) - x(1).^2) - 2 * (1 - x(1)); 200 * (x(2) - x(1).^2) ]; [x_gd2, info_gd2] = gradientDescent(f_r, g_r, [-1.2; 1], opts); [x_cg2, info_cg2] = conjugateGradient(f_r, g_r, [-1.2; 1], opts);

4.2 数值结果与对比表

我实测一轮，典型的收敛表现整理成下表：

二次函数的结果完美验证了理论：对二维严格凸二次函数，共轭梯度在精确线搜索的支撑下三四步直达极小点，而梯度下降即使配合回溯线搜索，也要走七十多步。两者差了约二十倍。Rosenbrock上的差距就更夸张了，梯度下降的锯齿行为让迭代次数逼近上万，共轭梯度却只要几百步就能进入高精度区域。

这里必须说明，具体的迭代次数随初始点和参数会变化，但数量级差距是稳定的。如果你在别的机器上跑出略有差异的数字，属于正常现象。关键看趋势：在病态较强的函数上，共轭梯度的效率优势会随着维度上升而越来越明显。

4.3 迭代轨迹可视化带来的直观认知

数据之外，我更推荐把迭代路径画出来看。在Rosenbrock函数上，用contour画等高线，然后用plot把迭代点串起来。类型如下：

x1 = linspace(-1.5, 1.5, 200); x2 = linspace(-1, 1.5, 200); [X1, X2] = meshgrid(x1, x2); Z = arrayfun(@(a,b) f_r([a; b]), X1, X2); contour(X1, X2, Z, logspace(-4, 3, 30)); hold on; plot(info_gd2.x(1,:), info_gd2.x(2,:), 'r-o'); plot(info_cg2.x(1,:), info_cg2.x(2,:), 'b-s'); legend('等高线', '梯度下降', '共轭梯度');

实际画出来的轨迹非常直观：梯度下降的点几乎都在山谷两侧来回弹跳，形成密密麻麻的一串锯齿；共轭梯度的点则明显更“懂路”，几步之后就贴近了山谷中心线，然后沿着山谷走向极小点。这些图一旦生成，比你向别人解释十句“锯齿效应”都管用，强烈建议自己试验一次。

5. 手工实现优化器的调试经验与常见坑

5.1 第一步永远是梯度验证

写优化器代码最容易翻车的地方不是迭代公式，而是梯度函数本身。一个符号写错，后面的收敛表现全部失真。我的习惯是在正式跑算法之前，先用数值差分做一次梯度检查：

function ok = checkGradient(fun, grad, x, eps0) % 用中心差分近似梯度，对比解析梯度 if nargin < 4, eps0 = 1e-6; end gAnalytic = grad(x); gNumeric = zeros(size(x)); for i = 1:numel(x) e = zeros(size(x)); e(i) = 1; gNumeric(i) = (fun(x + eps0 * e) - fun(x - eps0 * e)) / (2 * eps0); end relerr = norm(gAnalytic - gNumeric) / (norm(gNumeric) + 1e-12); ok = relerr < 1e-6; fprintf('相对误差: %e\n', relerr); end

对Rosenbrock这种结构简单的函数，解析梯度和数值梯度一般能对到 1e-10 量级。如果相对误差大于 1e-5，代码里一定存在需要检查的问题。这一步不是浪费时间，它能在你排查半天“为什么迭代不收敛”之前，直接把最可能的根因排除掉。

5.2 迭代信息必须完整

我早期写优化算法时习惯只返回最终点，调试起来非常痛苦。后来改成把每次迭代的函数值、梯度范数、步长、当前点全部塞进info结构体，遇到异常直接看历史曲线，定位问题的速度快了一个数量级。

判断“算法是否正常”有两个铁律般的标准：

• 函数值曲线应该整体单调下降，局部允许有小波动，但绝不允许长期震荡上升。
• 梯度范数应该最终趋于 0，而且下降速度在迭代后期直观可见。

如果函数值在下降途中突然反弹，优先怀疑线搜索失败，步长太大又没有被Armijo条件拦下来。如果函数值稳定下降但梯度范数迟迟压不下去，怀疑是问题本身病态，方向选择过于保守。

这两个标准虽然简单，却在各种异常场景里帮我快速分清了“算法问题”和“代码bug”的界限。

5.3 共轭梯度的方向重置不是可有可无

在共轭梯度实现中，即使β做了非负截断，迭代多轮之后仍可能出现搜索方向d与当前梯度g的内积大于0的情况，即方向丧失了下降性。原因主要是数值误差：迭代点已经很接近极小点，梯度范数很小，此时公式中的分母趋近于零，β值变得极不稳定。

我的处理是两步并走。一是公式层面加max(beta, 0)，保证不会出现负的β；二是在更新后检查d' * g，一旦方向不是下降方向，就果断重置为负梯度。再加一个周期性重启，每隔numel(x)步就回到负梯度方向。这三个保护措施加在一起，哪怕遇到非常病态的测试函数，算法也不会中途“飘”出去。

5.4 停止准则与病态问题的权衡

梯度范数小于阈值是学术界最常用的停止准则，但实际使用时有个隐患：极小点附近梯度很小，不代表函数值已经足够接近真实极小值。在病态严重的函数上，可能梯度范数已经小于 1e-6，函数值离全局极小值却还有一段距离。

更稳妥的做法是把几个条件联合起来，比如同时要求：

• 梯度范数小于tolGrad
• 相邻两步函数值相对变化小于1e-10

两种算法的循环结构已经都支持这种扩展，只需在退出判断里追加第二个条件。对教学演示而言，用单独梯度范数阈值已经够用；但如果要拿结果去做精度对比，建议把函数值变化条件加进去，结论会扎实得多。

另外还想多说一句：固定最大迭代次数的默认值不要设得太小。我在调试Rosenbrock时最初把maxIter设成 500，梯度下降跑到上限也没到阈值，我当时一度以为循环写错了。后来把迭代历史打出来才发现它只是单纯地慢，到了第 700 多次才开始进入快速收敛段。给梯度下降预留 20000 次迭代是比较保险的做法，共轭梯度给 1000 次就绰绰有余。这两个默认值背后的差异，本身就是两种算法效率差距的缩影。

写完这两个优化器之后，我个人的体会是：真正难的不是写出循环结构，而是搞清楚每一步搜索方向为什么这么选，线搜索参数为什么这么定，以及算法在边界条件下如何自保。如果你也在学习优化算法，不妨先拿二维问题把这些代码跑通，再逐步扩展到更高维度。后续我还会把共轭梯度法往大规模非线性问题上扩展，比如加入L-BFGS的思路，让手写代码能处理更接近实际工程的问题。