从实验数据进行建模与参数估计

课程目标与建模的本质

这门课程的核心是从实验数据中建立数学模型并估计参数。模型本质上是对真实或人工过程的数学描述，但必须针对特定目标来构建。同一个物理系统，针对不同的目的可能需要完全不同的模型。

建模目标大致分为几类: 分析类(物理、化学、生物现象分析，软件传感器实现间接测量)，诊断与验证类(假设检验)，教学类(通过模拟器理解系统行为)，预测与控制类(短期预测用于自适应控制，长期预测用于规划，过程控制中的调节和最优控制)，以及信号处理类(如卡尔曼滤波)。建模过程必须始终围绕最终用途来设计模型的复杂度和侧重点。

系统的定义与组成

系统是我们从整个宇宙中划分出来并与之交互的一部分。实际操作中就是明确研究对象的边界。系统与外界有四种关键的交互方式。

首先是输出 y，这是能够直接观测和测量的量。但我们真正关心的往往是另一个量 z，它代表系统内部感兴趣的演化过程，z 可能与 y 不同。能测量的不一定是我们最关心的，最关心的不一定能直接测量。

其次是输入 u，这是主动施加到系统上的作用，用于控制或影响系统行为。最后是扰动 b，它可能作用在系统本身上，也可能作用在测量环节，代表噪声或不确定性。这四个量构成了系统建模的基本框架: 通过 u 驱动系统，系统内部状态演化为 z，我们测量到的是带有扰动 b 的输出 y。

电机系统实例

以电机为例。输出 y 是电机轴的转速，可以直接测量。但我们真正关心的量 z 是电机的扭矩，它决定了负载能力却无法直接测量。输入 u 是供给电机的电流和电压，这是我们可以控制的。扰动 b 体现在 u 和 y 的测量不确定性上。

这里会遇到一个实际问题: 有时识别系统的输入输出并不容易。第一个困难是在某些领域(如生物系统)或某些工况(如不稳定的开环系统)中，无法将系统与环境完全隔离，边界模糊导致输入输出定义困难。第二个困难是有时很难区分哪些是原因(输入)哪些是结果(输出)，因果关系可能相互耦合或循环影响。

模型的定义与分类

模型是一套规则(可以是方程组，也可以是程序代码)，它的作用是: 已知或测量了某些量(y 和/或 u)，通过这套规则计算出其他量，这些计算出的量应该能复现真实系统 S 中对应量的行为。所谓复现就是模型输出应该接近真实系统的对应量，因为模型总是简化的。

模型按照输入输出关系可以分为两大类。

并联模型

在这种模型中，输入是 u，模型计算输出 y_m，希望 y_m 接近真实输出 y。如果我们关心的量 z 与可测量的输出 y 不同，并联模型还可以同时计算 z_m，希望 z_m 接近真实的 z。这种模型与真实系统"并联"运行: 系统和模型都接收相同的输入 u，产生各自的输出，我们对比两者的输出来验证模型。

逆模型

也称为串联模型。在这种模型中，输入是输出 y，模型反向计算出 u_m，希望 u_m 接近真实输入 u。这相当于求解"什么输入会导致观测到的输出"，是一个逆向推断问题。典型应用是控制器设计: 知道期望的输出轨迹，反推需要施加的控制输入。

实际应用中，模型可以同时包含串联和并联部分，形成更复杂的结构。

参数化模型与估计问题

参数化模型的核心特征是模型中包含一些未知的量，这些量需要从收集到的数据 y 和 u 中估计出来，称为参数，记为 p。参数化建模过程需要区分两个概念: 模型的结构 M(.) 和给定参数后的具体模型 M(p)。

模型结构 M(.) 指的是方程的形式、变量之间的关系类型、微分方程的阶数等框架性内容，但其中的系数是未知的。例如，决定用二阶微分方程还是一阶差分方程来描述系统，这是结构选择问题。而给定模型 M(p) 是在结构确定后，用具体的参数值 p 填充进去，得到可以用于计算的完整模型。

特征化

特征化是指选择模型结构 M(.) 的过程，这是重要且困难的步骤。困难在于: 结构选得太简单，模型无法捕捉系统的关键动态;结构选得太复杂，参数无法从有限数据中可靠估计，且模型可能过拟合。特征化的难点在于在没有完全了解系统的情况下，如何选择一个既能描述主要行为、又能从数据中识别参数的恰当结构。

参数估计的具体实施流程

参数调整与尺度问题

一旦模型结构 M(.) 确定后，接下来的任务是尽可能好地调整其中的参数。这里"尽可能好"需要一个明确的定义，这就涉及参数值尺度的问题。在实际参数估计中，不同参数可能有完全不同的数量级和物理意义，比如一个参数可能是时间常数(秒级)，另一个可能是增益系数(可能很大或很小)。如果不对参数进行适当的归一化或缩放，优化算法可能在某些参数方向上步长过大而在其他方向上步长过小，导致收敛困难。

为了度量模型的好坏，引入一个标量函数 j(.)，它是参数的函数(可能还依赖于模型结构)。这个函数被称为准则或代价函数。准则函数的作用是将模型性能映射为一个数值，数值越小表示模型越好。这样就把"调整参数使模型尽可能好"转化为"找到使 j(.) 最小的参数值"这个明确的数学优化问题。

当存在多个候选模型结构时，不仅要为每个结构找到最优参数，还要在不同结构之间进行比较。这个比较同样使用准则函数: 计算每个结构在其最优参数下的准则值，选择准则值最小的那个结构。这里隐含一个关键逻辑: 最佳结构和该结构的最佳参数值都依赖于所选择的准则。如果换一个准则函数，可能会得到完全不同的最优结构和参数。

准则函数的定义与含义

一般情况下，假设准则函数 j(.) 需要被最小化。基于这个约定，可以明确定义模型之间的优劣关系。如果有两个模型 M_1(p_1) 和 M_2(p_2)，我们说 M_1(p_1) 优于 M_2(p_2)，当且仅当

这个不等式给出了模型比较的数学判据。

对于并联模型，最自然的准则是基于输出误差。输出误差定义为真实系统输出与模型输出之差，记为

这里 t 是时间变量，p 是参数向量。这个误差的物理意义是: 在时刻 t，给定参数 p 时，模型预测的输出与实际观测到的输出之间的偏差。真实系统 S 和模型 M(p) 都接收相同的输入 u(t) 和扰动 b(t)，系统产生输出 y(t)，模型产生输出 y_m(t， p)，两者相减得到误差 e_s(t， p)。

准则的具体形式与建模目标的关系

仅仅定义误差还不够，我们需要将所有时刻的误差综合成一个标量准则。一个常用的做法是最小化输出误差的平方和。为什么用平方和?因为误差可能为正或负，直接求和会相互抵消，而平方和既避免了正负抵消，又对大误差给予了更大的惩罚权重。准则函数可以表示为

但这只是众多可能准则中的一种。准则的选择应该反映建模目标。不同的应用场景需要不同的准则。例如，如果建模目标是预测，可能更关心未来时刻的误差;如果建模目标是控制，可能更关心稳态误差或峰值误差;如果建模目标是故障诊断，可能更关心模型对异常数据的敏感性。选择准则时必须明确: 我们希望模型在哪些方面表现最好?

无论选择了哪种准则，接下来都必须对这个准则进行优化。准则不是用来评价模型的静态工具，而是驱动参数搜索的动态目标。选定准则后，整个参数估计问题转化为一个优化问题。

优化算法的工作原理

优化算法的任务是根据所有可用信息计算出参数估计值 \hat{p}，使得准则 j(.) 达到最小值。这里"所有可用信息"包括测量数据、模型结构、准则函数形式等。在某些情况下，优化算法不仅要找到最优参数 \hat{p}，还可能同时选择最优模型结构 \hat{M}，因此输出可能是 (\hat{p}， \hat{M}) 的组合。

优化过程是一个闭环结构。真实系统 S 在扰动 b(t) 和输入 u(t) 作用下产生输出 y(t)。模型 M(p) 接收同样的输入，产生模型输出 y_m(t， p) 和 z_m(t， p)。这些输出进入准则计算模块，该模块根据准则定义计算当前准则值 j(p)。然后这个准则值反馈给优化器，优化器根据准则值调整参数 p，形成闭环。这个闭环不断迭代，直到找到使准则最小的参数。

迭代优化的实施细节

在绝大多数实际情况下，优化是以迭代方式执行的，而非一步求解。迭代优化的基本逻辑是: 从参数的第 k 次估计 \hat{p}_k 出发，计算出下一次估计 \hat{p}_{k+1}，要求新的估计使准则值下降，即

通过不断重复这个过程，准则值逐步减小，参数逐渐逼近最优值。这种迭代策略的优势在于每一步只需要做局部改进，而不需要全局搜索整个参数空间。

迭代优化面临两个核心困难。第一个是初始化问题，即如何选择初始参数 \hat{p}_0。初始值的选择可能显著影响优化结果: 如果初始值离最优解太远，优化可能需要很多迭代步;如果初始值选在某些特殊点(如局部极小点附近或鞍点附近)，优化可能陷入局部最优或无法收敛。第二个是终止条件，即何时停止迭代。常见的终止准则包括: 准则值的变化量小于阈值、参数的变化量小于阈值、迭代次数达到上限、或者梯度范数小于阈值。

根据优化算法执行的时机，可以分为两类。在线算法是指在系统运行过程中实时执行参数估计，可能需要在每个采样时刻或短时间窗口内更新参数，要求计算速度快，适用于参数时变的系统或需要自适应控制的场景。离线算法是指在数据收集完成后，容许一定时间延迟来进行参数估计。离线算法可以使用更复杂的优化方法、更多的迭代次数，追求更高的精度，适用于参数固定或更新周期较长的系统。

估计不确定性与模型集合

参数估计完成后，不应该认为 \hat{p} 对应的模型是唯一可接受的模型。参数估计的结果不是一个确定的模型，而是一个可接受模型的集合。为什么会有多个可接受的模型?因为测量数据本身存在不确定性(噪声、量化误差、传感器误差等)，这些不确定性传播到参数估计中，导致存在多组参数值都能较好地拟合数据，它们的准则值可能非常接近，在不确定性范围内无法区分优劣。

这个可接受模型集合的大小和形状取决于数据采集实验的质量，具体包括以下因素。第一是输入信号的形式选择: 输入信号的频率内容、幅值、持续时间直接影响系统的激励充分性，输入设计不当会导致某些参数不可观测或弱可观测，使得这些参数的不确定性很大。第二是传感器和执行器的位置与类型: 测量什么量、在哪里测量、用什么精度的传感器，决定了数据中包含多少关于系统的信息。第三是测量时间: 采样频率和测量时长影响数据的时间分辨率和统计特性。第四是数字数据处理: 滤波、降采样、预处理等操作可能引入额外的信息损失或失真。

这些因素共同构成了"实验"的概念。为了使可接受模型集合尽可能小(即参数估计的不确定性尽可能低)，需要进行实验设计。实验设计的目标是在给定约束(如时间、成本、安全性)下，优化上述因素的配置，以收集到最具信息量的数据。

结果的批判性分析

模型构建的整个过程涉及许多选择，这些选择或多或少都是基于某些假设或直觉，可能缺乏严格的理论依据。例如选择什么样的模型结构、选择什么准则函数、选择什么初始参数、选择什么终止条件等。因此在得到最终模型后，质疑这些选择是很重要的。

具体的质疑方式是通过一系列测试来挑战模型。这些测试的目的是证伪，或者更严格地说，是反驳。证伪的思想来自科学哲学: 我们无法证明一个理论是对的，但可以通过反例证明它是错的。同样，我们无法证明一个模型在所有情况下都是好的，但可以通过设计特定测试来暴露模型的不足。

不幸的是，没有这样的测试能够证明一个模型是好的(那样的测试被称为验证)。换句话说，验证(确认模型是正确的)在严格意义上是不可能的，因为我们永远无法穷尽所有可能的使用场景。我们能做的只是尽可能多地进行证伪测试，如果模型通过了这些测试，我们增加对它的信心，但仍不能说它已经被"验证"了。

建模流程的系统总结

综合前面所有内容，建模过程本质上包含六个核心步骤，它们构成了一个完整的工作流程。

第一步是数据收集。这是整个流程的基础，数据的质量和信息量直接决定了后续步骤的效果。数据收集不是简单的测量记录，而是涉及实验设计的系统工程。

第二步是选择一个或多个模型结构。这一步对应前面讨论的特征化过程，需要根据系统的物理机制、先验知识、以及建模目标来确定方程形式。之所以说"一个或多个"，是因为实践中可能同时考察多个候选结构，通过后续步骤比较它们的性能。

第三步是选择代价函数。代价函数的选择必须与建模目标一致，不同的目标需要不同的准则。

第四步是优化准则以获得参数的最优数值，并可能选择最合适的结构。这一步执行实际的参数搜索，通过迭代算法最小化代价函数。如果有多个候选结构，这一步还包括结构选择。

第五步是评估估计参数的不确定性。这一步对应前面讨论的可接受模型集合的特征化，量化参数估计的置信区间或置信域，告诉我们参数估计有多可靠。

第六步是结果的批判性分析。这一步通过各种测试挑战模型，寻找可能的缺陷或改进空间。

这六个步骤之间可能需要多次循环。建模不是线性的一次性过程，而是迭代改进的循环过程。例如，如果批判性分析发现模型不够好，可能需要返回第二步重新选择结构，或返回第一步收集更多数据，或返回第三步重新定义准则。每一轮循环都是对模型的一次refinement，直到达到满意的性能或资源耗尽。

模型结构的分类与选择

课程大纲

这一讲的核心主题是模型结构的分类与选择。内容涵盖: 知识驱动模型和行为模型的区分，线性与非线性模型的差异，连续时间与离散时间建模方式，模型中扰动的处理，模型复杂度的选择，以及模型的结构性质(可辨识性和可区分性及其关系)。

结构选择的本质含义

结构选择就是特征化过程。假设我们选择了如下微分方程结构来描述模型输出 y_m 随时间的演化:

这里 u 是输入，p_1 和 p_2 是待估计参数。这个结构选择隐含了两层信息。

第一层是关于系统行为特征的假设: 要求 p_1 和 p_2 为正值。这个正值约束基于物理直觉。如果 p_1 > 0，则 -p_1 y_m 项提供负反馈，使系统稳定;如果 p_2 > 0，则输入 u 正向驱动输出。正值约束定义了一类可能的行为模式: 系统不会产生振荡，是稳定的，并且是线性的。结构选择实际上是在划定"系统可能表现出哪些行为"的边界。

第二层信息是关于参数空间的先验约束: 参数向量 p 必须属于一个先验可允许参数集合。这个集合的定义来自于我们对系统的先验知识或物理合理性的要求。例如，如果从物理意义上某个参数必须为正(如质量、电阻、反应速率常数)，则参数空间就被限制在正实数范围内。这种先验约束减少了参数估计的搜索空间，提高了估计效率，同时也保证了估计结果的物理可解释性。

符号约定与基本概念

在接下来的讨论中，需要明确几个核心符号的含义。M(.) 表示模型结构，它是一个框架或模板，其中包含待定参数。M(p) 表示给定参数向量 p 后的具体模型，此时所有系数都已确定，模型可以用于计算和预测。\mathbb{P} 表示先验可允许参数向量的集合，它定义了参数的可行域。

选择 M(.) 和 \mathbb{P} 的过程具有双重性质。一方面，这很困难，因为它需要领域专业知识、物理直觉和经验。对于同一个系统，不同的专家可能提出不同的结构，没有绝对的正确答案，只能通过后验检验来评估。另一方面，这是根本性的，因为结构选择限制了模型可能产生的输出类型。一旦结构确定，模型的表达能力就被框定了。如果结构本身无法描述系统的某些关键行为(例如用线性结构去描述本质上非线性的现象)，那么无论如何优化参数都无法得到好的模型。

模型类型的范围界定

为了简化讨论，本课程只考虑具有局部化参数的模型，即参数的数量是有限的，并且可以用常数向量表示。这些模型可以用三种数学形式描述。第一种是代数方程，直接描述输入输出之间的静态关系或约束关系，不涉及时间导数或差分。第二种是常微分方程(ODE)，描述系统状态随时间连续变化的动态过程。第三种是有限差分方程，描述离散时间系统的递推关系。

这里介绍的方法原则上可以推广到偏微分方程(PDEs)。偏微分方程描述时间和空间中的分布参数系统，例如热传导、流体动力学等。但偏微分方程的参数估计涉及更复杂的数学工具(如泛函分析、逆问题理论)，并且计算成本高昂，因此这里不展开讨论。

结构选择中的多维权衡

模型结构的选择需要综合考虑多个维度的权衡。

第一个维度是所追求的目标。如果目标是实时控制，模型必须计算快速;如果目标是长期预测，模型必须具有良好的泛化能力;如果目标是物理机制理解，模型必须具有明确的物理意义。

第二个维度是模型将被使用的条件。这包括工作范围、输入信号的性质、以及与控制系统其他元素的通信方式等。例如，如果模型需要在宽广的工作范围内有效，就不能采用局部线性化;如果输入信号可能包含高频成分，模型必须具有足够的时间分辨率;如果模型需要嵌入到分布式控制系统中，就必须考虑通信延迟和数据格式兼容性。

第三个维度是构建模型的成本，包括时间成本和资源成本。高保真度的复杂模型可能需要数月的建模工作和大量的计算资源，这在实际工程中可能不可行。

第四个维度是研究的可能性。即使理论上某个复杂模型是最优的，但如果我们无法获得足够详细的数据来估计其中的大量参数，那么这个模型在实践中就是不可用的。因此结构选择必须与数据的可获得性相匹配。开发一个参数很多的非常复杂的模型是没有意义的，如果无法获得详细数据来估计这些参数。

最后有一个重要的假设: 大多数模型假设如果初始系统处于静止状态并施加一个恒为零的输入 u，则输出 y 在没有任何扰动的情况下也将恒为零。这意味着系统没有自发的内部动力学，其演化完全由外部输入驱动。这一假设简化了分析，但要求 u 和 y 以适当的参考坐标系表示，即它们应该表示相对于某个平衡点的偏差，而不是绝对值。

知识驱动模型的构建原理

知识驱动模型是基于物理学基本原理构建的模型。其构建过程依赖于两类核心方程。第一类是守恒方程，包括质量守恒、动量守恒、能量守恒等。这些守恒定律是物理世界的基本约束，确保模型在宏观层面上满足热力学和力学的基本法则。第二类是平衡方程，描述系统内部不同组分或不同力之间的平衡关系，例如化学平衡、热平衡、力的平衡等。

知识驱动模型适合利用先验可用信息。如果我们对系统的物理机制有深入理解，知道哪些守恒律适用、哪些物理过程主导系统行为，就可以从第一性原理出发推导模型方程，而不需要完全依赖数据。同时，它适合对估计参数的数量级进行后验控制。知识驱动模型的参数通常具有明确的物理意义(如反应速率常数、传热系数、摩擦系数等)，因此我们可以根据物理直觉或文献值判断估计出的参数是否合理。如果某个参数的估计值与物理预期严重偏离，这可能提示模型结构有问题或数据质量不佳。

化学反应系统的建模实例

考虑一个化学反应系统。反应机理为: 物质 A 与物质 B 之间以速率常数 p_1 正向反应和速率常数 p_2 逆向反应达到可逆平衡，同时 B 以速率常数 p_3 不可逆地转化为 C。反应图式为:

建模时做出两个简化假设。第一，所有反应动力学都是一阶的，即反应速率与反应物浓度成正比。第二，反应器是等温的且搅拌完全，即温度和浓度在整个反应器中均匀分布。基于这些假设和质量平衡原理，可以写出三个物质浓度的微分方程:

这三个方程的推导逻辑如下。对于 A，正向反应消耗 A(速率为 p_1[A])，逆向反应生成 A(速率为 p_2[B])，因此 A 的净变化率为两者之差。对于 B，正向反应生成 B(速率为 p_1[A])，逆向反应消耗 B(速率为 p_2[B])，向 C 的反应也消耗 B(速率为 p_3[B])，因此 B 的净变化率为生成项减去两个消耗项。对于 C，只有从 B 的转化生成 C(速率为 p_3[B])，没有消耗，因此 C 单调增加。

这个例子展示了知识驱动模型的核心特征: 模型结构由对系统的先验知识强加。我们不是从数据中"学习"方程形式，而是根据反应机理直接写出方程。参数和状态变量具有具体的物理或化学意义: 参数 p_i 是反应速率常数，状态变量 [A]， [B]， [C] 是反应物浓度。这种具体意义使得模型具有可解释性和可验证性。

知识驱动模型的复杂性挑战

当处理复杂的实际过程时，知识驱动模型可能面临一系列挑战。

第一，模型可能具有大维度。例如，如果从偏微分方程出发对空间分布系统建模，需要对空间进行离散化，导致状态变量数量大幅增加。一个三维流体动力学问题可能涉及数百万个网格点，每个点都有多个状态变量。

第二，模型可能由许多方程组成。多组分化学反应、多物理场耦合系统往往需要数十甚至上百个耦合方程来描述。第三，模型往往是非线性的。真实物理系统很少是完全线性的，非线性来源于化学反应的质量作用定律、材料的非线性本构关系、流体的对流项等。

第四，这些模型有时难以仿真。非线性、高维、刚性(stiff)方程的数值求解可能面临收敛困难、计算时间长、数值不稳定等问题。第五，这些模型很少能直接用于设计控制律。控制器设计通常需要模型的解析形式或者快速可计算的输入输出关系，而复杂的知识驱动模型往往无法满足实时性要求。

但知识驱动模型有一个重要优势: 它非常适合进行精确的长期仿真。因为模型基于物理守恒律，不会随时间累积系统性偏差，能够可靠地预测系统在长时间尺度上的演化。

行为模型的本质特征

行为模型采用完全不同的建模哲学。行为模型的目标很直接: 简单地复现观测到的行为。这种模型有三个根本性的特征。

第一个特征是: 你甚至不需要知道输出是什么。这听起来有些反直觉，但实际意思是你不需要理解输出的物理意义。例如，在黑箱建模中，你可能只知道某个传感器读数与系统状态相关，但不清楚这个读数的物理解释(它可能是多个物理量的复杂混合)。行为模型不关心这些，只要能拟合数据就行。

第二个特征是: 模型结构不对应于过程的物理结构。知识驱动模型的方程结构反映物理机制，而行为模型的结构是数学上方便拟合的形式，与物理过程的内部结构无关。例如，一个多项式模型或神经网络模型的"层"或"项"没有物理对应物。

第三个特征是: 模型参数因此与所研究的过程没有清晰的关联。行为模型的参数是纯数学系数，没有物理单位或物理意义。你无法说"这个参数代表热容"或"那个参数代表阻尼系数"。参数只是调节拟合曲线形状的旋钮。

多项式模型的万能逼近性

多项式模型是行为模型的典型代表。一个关于时间 t 的多项式模型可以写为:

这个模型能够以任意精度复现任何实验数据集。这是基于魏尔斯特拉斯逼近定理: 任何在闭区间上的连续函数都可以用多项式任意精度地逼近，只要多项式的阶数足够高。

这个性质看似强大，但实际上暴露了行为模型的根本问题: 拟合能力不等于预测能力。虽然高阶多项式可以完美拟合训练数据，但它可能在训练区间之外表现极差(外推性差)，并且高阶项导致过拟合，捕捉噪声而非真实规律。因此，行为模型需要通过正则化、交叉验证等技术来控制复杂度。

除了简单的多项式，还有更复杂的行为模型工具，包括深度神经网络、样条函数、支持向量机。这些方法在机器学习中广泛应用，它们的共同特点是具有强大的拟合能力但缺乏可解释性。

知识驱动与行为模型的系统对比

两类模型有明显的差异。

参数的意义

知识驱动模型的参数具有具体的物理或化学含义，可以与文献值或理论计算对比验证。行为模型的参数没有具体含义，只是拟合系数。

仿真难度

知识驱动模型的仿真往往困难，因为涉及非线性微分方程组的数值求解，可能存在刚性、多尺度等数值挑战。行为模型的仿真容易，通常只是简单的代数计算或函数评估，计算成本低。

先验信息的利用

知识驱动模型可以纳入先验信息，例如守恒律、对称性、参数的物理范围等。行为模型忽略先验信息，完全依赖数据驱动，不考虑物理约束。

有效性域

知识驱动模型的有效性域大，前提是模型正确。因为基于物理原理，模型在理论上适用于所有满足这些原理的工况，包括未在训练数据中出现的工况。行为模型的有效性域受限，只能可靠地预测训练数据覆盖的区域附近的行为，外推能力弱。

这个对比揭示了一个核心权衡: 知识驱动模型以建模成本和计算成本换取泛化能力和可解释性，行为模型以有效性域和可解释性换取建模简便性和计算效率。

线性模型的双重含义

线性性有两种不同的定义，它们对应不同的数学和物理意义，必须加以区分。

关于输入的线性性(LI)

这种线性性意味着叠加原理适用于输入 u。数学表达为: 对于任意标量 \lambda， \mu \in \mathbb{R}^2 和任意时间 t \in \mathbb{R}^+，有

这个性质的物理意义是: 系统对两个输入的线性组合的响应等于对各个输入响应的相同线性组合。这是控制理论意义上的线性性。这种线性性还隐含假设了平稳性，即系统对时间平移不变。如果在时刻 t_0 施加输入产生某个输出，则在时刻 t_0 + \tau 施加相同输入将产生时间平移 \tau 后的相同输出。这对应时不变系统的概念。

关于参数的线性性(LP)

这种线性性意味着叠加原理适用于参数向量 p。数学表达为: 对于任意标量 \lambda， \mu \in \mathbb{R}^2 和任意时间 t \in \mathbb{R}^+，有

这个性质是统计意义上的线性性。它的重要性在于参数估计的计算复杂度。如果模型关于参数线性，参数估计问题通常可以转化为线性最小二乘问题，存在解析解或高效的数值算法。如果模型关于参数非线性，参数估计变成非凸优化问题，需要迭代算法，可能陷入局部最优，初始值选择变得关键。

需要强调的是，这两种线性性是相互独立的。一个模型可以关于输入线性但关于参数非线性(例如 y = p_1 \sin(p_2 u))，也可以关于输入非线性但关于参数线性(例如 y = p_1 u^2 + p_2 u^3)，或者两者都线性(线性系统的线性参数化)，或者两者都非线性(一般非线性系统)。

线性性判断的实际应用

在实际建模中，必须判断一个给定的模型结构 M(.) 是否满足关于参数的线性性(LP)或关于输入的线性性(LI)。这个判断对估计算法的选择有重要影响。

例子1: y_m(t， p) = pu(t)

这个模型关于参数 p 是线性的，因为输出与 p 成正比关系;同时关于输入 u 也是线性的，因为输出与 u 成正比。因此这是一个既是LP又是LI的模型。

例子2: y_m(t， p) = p_1 y_m(t-1， p) + p_2 u(t)

这个模型是一个自回归模型。检查关于参数的线性性: 由于 y_m(t-1， p) 本身依赖于 p，形成嵌套依赖，因此不是LP。检查关于输入的线性性: 由于存在自回归项 y_m(t-1， p)，而这一项在时间上累积了历史输入的影响，叠加原理仍然成立，因此是LI。

例子3: y_m(t， p) = pu^2(t)

关于参数 p 仍然是线性的，因为 p 以线性方式出现。但关于输入 u 是非线性的，因为输出与 u 的平方成正比，不满足叠加原理。这是LP但非LI的典型例子。

例子4: y_m(t， p) = p_1 y_m^2(t-1， p) + p_2 u(t)

这里自回归项变成了平方项 y_m^2(t-1， p)。这导致关于参数非线性(因为 y_m^2 对 p 的依赖是非线性的)且关于输入也非线性(因为平方项破坏了叠加原理)。因此这既不是LP也不是LI。

例子5: y_m(t， p) = p_1 y(t) + p_2 u(t)

这里的关键是输出表达式中包含了真实系统的输出 y(t)(而不是模型输出 y_m)。关于参数 p 是线性的，这很明显。但关于输入 u 的线性性无法确定。原因是: 真实输出 y(t) 与输入 u 之间的关系由真实系统决定，而我们不知道真实系统是否满足叠加原理。如果真实系统是线性的，则这个模型关于输入线性;如果真实系统是非线性的，则这个模型关于输入也非线性。

线性模型的优先性及其权衡

在可能的情况下，LP和LI模型是优先选择的。

对于LI模型，其优势在于存在强大的数学理论支撑。线性系统理论提供了关于稳定性分析、最优控制设计、扰动影响分析等方面的完整工具箱。但LI模型的局限性也很明显: 有效性域受限。所有真实系统在大输入下都会表现出非线性行为，线性模型只能在输入幅值较小的工作范围内有效。

对于LP模型，其关键优势是关于参数的线性性极大地简化了参数估计的三个方面。第一，参数估计本身变得简单，可以用线性最小二乘等解析方法求解，不需要迭代优化。第二，确定这些参数的不确定性变得容易，因为线性模型的参数估计统计性质有明确的理论结果(如参数估计的方差协方差矩阵有解析表达式)。第三，实验设计也变得更直接，因为可以用线性代数工具分析哪些输入配置能最好地激励系统以估计参数。

但LP模型有一个根本缺陷: LP模型的参数往往没有物理意义。这是因为为了获得参数线性性，往往需要选择数学上方便但物理上不自然的参数化方式。例如后面将看到的变量变换方法，虽然实现了线性化，但新参数与原始物理参数的关系可能很复杂，失去了直接的物理解释。

通过变量变换实现线性化

有时可以通过变量变换将非线性参数模型(non-LP)转化为线性参数模型(LP)。

考虑一个非LP模型:

这个模型关于参数 p_1 是线性的，但关于参数 p_2 是非线性的(因为 p_2 出现在指数函数的指数位置)。直接估计这个模型需要非线性优化。

但如果对输出取对数，得到:

