首页 / 资讯中心 / 文章详情

迭代算法核心原理与实战:从收敛性到工程实现

迭代算法核心原理与实战:从收敛性到工程实现 ★ FEATURED ARTICLE
1. 迭代算法到底在解决什么问题第一次接触“迭代算法”这个词很多人会以为它是一门具体的算法比如快速排序、二分查找那种。其实不是。迭代算法是一大类求解思路的统称核心逻辑就一句话从一个不太准的初始答案出发反复用同一套规则去修正它让每一次修正后的结果都比上一次更接近真实解直到误差小到可以接受为止。这个思路听起来朴素但它撑起了数值计算、机器学习、运筹优化、图像处理、工程仿真等一大片领域。你平时用的地图路径规划、推荐系统里的矩阵分解、训练神经网络时的梯度下降、求解电路方程时的牛顿法背后都是迭代算法在干活。为什么非要“迭代”因为很多问题的精确解根本写不出解析表达式。比如解一个五次以上的多项式方程数学上已经证明没有通用的求根公式再比如求解一个包含几百万个未知数的稀疏线性方程组直接求逆矩阵的计算量是天文数字。这时候迭代法就派上用场了——我不追求一步到位我一步步逼近用可控的计算量换一个足够好的近似解。这篇文章适合谁看如果你是刚接触数值计算的学生或者在工作中需要手写优化逻辑的开发者又或者你只是好奇“为什么机器学习能训练出来”这件事那接下来的内容会帮你把迭代算法的骨架和血肉都理清楚。我会从设计思路讲到具体实现再把我自己踩过的坑和排查经验一并倒出来。2. 迭代算法的整体设计思路拆解2.1 迭代法的基本骨架不动点思想所有迭代算法都可以抽象成同一个数学模型。假设我们要求解方程 ( f(x) 0 )直接解解不出来那就把它改写成等价形式 ( x g(x) )。然后从一个初始猜测值 ( x_0 ) 出发反复计算[ x_{k1} g(x_k) ]当序列 ( x_0, x_1, x_2, \ldots ) 收敛时极限值 ( x^* ) 就满足 ( x^* g(x^*) )也就是原方程的解。这就是不动点迭代的基本框架。举个生活化的例子。你想找一个数的平方根比如求 ( \sqrt{2} )。你可以猜一个值比如1.5然后用公式 ( x_{k1} \frac{1}{2}(x_k \frac{2}{x_k}) ) 去修正。第一次算出来是1.4167第二次是1.41422第三次就非常接近1.41421了。这个公式就是牛顿法求平方根的特例本质上就是一个不动点迭代。为什么这个框架这么重要因为它把“求解”这个动作转化成了“反复执行一个简单计算”的过程。简单计算意味着可以用计算机高效执行而反复执行意味着精度可以通过增加迭代次数来提升。这种“用时间换精度”的策略是迭代算法最根本的设计哲学。2.2 收敛性迭代算法的生命线迭代算法最怕的事情就是不收敛。你反复算了几万次结果不但没逼近答案反而越跑越远或者在一个区间里来回震荡。所以判断一个迭代算法能不能用第一件事就是看它的收敛性。收敛性分析的核心工具是压缩映射原理。简单说如果迭代函数 ( g(x) ) 满足一个条件存在一个常数 ( L 1 )使得对任意 ( x, y ) 都有 ( |g(x) - g(y)| \leq L|x - y| )那么无论从哪个初始值出发迭代都一定收敛到唯一的不动点。这个 ( L ) 叫做利普希茨常数它衡量的是迭代函数“压缩距离”的能力。用大白话解释如果每次迭代都能把当前答案和真实答案之间的距离缩小到原来的 ( L ) 倍( L ) 小于1那距离就会像滚雪球一样越来越小最终趋近于零。反过来如果 ( L \geq 1 )距离不会缩小迭代就可能发散。在实际操作中我们通常不会去严格验证利普希茨条件而是用收敛阶来衡量迭代速度。收敛阶 ( p ) 的定义是当迭代接近收敛时误差 ( e_{k1} ) 与 ( e_k^p ) 成正比。( p1 ) 叫线性收敛( p2 ) 叫二次收敛。牛顿法在单根附近就是二次收敛意思是每次迭代有效数字大约翻倍收敛非常快。而普通的梯度下降通常是线性收敛速度慢一些但每步计算便宜。注意收敛阶描述的是“接近收敛时”的行为。如果初始值离解很远牛顿法也可能发散。所以实际使用中初始值的选择和收敛判据的设计同样关键。2.3 迭代终止条件什么时候该停下来迭代不能无限跑下去必须有一个停止规则。常见的终止条件有三类第一类是残差判据。比如求解 ( f(x)0 )当 ( |f(x_k)| \epsilon ) 时就停止其中 ( \epsilon ) 是预设的容差。这个判据直接衡量当前答案离满足方程还有多远比较直观。第二类是步长判据。当 ( |x_{k1} - x_k| \epsilon ) 时停止意思是两次迭代的结果已经非常接近再算下去变化也不大了。第三类是最大迭代次数。无论收敛与否迭代到一定次数就强制停止防止死循环。这个条件通常作为兜底和前面两个条件配合使用。选择哪个判据取决于具体问题。如果函数值的量级和自变量的量级差异很大残差判据可能过早停止或过晚停止。这时候需要结合问题的物理意义来设定容差。我个人的经验是至少同时使用残差判据和最大迭代次数步长判据作为辅助参考。因为残差判据直接反映方程满足程度而最大迭代次数防止意外情况。容差 ( \epsilon ) 的设定也有讲究。设得太小迭代次数暴增计算时间不可接受设得太大结果精度不够。一般来说如果数据本身只有单精度约7位有效数字容差设到 ( 10^{-6} ) 就够了如果是双精度约16位有效数字可以设到 ( 10^{-12} ) 甚至更小。但也要看问题的条件数——病态问题即使残差很小解的误差也可能很大。3. 核心细节解析与实操要点3.1 初始值的选择策略迭代算法对初始值敏感吗答案是看算法。线性收敛的算法通常对初始值不太敏感只要在收敛域内从哪出发都能慢慢爬到解。但二次收敛的算法如牛顿法对初始值非常敏感初始值选得不好轻则收敛慢重则直接发散。我试过一个经典的例子求解 ( x^3 - 2x - 5 0 )。从 ( x_0 2 ) 出发牛顿法三次迭代就得到高精度解。但从 ( x_0 0 ) 出发第一次迭代就跳到 ( -2.5 )然后跑到 ( -0.5 ) 附近震荡最后收敛到一个完全不同的根。这就是初始值的影响。那怎么选初始值几个实用策略利用物理背景如果问题来自实际工程初始值通常有物理意义。比如求解温度分布初始值可以设为环境温度求解结构变形初始值可以设为零位移。网格搜索在可能的解区间内均匀取几个点分别迭代几步看哪个收敛最快就用哪个作为正式初始值。降阶近似先忽略高阶项或小参数求解简化后的问题把简化解作为原问题的初始值。前一步结果如果是随时间演化的问题用上一时刻的解作为当前时刻的初始值这叫“热启动”通常非常有效。实操心得对于牛顿法这类敏感算法我习惯先跑几步最速下降法梯度下降来“预热”等接近解附近再切换到牛顿法。这样既有全局收敛的保障又有局部快速收敛的优势。3.2 迭代矩阵与谱半径对于线性方程组 ( Ax b )迭代法可以写成 ( x_{k1} M x_k c ) 的形式其中 ( M ) 叫迭代矩阵。收敛的充要条件是 ( M ) 的谱半径 ( \rho(M) 1 )。谱半径定义为矩阵所有特征值模的最大值。这个条件比利普希茨条件更具体因为它直接和矩阵的特征值挂钩。谱半径越小收敛越快。所以设计迭代法时一个核心目标就是构造一个谱半径尽可能小的迭代矩阵。以雅可比迭代和高斯-赛德尔迭代为例。雅可比迭代把矩阵 ( A ) 拆成对角部分 ( D ) 和非对角部分 ( LU )迭代矩阵是 ( -D^{-1}(LU) )。高斯-赛德尔迭代则利用了最新计算出的分量迭代矩阵是 ( -(DL)^{-1}U )。理论可以证明如果 ( A ) 是对称正定矩阵高斯-赛德尔迭代的谱半径是雅可比迭代谱半径的平方所以收敛速度大约快一倍。但高斯-赛德尔也不是万能的。如果矩阵不是对角占优的两种方法都可能发散。这时候就需要更高级的方法比如超松弛迭代SOR通过引入松弛因子 ( \omega ) 来调节收敛速度。( \omega 1 ) 退化为高斯-赛德尔( \omega 1 ) 是超松弛( \omega 1 ) 是低松弛。最优松弛因子 ( \omega_{opt} ) 可以通过雅可比矩阵的谱半径计算出来[ \omega_{opt} \frac{2}{1 \sqrt{1 - \rho^2}} ]这个公式在求解大型稀疏线性方程组时非常有用能把收敛速度提升一个量级。3.3 数值稳定性与误差累积迭代算法在计算机上执行时浮点数的舍入误差会逐步累积。如果算法本身是数值不稳定的舍入误差会被放大最终淹没真实解。一个经典的稳定性问题是用迭代法求解递推关系时如果迭代函数在某处的导数绝对值大于1误差就会被放大。即使理论收敛实际计算也可能因为舍入误差而偏离。怎么判断数值稳定性一个实用的方法是条件数分析。对于线性方程组条件数 ( \kappa(A) |A| \cdot |A^{-1}| ) 衡量了输入扰动对解的影响。条件数越大问题越病态迭代法越难收敛到高精度。如果条件数达到 ( 10^{12} ) 以上双精度浮点数可能只能保证几位有效数字。应对策略包括预处理对原问题做变换降低条件数。比如用对角矩阵缩放或者用不完全分解作为预条件子。高精度计算在关键步骤使用扩展精度或任意精度算术减少舍入误差。残差修正每迭代若干步后计算残差并求解修正方程把累积误差“洗掉”。注意预处理虽然能改善条件数但预处理本身也有计算成本。如果预处理太复杂省下来的迭代次数可能抵不过预处理的开销。所以预条件子的选择要在效果和成本之间权衡。4. 实操过程与核心环节实现4.1 手写一个通用的迭代求解框架下面我用 Python 写一个通用的迭代求解框架以求解非线性方程 ( f(x)0 ) 为例。这个框架包含牛顿法、割线法和不动点迭代三种模式方便对比效果。import math def newton_method(f, df, x0, tol1e-10, max_iter100): 牛顿法求解 f(x)0 f: 目标函数 df: 导函数 x0: 初始猜测 tol: 容差 max_iter: 最大迭代次数 x x0 history [] for k in range(max_iter): fx f(x) dfx df(x) if abs(dfx) 1e-15: print(f第{k}次迭代导数接近零迭代终止) break x_new x - fx / dfx residual abs(f(x_new)) step abs(x_new - x) history.append((k, x_new, residual, step)) if residual tol or step tol: print(f第{k1}次迭代收敛解为 {x_new:.12f}) return x_new, history x x_new print(f达到最大迭代次数 {max_iter}当前解为 {x:.12f}) return x, history def secant_method(f, x0, x1, tol1e-10, max_iter100): 割线法求解 f(x)0不需要导数 x_prev, x_curr x0, x1 history [] for k in range(max_iter): f_prev f(x_prev) f_curr f(x_curr) if abs(f_curr - f_prev) 1e-15: print(f第{k}次迭代函数值差接近零迭代终止) break x_new x_curr - f_curr * (x_curr - x_prev) / (f_curr - f_prev) residual abs(f(x_new)) step abs(x_new - x_curr) history.append((k, x_new, residual, step)) if residual tol or step tol: print(f第{k1}次迭代收敛解为 {x_new:.12f}) return x_new, history x_prev, x_curr x_curr, x_new print(f达到最大迭代次数 {max_iter}当前解为 {x_curr:.12f}) return x_curr, history def fixed_point_iteration(g, x0, tol1e-10, max_iter100): 不动点迭代求解 x g(x) x x0 history [] for k in range(max_iter): x_new g(x) residual abs(x_new - x) history.append((k, x_new, residual)) if residual tol: print(f第{k1}次迭代收敛解为 {x_new:.12f}) return x_new, history x x_new print(f达到最大迭代次数 {max_iter}当前解为 {x:.12f}) return x, history这个框架可以直接拿来测试不同问题。比如求解 ( f(x) x^3 - 2x - 5 0 )导函数是 ( f(x) 3x^2 - 2 )。从 ( x_0 2 ) 出发牛顿法大约4次迭代就能达到 ( 10^{-10} ) 的精度。4.2 参数选择与计算过程实录以牛顿法求解 ( x^3 - 2x - 5 0 ) 为例我记录一下实际迭代过程迭代次数当前 xf(x)步长02.0000000000-1.0000000000-12.10000000000.06100000000.100000000022.09456812110.00018570000.005431878932.09455148170.00000000180.000016639442.09455148150.00000000000.0000000002可以看到第一次迭代后残差从1降到0.06第二次降到0.00018第三次降到十亿分之一级别。这就是二次收敛的威力——有效数字大约每步翻倍。初始值的选择在这里很关键。如果从 ( x_0 0 ) 出发第一次迭代得到 ( x_1 -2.5 )然后 ( x_2 -0.5 )接着在 ( -0.5 ) 附近震荡最终收敛到另一个根。所以我在代码里加了一个判断如果迭代过程中步长突然变得很大就触发“阻尼”机制把步长乘以0.5再试。阻尼牛顿法的核心改动就一行x_new x - alpha * fx / dfx # alpha 初始为1如果残差没下降就减半这个小小的改动能把牛顿法的收敛域扩大很多。代价是接近解时收敛速度略慢于纯牛顿法但换来的是更强的鲁棒性。4.3 收敛判据的代码实现细节在实际代码中收敛判据不能只写一个abs(f(x)) tol。因为如果 ( f(x) ) 的量级本身很小比如 ( f(x) ) 的值在 ( 10^{-15} ) 附近那残差判据可能永远无法满足。这时候需要结合相对残差residual abs(f(x_new)) relative_residual residual / max(abs(f(x0)), 1e-15) if residual tol or relative_residual tol: break另外步长判据也要考虑相对步长step abs(x_new - x) relative_step step / max(abs(x_new), 1e-15) if relative_step tol: break我通常把绝对容差和相对容差都设上取“或”的关系。这样无论解的量级是大是小都能合理终止。还有一个细节迭代历史记录。把每次迭代的 x、残差、步长都存下来不仅方便调试还能在收敛异常时快速定位问题。比如如果残差在下降但步长在增大说明迭代可能在“ overshoot”需要加阻尼。如果残差震荡不下降说明迭代矩阵谱半径可能大于1需要换方法。5. 常见问题与排查技巧实录5.1 迭代不收敛的典型原因迭代算法不收敛原因通常出在以下几个方面。我整理了一个速查表方便对照排查现象可能原因排查方法解决方案残差震荡不下降迭代矩阵谱半径≥1计算迭代矩阵特征值改用收敛的迭代格式如从雅可比换高斯-赛德尔残差先降后升初始值离解太远牛顿法过冲打印每步的 x 和 f(x)加阻尼或先用梯度下降预热迭代次数异常多收敛阶低或容差设太小检查收敛阶检查容差换高阶方法或放宽容差结果精度不够问题病态条件数大计算条件数预处理或提高计算精度迭代直接发散初始值不在收敛域内网格搜索初始值换初始值或换全局收敛方法残差降到一定程度不再下降舍入误差主导检查残差与机器精度关系接受当前精度或换高精度计算这张表是我自己排查问题时总结的基本上覆盖了八成以上的不收敛情况。其中“残差先降后升”是最常见的尤其是用牛顿法解非线性方程时。这时候不要急着换算法先试试阻尼——把步长乘以一个小于1的因子往往就能救回来。5.2 收敛速度慢的优化技巧如果迭代能收敛但速度太慢有几个立竿见影的优化手段第一换收敛阶更高的方法。比如把线性收敛的不动点迭代换成二次收敛的牛顿法。但要注意高阶方法每步计算量更大如果函数和导数计算很贵可能得不偿失。第二使用松弛因子。对于线性方程组迭代超松弛SOR能把收敛速度提升数倍。最优松弛因子可以通过谱半径公式计算。对于非线性问题也可以在迭代步上乘以一个松弛因子 ( \omega )通过线搜索确定最优值。第三预处理。对于病态问题预处理能显著改善条件数。最简单的预处理是对角缩放把方程组两边同时乘以对角矩阵的逆。复杂一点的用不完全 LU 分解作为预条件子。第四热启动。如果是系列问题比如随时间演化用上一步的解作为当前步的初始值通常能减少一半以上的迭代次数。第五并行化。有些迭代算法天然适合并行比如雅可比迭代的每个分量可以同时更新。虽然收敛速度不变但每步计算时间大幅缩短。实操心得我做过一个测试同一个线性方程组雅可比迭代需要2000次高斯-赛德尔需要1000次SOR最优松弛因子只需要200次。但SOR需要预先估计谱半径这个估计本身有成本。所以如果只解一次高斯-赛德尔可能更划算如果要解很多次相似问题SOR的优势就体现出来了。5.3 数值精度问题的处理经验浮点数精度是迭代算法绕不开的坎。双精度浮点数有大约16位有效数字但迭代过程中舍入误差会累积。如果问题条件数很大最终结果可能只有几位有效数字。我处理精度问题的经验是先估计条件数。如果条件数在 ( 10^6 ) 以下双精度通常够用如果在 ( 10^{12} ) 以上就要考虑高精度计算或预处理。监控残差下降曲线。如果残差降到某个值后不再下降说明已经到达精度极限继续迭代没有意义。使用补偿求和。在累加大量浮点数时用 Kahan 求和算法可以减少舍入误差。必要时用任意精度库。Python 的decimal模块或mpmath库可以提供任意精度算术代价是速度慢很多。还有一个容易被忽视的点迭代公式的代数等价变换。有时候两个数学上等价的公式在浮点运算下精度差异很大。比如计算 ( \sqrt{x^21} - 1 )当 ( x ) 很小时直接算会损失精度改成 ( \frac{x^2}{\sqrt{x^21}1} ) 就稳定得多。迭代算法中的每一步计算都要注意这类数值稳定性问题。6. 迭代算法的扩展与变体6.1 从定常迭代到非定常迭代前面讨论的迭代算法迭代规则是固定的叫定常迭代。还有一类非定常迭代每次迭代的规则会变化。最典型的是共轭梯度法它在线性方程组求解中非常流行。共轭梯度法的核心思想是每次迭代沿着一个“共轭方向”搜索而不是像最速下降法那样沿着梯度方向。这样能保证在 ( n ) 步内收敛到 ( n ) 维线性方程组的精确解忽略舍入误差。对于大型稀疏矩阵共轭梯度法通常比高斯-赛德尔快几个数量级。共轭梯度法的迭代公式稍微复杂一些但代码实现并不长def conjugate_gradient(A, b, x0, tol1e-10, max_iter1000): 共轭梯度法求解 Ax bA 对称正定 x x0.copy() r b - A x p r.copy() rs_old r r for k in range(max_iter): Ap A p alpha rs_old / (p Ap) x x alpha * p r r - alpha * Ap rs_new r r if math.sqrt(rs_new) tol: print(f第{k1}次迭代收敛) return x beta rs_new / rs_old p r beta * p rs_old rs_new print(f达到最大迭代次数 {max_iter}) return x这个算法在求解偏微分方程离散化后的大型稀疏系统时非常有用。我试过一个 10000 维的稀疏系统高斯-赛德尔需要上万次迭代共轭梯度法几百次就收敛了。6.2 随机迭代与蒙特卡洛方法还有一类迭代算法引入了随机性比如随机梯度下降SGD。每次迭代不用全部数据而是随机抽取一个或一小批样本计算梯度。这样做的好处是每步计算量小适合大规模机器学习。SGD 的收敛性和批量梯度下降不同。因为梯度有噪声SGD 不会精确收敛到最小值而是在最小值附近震荡。所以学习率通常要随着迭代逐渐减小这叫“学习率衰减”。随机迭代的收敛分析需要用到随机逼近理论。核心结论是如果学习率满足一定条件比如 ( \sum \eta_k \infty ) 且 ( \sum \eta_k^2 \infty )SGD 几乎必然收敛到局部最小值。这个条件在实践中通常通过分段常数学习率或指数衰减来满足。6.3 迭代算法在机器学习中的角色机器学习几乎就是迭代算法的天下。神经网络的训练过程本质上是求解一个非凸优化问题用的就是迭代法。前向传播计算损失反向传播计算梯度然后用梯度下降更新参数循环往复。为什么不用直接法因为神经网络的参数量动辄百万千万直接求解需要计算 Hessian 矩阵的逆计算量和存储量都不可接受。迭代法每次只用一阶梯度信息计算便宜虽然收敛慢但胜在可扩展。深度学习中的迭代算法还有一些特殊技巧动量法在梯度更新中引入动量项加速收敛并减少震荡。自适应学习率Adam、RMSProp 等算法根据梯度历史自动调整每个参数的学习率。批量归一化在每层输入上做归一化改善优化 landscape让迭代更稳定。这些技巧本质上都是在改善迭代算法的收敛性和稳定性。理解了迭代算法的基本原理再看这些技巧就不会觉得神秘了。7. 我个人的实操体会迭代算法这东西理论分析是一回事实际写代码是另一回事。我刚开始用牛顿法的时候觉得二次收敛太美了什么间题都想往上套。结果有一次求解一个化学反应平衡方程初始值选得不好迭代直接跑到负浓度去了物理上完全不可接受。后来加了边界约束和阻尼才稳定下来。我的体会是永远不要相信“理论收敛”就等于“实际能收敛”。理论分析假设精确算术实际计算有舍入误差理论分析假设初始值在收敛域内实际中你根本不知道收敛域在哪。所以写迭代算法一定要加保护措施最大迭代次数、步长限制、残差监控、异常值检测一个都不能少。另一个体会是收敛判据要结合问题背景来设。我见过有人把容差设到 ( 10^{-15} )结果迭代了几万次还没停因为问题条件数太大残差根本降不到那么低。后来把容差改成 ( 10^{-8} )迭代几十次就停了结果精度完全够用。容差不是越小越好够用就行。最后分享一个小技巧把迭代历史画成图。残差随迭代次数的变化曲线能直观反映收敛速度。如果曲线在对数坐标下是一条直线说明线性收敛如果越来越陡说明超线性收敛如果震荡说明有问题。这个图我每次调试迭代算法都会画比看数字快多了。
阅读完成 · 觉得有帮助?
咨询建站