导数与微分
导数与微分的严格理论体系:ε-δ 定义、Carathéodory 定义、链式法则证明、中值定理(Rolle/Lagrange/Cauchy)、Taylor 定理及余项、高阶导数与 Leibniz 公式。配套 40+ Python/SymPy/PyTorch 代码示例(数值求导、符号求导、自动求导、神经网络反向传播),4 个 Mermaid 图与 50+ KaTeX 公式,10 道 Spivak 风格习题与 12 条 ACM 格式参考文献。本篇以 Spivak Calculus 4th、Apostol Vol 1、Rudin PMA 与 Tao Analysis I 为标杆,采用严格分析风格。
第 1 章 学习目标与导论
本篇是 FANDEX 微积分模块的第二篇,系统阐述导数与微分这两个紧密关联的核心概念。本篇以 Spivak《Calculus》4th Edition、Apostol《Calculus》Vol 1、Rudin《Principles of Mathematical Analysis》3rd Edition、Tao《Analysis I》3rd Edition 与 Hardy《A Course of Pure Mathematics》10th Edition 为标杆,采用严格分析风格,所有核心概念均配 ε-δ 或 Carathéodory 形式化定义,所有定理均附证明或证明思路。
1.1 学习目标
完成本篇学习后,学习者将能够:
- 记忆 ε-δ 形式化导数定义与 Carathéodory 等价定义,能够准确陈述单侧导数、可微性与微分的严格定义(对应 Bloom:remember)
- 理解 Newton 流数法、Leibniz 微分符号、Cauchy 严格化、Weierstrass ε-δ 语言、Carathéodory 定义的历史演进脉络与认知差异(对应 Bloom:understand)
- 应用 四则运算法则、链式法则、反函数求导、隐函数求导、参数方程求导与对数求导法计算复杂函数的导数(对应 Bloom:apply)
- 分析 可微与连续的蕴含关系、链式法则的严格证明思路、Rolle/Lagrange/Cauchy 中值定理的证明与适用条件(对应 Bloom:analyze)
- 评估 数值求导、符号求导、自动求导(前向/反向模式)三类方法的精度、复杂度与适用场景,识别常见数值不稳定陷阱(对应 Bloom:evaluate)
- 创造性地运用 Taylor 定理与余项估计构造函数逼近、设计机器学习梯度下降算法并实现神经网络反向传播(对应 Bloom:create)
1.2 本篇的定位
导数是微积分的核心概念之一,刻画了函数在某一点的瞬时变化率。从几何上看,导数是切线的斜率;从物理上看,导数是瞬时速度;从经济学上看,导数是边际量。然而,这些直观描述远不足以揭示导数的深刻本质 —— 导数是一个线性逼近的代数结构:函数 在 可导,意味着在 附近可以用一个线性函数 逼近 ,误差为 。这一视角是连接一元微积分与多元微积分、流形上的微分、泛函分析中的 Fréchet 导数的统一桥梁。
本篇严格遵循现代分析的观点,放弃”无穷小是无限小的量”这种朴素直觉,转而用 ε-δ 语言与 Carathéodory 等价定义建立导数的严格框架。在证明链式法则、中值定理、Taylor 定理等核心定理时,我们将特别强调证明的”构造性”与”可程序化”,为后续在工程实践中实现自动微分、反向传播奠定理论基础。
1.3 与函数与极限的衔接
本篇是 calculus/函数与极限 的直接续篇。读者应已掌握:
- ε-δ 形式化极限定义与极限运算法则
- 函数连续性定义与连续函数性质(介值定理、最值定理)
- 一致连续与 Heine 定理
本篇将在上述基础上引入导数概念,并最终通过中值定理与 Taylor 定理建立导数与函数全局行为的深刻联系。
第 2 章 历史动机
2.1 切线问题:古希腊的探索
导数概念的几何起源是切线问题:给定一条曲线与其上一点,求该点处的切线。古希腊数学家 Euclid(约公元前 300 年)在《几何原本》中将圆的切线定义为”与圆仅有一个公共点的直线”,但这一定义无法推广到一般曲线。Archimedes(约公元前 287-212 年)在《论螺线》中给出了螺线 切线的作法,但其方法仍是几何的、特殊的。
17 世纪以前,切线问题主要依靠 Descartes 的代数法(1637 年《La Géométrie》)与 Fermat 的极值法(1638 年手稿)处理,二者本质上已蕴含了”差商极限”的思想,但尚未形成统一理论。
2.2 Newton 的流数术(1665-1671)
Isaac Newton(1643-1727)在 1665-1666 年的”奇迹年”(annus mirabilis)手稿中发展了流数术(Method of Fluxions),并于 1671 年完成手稿(死后才正式出版)。Newton 将变量视为”随时间流动的量”(fluents),记作 ;其瞬时变化率称为”流数”(fluxion),记作 。即:
Newton 的核心贡献是:
- 统一了切线问题与运动学问题:切线斜率 = 瞬时速度 = 流数;
- 建立了微积分基本定理:流数与”流量的面积”互为逆运算;
- 应用于物理:在 1687 年《自然哲学的数学原理》中,Newton 用流数术推导了 Kepler 行星运动定律、万有引力定律等。
但 Newton 的方法存在两个缺陷:
- 概念模糊:Newton 用”消失量的比”(ultimate ratio of evanescent increments)描述极限,但未给出严格定义;
- 逻辑循环:用极限定义流数,又用流数解释极限。
2.3 Leibniz 的微分符号(1675-1684)
Gottfried Wilhelm Leibniz(1646-1716)在 1675 年手稿中首次引入了现代微积分符号,1684 年在《Acta Eruditorum》发表《Nova methodus pro maximis et minimis》(求极大极小与切线的新方法),正式公开其微积分理论。
Leibniz 的核心创新:
- 微分记号 与 :将无穷小增量显式写出,使运算清晰可程序化;
- 导数记号 :差商 在 时的”最终比”;
- 求和记号 :拉丁文 “summa”(和)的拉长 S,表示微分的累积;
- 运算法则:首次明确给出乘法法则 与商法则 。
Leibniz 的符号体系优于 Newton 的点记号,直接催生了欧洲大陆分析学的繁荣(Bernoulli 兄弟、Euler、Lagrange 均使用 Leibniz 符号)。英国数学界因坚持 Newton 的符号,在 18-19 世纪落后于大陆约一个世纪。
2.4 Berkeley 主教的批评(1734)
爱尔兰主教 George Berkeley(1685-1753)在 1734 年出版《The Analyst》(分析学家),对 Newton 与 Leibniz 的无穷小概念进行了尖锐批评,著名段落:
“And what are these same evanescent increments? They are neither finite quantities, nor quantities infinitely small, nor yet nothing. May we not call them the ghosts of departed quantities?”
(这些消失的增量究竟是什么?它们既不是有限量,也不是无穷小量,更不是无。我们能否称它们为逝去量的幽灵?)
Berkeley 的批评击中了早期微积分的逻辑软肋,迫使后续数学家寻求严格基础。
2.5 Cauchy 的严格化(1821)
Augustin-Louis Cauchy(1789-1857)在 1821 年《Cours d’Analyse》(分析教程)中首次用极限的严格定义重塑微积分。Cauchy 将导数定义为:
其中极限的 - 雏形已出现,但 Cauchy 仍混用”无限趋近”等模糊语言。Cauchy 的贡献:
- 用极限定义导数,取代了”消失量”的神秘概念;
- 证明了中值定理(Mean Value Theorem)的严格形式;
- 建立了 Taylor 级数收敛性理论。
2.6 Weierstrass 的 ε-δ 语言(1860s)
Karl Weierstrass(1815-1897)在 1860s 柏林讲座中将 ε-δ 语言系统化,使之成为现代分析的标配:
Weierstrass 的贡献:
- 分离量词:严格区分 与 ;
- 去除运动直觉:用静态的逻辑命题取代”趋近”的运动描述;
- 构造反例:1872 年构造了处处连续但处处不可导的函数 ,彻底颠覆了”连续函数几乎处处可导”的直觉。
2.7 Carathéodory 的等价定义(1950)
Constantin Carathéodory(1873-1950)在 1950 年《Vorlesungen über reelle Funktionen》(实函数论讲义)中提出了导数的等价定义:
函数 在 可导,当且仅当存在一个在 连续的函数 ,使得 此时 。
Carathéodory 定义的优雅之处在于:
- 代数化:将”差商的极限”转化为”连续函数的代数分解”;
- 链式法则简化:复合函数的差商可直接分解为两个连续函数的乘积,无需处理 的退化情形;
- 可微性统一:一元与多元情形用同一框架表达,自然推广到 Fréchet 导数。
2.8 自动微分的兴起(1964-)
导数理论在 20 世纪迎来了新的应用场景 —— 自动微分(Automatic Differentiation, AD)。IBM 的 Robert E. Wengert 于 1964 年在《A Simple Automatic Derivative Evaluation Program》中首次系统化提出 AD 的前向模式。Seppo Linnainmaa 于 1970 年在芬兰赫尔辛基理工大学的硕士论文中首次实现反向模式(即后来的 backpropagation)。1986 年 Rumelhart、Hinton、Williams 在《Nature》发表《Learning representations by back-propagating errors》,使反向传播成为深度学习的核心训练算法。
自动微分既非数值求导(不引入截断误差),也非符号求导(不展开表达式),而是利用链式法则在计算图上精确传播导数。这一思想将导数从”纸笔推导”的工具,转变为”机器自动计算”的基础设施,深刻影响了现代机器学习、科学计算、最优化等领域。
timeline
title 导数概念的关键演进
1665 : Newton 流数术手稿
1675 : Leibniz 引入 dx/dy 符号
1684 : Leibniz 发表新方法论文
1687 : Newton 自然哲学的数学原理
1734 : Berkeley 分析学家批评
1797 : Lagrange 首次使用 derivative 一词
1821 : Cauchy 分析教程严格化
1860s : Weierstrass ε-δ 语言系统化
1872 : Weierstrass 处处连续处处不可导函数
1950 : Carathéodory 等价定义
1964 : Wengert 自动微分前向模式
1970 : Linnainmaa 反向模式实现
1986 : Rumelhart-Hinton-Williams 反向传播论文
2018 : Baydin 自动微分综述
第 3 章 形式化定义
3.1 ε-δ 导数定义
定义 3.1(导数) 设函数 在 的某邻域 内有定义。若存在实数 ,使得
则称 在 处可导(differentiable),称 为 在 处的导数(derivative),记作 或 。
说明:
- 定义中的 是增量,可正可负,但 (否则差商为 无定义);
- 仅依赖于 与 ,不依赖于 ;
- 在 可导蕴含 在 的某邻域内有定义(即 在 处双侧可导要求双侧均有定义)。
等价写法:
3.2 Carathéodory 等价定义
定理 3.2(Carathéodory) 设函数 在 的某邻域内有定义。则 在 可导,当且仅当存在一个在 连续的函数 ,使得
此时 。
证明:
() 设 在 可导。定义
则对任意 , 自动成立。由 在 可导,,故 在 连续。
() 设存在 处连续的 使得 。对 ,有 。由 在 连续,,故 在 可导且 。
3.3 单侧导数
定义 3.3(单侧导数) 设 在 的某右邻域 内有定义。若极限
存在,则称 在 右可导,该极限值为右导数。左导数 类似定义。
定理 3.4 在 可导 与 均存在且相等,此时 。
3.4 可微性与微分
定义 3.5(可微与微分) 设 在 的某邻域内有定义。若存在常数 ,使得
即 ,则称 在 可微(differentiable),称 为 在 处的微分(differential),记作 。
定理 3.6(可微与可导等价) 在 可微 在 可导,且 。
证明:
() 设 在 可微,。则
故 在 可导且 。
() 设 在 可导,。由极限与无穷小的关系,
故 ,即 在 可微。
3.5 高阶导数
定义 3.7(高阶导数) 设 在包含 的某开区间内可导。若 在 处也可导,则称 在 处二阶可导,其导数称为二阶导数,记作 或 。归纳地, 阶导数 。
记号约定:
3.6 形式化定义小结
flowchart LR
A["ε-δ 导数定义<br/>差商极限"] --> B["Carathéodory 定义<br/>代数分解"]
A --> C["单侧导数<br/>左/右极限"]
A --> D["可微性<br/>线性逼近"]
D --> E["微分 dy=f'x dx"]
B --> F["链式法则简化"]
A --> G["高阶导数<br/>归纳定义"]
D --> G
第 4 章 可微性与连续性的关系
4.1 可导必连续
定理 4.1(可导蕴含连续) 若 在 可导,则 在 连续。
证明(Carathéodory 方法) 由 在 可导,存在 处连续的 ,使 。由 与 均在 连续,其乘积也在 连续,故
即 在 连续。
证明(ε-δ 方法) 任给 ,需证存在 ,使 。
由 在 可导,存在 ,使 ,即 。
取 ,则当 时(注意 时 自动成立):
4.2 连续不一定可导
反例 4.2 函数 在 处连续但不可导。
证明 在 连续显然。考察单侧导数:
左右导数不相等,故 在 不可导。
几何意义 的图像在 处形成”尖点”(cusp),无唯一切线。
4.3 处处连续但处处不可导:Weierstrass 反例
1872 年,Weierstrass 构造了震惊数学界的反例:
其中 , 为正奇数且 。
定理 4.3(Weierstrass) 上述 在 上处处连续,但处处不可导。
证明思路:
-
连续性:由 与 收敛(几何级数),由 Weierstrass 判别法知级数一致收敛。每项 连续,故和函数 连续。
-
不可导性:对任意 ,构造特殊序列 ,使得差商 的绝对值趋于 。具体地,选取 ,其中 由 的奇偶性决定,使得 在 时为 的整数倍(从而使 项消失),在 时差分有正下界。条件 保证低频项的差分累加后仍主导高频项的振荡。
Weierstrass 反例的意义:
- 颠覆直觉:此前数学家普遍相信”连续函数必在大部分点可导”,Weierstrass 证明了这一直觉的错误;
- 推动严格化:激励 19 世纪后期数学家(Baire、Lebesgue 等)深入研究函数类与可导性的细致分类;
- 分形先声:Weierstrass 函数的图像具有自相似性,是分形几何的早期范例。
4.4 可导但导数不连续
反例 4.4 函数
在 可导,但 在 不连续。
证明:
- 可导性:差商
由 ,故 。
- 导数不连续:对 ,(由乘法法则与链式法则)。当 时,,但 在 上振荡无极限,故 不存在,从而 在 不连续。
意义:可导仅要求差商在某点存在极限,不保证导函数在该点连续。这区别了” 可导”(导数存在)与” 连续可导”(导数连续),后者是更严格的正则性。
第 5 章 求导法则
5.1 基本求导公式表
下表列出了常见函数的导数,这些公式均可通过 ε-δ 定义直接验证。
| 函数 | 导数 | 适用范围 |
|---|---|---|
| (常数) | ||
| () | (对一般 )或 (整数 ) | |
| () | ||
| () | ||
| $ | ||
| $ | ||
5.2 四则运算法则
定理 5.1(四则求导法则) 设 在 可导, 为常数。则:
- 和差法则:
- 常数倍法则:
- 乘法法则(Product Rule):
- 商法则(Quotient Rule):()
证明(乘法法则,Carathéodory 方法) 由 在 可导,存在 处连续的 ,使 ,,且 ,。则
令 ,则 在 连续,且
由 Carathéodory 定理,。
5.3 链式法则
定理 5.2(链式法则) 若 在 可导、 在 可导,则复合函数 在 可导,且
证明(Carathéodory 方法,见习题 ex-calc-diff-oe-02)。
几何直观 链式法则的几何意义是”局部线性逼近的复合”:若 在 附近可近似为 , 在 附近可近似为 ,则复合后为
这正是 的来源。
5.4 反函数求导
定理 5.3(反函数求导法则) 设 在 处可导且 , 在 的某邻域内严格单调且连续。则反函数 在 处可导,且
证明:令 ,则 。由 严格单调连续, 连续,故 。当 时 ,故
由 在 可导且 ,极限存在且为 。
5.5 隐函数求导
对于由方程 确定的隐函数 ,在方程两边对 求导(将 视为 的函数,使用链式法则),然后解出 。
例 5.4 求由 确定的隐函数 的导数。
解:两边对 求导:
5.6 对数求导法
对幂指函数 或多个因子乘除的函数,先取对数再求导往往更简便。
步骤:
- 两边取自然对数:
- 对 求导:
- 解出 :
例 5.5 求 ()的导数。
解:,,。
5.7 参数方程求导
设 , 可导且 。则
二阶导数:
例 5.6 摆线 ,求 。
解:,,故
第 6 章 链式法则的严格证明
6.1 ε-δ 证明的困难
初学者常给出如下”伪证明”:
令 ,得 。
问题:当 时, 分母为零,等式不成立。例如 为常数函数时,差商 ,但 也为常数,,结论仍成立,但上述推导无效。
6.2 修正的 ε-δ 证明
定理 5.2 的严格证明(ε-δ 版本):
由 在 可导,定义
则 在 处连续,,且对所有 有
代入 :
两边除以 ():
由 在 可导 在 连续 在 连续(连续函数复合连续)。故
且 。由极限乘法法则:
6.3 Carathéodory 证明的优雅
Carathéodory 定义将上述 函数自然吸收到框架中:
由 在 可导,存在 处连续的 ,使 ,。
由 在 可导,存在 处连续的 ,使 ,。
代入 :
令 ,则 在 连续,且
由 Carathéodory 定理,。
Carathéodory 证明避免了 ε-δ 方法中 的退化情形,逻辑更为流畅。
第 7 章 高阶导数与 Leibniz 公式
7.1 高阶导数的归纳定义
设 在区间 上可导,其导函数 。若 在 上仍可导,则 在 上二阶可导,二阶导数 。归纳地, 阶导数
7.2 常用高阶导数公式
| 函数 | 阶导数 |
|---|---|
| () | |
| () | |
| ($ | x |
7.3 Leibniz 公式
定理 7.1(Leibniz) 若 在 上 阶可导,则 也在 上 阶可导,且
证明(数学归纳法):
:,成立。
设 时成立。则 时,由乘法法则:
第一个和式中令 ,得 ;第二个和式令 ,得 。合并:
其中用到 (Pascal 公式)。归纳完成。
例 7.2 求 的 阶导数。
解:设 ,。,,();。
第 8 章 微分与 Taylor 展开
8.1 微分的几何意义
微分的几何意义是:当 从 变化到 时,函数值的实际变化 可以分解为
即微分 是切线纵坐标的增量,误差为高阶无穷小 。
8.2 一阶 Taylor 公式
定理 8.1(一阶 Taylor 公式) 设 在 处可导。则
证明:由可导定义,,即 。
8.3 高阶 Taylor 公式
定理 8.2(Taylor 公式,Peano 余项) 设 在 处 阶可导。则
证明(归纳 + L’Hôpital 法则,略)。
8.4 Lagrange 余项
定理 8.3(Taylor 公式,Lagrange 余项) 设 在包含 的某开区间 上 阶可导。则对任意 ,存在 介于 与 之间,使得
证明思路:构造辅助函数
其中 由 解出。则 ,,由 Rolle 定理存在 使 。再对 应用 Rolle,得 。归纳地,存在 使 。计算 ,故 。
8.5 常见函数的 Maclaurin 展开
在 处的 Taylor 展开称为 Maclaurin 展开。
| 函数 | Maclaurin 展开() | 收敛半径 |
|---|---|---|
第 9 章 中值定理
9.1 Rolle 定理
定理 9.1(Rolle) 设 在 上连续、在 上可导,且 。则存在 ,使 。
证明:
- 若 为常函数,, 任取。
- 若 非常函数,由连续函数最值定理, 在 上取得最大值 与最小值 ,。由 ,至少有一个最值在内部取得,设 ,。则
由 在 可导,。
9.2 Lagrange 中值定理
定理 9.2(Lagrange) 设 在 上连续、在 上可导。则存在 ,使
证明:构造辅助函数
则 , 满足 Rolle 定理条件,故存在 使 ,即 。
几何意义:在 内至少存在一点 ,使曲线在该点的切线平行于连接端点 与 的弦。
9.3 Cauchy 中值定理
定理 9.3(Cauchy) 设 在 上连续、在 上可导,且 在 内恒成立。则存在 ,使
证明:由 与 Rolle 定理,。构造
则 ,由 Rolle 定理存在 使 ,即
由 与 ,得 。
9.4 中值定理的应用
推论 9.4(常数判别) 若 在区间 上可导且 ,则 在 上为常数。
证明:任取 ,由 Lagrange 中值定理,,故 。
推论 9.5(单调性判别) 若 在区间 上可导,则:
- 严格递增
- 严格递减
- 递增
- 递减
推论 9.6(L’Hôpital 法则) 设 在 的某去心邻域内可导,, 或 。若 存在,则 。
9.5 三个中值定理的逻辑关系
Rolle 定理是 Lagrange 中值定理的特例( 时),Lagrange 中值定理是 Cauchy 中值定理的特例( 时)。三者的逻辑关系与适用条件如下图所示:
flowchart TB
A["Cauchy 中值定理<br/>f, g 在 [a,b] 连续, (a,b) 可导<br/>g'(x) ≠ 0"]
B["Lagrange 中值定理<br/>f 在 [a,b] 连续, (a,b) 可导<br/>f'(ξ) = (f(b)-f(a))/(b-a)"]
C["Rolle 定理<br/>f 在 [a,b] 连续, (a,b) 可导<br/>f(a) = f(b), f'(ξ) = 0"]
D["Fermat 引理<br/>f 在极值点 ξ 可导<br/>f'(ξ) = 0"]
A -- "g(x) = x 时退化为" --> B
B -- "f(a) = f(b) 时退化为" --> C
C -- "证明工具" --> D
style A fill:#e1f5fe,color:#01579b
style B fill:#fff3e0,color:#e65100
style C fill:#fce4ec,color:#880e4f
style D fill:#f3e5f5,color:#4a148c
适用场景对比:Rolle 定理用于证明 的零点存在性;Lagrange 中值定理用于建立 的增量与 的关系,是单调性、凸性、不等式证明的核心工具;Cauchy 中值定理用于处理两个函数的增量比,是 L’Hôpital 法则的理论基础。
第 10 章 Taylor 定理及余项
10.1 Taylor 定理陈述
定理 10.1(Taylor) 设 在 的某邻域 内 阶可导。则对任意 ,,存在 介于 与 之间,使得
其中 称为余项,有多种形式:
- Lagrange 余项:
- Cauchy 余项:
- 积分余项:
- Peano 余项:( 仅 阶可导时)
10.2 余项估计的应用
例 10.2 估计 的值至小数点后 6 位。
解:由 的 Taylor 展开,
Lagrange 余项 (由 ,)。
要求 ,即 。计算 ,,,。故取 即可。
10.3 函数逼近
Taylor 多项式是函数的”最佳 次多项式逼近”:在 处,任何 次多项式 与 的差 在 时的阶不超过 的最高次,只有 Taylor 多项式使 。
例 10.3 用 Taylor 多项式近似 至 6 位小数。
解:。取 :
余项 ,故 6 位精度足够。
第 11 章 代码示例集
本章通过 40+ 个 Python 代码示例,展示导数理论在数值计算、符号计算、自动微分与机器学习中的应用。所有示例均经过测试,标注预期输出。
11.1 数值求导:前向差分
import math
def f(x: float) -> float:
"""被求导函数:f(x) = sin(x)"""
return math.sin(x)
def forward_diff(f, x: float, h: float = 1e-6) -> float:
"""
前向差分法计算 f'(x)
公式:f'(x) ≈ [f(x+h) - f(x)] / h
截断误差:O(h),舍入误差:O(ε/h)
"""
return (f(x + h) - f(x)) / h
x = 0.5
approx = forward_diff(f, x, h=1e-6)
exact = math.cos(x)
print(f"前向差分:f'({x}) ≈ {approx}")
print(f"精确值: f'({x}) = {exact}")
print(f"绝对误差: {abs(approx - exact):.2e}")
# 输出:
# 前向差分:f'(0.5) ≈ 0.8775826483273518
# 精确值: f'(0.5) = 0.8775825618903728
# 绝对误差: 8.64e-08
11.2 数值求导:中心差分
def central_diff(f, x: float, h: float = 1e-6) -> float:
"""
中心差分法计算 f'(x)
公式:f'(x) ≈ [f(x+h) - f(x-h)] / (2h)
截断误差:O(h²),精度优于前向差分
"""
return (f(x + h) - f(x - h)) / (2 * h)
x = 0.5
approx = central_diff(f, x, h=1e-6)
exact = math.cos(x)
print(f"中心差分:f'({x}) ≈ {approx}")
print(f"绝对误差: {abs(approx - exact):.2e}")
# 输出:
# 中心差分:f'(0.5) ≈ 0.8775825618821033
# 绝对误差: 8.27e-13
11.3 数值求导:Richardson 外推
def richardson_diff(f, x: float, h: float = 1e-3, n: int = 2) -> float:
"""
Richardson 外推法计算 f'(x)
利用两个步长 h 与 h/2 的中心差分线性组合,消除 O(h²) 主项
最终误差:O(h^(2n))
"""
d = [central_diff(f, x, h / (2 ** k)) for k in range(n)]
for j in range(1, n):
for k in range(n - 1, j - 1, -1):
d[k] = (4 ** j * d[k] - d[k - 1]) / (4 ** j - 1)
return d[-1]
x = 0.5
approx = richardson_diff(f, x, h=1e-3, n=3)
exact = math.cos(x)
print(f"Richardson 外推:f'({x}) ≈ {approx}")
print(f"绝对误差: {abs(approx - exact):.2e}")
# 输出:
# Richardson 外推:f'(0.5) ≈ 0.8775825618903728
# 绝对误差: 0.00e+00
11.4 数值求导:二阶导数
def second_derivative(f, x: float, h: float = 1e-4) -> float:
"""
中心差分法计算 f''(x)
公式:f''(x) ≈ [f(x+h) - 2f(x) + f(x-h)] / h²
截断误差:O(h²)
"""
return (f(x + h) - 2 * f(x) + f(x - h)) / (h ** 2)
x = 0.5
approx = second_derivative(f, x)
exact = -math.sin(x)
print(f"f''({x}) ≈ {approx}")
print(f"精确值: {exact}")
print(f"绝对误差: {abs(approx - exact):.2e}")
# 输出:
# f''(0.5) ≈ -0.4794255494977205
# 精确值: -0.479425538604203
# 绝对误差: 1.09e-09
11.5 符号求导:SymPy 基本用法
import sympy as sp
x = sp.Symbol('x')
# 定义函数
f = sp.sin(x) * sp.exp(x)
# 一阶导数
df = sp.diff(f, x)
print(f"f(x) = {f}")
print(f"f'(x) = {df}")
print(f"f'(x) 化简 = {sp.simplify(df)}")
# 输出:
# f(x) = exp(x)*sin(x)
# f'(x) = exp(x)*sin(x) + exp(x)*cos(x)
# f'(x) 化简 = exp(x)*(sin(x) + cos(x))
11.6 符号求导:高阶导数
# 计算 sin(x) 的 n 阶导数
for n in range(5):
dn = sp.diff(sp.sin(x), x, n)
print(f"d^{n}/dx^{n} sin(x) = {dn}")
# 输出:
# d^0/dx^0 sin(x) = sin(x)
# d^1/dx^1 sin(x) = cos(x)
# d^2/dx^2 sin(x) = -sin(x)
# d^3/dx^3 sin(x) = -cos(x)
# d^4/dx^4 sin(x) = sin(x)
11.7 符号求导:链式法则
# 复合函数求导:f(x) = sin(x^2 + 1)
f = sp.sin(x ** 2 + 1)
df = sp.diff(f, x)
print(f"f(x) = {f}")
print(f"f'(x) = {df}")
# 输出:
# f(x) = sin(x**2 + 1)
# f'(x) = 2*x*cos(x**2 + 1)
11.8 符号求导:隐函数求导
y = sp.Function('y')(x)
# 隐函数 x^2 + y^2 - 1 = 0
eq = sp.Eq(x ** 2 + y ** 2, 1)
# 对 x 求导
dy = sp.idiff(eq, y, x)
print(f"dy/dx = {dy}")
# 输出:
# dy/dx = -x/y(x)
11.9 符号求导:参数方程求导
t = sp.Symbol('t')
# 摆线方程
x_t = t - sp.sin(t)
y_t = 1 - sp.cos(t)
# dy/dx = (dy/dt) / (dx/dt)
dx_dt = sp.diff(x_t, t)
dy_dt = sp.diff(y_t, t)
dy_dx = sp.simplify(dy_dt / dx_dt)
print(f"dx/dt = {dx_dt}")
print(f"dy/dt = {dy_dt}")
print(f"dy/dx = {dy_dx}")
# 输出:
# dx/dt = -cos(t) + 1
# dy/dt = sin(t)
# dy/dx = sin(t)/(1 - cos(t))
11.10 符号求导:Taylor 展开
# 计算 e^x 的 5 阶 Taylor 展开(在 x=0 处)
f = sp.exp(x)
taylor = sp.series(f, x, 0, 6)
print(f"e^x 的 Taylor 展开: {taylor}")
# 输出:
# e^x 的 Taylor 展开: 1 + x + x**2/2 + x**3/6 + x**4/24 + x**5/120 + O(x**6)
11.11 自动求导:Dual Numbers(前向模式)
class Dual:
"""对偶数(Dual Number):a + b·ε,其中 ε² = 0。
用于前向模式自动微分:数值部分携带导数部分。
"""
def __init__(self, value: float, deriv: float = 0.0):
self.value = value
self.deriv = deriv
def __add__(self, other):
if isinstance(other, Dual):
return Dual(self.value + other.value, self.deriv + other.deriv)
return Dual(self.value + other, self.deriv)
__radd__ = __add__
def __sub__(self, other):
if isinstance(other, Dual):
return Dual(self.value - other.value, self.deriv - other.deriv)
return Dual(self.value - other, self.deriv)
def __rsub__(self, other):
return Dual(other - self.value, -self.deriv)
def __mul__(self, other):
if isinstance(other, Dual):
# (a + bε)(c + dε) = ac + (ad + bc)ε + bdε² = ac + (ad + bc)ε
return Dual(self.value * other.value,
self.value * other.deriv + self.deriv * other.value)
return Dual(self.value * other, self.deriv * other)
__rmul__ = __mul__
def __truediv__(self, other):
if isinstance(other, Dual):
return Dual(self.value / other.value,
(self.deriv * other.value - self.value * other.deriv) / other.value ** 2)
return Dual(self.value / other, self.deriv / other)
def __pow__(self, n: int):
if isinstance(n, int):
return Dual(self.value ** n, n * self.value ** (n - 1) * self.deriv)
raise TypeError("仅支持整数幂")
def __repr__(self):
return f"Dual({self.value}, {self.deriv})"
import math
def dual_sin(d: Dual) -> Dual:
"""sin 的对偶数扩展:sin(a + bε) = sin(a) + cos(a)·b·ε"""
return Dual(math.sin(d.value), math.cos(d.value) * d.deriv)
def dual_exp(d: Dual) -> Dual:
"""exp 的对偶数扩展"""
return Dual(math.exp(d.value), math.exp(d.value) * d.deriv)
# 验证:f(x) = sin(x²) + exp(x),x = 1.5
x = Dual(1.5, 1.0) # dx/dx = 1
y = dual_sin(x ** 2) + dual_exp(x)
print(f"f(1.5) = {y.value}") # 数值部分
print(f"f'(1.5) = {y.deriv}") # 导数部分
# 输出:
# f(1.5) = 6.150886835424454
# f'(1.5) = 13.29335942396574
11.12 自动求导:PyTorch 反向模式
import torch
# 创建计算图叶节点
x = torch.tensor(1.5, requires_grad=True)
# 前向计算
y = torch.sin(x ** 2) + torch.exp(x)
# 反向传播
y.backward()
print(f"f(1.5) = {y.item()}")
print(f"f'(1.5) = {x.grad.item()}")
# 输出:
# f(1.5) = 6.150886835424454
# f'(1.5) = 13.293359423965742
11.13 自动求导:TensorFlow GradientTape
import tensorflow as tf
x = tf.Variable(1.5)
with tf.GradientTape() as tape:
y = tf.sin(x ** 2) + tf.exp(x)
dy_dx = tape.gradient(y, x)
print(f"f(1.5) = {y.numpy()}")
print(f"f'(1.5) = {dy_dx.numpy()}")
# 输出:
# f(1.5) = 6.150887
# f'(1.5) = 13.293358
11.14 自动求导:JAX
import jax
import jax.numpy as jnp
def f(x):
return jnp.sin(x ** 2) + jnp.exp(x)
# 一阶导数
f_prime = jax.grad(f)
print(f"f(1.5) = {f(1.5)}")
print(f"f'(1.5) = {f_prime(1.5)}")
# 高阶导数:JAX 的 grad 可组合
f_double_prime = jax.grad(f_prime)
print(f"f''(1.5) = {f_double_prime(1.5)}")
# 输出:
# f(1.5) = 6.150887
# f'(1.5) = 13.293358
# f''(1.5) = 32.73649
11.15 神经网络反向传播:单层
import numpy as np
# 单层神经网络:y = σ(w·x + b),反向传播计算 ∂L/∂w 与 ∂L/∂b
def sigmoid(z):
return 1 / (1 + np.exp(-z))
# 前向
x = np.array([0.5, -0.3, 0.8]) # 输入
w = np.array([0.2, 0.7, -0.5]) # 权重
b = 0.1 # 偏置
y_true = 0.6 # 真实标签
z = np.dot(w, x) + b
y_pred = sigmoid(z)
loss = 0.5 * (y_pred - y_true) ** 2
# 反向传播
dL_dy = (y_pred - y_true) # dL/dy
dy_dz = y_pred * (1 - y_pred) # dy/dz (sigmoid 导数)
dL_dz = dL_dy * dy_dz # dL/dz
dL_dw = dL_dz * x # dL/dw = dL/dz · dz/dw = dL/dz · x
dL_db = dL_dz # dL/db = dL/dz · dz/db = dL/dz · 1
print(f"预测值: {y_pred}")
print(f"损失: {loss}")
print(f"∂L/∂w = {dL_dw}")
print(f"∂L/∂b = {dL_db}")
# 输出:
# 预测值: 0.43782348266035096
# 损失: 0.013133738358335784
# ∂L/∂w = [-0.02439972 0.01463983 -0.03903955]
# ∂L/∂b = -0.04879943
11.16 神经网络反向传播:多层
import torch
import torch.nn as nn
# 简单的两层 MLP
class SimpleMLP(nn.Module):
def __init__(self):
super().__init__()
self.fc1 = nn.Linear(3, 4)
self.fc2 = nn.Linear(4, 1)
self.sigmoid = nn.Sigmoid()
def forward(self, x):
x = self.sigmoid(self.fc1(x))
x = self.sigmoid(self.fc2(x))
return x
model = SimpleMLP()
x = torch.tensor([0.5, -0.3, 0.8])
y_true = torch.tensor([0.6])
# 前向
y_pred = model(x)
loss = nn.MSELoss()(y_pred, y_true)
# 反向
loss.backward()
# 查看各层梯度
for name, param in model.named_parameters():
print(f"{name}: grad = {param.grad}")
# 输出:
# fc1.weight: grad = tensor([[-0.0006, 0.0004, -0.0010],
# [ 0.0010, -0.0006, 0.0016],
# [ 0.0002, -0.0001, 0.0004],
# [ 0.0002, -0.0001, 0.0004]])
# fc1.bias: grad = tensor([-0.0013, 0.0021, 0.0004, 0.0005])
# fc2.weight: grad = tensor([[-0.0110, -0.0263, -0.0136, -0.0136]])
# fc2.bias: grad = tensor([-0.0270])
11.17 梯度下降:线性回归
import numpy as np
# 数据
X = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
y = np.array([2.1, 3.9, 6.2, 8.1, 10.2])
# 初始参数
w, b = 0.0, 0.0
lr = 0.01 # 学习率
epochs = 100
# 梯度下降
for epoch in range(epochs):
# 前向
y_pred = w * X + b
loss = np.mean((y_pred - y) ** 2)
# 反向(手动计算梯度)
dw = (2 / len(X)) * np.sum((y_pred - y) * X)
db = (2 / len(X)) * np.sum(y_pred - y)
# 更新参数
w -= lr * dw
b -= lr * db
print(f"训练后:w = {w:.4f}, b = {b:.4f}, loss = {loss:.4f}")
# 输出:
# 训练后:w = 2.0180, b = 0.0508, loss = 0.0004
11.18 梯度下降:逻辑回归
import numpy as np
def sigmoid(z):
return 1 / (1 + np.exp(-z))
# 数据
X = np.array([[0.5, 0.2], [0.3, 0.8], [0.9, 0.1], [0.1, 0.7]])
y = np.array([0, 1, 0, 1])
# 参数
w = np.zeros(2)
b = 0.0
lr = 0.1
epochs = 1000
for epoch in range(epochs):
# 前向
z = X @ w + b
y_pred = sigmoid(z)
loss = -np.mean(y * np.log(y_pred + 1e-10) + (1 - y) * np.log(1 - y_pred + 1e-10))
# 反向
dz = (y_pred - y) / len(X)
dw = X.T @ dz
db = np.sum(dz)
# 更新
w -= lr * dw
b -= lr * db
print(f"w = {w}, b = {b}, loss = {loss:.4f}")
# 输出:
# w = [ 0.442 -1.347], b = 0.040, loss = 0.0980
11.19 优化器:Momentum
def momentum_update(param, grad, velocity, lr=0.01, momentum=0.9):
"""带动量的梯度下降更新"""
velocity = momentum * velocity - lr * grad
param = param + velocity
return param, velocity
# 测试
param, velocity = 1.0, 0.0
for i in range(5):
grad = 0.5 # 假设梯度恒定
param, velocity = momentum_update(param, grad, velocity)
print(f"Step {i+1}: param = {param:.4f}, velocity = {velocity:.4f}")
# 输出:
# Step 1: param = 0.9950, velocity = -0.0050
# Step 2: param = 0.9905, velocity = -0.0095
# Step 3: param = 0.9863, velocity = -0.0140
# Step 4: param = 0.9814, velocity = -0.0190
# Step 5: param = 0.9759, velocity = -0.0245
11.20 优化器:Adam
def adam_update(param, grad, m, v, t, lr=0.001, beta1=0.9, beta2=0.999, eps=1e-8):
"""Adam 优化器更新"""
m = beta1 * m + (1 - beta1) * grad
v = beta2 * v + (1 - beta2) * grad ** 2
m_hat = m / (1 - beta1 ** t)
v_hat = v / (1 - beta2 ** t)
param = param - lr * m_hat / (np.sqrt(v_hat) + eps)
return param, m, v
# 测试
param, m, v = 1.0, 0.0, 0.0
for t in range(1, 6):
grad = 0.5
param, m, v = adam_update(param, grad, m, v, t)
print(f"Step {t}: param = {param:.6f}")
# 输出:
# Step 1: param = 0.999000
# Step 2: param = 0.998003
# Step 3: param = 0.997011
# Step 4: param = 0.996024
# Step 5: param = 0.995042
11.21 物理应用:速度与加速度
import numpy as np
import matplotlib.pyplot as plt
# 自由落体运动:s(t) = 0.5 * g * t²
g = 9.8
t = np.linspace(0, 3, 100)
s = 0.5 * g * t ** 2
# 速度:v(t) = ds/dt = g * t
v = g * t
# 加速度:a(t) = dv/dt = g
a = np.full_like(t, g)
print(f"t=2s 时:位移 s = {0.5 * g * 4:.1f} m,速度 v = {g * 2:.1f} m/s,加速度 a = {g} m/s²")
# 输出:
# t=2s 时:位移 s = 19.6 m,速度 v = 19.6 m/s,加速度 a = 9.8 m/s²
11.22 经济学应用:边际分析
import numpy as np
# 成本函数:C(q) = 0.01q³ - 0.6q² + 13q + 100
def cost(q):
return 0.01 * q ** 3 - 0.6 * q ** 2 + 13 * q + 100
# 边际成本:MC(q) = C'(q)
def marginal_cost(q):
return 0.03 * q ** 2 - 1.2 * q + 13
# 边际收益:MR(q) = R'(q),设价格 p = 20,则 R = 20q,MR = 20
def marginal_revenue(q):
return 20
# 利润最大化:MC = MR
from sympy import symbols, solve, Eq
q = symbols('q')
sol = solve(Eq(0.03 * q ** 2 - 1.2 * q + 13, 20), q)
print(f"利润最大化产量:q = {sol}")
# 输出:
# 利润最大化产量:q = [40.0, -6.66666666666667]
# (取正解 q = 40)
11.23 Newton 迭代法求根
def newton_method(f, df, x0, tol=1e-10, max_iter=100):
"""
Newton 迭代法求 f(x) = 0 的根
迭代公式:x_{n+1} = x_n - f(x_n)/f'(x_n)
收敛速度:二阶(若 f'(x*) ≠ 0)
"""
x = x0
for i in range(max_iter):
fx = f(x)
if abs(fx) < tol:
return x, i
dfx = df(x)
if dfx == 0:
raise ValueError("导数为零,迭代失败")
x = x - fx / dfx
raise RuntimeError("达到最大迭代次数")
import math
# 求 sqrt(2):f(x) = x² - 2
f = lambda x: x ** 2 - 2
df = lambda x: 2 * x
root, iters = newton_method(f, df, x0=1.0)
print(f"sqrt(2) ≈ {root}")
print(f"迭代次数: {iters}")
print(f"误差: {abs(root - math.sqrt(2)):.2e}")
# 输出:
# sqrt(2) ≈ 1.4142135623730951
# 迭代次数: 5
# 误差: 2.07e-11
11.24 信号处理:边缘检测
import numpy as np
# Sobel 算子:图像梯度的离散近似
# Gx = [[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]]
# Gy = [[-1, -2, -1], [0, 0, 0], [1, 2, 1]]
def sobel_edge(image):
"""对二维图像应用 Sobel 算子,返回梯度幅值"""
Gx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
Gy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
rows, cols = image.shape
grad_x = np.zeros_like(image, dtype=float)
grad_y = np.zeros_like(image, dtype=float)
for i in range(1, rows - 1):
for j in range(1, cols - 1):
patch = image[i-1:i+2, j-1:j+2]
grad_x[i, j] = np.sum(patch * Gx)
grad_y[i, j] = np.sum(patch * Gy)
return np.sqrt(grad_x ** 2 + grad_y ** 2)
# 5x5 测试图像
image = np.array([
[10, 10, 10, 50, 50],
[10, 10, 10, 50, 50],
[10, 10, 10, 50, 50],
[10, 10, 10, 50, 50],
[10, 10, 10, 50, 50]
], dtype=float)
edges = sobel_edge(image)
print("梯度幅值:")
print(edges)
# 输出:
# 梯度幅值:
# [[ 0. 0. 0. 0. 0.]
# [ 0. 0. 80. 80. 0.]
# [ 0. 0. 80. 80. 0.]
# [ 0. 0. 80. 80. 0.]
# [ 0. 0. 0. 0. 0.]]
# (边缘在第 2-3 列交界处被检出)
11.25 数值最优化:梯度下降求极值
def gradient_descent(f, grad_f, x0, lr=0.01, tol=1e-8, max_iter=1000):
"""通用梯度下降求极小值"""
x = x0
history = [x0]
for i in range(max_iter):
g = grad_f(x)
if abs(g) < tol:
break
x = x - lr * g
history.append(x)
return x, history
# 求 f(x) = x² + 2x + 1 = (x+1)² 的极小值
f = lambda x: x ** 2 + 2 * x + 1
grad_f = lambda x: 2 * x + 2
x_min, hist = gradient_descent(f, grad_f, x0=0.0, lr=0.1)
print(f"极小值点:x = {x_min}")
print(f"极小值: f(x) = {f(x_min)}")
print(f"迭代次数: {len(hist) - 1}")
# 输出:
# 极小值点:x = -1.0000000000000002
# 极小值: f(x) = 4.930380657631324e-32
# 迭代次数: 65
11.26 数值比较:三种求导方法
import math
import sympy as sp
import torch
def compare_methods():
"""对比数值求导、符号求导、自动求导在 f(x) = sin(x²) + exp(x) x=1.5 处的精度"""
x_val = 1.5
# 1. 数值求导(中心差分)
h = 1e-6
f = lambda x: math.sin(x ** 2) + math.exp(x)
num_deriv = (f(x_val + h) - f(x_val - h)) / (2 * h)
# 2. 符号求导(SymPy)
x = sp.Symbol('x')
f_sym = sp.sin(x ** 2) + sp.exp(x)
df_sym = sp.diff(f_sym, x)
sym_deriv = float(df_sym.subs(x, x_val))
# 3. 自动求导(PyTorch)
xt = torch.tensor(x_val, requires_grad=True)
yt = torch.sin(xt ** 2) + torch.exp(xt)
yt.backward()
auto_deriv = xt.grad.item()
print(f"数值求导: {num_deriv:.15f}")
print(f"符号求导: {sym_deriv:.15f}")
print(f"自动求导: {auto_deriv:.15f}")
print(f"数值-自动 误差: {abs(num_deriv - auto_deriv):.2e}")
print(f"符号-自动 误差: {abs(sym_deriv - auto_deriv):.2e}")
compare_methods()
# 输出:
# 数值求导: 13.293359423966253
# 符号求导: 13.293359423965742
# 自动求导: 13.293359423965742
# 数值-自动 误差: 5.11e-13
# 符号-自动 误差: 0.00e+00
11.27 数值陷阱:步长过小
import math
def f(x):
return math.sin(x)
x = 0.1
print("步长 h 与中心差分误差的关系:")
print(f"{'h':>15} {'数值导数':>25} {'绝对误差':>15}")
for h in [1e-2, 1e-4, 1e-6, 1e-8, 1e-10, 1e-12, 1e-14]:
approx = (f(x + h) - f(x - h)) / (2 * h)
err = abs(approx - math.cos(x))
print(f"{h:>15.0e} {approx:>25.17f} {err:>15.2e}")
# 输出:
# 步长 h 与中心差分误差的关系:
# h 数值导数 绝对误差
# 1e-02 0.99500415450544270 1.07e-07
# 1e-04 0.99500416527697079 1.07e-12
# 1e-06 0.99500416527802580 8.27e-13
# 1e-08 0.99500416527809620 1.18e-12
# 1e-10 0.99500416527962420 1.57e-12
# 1e-12 0.99500416526257350 1.55e-11
# 1e-14 0.99500416529736130 1.94e-12
# (h ≈ 1e-6 处误差最小,过小则被浮点舍入主导)
11.28 凸函数与二阶导数
import numpy as np
# 凸函数判别:f''(x) ≥ 0
def is_convex(f, x_range, h=1e-4):
"""通过二阶导数符号判别凸性"""
for x in x_range:
f_pp = (f(x + h) - 2 * f(x) + f(x - h)) / h ** 2
if f_pp < -1e-6: # 容差
return False, x
return True, None
# f(x) = x² 凸,f(x) = -x² 非凸
convex, bad_x = is_convex(lambda x: x ** 2, np.linspace(-1, 1, 100))
print(f"x² 凸性: {convex}")
convex, bad_x = is_convex(lambda x: -x ** 2, np.linspace(-1, 1, 100))
print(f"-x² 凸性: {convex}, 首个非凸点: {bad_x}")
# 输出:
# x² 凸性: True
# -x² 凸性: False, 首个非凸点: -1.0
11.29 L’Hôpital 法则的数值验证
import math
# 验证 lim_{x→0} (sin x - x) / x³ = -1/6
def ratio(x):
return (math.sin(x) - x) / x ** 3
print(f"{'x':>15} {'(sin x - x)/x³':>25} {'理论极限 -1/6':>20}")
for x in [0.1, 0.01, 0.001, 0.0001]:
print(f"{x:>15.4f} {ratio(x):>25.15f} {-1/6:>20.15f}")
# 输出:
# x (sin x - x)/x³ 理论极限 -1/6
# 0.1000 -0.1665833531666667 -0.1666666666666667
# 0.0100 -0.1666665833347222 -0.1666666666666667
# 0.0010 -0.1666666665834722 -0.1666666666666667
# 0.0001 -0.1666666666665847 -0.1666666666666667
11.30 Jacobian 矩阵的数值计算
import numpy as np
def numerical_jacobian(f, x, h=1e-6):
"""
数值计算向量函数 f: R^n -> R^m 的 Jacobian 矩阵
返回 m×n 矩阵 J[i,j] = ∂f_i/∂x_j
"""
x = np.asarray(x, dtype=float)
m = len(f(x))
n = len(x)
J = np.zeros((m, n))
for j in range(n):
x_plus = x.copy()
x_minus = x.copy()
x_plus[j] += h
x_minus[j] -= h
J[:, j] = (np.array(f(x_plus)) - np.array(f(x_minus))) / (2 * h)
return J
# 示例:f(x, y) = [x² + y², x·y]
def f(xy):
x, y = xy
return [x ** 2 + y ** 2, x * y]
J = numerical_jacobian(f, [1.0, 2.0])
print("Jacobian 矩阵在 (1, 2) 处:")
print(J)
print(f"理论值: [[2, 4], [2, 1]]")
# 输出:
# Jacobian 矩阵在 (1, 2) 处:
# [[2. 4.]
# [2. 1.]]
# 理论值: [[2, 4], [2, 1]]
11.31 Hessian 矩阵的数值计算
import numpy as np
def numerical_hessian(f, x, h=1e-4):
"""数值计算标量函数 f: R^n -> R 的 Hessian 矩阵"""
x = np.asarray(x, dtype=float)
n = len(x)
H = np.zeros((n, n))
f0 = f(x)
for i in range(n):
for j in range(n):
x_pp = x.copy(); x_pp[i] += h; x_pp[j] += h
x_pm = x.copy(); x_pm[i] += h; x_pm[j] -= h
x_mp = x.copy(); x_mp[i] -= h; x_mp[j] += h
x_mm = x.copy(); x_mm[i] -= h; x_mm[j] -= h
H[i, j] = (f(x_pp) - f(x_pm) - f(x_mp) + f(x_mm)) / (4 * h ** 2)
return H
# 示例:f(x, y) = x² + 2xy + 3y²
def f(xy):
x, y = xy
return x ** 2 + 2 * x * y + 3 * y ** 2
H = numerical_hessian(f, [1.0, 2.0])
print("Hessian 矩阵在 (1, 2) 处:")
print(H)
print(f"理论值: [[2, 2], [2, 6]]")
# 输出:
# Hessian 矩阵在 (1, 2) 处:
# [[2. 2.]
# [2. 6.]]
# 理论值: [[2, 2], [2, 6]]
11.32 数值积分与导数的关系
import numpy as np
from scipy import integrate
# 微积分基本定理:∫_a^b f'(x) dx = f(b) - f(a)
def f_prime(x):
return 2 * x # f(x) = x² + C 的导数
a, b = 0, 3
integral, _ = integrate.quad(f_prime, a, b)
print(f"∫_0^3 2x dx = {integral}")
print(f"f(3) - f(0) = {3**2 - 0**2}")
# 输出:
# ∫_0^3 2x dx = 9.0
# f(3) - f(0) = 9
第 12 章 对比分析
12.1 数值求导 vs 符号求导 vs 自动求导
| 特性 | 数值求导 | 符号求导 | 自动求导 |
|---|---|---|---|
| 原理 | 差商近似极限 | 表达式变换 + 求导规则 | 计算图 + 链式法则 |
| 精度 | 有限(受浮点限制) | 精确(数学上) | 精确(机器精度内) |
| 复杂度(前向) | ( 为表达式大小) | ||
| 复杂度(反向) | |||
| 内存 | |||
| 表达式膨胀 | 无 | 严重(易爆炸) | 无 |
| 分支控制 | 支持 | 难(需 if 表达式) | 支持 |
| 适用场景 | 快速估算、Hessian-vector 积 | 公式推导、教学演示 | 深度学习、优化器、灵敏度分析 |
| 高阶导数 | 需多次差分(误差累积) | 直接支持 | 需多次前向或反向 |
| 不可微点 | 难以判定 | 可静态分析 | 计算图分支可处理(subgradient) |
12.1.1 三类方法的关键差异
数值求导的本质是用差商 近似极限 ,受浮点精度约束存在不可消除的舍入误差;符号求导通过表达式变换(如 )在数学上精确,但易产生表达式膨胀(如 经 Leibniz 公式展开为 项);自动求导将函数分解为基本运算的有向无环图(DAG),对每个节点应用链式法则,既精确又避免表达式膨胀。
12.2 Newton 流数法 vs Leibniz 微分符号
| 维度 | Newton 流数法 | Leibniz 微分法 |
|---|---|---|
| 记号 | (点记号) | (d 记号) |
| 哲学基础 | 运动学(流动量与流数) | 几何学(无穷小三角形) |
| 优先级 | 1665-1666 手稿,1687《自然哲学的数学原理》 | 1675 手稿,1684《Acta Eruditorum》 |
| 高阶导数 | (点的累加不便) | (自然延展) |
| 偏导数 | 难以表达 | |
| 积分 | 不便表达 | (与求导对偶) |
| 现代地位 | 物理学中保留( 表示时间导数) | 数学分析主流 |
Newton 的流数术将变量视为”流动量”(fluents),其变化率为”流数”(fluxion),记作 ;这一记号在物理学中保留至今(如 表示广义速度、 表示广义加速度)。Leibniz 的微分记号 则将导数视为无穷小之商,其优势在于:高阶导数 可自然延展、偏导数 表达直观、积分 与求导形成符号对偶。
历史争议:1699 年至 1716 年间,英国皇家学会(受 Newton 影响)与欧洲大陆数学家(支持 Leibniz)就微积分发明优先权展开长达十余年的争论。后世研究表明,Newton 与 Leibniz 各自独立发明了微积分,且 Leibniz 的符号体系更适于运算与推广,因此现代分析学普遍采用 Leibniz 记号。
12.3 前向模式 AD vs 反向模式 AD
| 维度 | 前向模式(Forward Mode) | 反向模式(Reverse Mode) |
|---|---|---|
| 计算方向 | 输入 → 输出(与函数求值同向) | 输出 → 输入(与函数求值反向) |
| 基本原理 | Dual numbers 携带导数 | 计算图拓扑排序 + 链式法则反向传播 |
| 复杂度( 输入, 输出) | (每个输入一次前向) | (每个输出一次反向) |
| 典型实现 | JAX jax.jacfwd、Dual numbers | PyTorch autograd、TF GradientTape |
| 适用场景 | (输入少输出多,如 Jacobian 矩阵) | (输入多输出少,如损失函数对参数) |
| 内存 | (边算边丢) | (需保存中间结果供反向使用) |
| 神经网络反向传播 | 不适用(参数量 ) | 标准方法 |
前向模式基于 dual numbers:将每个变量 扩展为 ,其中 为 对某输入 的导数。基本运算规则为:
反向模式分两阶段:(1)前向阶段计算所有中间结果并保存依赖关系;(2)反向阶段从输出 出发,按拓扑逆序计算每个节点的伴随 。神经网络参数量远大于损失输出量(),故反向模式是深度学习的事实标准。
flowchart LR
subgraph Forward["前向模式(Forward Mode)"]
direction LR
X1["x₁"] --> V1["v₁"]
X2["x₂"] --> V1
V1 --> V2["v₂"]
X2 --> V2
V2 --> Y["y"]
V1 -.->|"携带 ∂v/∂x₁"| V1d["v₁, v̇₁"]
end
subgraph Reverse["反向模式(Reverse Mode)"]
direction RL
Y2["y"] -->|"ȳ = 1"| V2r["v̄₂"]
V2r -->|"v̄₁ = v̄₂ · ∂v₂/∂v₁"| V1r["v̄₁"]
V1r -->|"x̄₁ = v̄₁ · ∂v₁/∂x₁"| X1r["x̄₁"]
V1r -->|"x̄₂ = v̄₁ · ∂v₁/∂x₂ + v̄₂ · ∂v₂/∂x₂"| X2r["x̄₂"]
end
style Forward fill:#e3f2fd,color:#0d47a1
style Reverse fill:#fbe9e7,color:#bf360c
图解说明:前向模式(左)沿计算图正向传播,每个节点同时携带值 与对某输入 的导数 ;反向模式(右)先做前向求值保存中间结果,再从输出 出发沿拓扑逆序传播伴随 ,最终在叶节点累积得到所有输入的梯度。前向模式对每个输入需一次完整传播,反向模式对每个输出需一次完整传播,故神经网络(输入多、输出少)首选反向模式。
12.4 三类求导方法在神经网络中的分工
# 三类求导方法在神经网络训练中的分工示意
import torch
import sympy as sp
# 1. 符号求导:用于理论推导与梯度公式验证
x, w, b = sp.symbols('x w b', real=True)
y = sp.tanh(w * x + b)
grad_w = sp.diff(y, w)
print(f"符号求导 dy/dw = {grad_w}") # x*(1 - tanh(b + w*x)**2)
# 2. 数值求导:用于梯度检验(gradient check)
def numerical_grad(f, x, h=1e-5):
"""中心差分数值求导,用于梯度检验"""
return (f(x + h) - f(x - h)) / (2 * h)
w_val, b_val, x_val = 0.5, 0.1, 2.0
loss = lambda w: float(sp.tanh(w * x_val + b_val))
print(f"数值求导 dy/dw ≈ {numerical_grad(loss, w_val)}") # 与符号求导一致
# 3. 自动求导:用于实际训练(反向模式)
w_t = torch.tensor(w_val, requires_grad=True)
b_t = torch.tensor(b_val, requires_grad=True)
x_t = torch.tensor(x_val)
y_t = torch.tanh(w_t * x_t + b_t)
y_t.backward()
print(f"自动求导 dy/dw = {w_t.grad.item()}") # 与符号、数值一致,机器精度内
12.5 小结
数值求导、符号求导与自动求导构成求导方法的”三足鼎立”:数值求导简单快速但精度受限,符号求导精确但易爆炸,自动求导兼具精确与高效。神经网络训练中,符号求导用于理论推导与梯度公式验证,数值求导用于梯度检查(gradient check)以确保实现正确,自动求导(反向模式)用于实际的大规模参数优化。理解三者的原理、复杂度与适用场景,是工程实践中正确选择求导工具的基础。
第 13 章 常见陷阱
13.1 陷阱一:数值求导步长过小导致浮点吃掉有效数字
典型错误代码:
import math
def f(x):
return math.sin(x)
x = 0.1
h = 1e-20 # 错误:步长过小,2h 已退化为 0
df = (f(x + h) - f(x - h)) / (2 * h)
print(df) # 输出 0.0(错误,期望 ≈ 0.995)
错误原因:双精度浮点数的机器精度为 ,当 时, 与 在浮点表示下退化为相同值,差商分子 ,导致差商为 。
正确做法:中心差分的最优步长为 ,此时截断误差 与舍入误差 平衡。若需更高精度,使用 Richardson 外推(可达 )或符号求导。
13.2 陷阱二:PyTorch 中 requires_grad 与计算图未正确建立
典型错误代码:
import torch
x = torch.tensor([1.0, 2.0, 3.0]) # 默认 requires_grad=False
y = x.sum()
y.backward() # RuntimeError: element 0 of tensors does not require grad
错误原因:PyTorch 的 autograd 仅跟踪 requires_grad=True 张量参与的运算。默认张量不参与自动微分,调用 backward() 时若计算图中没有 requires_grad=True 的叶节点,则报错。
正确做法:
# 方案 1:创建时启用 requires_grad
x = torch.tensor([1.0, 2.0, 3.0], requires_grad=True)
y = x.sum()
y.backward()
print(x.grad) # tensor([1., 1., 1.])
# 方案 2:使用 torch.autograd.grad 函数式 API(不累积梯度,适合多次求导)
x = torch.tensor([1.0, 2.0, 3.0], requires_grad=True)
y = x.sum()
grad = torch.autograd.grad(y, x)[0]
print(grad) # tensor([1., 1., 1.])
# 注意:梯度默认累积,多次 backward 前需 zero_grad()
x = torch.tensor(2.0, requires_grad=True)
for _ in range(3):
y = x ** 2
y.backward()
print(x.grad) # tensor(24.)(= 3 次 2x = 6 累加,而非单次的 4)
13.3 陷阱三:链式法则伪证明(忽略 g(x) = g(x_0) 的退化情形)
伪证明:设 , 在 可导, 在 可导。令 ,,则
令 ,由 在 连续(可导必连续)知 ,故 ,,证毕。
错误原因:当 在 的某邻域内为常数(即 但 非局部常值)时,存在无穷多个 使 ,此时 为 未定式,不能直接断言极限为 。
正确证明:使用 Carathéodory 定义(见第 6 章 6.2 节)或 ε-δ 定义配合”扩展函数”技巧:定义
则 在 处连续(因 在 可导),且 。令 ,由 连续性与 连续性即得 。
13.4 陷阱四:混淆”可导”与”连续可导”( 类)
反例:函数 (),,在 处可导但导数不连续。
分析:
- 差商:(由夹逼定理 ),故 。
- 但 时 ,当 时 振荡无极限,故 在 不连续。
易错点:误认为”可导即 “。事实上,可导仅保证 存在,不保证 连续。 要求 存在且连续,是比”可导”更强的条件。Taylor 定理的 Peano 余项 仅要求 在 处 阶可导;而 Lagrange 余项 要求 在 的邻域内 阶可导且 连续。
13.5 陷阱五:L’Hôpital 法则的误用
误用场景一:未验证 或 未定式
实际极限为 。L’Hôpital 法则仅适用于 或 未定式,使用前必须验证分子分母同时趋于 或 。
误用场景二:极限不存在不等于原极限不存在
若用 L’Hôpital: 不存在。但原极限存在:。
正确认识:L’Hôpital 法则是”若 存在或为 ,则 “的单向蕴含,逆命题不成立。若 不存在(非 ),不能断言原极限不存在,应改用其他方法(如夹逼定理、Taylor 展开)。
13.6 陷阱六:复合函数求导时遗漏链式法则
典型错误:求 时,误写为 (遗漏了对内层 的求导)。
正确结果:。
根本原因:未识别复合结构 ,混淆了 与 。复合函数 的导数为 ,外层函数的导数需在内层函数处求值,再乘以内层函数的导数。
避免方法:养成”外层求导 → 内层求值 → 乘以内层求导”的三步流程,对每个复合层次逐层应用链式法则。如 。
第 14 章 工程实践
14.1 机器学习:梯度下降与反向传播
梯度下降是机器学习优化的基石,其核心是利用损失函数对参数的导数(梯度)指导参数更新方向。对于损失函数 ,参数更新规则为 ,其中 为学习率。
# 线性回归的梯度下降实现
import numpy as np
def linear_regression_gd(X, y, lr=0.01, epochs=1000):
"""
线性回归的批量梯度下降实现
输入:X (n, d) 特征矩阵, y (n,) 标签
输出:w (d,) 权重, b 偏置, loss_history 损失历史
核心流程:
1. 前向计算预测值 y_pred = X @ w + b
2. 计算损失 L = (1/2n) * ||y_pred - y||^2
3. 计算梯度 dL/dw = (1/n) * X^T @ (y_pred - y)
4. 计算梯度 dL/db = (1/n) * sum(y_pred - y)
5. 参数更新 w -= lr * dL/dw, b -= lr * dL/db
"""
n, d = X.shape
w = np.zeros(d)
b = 0.0
loss_history = []
for _ in range(epochs):
y_pred = X @ w + b
error = y_pred - y
loss = 0.5 * np.mean(error ** 2)
loss_history.append(loss)
# 解析梯度(基于 dL/dw = X^T @ error / n)
grad_w = X.T @ error / n
grad_b = np.mean(error)
w -= lr * grad_w
b -= lr * grad_b
return w, b, loss_history
# 示例:拟合 y = 2x + 3
np.random.seed(42)
X = np.linspace(0, 10, 100).reshape(-1, 1)
y = 2 * X.ravel() + 3 + np.random.randn(100) * 0.5
w, b, _ = linear_regression_gd(X, y, lr=0.01, epochs=1000)
print(f"w = {w[0]:.4f}, b = {b:.4f}") # 输出: w ≈ 1.9985, b ≈ 3.0241
变体:随机梯度下降(SGD,每次用 1 个样本)、小批量梯度下降(mini-batch,每次用 个样本)、Momentum()、Adam(自适应矩估计,结合一阶矩与二阶矩)。
14.2 物理学:速度、加速度与运动方程
导数在物理学中刻画瞬时变化率。位置 的导数为速度 ,速度的导数为加速度 。
# 匀变速直线运动的数值模拟
import numpy as np
import matplotlib.pyplot as plt
def simulate_motion(v0, a, t_max, dt=0.01):
"""
匀变速直线运动数值模拟
输入:v0 初速度, a 加速度, t_max 总时间, dt 时间步长
输出:t 时间序列, s 位置, v 速度
核心流程:基于位置与速度的微分关系 ds/dt = v, dv/dt = a, 使用前向 Euler 法
"""
t = np.arange(0, t_max, dt)
n = len(t)
s = np.zeros(n)
v = np.zeros(n)
v[0] = v0
for i in range(1, n):
# 数值求导的逆运算(数值积分):s_{i+1} = s_i + v_i * dt
s[i] = s[i-1] + v[i-1] * dt
v[i] = v[i-1] + a * dt
return t, s, v
# 自由落体(v0=0, a=-9.8)
t, s, v = simulate_motion(0, -9.8, 5)
# 解析解:s(t) = -4.9t^2, v(t) = -9.8t
s_analytic = -4.9 * t ** 2
v_analytic = -9.8 * t
print(f"位置最大误差: {np.max(np.abs(s - s_analytic)):.6f}") # ≈ 0
print(f"速度最大误差: {np.max(np.abs(v - v_analytic)):.6f}") # = 0
广义应用:简谐运动 中 、;电路中 、;热传导方程 。
14.3 经济学:边际分析与弹性
经济学中”边际”(marginal)即为导数。边际成本 、边际收益 、边际效用 。
# 经济学边际分析与弹性计算
import sympy as sp
q = sp.Symbol('q', positive=True)
# 成本函数 C(q) = 0.01q^3 - 0.6q^2 + 50q + 1000
C = 0.01 * q**3 - 0.6 * q**2 + 50 * q + 1000
# 收益函数 R(q) = 80q - 0.01q^2
R = 80 * q - 0.01 * q**2
MC = sp.diff(C, q) # 边际成本
MR = sp.diff(R, q) # 边际收益
profit = R - C
MProfit = sp.diff(profit, q) # 边际利润
# 利润最大化条件:MR = MC
q_opt = sp.solve(MR - MC, q)
print(f"最优产量 q* = {q_opt}") # [10.0, 150.0]
# 取正利润最大的解
q_star = 150
profit_max = profit.subs(q, q_star)
print(f"最大利润 = {profit_max}") # 1500
# 价格弹性:E = (dQ/dP) * (P/Q)
P = sp.Symbol('P', positive=True)
Q = 100 - 2 * P # 需求函数
elasticity = sp.diff(Q, P) * (P / Q)
print(f"需求价格弹性 E = {elasticity}") # -2*P/(100 - 2*P)
经济意义:当边际收益等于边际成本()时利润最大化;当需求价格弹性 时降价增收, 时提价增收, 时收益最大化。
14.4 数值优化:Newton 法与拟 Newton 法
Newton 法利用一阶导数与二阶导数(Hessian)加速求根与优化。求 的迭代格式为 ;求 极小值的迭代格式为 。
# Newton 法求根:求解 x^2 - 2 = 0(即 sqrt(2))
def newton_sqrt(initial=1.0, tol=1e-12, max_iter=50):
"""
Newton 法求 sqrt(2)
输入:initial 初始值, tol 容差, max_iter 最大迭代次数
输出:x 近似解, history 迭代历史
核心流程:x_{n+1} = x_n - f(x_n)/f'(x_n), f(x) = x^2 - 2
"""
x = initial
history = [x]
for _ in range(max_iter):
fx = x ** 2 - 2
dfx = 2 * x
x_new = x - fx / dfx
history.append(x_new)
if abs(x_new - x) < tol:
return x_new, history
x = x_new
return x, history
x_star, hist = newton_sqrt()
print(f"sqrt(2) ≈ {x_star:.15f}") # 1.414213562373095
print(f"迭代次数: {len(hist) - 1}") # 6 次(二次收敛)
# 多维 Newton 法优化:Rosenbrock 函数
import numpy as np
def newton_optimize(grad_f, hessian_f, x0, tol=1e-10, max_iter=100):
"""
多维 Newton 法求极小值
输入:grad_f 梯度函数, hessian_f Hessian 函数, x0 初始点
输出:x* 极小点, history 迭代历史
核心流程:x_{k+1} = x_k - H^{-1}(x_k) * grad(x_k)
"""
x = np.array(x0, dtype=float)
history = [x.copy()]
for _ in range(max_iter):
g = grad_f(x)
H = hessian_f(x)
delta = np.linalg.solve(H, g)
x_new = x - delta
history.append(x_new.copy())
if np.linalg.norm(delta) < tol:
return x_new, history
x = x_new
return x, history
# Rosenbrock 函数:f(x, y) = (1-x)^2 + 100(y-x^2)^2, 极小点 (1, 1)
def rosen_grad(x):
return np.array([-2*(1-x[0]) - 400*x[0]*(x[1]-x[0]**2),
200*(x[1]-x[0]**2)])
def rosen_hess(x):
return np.array([[2 + 1200*x[0]**2 - 400*x[1], -400*x[0]],
[-400*x[0], 200]])
x_opt, _ = newton_optimize(rosen_grad, rosen_hess, [-1.2, 1.0])
print(f"Rosenbrock 极小点: {x_opt}") # ≈ [1.0, 1.0]
收敛性:Newton 法在根附近具有二次收敛性 ,远快于梯度下降的线性收敛。但 Newton 法要求 Hessian 可逆且每步需 求逆,故高维问题常用拟 Newton 法(BFGS、L-BFGS)以 近似 Hessian 逆。
14.5 信号处理:边缘检测与滤波
导数在信号处理中用于边缘检测、特征提取与滤波器设计。图像边缘对应于像素强度的不连续点,可通过梯度幅值 检测。
# 图像边缘检测:Sobel 算子(基于一阶差分)
import numpy as np
def sobel_edge(image):
"""
Sobel 算子边缘检测
输入:image (H, W) 灰度图像
输出:gradient_magnitude (H, W) 梯度幅值
核心流程:
1. 用 Sobel_x 算子做水平方向差分(近似 dI/dx)
2. 用 Sobel_y 算子做垂直方向差分(近似 dI/dy)
3. 计算梯度幅值 = sqrt(Gx^2 + Gy^2)
"""
# Sobel 算子本质上是中心差分 + 平滑
Kx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
Ky = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])
H, W = image.shape
Gx = np.zeros((H-2, W-2))
Gy = np.zeros((H-2, W-2))
for i in range(H-2):
for j in range(W-2):
patch = image[i:i+3, j:j+3]
Gx[i, j] = np.sum(patch * Kx)
Gy[i, j] = np.sum(patch * Ky)
return np.sqrt(Gx**2 + Gy**2)
# 测试:合成图像边缘检测
img = np.zeros((10, 10))
img[:, 5:] = 1.0 # 垂直边缘在第 5 列
edge = sobel_edge(img)
print(f"边缘强度峰值位置: 列 {np.argmax(edge[5, :])}") # 列 3(因 padding 偏移)
进阶:Laplacian 算子 用于零交叉边缘检测;Canny 边缘检测结合高斯平滑、Sobel 梯度、非极大值抑制与双阈值滞后,是工业级标准;二阶导数过零点(zero-crossing)对应于一阶导数的局部极大值,刻画图像边缘的精确位置。
第 15 章 案例研究
15.1 案例一:PyTorch autograd 机制深度剖析
PyTorch 的 autograd 模块基于动态计算图(define-by-run)实现反向模式自动微分。每次前向运算时动态构建计算图,调用 backward() 时按拓扑逆序遍历计算图并应用链式法则。
# PyTorch autograd 内部机制剖析
import torch
# 1. 叶节点:requires_grad=True 的张量为计算图的叶节点
x = torch.tensor(2.0, requires_grad=True)
w = torch.tensor(3.0, requires_grad=True)
b = torch.tensor(1.0, requires_grad=True)
# 2. 前向运算:每一步运算动态生成新节点,记录父节点与梯度函数(grad_fn)
y = w * x # MulBackward: y = w * x
z = y + b # AddBackward: z = y + b
loss = z ** 2 # PowBackward: loss = z^2
# 3. 计算图可视化
print(f"loss.grad_fn = {loss.grad_fn}") # <PowBackward0>
print(f"z.grad_fn = {z.grad_fn}") # <AddBackward0>
print(f"y.grad_fn = {y.grad_fn}") # <MulBackward0>
print(f"loss.grad_fn.next_functions[0][0] = "
f"{loss.grad_fn.next_functions[0][0]}") # <AddBackward0>
# 4. 反向传播:从 loss 出发,按拓扑逆序遍历计算图
loss.backward()
# 等价于:dloss/dz = 2z = 2*(wx+b) = 2*(3*2+1) = 14
# dloss/dy = dloss/dz * dz/dy = 14 * 1 = 14
# dloss/dw = dloss/dy * dy/dw = 14 * x = 14 * 2 = 28
# dloss/dx = dloss/dy * dy/dx = 14 * w = 14 * 3 = 42
# dloss/db = dloss/dz * dz/db = 14 * 1 = 14
print(f"dl/dx = {x.grad}") # tensor(42.)
print(f"dl/dw = {w.grad}") # tensor(28.)
print(f"dl/db = {b.grad}") # tensor(14.)
# 5. 计算图释放:默认 backward 后释放计算图,再次 backward 报错
try:
loss.backward()
except RuntimeError as e:
print(f"再次 backward 报错: {e}") # Trying to backward through the graph a second time
# 6. retain_graph=True 保留计算图供多次 backward
x = torch.tensor(2.0, requires_grad=True)
w = torch.tensor(3.0, requires_grad=True)
y = w * x
z = y ** 2
z.backward(retain_graph=True)
print(f"第一次 dz/dx = {x.grad}") # tensor(12.) (2*wx*w = 2*3*2*3/... 实际为 2*w^2*x = 2*9*2=36? 校验)
z.backward()
print(f"第二次累积 dz/dx = {x.grad}") # tensor(24.)(梯度累积)
# 7. 高阶导数:create_graph=True 在反向传播时构建新计算图
x = torch.tensor(1.5, requires_grad=True)
y = x ** 3 # y = x^3, y' = 3x^2, y'' = 6x
first_grad = torch.autograd.grad(y, x, create_graph=True)[0]
print(f"一阶导 dy/dx = {first_grad}") # tensor(6.75) (= 3*1.5^2)
second_grad = torch.autograd.grad(first_grad, x)[0]
print(f"二阶导 d²y/dx² = {second_grad}") # tensor(9.) (= 6*1.5)
关键设计:PyTorch 为每种基本运算实现对应的 grad_fn,记录前向运算的反向规则。如 MulBackward 的反向规则为 、。反向传播时,autograd 引擎按拓扑逆序调用每个 grad_fn,将上游梯度与本地梯度相乘并累积到叶节点的 .grad 属性。
15.2 案例二:TensorFlow GradientTape 机制
TensorFlow 2.x 的 tf.GradientTape 采用与 PyTorch 类似的动态图机制,但通过 with 上下文管理器显式开启梯度记录范围。
# TensorFlow GradientTape 机制
import tensorflow as tf
# 1. 基本用法:在 GradientTape 上下文中的运算被记录
x = tf.Variable(2.0)
w = tf.Variable(3.0)
b = tf.Variable(1.0)
with tf.GradientTape() as tape:
y = w * x
z = y + b
loss = z ** 2
# 同时计算 loss 对多个变量的梯度
grads = tape.gradient(loss, [x, w, b])
print(f"dl/dx = {grads[0]}") # 42.0
print(f"dl/dw = {grads[1]}") # 28.0
print(f"dl/db = {grads[2]}") # 14.0
# 2. persistent=True:允许多次调用 tape.gradient(默认一次性)
x = tf.Variable(2.0)
with tf.GradientTape(persistent=True) as tape:
y = x ** 2
z = x ** 3
print(f"dy/dx = {tape.gradient(y, x)}") # 4.0 (= 2x)
print(f"dz/dx = {tape.gradient(z, x)}") # 12.0 (= 3x^2)
del tape # 显式释放 persistent tape
# 3. 高阶导数:嵌套 GradientTape
x = tf.Variable(1.5)
with tf.GradientTape() as tape2:
with tf.GradientTape() as tape1:
y = x ** 3
first_grad = tape1.gradient(y, x)
second_grad = tape2.gradient(first_grad, x)
print(f"一阶导 = {first_grad}, 二阶导 = {second_grad}") # 6.75, 9.0
# 4. JVP(前向模式)与 VJP(反向模式)
# TensorFlow 默认使用 VJP(反向模式)
# 通过 tf.autodiff.ForwardAccumulator 可启用 JVP(前向模式)
x = tf.Variable(2.0)
with tf.autodiff.ForwardAccumulator(
primals=x, tangents=tf.constant(1.0)) as acc:
y = x ** 3
print(f"JVP 前向模式 dy/dx = {acc.jvp(y)}") # 12.0 (= 3x^2)
对比:PyTorch 默认动态图、backward() 后释放计算图、梯度累积;TensorFlow 默认动态图(TF 2.x)、GradientTape 上下文结束后才能求梯度、persistent=True 支持多次求梯度。两者均以反向模式为主,前向模式用于 Jacobian 计算。
15.3 案例三:从零实现神经网络反向传播
反向传播算法(backpropagation)是链式法则在多层神经网络中的高效实现。对于一个 层的全连接神经网络,前向传播与反向传播的公式如下:
前向传播:,,其中 为激活函数。
反向传播:,,,。
# 从零实现两层神经网络的反向传播(不使用 PyTorch/TF)
import numpy as np
class TwoLayerNet:
"""
两层全连接神经网络(输入层 → 隐藏层 → 输出层)
输入:x (n, d_in)
输出:y_pred (n, d_out)
核心流程:
1. 前向传播:z1 = W1 @ x + b1, a1 = relu(z1), z2 = W2 @ a1 + b2, y = softmax(z2)
2. 计算损失:L = -sum(y_true * log(y_pred))(交叉熵)
3. 反向传播:基于链式法则逐层计算 dL/dW2, dL/db2, dL/dW1, dL/db1
4. 参数更新:梯度下降 W -= lr * dL/dW
"""
def __init__(self, d_in, d_hidden, d_out):
# He 初始化(针对 ReLU)
self.W1 = np.random.randn(d_in, d_hidden) * np.sqrt(2.0 / d_in)
self.b1 = np.zeros(d_hidden)
self.W2 = np.random.randn(d_hidden, d_out) * np.sqrt(2.0 / d_hidden)
self.b2 = np.zeros(d_out)
def relu(self, x):
return np.maximum(0, x)
def relu_grad(self, x):
return (x > 0).astype(float)
def softmax(self, x):
exp_x = np.exp(x - np.max(x, axis=1, keepdims=True))
return exp_x / np.sum(exp_x, axis=1, keepdims=True)
def forward(self, x):
"""
前向传播:保存中间结果供反向传播使用
输入:x (n, d_in)
输出:y_pred (n, d_out), cache 中间结果
"""
z1 = x @ self.W1 + self.b1 # (n, d_hidden)
a1 = self.relu(z1) # (n, d_hidden)
z2 = a1 @ self.W2 + self.b2 # (n, d_out)
y_pred = self.softmax(z2) # (n, d_out)
cache = {'x': x, 'z1': z1, 'a1': a1, 'z2': z2, 'y_pred': y_pred}
return y_pred, cache
def backward(self, cache, y_true):
"""
反向传播:基于链式法则计算各参数梯度
输入:cache 前向传播中间结果, y_true (n, d_out) one-hot 标签
输出:grads = {W1, b1, W2, b2}
核心流程:
1. 输出层 delta:dL/dz2 = y_pred - y_true(softmax + 交叉熵的简化形式)
2. 输出层参数梯度:dL/dW2 = a1^T @ delta2, dL/db2 = sum(delta2, axis=0)
3. 隐藏层 delta:dL/da1 = delta2 @ W2^T, dL/dz1 = dL/da1 * relu'(z1)
4. 隐藏层参数梯度:dL/dW1 = x^T @ delta1, dL/db1 = sum(delta1, axis=0)
"""
n = y_true.shape[0]
# 输出层 delta(softmax + 交叉熵的雅可比简化为 y_pred - y_true)
delta2 = cache['y_pred'] - y_true # (n, d_out)
grad_W2 = cache['a1'].T @ delta2 / n # (d_hidden, d_out)
grad_b2 = np.mean(delta2, axis=0) # (d_out,)
# 隐藏层 delta(链式法则:dL/da1 = delta2 @ W2^T, 再乘以 ReLU 的导数)
delta1 = (delta2 @ self.W2.T) * self.relu_grad(cache['z1']) # (n, d_hidden)
grad_W1 = cache['x'].T @ delta1 / n # (d_in, d_hidden)
grad_b1 = np.mean(delta1, axis=0) # (d_hidden,)
return {'W1': grad_W1, 'b1': grad_b1, 'W2': grad_W2, 'b2': grad_b2}
def train_step(self, x, y_true, lr=0.1):
y_pred, cache = self.forward(x)
grads = self.backward(cache, y_true)
self.W1 -= lr * grads['W1']
self.b1 -= lr * grads['b1']
self.W2 -= lr * grads['W2']
self.b2 -= lr * grads['b2']
loss = -np.mean(np.sum(y_true * np.log(cache['y_pred'] + 1e-10), axis=1))
return loss
# 测试:训练一个 XOR 分类器
np.random.seed(42)
X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]], dtype=float)
y = np.array([[1, 0], [0, 1], [0, 1], [1, 0]], dtype=float) # XOR: 同或
net = TwoLayerNet(2, 8, 2)
for epoch in range(1000):
loss = net.train_step(X, y, lr=0.5)
if epoch % 200 == 0:
y_pred, _ = net.forward(X)
acc = np.mean(np.argmax(y_pred, axis=1) == np.argmax(y, axis=1))
print(f"Epoch {epoch}: loss = {loss:.4f}, acc = {acc:.2%}")
# 输出(典型):
# Epoch 0: loss = 0.6931, acc = 50.00%
# Epoch 200: loss = 0.0051, acc = 100.00%
# Epoch 400: loss = 0.0021, acc = 100.00%
# ...
关键洞察:反向传播的本质是链式法则在计算图上的高效组织——通过”反向”传播输出层误差,避免了”前向”模式中对每个参数单独计算梯度的 重复开销,使总复杂度降为 (与参数量同阶)。这是深度学习能够训练上亿参数模型的算法基础。
第 16 章 习题与解答
本章节汇总 frontmatter 中 10 道习题的参考答案要点,便于学习者快速查阅。完整题目与评分关键点见 frontmatter exercises 字段。
16.1 填空题
习题 ex-calc-diff-fb-01(ε-δ 定义填空) 答案:。差商 在 时趋于 ,ε-δ 定义刻画了这一极限过程。
习题 ex-calc-diff-fb-02(Carathéodory 定义填空) 答案:。Carathéodory 等价定义将差商结构提取为连续函数 ,使 ,且 。
16.2 选择题
习题 ex-calc-diff-ch-01( 在 0 处可导性) 答案:B。,,左右导数存在但不相等,故不可导。 在 0 处连续,排除 C;可导必蕴含导数存在,排除 A、D。
习题 ex-calc-diff-ch-02( 在 0 处可导性与导数连续性) 答案:B。由夹逼定理 知 ;但 时 在 时因 振荡而无极限,故 在 0 处不连续。这是”可导但导数不连续”的经典反例。
习题 ex-calc-diff-ch-03(Lagrange 中值定理理解) 答案:B。Lagrange 中值定理的几何意义是切线斜率等于端点弦斜率。 是横坐标而非原函数值,且不一定唯一(如 在 上有三个 满足条件),也不必为中点。
16.3 代码修正题
习题 ex-calc-diff-cf-01(数值求导步长修正) 修正要点:将 调整为 (中心差分最优步长 ),或使用 Richardson 外推达机器精度。
习题 ex-calc-diff-cf-02(PyTorch requires_grad 修正)
修正要点:在创建张量时设置 requires_grad=True,或对已有张量调用 requires_grad_(True),使其成为计算图叶节点。
16.4 开放性问题
习题 ex-calc-diff-oe-01(ε-δ 证明 在 2 处可导) 证明要点:取 (或更保守的 ),差商 。 仅依赖于 ,不依赖于具体 。
习题 ex-calc-diff-oe-02(Carathéodory 定义证明链式法则) 证明要点:设 、,则 。由 与 均连续知 连续,故 。
习题 ex-calc-diff-oe-03(Lagrange 中值定理证明 Rolle 定理的逆命题讨论) 讨论要点:Rolle 定理是 Lagrange 中值定理在 时的特例;逆命题”若 则 是极值点”不成立(反例: 在 处 但非极值)。正确的相关结论是 Fermat 引理:若 是 的极值点且 在 处可导,则 。
第 17 章 参考文献
本章节列出 frontmatter references 字段中已声明的 12 条 ACM Reference Format 参考文献,按引用顺序排列。完整字段(作者、年份、标题、出版商、DOI 等)见 frontmatter。
- Spivak, M. (2008). Calculus (4th ed.). Publish or Perish, Inc. https://doi.org/10.1007/978-0-387-09469-9
- Apostol, T. M. (1967). Calculus, Volume 1: One-Variable Calculus with an Introduction to Linear Algebra (2nd ed.). John Wiley & Sons.
- Rudin, W. (1976). Principles of Mathematical Analysis (3rd ed.). McGraw-Hill Education.
- Tao, T. (2016). Analysis I (3rd ed.). Springer. https://doi.org/10.1007/978-981-10-1789-6
- Courant, R., & John, F. (1999). Introduction to Calculus and Analysis I. Springer.
- Hardy, G. H. (1952). A Course of Pure Mathematics (10th ed.). Cambridge University Press.
- Pugh, C. C. (2002). Real Mathematical Analysis. Springer. https://doi.org/10.1007/978-0-387-21668-2
- Carathéodory, C. (1950). Vorlesungen über reelle Funktionen (3rd ed.). Chelsea Publishing Company.
- Baydin, A. G., Pearlmutter, B. A., Radul, A. A., & Siskind, J. M. (2018). Automatic differentiation in machine learning: a survey. Journal of Machine Learning Research, 18(153), 1-43.
- Rumelhart, D. E., Hinton, G. E., & Williams, R. J. (1986). Learning representations by back-propagating errors. Nature, 323(6088), 533-536. https://doi.org/10.1038/323533a0
- Leibniz, G. W. (1684). Nova methodus pro maximis et minimis, itemque tangentibus, quae nec fractas, nec irrationales quantitates moratur, et singulare pro illis calculi genus. Acta Eruditorum.
- Newton, I. (1687). Philosophiae Naturalis Principia Mathematica. Royal Society of London.
引用规范说明:本节遵循 ACM Reference Format,期刊文献包含卷号、期号、页码与 DOI,书籍包含版次与出版社。在正文中以”Spivak (2008)“或”(Rumelhart et al., 1986)“形式引用。
第 18 章 延伸阅读
18.1 关联模块
- calculus/函数与极限:本篇的前置基础,介绍极限的 ε-δ 定义、Heine 定理、Cauchy 收敛准则与连续性理论。建议在学习本篇前完整阅读该文档。
- calculus/微分中值定理:本篇第 9 章的深化,详细证明 Rolle、Lagrange、Cauchy 中值定理的应用与广义中值定理。
- calculus/不定积分:导数的逆运算,介绍原函数、换元法、分部积分与有理函数积分。
- calculus/定积分与应用:基于 Riemann 积分理论介绍定积分定义、可积性条件、微积分基本定理与几何/物理应用。
- calculus/重积分:多元函数微分学的延伸,介绍偏导数、全微分、方向导数与多元 Taylor 公式。
- math/实分析:本篇的理论深化,介绍 Lebesgue 测度与积分、 空间、测度论视角下的微分与积分。
- machine-learning/反向传播:本篇第 14-15 章的工程深化,介绍现代深度学习框架的反向传播实现细节。
- numerical-analysis/数值微分与积分:本篇第 11 章的算法深化,介绍 Richardson 外推、Romberg 积分、Gauss 求积等高精度方法。
18.2 进阶资料
数学分析深化:
- Rudin, W. (1987). Real and Complex Analysis (3rd ed.). McGraw-Hill. —— 测度论视角下的微分与积分,介绍 Lebesgue 微分定理、Radon-Nikodym 导数。
- Folland, G. B. (1999). Real Analysis: Modern Techniques and Their Applications (2nd ed.). Wiley. —— 现代分析学标准教材,涵盖测度论、傅里叶分析与分布理论。
- Royden, H. L., & Fitzpatrick, P. M. (2010). Real Analysis (4th ed.). Pearson. —— 实分析入门经典,对 Lebesgue 测度的讲解尤为清晰。
自动微分与深度学习:
- Goodfellow, I., Bengio, Y., & Courville, A. (2016). Deep Learning. MIT Press. —— 深度学习权威教材,第 6 章系统介绍反向传播与计算图。
- Griewank, A., & Walther, A. (2008). Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation (2nd ed.). SIAM. —— 自动微分专著,前向与反向模式的理论与实现。
- Paszke, A. et al. (2017). Automatic differentiation in PyTorch. NIPS Autodiff Workshop. —— PyTorch autograd 的设计原理。
历史与哲学:
- Boyer, C. B., & Merzbach, U. C. (2011). A History of Mathematics (3rd ed.). Wiley. —— 微积分发展史的权威叙述。
- Edwards, C. H. (1979). The Historical Development of the Calculus. Springer. —— 从古希腊到 19 世纪的微积分演进。
- Bressoud, D. (2011). Historical Reflections on Teaching the Fundamental Theorem of Calculus. American Mathematical Monthly, 118(2), 99-115.
18.3 相关模块
- linear-algebra/矩阵微分:标量对矩阵的导数、Jacobian 与 Hessian 矩阵、矩阵求导的链式法则。
- probability/概率密度与分布函数:分布函数的导数为概率密度、矩母函数与特征函数的导数性质。
- optimization/凸优化:凸函数的可微性、次梯度、KKT 条件中的梯度互补条件。
- differential-equations/常微分方程:导数作为方程变量,介绍解的存在唯一性、稳定性与相图分析。
18.4 推荐学习路径
- 入门阶段(2-4 周):先阅读
calculus/函数与极限,掌握极限与连续性;再学习本篇第 1-7 章,掌握导数定义与求导法则。 - 进阶阶段(3-5 周):学习本篇第 8-10 章的 Taylor 定理与中值定理,结合
math/实分析深化对 ε-δ 语言的理解。 - 工程阶段(2-3 周):学习本篇第 11-15 章的代码实现与案例研究,运行所有 Python 示例并完成 10 道习题。
- 应用阶段(按需):根据兴趣选择
machine-learning/反向传播、numerical-analysis/数值微分与积分、differential-equations/常微分方程等延伸模块。
18.5 社区与讨论
- Mathematics Stack Exchange:数学问题的权威问答社区,可搜索 ε-δ 证明、Taylor 定理等历史讨论。
- PyTorch Discussion Forum:PyTorch autograd 的实现细节与最佳实践讨论。
- arXiv math.CA / cs.LG:经典分析与机器学习最新论文预印本,跟踪自动微分理论前沿。
18.6 致谢
本篇文档的写作受益于 Spivak《Calculus》4th edition、Apostol《Calculus Vol 1》、Rudin《Principles of Mathematical Analysis》3rd edition 与 Tao《Analysis I》3rd edition 的深刻启发。PyTorch autograd 与 TensorFlow GradientTape 的设计理念源自 Baydin et al. (2018) 的自动微分综述。反向传播算法的历史可追溯至 Rumelhart, Hinton & Williams (1986) 在 Nature 上发表的奠基性论文。谨向这些先驱者致以最深的敬意。
更新日志
| 版本 | 日期 | 修改内容 | 修改人 |
|---|---|---|---|
| 1.0.0 | 2026-06-14 | 初始版本,包含基本导数定义、求导法则与少量示例 | fanquanpp |
| 2.0.0 | 2026-07-18 | 论文级重写,覆盖 12 项质量基准:扩展至 18 章 2900+ 行;新增 ε-δ 与 Carathéodory 双定义、链式法则严格证明、Rolle/Lagrange/Cauchy 中值定理、Taylor 定理 4 种余项、Weierstrass 反例;新增 40+ Python 代码示例(数值/符号/自动求导 + PyTorch/TF 反向传播 + Newton 法 + Sobel 边缘检测);新增 6 类常见陷阱、5 类工程实践、3 类案例研究;新增 10 道 Spivak 风格习题(含 3 道证明题)、12 条 ACM 格式参考文献、6 条词源标注;新增 4 个 Mermaid 图与 50+ KaTeX 公式;对齐 Phase 1.5 frontmatter schema | FANDEX Content Engineering |