其中定义了新参数 q_1 = \ln p_1 和 q_2 = p_2。现在对数输出 \ln y_m 关于新参数 q = [q_1， q_2]^T 是线性的。这样就可以用线性回归方法估计 q，然后通过反变换 p_1 = \exp(q_1) 和 p_2 = q_2 得到原始参数。

但这种变换对估计值有影响。在存在噪声的情况下，直接估计 p 与估计 q 后反变换得到的 p 的统计性质不同。原因是变换改变了噪声的分布特性。如果原始模型中噪声是加性高斯的(y = y_m + \epsilon，其中 \epsilon 服从高斯分布)，那么对数变换后噪声变成乘性的(\ln y = \ln y_m + \ln(1 + \epsilon/y_m) \approx \ln y_m + \epsilon/y_m)，不再是简单的加性高斯噪声，这会影响估计的统计效率和最优性。

尽管如此，这种线性化技巧仍然很有价值。它可以用来初始化非LP模型的迭代搜索。具体做法是: 先用变换后的LP模型快速得到一个参数估计，然后用这个估计作为原始非LP模型迭代优化的初始值。这样可以避免初始值选择不当导致的收敛问题，提高非线性优化的成功率。

连续时间模型的基本形式

自然系统通常在连续时间中演化，因此用常微分方程或偏微分方程描述是很自然的。连续时间模型的标准形式是状态空间表示，由两个方程组成:

第一个方程是状态方程，描述状态向量 \mathbf{x}(t) 随时间的演化规律。状态是系统的内部变量，包含了系统在当前时刻的全部"记忆"。状态导数由函数 \mathbf{f} 确定，该函数依赖于当前状态 \mathbf{x}、参数 \mathbf{p}、输入 \mathbf{u} 和时间 t。初始条件 \mathbf{x}(0) = \mathbf{x}_0(\mathbf{p}) 可能依赖于参数，因为初始状态本身可能是待估计的量。

第二个方程是输出方程，描述可观测输出 \mathbf{y}_m(t) 与状态、输入和参数的关系。输出方程建立了内部状态与外部观测之间的联系。在很多情况下，状态本身不能直接测量，只能通过输出间接推断。

平稳线性连续时间模型

对于平稳(时不变)线性系统，状态空间方程具有特殊的矩阵形式:

这里 \mathbf{A}(\mathbf{p}) 是系统矩阵(描述状态自演化)，\mathbf{B}(\mathbf{p}) 是输入矩阵(描述输入如何驱动状态)，\mathbf{C}(\mathbf{p}) 是输出矩阵(描述状态如何映射到输出)，\mathbf{D}(\mathbf{p}) 是直通矩阵(描述输入对输出的直接影响)。

如果初始状态为零(\mathbf{x}(0) = 0)，可以考虑传递函数表示:

其中 s 是拉普拉斯变量，\mathbf{H}(s， \mathbf{p}) 是传递函数矩阵，定义为:

这个公式的推导来自对状态方程进行拉普拉斯变换。传递函数表示的优势是直接给出输入输出关系，不涉及内部状态，便于频域分析。

在时域中可以得到连接输入到输出的常微分方程组，形式为:

这是输入输出微分方程表示，\mathbf{P}_i(\mathbf{p}) 和 \mathbf{Q}_j(\mathbf{p}) 是参数矩阵。这种表示的物理意义是: 输出及其各阶导数的线性组合等于输入及其各阶导数的线性组合。这对应于传递函数的有理分式形式。

在不同表示之间切换是容易的，但对于非LI模型这很困难。线性系统的各种表示(状态空间、传递函数、微分方程)之间存在成熟的转换算法，可以用线性代数工具实现。但对于非线性系统，这些表示之间的等价性不复存在，转换变得复杂甚至不可能。

离散时间模型的动机与形式

离散时间模型的数值仿真具有两个显著优势: 更简单和更快。连续时间微分方程的数值求解需要选择积分器(如Runge-Kutta方法)、控制时间步长、处理刚性等问题，计算复杂且耗时。而离散时间模型只涉及代数递推，每一步计算是直接的函数评估，无需迭代求解微分方程。

这就是为什么离散时间模型被广泛用于实时过程控制。在实时控制中，控制器必须在每个采样周期内完成计算并输出控制信号，计算速度至关重要。离散时间模型的快速计算特性使其成为模型预测控制等先进控制策略的基础。

但离散时间建模有一个根本代价: 丧失了关于底层连续时间模型行为的信息。采样过程本质上是信息损失的过程。如果系统在采样间隔内存在快速振荡或瞬态响应，这些信息在采样点之间会丢失。此外，采样可能引入混叠现象，高频动态被误认为低频动态。

离散时间状态空间模型

类似于连续时间情况，离散时间模型也可以用状态方程表示:

这里时间 t 是离散的整数索引，t = 0， 1， 2， \ldots。状态方程变成差分方程: 下一时刻的状态 \mathbf{x}(t+1) 由当前时刻的状态、输入和参数通过函数 \mathbf{f} 递推得到。

对于线性情况，离散时间状态空间模型为:

这与连续时间线性模型形式上非常相似，只是微分变成了差分。

如果 \mathbf{x}(0) = 0，可以考虑传递函数表示:

其中 z 是Z变换变量(离散时间的频域变量)，传递函数为:

这与连续时间的拉普拉斯传递函数形式平行，只是 s 换成了 z。

离散时间差分方程与延迟算子

在时域中，离散时间线性系统可以表示为差分方程:

这个方程说明当前输出是过去输出和过去输入的线性组合。为了简化这种涉及时间延迟项的表达式，引入单位延迟算子 q^{-1}，定义为:

算子 q^{-1} 作用在信号上，将信号向后平移一个时间步。这样，q^{-2}x(t) = x(t-2)，依此类推。

使用延迟算子，差分方程可以重写为算子多项式形式。例如，方程

可以改写为:

其中 A(q^{-1}) = 1 + a_1 q^{-1} + a_2 q^{-2} 和 B(q^{-1}) = b_1 q^{-1} + b_2 q^{-2}。这种算子表示法使得差分方程的代数操作变得类似于多项式运算，简化了理论分析和控制器设计。

离散时间模型的固有问题

离散时间模型存在三个本质问题。第一，这些模型对测量时间点施加了约束。一旦选择了采样周期，所有测量和控制动作都必须在这些离散时刻进行，失去了连续时间的灵活性。如果需要改变采样率，模型参数也会改变，因为离散时间模型的参数依赖于采样周期。

第二，离散时间模型可能隐藏相关连续系统的振荡特性。这是采样的本质限制。如果系统在采样间隔内发生振荡，但振荡频率高于奈奎斯特频率(采样频率的一半)，这些振荡在采样序列中不可见。例如，如果连续系统以某个频率振荡，但每次采样恰好在振荡的同一相位点，采样序列看起来是常数，完全掩盖了振荡。

第三，离散时间模型的参数没有清晰的物理意义，并且其数值依赖于离散化步长。这是因为离散化过程引入了与采样周期相关的人为参数。例如，连续时间系统的时间常数在离散化后变成与采样周期相关的复杂组合，失去了原有的物理解释。改变采样周期会改变离散模型的参数值，即使底层物理系统没有变化。

还有一个重要的事实: 离散和连续模型的性质并不总是类似的。例如离散时间振荡的一阶线性模型:

这个模型在 0 < p_1 < 1 时会产生振荡行为(因为负系数导致符号交替)。但对于连续时间的一阶线性模型，振荡是不可能的。一阶连续系统的特征方程只有一个实根，不可能产生复数极点，因此不能产生振荡。这说明离散化可能引入连续系统中不存在的行为模式，使得通过离散模型反推连续系统变得复杂。

连续与离散时间模型的系统对比

连续时间和离散时间模型有几个关键差异。

参数

连续时间模型的参数独立于采样时间，因为模型本身不涉及采样。离散时间模型的参数依赖于采样时间，改变采样周期需要重新计算参数。

仿真难度

连续时间模型的仿真复杂，需要数值积分方法求解微分方程。离散时间模型的仿真容易，只需要简单的递推计算。

先验信息的利用

连续时间模型可以利用先验信息，因为物理知识通常以连续时间微分方程的形式表达(如牛顿定律、能量守恒等)，这些先验知识可以直接融入连续时间模型。离散时间模型不能利用先验信息，因为物理定律不是以差分方程形式给出的，将连续物理定律转化为离散形式会丢失信息和直观性。

测量时间点的要求

连续时间模型对测量时间点没有限制，原则上可以在任意时刻进行测量，模型可以预测任意时刻的输出。离散时间模型强制要求特定的测量时间点，只能在采样时刻进行测量和预测。

扰动建模的基本框架

扰动通常用随机过程的实现来描述，即随机变量序列。在这一部分，考虑具有简化结构的系统: 标量输入 u(.)、标量输出 y(.)、标量扰动 b(.)，并且用差分方程描述。

这种简化设定的目的是聚焦于扰动建模的核心概念，而不被高维系统的复杂性分散注意力。标量情况下的原理可以直接推广到向量情况。选择差分方程形式是因为离散时间随机过程的理论更为成熟和直观，概率分布、期望、方差等概念在离散时间中更容易定义和计算。

扰动建模的核心挑战是: 扰动是不可测的随机过程，但它影响系统的输入输出行为。我们需要在模型中以某种方式表征扰动的统计特性(如均值、方差、自相关函数等)，使得模型预测能够反映这种随机性。

ARX结构：最简单的扰动建模

在差分方程中引入扰动的最简单方法是直接在递归方程中添加噪声项。考虑标准的自回归方程：

这个方程的结构可以分解为三部分理解。等号左边是输出的自回归部分（auto-regressive），包含当前输出和过去 n_a 个时刻的输出。等号右边第一项是外生输入部分（exogenous），包含输入信号及其延迟项，其中 n_r 是输入输出之间的纯延迟（时滞）。等号右边第二项 b(t) 是噪声项。

噪声项 b(t) 在时刻 t = 1, 2, \ldots 是独立同分布随机变量 \varepsilon(t) 的实现。独立同分布（iid, independent and identically distributed）意味着每个时刻的噪声相互独立，且都服从相同的概率分布（通常假设为零均值、方差为 \sigma^2 的正态分布）。

PPT中有一个重要的括号注记：对于离散化的连续时间系统，输入输出延迟 n_r 至少为1。这是因为在因果系统中，当前时刻的输出不能依赖于当前时刻的输入（需要时间来传播），因此在离散时间表示中，输入至少要延迟一个采样周期才能影响输出。

一旦确定了阶数 n_a、n_b 和延迟 n_r，参数向量为：

通常选择 n_a = n_b，这简化了结构选择过程。这种结构被称为ARX结构，其中AR代表自回归（Auto-Regressive），X代表外生变量（eXogeneous variable）。ARX是最基础的线性动态系统建模框架。

ARMAX结构：引入移动平均噪声

ARX结构的一个限制是噪声灵活性不足。ARX假设递归方程中的噪声是独立同分布的白噪声，但实际系统中的扰动往往具有时间相关性。为了解决这个问题，可以用独立同分布噪声 \varepsilon(t) 的线性组合来替代 b(t)，这种线性组合称为移动平均（moving average）：

这里噪声部分变成了 c_1 \varepsilon(t-1) + \cdots + c_{n_c} \varepsilon(t-n_c)，它是过去 n_c 个独立噪声的加权和。移动平均的物理意义是：当前观测到的扰动不仅受当前白噪声影响，还受过去若干时刻白噪声的累积影响。这种结构能够描述具有短期记忆的随机过程。

确定阶数 n_a、n_b、n_c 和延迟 n_r 后，参数向量扩展为：

这种结构称为ARMAX，其中MA代表移动平均（Moving Average）。ARMAX在ARX基础上增加了噪声建模的灵活性，能够更好地拟合具有相关扰动的系统。

ARARX结构：自回归噪声模型

除了移动平均，还有另一种增强噪声灵活性的方法：将ARX结构中的白噪声 b(t) 替换为自回归部分 \eta(t)。自回归噪声 \eta(t) 本身由另一个差分方程定义：

这个方程说明当前噪声 \eta(t) 依赖于过去的噪声值和当前的白噪声激励 \varepsilon(t)。将这个自回归噪声代入主方程，得到：

确定阶数后，参数向量为：

这种结构称为ARARX。相比ARMAX用移动平均描述噪声，ARARX用自回归描述噪声。自回归噪声的特点是具有长期记忆：过去的噪声会通过递归关系持续影响未来，衰减速度由参数 d_i 控制。这适合建模具有慢变趋势或低频扰动的系统。

ARARMAX结构：最通用的扰动模型

最通用的结构是将自回归噪声和移动平均噪声结合，得到ARARMAX结构。噪声部分 \eta(t) 由ARMA模型描述：

代入主方程：

参数向量包含所有系数：

ARARMAX结构最灵活，能够描述各种复杂的噪声特性。但灵活性的代价是参数数量增加，估计难度和计算成本上升。

统一的算子多项式表示

为了统一表示这些结构，引入延迟算子多项式。定义四个算子多项式：

这里 A(q, p) 对应输出的自回归部分，B(q, p) 对应输入的外生部分（注意 q^{1-n_r} 项引入了时滞），C(q, p) 对应移动平均噪声部分，D(q, p) 对应自回归噪声部分。

使用这些算子多项式，ARARMAX结构可以紧凑地写为：

这里 p^* 表示参数的真实值。这个公式揭示了系统的完整传递函数结构：输入到输出的传递函数为 B(q, p)/A(q, p)，噪声到输出的传递函数为 C(q, p)/[A(q, p)D(q, p)]。

PPT还给出了框图表示，清晰展示了信号流。白噪声 \varepsilon(t) 通过 C(q, p)/D(q, p) 滤波器产生相关噪声，输入 u(t) 通过 B(q, p) 滤波器产生确定性输出，两者相加后通过 1/A(q, p) 动态系统产生最终输出 y(t)。

结构之间的包含关系

ARARMAX结构包含了ARX、ARARX和ARMAX作为特例。ARX对应 C(q, p) = 1 且 D(q, p) = 1，此时噪声传递函数退化为 1/A(q, p)，噪声直接加在递归方程中。ARMAX对应 D(q, p) = 1，噪声传递函数为 C(q, p)/A(q, p)，具有移动平均特性。ARARX对应 C(q, p) = 1，噪声传递函数为 1/[A(q, p)D(q, p)]，具有自回归特性。

PPT的框图展示了这些特例。ARX最简单：噪声 \varepsilon(t) 和输入滤波输出直接相加，经过 1/A(q, p) 产生输出。ARMAX在噪声通路增加 C(q, p) 滤波器。ARARX在噪声通路增加 1/D(q, p) 滤波器。

PPT最后指出一个重要性质：所有这些结构都强制噪声和输入到输出的传递函数共享分母的一部分。具体来说，从真实参数 p^* 的角度，输出可以写为：

注意到 A(q, p^*) 同时出现在输入传递函数和噪声传递函数的分母中。这意味着系统的极点（由 A(q, p^*) 的根决定）同时影响确定性动态和随机动态。这是一个结构约束，在某些情况下可能限制模型的灵活性，但也大大减少了需要估计的参数数量。

模型复杂度选择的基本原则

考虑两个模型结构 M_1(.) 和 M_2(.)，其中 M_2(.) 是 M_1(.) 的超集（super-set）。超集的含义是 M_2(.) 包含 M_1(.) 的所有可能模型，并且还包含更多可能性。因此 M_2(.) 比 M_1(.) 具有更多的自由度。

PPT给出三个例子来说明超集关系。例子1：两个传递函数

H_2 的分母包含二阶项，因此可以描述更复杂的动态（如欠阻尼振荡）。通过设置 p_3 = 0，H_2 退化为 H_1，因此 H_2 是 H_1 的超集。

例子2：M_1(.) 和 M_2(.) 对应阶数递增的ARMAX结构。更高阶的ARMAX可以通过将额外的参数设为零退化为低阶ARMAX，因此高阶模型是低阶模型的超集。

例子3：阶数递增的多项式模型。高阶多项式可以通过将高阶系数设为零退化为低阶多项式。

复杂度与拟合能力的权衡

由于 M_2(.) 具有更多自由度，它能够更精细地描述给定的实验数据（can describe more finement of the given experimental results than M_1(.)）。这在数学上是显然的：更多的参数提供了更大的搜索空间，因此最小化拟合误差时，M_2 的最优误差不会大于 M_1 的最优误差。

因此，人们可能认为唯一需要解决的问题是复杂度与性能之间的权衡（complexity/performance trade-off）。更复杂的模型拟合更好但计算成本更高、参数估计更困难。如果只考虑这一点，选择复杂度就变成了在计算资源约束下尽可能选择复杂模型。

但PPT强调：实际上，问题并不那么简单（In fact, the problem is not so simple）。这是一个关键警示，预示着接下来将揭示复杂度选择的更深层困难。

过拟合问题：噪声数据的陷阱

考虑一个真实结构为 M_1(.) 的过程，但实验数据非常嘈杂。此时会出现一个悖论性的现象。结构 M_1(.) 的最优模型只能近似地表现得像真实过程（only very approximately behaves like the process），因为数据噪声掩盖了真实信号，使得即使用正确的结构也无法完美拟合。

在这种情况下，M_2(.) 结构看似能够获得更好的结果。它确实可以将拟合误差降得更低，看起来更"准确"。但这是一个误导性的表象。M_2(.) 的额外自由度在这里仅仅用于建模噪声的特定实现（only serve here to model a particular realization of the noise），而不是捕捉系统的真实动态。

这里的关键洞察是：噪声是随机的。如果我们重复实验，会得到不同的噪声实现。用 M_2(.) 在第一组数据上估计的模型，其额外参数是为第一次实验的噪声"量身定制"的。当用这个模型预测新数据时，新数据包含不同的噪声实现，那些为旧噪声定制的参数现在变成了累赘，导致预测性能变差。因此 M_2(.) 的最优模型会提供比 M_1(.) 的最优模型更差的预测（will provide a worse prediction than the best model of structure M_1(.)）。

这个现象称为过拟合（overfitting）：模型过度拟合了训练数据的噪声，损失了泛化能力。这揭示了模型复杂度选择的核心困境：在训练数据上表现最好的模型，在新数据上可能表现最差。

过拟合实例分析

PPT给出一个数值例子来具体说明过拟合。真实系统生成的数据遵循规律：

这是一个线性趋势加上随机噪声。选择采样间隔 \Delta = 0.5 秒，使用 i = 0, \ldots, 3 的数据作为训练集，i = 4, \ldots, 8 的数据作为测试集。

考虑三个多项式模型：

对这三个模型进行最小二乘参数估计。第一个图展示了在训练数据上的拟合效果，结果显示最复杂的三阶模型（ordre 3）使输出误差最小化。这符合前面的分析：更多参数可以更好地拟合训练数据。

但第二个图展示了在测试数据上的预测能力，结果完全逆转：最简单的一阶模型（ordre 1）具有最佳的预测能力。三阶模型在训练数据上拟合最好，但外推到新数据时发散得最厉害，预测误差巨大。这是因为三阶模型的高阶项拟合了训练数据中的噪声波动，这些波动在新数据中不复现，导致外推失败。

这个例子直观地证实了理论分析：复杂模型虽然在训练集上表现更好，但泛化能力可能更差。模型复杂度必须与真实系统的复杂度相匹配，过度复杂会导致过拟合。

复杂度选择的统计准则与实际约束

为了避免过拟合，需要使用统计准则（Statistical criteria）来决定何时不值得进一步增加复杂度。这类准则通常基于信息论（如AIC、BIC准则）或交叉验证，平衡拟合精度和模型复杂度。这些方法在统计学习理论中有严格的理论基础。

但PPT指出，统计准则不是唯一的考虑。先验知识（物理学）可能指向更复杂的模型（a priori knowledge (physics) may lead to a more complex model）。如果我们从物理原理知道系统应该包含某些动态特性（如二阶惯性、时滞等），即使简单模型在统计上"足够好"，我们仍可能选择更复杂但物理上合理的模型，以保证模型的外推能力和可解释性。

另一方面，应用约束可能强制使用更简单的模型（application constraints may impose a simpler model）。典型的例子是实时控制：控制器需要在每个采样周期内计算控制律，如果模型过于复杂，计算时间超过采样周期，就无法实时运行。此时必须牺牲一些精度，选择计算更快的简化模型。

频带选择性建模策略

PPT提出一个实用的建模策略：可以寻找一个非常简单的模型，只在给定频率带内（接近临界频率附近）描述过程行为。这种策略的动机是许多控制问题只关心系统在特定频率范围内的响应。例如，在回路整形控制设计中，只需要在交叉频率附近有准确的模型。

为了实现这一点，可以用带通滤波器预处理数据，消除无关信息和减少由偏移和高频噪声引起的误差（eliminate irrelevant information, reduce errors induced by offsets and high-frequency noise）。滤波的两个好处是：第一，去除低频偏移（直流分量）可以避免模型被迫拟合这些与动态无关的常数项；第二，去除高频噪声减少了参数估计的方差，提高了估计精度。

这种方法的哲学是"less is more"：通过限制建模目标（只在特定频带内准确），可以用更简单的模型达到目的。这避免了用复杂模型在全频带建模的困难，同时满足了应用需求。

最后强调，这种频带选择性建模需要实验设计工具（Experiment design tools）的支持。实验设计需要确定：应该在哪些频率激励系统（输入信号的频率内容）、如何设置滤波器参数、需要多长的数据以获得足够的统计精度等。这些问题的答案依赖于对系统的先验知识和建模目标的明确定义。因此实验设计不是建模的附属步骤，而是贯穿整个建模流程的核心环节，从结构选择到数据采集都需要系统化的设计。

模型的结构性质：概念框架

假设我们已经选择了一个模型结构（可能是多个候选结构，需要在其中决定）。在收集数据之前，研究模型结构的性质是很重要的。这里的核心思想是：应该尽可能独立于参数所取的值来研究这些性质（as independently as possible from the values taken by the parameters）。为什么要强调独立于参数值？因为在数据收集之前，我们还不知道参数的真实值，如果模型性质强烈依赖于参数的具体数值，我们就无法在实验设计阶段预测和避免潜在问题。因此需要在数据收集前就检测可能存在的问题。

一个性质被称为结构性质（structural），当且仅当它对几乎所有的参数值都成立，并且可能仅在参数空间的零测度子空间上为假。零测度子空间是测度论中的概念，直观理解就是在参数空间中"极其稀少"的点集合，如果随机选择参数，落在这个集合上的概率为零。因此结构性质是"通用"的性质，不依赖于参数的特殊取值。这样定义的好处是：如果随机抽取参数向量 p，我们几乎必然不会遇到反例（null probability of finding an atypical value）。

在建模背景下，有两个结构性质特别重要：可辨识性（identifiability）和可区分性（distinguishability）。这两个概念将在接下来详细展开，它们决定了参数估计问题是否有意义以及不同模型结构是否能从数据中区分开来。

可辨识性问题的提出

面对参数估计任务，一个自然的问题是：参数 p 的估计是否有成功的机会？这个问题的直观含义是：给定数据，我们能否唯一地确定参数值，或者至少将参数限制在一个小范围内？

但PPT立即指出：如果以这样模糊的方式提问，这个问题没有答案（In such vague terms, the question has no answer）。为什么没有答案？因为"成功的机会"取决于太多因素：数据的质量、数据的数量、噪声的大小、输入信号的选择、模型结构与真实系统的匹配程度等。在没有明确这些条件的情况下，无法给出确定的回答。

为了使问题有意义，需要构建一个理想化框架（idealized framework），在这个框架中去除所有不确定性因素，只关注结构本身的可辨识性。理想化框架包含三个假设。第一，过程和模型具有相同的结构，即不存在特征化误差（no characterization error）。这意味着我们已经选对了模型结构的形式，只是不知道参数值。第二，数据是无噪声的（the data is noiseless）。这去除了随机扰动的影响，使问题变成确定性的。第三，应用的输入和测量时间可以自由选择（the applied input and the measuring times can be freely selected）。这意味着我们可以设计最优的实验，不受采样率、输入幅值、测量时刻等实际约束的限制。

在这个理想化框架下，如果连参数都无法唯一确定，那么在实际情况下（有噪声、有模型误差、输入受限）就更不可能唯一确定了。因此理想化框架给出了可辨识性的一个必要条件。如果理想化框架下可辨识，实际中仍可能因为数据不足或噪声过大而无法准确估计；但如果理想化框架下就不可辨识，那么无论数据多么好都无法唯一确定参数。

在这个理想化框架中，总是可以设置模型参数 \hat{p}，使得模型与过程具有相同的输入输出行为（same input-output behavior），即 M(p) = M(\hat{p})。这里 p 是真实系统的参数（在理想化框架中系统也被假设为同结构模型），\hat{p} 是我们为模型选择的参数。这个等式的含义是：对于所有可能的输入和所有时刻，两个模型产生完全相同的输出。

结构可辨识性的定义与图示

PPT给出了研究结构可辨识性的理想化框架的框图。框图显示输入 u(t) 同时馈入两个模型：一个是参数为真值 p^* 的模型 M(p^*)，另一个是参数为估计值 \hat{p} 的模型 M(\hat{p})。两个模型分别产生输出 y_m(t, p^*, u) 和 y_m(t, \hat{p}, u)，以及感兴趣的量 z_m(t, p^*, u) 和 z_m(t, \hat{p}, u)。两个输出相减得到输出误差 e_s = 0，这是理想化框架的体现：我们希望找到 \hat{p} 使得误差为零。

方程 M(\hat{p}) = M(p^*) 总有解 \hat{p} = p^*（has always the solution \hat{p} = p^*）。这是显然的：如果我们把估计值设为真值，模型自然与真实系统一致。关键问题是：我们希望没有其他解（We hope that there is no other solution）。如果存在其他解 \hat{p} \neq p^* 也满足 M(\hat{p}) = M(p^*)，那么从输入输出数据无法区分 \hat{p} 和 p^*，参数估计就不唯一。

基于这个思想，PPT给出结构全局可辨识性的严格定义。参数 p_i 是结构全局可辨识的（structurally globally identifiable, s.g.i.），如果对于几乎所有 p^* \in \mathbb{P}，方程 M(\hat{p}) = M(p^*) 成立蕴含 \hat{p}_i = p_i^*。这里"几乎所有"再次强调零测度例外的允许。这个定义的含义是：只要知道模型的输入输出行为与真实系统一致，就能唯一确定参数 p_i 的值（至多除了零测度的特殊情况）。

如果模型的所有参数都是s.g.i.，则模型结构 M(.) 本身是s.g.i.。这意味着整个参数向量可以从输入输出行为唯一确定，不存在不同的参数组合产生相同的输入输出关系。

局部可辨识性：弱化的概念

PPT指出一个不幸的现实：M(.) 并不总是s.g.i.。许多实际模型结构存在全局不可辨识性，不同的参数值可能产生相同或几乎相同的输入输出行为。更糟糕的是，即使 M(.) 确实是s.g.i.，证明这一点也可能很困难（even if it is, it may be difficult to prove this）。全局可辨识性的检验通常需要求解非线性方程组，分析解的唯一性，这在高维参数空间中是非常困难的数学问题。

为了应对这些困难，引入一个弱化的概念：结构局部可辨识性（structurally locally identifiable, s.l.i.）。参数 p_i 是s.l.i.，如果对于几乎所有 p^* \in \mathbb{P}，存在 p^* 的某个邻域 \mathbb{N}(p^*)，使得对于该邻域内的 \hat{p} \in \mathbb{N}(p^*)，若 M(\hat{p}) = M(p^*) 则 \hat{p}_i = p_i^*。

局部可辨识性的含义是：参数在真值附近是唯一的，但在远离真值的区域可能存在其他解。换句话说，如果我们的参数估计已经比较接近真值，那么进一步精确化是可能的，不会陷入歧义；但如果初始估计离真值很远，可能存在其他局部极小点使得参数估计收敛到错误的值。

PPT强调：局部可辨识性是全局可辨识性的必要条件（Local identifiability is a necessary condition for global identifiability）。这个逻辑关系很清楚：如果全局可辨识，当然在任何局部邻域内也可辨识；反之，如果连局部都不可辨识，全局肯定也不可辨识。因此局部可辨识性是一个较弱但更容易检验的条件，可以作为初步筛选工具。

如果模型的所有参数都是s.l.i.，则称 M(.) 是s.l.i.。这给出了模型结构可辨识性的层次：s.g.i. > s.l.i. > s.n.i.（结构不可辨识）。

可辨识性的具体例子

PPT通过欧姆定律给出最简单的例子。测量电压 U 和电流 I，模型为 U = RI，参数 p = R。这显然是s.g.i.：给定任意的 (U, I) 数据对（I \neq 0），可以唯一确定 R = U/I。不存在其他 R 值能产生相同的 U 与 I 关系。

但如果选择不同的参数化，情况会改变。假设模型写成 U = \cos(p)I，参数 p 是s.l.i.但不是s.g.i.。为什么？因为余弦函数是周期的，\cos(p) 和 \cos(p + 2\pi k)（k 为任意整数）产生相同的值。因此从数据只能确定 p 模 2\pi 的值，存在无穷多个全局解。但在局部（例如 p \in [0, \pi)），解是唯一的，因此是s.l.i.。这个例子说明可辨识性依赖于参数化方式。

更有趣的是第三个例子。如果 U 被测量（而不是电流 I），并且模型为 U = p_1 p_2，那么参数向量 (p_1, p_2)^T 是s.n.i.（结构不可辨识）。原因很明显：任何满足 p_1 p_2 = U 的参数对都产生相同的输出。存在不可数无穷多个解，例如 (2U, 0.5)、(U, 1)、(0.5U, 2) 等都等价。这个例子说明：如果只测量输出而不测量中间变量，某些参数组合可能无法分离。

这三个例子揭示了可辨识性分析的必要性：仅仅写出模型方程是不够的，必须分析在给定测量配置下参数是否可以唯一确定。最后一句话总结：需要测试方法（Test methods are now required）来系统化地检验可辨识性。

结构不可辨识的定义

在给出测试方法之前，PPT先完整定义结构不可辨识性。参数 p_i 是结构不可辨识的（structurally non-identifiable, s.n.i.），如果对于几乎所有 p^* \in \mathbb{P}，存在不可数无穷多个 \hat{p}_i 值使得 M(\hat{p}) = M(p^*)。不可数无穷的含义是：满足条件的参数值形成一个连续集合（如一个区间或流形），而不仅仅是离散的点。

如果模型至少有一个参数是s.n.i.，则模型 M(.) 是s.n.i.。注意这个定义的不对称性：模型是s.g.i.需要所有参数都是s.g.i.，但模型是s.n.i.只需要至少一个参数是s.n.i.。这种不对称性是合理的，因为任何一个参数的不可辨识都会导致参数估计问题没有唯一解。

可辨识性判断的复杂性

PPT指出一个重要的现实：对于一个s.n.i.的模型，某些参数可能是s.l.i.甚至s.g.i.。这意味着可辨识性是按参数分析的，不同参数可能有不同的可辨识性质。例如在 U = p_1 p_2 的例子中，虽然 (p_1, p_2) 作为向量是s.n.i.，但如果我们引入约束（如 p_1 = 1），则 p_2 变成s.g.i.。

更微妙的情况是：可辨识性有时依赖于参数的数值取值，而非典型区域的测度为零（The identifiability sometimes depends on the numerical value taken by the parameters without the atypical region being of null measurement）。这句话的含义是：存在某些特殊的参数值点或曲线，在这些点上可辨识性丧失，但这些点构成的集合虽然测度为零，却不是"完全可以忽略"的。例如在某个传递函数中，当某个参数恰好使得零点和极点相消时，系统的阶数降低，此时可辨识性分析会改变。

因此，在这种情况下无法给出结构性结论（It is then impossible to give a structural conclusion）。结构性质要求对几乎所有参数值都成立，但如果可辨识性敏感地依赖于参数的数值（例如依赖于某个参数是否为零），那么我们无法在不知道参数值的情况下判断可辨识性。这种依赖于参数值的可辨识性问题在实践中很常见，提醒我们在应用理论结果时需要谨慎。

平稳LI模型的可辨识性测试方法

为了系统化地检验可辨识性，PPT介绍了一个适用于平稳线性输入（stationary LI models）模型的测试方法。这个方法基于传递函数分析，利用线性系统理论的工具。

测试方法包含四个步骤。第一步是评估相关的传递矩阵（Evaluate the associated transfer matrix）：

这是从状态空间表示推导出的输入输出传递函数。传递矩阵 \mathbf{H}(s, p) 完全刻画了系统在频域中的输入输出关系，包含了所有可以从外部观测获得的信息。

第二步是将 \mathbf{H}(s, p) 写成规范形式（canonical form） \mathbf{H}_c(s, p)。规范形式是指一种唯一的、标准化的表示方式。具体要求是：首先通过分子和分母的最大公因子（GCD）简化 \mathbf{H}(s, p) 的每个元素，消除零极点对消；然后将分子和分母按 s 的幂次递增展开；最后将每个分母的最高次项系数归一化为1。这些操作确保了传递函数的表示是最简且唯一的。

第三步是对 \mathbf{H}_c(s, \hat{p}) 和 \mathbf{H}_c(s, p^*) 的系数进行逐项识别（Perform a term-by-term identification）。如果 M(\hat{p}) = M(p^*)，则它们的传递函数必须相等，即 \mathbf{H}_c(s, \hat{p}) = \mathbf{H}_c(s, p^*) 对所有 s 成立。由于是规范形式，这等价于对应次数的系数相等。通过对比每一项的系数，可以得到一组关于 \hat{p} 和 p^* 的代数方程。

第四步是求解得到的方程组（Solve the resulting equations in \hat{p}）。这些方程通常是非线性的，但由于来自传递函数系数匹配，具有特定的结构。求解这些方程，分析解的情况，就可以判断可辨识性。

可辨识性的判断准则

根据方程组解的情况，可以判断模型的可辨识性类型。PPT给出了三种情况的判断准则。

如果对于几乎所有 p^*，\hat{p} 有唯一解且等于 p^*（there is a unique solution for \hat{p}, equal to p^*），则 M(.) 是s.g.i.。这是理想情况：从输入输出行为可以唯一反推参数，不存在歧义。

如果解的集合是可数的（通常是有限的）（the set of solutions is countable (usually finite)），则 M(.) 是s.l.i.。这意味着存在有限个离散的参数值产生相同的输入输出行为。在参数真值的局部邻域内解是唯一的，但在全局存在其他孤立的等价参数值。例如前面提到的 \cos(p) 的例子，解集合是 \{p^* + 2\pi k : k \in \mathbb{Z}\}，这是可数无穷集合。

如果解的集合是不可数的（the set of solutions is uncountable），则 M(.) 是s.n.i.。这表示存在连续的参数自由度无法从数据中确定。例如 U = p_1 p_2 的情况，满足 \hat{p}_1 \hat{p}_2 = p_1^* p_2^* 的参数对 (\hat{p}_1, \hat{p}_2) 形成一个双曲线，是不可数集合。

这三种情况的判断准则提供了清晰的决策规则，使得可辨识性分析变得可操作。

参数级别的可辨识性分析

前面的判断准则针对模型整体，但PPT进一步细化到参数级别。一个参数如果在所有解中总是取相同的值，则是s.g.i.（A parameter which takes always the same value in all solutions is s.g.i.）。这意味着虽然参数向量可能有多个解，但某个特定参数在所有解中都一致，因此该参数是全局可辨识的。

一个参数如果其值在可数集合中取值，则是s.l.i.（A parameter taking its value in a countable set is s.l.i.）。例如，如果所有解中某参数只能取 \{1, -1\} 两个值，则该参数是局部可辨识的。

对于给定的结构 M(.)，不同参数可能具有不同的可辨识性质：某些参数可能是s.g.i.，某些是s.l.i.，某些是s.n.i.。这种混合情况在实践中很常见。可辨识的参数可以从数据中可靠估计，而不可辨识的参数需要通过其他手段确定（如固定为已知值或引入额外的测量）。

PPT特别强调：规范形式是强制性的（The canonical form is mandatory）。这个要求的重要性通过下面的例子说明。

规范形式的必要性：一个反例

考虑传递函数：

如果过快地得出结论，可能会说：\mathbf{H}(s, \hat{p}) = \mathbf{H}(s, p^*) 意味着 \hat{p} = p^*，因此模型是s.g.i.。但这个结论是错误的。

问题在于这不是规范形式。注意到分母没有归一化：最高次项 s 的系数是 p_3 而不是1。将传递函数改写为规范形式：

现在传递函数完全由两个参数组合 p_1/p_2 和 p_3/p_2 决定。因此条件 M(\hat{p}) = M(p^*) 等价于：

这两个方程只能确定两个比值，而不能确定三个参数的绝对值。事实上，存在无穷多个解 \hat{p} 满足这些方程。例如，对于任意非零常数 \lambda，参数 (\lambda p_1^*, \lambda p_2^*, \lambda p_3^*) 都满足条件。这些解构成一条通过原点的射线，是不可数集合，因此 M(.) 是s.n.i.。

这个例子清晰地展示了为什么规范形式是强制性的。如果不先化简和归一化，直接比较系数会遗漏参数之间的内在约束关系，导致错误的可辨识性结论。规范形式通过消除冗余自由度，揭示了传递函数中真正独立的参数组合，从而正确反映可辨识性。

结构可辨识性与可区分性

结构可辨识性 - 测试方法

我们从一个线性时不变但非线性参数化（LI, non-LP）的模型开始分析。这个模型由两个状态变量 x_1 和 x_2 描述，其动态方程为

这里 p_1, p_2, p_3 是待辨识的参数，u 是输入，y_m 是可测量的输出。初始条件都设为零。

为了分析这个系统的可辨识性，我们需要从输入输出关系中提取参数信息。最直接的方法是对系统进行拉普拉斯变换。变换后得到

这里所有变量都是在 s 域中的。现在我们通过消元法来建立输入输出之间的直接关系。从第二个方程可以解出 x_1(s) = \frac{s + p_3}{p_1}x_2(s)，代入第一个方程，经过整理可以消去 x_1 和 x_2，最终得到

这个等式描述了输出 y_m 与输入 u 之间的关系。将其改写为传递函数的标准形式，即输出与输入之比，就得到了典型传递函数

这个传递函数是从输入到输出的完整描述，包含了所有参数信息。

传递函数的参数约束分析

现在考察传递函数的一个关键性质。对于任何传递函数 \mathbf{H}_c(s,\mathbf{p})，如果它对应一个真实的物理系统，那么它必须满足共轭对称性：\mathbf{H}_c(s,\bar{\mathbf{p}}) = \mathbf{H}_c(s,\mathbf{p}^*)。这里 \bar{\mathbf{p}} 表示参数向量 \mathbf{p} 的某种变换，而 \mathbf{p}^* 表示其复共轭。这个关系实际上反映了传递函数在实数域上的对称性要求。

将这个关系应用到我们的传递函数上，通过比较分子和分母的系数，可以得到一组约束方程

第一个方程直接告诉我们 \hat{p}_1 = p_1^*，这意味着 p_1 必须是实数。后两个方程涉及 p_2 和 p_3，它们构成了一个关于 \hat{p}_2 和 \hat{p}_3 的二次方程组。通过韦达定理，如果知道两个数的和与积，这两个数是某个二次方程的两个根。因此，这个方程组在参数空间 \hat{\mathbf{p}} 中有两个解

这两个解的物理意义是：第一个解就是原参数向量本身，而第二个解是将 p_2 和 p_3 互换位置得到的。

这个结果揭示了参数可辨识性的层次：p_1 是结构全局可辨识的（structurally globally identifiable, s.g.i.），因为无论数据如何，它都只有唯一解。而 p_2 和 p_3 虽然可能取两个不同的值，但它们只是结构局部可辨识的（structurally locally identifiable, s.l.i.）。这意味着从无噪声数据中，我们能够确定 p_1 的真实值，但对于 p_2 和 p_3，我们会得到两个可能的值，并且无法仅从输入输出数据来判断哪个是正确的。

结构可辨识性的实际含义

上述分析引出几个关键观察。在无噪声的理想情况下，如果使用迭代算法来辨识参数，算法会收敛到两个解中的某一个。具体收敛到哪个解取决于初始值的选择和算法的路径，但一旦收敛，就会稳定在该解上。

当存在噪声时，问题并不会消失。我们仍然有两个参数向量与相同的输入输出行为相关联。如果知道这两个解中的一个，我们现在就能够推导出另一个解——它就是将 p_2 和 p_3 对调。

这个例子说明了为什么在定义中需要"几乎任何值"（for almost any value）这样的表述。如果参数向量 \mathbf{p}^* 恰好落在某个特殊的参数空间平面上，比如由 p_1 = 0 定义的平面，那么过程的输出会恒等于零，无法提供任何信息。在这种情况下，p_2 和 p_3 会变得完全不可辨识。然而，我们仍然说它们是结构局部可辨识的，因为 p_1 = 0 这种情况是非典型的（atypical），在实际应用中几乎不会遇到。

模型不是结构全局可辨识的这一事实，意味着我们无法从输入输出行为的单独知识中以唯一的方式重构状态 \mathbf{x}。根据选择的模型（即选择了 \hat{\mathbf{p}}_1 还是 \hat{\mathbf{p}}_2），我们会得到 x_1 的两个可能值。

另一方面，如果考虑的模型仅用于预测过程输出，那么可辨识性的概念就变得不相关了。因为两个参数向量产生完全相同的输出预测，选择哪个在预测意义上没有区别。

实用的可辨识性测试方法

测试可辨识性通常需要进行繁琐的代数运算，这些运算在复杂模型中很难手工完成。可以借助计算机代数软件（如MAPLE）来简化这些操作。

然而，有些计算即使在最强大的计算机上也过于复杂。对于这种情况，存在一个实用的数值测试方法来判断模型 M(.) 是否至少是局部可辨识的。这个方法包含以下步骤：首先，在参数空间中随机选择一个名义参数值 \mathbf{p}_0（这里"随机"意味着在参数空间 \mathbb{P} 中任意选取）。然后，用这个参数值 \mathbf{p}_0 对模型 M(\mathbf{p}_0) 进行高精度仿真，生成虚拟测量数据 \mathbf{y}_f。接下来，将这个虚拟数据 \mathbf{y}_f 作为"真实"观测，使用某种估计方法（如牛顿法或高斯-牛顿法）通过最小化输出误差的二次准则来估计参数 \mathbf{p}。估计算法从初始值 \hat{\mathbf{p}}_0 = \mathbf{p}_0 开始迭代。如果估计出的参数 \hat{\mathbf{p}}_k 在 \mathbf{p}_0 处保持稳定（即收敛回 \mathbf{p}_0），则可以认为 M(.) 是结构局部可辨识的。

这个方法的逻辑是：如果模型在某个参数值附近是局部可辨识的，那么从该参数值产生的数据应该能够唯一地反推回这个参数值。如果估计过程不稳定或收敛到其他值，则说明存在局部不可辨识性。

可区分性问题的引入

在建模过程中，我们可能在多个候选模型结构之间犹豫不决。这时自然要问：我们正在考虑的测量数据是否能够帮助我们判断哪个模型结构是最优的？

这就是结构可区分性（distinguishability of structures）问题。这个问题与可辨识性研究在同样的理想化框架下得到部分解答。框架的假设是：真实"过程"是某个结构 M(.) 的一个模型；而我们考虑的"候选模型"具有与 M(.) 不同的结构，记为 \hat{M}(.)。

我们需要引入明确的记号。真实过程结构 M(.) 的参数向量（不一定与候选模型的参数向量维数相同）记为 \mathbf{p}；候选模型结构 \hat{M}(.) 的参数向量记为 \hat{\mathbf{p}}。

可区分性研究的理想化框架可以用框图表示：相同的输入 u(t) 同时送入两个不同结构的系统 M(\mathbf{p}) 和 \hat{M}(\hat{\mathbf{p}})，它们分别产生输出 y_m(t, \mathbf{p}, u) 和 \hat{y}_m(t, \hat{\mathbf{p}}, u)。两个输出的差 e_s 被计算出来。

我们希望找到这样的情况：不存在参数 \hat{\mathbf{p}} 使得 \hat{M}(\hat{\mathbf{p}}) = M(\mathbf{p})。如果能证明这种不可能性，那么就能帮助我们排除候选模型 \hat{M}(.)，因为它永远无法复制真实过程的行为。

结构可区分性的严格定义

更精确地说，模型 \hat{M}(.) 被称为与 M(.) 结构可区分的（structurally distinguishable, s.d.），如果对于几乎所有的参数 \mathbf{p}，都不存在 \hat{\mathbf{p}} 使得

这里的等号意味着两个模型对所有可能的输入产生完全相同的输出。

需要注意一个重要的非对称性：\hat{M}(.) 与 M(.) 结构可区分，并不意味着 M(.) 也与 \hat{M}(.) 结构可区分。换句话说，可能存在这样的情况：第一个模型的类可能是第二个模型类的真子集。如果是这样，那么 \hat{M}(.) 永远无法模拟 M(.) 的所有可能行为（所以它们可区分），但 M(.) 可能包含能够复制 \hat{M}(.) 的所有行为的特殊参数值（所以反向不可区分）。

当 \hat{M}(.) 与 M(.) 结构可区分，并且 M(.) 也与 \hat{M}(.) 结构可区分时，我们说 M(.) 和 \hat{M}(.) 互相结构可区分（mutually s.d.）。

结构可区分性的测试技术

结构可区分性的测试技术与可辨识性测试技术非常相似，但存在两个显著区别。

首先，我们现在考虑的是具有不同结构的两个模型，而不是同一结构下不同参数值的模型。这意味着我们需要比较不同形式的方程或传递函数。

其次，目标发生了变化：在可辨识性研究中，我们希望证明解的唯一性；而在可区分性研究中，我们希望证明方程组 \hat{M}(\hat{\mathbf{p}}) = M(\mathbf{p}) 在 \hat{\mathbf{p}} 中没有解，即证明解的不存在性。

结构可区分性的例子

考虑由以下方程定义的模型 M(.)

以及由下式定义的候选模型 \hat{M}(.)

注意这两个模型的输出都是 x_1，但它们的内部结构不同。

通过拉普拉斯变换和消元，可以得到两个模型的典型传递函数形式。模型 M(.) 的传递函数为

候选模型 \hat{M}(.) 的传递函数为

现在判断这两个模型是否可区分。条件 \hat{M}(\hat{\mathbf{p}}) = M(\mathbf{p}) 等价于两个传递函数相等，即它们的分子和分母系数分别对应相等。这给出方程组

从第一个方程得 \hat{p}_1 = p_3 - \hat{p}_3，代入第二个方程得 \hat{p}_2 = p_1 + p_2。然后代入第三个方程进行验证：(p_3 - \hat{p}_3)(p_1 + p_2) = p_1p_3，整理得 \hat{p}_3(p_1 + p_2) = p_2p_3，即 \hat{p}_3 = \frac{p_2p_3}{p_1 + p_2}。

这意味着对于任何给定的 \mathbf{p}，我们总能找到相应的 \hat{\mathbf{p}} 满足上述方程组。因此 M(.) 和 \hat{M}(.) 是结构不可区分的（structurally undistinguishable）。

实际上，这两个模型虽然内部结构看起来不同，但它们能够产生完全相同的输入输出行为。要消除这种不确定性，一种方法是对 x_2 也进行测量，这样就能从内部状态的差异来区分两个模型。

可区分性分析的注意事项

与可辨识性类似，在某些情况下无法得出明确的结构结论。考虑两个具有传递函数的模型

其中 \hat{\mathbf{p}} 和 \mathbf{p} 都属于 \mathbb{R}^2。

如果 \hat{\mathbf{p}} 的取值使得 \hat{\mathbf{H}}_c(s,\hat{\mathbf{p}}) 有两个实数极点，那么通过选择合适的 \mathbf{p} 可以使 \mathbf{H}_c(s,\mathbf{p}) 产生完全相同的极点配置，因此 \mathbf{H}_c(s,\mathbf{p}) 与 \hat{\mathbf{H}}_c(s,\hat{\mathbf{p}}) 不可区分。

反之，如果 \hat{\mathbf{H}}_c(s,\hat{\mathbf{p}}) 的极点是复共轭对，那么无论如何选择 \mathbf{p}，\mathbf{H}_c(s,\mathbf{p}) 都只能产生实数极点，因此两者可区分。

问题在于，这两种情况都不能被认为是典型的。参数空间被分成了两个区域，在一个区域内模型可区分，在另一个区域内不可区分，而这两个区域的边界（即 \hat{\mathbf{H}}_c 有重实根的情况）是一个测度为零的集合。这种情况下，我们无法对模型做出一般性的可区分性判断。

可辨识性与可区分性的关系

两个结构的可辨识性既不是它们可区分性的必要条件，也不是充分条件。一个模型可能是全局可辨识的但与另一个模型不可区分；反过来，两个模型可能互相可区分但各自都不是全局可辨识的。

尽管如此，用于测试这两类性质的技术相当相似。它们都基于传递函数或等价的输入输出关系的代数操作，都需要分析方程组的解的存在性和唯一性，都在同样的理想化框架（无噪声、充分激励等）下进行。

总结性思考

关于模型结构选择，我们需要认识到这不是一门精确的科学。它涉及或多或少的任意决策，但这些决策会带来重要的后果——最佳的可能模型必须在所考虑的模型类别中寻找。因此，任何选定的模型结构都必须被视为临时性的、可质疑的。

一个模型结构可以在进行任何实际测量之前就受到质疑，方法就是通过可辨识性或可区分性研究。这些研究能够在数据收集之前就揭示模型的内在局限性。

当参数或状态变量具有物理意义时（即基于知识的模型），可辨识性和可区分性的概念是根本性的。因为在这种情况下，我们关心的不仅是输入输出行为，还关心内部机制的正确表示和参数的物理解释。

相反，当我们只对再现输入输出行为感兴趣时（即行为模型），这些概念的重要性就变得边缘化了。在纯粹的预测场景中，只要模型能准确预测输出，其内部参数是否唯一确定或模型结构是否能与其他结构区分开来，都不是关键问题。

参数估计的代价函数

课程概览

本课程主要讨论参数估计中使用的各种代价函数。课程内容包括：最小二乘法、最小绝对值准则（L_1 范数）、极小极大准则（L_\infty 范数）以及最大似然估计。这些不同的代价函数对应着不同的优化目标和数学性质，适用于不同的实际问题场景。

参数估计问题的基本框架

参数估计的核心问题可以用框图来理解：真实系统 S 受到输入 u(t) 和扰动 b(t) 的共同作用，产生输出 y(t)。我们建立一个参数化模型 M(\mathbf{p})，它接收相同的输入 u(t)，产生模型输出 y_m(\mathbf{p};t)。将真实输出 y(t) 与模型输出 y_m(\mathbf{p};t) 进行比较，得到误差信号，这个误差被送入某个准则（代价函数）进行评估。优化器根据代价函数的值 z(\mathbf{p};t) 来调整参数 \mathbf{p}，目标是使代价函数最小化。

问题的一般性描述

设 \mathbf{y} 是一个向量，它收集了对真实系统 S 的所有观测数据。\mathbf{y}_m(\mathbf{p}) 是当参数取值为 \mathbf{p} 时，模型对这些观测条件下的响应预测。参数估计的目标是选择参数 \mathbf{p}，使得 \mathbf{y}_m(\mathbf{p}) 在某种意义下"看起来像" \mathbf{y}。

这个"看起来像"的具体含义需要数学化，即最小化 j(\mathbf{p})，这个 j(\mathbf{p}) 衡量的是 \mathbf{y}_m(\mathbf{p}) 与 \mathbf{y} 之间的距离。用符号表示就是：选择 \mathbf{p} 使得估计出的参数记为 \hat{\mathbf{p}}(\mathbf{y})。

这里隐含着两个关键问题。第一个问题（Rq1）：由于 \mathbf{y} 是随机的（包含测量噪声、系统扰动等随机因素），因此 \hat{\mathbf{p}}(\mathbf{y}) 也是随机的。这意味着估计出的参数本身具有统计特性，包括偏差（bias）和方差（variance）等。第二个问题（Rq2）：如何优化函数 j(.)？这涉及优化算法的选择和实现，后续会讨论。

最小二乘准则

最小二乘法使用的是二次距离，也称为 L_2 范数。代价函数定义为

这里 i 是测量的索引。如果系统输出是以周期 T 进行采样的，那么 i 对应观测时刻 t_i，其中 t_i = Ti。更一般地，i 可以对应不同的实验条件 \xi_i，比如不同的输入序列、不同的初始条件等。

这个准则的物理意义是：我们希望模型预测值与实际观测值的偏差平方和最小。平方的使用有几个原因：它惩罚大的偏差比小偏差更严重（因为是二次增长），同时数学上可微，便于优化。

加权最小二乘

在某些情况下，不同的观测数据点可能具有不同的重要性或可靠性。这时可以引入权重矩阵 \mathbf{W}，得到加权最小二乘准则

权重矩阵是对角矩阵 \mathbf{W} = \text{diag}(w_i, i = 1,\ldots,N)，其中每个权重满足 w_i \geq 0，对所有 i 成立。这里的一个实际问题是：如何选择这些权重 w_i？权重的选择通常基于对测量精度的先验知识。如果某个测量点的噪声较大，我们会给它分配较小的权重，反之亦然。

最小二乘的几何解释 - 局部极小问题

从几何角度看，最小二乘估计 \hat{\mathbf{p}} 的含义是：\mathbf{y}_m(\hat{\mathbf{p}}) 是观测向量 \mathbf{y} 在模型响应曲面上的正交投影。这个曲面定义为 \mathcal{S}_r = \{\mathbf{y}_m(\mathbf{p}), \mathbf{p} \in R^p\}，即所有可能的参数值对应的模型响应的集合。

在理想情况下，如果这个曲面 \mathcal{S}_r 是平的（或者说是凸的），那么最小二乘问题只有一个全局最优解。从图中可以看到，在二维输出空间 (y_1, y_2) 中，观测点 \mathbf{y} 到曲面的最短距离点就是 \mathbf{y}_m(\hat{\mathbf{p}})，这个点通过垂直线连接。

但实际问题中，当曲面 \mathcal{S}_r 是弯曲的（非凸的），就可能出现局部极小点。下面的图展示了这种情况：存在另一个参数值 \hat{\mathbf{p}}'，使得 \mathbf{y}_m(\hat{\mathbf{p}}') 也是 \mathbf{y} 到曲面的一个局部最短距离点。虽然 \mathbf{y}_m(\hat{\mathbf{p}}) 到 \mathbf{y} 的距离可能比 \mathbf{y}_m(\hat{\mathbf{p}}') 到 \mathbf{y} 的距离更短（全局最优），但如果优化算法从某个初始值开始迭代，可能会收敛到 \hat{\mathbf{p}}' 这个局部最优点而无法跳出。这是非线性参数估计问题的一个根本困难。

最小绝对值准则（L_1 范数）

与最小二乘不同，最小绝对值准则使用的是 L_1 范数

加权版本为

这个准则有两个重要的特性需要注意。第一，它在任何地方都不可微。具体来说，绝对值函数 |\cdot| 在零点处不可微，这意味着传统的基于梯度的优化方法（如梯度下降法）无法直接应用。第二，即使模型是结构全局可辨识的（s.g.i.），最优解也不一定唯一。

通过一个例子可以看出非唯一性。假设 y_m(\mathbf{p}, 1) = p_1 + p_2^2，y_m(\mathbf{p}, 2) = y_m(\mathbf{p}, 3) = p_1 + p_2，观测值为 \mathbf{y} = (5, 2, 4)。那么 j_{LAV}(\mathbf{p}) = 2 对所有满足 p_1 + p_2^2 = 5 且 2 \leq p_1 + p_2 \leq 4 的参数 \mathbf{p} 都成立。这意味着存在无穷多个参数组合都能达到相同的最小代价值。

因此，虽然 L_1 准则在某些应用中很有用（特别是对离群值具有鲁棒性），但它不能直接使用，需要特殊的处理方法。

L_1 准则的线性规划求解

当模型是线性参数化的，即 y_m(\mathbf{p}, i) = \mathbf{r}_i^T\mathbf{p} + c_i 时，最小绝对值问题可以转化为线性规划问题。目标是最小化 \sum_{i=1}^{N}\alpha_i，约束条件为

这里引入了 N 个辅助变量 \alpha_i，它们代表了每个残差的绝对值。这样，原问题变成了一个具有 N + p 个变量（\alpha 和 \mathbf{p}）以及 2N 个约束的线性规划问题。

关键观察是：约束条件和代价函数都是关于变量线性的。这使得可以使用成熟的线性规划算法（如单纯形法）来高效求解。虽然引入了额外的变量，但线性规划的计算复杂度相对可控，特别是对于大规模问题。

极小极大准则（L_\infty 范数）

极小极大准则关注的是最大误差，即所有残差中的最坏情况。无权重版本定义为

这个准则寻找的是使得最大残差最小的参数。换句话说，它不关心残差的平均水平或总和，而是要控制最坏的那个数据点的偏差。加权版本为

当模型是线性参数化的，即 y_m(\mathbf{p}, i) = \mathbf{r}_i^T\mathbf{p} + c_i 时，极小极大问题也可以转化为线性规划。目标是最小化单个变量 \alpha，约束条件为

这里只引入了一个辅助变量 \alpha，它的值必须大于等于所有残差的绝对值，因此 \alpha 最终等于最大残差。整个问题变成了具有 1 + p 个变量（\alpha 和 \mathbf{p}）以及 2N 个约束的线性规划。相比 L_1 准则，这里的变量数更少，但约束数量相同。

极小极大准则在对最坏情况敏感的应用中很有用，比如当我们需要保证所有数据点的拟合质量都在某个界限内时。

最大似然估计的基本概念

最大似然估计从概率的角度来看待参数估计问题。我们假设观测变量 \mathbf{y} 是一个随机变量，它由确定性模型加上参数 \bar{\mathbf{p}}（真实但未知的参数）以及一个扰动序列 (\epsilon_k)_k 生成。这个扰动序列是独立同分布的随机变量序列，可能经过了系统的某个部分的滤波。

对于给定的观测数据 \mathbf{y}，不同的参数值 \mathbf{p} 会产生不同的观测概率。\mathbf{y} 的似然函数定义为 \pi(\mathbf{y}|\mathbf{p})，它表示在参数为 \mathbf{p} 的条件下观测到 \mathbf{y} 的概率密度。注意这里的符号：\mathbf{y} 是固定的（已经观测到的数据），而 \mathbf{p} 是变量。

最大似然估计的思想是：既然我们观测到了 \mathbf{y}，那么应该选择使得这个观测结果出现概率最大的参数值。因此，最大似然估计 \hat{\mathbf{p}}_{ML} 定义为使 \pi(\mathbf{y}|\mathbf{p}) 最大的参数值。

在实际计算中，通常不直接最大化似然函数，而是最大化对数似然函数 \log\pi(\mathbf{y}|\mathbf{p})。这样做有几个原因：对数函数是单调递增的，所以最大化 \log\pi 等价于最大化 \pi；对数可以将概率的乘积变成和，简化计算；对数函数在数值计算上更稳定，避免了概率值非常小时的数值下溢问题。

似然函数的图示理解

通过一个图示可以更直观地理解似然函数。假设我们有两个候选参数值 \mathbf{p}_1 和 \mathbf{p}_2，对应的似然函数分别是 \pi(\mathbf{y}|\mathbf{p}_1) 和 \pi(\mathbf{y}|\mathbf{p}_2)。横轴代表观测值 y，纵轴代表概率密度 \pi(\mathbf{y}|\mathbf{p})。

图中显示了两条概率密度曲线，一条对应 \mathbf{p}_1，另一条对应 \mathbf{p}_2。假设实际观测值为"y observé"（观测到的y）。在这个观测点处，\pi(\mathbf{y}|\mathbf{p}_1) 的值高于 \pi(\mathbf{y}|\mathbf{p}_2) 的值，这意味着参数 \mathbf{p}_1 能够更好地解释观测数据，因此 \mathbf{p}_1 更有可能是真实参数。最大似然估计就是在所有可能的参数值中，找到使得似然函数在观测点处取得最大值的那个参数。

最大似然估计的例子

第一个例子（Ex1）考虑了一个简单但基础的情况。假设有 N 个观测值可用，记为 i = 1,\ldots,N。观测模型为 y(i) = \epsilon_i \sim \mathcal{N}(\mu, \sigma^2)，即每个观测值都是均值为 \mu、方差为 \sigma^2 的高斯随机变量。参数向量为 \mathbf{p} = (\mu, \sigma^2)。

这里提出了几个问题。首先，最大似然估计 \hat{\mathbf{p}}_{ML} 是什么？其次，这个估计的偏差（Bias）是多少，即 \mathbb{E}[\hat{\mathbf{p}}(\mathbf{y})|\mathbf{p}] - \mathbf{p} 等于多少？第三，方差是多少，即 \mathbb{E}\{[\hat{\mathbf{p}}(\mathbf{y}) - \mathbb{E}\{\hat{\mathbf{p}}(\mathbf{y})|\mathbf{p}\}][\hat{\mathbf{p}}(\mathbf{y}) - \mathbb{E}\{\hat{\mathbf{p}}(\mathbf{y})|\mathbf{p}\}]^T|\mathbf{p}\} 等于多少？

虽然后续的例子（Ex2和Ex3）内容显示为灰色，但它们涉及更复杂的情况。Ex2考虑 y(i) = y_m(\mathbf{p}, i) + \epsilon_i，其中 (\epsilon_i)_i 独立同分布服从 \mathcal{N}(0, \sigma_i^2)。这涉及三种子情况：方差 \sigma_i^2 已知、所有方差相等但未知、方差未知。Ex3则考虑多维输出的情况，y(i) = y_m(\mathbf{p}, i) + \epsilon_i，其中 \mathbf{y}_i \in \mathbb{R}^n，(\epsilon_i)_i 独立同分布服从 \mathcal{N}(\mathbf{0}, \Sigma)。

模型选择与信息准则

当我们面对多个竞争的模型结构 \mathcal{M}_i 时（例如参数数量 n_p 递增的一系列模型），需要同时解决两个问题：确定最佳结构，以及确定该结构的参数值。

针对这个问题，赤池信息准则（Akaike Information Criterion, AIC）定义为

贝叶斯信息准则（Bayesian Information Criterion, BIC）定义为

这两个准则都需要最小化。它们的基本思想是在模型拟合度（由负对数似然项衡量）和模型复杂度（由参数数量项衡量）之间寻找平衡。具体做法是：对每个结构 \mathcal{M}_i 先进行最大似然估计得到 \hat{\mathbf{p}}_i，然后选择使 j_{AIC}(\hat{\mathbf{p}}_i, i) 或 j_{BIC}(\hat{\mathbf{p}}_i, i) 最小的模型。

这两个准则的第二项是惩罚项。AIC使用线性惩罚 n_p，而BIC使用对数惩罚 n_p\log(N)。当样本量 N 较大时，BIC对复杂模型的惩罚更重，倾向于选择更简单的模型。其他代价函数也被提出，但核心思想都是对复杂结构进行惩罚。

最大似然估计的渐近性质

最大似然估计为什么值得关注？这需要从其统计性质来回答。在一系列假设条件下，最大似然估计具有优良的渐近性质。

这些假设包括：数据由真实参数 \mathbf{p}^* 的模型 \mathcal{M}(\mathbf{p}^*) 生成；扰动 (\epsilon_i)_i 是独立同分布的；似然函数关于 \mathbf{y} 的导数 \{\mathbf{y}/\pi(\mathbf{y}|\mathbf{p}) > 0\} 与参数 \mathbf{p} 独立；对数似然的二阶偏导数 \frac{\partial^2\log\pi(\mathbf{y}|\mathbf{p})}{\partial p\partial p} 存在，关于 \mathbf{p} 连续，关于 \mathbf{y} 一致连续；一阶导数的期望平方有界，即 \mathbb{E}\left\{\left|\frac{\partial\log\pi(\mathbf{y}|\mathbf{p})}{\partial p_i}\right|^2\bigg|\mathbf{p}\right\} < \infty 对所有 i = 1,\ldots,p 和所有 \mathbf{p} 成立；二阶混合偏导数的期望也有界，即 \mathbb{E}\left\{\left|\frac{\partial^2\log\pi(\mathbf{y}|\mathbf{p})}{\partial p_i\partial p_j}\right|\bigg|\mathbf{p}\right\} < \infty 对所有 (i,j) 和所有 \mathbf{p} 成立。

在这些条件下，当样本量 N 趋于无穷时，最大似然估计 \hat{\mathbf{p}}_{ML}^N 的分布渐近服从

其中 \mathbf{M}(\mathbf{p}^*, \xi) 是费舍尔信息矩阵（Fisher information matrix），它反映了观测数据中包含的关于参数的信息量。这个矩阵定义为

或者等价地

这里 \xi 刻画了所考虑的实验条件。这个结果表明，在大样本情况下，最大似然估计是渐近无偏的（asymptotically without bias）且渐近有效的（asymptotically efficient）。

费舍尔信息矩阵的推导

为了理解费舍尔信息矩阵的两个等价形式，考虑模型 y(i) = y_m(\mathbf{p}, i) + \epsilon_i，其中 (\epsilon_i)_i 是独立同分布的。如果不满足这个条件，可以考虑预测误差及其分布 \pi[\epsilon(\mathbf{p}, i)|\mathbf{p}]。

对数似然关于参数的梯度可以写为

在费舍尔信息矩阵的定义中，梯度的外积的期望为

关键观察是，梯度的期望等于零。这可以通过以下推导看出：

最后一步成立是因为概率密度函数的积分恒等于1。因此信息矩阵的第一个形式简化，可以得到

Cramér-Rao不等式

Cramér-Rao不等式是参数估计理论中的一个基本结果，它给出了任何无偏估计量方差的下界。对于任何无偏估计器 \hat{\mathbf{p}}（即满足 \mathbb{E}\{\hat{\mathbf{p}}(\mathbf{y})|\mathbf{p}\} = \mathbf{p}），其协方差矩阵满足

这里的符号 \mathbf{A} \succeq \mathbf{B} 表示 \mathbf{A} - \mathbf{B} 是非负定矩阵，即对于任意向量 \mathbf{v}，都有 \mathbf{v}^T(\mathbf{A} - \mathbf{B})\mathbf{v} \geq 0。

这个不等式的重要意义在于：它告诉我们无论使用什么样的无偏估计方法，估计量的方差都不可能小于费舍尔信息矩阵的逆。这个下界只依赖于数据的统计模型和样本量，与具体的估计方法无关。前面我们证明的最大似然估计的渐近分布恰好达到了这个下界，这说明最大似然估计在大样本情况下是最优的——它达到了理论上的最小方差。

Cramér-Rao不等式的证明要素

证明这个不等式需要构造特定的矩阵并证明其非负定性。首先定义向量函数

这个向量包含两部分：估计误差和对数似然的梯度。然后计算这个向量的协方差矩阵

这个分块矩阵的右上角是单位矩阵 \mathbf{I}_p，这来自于估计量无偏性的条件。由于任何协方差矩阵都是半正定的，所以 \mathbf{A} \succeq \mathbf{O}。

接下来的关键步骤是计算另一个矩阵

展开这个矩阵乘法，可以得到 \text{Cov}(\hat{\mathbf{p}}) - \mathbf{M}^{-1}(\mathbf{p}, \xi)/N \succeq \mathbf{O}，这正是Cramér-Rao不等式。证明的核心思想是利用了二次型的非负性以及矩阵的合同变换保持半正定性这一性质。

积分内求导的数学工具

证明过程中需要交换积分和求导的顺序。这里给出了Leibniz积分法则的一个版本。如果函数 f:[a,b] \times \mathcal{S} \to \mathbb{R}（其中 -\infty < a < b \leq \infty，\mathcal{S} \subseteq \mathbb{R}）满足以下条件：对所有 x \in [a,b] 和 y \in \mathcal{S}，偏导数 \frac{\partial f(x,y)}{\partial y} 存在；对所有 y \in \mathcal{S}，函数 f(x,y) 和 \frac{\partial f(x,y)}{\partial y} 关于 x 在 [a,b] 上可积；对所有 x \in [a,b] 和 y \in \mathcal{S}，存在可积函数 g(x) 使得 \left|\frac{\partial f(x,y)}{\partial y}\right| \leq g(x)。

那么可以交换积分与求导的顺序

这个定理确保了我们在推导费舍尔信息矩阵时进行的积分求导操作是合法的。

贝叶斯估计的基本框架

贝叶斯估计与最大似然估计的根本区别在于对参数 \mathbf{p} 的认识论。在最大似然框架中，参数被视为未知但固定的常数；而在贝叶斯框架中，参数本身被视为随机变量，具有概率分布。

在贝叶斯方法中，我们对参数有一个先验分布（a priori pdf）\pi(\mathbf{p})，它反映了在看到数据之前对参数的认识或信念。观测到数据 \mathbf{y} 之后，我们得到参数的后验分布（a posteriori probability）\pi(\mathbf{p}|\mathbf{y})，它综合了先验信息和数据信息。

图中展示了这个过程：横轴是参数空间，纵轴是概率密度。先验分布 \pi(\mathbf{p}) 是一条较为平坦的曲线，表示我们初始时对参数的不确定性。参数的真实值可能在 \hat{\mathbf{p}}_0 附近，有一定的不确定性范围 \sigma_0。在观测到数据后，后验分布 \pi(\mathbf{p}|\mathbf{y}) 变成了一条尖锐的峰，集中在 \hat{\mathbf{p}}_{MAP} 附近。这个峰的位置综合了先验信息和似然信息，通常会介于先验均值和最大似然估计之间。

贝叶斯规则与MAP估计

连接先验、似然和后验的是贝叶斯规则

分子是联合分布，可以分解为似然函数乘以先验。分母 \pi(\mathbf{y}) 是边缘分布，它作为归一化常数保证后验分布积分为1。

最大后验估计器（Maximum a posteriori estimator, MAP）定义为使后验分布最大的参数值 \hat{\mathbf{p}}_{MAP}。由于 \pi(\mathbf{y}) 不依赖于 \mathbf{p}，最大化 \pi(\mathbf{p}|\mathbf{y}) 等价于最大化

这个表达式包含两项：第一项是对数似然，第二项是对数先验。因此MAP估计可以看作是在最大似然估计的基础上加入了先验信息的正则化项。当先验是均匀分布（无信息先验）时，MAP估计退化为最大似然估计。

贝叶斯估计的渐近性质

在与最大似然估计相同的假设条件下，如果先验分布是正态分布 \mathcal{N}(\bar{\mathbf{p}}^0, \Omega)，则MAP估计具有与最大似然估计相同的渐近性质。具体地说，\hat{\mathbf{p}}_{MAP} 最大化

第二项是正态先验对数密度的二次形式。

当模型是线性参数化的，即 \mathbf{y}(\mathbf{p}) = \mathbf{R}\mathbf{p} + \epsilon，其中 \epsilon \sim \mathcal{N}(\mathbf{0}, \Sigma)，那么MAP估计有显式解

这个公式展示了数据信息（通过 \mathbf{R}^T\Sigma^{-1}\mathbf{y}）和先验信息（通过 \Omega^{-1}\bar{\mathbf{p}}^0）如何被加权组合。数据的权重是 \mathbf{R}^T\Sigma^{-1}\mathbf{R}（与观测精度相关），先验的权重是 \Omega^{-1}（与先验不确定性的逆相关）。当观测数据很多且精确时，第一项主导；当先验很强（\Omega 很小）时，第二项主导。

最小风险估计

贝叶斯框架还允许我们定义更一般的估计准则。最小风险估计器（Minimum risk estimator）的思想是：定义一个代价函数 j(\mathbf{p}|\mathbf{p}^*)，它衡量当真实参数是 \mathbf{p}^* 时用 \hat{\mathbf{p}} 估计带来的损失。由于 \mathbf{p}^* 本身是随机的，我们对这个代价取期望

这可以进一步改写为

这个积分综合了三个因素：估计误差的代价、数据的似然以及参数的先验分布。

二次代价函数下的最小风险估计

当代价函数是二次形式，即 j(\mathbf{p}|\mathbf{p}^*) = (\mathbf{p} - \mathbf{p}^*)^T\mathbf{Q}(\mathbf{p} - \mathbf{p}^*)，其中 \mathbf{Q} 是对称正定矩阵，那么可以通过对 j_{RM}(\mathbf{p}) 关于 \mathbf{p} 求导并令其为零来找到最优估计。驻点条件为

由于 \mathbf{Q} 可逆，这等价于

这个结果表明：在二次代价函数下，最小风险估计就是后验分布的均值。这是一个重要的结论，因为后验均值在计算上通常比后验众数（MAP）更稳定。如果后验分布是对称的，那么后验均值等于MAP估计，即 \hat{\mathbf{p}}_{RM} = \mathbb{E}\{\mathbf{p}|\mathbf{y}\} = \hat{\mathbf{p}}_{MAP}。

鲁棒估计的动机

在实际数据中，离群值（outliers）是一个常见且严重的问题。通过一个具体例子可以看出其影响：考虑观测模型 y(i) = p_1\exp(-p_2t_i) + \epsilon_i，其中真实参数为 \mathbf{p} = (1, 2)^T，噪声 (\epsilon_i)_i 是独立同分布的高斯噪声，标准差为 \sigma = 5 \times 10^{-3}。在时间区间 [0,1] 内采集了200个数据点。

但是数据中存在离群值：从 i = 11 到 i = 30 以及从 i = 50 到 i = 64 的观测值被人为设置为零。图中展示了两条曲线：虚线表示最大似然估计（实际上就是最小二乘估计）得到的模型响应 y_m(\hat{\mathbf{p}}_{MC}, i)，实线表示使用鲁棒方法（MVA方法）得到的模型响应 y_m(\hat{\mathbf{p}}_{MVA}, i)。

可以看到，最小二乘估计受到离群值的严重影响，拟合曲线被拉向零值点，偏离了真实的指数衰减趋势。而鲁棒估计方法能够识别并降低离群值的影响，得到的曲线更接近真实模型。这说明当数据质量不理想时，标准的最大似然估计可能失效，需要使用对离群值不敏感的估计方法。

M-估计器的一般框架

最小二乘和最小绝对值准则实际上是一类更广泛的估计方法的特例，这类方法被称为M-估计器（M-estimators）。M-估计器的一般形式为

其中 e(\mathbf{p}, i) = y(i) - y_m(\mathbf{p}, i) 是输出误差（如果模型不同，也可以是其他类型的残差）。关键在于函数 \rho(\cdot) 的选择，不同的 \rho 函数对应不同的估计器，也决定了估计器对离群值的敏感程度。

两个具有代表性的M-估计器是Huber估计器和Tukey估计器。它们的 \rho 函数都包含一个阈值参数 \delta，用于区分正常数据和可能的离群值。

左图展示了 \rho_1(e)（实线）及其导数 |d\rho_1(e)/de|（虚线）的形状。可以看到，在误差较小时，\rho_1 是二次的（类似最小二乘），但当误差超过阈值后，增长速度变慢。右图展示了 \rho_2(e) 及其导数。\rho_2 在小误差时也是二次的，但在大误差时会达到一个上界并保持平坦，这意味着极端离群值对代价函数的贡献被限制住了。

Huber估计器

Huber估计器的 \rho 函数定义为

这个函数的设计思想是在小误差范围内保持二次性质（即保持高斯噪声下的最优性），而在大误差时切换到线性增长。为什么这样设计？因为二次函数对离群值过于敏感——如果一个残差是 10\sigma（10倍标准差），它对代价函数的贡献是 100\sigma^2，这会严重扭曲估计结果。而线性增长意味着极端值的影响被大大削弱。

阈值 \delta 的选择需要权衡：\delta 太小会把正常的大噪声点也当作离群值处理，损失效率；\delta 太大则无法有效抑制离群值。通常 \delta 设置为噪声标准差的1.5到3倍。

Tukey估计器

Tukey估计器（也称为双权重或biweight估计器）采用更激进的策略

这个函数在 |e| > \delta 时变为常数，意味着超过阈值的残差对代价函数的贡献完全固定，进一步的增大不再有任何影响。这相当于完全忽略了离群值。在 |e| \leq \delta 范围内，函数形状经过精心设计，在零点附近近似二次，但随着误差增大逐渐过渡到平坦。

从导数的角度看更清楚：Huber的导数在 |e| > \delta 时仍然是常数（线性增长的导数），而Tukey的导数在超过阈值后直接变为零，这意味着离群值在优化过程中不再提供任何梯度信息，被完全屏蔽了。

崩溃点：离群值何时压倒估计

M-估计器虽然对离群值有一定抵抗力，但这种抵抗力是有限的。崩溃点（breakdown point）的概念量化了这个限度。它定义为：当数据中有多大比例的离群值时，估计器会完全失效（产生任意大的偏差）。

设原始数据为 \mathbf{y}，被污染的数据为 \mathbf{y}'(\alpha)，其中包含 \alpha\% 的离群值。估计器 \hat{\mathbf{p}} 基于数据的偏差定义为

这个定义考虑了最坏情况：在所有包含 \alpha\% 离群值的数据集中，哪个导致的估计偏差最大。

崩溃点则定义为

即使得偏差变为无穷大的最小离群值比例。

图中展示了这个概念：纵轴是观测值 y(i)，横轴是索引或时间。大部分数据点（标记为+）沿着一条平滑的曲线 y_m(\hat{\mathbf{p}}(\mathbf{y})) 分布，但有一个点 y(1/7) 远离曲线。如果离群值的数量继续增加，最终会达到一个临界点，使得估计器无法区分哪些是真实数据，哪些是离群值。

对于标准的最小二乘估计，\text{BP}[\mathbf{y}, \hat{\mathbf{p}}_{MC}(\cdot)] = \text{BP}[\mathbf{y}, \hat{\mathbf{p}}_{MVA}(\cdot)] = 1/N \to 0 当 N \to \infty。这意味着即使只有一个离群值（比例 1/N），在样本量很大时，就可能导致最小二乘估计严重偏离。这解释了为什么最小二乘法对离群值如此脆弱。

误差平方的中位数

一个具有高崩溃点的鲁棒方法是使用误差平方的中位数（median of squared errors）作为代价函数

这里"med"表示中位数操作。中位数的定义是：将 N 个平方误差按从小到大排序，如果将它们标记为 e_s^2(\mathbf{p}, 1) \leq e_s^2(\mathbf{p}, 2) \leq \cdots \leq e_s^2(\mathbf{p}, N)，那么中位数是排序后序列的中间值。

图示展示了这个概念：横轴代表误差平方 e^2(\mathbf{p}, i)，每个×号代表一个数据点的误差平方值。j_{MED}(\mathbf{p}) 对应的位置在中间，意味着有一半的误差平方比它小，一半比它大。

具体地，如果 h 个最小的平方误差满足 e_s^2(\mathbf{p}, 1) \leq e_s^2(\mathbf{p}, 2) \leq \cdots \leq e_s^2(\mathbf{p}, N)，那么 j_{sh}(\mathbf{p}) = \sum_{i=1}^{h}e_s^2(\mathbf{p}, i)。当 h 选择为 [N/2]（N 的一半，向下取整）时，这种方法称为最小中位数平方（Least Median of Squares, LMS）。

虽然PPT中有些内容标注为灰色，但可以推断出中位数方法的崩溃点接近50%。这是因为只有当超过一半的数据是离群值时，中位数才会被离群值主导。相比之下，均值（最小二乘使用的）只要有少量大的离群值就会被显著影响。这解释了中位数方法的鲁棒性来源。

图像处理中的应用

这些鲁棒估计方法在图像处理中有直接的应用，特别是图像配准（image registration）问题。假设有一个变换 T_{\mathbf{p}}，它将图像 A 映射到图像 T_{\mathbf{p}}(A)，其中 \mathbf{p} 代表变换参数，可能包括平移、旋转、缩放、伽马校正等。

给定两幅图像 A 和 B，自动配准的目标是找到变换参数 \mathbf{p}，使得变换后的 B 看起来像 A（这在医学成像中很常见，比如对齐不同时间或不同模态的扫描图像）。传统方法是将 B 视为 T_{\mathbf{p}}(A) 加上噪声，然后估计 \mathbf{p} 使得 T_{\mathbf{p}}(A) 与 B 之间的相似度最大化。

但问题在于：两幅图像之间可能存在显著的差异，这些差异可能无法用变换 T_{\mathbf{p}} 解释，它们相当于离群值。比如，图像中可能有某些区域的内容完全不同，或者有遮挡、阴影等。如果使用标准的最小二乘方法（最大化相似度通常等价于最小化像素值差的平方和），这些离群区域会严重干扰参数估计。

一个鲁棒的解决方案是：将图像 A 和 B 转换为一维信号（例如通过光栅扫描），然后估计参数 \hat{\mathbf{p}} 使得符号变化次数最大

这里"符号变化"是指相邻像素的差值从正变负或从负变正的次数。图示说明了这个想法：下方的曲线显示了两种情况。当配准不佳时（1个符号变化），差值信号几乎单调，符号变化很少。当配准良好时（6个符号变化），差值信号在零附近快速振荡，符号频繁变化。

为什么这个方法鲁棒？因为它只关注符号，而不关心差值的大小。即使某些区域的差值非常大（离群值），只要大部分区域配准良好，符号变化次数仍然会很高。这个方法易于计算，并且对离群值具有天然的鲁棒性，使其在实际图像处理任务中非常实用。

Ex3：多维高斯噪声与未知协方差矩阵

这个例子扩展到了多维输出的情况，并且噪声的协方差矩阵是未知的。具体设定是：在同一个过程上同时观测多个输出变量。观测模型为

这里 \mathbf{y}、\mathbf{y}_m 和 \mathbf{b}_i 都是大小为 n_y 的向量，即每个时刻观测 n_y 个变量。噪声 \mathbf{b}_i 是独立的零均值随机向量，但其协方差矩阵 \Sigma 是未知的。这与之前只考虑标量输出且方差已知的情况有本质区别。

单个观测向量 \mathbf{b}_i 的概率密度函数为

这是 n_y 维多元正态分布的标准形式。指数项中的二次型 \mathbf{b}_i^T\Sigma^{-1}\mathbf{b}_i 是马氏距离的平方，它衡量了向量 \mathbf{b}_i 相对于协方差结构 \Sigma 的"大小"。

整个观测序列的似然函数

由于各个时刻的噪声 \mathbf{b}_i 相互独立，整个观测序列 \mathbf{y}^S（包含所有 n_t 个时刻的观测）的似然函数是各时刻概率密度的乘积

其中残差向量定义为

这个似然函数现在依赖于两组未知量：模型参数 \mathbf{p} 和协方差矩阵 \Sigma。在标量情况下，如果方差已知，我们只需要估计 \mathbf{p}；但现在 \Sigma 本身包含 n_y(n_y+1)/2 个独立参数（因为协方差矩阵是对称的），也需要从数据中估计。

对数似然函数的展开

为了便于优化，考虑对数似然函数。将上述似然函数取对数并整理，得到

这个表达式包含三项。第一项是常数项，不影响优化。第二项涉及 \Sigma 的行列式，它体现了协方差矩阵的"体积"。第三项是所有残差的加权平方和，权重由 \Sigma^{-1} 给出。

由于 \Sigma 未知，我们需要同时求解 \hat{\mathbf{p}}_{\text{mv}} 和 \hat{\Sigma}_{\text{ML}}。这通过对对数似然分别关于 \Sigma 和 \mathbf{p} 求偏导并令其为零来实现

这是一个联合优化问题，两个方程需要同时满足。

矩阵求导规则

为了求解上述方程，需要用到矩阵微积分的一些基本规则。这里给出两个关键的求导公式。

第一个是关于行列式对数的导数：对于正定矩阵 \Sigma，有

这个结果可以从行列式的微分公式推导出来：d(\ln\det\Sigma) = \text{tr}(\Sigma^{-1}d\Sigma)，因此梯度是 \Sigma^{-1}。

第二个是关于二次型的导数：对于矩阵函数 f(\Sigma) = \mathbf{A}\Sigma^{-1}\mathbf{B}（其中 \mathbf{A} 和 \mathbf{B} 不依赖于 \Sigma），有

这个公式反映了逆矩阵求导的链式法则。负号来自于逆矩阵的导数 d(\Sigma^{-1}) = -\Sigma^{-1}(d\Sigma)\Sigma^{-1}。

参数估计器的推导

利用上述求导规则，对 \Sigma 求偏导并令其为零。从对数似然的第二项得到 -\frac{n_t}{2}\Sigma^{-1}，从第三项得到 \frac{1}{2}\sum_{i=1}^{n_t}\Sigma^{-1}\mathbf{e}_S(t_i, \mathbf{p})\mathbf{e}_S(t_i, \mathbf{p})^T\Sigma^{-1}。令总和为零并消去公共因子 \Sigma^{-1}，得到协方差矩阵的估计

这是残差向量的样本协方差矩阵。注意这里使用了 \hat{\mathbf{p}}_{\text{mv}}（即参数的最大似然估计），说明协方差估计依赖于参数估计。

对 \mathbf{p} 求偏导的结果显示，\hat{\mathbf{p}}_{\text{ML}} 可以通过最小化以下加权代价函数得到

这个代价函数的几何意义是：它最小化残差向量张成的椭球体的"体积"（由行列式衡量）。与标量情况下的简单平方和不同，多维情况需要考虑各输出变量之间的相关性结构，因此使用矩阵的行列式作为优化目标。实际计算中，这个问题通常通过迭代算法求解：先固定 \Sigma 估计 \mathbf{p}，再固定 \mathbf{p} 更新 \Sigma，如此交替进行直到收敛。

优化（Optimization）课程笔记

课程概述

本课程主要讲解优化理论的基础内容，包括无约束最小化问题、最优性条件、二次代价函数的特殊情况以及下降算法。课程的核心目标是建立从理论定义到实际算法的完整框架。

第一部分：无约束最小化（Unconstrained Minimization）

问题设定

我们考虑一个定义在 \mathcal{D} \subset \mathbb{R}^n 上的代价函数 f(\mathbf{x})。这里的定义域 \mathcal{D} 是 n 维实数空间的一个子集，代价函数将这个子集中的每个点映射到一个实数值。

严格局部最小化器

点 \mathbf{x}^* 被称为函数 f(\mathbf{x}) 的严格局部最小化器，需要满足：存在某个正数 \delta > 0，使得在 \mathbf{x}^* 的邻域内，函数行为满足特定条件。

具体来说，首先函数 f(\mathbf{x}) 必须在邻域

上有定义。这个邻域包含了所有与 \mathbf{x}^* 距离小于 \delta 的点。其次，对于这个邻域内除了 \mathbf{x}^* 本身之外的所有点 \mathbf{y}，都必须有 f(\mathbf{x}^*) < f(\mathbf{y})。这意味着 \mathbf{x}^* 的函数值严格小于其邻域内任何其他点的函数值，这就是"严格"的含义所在。

局部最小化器

相比严格局部最小化器，局部最小化器的定义放宽了条件。点 \mathbf{x}^* 是函数 f(\mathbf{x}) 的局部最小化器，同样要求存在 \delta > 0 使得函数在邻域 N(\mathbf{x}^*, \delta) 上有定义，但此时只需要 f(\mathbf{x}^*) \leq f(\mathbf{y}) 对所有 \mathbf{y} \in N(\mathbf{x}^*, \delta)、\mathbf{y} \neq \mathbf{x}^* 成立即可。注意这里用的是小于等于而非严格小于，这表示邻域内可能存在其他点与 \mathbf{x}^* 的函数值相同。同时，定义明确指出 \mathbf{x}^* 不是严格局部最小化器。

无论是哪种情况，f(\mathbf{x}^*) 都被称为 f(\mathbf{x}) 的局部最小值。

不同类型最小值的几何理解

通过函数图像可以直观理解这些概念。一个函数可能有多个局部最小值点（标注为"Minima locaux"），这些点在各自的局部邻域内是最小的。其中某些点可能是严格局部最小值（标注为"Minimum local strict"），而整个定义域上函数值最小的点称为全局最小值（标注为"Minimum global"）。这个几何图景说明了局部性质与全局性质的区别。

从定义到算法的过渡

前面给出的严格局部最小化器和局部最小化器的定义，虽然在数学上是严格的，但它们本身并不能直接告诉我们如何去寻找这些最小化器。这些定义只是描述了最小化器应该满足什么性质，却没有提供构造性的方法来定位它们。

为了发展出实际可用的搜索算法，我们需要更强的条件。当代价函数具有足够的光滑性时，我们可以推导出构造性的最优性条件。这些条件将利用函数的导数信息，从而为算法设计提供理论基础。

第二部分：最优性条件（Optimality Conditions）

一维情况的基本定理

我们要解决的问题是：在实数域 \mathbb{R} 上最小化函数 f(x)。

当函数满足适当的光滑性条件时，可以建立一个基础性定理。假设函数 f 属于 C^2 类（即二阶连续可微），并且其二阶导数 f'' 满足Lipschitz条件。Lipschitz条件保证了二阶导数的变化是受控的，不会出现突变。在这些前提下，如果 x^* 是函数 f(x) 的局部最小化器，那么必然有 f'(x^*) = 0。

这个结论的直观理解是：在局部最小值点处，函数的切线必须是水平的。如果导数不为零，那么沿着导数的反方向移动一小步，函数值会下降，这与 x^* 是局部最小化器矛盾。因此，一阶导数为零是局部最小化器的必要条件。

定理的反证法证明

第一个必要条件的证明

前面提到的定理（如果 x^* 是局部最小化器且 f \in C^2 满足 f'' Lipschitz，则 f'(x^*) = 0）可以通过反证法来证明。我们假设结论不成立，即假设 f'(x^*) \neq 0，然后推导出矛盾。

在点 x^* + h 处对函数 f 进行Taylor展开，存在某个 \theta \in [0,1] 使得

这个展开式表达了函数在 x^* 附近的行为。由于 f'' 满足Lipschitz条件，二阶导数 f''(x^* + h\theta) 在邻域 N(x^*, \delta) 内被某个常数 M 所界定。现在选择特殊的步长 h = -\alpha f'(x^*)，其中 \alpha 是一个待定的正数。将这个选择代入Taylor展开式，得到

观察这个不等式右边的结构：第二项 -\alpha(f'(x^*))^2 是负的（因为导数不为零），而第三项 \frac{\alpha^2 M^2}{2}(f'(x^*))^2 是正的但包含 \alpha^2。当 \alpha 足够小时，线性的负项会主导二次的正项，因此存在充分小的 \alpha 使得

这意味着我们找到了一个点 x^* - \alpha f'(x^*) 在 x^* 的邻域内，其函数值严格小于 f(x^*)，这与 x^* 是局部最小化器矛盾。因此原假设 f'(x^*) \neq 0 不能成立。

二阶条件的证明

有了一阶必要条件 f'(x^*) = 0 之后，我们需要进一步探讨二阶导数的性质。这次假设函数 f \in C^3，并且三阶导数 f^{(3)} 满足Lipschitz条件。假设 x^* 是局部最小化器。

第一个定理断言：在上述条件下，不仅有 f'(x^*) = 0，还必须有 f''(x^*) \geqslant 0。证明同样采用反证法，假设 f''(x^*) < 0。在 x^* + h 处进行Taylor展开，由于已知 f'(x^*) = 0，展开式变为

其中 \theta \in [0,1]。由于 f^{(3)} 满足Lipschitz条件，三阶导数项在邻域 N(x^*, \delta) 内被常数 M 界定。选择 h = \alpha f''(x^*)，代入后得到

当 \alpha 充分小时，右边的 \frac{\alpha^2}{2}(f''(x^*))^3 + \frac{\alpha^3 M^3}{6}(f''(x^*))^3 < 0（因为假设 f''(x^*) < 0，所以三次方也是负的），这导致

这再次与 x^* 是局部最小化器矛盾，因此 f''(x^*) \geqslant 0 必须成立。

充分条件：从局部到严格局部

前面讨论的都是必要条件，即"如果是最小化器，那么必须满足什么"。现在考虑一个反向的问题：如果条件更强，能否保证 x^* 是严格局部最小化器？

定理指出：假设 f \in C^2，f^{(2)} 满足Lipschitz条件，如果 x^* 是局部最小化器，并且 f'(x^*) = 0 且 f''(x^*) > 0（注意这里是严格大于零），那么 x^* 必然是严格局部最小化器。

证明的思路是直接验证严格局部最小化器的定义。考虑 x^* 的一个邻域 N(x^*, \delta)，在这个邻域内 f''(x^*) > 0。对于任何 x^* + h \in N(x^*, \delta)，Taylor展开给出

这里用到了 f'(x^*) = 0 以及二阶导数在邻域内保持正性。这个严格不等式对所有 h \neq 0 成立，因此 x^* 是严格局部最小化器。

多维情况的推广

从一维到多维

现在要解决的问题变成了：在 \mathbb{R}^n 空间中最小化函数 f(\mathbf{x})。一维情况的所有结论都可以推广到多维，但符号表示需要相应调整。

定理陈述如下：假设 f \in C^2，\mathbf{x}^* 是函数 f(\mathbf{x}) 的局部最小化器（或者更强地，是严格局部最小化器），那么必须满足梯度条件

其中 \mathbf{g}(\mathbf{x}) 是 f 在点 \mathbf{x} 的梯度向量，它是所有偏导数组成的向量。梯度为零意味着函数在 \mathbf{x}^* 处沿任何方向的方向导数都为零，这是一维情况下 f'(x^*) = 0 的自然推广。

除了一阶条件，还有二阶条件。\mathbf{H}(\mathbf{x}^*) 表示 f 在 \mathbf{x}^* 处的Hessian矩阵（即所有二阶偏导数组成的矩阵），它必须是半正定的，即对于任意向量 \mathbf{s} \in \mathbb{R}^n 都有 \mathbf{s}^T \mathbf{H}(\mathbf{x}^*) \mathbf{s} \geqslant 0。如果 \mathbf{x}^* 是严格局部最小化器，那么Hessian矩阵是正定的，即 \mathbf{s}^T \mathbf{H}(\mathbf{x}^*) \mathbf{s} > 0 对所有非零向量 \mathbf{s} 成立。

这些多维定理的证明思路与一维情况完全类似，都是通过Taylor展开和反证法进行。Hessian矩阵在多维优化中扮演的角色，正如二阶导数在一维中的角色一样，刻画了函数的局部曲率。

第三部分：二次代价函数（Quadratic Cost Functions）

二次函数的重要性

二次代价函数构成了优化理论中的一个特殊但极其重要的类别。这类函数之所以值得单独讨论，有两个主要原因。首先，许多实际问题本身就直接涉及二次函数的最小化，例如最小二乘问题和某些控制问题。其次，即使原问题不是二次的，二次近似也常常被用作优化算法的基础。许多数值算法的核心思想就是在每一步迭代中用二次函数近似原函数，然后最小化这个二次近似，这样可以利用二次函数的良好性质。

二次代价函数的实例：数字通信系统中的信道估计

考虑一个数字通信系统场景。发射端发送一个符号序列 \{x_n\}，其中每个符号 x_n = \pm 1，即只能取两个值。这些符号通过信道传输后，接收端得到的序列为 \{y_n\}，这里 y_n \in \mathbb{R} 是实数值。信道的作用通过输入输出关系来建模，具体表现为

这个公式描述了线性时不变信道的卷积模型。接收到的每个输出符号 y_n 是由当前及过去 L-1 个输入符号的线性组合形成的，权重系数 h_k（k = 0, \ldots, L-1）就是信道参数。这些参数完全刻画了信道的特性，但它们是未知的。

为了估计这些未知的信道参数，需要采用训练过程。具体做法是：发送一个已知的训练序列 \{x_n\}_{n=1}^N，这个序列在发射端和接收端都是已知的。接收端收到对应的输出 \{y_n\} 后，就可以利用输入输出关系来反推信道参数。

从卷积关系到矩阵方程

根据信道的输入输出关系，可以写出一系列方程。从时刻 L 开始（之前的时刻没有足够的历史输入），有

这些方程可以用矩阵形式紧凑地表示为 \mathbf{y} = \mathbf{X}\mathbf{h}。这里输出向量为 \mathbf{y} = (y_L, y_{L+1}, \ldots, y_N)^T，未知参数向量为 \mathbf{h} = (h_0, h_1, \ldots, h_{L-1})^T。输入矩阵 \mathbf{X} 是一个 (N-L+1) \times L 的矩阵，其结构为

每一行对应一个时刻的方程，包含了产生该时刻输出所需的 L 个历史输入。

最小二乘估计与二次代价函数

一般情况下，矩阵 \mathbf{X} 不是方阵，方程组 \mathbf{y} = \mathbf{X}\mathbf{h} 可能没有精确解（特别是当存在噪声时）。因此需要寻找最小二乘意义下的最优估计 \hat{\mathbf{h}}_{LS}，即最小化残差的二范数平方

将这个目标函数展开。首先注意到

展开这个乘积得到

利用对称性 \mathbf{h}^T \mathbf{X}^T \mathbf{y} = \mathbf{y}^T \mathbf{X}\mathbf{h}（两者都是标量且互为转置），可以简化为

因此代价函数表示为

这正是二次函数的标准形式。最后一项 \mathbf{y}^T \mathbf{y} 与优化变量 \mathbf{h} 无关，不影响最优解的位置。这个例子展示了实际问题如何自然地导致二次代价函数。

一般二次代价函数的结构

现在考虑一般形式的二次代价函数

这里 \alpha 是标量常数，\mathbf{c} 是 n 维向量，\mathbf{H} 是 n \times n 对称矩阵。对称性是自然的假设，因为任何矩阵都可以分解为对称和反对称部分，而在二次型 \mathbf{x}^T \mathbf{H}\mathbf{x} 中，反对称部分的贡献为零。

计算这个函数的梯度。对于二次项，利用 \nabla(\mathbf{x}^T \mathbf{H}\mathbf{x}) = (\mathbf{H} + \mathbf{H}^T)\mathbf{x}，当 \mathbf{H} 对称时简化为 2\mathbf{H}\mathbf{x}。对于线性项，\nabla(\mathbf{c}^T \mathbf{x}) = \mathbf{c}。因此梯度为

注意这里系数 \frac{1}{2} 的作用：它使得梯度公式更简洁。二阶导数即Hessian矩阵为

这说明二次函数的Hessian矩阵是常数，不依赖于位置 \mathbf{x}。根据前面建立的最优性条件，函数 \Phi(\mathbf{x}) 有最小化器 \mathbf{x}^* 的充要条件是梯度为零，即

这是一个线性方程组。如果 \mathbf{H} 可逆，最小化器唯一且为 \mathbf{x}^* = -\mathbf{H}^{-1}\mathbf{c}。

从最小化器出发的函数值分析

假设 \mathbf{x}^* 满足 \mathbf{H}\mathbf{x}^* = -\mathbf{c}，考虑任意其他点 \mathbf{x}^* + a\mathbf{p}，其中 a 是标量，\mathbf{p} 是任意方向向量。将这个点代入代价函数

展开得到

利用条件 \mathbf{H}\mathbf{x}^* = -\mathbf{c}，有 \mathbf{x}^{*T}\mathbf{H}\mathbf{p} = -\mathbf{c}^T\mathbf{p}，这两项相互抵消。因此

这个简洁的公式表明：从 \mathbf{x}^* 出发沿任何方向 \mathbf{p} 移动距离 a，函数值的变化由二次项 \frac{1}{2}a^2\mathbf{p}^T\mathbf{H}\mathbf{p} 决定。

特征值对函数行为的刻画

为了理解这个二次项的符号，需要考察 \mathbf{H} 的特征值。这里简要回顾特征值的定义：\lambda 是矩阵 \mathbf{H} 的特征值，如果存在非零向量 \mathbf{u} 满足 \|\mathbf{u}\|^2 = 1 和 \mathbf{H}\mathbf{u} = \lambda\mathbf{u}。一般的 n \times n 矩阵有 n 个特征值（可能是复数，可能重复）及对应的特征向量。但对于对称矩阵，有重要的特殊性质：所有特征值都是实数，并且存在 n 个对应的特征向量。

设 \lambda_i 和 \mathbf{u}_i（i = 1, \ldots, n）分别表示 \mathbf{H} 的特征值和对应的标准正交特征向量。如果选择方向 \mathbf{p} = \mathbf{u}_i，代入前面的公式得到

利用特征值定义 \mathbf{H}\mathbf{u}_i = \lambda_i\mathbf{u}_i 以及 \mathbf{u}_i^T\mathbf{u}_i = 1，得到

这个结果揭示了特征值的几何意义：当沿着特征向量 \mathbf{u}_i 的方向从 \mathbf{x}^* 移动时，函数值的变化率由对应的特征值 \lambda_i 控制。具体而言，如果 \lambda_i > 0，那么 \Phi(\mathbf{x}^* + a\mathbf{u}_i) 会增加，说明 \mathbf{x}^* 在这个方向上是局部最小；如果 \lambda_i = 0，函数值保持不变，\mathbf{x}^* 在这个方向上既不是最小也不是最大；如果 \lambda_i < 0，函数值会减小，说明可以通过沿这个方向移动找到更小的函数值，因此 \mathbf{x}^* 不是最小化器。

二次代价函数的几何形态

PPT展示了三个二维二次函数的例子，通过三维曲面图和等高线图来直观理解它们的几何特性。

左图展示了一个具有唯一最小值点的情况，曲面呈现出一个明确的"碗"状，等高线是椭圆形的闭合曲线，中心点就是全局最小化器。中间图的情况类似，但等高线的椭圆更为狭长，说明在不同方向上函数的曲率差异较大。右图则展示了一个鞍点的情况，曲面既不完全向上弯曲也不完全向下弯曲，等高线呈现出双曲线形状。从三维看，某些方向上函数值上升（蓝色到黄色），而在另一些方向上函数值下降，这对应于Hessian矩阵既有正特征值又有负特征值的情形。

二次代价函数最小化器的完整刻画

对于一般形式的二次代价函数

其最小化器的存在性和唯一性完全由Hessian矩阵 \mathbf{H} 的性质决定。第一种情况：如果线性方程组 \mathbf{H}\mathbf{x} = -\mathbf{c} 有唯一解 \mathbf{x}^*，并且 \mathbf{H} 是正定的（即所有特征值都严格大于零），那么 \mathbf{x}^* 是唯一的全局最小化器（孤立的全局最小化器）。正定性保证了从 \mathbf{x}^* 出发沿任何非零方向移动，函数值都会严格增加，这正是唯一全局最小值点的特征。

第二种情况：如果方程 \mathbf{H}\mathbf{x} = -\mathbf{c} 有一组解（即解集合不止一个点），并且 \mathbf{H} 是半正定的（即特征值为正或零），那么存在一个全局最小化器的集合。半正定意味着某些方向上特征值为零，函数在这些方向上是平坦的，因此沿着这些方向的所有点都具有相同的最小函数值。这种情况下优化问题虽然有解，但解不唯一。

二次代价函数的最小化算法步骤

对于二次代价函数的最小化，有明确的算法步骤。第一步是求解线性方程组 \mathbf{H}\mathbf{x} = -\mathbf{c}。这个方程的解就是候选的最小化器。第二步是验证矩阵 \mathbf{H} 是否为半正定。只有当 \mathbf{H} 至少是半正定时，第一步得到的解才真正是最小化器。

关于第一步求解线性方程组，PPT提供了两个技术性备注。第一个备注指出，方程 \mathbf{H}\mathbf{x} = -\mathbf{c} 的解可以通过对 \mathbf{H} 进行QR分解来获得。QR分解是一种将矩阵分解为正交矩阵和上三角矩阵乘积的方法，这种分解使得线性方程组的求解在数值上更加稳定。或者，当 \mathbf{H} 可以表示为 \mathbf{H} = \mathbf{A}^T\mathbf{A} 时（这种情况在最小二乘问题中很常见），可以通过对矩阵 \mathbf{A} 进行QR分解来间接求解。第二个备注针对大规模问题：当变量 \mathbf{x} 的维度非常大时（例如在图像处理应用中），直接求解线性方程组可能计算量过大，此时可以使用迭代下降算法作为替代。

第四部分：下降算法（Descent Algorithms）

下降算法的基本框架

现在考虑更一般的优化问题：在 \mathbb{R}^n 空间中最小化函数 f(\mathbf{x})。对于这个问题，通常假设函数 f 属于 C^1 类（一阶连续可微），有时也假设属于 C^2 类（二阶连续可微）。与二次函数不同，一般函数可能没有显式解，需要通过迭代算法来逼近最优解。

下降算法的核心思想是通过迭代方式搜索函数 f(\mathbf{x}) 的局部最小化器 \mathbf{x}^*。算法从某个初始点 \mathbf{x}_0 开始，构造一个序列 \{\mathbf{x}_k\}（k = 1, 2, \ldots），使得这个序列满足三个基本性质。第一，函数值严格下降，即 f(\mathbf{x}_{k+1}) < f(\mathbf{x}_k)，这是"下降"的直接含义。第二，梯度序列 \{\nabla f(\mathbf{x}_k)\} 收敛到零向量。根据前面建立的最优性条件，梯度为零是局部最小化器的必要条件，因此这个性质保证了序列向着满足必要条件的点收敛。

在后续讨论中，为了简化符号，记 \nabla f(\mathbf{x}_k) = \mathbf{g}_k 表示在第 k 次迭代点的梯度，\nabla^2 f(\mathbf{x}_k) = \mathbf{H}_k 表示在该点的Hessian矩阵。

下降算法的迭代机制

右侧的图示直观展示了下降算法的工作方式。左边的三维曲面图显示了函数 f(\mathbf{x}) 的形状，右边的等高线图显示了算法迭代路径在二维平面上的投影。从初始点 \mathbf{x}_0 开始，算法在每次迭代中都会产生一个新的点，这些点连成的轨迹逐步向着函数的最小值点移动，等高线的密集程度反映了函数梯度的大小。

每次迭代包含两个关键步骤。第一步是寻找下降方向 \mathbf{p}_k：从当前点 \mathbf{x}_k 出发，需要找到一个方向 \mathbf{p}_k 使得沿这个方向移动会使函数值下降。数学上的判断准则是方向导数必须为负，即 \mathbf{p}_k^T\mathbf{g}_k < 0（这里假设梯度 \mathbf{g}_k \neq \mathbf{0}，因为如果梯度已经为零就没必要继续迭代了）。这个不等式的几何意义是：方向 \mathbf{p}_k 与梯度 \mathbf{g}_k 的夹角必须大于90度，即两者大致指向相反的方向。

为什么这个条件能保证函数下降？这可以通过Taylor定理来理解。在点 \mathbf{x}_k 附近，对于足够小的步长，函数可以近似为

如果 \mathbf{p}_k^T\mathbf{g}_k < 0，那么当 \alpha 充分小时，右边的线性项是负的，从而 f(\mathbf{x}_k + \alpha\mathbf{p}_k) < f(\mathbf{x}_k)，这保证了函数值的下降。

第二步是寻找步长 \alpha_k：确定了下降方向后，需要决定沿这个方向移动多远。步长 \alpha_k 必须选择得使得 f(\mathbf{x}_k + \alpha_k\mathbf{p}_k) < f(\mathbf{x}_k)。步长的选择需要平衡两个因素：如果步长过小，虽然能保证下降但收敛速度会很慢；如果步长过大，可能会越过最小值点或者函数的线性近似失效导致函数值反而上升。确定步长后，下一个迭代点就是 \mathbf{x}_{k+1} = \mathbf{x}_k + \alpha_k\mathbf{p}_k。

通用下降算法的伪代码

PPT给出了下降算法的完整伪代码结构，这个框架是所有具体下降算法的通用模板。算法从初始化开始：设置迭代计数器 k := 0 和完成标志 finished := false。然后进入主循环，执行以下步骤。

首先调用 \mathbf{p}_k := \text{descent\_direction}(\mathbf{x}_k) 来计算当前点 \mathbf{x}_k 的下降方向。这个函数是算法的核心部分之一，不同的下降算法主要区别就在于如何选择这个方向。如果函数返回的下降方向不存在（即 \mathbf{p}_k does not exist），说明当前点已经满足某种终止条件（比如梯度已经足够接近零），此时设置 finished := true 标记算法完成。

如果下降方向存在，则继续第二步：通过 \alpha_k := \text{line\_search}(\mathbf{x}_k, \mathbf{p}_k) 来确定沿着方向 \mathbf{p}_k 应该移动多远。线搜索函数接收当前点和下降方向作为输入，返回最优的步长。得到步长后，更新迭代点 \mathbf{x}_{k+1} := \mathbf{x}_k + \alpha_k\mathbf{p}_k，并将迭代计数器增加 k := k + 1。

最后检查 \text{stopping\_condition}(finished, \mathbf{p}_k, \mathbf{x}_k, \mathbf{x}_{k+1}) 是否为假。这个停止条件函数综合考虑多个因素来判断是否应该终止迭代，包括是否已经标记完成、当前的下降方向、连续两次迭代点的位置等。只要停止条件返回假，循环就继续执行。

步长选择不当的后果

为了理解步长选择的重要性，PPT展示了一个简单但富有启发性的例子。考虑代价函数 f(x) = x^2，这是最简单的凸二次函数，其最小值显然在 x = 0 处。选择下降方向为 p_k = (-1)^k，即方向在每次迭代中左右交替。步长公式为 \alpha_k = 2 + 3/2^{k+1}，初始点设为 x_0 = 2。

第一个图展示了步长过大的情况。从 x_0 = 2 开始，第一步的方向是 p_0 = 1（因为 (-1)^0 = 1），实际上这个方向是错误的，它指向远离最小值的方向。但更严重的问题是步长：\alpha_0 = 2 + 3/2 = 3.5，这导致第一步就移动到 x_1 = 2 + 3.5 \times 1 = 5.5，不但没有靠近最小值点0，反而离得更远了。后续的迭代虽然方向会改变，但由于步长设计不合理，迭代点在抛物线两侧大幅震荡，收敛极慢甚至发散。图中的虚线箭头显示了这种跳跃式的、不稳定的迭代轨迹。

第二个图展示了相反的极端：步长过小。这次下降方向设为 p_k = -1（固定指向左侧），步长为 \alpha_k = 1/2^{k+1}。从 x_0 = 2 开始，每步移动的距离是 \alpha_0 = 1/2, \alpha_1 = 1/4, \alpha_2 = 1/8, \ldots，呈指数衰减。虽然每一步都确实在向最小值点靠近，但由于步长太小，进展非常缓慢。图中可以看到迭代点 x_1, x_2, x_3, x_4, x_5 密集地挤在一起，虽然理论上会收敛到0，但需要无穷多步才能到达，实际计算中会因为收敛速度过慢而效率低下。

这两个例子说明：步长既不能太大（否则会震荡或发散），也不能太小（否则收敛过慢）。因此需要一个系统的方法来选择合适的步长，这就是线搜索要解决的问题。

线搜索：一维优化子问题

当在第 k 次迭代中确定了起始点 \mathbf{x}_k 和下降方向 \mathbf{p}_k 后，线搜索本质上变成了一个一维优化问题。我们要寻找的是最优步长

这个公式的含义是：把 \alpha 作为唯一变量，在保持 \mathbf{x}_k 和 \mathbf{p}_k 固定的情况下，寻找使得函数值下降最多的步长。注意这里是在优化一个单变量函数，因此称为"一维"优化。原本 n 维空间中的优化问题，在确定了方向后，就简化为沿着这条直线的优化问题。

对于这个一维优化，原则上可以采用两种策略。第一种是精确线搜索：求解上述最小化问题得到真正的最优 \alpha_k。然而PPT明确指出，精确线搜索通常不被采用，原因是计算效率比不佳。精确求解一维优化本身需要多次函数求值或导数计算，这些额外开销可能抵消掉找到更优步长带来的收益。特别是当函数 f 本身计算代价很高时，每次都精确求最优步长是不经济的。

第二种是近似线搜索：不寻求绝对最优，而是找一个"足够好"的步长。这是实际算法中通常采用的策略。近似的标准需要在两个目标间平衡：步长要足够大以保证每步有实质性进展，但又不能大到违反下降性质或使近似失效。

回溯算法：一种实用的近似线搜索

回溯算法（Backtracking algorithm）是实现近似线搜索的一种简单而有效的方法。算法的初始化阶段选择一个初始步长 \alpha_{\text{init}} > 0（例如可以简单地设为1），并设置 \alpha^{(0)} := \alpha_{\text{init}}，内部迭代计数器 \ell := 0。

然后进入while循环，检查条件 f(\mathbf{x}_k + \alpha^{(\ell)}\mathbf{p}_k) \geqslant f(\mathbf{x}_k)。这个条件判断当前步长 \alpha^{(\ell)} 是否导致函数值上升或不变。如果不等式成立，说明步长太大了，需要缩减。缩减的方式是 \alpha^{(\ell+1)} := \tau\alpha^{(\ell)}，其中 \tau \in (0,1) 是收缩因子（例如 \tau = 0.5），然后增加计数器 \ell := \ell + 1。这个过程不断重复，每次将步长乘以 \tau 使其变小，直到找到一个能使函数值下降的步长为止。最终接受的步长是 \alpha_k := \alpha^{(\ell)}。

PPT指出了这个回溯方法的一个特点：它能避免步长过小（因为从一个较大的初始值开始，只在必要时才缩减），但无法避免步长过大的情况（如果初始步长或收缩不够，最终接受的步长仍可能过大）。这就需要更精细的准则。

Armijo条件：充分下降准则

Armijo条件提供了一个更精细的判断准则，不仅要求函数值下降，还要求下降的"充分性"。条件表述为

这里 \beta \in (0,1)（典型值如 \beta = 0.1 或 0.001）。为了理解这个条件的含义，先看右边的项 \alpha_k\beta\mathbf{p}_k^T\mathbf{g}_k。回忆 \mathbf{p}_k^T\mathbf{g}_k < 0（下降方向的定义），因此这一项是负的，它代表了函数在点 \mathbf{x}_k 处沿方向 \mathbf{p}_k 的线性近似下降量的 \beta 倍。

Armijo条件要求：实际函数值的下降（左边）必须至少达到线性预测下降的 \beta 倍（右边）。由于 \beta 通常取很小的值（如0.1），这个条件其实是相当宽松的——只要求实际下降量达到理论预测的10%即可。这种"充分下降"的要求既保证了函数值确实在减小，又不至于要求过严使得很难找到满足条件的步长。结合前面的回溯框架，可以将while循环的条件改为检查Armijo条件是否被违反，从而得到一个既避免步长过小、又确保充分下降的实用算法。

Armijo条件的几何解释

图中展示了Armijo条件的几何意义。横轴表示步长 \alpha，纵轴表示沿下降方向的函数值。蓝色实线曲线 f(\mathbf{x}_k + \alpha\mathbf{p}) 表示从点 \mathbf{x}_k 出发沿方向 \mathbf{p} 移动距离 \alpha 后的实际函数值。这条曲线在 \alpha = 0 处的切线斜率为 \mathbf{p}^T\mathbf{g}_k（一阶导数），由于 \mathbf{p} 是下降方向，这个斜率是负的，因此曲线在原点附近向下倾斜。

红色虚线 f(\mathbf{x}_k) + \alpha\beta\mathbf{p}^T\mathbf{g}_k 表示Armijo条件的右侧界限。这是一条过点 (0, f(\mathbf{x}_k)) 的直线，其斜率为 \beta\mathbf{p}^T\mathbf{g}_k。由于 \beta \in (0,1)，这条直线的斜率比实际函数在原点处的切线斜率更平缓（绝对值更小）。Armijo条件要求蓝色曲线必须位于红色直线下方或与其相交。

上方的红色虚线 f(\mathbf{x}_k) + \alpha\mathbf{p}^T\mathbf{g}_k 代表一阶Taylor近似（\beta = 1 的情况）。从图中可以看出，在较小的 \alpha 值范围内，实际函数曲线确实位于Armijo界限下方，说明这些步长都满足充分下降条件。图的说明指出：所有 \alpha < 1.8 的值都是可能的，即在这个范围内的步长都满足Armijo条件。这个范围的存在为线搜索算法提供了灵活性。

PPT还提到存在许多其他线搜索方法，Armijo条件只是其中一种常用的准则。

满足Armijo条件的步长存在性定理

一个自然的问题是：是否总能找到满足Armijo条件的步长 \alpha_k？定理回答了这个问题。假设函数 f \in C^1，梯度 \mathbf{g}_k 在点 \mathbf{x}_k 处满足Lipschitz条件（Lipschitz常数为 \gamma(\mathbf{x}_k)），\beta \in (0,1)，并且 \mathbf{p}_k 是 \mathbf{x}_k 处的下降方向。在这些条件下，Armijo条件

对所有 \alpha \in [0, \alpha_{\max}(\mathbf{x}_k, \mathbf{p}_k)] 成立，其中

这个公式给出了满足Armijo条件的最大步长。注意到分子中 (\beta - 1)\mathbf{p}_k^T\mathbf{g}_k 是正的（因为 \beta < 1 使得 \beta - 1 < 0，而 \mathbf{p}_k^T\mathbf{g}_k < 0，负负得正），因此 \alpha_{\max} 确实是正数。这个结果保证了满足条件的步长区间非空，线搜索总能找到可行的步长。

更进一步，定理指出回溯-Armijo算法（结合了回溯策略和Armijo条件）产生的步长 \alpha_k 满足

这个下界保证了算法不会选择过小的步长。如果初始步长 \alpha_{\text{init}} 足够小，或者梯度和方向的几何关系使得右边第二项较大，那么算法会接受一个合理大小的步长。这个结果说明回溯-Armijo算法既避免步长过小（有下界保护），又保证充分下降（通过Armijo条件）。

通用下降算法的收敛性

现在考察通用下降算法的收敛行为。假设函数 f \in C^1，梯度 \mathbf{g}_k 在整个 \mathbb{R}^n 上满足Lipschitz条件。定理指出：通用下降算法的连续迭代序列必然满足以下三种情况之一。

第一种情况：存在某个有限的迭代次数 \ell > 0 使得 \mathbf{g}_\ell = \mathbf{0}。这意味着算法在有限步内找到了一个驻点（梯度为零的点）。这种情况虽然理论上可能，但在实际计算中由于数值误差很少恰好发生。

第二种情况：当 k \to +\infty 时，f_k = f(\mathbf{x}_k) \to -\infty。这表示函数值无界下降。这种情况发生在函数没有下界时，算法会持续找到更低的函数值点。虽然这在数学上是一种可能性，但在实际问题中通常意味着问题设定有误或需要添加约束。

第三种情况：

这个极限条件比较复杂，它涉及两个量的最小值。第一个量 |\mathbf{p}_k^T\mathbf{g}_k| 是下降方向与梯度的内积的绝对值，反映了沿当前方向的方向导数大小。第二个量 \frac{|\mathbf{p}_k^T\mathbf{g}_k|}{\|\mathbf{p}_k\|^2} 是对方向长度归一化后的结果。

PPT特别指出：前两种情况是令人满意的终止状态（第一种找到了精确驻点，第二种说明函数无下界），但第三种情况不那么理想。第三种情况可能发生在下降方向逐渐变得与梯度正交的时候。如果 \mathbf{p}_k^T\mathbf{g}_k \to 0 而 \mathbf{g}_k 本身不趋于零，这意味着算法选择的方向越来越不是好的下降方向。这种现象在某些特殊情况下会出现，因此需要对下降方向的选择方式进行更仔细的设计。

停止条件的实际考量

由于算法在理论上可能需要无限次迭代才能严格收敛，实际应用中必须设定停止条件以在有限时间内给出解。算法通常在以下情况之一发生时停止。

第一个标准是达到最大迭代次数：k > k_{\max}。这是一个简单的保护性措施，防止算法运行时间过长。当计算资源有限或需要实时响应时，这个条件确保算法在预定时间内终止。

第二个标准是梯度变得足够小：\|\mathbf{g}_k\| < \varepsilon_{\mathbf{g}}。由于梯度为零是最优性的必要条件，当梯度范数小于某个阈值 \varepsilon_{\mathbf{g}} 时，可以认为已经足够接近驻点。这是最常用和最直观的停止准则。

第三个标准基于函数值的变化：代价函数下降不够显著。有两种具体形式。绝对变化准则：|f(\mathbf{x}_{k+1}) - f(\mathbf{x}_k)| < \varepsilon_f，即连续两次迭代的函数值差异小于阈值。相对变化准则：|f(\mathbf{x}_{k+1}) - f(\mathbf{x}_k)| < \varepsilon_f \|f(\mathbf{x}_k)\|，考虑了函数值量级的影响，当函数值本身很大时要求更大的绝对下降量。

PPT特别警告：不要将 \varepsilon_{\mathbf{g}} 或 \varepsilon_f 设置得过小。虽然理论上更小的阈值意味着更高的精度，但实际计算中会遇到计算精度问题。由于浮点运算的舍入误差，当梯度或函数值变化已经接近机器精度时，继续迭代不仅无法提高解的质量，反而可能导致数值不稳定。因此阈值应该根据问题的数值特性和所需精度合理选择。

梯度下降法

梯度下降法是最基本也最直观的下降算法。它的核心思想是选择使代价函数局部下降最快的方向作为搜索方向。

为了确定这个"最快下降"的方向，考虑函数在点 \mathbf{x}_k 处的线性模型。通过一阶Taylor展开，函数在 \mathbf{x}_k + \mathbf{p} 附近可以近似为

这个线性模型忽略了高阶项，只保留了常数项和一阶导数项。在所有可能的方向 \mathbf{p} \in \mathbb{R}^n 中，梯度下降法选择使这个线性模型下降最快的方向。具体来说，求解优化问题

这里添加了一个约束：搜索方向的范数等于梯度的范数。这个约束是必要的，因为没有长度限制的话，最小化问题没有有限解（可以无限延长方向使线性项任意负）。通过这个约束，问题变成在固定长度的所有方向中找最优方向。

这个约束优化问题的解非常简洁：

即梯度下降方向就是负梯度方向。这个结果的直观理解是：在线性近似下，函数沿梯度方向增长最快，因此沿负梯度方向下降最快。这也解释了为什么这个方法叫做"梯度下降"——每一步都沿着负梯度方向前进。

梯度下降法的收敛性

对于采用梯度下降方向的算法，可以获得更强的收敛性结果。假设函数 f \in C^1，梯度 \mathbf{g}_k 在 \mathbb{R}^n 上满足Lipschitz条件。使用梯度下降方向的通用下降算法的连续迭代满足以下三种情况之一。

第一种情况：存在某个 \ell > 0 使得 \mathbf{g}_\ell = \mathbf{0}，即在有限步找到驻点。第二种情况：\lim_{k \to +\infty} f_k = -\infty，函数值无界下降。第三种情况：\lim_{k \to +\infty} \mathbf{g}_k = \mathbf{0}，梯度序列收敛到零向量。

对比前面通用下降算法的收敛定理，梯度下降法的第三种情况有本质改进。在通用情况下，第三种可能是下降方向与梯度逐渐正交，这并不保证梯度趋于零。但当使用梯度下降方向 \mathbf{p}_k = -\mathbf{g}_k 时，第三种情况直接保证了 \mathbf{g}_k \to \mathbf{0}。这是因为梯度方向始终是最陡下降方向，不会出现与梯度正交的病态情况。因此，梯度下降法从理论上保证了收敛到驻点（在函数有下界的情况下）。

不过PPT在最后加了一句现实的评论：理论上梯度下降法收敛，但实践中并非那么容易。这暗示了虽然理论保证存在，但收敛速度、数值稳定性、局部最小值陷阱等实际问题仍然存在，需要在具体应用中仔细处理。

梯度下降法的实例：Rosenbrock函数

为了直观理解梯度下降法的实际表现，考虑一个经典的测试函数——Rosenbrock函数

这个函数的最小值点在 (1, 1) 处，但其等高线形状呈现出狭长的弯曲山谷。图中的等高线密集程度反映了函数值的变化率：等高线密集的区域梯度大，稀疏的区域梯度小。红色曲线显示了从某个初始点出发的梯度下降迭代轨迹。

观察这条轨迹可以发现一个重要现象：虽然每一步都选择了当前位置的最陡下降方向（负梯度方向），但整体收敛过程非常缓慢。迭代路径呈现锯齿状，在山谷中来回震荡，缓慢地向最小值点逼近。这种行为的原因在于：在狭长山谷中，梯度方向主要指向谷底（垂直于山谷走向），而不是指向真正的最小值点（沿着山谷方向）。每次迭代都在试图快速到达谷底，但到达谷底后又发现偏离了最优点，需要重新调整方向。

PPT的结论指出：即使下降方向在局部是最优的（最陡下降），收敛速度仍然可能非常慢。这揭示了梯度下降法的一个根本局限性——它只利用一阶导数信息，无法感知函数的曲率结构，因此在某些函数形状下效率低下。

梯度下降对变量单位的依赖性

梯度下降法还有一个微妙但重要的问题：它的行为依赖于变量的单位选择。左图显示了一个接近圆形的等高线，中心在原点。在这种情况下，从外围任意点 \mathbf{x}_k 出发，负梯度方向几乎直接指向中心，因此梯度下降会沿着大致直线的路径快速收敛到最小值点。

右图展示了同一个函数，但等高线变成了狭长的椭圆。这并不是因为函数本身改变了，而可能是因为变量的尺度不同。例如，如果第一个变量以米为单位，第二个变量以千米为单位，即使它们在物理意义上等价，数值上的单位差异也会导致等高线拉伸。在这种情况下，从点 \mathbf{x}_k 出发的负梯度方向不再指向中心，而是大致垂直于椭圆的长轴。结果是梯度下降轨迹会像第一个Rosenbrock例子那样，在最优点附近来回震荡，收敛速度大幅下降。

这个现象说明：梯度下降法对变量的尺度非常敏感。在实际应用中，通常需要对变量进行归一化或标准化预处理，使得各个变量的数值范围相当，否则算法的效率会受到严重影响。这种对坐标系选择的敏感性是梯度下降法的一个固有缺陷。

牛顿法及其变体

基于正定矩阵的下降方向

梯度下降法的问题促使我们寻找更好的方向选择策略。考虑一个对称正定矩阵 \mathbf{B}_k，如果向量 \mathbf{p}_k 满足线性方程

那么 \mathbf{p}_k 是一个下降方向。为了验证这一点，计算方向导数

这里利用了 \mathbf{B}_k\mathbf{p}_k = -\mathbf{g}_k 的条件。由于 \mathbf{B}_k 是正定的，对任何非零向量 \mathbf{p}_k 都有 \mathbf{p}_k^T\mathbf{B}_k\mathbf{p}_k > 0。因此当梯度 \mathbf{g}_k \neq \mathbf{0} 时，\mathbf{p}_k^T\mathbf{g}_k < 0 成立，这正是下降方向的定义条件。

这个构造的关键在于正定矩阵 \mathbf{B}_k 的引入。不同的 \mathbf{B}_k 选择会导致不同的下降方向，从而产生不同的优化算法。梯度下降法实际上对应于 \mathbf{B}_k = \mathbf{I}（单位矩阵）的特殊情况。通过选择更好的 \mathbf{B}_k，我们可以克服梯度下降法对坐标系敏感的问题。

二次近似与牛顿方向

现在考虑函数在点 \mathbf{x}_k 处的二次近似。通过二阶Taylor展开，函数在 \mathbf{x}_k + \mathbf{p} 附近可以近似为

这里 \mathbf{B}_k 是某个对称矩阵。与线性近似相比，二次近似多了一个二次项 \frac{1}{2}\mathbf{p}^T\mathbf{B}_k\mathbf{p}，这个项捕捉了函数的曲率信息。现在的策略是：不再最小化线性模型，而是最小化二次模型，即求解

对二次函数求最小值，根据前面二次代价函数的理论，最优解满足梯度为零的条件。计算 m_k^Q 关于 \mathbf{p} 的梯度并令其为零，得到

这正是前面提到的线性方程。当矩阵 \mathbf{B}_k 选择为函数在 \mathbf{x}_k 处的Hessian矩阵 \mathbf{H}_k 时，即 \mathbf{B}_k = \mathbf{H}_k，方程变为

这就是牛顿下降方向的定义。牛顿方向的几何意义是：它是使得二阶Taylor近似最小的方向。与梯度方向（使一阶近似下降最快）相比，牛顿方向利用了函数的曲率信息（通过Hessian矩阵），因此能更准确地预测最优移动方向。

PPT特别强调：牛顿方向只有在 \mathbf{H}_k 正定时才是下降方向。如果Hessian矩阵不是正定的（比如在鞍点附近或接近最大值点时），牛顿方向可能不满足 \mathbf{p}_k^T\mathbf{g}_k < 0，甚至可能指向函数值上升的方向。这是牛顿法的一个重要限制，在实际应用中需要对Hessian矩阵进行修正以保证正定性。

牛顿法及其变体的收敛性

对于使用牛顿方向或其变体的通用下降算法，可以建立收敛性定理。假设函数 f \in C^1，梯度 \mathbf{g}_k 在 \mathbb{R}^n 上满足Lipschitz条件。算法使用牛顿下降方向 \mathbf{H}_k\mathbf{p}_k = -\mathbf{g}_k 或其变体 \mathbf{B}_k\mathbf{p}_k = -\mathbf{g}_k。在这些条件下，连续迭代满足以下三种情况之一。

第一种情况：存在某个 \ell > 0 使得 \mathbf{g}_\ell = \mathbf{0}，即有限步找到驻点。第二种情况：\lim_{k \to +\infty} f_k = -\infty，函数值无界下降。第三种情况：\lim_{k \to +\infty} \mathbf{g}_k = \mathbf{0}，梯度序列收敛到零。

但这个结果需要一个关键的附加条件：矩阵 \mathbf{B}_k 的特征值必须一致有界且充分远离零。具体来说，存在常数使得所有迭代中 \mathbf{B}_k 的最小特征值和最大特征值都在某个固定范围内，并且最小特征值不能太接近零。这个条件保证了矩阵 \mathbf{B}_k 既不会变得病态（条件数过大），也不会接近奇异。

在标准牛顿法中，\mathbf{B}_k = \mathbf{H}_k 是函数的真实Hessian矩阵，其特征值随着迭代点变化。如果函数的Hessian矩阵在迭代路径上保持良好的正定性（特征值始终在合理范围内），那么收敛性就能得到保证。但如果Hessian矩阵的条件数很大或接近奇异，收敛性可能无法保证，这时需要对算法进行修正。

牛顿法的收敛速度：二次收敛

最后一张图再次展示了Rosenbrock函数，但这次使用牛顿法进行优化。与梯度下降的锯齿状轨迹形成鲜明对比，牛顿法的迭代路径要平滑和直接得多。图像显示迭代点快速且稳定地向最优点收敛，没有来回震荡的现象。

PPT指出：当牛顿法收敛时，其收敛速度是二次的。这是一个非常强的性质。所谓二次收敛是指：设迭代序列 \{\mathbf{x}_k\} 收敛到最优点 \mathbf{x}^*，误差为 \mathbf{e}_k = \mathbf{x}_k - \mathbf{x}^*，那么存在常数 C 使得

这意味着每次迭代后，误差的平方会进一步减小。二次收敛的实际含义是：一旦迭代点接近最优点，收敛速度会变得极快。如果某次迭代的误差是 10^{-2}，下次迭代误差可能就变成 10^{-4}，再下次变成 10^{-8}，有效数字几乎每次迭代翻倍。这种快速收敛是牛顿法相对于梯度下降法的主要优势。

然而需要注意的是，二次收敛只在迭代点已经足够接近最优点、且Hessian矩阵保持正定的情况下才能保证。在迭代初期或远离最优点时，牛顿法的全局收敛性较弱，可能需要结合线搜索等技术来稳定算法。此外，每次迭代需要计算和求解涉及Hessian矩阵的线性方程组，计算代价远高于梯度下降法，这在高维问题中可能成为瓶颈。

修正牛顿法（Modified Newton's Method）

处理非正定Hessian矩阵的问题

标准牛顿法要求Hessian矩阵 \mathbf{H}_k 是正定的，才能保证牛顿方向 \mathbf{p}_k 满足 \mathbf{H}_k\mathbf{p}_k = -\mathbf{g}_k 是下降方向。然而在实际优化过程中，特别是在远离最优点的区域、接近鞍点的位置、或者当前点是局部最大值点附近时，Hessian矩阵可能不是正定的（存在零或负特征值），甚至可能是未定义的（不可逆）。此时直接使用方程 \mathbf{H}_k\mathbf{p}_k = -\mathbf{g}_k 求得的 \mathbf{p}_k 不一定是下降方向，可能导致算法失败或朝错误方向移动。

修正牛顿法的核心思想是：通过添加一个修正矩阵 \mathbf{M}_k 来"修复"Hessian矩阵，使得修正后的矩阵成为正定的。具体做法是选择搜索方向 \mathbf{p}_k 满足

其中修正矩阵 \mathbf{M}_k 的选择必须保证 \mathbf{H}_k + \mathbf{M}_k 是正定的。这样一来，即使原始Hessian矩阵存在问题，修正后的系统仍然能够给出有效的下降方向。修正矩阵的选择有多种策略，需要在保证正定性和保持原Hessian信息之间取得平衡。

修正策略一：基于谱分解的修正

第一种修正策略基于Hessian矩阵的谱分解（特征值分解）。假设 \mathbf{H}_k 可以进行特征值分解为

这里 \mathbf{Q}_k 是正交矩阵，其列向量是 \mathbf{H}_k 的特征向量，\mathbf{D}_k 是对角矩阵，对角元素是对应的特征值。如果 \mathbf{H}_k 不是正定的，那么 \mathbf{D}_k 中会有零或负的对角元素。修正的方法是将这些问题特征值替换为正值，具体选择修正矩阵 \mathbf{M}_k 使得

这里 \varepsilon 是一个小的正数（如 10^{-6}），\max(\varepsilon\mathbf{I}, |\mathbf{D}_k|) 表示对 \mathbf{D}_k 的每个对角元素取绝对值，然后与 \varepsilon 比较取较大者。这个操作确保了所有特征值都是正的且有一个下界，从而保证修正后的矩阵是正定的。这种方法的优点是精确且理论上完美，但PPT指出分解的计算代价很高。对于 n \times n 矩阵，特征值分解需要 O(n^3) 的运算量，这在高维问题中可能难以承受。

修正策略二：基于最小特征值估计的修正

第二种策略更加实用且计算效率更高。基本思路是：如果能估计出Hessian矩阵的最小特征值 \lambda_{\min}(\mathbf{H}_k)，就可以直接构造一个对角修正矩阵。具体做法是选择

这个公式的含义是：如果最小特征值已经大于阈值 \varepsilon（即 \mathbf{H}_k 已经充分正定），那么 \mathbf{M}_k = \mathbf{0}，不需要修正；如果最小特征值小于 \varepsilon（包括负数情况），就添加一个适当大小的单位矩阵倍数，将最小特征值提升到 \varepsilon。修正后的矩阵 \mathbf{H}_k + \mathbf{M}_k 的最小特征值至少为 \varepsilon，因此是充分正定的。

估计最小特征值有多种方法，比如Gershgorin圆盘定理、幂迭代法的变体等，这些方法的计算代价远低于完整的特征值分解。这种策略在实际应用中更受欢迎，因为它在保证正定性的同时，计算效率较高，且当Hessian矩阵本身接近正定时，对其改动最小。

高斯-牛顿法（Gauss-Newton Method）

最小二乘问题的特殊结构

高斯-牛顿法是专门为一类特殊的优化问题设计的，这类问题的代价函数具有最小二乘的形式。考虑由 m 个分量函数组成的向量值函数 \mathbf{f}(\mathbf{x}) = (f_1(\mathbf{x}), \ldots, f_m(\mathbf{x}))^T，优化目标是

这种形式在系统辨识、参数估计、数据拟合等问题中频繁出现。PPT特别指出这是辨识问题中常遇到的代价函数形式，比如前面讨论过的信道参数估计例子就属于这一类。最小二乘问题的特殊性在于：代价函数本身是向量函数的范数平方，这个结构可以被利用来设计更高效的算法。

高斯-牛顿法的主要思想包含两个步骤：首先对代价函数进行二次近似，然后最小化这个二次近似。与标准牛顿法不同，高斯-牛顿法不直接计算代价函数 f(\mathbf{x}) 的Hessian矩阵，而是通过向量函数 \mathbf{f}(\mathbf{x}) 的Jacobian矩阵来构造二次近似。

通过线性化构造二次近似

方法从向量函数 \mathbf{f}(\mathbf{x}) 在点 \mathbf{x} 处的Taylor展开开始。一阶Taylor展开给出

这里 \mathbf{J}(\mathbf{x}) = \nabla\mathbf{f}(\mathbf{x}) 是向量函数的Jacobian矩阵，它是一个 m \times n 矩阵，第 i 行是标量函数 f_i 的梯度。忽略高阶项后，得到 \mathbf{f}(\mathbf{x} + \mathbf{p}) 的线性近似

现在利用代价函数的特殊形式 f(\mathbf{x} + \mathbf{p}) = \frac{1}{2}\mathbf{f}(\mathbf{x} + \mathbf{p})^T\mathbf{f}(\mathbf{x} + \mathbf{p})。用线性近似 \boldsymbol{\ell}(\mathbf{p}) 替换 \mathbf{f}(\mathbf{x} + \mathbf{p})，得到 f(\mathbf{x} + \mathbf{p}) 的二次近似

展开这个表达式。首先

展开乘积得到

利用对称性 \mathbf{f}^T\mathbf{J}\mathbf{p} = \mathbf{p}^T\mathbf{J}^T\mathbf{f}（两者都是标量且互为转置），合并同类项

因此二次近似为

注意到第一项 \frac{1}{2}\mathbf{f}^T\mathbf{f} = f(\mathbf{x})，因此可以重写为

这正是标准二次函数的形式，其中线性项系数是 \mathbf{J}^T\mathbf{f}，二次项系数矩阵是 \mathbf{J}^T\mathbf{J}。

高斯-牛顿方向的推导

高斯-牛顿法的目标是找到使二次近似最小的方向，即求解

对二次函数求最小值，根据前面建立的理论，需要令其梯度为零。计算 m^Q 关于 \mathbf{p} 的梯度

令梯度为零得到

在迭代点 \mathbf{x}_k 处，记 \mathbf{J}_k = \mathbf{J}(\mathbf{x}_k) 和 \mathbf{f}_k = \mathbf{f}(\mathbf{x}_k)，高斯-牛顿方向 \mathbf{p}_k 满足方程

这个方程称为正规方程（normal equations）。注意到代价函数 f(\mathbf{x}) 的梯度是 \mathbf{g}(\mathbf{x}) = \mathbf{J}(\mathbf{x})^T\mathbf{f}(\mathbf{x})（通过链式法则可以验证），因此正规方程也可以写成

这与牛顿方程 \mathbf{H}_k\mathbf{p}_k = -\mathbf{g}_k 的形式一致，但这里用 \mathbf{J}_k^T\mathbf{J}_k 替代了真实的Hessian矩阵 \mathbf{H}_k。

高斯-牛顿法的正定性和唯一性

二次近似 m^Q(\mathbf{x} + \mathbf{p}) 的Hessian矩阵是 \nabla^2 m^Q(\mathbf{x} + \mathbf{p}) = \mathbf{J}^T(\mathbf{x})\mathbf{J}(\mathbf{x})。这个矩阵是对称的（显然 (\mathbf{J}^T\mathbf{J})^T = \mathbf{J}^T\mathbf{J}）。关键问题是它是否正定。如果Jacobian矩阵 \mathbf{J} 是满秩的（即秩为 n），那么 \mathbf{J}^T\mathbf{J} 是正定的。这是因为对任何非零向量 \mathbf{p}，有 \mathbf{p}^T\mathbf{J}^T\mathbf{J}\mathbf{p} = \|\mathbf{J}\mathbf{p}\|^2 > 0（最后的不等号成立因为 \mathbf{J} 满秩意味着 \mathbf{J}\mathbf{p} \neq \mathbf{0}）。

当 \mathbf{J}^T\mathbf{J} 正定时，二次近似 m^Q(\mathbf{x} + \mathbf{p}) 有唯一的全局最小值，通过正规方程 (\mathbf{J}^T\mathbf{J})\mathbf{p} = -\mathbf{J}^T\mathbf{f} 求得。得到高斯-牛顿方向后，更新迭代点为

PPT指出步长 \alpha 的默认值是1（即完全步长），但也可以通过线搜索来确定更优的步长。使用完全步长意味着完全信任二次近似在整个搜索方向上的准确性，这在接近最优点时通常是合理的，但在远离最优点时可能需要通过线搜索来保证稳定性。

高斯-牛顿法的计算效率优势

PPT特别强调高斯-牛顿法的一个重要优势：矩阵 \mathbf{J}^T\mathbf{J} 是函数 f(\mathbf{x}) 的Hessian矩阵的近似，但计算它比直接计算 \nabla^2 f(\mathbf{x}) 容易得多。为什么？

真实的Hessian矩阵 \mathbf{H}(\mathbf{x}) = \nabla^2 f(\mathbf{x}) 需要计算所有二阶偏导数。对于最小二乘形式的 f(\mathbf{x}) = \frac{1}{2}\sum_{i=1}^m f_i(\mathbf{x})^2，通过链式法则可以得到

第一项正是高斯-牛顿法使用的近似，第二项是额外的二阶项。高斯-牛顿法忽略了第二项，相当于假设残差 f_i(\mathbf{x}) 很小，或者 f_i 的二阶导数很小。在很多实际问题中，特别是当优化接近收敛时，残差确实变得很小，因此这个近似是合理的。

计算 \mathbf{J}^T\mathbf{J} 只需要一阶导数信息（Jacobian矩阵），然后进行矩阵乘法，而直接计算完整的Hessian需要所有二阶导数。对于大规模问题，避免二阶导数的计算可以显著降低计算成本。这使得高斯-牛顿法在最小二乘问题中成为非常实用的选择，它在计算效率和收敛速度之间取得了良好的平衡。

课程总结

通过本课程的学习，应当掌握以下核心内容：

首先是二次代价函数的最小化方法。理解二次函数的结构、最优性条件以及如何通过求解线性方程组 \mathbf{H}\mathbf{x} = -\mathbf{c} 来找到最小化器。掌握判断最小化器存在性和唯一性的条件，即Hessian矩阵的正定性或半正定性的作用。

其次是通用下降算法的框架结构。理解下降算法的基本组成部分：初始化、选择下降方向、执行线搜索、更新迭代点、检查停止条件。这个通用框架是所有具体下降算法的基础，不同算法的区别主要在于如何选择下降方向。

第三是近似线搜索的实现方法。掌握回溯算法的基本流程，理解Armijo条件的含义和作用。知道如何在保证函数值充分下降和避免步长过小之间取得平衡。理解为什么精确线搜索通常不被采用，以及近似线搜索在实际中的实用性。

第四是各种下降方向的计算和评估。理解梯度下降方向的定义及其局限性（对坐标系敏感、收敛可能很慢）。掌握牛顿方向的推导、其利用二阶信息的优势以及二次收敛的性质。理解牛顿方向要求Hessian正定这一前提条件。

最后是处理Hessian矩阵非正定情况的方法。学会修正牛顿法的两种策略：基于谱分解的精确修正和基于最小特征值估计的实用修正。理解高斯-牛顿法如何利用最小二乘问题的特殊结构，通过 \mathbf{J}^T\mathbf{J} 近似Hessian矩阵，在避免直接计算二阶导数的同时保持良好的收敛性能。

这些知识构成了数值优化的基础工具箱，能够应对各种无约束优化问题，从理论理解到算法实现都建立了完整的框架。

梯度、雅可比矩阵和海森矩阵的计算方法

这门课程主要讲解如何计算梯度、雅可比矩阵和海森矩阵，内容分为四个部分：有限差分方法、灵敏度函数、伴随状态方法，以及伴随代码的实现。

有限差分方法

基本思想与标准有限差分

有限差分是计算梯度最直观的数值方法，其核心思想来自于导数的定义。对于一个代价函数 j(\mathbf{p})，其中 \mathbf{p} 是参数向量，我们想要计算它关于参数的梯度 \mathbf{g}(\mathbf{p})。标准的有限差分方法定义梯度的第 i 个分量为：

这个公式的本质是用前向差分来近似偏导数。具体来说，\mathbf{e}_i 是第 i 个标准基向量（只有第 i 个位置是1，其余位置都是0），\mu_i 是一个很小的扰动量。公式计算的是：在参数向量的第 i 个维度上移动一个小步长 \mu_i 后，函数值的变化除以步长，这正是导数定义式的离散化形式。

这个方法的问题在于计算成本。要得到完整的梯度向量，我们需要计算 n_p + 1 次函数 j(\cdot)：一次计算 j(\mathbf{p}) 本身，然后对每个参数维度各计算一次 j(\mathbf{p} + \mu_i \mathbf{e}_i)。当参数维度 n_p 很大时，这个计算量就变得难以承受。

随机方向的有限差分

为了解决高维参数空间中的计算成本问题，我们可以使用随机方向的有限差分。这个方法不是沿着每个坐标轴方向分别扰动，而是沿着一个随机方向同时扰动所有参数。梯度的近似变为：

这里 \mathbf{u} 是随机方向向量，\mu 是扰动的幅度，\phi(n_p) 是一个与参数维度相关的缩放因子。随机向量 \mathbf{u} 有两种常见的选择方式：第一种是从标准正态分布采样 \mathbf{u} \sim \mathcal{N}(\mathbf{0}, \mathbf{I})，此时 \phi(n_p) = 1；第二种是从单位球面上的均匀分布采样 \mathbf{u} \sim \mathcal{U}(\mathcal{S}(0,1))，此时 \phi(n_p) = n_p。

这个方法的巨大优势是只需要两次函数计算，无论参数维度有多高。但代价是引入了估计误差。误差的期望可以分解为两部分：

第一项 O(n_p) \|\mathbf{g}(\mathbf{p})\|_2^2 是随机采样带来的方差，它与参数维度 n_p 成正比，这意味着在高维空间中，单次随机采样的估计会很不准确。第二项是由有限差分步长 \mu 引入的偏差，包含 \mu^2 是因为我们使用的是一阶差分近似。当 \mu 太大时偏差大，\mu 太小时数值误差会放大，需要权衡选择。

小批量采样改进

为了降低随机采样的方差，自然的想法是使用多个随机方向并取平均。这就是小批量采样方法：

这里 b 是批量大小，我们采样 b 个独立的随机方向 \mathbf{u}_i。每个 \mathbf{u}_i 的分布与之前相同：可以是标准正态分布（\phi(n_p) = 1）或单位球面均匀分布（\phi(n_p) = n_p）。

这个方法需要 b + 1 次函数计算（一次计算 j(\mathbf{p})，b 次计算扰动后的函数值），是标准有限差分和单次随机采样之间的折中方案。误差分析显示：

与单次采样相比，方差项从 O(n_p) 降低到 O(n_p/b)，偏差的第一部分也从 O(\mu^2 n_p^3/\phi(n_p)) 降低到 O(\mu^2 n_p^3/[\phi(n_p) b])。这意味着通过增加批量大小 b，我们可以用线性增长的计算成本来显著降低估计误差。但最后一项 O(\mu^2 n_p/\phi(n_p)) 不随 b 变化，它代表了方法本身的固有偏差下界。

灵敏度函数方法

从梯度计算到灵敏度函数

在许多优化问题中，代价函数 j(\mathbf{p}) 实际上是通过输出误差 e(\mathbf{p}, i) 构造的。这里的输出误差定义为测量值与模型预测值之间的差异，例如 e(\mathbf{p}, i) = y(i) - y_m(\mathbf{p}, i)，其中 y(i) 是在第 i 个时刻的实际测量值，y_m(\mathbf{p}, i) 是参数为 \mathbf{p} 时模型的预测输出。当代价函数依赖于这种误差时，计算梯度 \mathbf{g}(\mathbf{p}) 就涉及到计算输出误差关于参数的导数 \partial[e(\mathbf{p}, i)]/\partial \mathbf{p}。

为了简化这个导数的表示和计算，我们引入灵敏度函数的概念。灵敏度函数定义为：

它刻画了输出误差对参数变化的敏感程度。如果某个参数 p_j 的微小变化会导致误差 e(\mathbf{p}, i) 产生很大变化，那么对应的灵敏度分量就会很大。

加权最小二乘例子中的应用

以加权最小二乘代价函数为例，假设我们要最小化：

这里 w_i 是第 i 个观测点的权重。要计算这个函数的梯度，我们对参数求导：

这个公式揭示了灵敏度函数的核心作用：梯度的计算被分解为两部分的乘积——当前的误差值 e(\mathbf{p}, i) 和灵敏度 \mathbf{s}_e(\mathbf{p}, i)。系数2来自对平方项求导，权重 w_i 反映不同数据点的重要性。这意味着只要我们能够计算出灵敏度函数，梯度的计算就变成了一个简单的加权求和问题。

从输出误差到模型输出的灵敏度

关键问题是如何计算灵敏度函数 \mathbf{s}_e(\mathbf{p}, i)。对于输出误差 e(\mathbf{p}, i) = y(i) - y_m(\mathbf{p}, i)，由于测量值 y(i) 不依赖于参数 \mathbf{p}（它是给定的观测数据），我们可以直接推导出：

其中 \mathbf{s}_y(\mathbf{p}, i) = \partial y_m(\mathbf{p}, i)/\partial \mathbf{p} 是模型输出关于参数的灵敏度。负号表明当参数变化使模型输出增大时，误差会相应减小。因此，计算输出误差的灵敏度问题转化为计算模型输出的灵敏度问题。

线性不变系统中的灵敏度计算

当模型是一个线性不变（LI）系统时，模型输出由 m 阶微分方程描述。要计算灵敏度函数 \mathbf{s}_y(\mathbf{p}, i)，对于 i = 1, \ldots, n_p 的每个参数分量，我们需要求解对应的灵敏度微分方程。假设我们忽略初始条件的影响（或初始条件不依赖于参数），那么为了得到模型输出 y_m 和所有参数的灵敏度 \mathbf{s}_y(\mathbf{p}, i)，我们需要模拟求解 2m 阶的微分方程系统——m 阶用于计算 y_m 本身，另外 m 阶用于每个参数的灵敏度。

具体例子：二阶系统的灵敏度

考虑一个由4个参数 p_1, p_2, p_3, p_4 确定的二阶线性微分方程：

其中 u(t) 是输入信号。系统的初始条件设为 y_m(\mathbf{p}, 0) = 0 和 \frac{dy_m(\mathbf{p}, t)}{dt}|_{t=0} = 0。

要计算关于每个参数 p_i 的偏导数（即灵敏度函数 s_i(\mathbf{p}, t)，i = 1, \ldots, 4），我们对原微分方程关于 p_i 求偏导。由于假设初始条件不依赖于参数，所有灵敏度函数的初始条件都是零。对原方程

关于各参数求偏导得到四个灵敏度微分方程：

这些方程的右侧来自对原方程中各项关于相应参数求导。例如，对 p_1 \frac{dy_m}{dt} 关于 p_1 求导时，产生项 \frac{dy_m}{dt}，因此第一个灵敏度方程右侧是 -\frac{dy_m}{dt}。

线性系统的计算优势

表面上看，我们有5个二阶微分方程（1个原方程 + 4个灵敏度方程），这似乎对应于一个10阶的微分方程组。但关键的观察是：这些灵敏度方程都是线性的，并且它们共享相同的左侧算子（即相同的齐次部分 \frac{d^2}{dt^2} + p_1 \frac{d}{dt} + p_2），只有右侧的驱动项不同。这意味着我们可以利用线性系统的叠加原理来高效求解。

具体计算流程是：首先模拟原方程以获得 y_m(\mathbf{p}, t)。在这个模拟过程中，我们通过对 y_m(\mathbf{p}, t) 求导得到 \frac{dy_m}{dt}。然后，我们可以利用这些已计算的量作为灵敏度方程的驱动项，通过求解四个共享相同左侧算子但右侧不同的线性方程来获得 s_1, s_2, s_3, s_4。最后，通过组合 y_m(\mathbf{p}, t) 与灵敏度函数 s_3(\mathbf{p}, t) 和 s_4(\mathbf{p}, t)，我们可以得到完整的输出灵敏度。PPT中的框图展示了这个计算流程的结构，其中积分器（1/s）、参数块（\theta 符号）以及反馈回路清晰地描绘了各量之间的依赖关系。

状态空间表示下的灵敏度计算

前面的二阶系统例子可以用状态空间形式统一描述。我们定义状态向量 \mathbf{x}(t) = [s_3, s_4, s_1, s_2]^{\top}，这个向量包含了4个积分器的输出。注意这里的顺序安排：先放灵敏度 s_3, s_4，后放 s_1, s_2。系统的动态方程写为：

其中系统矩阵和输入矩阵分别为：

矩阵 \mathbf{A} 的结构反映了系统的耦合关系：左上角的 2 \times 2 块对应 s_3, s_4 的动态，右下角的 2 \times 2 块反映了 s_1, s_2 如何依赖前面的状态，而左下角的块则体现了它们之间的耦合。输入矩阵 \mathbf{B} 表明只有第一个状态直接接收输入 u(t)。

我们关心的输出包括模型输出 y_m 和所有灵敏度函数，记为 \mathbf{z}(\mathbf{p}, t) = (y_m, s_1, s_2, s_3, s_4)。这些输出通过状态的线性组合得到：

矩阵 \mathbf{C} 的第一行给出如何从状态 s_3, s_4 构造出 y_m = p_3 s_3 + p_4 s_4，这对应于原微分方程中输入项的处理。后面几行则直接提取各个灵敏度分量。这种状态空间表示使得整个系统的模拟可以用标准的数值积分器高效完成。

非线性不变系统的一般情况

当模型不再是线性时，即模型为非线性不变（NLI）系统时，灵敏度函数的计算需要更通用的方法。一般的连续时间状态空间表示为：

这里函数 \mathbf{f} 和 \mathbf{h} 可以依赖于输入 u(t) 和时间 t。输出灵敏度函数的第 i 个分量通过链式法则展开：

右侧第一项是通过状态的变化间接影响输出，第二项是参数直接出现在输出函数中的贡献。为了计算状态关于参数的导数 \partial \mathbf{x}(t)/\partial p_i，我们定义状态灵敏度 [\mathbf{s}_x(\mathbf{p}, t)]_i = \partial \mathbf{x}(t)/\partial p_i。

对状态方程关于参数 p_i 求导，利用链式法则得到状态灵敏度的动态方程：

这是一个关于 [\mathbf{s}_x(\mathbf{p}, t)]_i 的线性时变微分方程。虽然原系统可能是非线性的，但灵敏度方程总是线性的——这是因为我们对参数的依赖关系求导，导数运算的线性性保证了结果的线性性。方程右侧第一项的雅可比矩阵 \partial \mathbf{f}/\partial \mathbf{x}^{\top} 和第二项 \partial \mathbf{f}/\partial p_i 都依赖于时间（通过 \mathbf{x}(t) 的变化），所以这是非平稳的。

灵敏度函数方法的计算流程

完整的计算流程是：首先进行一次模拟求解原系统以获得状态轨迹 \mathbf{x}(t)。然后对每个参数 p_i（i = 1, \ldots, n_p），我们需要额外进行一次模拟来求解对应的状态灵敏度微分方程，从而得到 [\mathbf{s}_x(\mathbf{p}, t)]_i。因此总共需要 n_p 次额外的模拟。

这个方法的关键特点是：虽然每个灵敏度微分方程都是线性的且非平稳的（因为雅可比矩阵随时间变化），但它的计算复杂度与有限差分方法相同——都需要 O(n_p) 次函数或微分方程的求解。然而，灵敏度函数方法的巨大优势在于它不存在有限差分中的数值近似误差。灵敏度方程是通过解析微分得到的精确方程，只要数值积分器的精度足够高，我们就能获得非常准确的梯度信息。相比之下，有限差分总是受到步长选择的困扰：步长太大导致截断误差，步长太小则放大舍入误差。

伴随状态方法

离散时间系统的基本设定

伴随状态方法处理的是离散时间状态空间表示的模型：

这里 k 是离散时间步，\mathbf{x}(k) 是第 k 步的状态向量。我们考虑一个可加形式的代价函数：

这个代价函数将每个时间步的即时代价 r_i 累加起来。这种可加结构在实际应用中很常见，例如轨迹跟踪问题中每个时刻的跟踪误差平方和，或者在强化学习中累积的奖励。需要注意这里的"可加"要求相对宽松——函数可以在每个时刻都不同，关键是不能有跨时刻的交叉项。

扩展状态的引入

为了用伴随方法处理这个问题，我们引入一个技巧：定义额外的标量状态 x^0(k) 来累积代价。这个辅助状态满足递推关系：

这样一来，原本分布在各个时刻的代价求和问题就被转化为终端状态 x^0(n+1) 的值，即 j(\mathbf{p}) = x^0(n+1)。我们现在定义扩展状态向量：

这个扩展状态有三个部分：累积代价 x^0(k)、原始状态 \mathbf{x}(k)、以及参数 \mathbf{p}。将参数也放入状态向量是一个关键技巧——虽然参数本身不随时间变化，但把它当作状态的一部分可以统一处理依赖关系。扩展状态的演化方程变为：

参数部分的演化方程就是 \mathbf{p} \to \mathbf{p}（恒等映射）。代价函数现在表达为 j(\mathbf{p}) = [1 \; 0 \; \ldots \; 0] \mathbf{x}_e(n+1)，即只提取扩展状态的第一个分量。

梯度的链式法则推导

现在我们要计算梯度 \mathbf{g}(\mathbf{p}) = \partial j(\mathbf{p})/\partial \mathbf{p}。由于 j(\mathbf{p}) = x^0(n+1) 只是终端扩展状态的一个分量，利用链式法则：

其中 \partial j/\partial \mathbf{x}_e(n+1) = [1 \; 0 \; \ldots \; 0]^{\top}。接下来需要计算扩展状态关于参数的雅可比矩阵 \partial \mathbf{x}_e^{\top}(n+1)/\partial \mathbf{p}。由于扩展状态在每一步都依赖于前一步，我们递归地应用链式法则。从 \mathbf{x}_e(n+1) = \mathbf{f}_e[\mathbf{x}_e(n)] 出发：

继续向前递推，我们得到从 k=n 一直到 k=1 的类似关系：

直到初始条件：

伴随状态的定义与递推

直接展开这个乘积链会导致计算复杂度随时间步数指数增长。伴随方法的巧妙之处在于反向累积这些雅可比矩阵。我们定义伴随状态向量 \mathbf{d}_e(k) 满足时间反向的递推关系：

终端条件（也就是递推的起点）设为：

这对应于 \partial j/\partial \mathbf{x}_e(n+1)。通过这个反向递推，当我们到达初始时刻时，梯度可以表达为：

初始扩展状态对参数的导数结构很清楚：第一行（累积代价）初始为0不依赖参数；第二部分是初始原始状态对参数的导数；第三部分是参数对自身的导数（单位矩阵 \mathbf{I}_p）。

伴随状态方法的计算效率

伴随状态方法的核心优势在于计算复杂度。整个梯度计算只需要：一次前向模拟计算状态轨迹 \mathbf{x}_e(k) 从 k=0 到 k=n+1（这给出代价 j(\mathbf{p})），然后一次反向模拟计算伴随状态 \mathbf{d}_e(k) 从 k=n+1 反向到 k=0（这给出梯度 \mathbf{g}(\mathbf{p})）。无论参数维度 \dim(\mathbf{p}) 有多大，我们都只需要这两次模拟就能得到完整的梯度向量，不存在任何近似误差。

相比之下，灵敏度函数方法需要 n_p 次额外的前向模拟（每个参数一次），有限差分方法也需要至少 n_p + 1 次函数计算。因此当参数维度很高时，伴随方法的优势极为显著。PPT最后提到的备注进一步说明了伴随状态的数学性质：伴随状态满足扩展状态演化的线性化（雅可比转置），并且沿着轨迹内积 \partial \mathbf{x}_e(k+1)^{\top} \mathbf{d}_e(k+1) = \partial \mathbf{x}_e(k)^{\top} \mathbf{d}_e(k) 保持不变，这可以用来验证计算的正确性。

伴随代码的实现

程序作为离散时间系统

伴随方法不仅适用于显式的状态空间模型，还可以直接应用于计算代价函数的程序本身。核心思想是将计算程序 j(\mathbf{p}) 看作一个离散时间系统，其中程序的每一行对应一个时间步。设 \mathbf{v} 为程序中所有变量的向量，程序的第 k 行更新某个单一变量（比如索引为 \mu(k) 的变量）：

这里 \mathcal{I}_k 是第 k 行需要用到的变量索引集合，\phi_k 是该行执行的运算（可能是加法、乘法、函数调用等）。这个更新可以用一个作用在整个变量向量上的变换 \Phi_k 来表示：

即除了被更新的变量 v_{\mu(k)}，其他所有变量保持不变。程序从初始到结束的全过程可以表示为：

其中 f 是程序的最后一行。这样，程序的执行就等价于一个 f 步的离散时间系统演化。

变量向量的约定排序

为了系统地应用伴随方法，我们对变量向量 \mathbf{v} 采用特定的排序约定。向量的前 n_p 个分量是参数 \mathbf{p}，这些是我们要对其求梯度的自变量。接下来是独立变量（对于计算 j(\mathbf{p}) 来说是给定的输入数据），例如测量值 y(i)（i = 1, \ldots, n）。然后是依赖变量，即程序中的中间计算结果。最后一个分量是代价函数的值 j(\mathbf{p}) 本身。按照这个约定：

即提取执行完所有 f 行代码后变量向量的最后一个分量。

链式法则与梯度的反向传播

要计算梯度 \mathbf{g}(\mathbf{p}) = \partial j/\partial \mathbf{p}，我们从终端开始应用链式法则。由于 j(\mathbf{p}) 依赖于 \mathbf{v}(f)，而 \mathbf{v}(f) 又依赖于之前所有步骤的变量，梯度可以展开为：

这是一个从 k=1 到 k=f 的雅可比矩阵乘积链。直接从左到右计算会非常低效，因为中间会产生大的矩阵。伴随代码的策略是从右到左反向累积。

对偶变量的引入与递推

我们定义对偶变量（伴随变量）\mathbf{d}(k) 来表示代价函数对第 k 步变量向量的梯度。终端条件初始化为：

然后从 k=f 反向递推到 k=1：

当反向传播完成后，梯度就是对偶变量向量的前 n_p 个分量（对应参数部分）：

因为初始变量向量的前 n_p 个分量就是参数自己（对自己求导是单位矩阵），后面的分量是独立变量不依赖参数（导数为零矩阵）。因此 \mathbf{g}(\mathbf{p}) 就是 \mathbf{d}(0) 的前 n_p 个分量，而 \mathbf{d}(0) 的其余分量表示 j(\cdot) 关于独立变量（如测量值 y(i)）的导数。

雅可比矩阵的稀疏结构

第 k 行代码只更新一个变量 v_{\mu(k)}，因此变换 \Phi_k 的雅可比矩阵 \partial \Phi_k^{\top}/\partial \mathbf{v}(k-1) 具有非常稀疏的结构。这个矩阵几乎是单位矩阵，只有第 \mu(k) 列不同，该列包含 \phi_k 关于其参数的偏导数。具体来说，这个雅可比矩阵的第 (i, j) 元素为：

PPT中的矩阵图示清晰地展示了这个结构：对角线上全是1（表示未被更新的变量保持不变），而第 \mu(k) 行（对应被更新的变量）包含了一组偏导数 \partial \phi_k/\partial \mathbf{v}(k-1)。

对偶变量的更新规则

利用这个稀疏结构，对偶变量的更新可以高效地逐分量计算。对于不是被更新变量的索引 i \neq \mu(k)：

这表明该变量的对偶值保持 d_i(k) 不变，但如果该变量参与了 \phi_k 的计算（即 i \in \mathcal{I}_k），还要加上一个反向传播来的贡献。对于被更新的变量索引 i = \mu(k)：

这个公式捕捉了如果 \phi_k 本身依赖于旧的 v_{\mu(k)} 值（例如累加操作 v \leftarrow v + \text{something}），梯度如何通过这种自依赖传播。最终梯度 \mathbf{g}(\mathbf{p}) 就是 \mathbf{d}(0) 的前 n_p 个分量，而接下来的分量给出了代价函数关于独立变量（如观测数据）的灵敏度。

伴随代码的内存优化

在实际实现伴随代码时，一个重要的优化是不需要保存对偶变量 \mathbf{d} 在每个时间步的完整历史。对于第 k 行代码的更新 v_{\mu(k)} = \phi_k(\mathbf{v}_i, i \in \mathcal{I}_k)，伴随代码中的对偶变量更新规则可以就地进行。对于所有参与计算的变量索引 i \in \mathcal{I}_k 但 i \neq \mu(k)：

对于被更新的变量本身：

这里有两种常见的特殊情况。第一种是初始化操作：如果第 k 行是 v_{\mu(k)} = \phi_k(\mathbf{v}_i, i \neq \mu(k))，即 \phi_k 不依赖于旧的 v_{\mu(k)} 值，那么对应的伴随代码中 d_{\mu(k)} = 0。第二种是增量操作：如果第 k 行是 v_{\mu(k)} = v_{\mu(k)} + \psi_k(\mathbf{v}_i, i \neq \mu(k))（即在原值基础上累加），那么对应的伴随代码是 d_{\mu(k)} = d_{\mu(k)}（保持不变），因为 \partial \phi_k/\partial v_{\mu(k)} = 1。

简单例子：平方误差的累积

考虑一个典型的代价累积代码行：

这行代码涉及3个变量：代价累积值 j、测量值 y(k) 和模型预测 y_m(k)。相应地，我们需要3个对偶变量：dj、dy(k) 和 dy_m(k)。为了推导伴随代码，我们需要计算这个操作关于各个输入的偏导数。记 e(k) = y(k) - y_m(k) 为误差，则 \partial[e^2]/\partial y(k) = 2e(k) 且 \partial[e^2]/\partial y_m(k) = -2e(k)。根据链式法则，伴随代码的更新为：

这里 dy(k) 和 dy_m(k) 累积了来自代价函数的梯度贡献。注意 dj 本身保持不变（增量操作的特性）。第三行的 dj = dj 实际上是无用的（标记为useless），因为它不改变任何值，只是形式上的完整性。

伴随代码的实现备注

PPT给出了实现伴随代码时的6个重要备注。第一，伴随代码的长度大约是直接代码的4倍，这是因为每个操作都需要相应的梯度计算和累积。第二，直接代码中的for循环在伴随代码中需要反向执行——如果直接代码是从 i=1 到 n 的循环，伴随代码要从 i=n 反向到 i=1。第三，分支条件（if语句）在伴随代码中必须保持相同，这是为了确保伴随代码沿着与直接代码完全相同的执行路径。第四，如果直接代码中有非线性表达式的变量值，这些值需要在伴随代码执行时保持可用（要么重新计算，要么从前向传播时保存）。第五，对于每个直接代码的子程序，应该构建对应的伴随子程序。第六，可以开发工具来并行执行或更新直接代码和伴随代码，或者使用自动微分库如JAX（https://github.com/google/jax）。最重要的结论是：梯度 \mathbf{g}(\mathbf{p}) 通过两次模拟获得（前向一次，反向一次），完全没有近似误差。

完整例子：递推系统的伴随代码

现在看一个更完整的例子。直接代码实现了一个简单的递推系统：

J = 0;
ym(0) = p2;
for i = 1 to nt do {
ym(i) = p1 * ym(i-1);
J = J + [y(i) - ym(i)]²;
}

这个代码有两个参数 p1 和 p2，其中 p2 是初始条件，p1 是递推系数。模型输出 ym(i) 通过递推生成，代价函数累积每个时刻的平方误差。对应的伴随代码结构如下：首先初始化所有对偶变量 dJ = 1（因为我们要计算 \partial J/\partial \mathbf{p}），以及对所有 i 设置 dy(i) = 0 和 dym(i) = 0。然后是关键的双层结构：第一，设置 dp1 = 0 和 dp2 = 0；第二，执行反向循环 i = nt 下降到 i = 1。在反向循环内部，我们对直接代码中第二条指令（代价累积）进行对偶化，得到：

接下来对偶化第一条指令（递推关系 ym(i) = p1 \cdot ym(i-1)）。这个操作依赖于 p1 和 ym(i-1)，因此：

完成反向循环后，我们还需要对偶化循环前的初始化 ym(0) = p2，这给出 dp2 = dp2 + dym(0)。最后，由于初始化时 dym(0) = 0 不再需要，我们可以直接设置 dym(0) = 0。最终，dp1 和 dp2 中存储的就是梯度 \partial J/\partial p1 和 \partial J/\partial p2。

这个例子清楚地展示了伴随代码的核心特征：反向执行、对偶变量的累积更新、以及如何处理循环结构。通过一次前向模拟（计算 J 和所有 ym(i)）和一次反向模拟（计算梯度），我们获得了完整且精确的梯度信息，无论参数维度有多少。

递归最小二乘与卡尔曼滤波器

本课程聚焦于递归最小二乘算法和卡尔曼滤波器，这两者是参数估计和状态估计中的核心方法。课程分为三个主要部分：首先介绍经典的最小二乘估计，然后推导递归最小二乘算法，最后引入卡尔曼滤波器及其各种应用场景。

最小二乘估计

参数线性模型的定义

我们考虑一类特殊但非常重要的模型——参数线性模型。这类模型的输出可以表示为参数向量的线性组合：

这个模型结构中，\mathbf{r}^T(t) 是已知的回归向量（regressor），它在时刻 t 的值完全已知，不依赖于待估计的参数。\mathbf{p} 是我们要估计的参数向量。y_m(t+1, \mathbf{p}) 是模型在时刻 t+1 的输出预测值。

这里需要强调的是，索引 t 不一定表示物理时间上的离散时间步。在许多情况下，t 可以理解为第 t 个实验条件或第 t 组测量数据。因此这个框架可以处理时间序列数据，也可以处理不同实验条件下的静态测量。

批量数据的矩阵形式

当我们收集多个时刻（或多组实验）的数据时，可以将所有模型输出组织成一个向量。对于 t = 1, \ldots, n_t 的数据，模型输出向量为：

这里定义的矩阵 \mathbf{R} 是回归矩阵，它的第 i 行是 \mathbf{r}^T(i-1)。这个矩阵形式揭示了整个估计问题的本质：我们有一个线性方程组 \mathbf{y}_m = \mathbf{R}\mathbf{p}，需要从观测数据 \mathbf{y} 来求解参数 \mathbf{p}。

双线性模型例子

一个重要的例子是双线性模型。这类模型在系统辨识中很常见，它同时包含输入和输出的历史信息，以及它们的乘积项：

这个模型包含三类项：第一类是输入 u 的滞后项（自回归部分），第二类是输出 y 的滞后项（移动平均部分），第三类是输入与输出滞后项的交叉乘积（双线性项）。

虽然双线性模型看起来复杂，但关键是它仍然满足参数线性的性质。我们可以构造回归向量 \mathbf{r}^T(t)，它包含所有必要的已知量：

这个回归向量将所有输入滞后值、输出滞后值，以及它们的乘积都排列在一起，使得模型输出可以写成 y_m(t) = \mathbf{r}^T(t) \mathbf{p} 的形式。

二次代价函数的构造

估计问题的目标是找到参数 \mathbf{p}，使得模型输出与实际测量值尽可能接近。我们采用二次代价函数：

这里 \mathbf{e}(\mathbf{p}) 是误差向量，定义为测量值与模型预测值之差：

权重矩阵 \mathbf{W} 是对称半正定的，它反映了我们对不同测量点误差的重视程度。基本假设是测量次数远多于参数个数，即 \dim \mathbf{y} \gg \dim \mathbf{p}，这保证了问题的超定性（方程多于未知数）。

权重矩阵的选择有统计含义。最简单的情况是 \mathbf{W} = \mathbf{I}，即所有测量点等权重。如果有先验信息表明测量噪声不是均匀的，可以根据噪声的协方差结构选择权重。一个特别重要的例子是：假设真实数据由真实参数 \mathbf{p}^* 加上噪声生成，即 \mathbf{y} = \mathbf{y}_m(\mathbf{p}^*) + \mathbf{b}，其中噪声 \mathbf{b} 是零均值随机向量，协方差矩阵为 \boldsymbol{\Sigma}。在这种情况下，最大似然估计理论告诉我们应该选择 \mathbf{W} = \boldsymbol{\Sigma}^{-1}。这个选择使得估计器具有最优统计性质。

最小二乘解的推导

将误差的定义代入代价函数，我们得到关于参数的显式二次型：

要找到最小化这个代价的参数估计值 \hat{\mathbf{p}}_{\text{LS}}（LS代表Least-Squares），我们对参数求导并令其为零。由于 \mathbf{W} 是对称矩阵，代价函数对参数的梯度为：

这个梯度为零的条件称为正规方程（normal equation）。系数-2可以消去，整理后得到：

如果矩阵 \mathbf{R}^T \mathbf{W} \mathbf{R} 是可逆的，我们就能解出参数估计的显式表达式：

这就是著名的加权最小二乘解。注意这个解是解析的，不需要任何迭代或初始化——只要给定数据 \mathbf{y} 和回归矩阵 \mathbf{R}，就能直接计算出估计值。

最小二乘估计的几何性质

最小二乘估计具有深刻的几何意义。第一个重要性质是，在最优估计处，残差与模型输出之间存在正交关系：

当权重矩阵取单位矩阵 \mathbf{W} = \mathbf{I} 时，这个关系特别直观——残差向量 \mathbf{e}(\hat{\mathbf{p}}_{\text{LS}}) 与模型输出向量 \mathbf{y}_m(\hat{\mathbf{p}}_{\text{LS}}) 在欧几里得空间中正交。PPT中的图示清晰地展示了这个几何关系：测量向量 \mathbf{y} 被投影到模型响应空间 \mathcal{S}_r 上，投影点就是 \mathbf{y}_m(\hat{\mathbf{p}}_{\text{LS}})，而残差 \mathbf{e}(\hat{\mathbf{p}}_{\text{LS}}) 垂直于这个空间。

第二个性质是关于投影的刻画。定义模型响应的超曲面为：

由于模型关于参数是线性的（即 \mathbf{y}_m(\mathbf{p}) = \mathbf{R}\mathbf{p}），这个"超曲面"实际上是一个超平面——它是 \mathbf{R} 的列空间。最小二乘解给出的 \mathbf{y}_m(\hat{\mathbf{p}}_{\text{LS}}) 正是测量向量 \mathbf{y} 在这个超平面上的正交投影。投影算子可以显式写为：

这个投影算子是幂等的（\boldsymbol{\Pi}^2 = \boldsymbol{\Pi}）且对称的，它将任意向量映射到 \mathbf{R} 的列空间中。

估计不确定性的信息

最小二乘解还有几个实用的性质。第一，矩阵 \mathbf{R}^T \mathbf{W} \mathbf{R} 是对称的，维度为 \dim \mathbf{p} \times \dim \mathbf{p}。在非病态情况下（即回归矩阵 \mathbf{R} 列满秩），这个矩阵是正定的，因此可逆。病态情况对应于参数之间存在共线性——某些参数组合无法从数据中区分出来。

第二，虽然最小二乘给出了闭式解，但在非线性优化问题中，\hat{\mathbf{p}}_{\text{LS}} 可以作为迭代算法的初始猜测。这是因为最小二乘解已经利用了所有数据，通常会给出一个相当合理的起点。

第三，也是最重要的一点，矩阵 (\mathbf{R}^T \mathbf{W} \mathbf{R})^{-1} 包含了参数估计不确定性的关键信息。在统计推断中，如果噪声 \mathbf{b} 是高斯分布且 \mathbf{W} = \boldsymbol{\Sigma}^{-1}，那么 (\mathbf{R}^T \mathbf{W} \mathbf{R})^{-1} 就是参数估计的协方差矩阵。它的对角元素给出各个参数估计的方差，非对角元素描述参数之间的相关性。这个矩阵的逆矩阵 \mathbf{R}^T \mathbf{W} \mathbf{R} 有时被称为信息矩阵或Fisher信息矩阵。

数值实现的技巧

在实际计算中，我们可以通过矩阵分解来简化问题。任何对称正定矩阵 \mathbf{W} 都可以用Cholesky分解表示为：

其中 \mathbf{U} 是可逆的上三角矩阵。这个分解在数值上非常稳定且高效。利用这个分解，代价函数可以改写为：

这表明，通过变量替换——用 \mathbf{U}\mathbf{y} 替换 \mathbf{y}，用 \mathbf{U}\mathbf{R} 替换 \mathbf{R}——我们可以将加权最小二乘问题转化为标准的（\mathbf{W} = \mathbf{I}）最小二乘问题。这种转化在数值实现中很有用，因为标准最小二乘问题有更成熟的数值算法，例如QR分解或奇异值分解（SVD）。通过这种预处理，我们既保留了加权的灵活性，又能利用标准算法的数值稳定性。

数值实现的稳定性问题

虽然最小二乘解的公式 \hat{\mathbf{p}}_{\text{LS}} = (\mathbf{R}^T \mathbf{W} \mathbf{R})^{-1} \mathbf{R}^T \mathbf{W} \mathbf{y} 看起来简洁明了，但直接按照这个公式计算在数值上并不高效。问题的根源在于需要先计算矩阵乘积 \mathbf{R}^T \mathbf{W} \mathbf{R}，然后求其逆矩阵，这两步都可能放大数值误差。

另一种方法是直接求解正规方程组：

这个方法计算量更小，但问题在于它可能数值不稳定。当矩阵 \mathbf{R}^T \mathbf{W} \mathbf{R} 接近奇异（即病态）时，求解线性方程组会遭遇严重的数值精度损失。

对于病态问题（\mathbf{R}^T \mathbf{W} \mathbf{R} 几乎奇异），奇异值分解（SVD）提供了更稳定的数值方案。

奇异值分解的应用

回归矩阵 \mathbf{R} 可以进行奇异值分解：

这个分解将矩阵分解为三部分。矩阵 \mathbf{S} 是对角矩阵，维度为 \dim \mathbf{p} \times \dim \mathbf{p}。矩阵 \mathbf{U} 的维度是 \dim \mathbf{y} \times \dim \mathbf{p}（与 \mathbf{R} 相同），满足正交性条件 \mathbf{U}^T \mathbf{U} = \mathbf{I}_{\dim \mathbf{p}}。矩阵 \mathbf{V} 的维度是 \dim \mathbf{p} \times \dim \mathbf{p}，满足 \mathbf{V}^T \mathbf{V} = \mathbf{I}_{\dim \mathbf{p}}。

对角矩阵 \mathbf{S} 的对角元素 s_{i,i} \geqslant 0 称为 \mathbf{R} 的奇异值，它们提供了关于矩阵数值条件的关键信息。矩阵 \mathbf{R}^T \mathbf{R} 是奇异的当且仅当 \mathbf{R} 有某些奇异值为零。更一般地，条件数（condition number）定义为最大奇异值与最小奇异值之比，它衡量了矩阵的数值稳定性。当条件数接近浮点数精度的倒数（双精度下约为 10^{16}）时，矩阵在数值计算中就被认为是病态的。

通过SVD改进计算

在标准最小二乘问题（\mathbf{W} = \mathbf{I}）中，需要求逆的矩阵是：

这里利用了 \mathbf{U}^T \mathbf{U} = \mathbf{I}。这个表达式揭示了两个重要事实：第一，\mathbf{R} 的奇异值是 \mathbf{R}^T \mathbf{R} 的特征值的平方根。具体来说，如果 s_i 是 \mathbf{R} 的第 i 个奇异值，那么 s_i^2 就是 \mathbf{R}^T \mathbf{R} 的第 i 个特征值。第二，\mathbf{R}^T \mathbf{R} 的条件数等于 \mathbf{R} 的条件数的平方。这意味着直接处理 \mathbf{R}^T \mathbf{R} 会使数值问题恶化。

为了避免这种条件数的平方恶化，可以利用SVD直接计算最小二乘解。从正规方程出发：

将SVD代入并逐步化简：

由于 \mathbf{V} 是正交矩阵，(\mathbf{V}^T)^{-1} = \mathbf{V} 且 \mathbf{V}^{-1} = \mathbf{V}^T，因此：

这个公式数值上更稳定，因为它只涉及一次矩阵求逆（\mathbf{S}^{-1}），而 \mathbf{S} 是对角矩阵，其逆只需要翻转对角元素。在正常情况下（\dim \mathbf{y} > \dim \mathbf{p}），所有奇异值都严格为正，存在唯一的最小二乘估计。

递归最小二乘

递归算法的动机

前面讨论的批量最小二乘算法一次性处理所有测量数据，这种方式称为批处理。然而在许多应用场景中，递归算法更有优势。递归算法在数据到达时逐个处理，这种方式称为在线处理。

递归算法的第一个优势是内存效率。当系统内存有限，无法存储所有历史数据时，我们仍然希望获得与批量最小二乘相同的估计结果，但以更少的内存代价。在这种情况下，我们假设真实参数 \mathbf{p}^* 是常数。

第二个优势是适应性。在某些应用中，需要根据观测数据实时做出决策，例如实时诊断或自适应控制。在这些场景中，参数 \mathbf{p}^* 可能随时间缓慢变化，递归算法能够跟踪这种变化。

参数恒定情况下的递归最小二乘

现在我们推导当真实参数 \mathbf{p}^* 保持恒定时的递归算法。假设我们已经利用前 t 个测量值计算出了参数估计 \hat{\mathbf{p}}_{\text{LS}}(t)，其表达式为：

这里的回归矩阵和测量向量分别是：

当第 t+1 个测量值到达时，回归矩阵和测量向量需要扩充新的一行：

我们的目标是在不存储完整的 \mathbf{R} 和 \mathbf{y} 的情况下计算 \hat{\mathbf{p}}_{\text{LS}}(t+1)，从而解决内存问题。

引入辅助变量的递推形式

将批量最小二乘公式展开为求和形式，矩阵乘积可以写为：

为了方便递推，我们引入两个辅助矩阵/向量：

矩阵 \mathbf{M}(t) 是累积的信息矩阵，它是对称正定的（在数据充分的情况下）。向量 \mathbf{v}(t) 是累积的加权测量。利用这两个量，参数估计可以简洁地写为：

这种表示法允许递归实现。当新的测量到达时，\mathbf{M} 和 \mathbf{v} 的更新规则为：

这两个更新式非常简单——每次只需要加上新的外积项和新的加权测量项。内存问题得到解决：我们只需要存储当前的 \mathbf{M}(t) 和 \mathbf{v}(t)，而不需要保留全部历史数据 \mathbf{R} 和 \mathbf{y}。

避免矩阵求逆的代数变换

虽然上面的递推形式解决了内存问题，但每一步都需要计算矩阵逆 \mathbf{M}^{-1}(t+1)，这在计算上仍然较昂贵。通过巧妙的代数变换，我们可以避免在每一步都进行完整的矩阵求逆。从递推更新出发：

由于 \hat{\mathbf{p}}_{\text{LS}}(t) = \mathbf{M}^{-1}(t) \mathbf{v}(t)，即 \mathbf{v}(t) = \mathbf{M}(t) \hat{\mathbf{p}}_{\text{LS}}(t)，代入上式：

现在利用 \mathbf{M}(t+1) = \mathbf{M}(t) + \mathbf{r}(t) \mathbf{r}^T(t)，我们可以将 \mathbf{M}(t) 表示为 \mathbf{M}(t+1) - \mathbf{r}(t) \mathbf{r}^T(t)：

展开并重新组合：

这个公式揭示了递归最小二乘的核心结构：新的估计等于旧的估计加上一个校正项。

预测-校正结构的解释

将上面推导出的递推公式进行重新组织，可以识别出三个关键成分。首先是一步超前预测：

这是用当前参数估计对下一时刻输出的预测。其次是预测误差：

这个误差衡量了实际测量值与模型预测值之间的差异。第三是校正增益向量：

这个增益向量在后续内容中会与卡尔曼增益联系起来。

利用这三个量，递归最小二乘的更新公式可以写为直观的预测-校正形式：

这个结构清楚地表明：参数估计的更新量正比于预测误差。当模型预测完全准确（e(t+1) = 0）时，参数估计保持不变，即 \hat{\mathbf{p}}_{\text{LS}}(t+1) = \hat{\mathbf{p}}_{\text{LS}}(t)。这符合直觉——如果模型已经能完美预测新数据，就没有必要调整参数。增益向量 \mathbf{k}(t+1) 则决定了误差如何转化为参数的校正，它依赖于累积的数据信息 \mathbf{M}(t+1) 和当前的回归向量 \mathbf{r}(t)。

矩阵求逆引理的应用

虽然递归形式 \hat{\mathbf{p}}_{\text{LS}}(t+1) = \hat{\mathbf{p}}_{\text{LS}}(t) + \mathbf{M}^{-1}(t+1) \mathbf{r}(t) (y(t+1) - \mathbf{r}^T(t) \hat{\mathbf{p}}_{\text{LS}}(t)) 避免了存储全部历史数据，但每一步都需要计算 \mathbf{M}^{-1}(t+1)。为了避免这个矩阵求逆操作，我们利用矩阵求逆引理（也称为Sherman-Morrison-Woodbury公式）：

这个引理的关键作用是将大矩阵的求逆转化为小矩阵的求逆。定义 \mathbf{P}(t+1) = \mathbf{M}^{-1}(t+1)，它就是我们要递推更新的矩阵。注意到：

将矩阵求逆引理应用于这个表达式，令 \mathbf{A} = \mathbf{M}^{-1}(t) = \mathbf{P}(t)，\mathbf{B} = \mathbf{r}(t)，\mathbf{C} = 1（标量），\mathbf{D} = \mathbf{r}^T(t)，我们得到：

这里矩阵求逆被一个标量除法替代，因为括号内的项 1 + \mathbf{r}^T(t)\mathbf{P}(t)\mathbf{r}(t) 是一个标量。最终的更新公式可以简化为：

这个公式的计算复杂度大大降低——从原本的 O(n_p^3)（矩阵求逆）降到 O(n_p^2)（矩阵-向量乘法和秩一更新）。

利用这个 \mathbf{P}(t+1) 的更新公式，校正增益 \mathbf{k}(t+1) 也可以得到更高效的计算形式。从定义出发：

代入 \mathbf{P}(t+1) 的表达式并简化：

这个简洁的公式表明，增益向量只需要一次矩阵-向量乘法和一个标量除法就能计算出来。

完整的递归算法流程

现在我们可以总结完整的递归最小二乘算法。每次迭代分为两个阶段：在新测量到达之前的准备阶段，以及测量到达后的更新阶段。

在测量 y(k+1) 到达之前，算法执行三个步骤。第一步是计算校正增益：

第二步是更新协方差矩阵（或称为估计误差协方差）：

这个形式比之前的公式更简洁，因为已经利用了增益向量的定义。第三步是预测下一时刻的输出：

一旦测量 y(k+1) 到达，算法立即执行参数更新：

这个更新使用了预先计算好的增益和预测值，因此非常高效。

算法的初始化策略

递归算法需要初始值 \hat{\mathbf{p}}_{\text{LS}}(0) 和 \mathbf{P}(0)。有两种常见的初始化方案。第一种方案是等待直到累积了足够的数据使得 \mathbf{M}(t) 变为可逆，然后用批量公式初始化：

从这个时刻开始，递归算法与批量算法在数值误差范围内完全等价。这种方案保证了算法的一致性，但需要等待初始阶段。

第二种方案是直接从零开始，设置 \hat{\mathbf{p}}_{\text{LS}}(0) = \mathbf{0}，同时设置 \mathbf{P}(0) 为一个很大的正定矩阵，例如 \mathbf{P}(0) = 10^{24}\mathbf{I}。这个选择的统计意义是：矩阵 \mathbf{P}(t) 表征参数估计的不确定性，一个很大的 \mathbf{P}(0) 反映了初始时我们对参数几乎一无所知。实际上，初值的精确选择并不太关键，因为 \mathbf{P}(t) 通常会在算法开始后迅速减小——随着数据的累积，估计不确定性快速下降。

递归最小二乘估计器的性质

递归最小二乘算法有几个重要的理论和实践性质。第一，算法中没有显式的矩阵求逆操作，可辨识性问题被隐藏起来了。如果回归向量存在共线性，问题会通过 \mathbf{P}(t) 的数值问题体现出来，而不是直接的求逆失败。

第二，矩阵 \mathbf{P} 的更新量具有明确的符号：

这个差值是对称且半负定的。这意味着（除非 \mathbf{r}(t) = \mathbf{0}）\mathbf{P}(t) 随时间单调递减，反映了参数估计不确定性的持续降低。相应地，校正增益 \mathbf{k}(t) 也会递减——这表明随着数据积累，新观测到的预测误差对参数更新的影响逐渐减弱。这符合直觉：早期的数据信息含量高，后期的数据对已经收敛的估计贡献较小。

第三，递归算法与批量最小二乘有相同的数值稳定性问题。当问题病态时（例如回归向量接近共线），可以采用类似的正则化方法，如向 \mathbf{P} 添加小的正定项 \epsilon\mathbf{I}（Levenberg-Marquart正则化）。另一个有效的数值技巧是对 \mathbf{P}(t) 使用Cholesky分解 \mathbf{P}(t) = \mathbf{U}(t)\mathbf{D}(t)\mathbf{U}^T(t)，其中 \mathbf{U}(t) 是单位上三角矩阵，\mathbf{D}(t) 是对角矩阵。通过直接更新 \mathbf{U}(t) 和 \mathbf{D}(t)，可以保持数值精度。

第四，由于 \mathbf{P}(t) 刻画了估计的不确定性（在统计框架下是协方差矩阵），我们可以跟踪其对角元素来监测各个参数估计的置信度。

参数时变情况：指数遗忘因子

前面的递归算法假设真实参数 \mathbf{p}^* 是恒定的。但在许多实际应用中，系统参数可能会随时间缓慢演变（例如由于老化、环境变化或工作点漂移）。在这种情况下，我们希望估计器 \hat{\mathbf{p}}_{\text{LS}}(t) 能够跟踪这种演变。

为了实现参数跟踪，关键思想是：旧的测量应该被赋予较小的权重，相当于让算法"遗忘"过去的数据。这通过在代价函数中引入指数遗忘因子来实现。具体来说，在时刻 t+1 计算代价时，我们给不同时刻的测量分配指数递减的权重：最新的测量 y(t+1) 权重为1，前一个测量 y(t) 权重为 \lambda，再前一个 y(t-1) 权重为 \lambda^2，以此类推。

参数 \lambda 称为遗忘因子，满足 0 < \lambda \leqslant 1。当 \lambda = 1 时，退化为标准的递归最小二乘（所有数据等权）。在实际应用中，\lambda 通常取接近1的值（如0.95到0.99），这样既能保留大部分历史信息的稳定性，又能对参数变化有足够的响应速度。\lambda 越小，算法对参数变化的响应越快，但估计的方差也会增大。

带遗忘因子的递归算法

引入遗忘因子后，辅助变量 \mathbf{M} 和 \mathbf{v} 的递推关系需要相应修改。新的累积信息矩阵更新为：

累积测量向量的更新为：

这两个公式的含义是：在添加新数据之前，先将所有历史累积量乘以 \lambda 进行折扣，然后再加上当前时刻的贡献。使用与 \lambda = 1 时完全相同的推导方法（应用矩阵求逆引理），我们得到带遗忘因子的完整递归算法。

校正增益的计算变为：

注意分母中出现了 \lambda 而不是1。协方差矩阵的更新为：

这里的关键是前面的因子 1/\lambda，它使得 \mathbf{P} 不再单调递减。输出预测和参数更新公式保持不变：

遗忘因子对算法性质的影响

带遗忘因子的递归最小二乘估计器有几个重要的性质变化。首先，当 \lambda = 1 时，算法退化为标准的递归最小二乘，这保证了兼容性。

其次，即使没有遗忘因子，如果 \mathbf{P} 变得奇异（接近不可逆），引入遗忘因子后 \mathbf{P} 仍然会保持奇异性。这是因为秩亏的问题源于回归向量的共线性，遗忘因子并不能改变这种结构性问题。

第三，当 \lambda < 1 时，\mathbf{P} 不再必然单调递减。具体来说，如果某个时刻回归向量为零（\mathbf{r}(t) = \mathbf{0}），此时没有新信息到来，但由于更新公式：

协方差矩阵会增大，反映了缺乏新数据时不确定性的增加。如果 \mathbf{r}(t) 长时间保持为零或接近零，\mathbf{P}(t) 可能会指数增长，这种现象称为"爆炸"。这在实际应用中可能导致数值不稳定。

第四，当 \lambda < 1 时，\mathbf{P}(t+1) 和 \mathbf{k}(t+1) 都不会收敛到零，而是会持续波动。这种持续的波动使得估计器能够对参数变化保持敏感，但代价是估计值 \hat{\mathbf{p}}_{\text{LS}} 的波动也会增大。遗忘因子越小（即 \lambda 远离1），增益 \mathbf{k} 越大，这意味着算法对参数变化的跟踪性能更好，但估计的波动性也更大。这是参数跟踪中响应速度与估计方差之间的经典权衡。

参数突变情况的处理

当系统参数发生突变（例如系统模式切换或故障）时，指数遗忘因子的表现并不理想。主要问题有两个：第一，估计值的行为可能变得不稳定，因为卡尔曼增益不会收敛到零，导致估计在突变前后震荡。第二，如果参数突变幅度较大，为了保证数值稳定性，需要将 \lambda 设置得接近1，但这会使得适应速度非常慢，无法快速捕捉突变。

对于参数突变场景，通常采用另一种策略：使用不带遗忘因子的标准递归最小二乘（\lambda = 1），但持续监测预测误差 e(t+1) = y(t+1) - y_m(t+1)。预测误差的方差在统计框架下可以表示为：

其中 \sigma^2 是测量噪声的方差。当观测到的预测误差远超过其预期方差时（例如通过卡方检验），这可能表明参数发生了突变。此时的应对策略是人工增大协方差矩阵 \mathbf{P}，例如将其重置为 \mathbf{P}(t) \leftarrow \alpha \mathbf{P}(t)，其中 \alpha \gg 1 是一个放大因子。增大 \mathbf{P} 的效果是增大卡尔曼增益 \mathbf{k}，使得算法暂时对新数据赋予更大的权重，从而快速适应参数突变。这种方法比固定使用小的 \lambda 更灵活，因为它只在检测到异常时才介入，正常运行时保持较小的估计方差。

卡尔曼滤波器

状态估计问题的引入

现在我们从参数估计转向状态估计问题。考虑一个由离散时间状态空间方程描述的系统：

这里 \mathbf{x}(t) 是系统的内部状态向量，\mathbf{u}(t) 是已知的输入，\mathbf{y}(t) 是可测量的输出。矩阵 \mathbf{A}^*、\mathbf{B}^* 和 \mathbf{C}^* 是系统的真实参数。我们的目标是从系统的输入-输出行为中估计状态 \mathbf{x}(t)。

一个直观的想法是：既然系统矩阵 \mathbf{A}^*、\mathbf{B}^* 和 \mathbf{C}^* 已知，我们可以构建一个状态模型，它由相同的输入 \mathbf{u}(t) 驱动。PPT中的框图展示了这个结构：模型接收与真实系统相同的输入，通过矩阵 \mathbf{A}、\mathbf{B} 和 \mathbf{C} 进行状态更新和输出预测。如果模型的初始状态与真实系统一致，并且模型参数准确，理论上模型状态 \mathbf{x}_m(t) 应该跟踪真实状态 \mathbf{x}(t)。

但这个简单方案有严重缺陷：在没有校正机制的情况下，即使模型输出 y_m(t) 与真实输出 y(t) 出现很大偏差，模型也不会自我纠正。初始状态误差、模型参数的微小不准确、或者未建模的扰动都会导致估计逐渐发散。

Kalman-Luenberger观测器

为了解决这个问题，Kalman和Luenberger提出了一个关键改进：利用输出误差 \mathbf{y}(t) - \mathbf{y}_m(t) 来校正模型的状态估计。具体来说，引入一个增益矩阵 \mathbf{K}，将输出预测误差反馈到状态更新中。修改后的框图显示，误差信号 e(t) = \mathbf{y}(t) - \mathbf{y}_m(t) 经过增益 \mathbf{K} 后被注入到状态更新环节。这种反馈机制使得当模型输出偏离真实输出时，状态估计会被推向正确的方向。

增益矩阵 \mathbf{K} 的选择决定了状态估计收敛到真实状态的速度。在确定性系统（无噪声）且模型完全准确的情况下，适当选择 \mathbf{K} 可以使估计误差按几何级数收敛。然而，这种确定性设计方法没有明确考虑噪声的影响。

卡尔曼滤波器采用了类似的结构，但其核心贡献是显式地对噪声进行建模，并通过统计最优准则自动确定增益矩阵。这使得卡尔曼滤波器能够在噪声环境中给出最优（在最小均方误差意义下）的状态估计。

卡尔曼滤波器的随机模型

卡尔曼滤波器的优势在于它能够自然地处理非平稳（时变）系统。考虑带有噪声的离散时间状态空间模型，状态方程为：

观测方程为：

这个模型与之前的确定性模型的关键区别在于引入了两类噪声。第一类是状态扰动 \mathbf{v}_t，它代表了建模误差、执行器的不完美、或过程噪声。状态扰动的存在意味着即使我们完全知道系统矩阵和初始状态，也无法精确预测未来的状态演化。第二类是测量噪声 \mathbf{w}_t，它代表传感器噪声，使得我们无法精确观测真实状态。

问题设定是：给定系统矩阵 \mathbf{A}_t、\mathbf{B}_t 和 \mathbf{C}_t 已知，以及输入 \mathbf{u} 和输出 \mathbf{y} 的历史观测值，估计当前状态 \mathbf{x}_t。

噪声的统计假设

卡尔曼滤波器对噪声做出了一系列统计假设。首先，状态扰动和测量噪声都假设为零均值随机变量：

零均值假设意味着噪声不会引入系统性偏差。其次，这些噪声序列对应于独立同分布的随机变量，它们的二阶统计特性由协方差矩阵刻画：

这里 \delta_{tk} 是克罗内克符号，定义为：当 t = k 时 \delta_{tk} = 1，否则 \delta_{tk} = 0。第一个等式表明状态扰动与测量噪声在任何时刻都不相关。第二个等式说明状态扰动在不同时刻之间是独立的（白噪声性质），且在时刻 t 的协方差矩阵为 \mathbf{V}_t。类似地，第三个等式表明测量噪声也是时间白色的，协方差为 \mathbf{W}_t。

矩阵 \mathbf{V}_t 和 \mathbf{W}_t 假设为已知的对称正定矩阵。这些协方差矩阵量化了噪声的强度和结构——对角元素越大表示对应分量的噪声越强，非对角元素则描述不同噪声分量之间的相关性。初始状态 \mathbf{x}_0 也被建模为一个已知均值 \mathbf{m}_0 和协方差矩阵 \mathbf{X}_0 的随机向量，并且假设它与所有时刻的 \mathbf{v}_t 独立。

从递归最小二乘到卡尔曼滤波器的路径

卡尔曼滤波器方程可以通过三个逐步推广的步骤从递归最小二乘方程推导出来。这种分步推导不仅提供了算法的数学严密性，也揭示了各个组件的物理意义。第一步是将标量输出的递归最小二乘扩展到向量输出的情况。在实际系统中，我们通常同时测量多个输出变量（例如多个传感器的读数），因此这个扩展是必要的。第二步考虑静态系统（状态不随时间演化）但存在测量噪声的情况。第三步引入状态扰动，从而得到完整的动态卡尔曼滤波器。

向量测量的递归最小二乘

在时刻 t，假设我们测量了 n_y 个标量输出，将它们组合成向量 \mathbf{y}_t。对于参数线性模型，输出预测变为：

这里的回归矩阵 \mathbf{R}_{t-1}^T 现在的维度是 n_y \times n_p（n_p 是参数的个数）。与标量情况不同，每一行对应一个输出分量的回归向量。加权最小二乘代价函数推广为：

其中 \mathbf{Q}_t 是正定半定的权重矩阵，维度为 n_y \times n_y。这个权重矩阵允许我们对不同输出分量赋予不同的重要性，也可以反映测量的相对精度。

最小化这个代价函数的估计器 \hat{\mathbf{p}}_{\text{LS}} 可以通过类似标量情况的迭代求解。递归公式包括四个步骤。第一步是计算卡尔曼增益：

注意这里需要求逆的矩阵维度是 n_y \times n_y（而不是参数空间的维度 n_p \times n_p）。当输出维度远小于参数维度时，这个求逆的计算成本更低。第二步是更新估计误差协方差矩阵：

第三步是预测下一时刻的输出向量：

第四步是参数更新：

矩阵 \mathbf{P}_t 在统计框架下可以解释为估计误差的协方差矩阵：

前提条件是 \mathbf{Q}_t 被取为测量噪声协方差矩阵的逆，即 \mathbf{Q}_t = \mathbf{W}_t^{-1}。这个选择使得算法在统计意义上最优。

静态系统情况

作为向卡尔曼滤波器过渡的中间步骤，我们考虑一个特殊的静态系统，其中状态不随时间变化但存在测量噪声：

在这个设定下，"参数" \mathbf{p} 实际上就是我们要估计的常数状态 \mathbf{x}，而回归矩阵 \mathbf{R}_{t-1}^T 对应于观测矩阵 \mathbf{C}_t。对状态 \mathbf{x}_t 和测量噪声 \mathbf{w}_t 施加前面提到的统计假设，我们通过最小化代价函数来估计状态：

其中 \mathbf{Q}_t = \mathbf{W}_t^{-1}。在高斯测量噪声的假设下，这个最小二乘估计器等价于最大似然估计器——它给出了在给定观测数据下最可能的状态值。

这个静态情况虽然简化，但它揭示了一个重要原理：当噪声的统计特性已知时，我们应该用噪声协方差矩阵的逆来加权代价函数。测量越精确（协方差越小），对应的权重越大。下一步是将这个静态结果推广到状态随时间演化且存在状态扰动的动态情况，从而得到完整的卡尔曼滤波器。

静态系统的卡尔曼滤波器推导

将向量测量的递归最小二乘算法应用于静态系统，并进行变量替换。用状态估计 \hat{\mathbf{x}}_t 替代参数估计 \hat{\mathbf{p}}_{\text{LS}}(t)，用测量噪声协方差的逆 \mathbf{W}_t^{-1} 替代权重矩阵 \mathbf{Q}_t，用观测矩阵 \mathbf{C}_{t+1} 替代回归矩阵 \mathbf{R}_t^T。经过这些替换，递归算法变为：

卡尔曼增益的计算为：

估计误差协方差矩阵的更新为：

输出预测为：

状态估计的更新为：

这就是静态情况下的卡尔曼滤波器。矩阵 \mathbf{P}_t 在这个框架下有明确的统计解释——它是估计误差的协方差矩阵：

这个解释的前提是权重矩阵 \mathbf{Q}_t 被正确选择为测量噪声协方差的逆。

动态系统的卡尔曼滤波器

现在考虑完整的动态系统，其中状态随时间演化且存在状态扰动：

系统的动态特性意味着状态在两次测量之间会发生演化，这个演化由系统矩阵 \mathbf{A}_t 和输入项 \mathbf{B}_t \mathbf{u}_t 确定。状态扰动 \mathbf{v}_t 的存在为这个演化过程带来了不确定性——即使我们完美地知道当前状态和输入，也无法精确预测下一时刻的状态。

这种动态性导致我们需要区分两类估计量。第一类是先验估计（a priori estimate），记为 \hat{\mathbf{x}}_{t+1|t} 和 \mathbf{P}_{t+1|t}，它们表示在 t+1 时刻测量到达之前，仅利用到时刻 t 为止的信息对状态和协方差的预测。第二类是后验估计（a posteriori estimate），记为 \hat{\mathbf{x}}_{t+1|t+1} 和 \mathbf{P}_{t+1|t+1}，它们表示在 t+1 时刻的测量到达之后，融合了新信息的更新估计。

状态和协方差的预测步骤

预测步骤的目的是从时刻 t 的后验估计 \hat{\mathbf{x}}_{t|t} 推算出时刻 t+1 的先验估计 \hat{\mathbf{x}}_{t+1|t}。关键观察是：状态扰动 \mathbf{v}_t 的过去值不提供关于其当前值的任何信息，因为我们假设它是时间白噪声（各时刻独立）。因此，预测器通过用状态扰动的期望值（即零）替换真实的 \mathbf{v}_t 来构造：

这个预测公式本质上是在无噪声假设下，从估计状态 \hat{\mathbf{x}}_{t|t} 出发对系统动态的确定性演化。它忽略了状态扰动的随机影响，因此会产生预测误差。

预测误差定义为真实状态与预测状态之差：

这个误差包含两部分来源：第一部分 \mathbf{A}_t(\mathbf{x}_t - \hat{\mathbf{x}}_{t|t}) 是由于前一时刻的估计误差通过系统动态传播而来，第二部分 \mathbf{v}_t 是当前时刻新引入的状态扰动。

如果后验估计是无偏的（即 \mathbb{E}(\mathbf{x}_t - \hat{\mathbf{x}}_{t|t}) = \mathbf{0}），那么先验预测也是无偏的（\mathbb{E}(\mathbf{x}_{t+1} - \hat{\mathbf{x}}_{t+1|t}) = \mathbf{0}）。在无偏假设下，预测误差的协方差矩阵可以推导如下：

将预测误差的表达式代入并展开：

由于估计误差 \mathbf{x}_t - \hat{\mathbf{x}}_{t|t} 与当前时刻的状态扰动 \mathbf{v}_t 独立（前者只依赖于到时刻 t 为止的历史，而 \mathbf{v}_t 与历史独立），交叉项的期望为零。因此：

这个公式揭示了不确定性如何在时间上传播：后验协方差 \mathbf{P}_{t|t} 通过系统动态矩阵 \mathbf{A}_t 的相似变换传播到下一时刻，同时叠加了过程噪声协方差 \mathbf{V}_t 的贡献。即使 \mathbf{P}_{t|t} 很小（估计很准确），状态扰动也会使先验不确定性 \mathbf{P}_{t+1|t} 增大。

状态和协方差的更新步骤

一旦时刻 t+1 的测量 \mathbf{y}_{t+1} 到达，我们需要利用这个新信息来更新先验估计。更新的目标是从 \hat{\mathbf{x}}_{t+1|t} 和 \mathbf{P}_{t+1|t} 得到 \hat{\mathbf{x}}_{t+1|t+1} 和 \mathbf{P}_{t+1|t+1}。

关键的观察是：在测量到达的这一瞬间，真实状态 \mathbf{x}_{t+1} 本身并不发生变化（状态演化发生在测量之间，而不是测量时刻）。变化的只是我们对状态的知识——新的测量提供了额外的信息来修正估计。

因此，更新步骤可以使用静态系统情况下推导的公式，只需要将时间索引适当调整。这是因为在固定时刻 t+1，问题退化为：给定一个状态的先验分布（均值 \hat{\mathbf{x}}_{t+1|t}，协方差 \mathbf{P}_{t+1|t}）和一个新的噪声测量 \mathbf{y}_{t+1}，求后验分布。这正是静态卡尔曼滤波器处理的问题。

完整的卡尔曼滤波器迭代算法

将预测和更新步骤组合起来，完整的卡尔曼滤波器在每个时刻执行以下四个步骤。给定时刻 t 的后验估计 (\hat{\mathbf{x}}_{t|t}, \mathbf{P}_{t|t})，计算时刻 t+1 的后验估计 (\hat{\mathbf{x}}_{t+1|t+1}, \mathbf{P}_{t+1|t+1})：

第一步，预测状态向量和协方差矩阵：

第二步，计算卡尔曼增益：

第三步，更新状态估计：

第四步，更新协方差矩阵估计：

卡尔曼滤波器的重要性质

关于算法的初始化，常见的做法是设置 \mathbf{x}_{0|0} = \mathbf{0} 和 \mathbf{P}_{0|0} = c\mathbf{I}，其中 c > 0 是一个较大的常数。这个选择反映了初始时对状态的完全不确定性。

协方差矩阵的演化有一个重要性质。与递归最小二乘类似，更新步骤导致的协方差变化是非正的：

这个差值是负半定矩阵，意味着新的测量 \mathbf{y}_{t+1} 不会增加不确定性——它只会减少或至多保持不确定性不变（当测量信息量为零时）。

卡尔曼滤波器还提供了一种跟踪时变参数的替代方法。考虑特殊的静态情况：设置 \mathbf{A}_t = \mathbf{I}（单位矩阵）且 \mathbf{B}_t = \mathbf{0}（无输入项），这时状态方程退化为 \mathbf{x}_{t+1} = \mathbf{x}_t + \mathbf{v}_t。在这个框架下，状态向量扮演参数的角色，而状态扰动解释了参数的时变性。通过调整状态扰动协方差 \mathbf{V}_t 的大小，我们可以控制每个"参数"允许变化的速度——\mathbf{V}_t 越大，算法对参数变化的响应越快。这种方法比递归最小二乘的遗忘因子更灵活，因为它允许对不同参数分量设置不同的变化速率。