1. 项目概述当方程无法显式求解时我们如何求导在工程优化、物理模拟和机器学习等领域我们常常会遇到这样的问题一个变量依赖于另一个变量但它们之间的关系并非一个简单的y f(x)公式而是被一个复杂的方程F(x, y) 0所隐含地定义。比如一个机械臂的末端位置(x, y)可能满足一个包含关节角度的复杂三角方程或者在经济学模型中均衡价格P和数量Q由市场供需方程组共同决定无法直接解出P f(Q)。这时候如果你想分析y如何随x变化即求dy/dx或者需要计算梯度进行数值优化直接求导的路就走不通了因为你连y的显式表达式都没有。隐函数定理Implicit Function Theorem就是解决这类问题的“尚方宝剑”。它不关心y具体等于什么而是直接告诉你在满足一定条件下y作为x的函数是存在的、可微的并且其导数可以通过原方程F的偏导数直接计算出来。这个项目就是要把这把理论上的“宝剑”打磨成一把可以在 C 代码中挥舞的实用“利刃”。我们将从最核心的一元隐函数求导公式dy/dx -F_x / F_y出发一直深入到多元方程组情形下的雅可比矩阵求逆运算并最终用 C 实现一个通用的、支持自动微分的隐函数偏导数计算模块。无论你是正在学习《高等数学》或《数学分析》的学生还是需要在仿真、优化算法中处理隐式约束的工程师这篇文章都将带你从理论到实践彻底掌握隐函数求导及其 C 实现的核心要领。2. 隐函数定理的核心思想与公式推导理解隐函数定理关键在于抓住其“局部存在性”和“导数可计算性”这两个核心。它不像代数求解那样追求全局的、显式的解而是保证在某个特定的点附近解是存在的、唯一的并且其微分行为完全由原函数在该点的局部性质决定。2.1 从直观例子理解定理条件让我们先看一个经典的例子单位圆方程x^2 y^2 - 1 0。我们想知道在点(0, 1)即圆的最顶端附近y能否表示为x的函数以及其导数。连续性条件C¹方程F(x, y) x^2 y^2 - 1本身及其偏导数F_x 2x,F_y 2y在(0, 1)附近都是连续且光滑的。这保证了函数行为是“良好”的没有奇点或突变。零点条件显然F(0, 1) 0点(0, 1)在圆上。非退化条件计算F_y(0, 1) 2 * 1 2 ≠ 0。这是最关键的条件。F_y ≠ 0意味着在(0, 1)点方程F(x, y)0在y方向上是“非奇异”的y的变化对F的值有显著影响。从几何上看圆在该点的切线不是竖直的。当这三个条件满足时隐函数定理断言在(0, 1)附近的一个小邻域内存在唯一的函数y φ(x)使得F(x, φ(x)) ≡ 0。并且这个隐函数φ(x)也是连续可微的。注意这个“局部”特性非常重要。对于整个圆来说y不能写成x的单一函数因为一个x对应两个y。但在(0, 1)这个点的附近比如右上半圆的一小段它是可以的。同样在(1, 0)点F_y 0定理条件不满足此时圆在该点的切线是竖直的y无法表示为x的函数但x可以表示为y的函数。2.2 核心公式的推导链式法则的巧妙应用定理不仅保证了隐函数存在还给出了计算其导数的公式。推导过程简洁而优美是链式法则的经典应用。假设由F(x, y) 0确定了隐函数y φ(x)且F和φ都可微。我们将方程F(x, φ(x)) 0的两边对x求导。这里F是一个二元函数x是直接变量y通过φ(x)也依赖于x。应用多元函数的链式法则d/dx [ F(x, φ(x)) ] (∂F/∂x) * (dx/dx) (∂F/∂y) * (dφ/dx) 0由于dx/dx 1上式简化为F_x(x, φ(x)) F_y(x, φ(x)) * φ(x) 0只要F_y ≠ 0我们就可以立即解出φ(x) - F_x(x, φ(x)) / F_y(x, φ(x))这就是一元隐函数求导的万能公式。你不需要知道φ(x)的具体形式只需要知道原函数F和该点的坐标(x, y)就能算出隐函数的导数。实操心得这个推导过程揭示了隐函数求导的本质——对方程两端同时关于自变量求全微分。对于更复杂的多元隐函数或方程组这个思想依然适用只是从标量方程升级到了矩阵方程。2.3 推广到多元隐函数与方程组实际问题往往是多维的。例如由方程组F(u, v, x, y) 0 G(u, v, x, y) 0确定了两个隐函数u u(x, y),v v(x, y)。我们想求∂u/∂x,∂v/∂x等。此时我们将所有方程对x求偏导视y为常数再次利用链式法则∂F/∂x (∂F/∂u)*(∂u/∂x) (∂F/∂v)*(∂v/∂x) (∂F/∂x)*(1) (∂F/∂y)*(0) 0 ∂G/∂x (∂G/∂u)*(∂u/∂x) (∂G/∂v)*(∂v/∂x) (∂G/∂x)*(1) (∂G/∂y)*(0) 0注意这里的∂F/∂x是指F对第一个显式变量x的偏导而∂u/∂x才是我们要求的隐函数偏导。整理后得到一个关于∂u/∂x和∂v/∂x的线性方程组[ ∂F/∂u ∂F/∂v ] [ ∂u/∂x ] [ -∂F/∂x ] [ ∂G/∂u ∂G/∂v ] [ ∂v/∂x ] [ -∂G/∂x ]左边的矩阵正是雅可比矩阵J_{(u,v)}F函数组(F, G)对变量组(u, v)的偏导数矩阵。只要这个雅可比矩阵在求解点可逆即行列式非零这是高维下的“非退化条件”我们就可以解出[ ∂u/∂x ] -1 [ -∂F/∂x ] [ ∂v/∂x ] [J_{(u,v)}F] * [ -∂G/∂x ]这就是隐函数定理的矩阵形式。它告诉我们隐函数组的导数等于原函数组对隐变量雅可比矩阵的逆乘以原函数组对显式变量的偏导数的负值。3. C实现的核心架构与设计选择将数学定理转化为代码我们需要做出几个关键的设计决策。我们的目标是构建一个灵活、高效且易于使用的ImplicitDifferentiator类。3.1 核心计算流程设计整个计算流程可以分解为以下几步这构成了我们类方法的骨架定义方程与变量用户提供方程F(x, y, ...) 0并指定哪些是自变量x哪些是因变量y即隐函数。计算偏导数计算F对所有变量的偏导数。这是最核心也最耗时的步骤。构造雅可比矩阵根据因变量和自变量的划分从偏导数集合中分别提取出J_yF对因变量的雅可比和J_xF对自变量的偏导向量/矩阵。求解线性系统计算- (J_yF)^{-1} * J_xF。对于单方程单隐函数这就是一个除法对于方程组这就是一个矩阵求逆或求解线性方程组的问题。输出结果返回隐函数关于自变量的偏导数标量、向量或矩阵。3.2 关键设计决策如何表示和计算偏导数这里有三个主流方案各有优劣方案一符号微分Symbolic Differentiation思路像 Mathematica 或 SymPy 一样解析地推导出偏导数的表达式。C实现需要构建表达式树实现符号化简和求导规则。可以使用像SymEngine这样的库。优点求导精确一次推导后可重复使用表达式对于复杂函数可能更快。缺点实现极其复杂表达式膨胀可能导致性能下降且难以处理控制流复杂的函数。结论对于专注于隐函数求导的专用工具来说过于重型不推荐。方案二数值微分Numerical Differentiation思路用有限差分近似偏导数例如∂F/∂x ≈ (F(xh, y) - F(x, y)) / h。C实现非常简单只需要能调用函数F即可。优点实现简单对函数F的形式没有任何要求即使是黑箱函数。缺点存在截断误差和舍入误差选择步长h需要技巧太小会放大舍入误差太大则截断误差大。对于条件数大的雅可比矩阵误差会被放大可能导致求解不稳定。结论适合快速原型验证或当F是黑箱函数时使用。但作为通用库的核心精度和稳定性是短板。方案三自动微分Automatic Differentiation, AD思路介于符号和数值之间。通过操作符重载在计算函数值的同时精确地计算其导数。它计算的是精确的导数值没有截断误差。C实现需要定义一个特殊的“双数”类型例如Dualdouble重载所有算术和初等函数操作。优点计算精度高机器精度效率通常优于符号微分且易于使用。缺点需要将函数F用重载了 AD 的类型重写。有轻微的计算和内存开销。结论这是本项目的最佳选择。它在精度、效率和易用性之间取得了最佳平衡。现代科学计算库如 Ceres Solver, Stan的核心都依赖于 AD。我们选择使用前向模式自动微分Forward Mode AD来实现。对于计算一个标量函数对多个变量的梯度前向模式非常直观高效。我们将实现一个简单的Dual类作为基础。3.3 类与接口设计基于以上分析我们设计以下核心类// 前向模式自动微分的双数类 templatetypename T class Dual { public: T value; // 函数值 T deriv; // 导数值 Dual(T v T(0), T d T(0)) : value(v), deriv(d) {} // 算术运算符重载 - * / 等 Dual operator(const Dual other) const { return Dual(value other.value, deriv other.deriv); } Dual operator*(const Dual other) const { return Dual(value * other.value, deriv * other.value value * other.deriv); } // ... 其他运算符和数学函数sin, cos, exp, log等 }; // 隐函数微分器主类 class ImplicitDifferentiator { public: // 构造函数可以接受一个代表方程的函数对象 templatetypename Func ImplicitDifferentiator(Func f) : func_(f) {} // 核心方法计算在点 p 处隐函数 y 关于自变量 x 的偏导数 // 输入 // - all_vars: 包含所有变量包括自变量和因变量取值的向量。 // - dep_indices: 因变量隐函数在 all_vars 中的索引。 // - indep_indices: 自变量在 all_vars 中的索引。 // 输出一个矩阵第 i 行对应第 i 个因变量第 j 列对应第 j 个自变量即 ∂y_i/∂x_j。 Eigen::MatrixXd computeDerivative(const Eigen::VectorXd all_vars, const std::vectorint dep_indices, const std::vectorint indep_indices); private: std::functiondouble(const Eigen::VectorXd) func_; // 方程 F // 内部会使用自动微分来计算雅可比矩阵 Eigen::MatrixXd computeJacobian(const Eigen::VectorXd point, const std::vectorint wrt_indices); };我们选择使用Eigen库来处理线性代数运算向量、矩阵、求逆因为它是 C 中事实上的标准线性代数库性能优异且接口友好。4. 基于自动微分的核心实现详解现在我们深入代码层面看看如何将理论公式一步步实现。4.1 实现前向模式自动微分Dual类自动微分的前向模式核心是二元数Dual Number算术。对于一个函数f(x)我们不仅计算f(x)还同时计算f(x)。规则来源于导数的基本运算法则。templatetypename T class Dual { public: T val; // 值部分 T grad; // 梯度导数部分 Dual(const T v T(0), const T g T(0)) : val(v), grad(g) {} // 基本算术运算 Dual operator(const Dual other) const { return Dual(val other.val, grad other.grad); } Dual operator-(const Dual other) const { return Dual(val - other.val, grad - other.grad); } Dual operator*(const Dual other) const { // 乘积法则: (uv) uv uv return Dual(val * other.val, grad * other.val val * other.grad); } Dual operator/(const Dual other) const { // 商法则: (u/v) (uv - uv) / v^2 T denom other.val * other.val; return Dual(val / other.val, (grad * other.val - val * other.grad) / denom); } // 与常数的运算 Dual operator(const T scalar) const { return Dual(val scalar, grad); } Dual operator-(const T scalar) const { return Dual(val - scalar, grad); } Dual operator*(const T scalar) const { return Dual(val * scalar, grad * scalar); } Dual operator/(const T scalar) const { return Dual(val / scalar, grad / scalar); } friend Dual operator(const T scalar, const Dual d) { return d scalar; } friend Dual operator*(const T scalar, const Dual d) { return d * scalar; } // ... 其他友元函数 // 初等函数重载 friend Dual sin(const Dual d) { return Dual(std::sin(d.val), std::cos(d.val) * d.grad); // sin(x)cos(x) } friend Dual cos(const Dual d) { return Dual(std::cos(d.val), -std::sin(d.val) * d.grad); // cos(x)-sin(x) } friend Dual exp(const Dual d) { T exp_val std::exp(d.val); return Dual(exp_val, exp_val * d.grad); // exp(x)exp(x) } friend Dual log(const Dual d) { return Dual(std::log(d.val), d.grad / d.val); // log(x)1/x } // ... sqrt, pow 等 };关键技巧如何计算多变量函数的梯度假设我们要计算f(x, y, z)在点(a, b, c)处关于x的偏导数。我们构造Dual变量x Dual(a, 1.0),y Dual(b, 0.0),z Dual(c, 0.0)。这里x.grad1表示我们对x求导y.grad0和z.grad0表示y, z视为常数。然后计算f(x, y, z)得到的Dual结果的grad成员就是∂f/∂x在(a,b,c)处的值。要计算关于y的偏导只需将y的grad设为 1其他设为 0。4.2 计算雅可比矩阵对于一个向量值函数F: R^n - R^m其雅可比矩阵J是一个m x n的矩阵其中J[i][j] ∂F_i / ∂x_j。利用前向模式 AD我们可以一次计算一列即关于一个变量x_j的所有偏导。Eigen::MatrixXd ImplicitDifferentiator::computeJacobian( const Eigen::VectorXd point, // 输入点长度为 n const std::vectorint wrt_indices) // 需要对其求导的变量索引 { int m 1; // 假设 F 是标量函数 (m1)。如果是向量函数需要知道其维度。 int n point.size(); int k wrt_indices.size(); // 实际求导的变量个数 Eigen::MatrixXd jacobian(m, k); std::vectorDualdouble dual_vars(n); // 首先将所有变量初始化为普通值导数为0 for (int i 0; i n; i) { dual_vars[i] Dualdouble(point(i), 0.0); } // 对每个需要求导的变量进行一次前向传播 for (int col 0; col k; col) { int var_idx wrt_indices[col]; // 将当前变量的导数分量设为1 dual_vars[var_idx].grad 1.0; // 使用 dual_vars 计算函数 F 的值现在是 Dual 类型 Dualdouble f_dual func_(dual_vars); // 需要 func_ 能处理 Dual 向量 // 雅可比矩阵的这一列就是 f_dual.grad (对于标量F) jacobian(0, col) f_dual.grad; // 重置当前变量的导数为0为下一个变量做准备 dual_vars[var_idx].grad 0.0; } return jacobian; }注意上面的代码假设func_是一个标量函数。如果F是向量函数m 1那么func_需要返回一个Eigen::VectorXDualdouble并且我们需要遍历这个向量的每个分量F_i取出其grad来填充雅可比矩阵的第i行。为了清晰我们先处理标量方程F0的情形。4.3 组装并求解隐函数导数有了计算雅可比矩阵的能力实现隐函数定理的公式就水到渠成了。Eigen::MatrixXd ImplicitDifferentiator::computeDerivative( const Eigen::VectorXd all_vars, const std::vectorint dep_indices, // 隐变量索引例如 [1] 表示 y const std::vectorint indep_indices) // 自变量索引例如 [0] 表示 x { // 1. 计算雅可比矩阵 J_Fy (F 对因变量 y 的偏导) // 这里 all_vars 的顺序是 [x, y, ...]我们需要计算在 all_vars 点处F 对 dep_indices 的偏导 Eigen::MatrixXd J_Fy computeJacobian(all_vars, dep_indices); // 对于标量FJ_Fy 是一个 1 x (dep_num) 的行向量。 // 但在多元方程组中F和y都是向量J_Fy是方阵。 // 2. 计算雅可比矩阵 J_Fx (F 对自变量 x 的偏导) Eigen::MatrixXd J_Fx computeJacobian(all_vars, indep_indices); // 对于标量FJ_Fx 是一个 1 x (indep_num) 的行向量。 // 3. 应用公式: dy/dx - (J_Fy)^{-1} * J_Fx // 对于标量F单隐函数J_Fy 是 1x1 矩阵即标量求逆就是取倒数。 // 对于标量F多隐函数不可能方程数应等于隐函数数。 // 我们假设是方程组即 F 和 y 维度相同。 int m dep_indices.size(); int n indep_indices.size(); // 检查 J_Fy 是否为方阵且可逆行列式不为零 if (J_Fy.rows() ! J_Fy.cols()) { throw std::invalid_argument(Number of equations must equal number of dependent variables for implicit function theorem.); } double det J_Fy.determinant(); if (std::fabs(det) 1e-10) { // 使用一个小的阈值 throw std::runtime_error(Jacobian w.r.t. dependent variables is singular (det ~ 0). Implicit function theorem condition violated.); } // 计算导数矩阵 Eigen::MatrixXd derivative -J_Fy.inverse() * J_Fx; // 维度: m x n return derivative; }4.4 完整示例计算一个具体隐函数的导数让我们用代码实现开头的单位圆例子。我们想求由F(x, y) x^2 y^2 - 1 0确定的隐函数y(x)在点(0.6, 0.8)处的导数dy/dx。解析解为dy/dx -x/y -0.6/0.8 -0.75。首先我们需要用Dual类型来定义函数F。#include iostream #include Eigen/Dense #include vector #include cmath // ... 此处插入 Dual 类和 ImplicitDifferentiator 类的实现 ... // 定义方程 F(x, y) x^2 y^2 - 1 templatetypename T T circle_equation(const Eigen::MatrixT, Eigen::Dynamic, 1 vars) { // vars[0] x, vars[1] y T x vars(0); T y vars(1); return x * x y * y - T(1.0); } int main() { // 1. 创建微分器传入方程函数 auto func [](const auto v) { return circle_equation(v); }; ImplicitDifferentiator diff(func); // 2. 定义点 (x0, y0) (0.6, 0.8) Eigen::VectorXd point(2); point 0.6, 0.8; // 3. 指定变量索引因变量 y 索引为 1自变量 x 索引为 0 std::vectorint dep_indices {1}; // y std::vectorint indep_indices {0}; // x // 4. 计算导数 dy/dx try { Eigen::MatrixXd dydx diff.computeDerivative(point, dep_indices, indep_indices); std::cout At point ( point(0) , point(1) )\n; std::cout dy/dx computed dydx(0,0) std::endl; std::cout dy/dx exact -point(0)/point(1) std::endl; } catch (const std::exception e) { std::cerr Error: e.what() std::endl; } return 0; }运行这段代码输出应该非常接近At point (0.6, 0.8) dy/dx computed -0.75 dy/dx exact -0.755. 高级应用、性能优化与边界情况处理一个基础的实现已经完成但要投入实际应用我们还需要考虑更多。5.1 处理多元方程组一个完整案例考虑一个更复杂的系统它定义了u(x, y)和v(x, y)F(u, v, x, y) u*x v*y - exp(x) 0 G(u, v, x, y) u*y - v*x - sin(y) 0我们想在点(x, y, u, v) (1.0, 2.0, ? , ?)附近计算∂u/∂x,∂v/∂x等。首先需要找到满足方程组的u, v可能需要用牛顿法求解假设我们已求得u≈0.8, v≈0.6。// 定义方程组 F_vec [F1, F2]^T templatetypename T Eigen::MatrixT, 2, 1 example_system(const Eigen::MatrixT, Eigen::Dynamic, 1 vars) { T u vars(0); T v vars(1); T x vars(2); T y vars(3); Eigen::MatrixT, 2, 1 result; result(0) u * x v * y - exp(x); result(1) u * y - v * x - sin(y); return result; } int main() { // 假设在 (x1.0, y2.0) 处解出 (u≈0.8, v≈0.6) Eigen::VectorXd point(4); point 0.8, 0.6, 1.0, 2.0; // 顺序: [u, v, x, y] auto func [](const auto v) { return example_system(v); }; ImplicitDifferentiator diff(func); // 因变量是 u, v (索引 0, 1)自变量是 x, y (索引 2, 3) std::vectorint dep_indices {0, 1}; std::vectorint indep_indices {2, 3}; Eigen::MatrixXd derivative diff.computeDerivative(point, dep_indices, indep_indices); // derivative 是一个 2x2 矩阵: // [ ∂u/∂x ∂u/∂y ] // [ ∂v/∂x ∂v/∂y ] std::cout Derivative matrix (dy/dx):\n derivative std::endl; return 0; }我们的computeJacobian和computeDerivative函数需要稍作修改以支持向量值函数F但核心逻辑不变分别计算J_Fy(2x2) 和J_Fx(2x2)然后计算-J_Fy^{-1} * J_Fx。5.2 性能优化策略避免矩阵求逆对于线性系统A * X B直接求逆X A^{-1} B在计算上不是最优的尤其是当A很大时。应使用矩阵分解法如 LU、QR、Cholesky。// 将 derivative -J_Fy.inverse() * J_Fx; // 替换为 Eigen::MatrixXd derivative -J_Fy.lu().solve(J_Fx); // 使用LU分解求解这通常更快、更数值稳定。稀疏性利用许多工程问题中的雅可比矩阵是稀疏的很多零元素。使用Eigen的稀疏矩阵模块SparseMatrix和对应的求解器可以极大节省内存和计算时间。自动微分模式选择对于计算一个标量函数对大量输入变量的梯度n很大前向模式需要n次计算。此时反向模式自动微分Reverse Mode AD通常更高效它只需要两次计算一次前向计算值一次反向传播梯度。但这需要记录计算图实现更复杂。可以考虑集成现有的 AD 库如CppAD或Adept。5.3 常见陷阱与排查指南“雅可比矩阵奇异”错误这是最常遇到的问题意味着det(J_Fy) ≈ 0。原因不满足隐函数定理的条件。在求解点附近隐函数可能不存在或不唯一如圆锥的顶点。排查检查你提供的点(x0, y0)是否确实满足方程F(x0, y0)0即使近似满足如果偏差大也可能在“错误”的点计算。检查方程和变量划分是否正确。你是否错误地将一个自变量指定为了因变量尝试对变量进行微小扰动或者换一个初始点重新计算。可能你刚好位于一个奇点。自动微分的精度问题确保所有运算都使用Dual类型。在定义方程F时必须使用重载过的数学函数如std::sin要替换为对Dual重载的sin。混用普通double会导致导数信息丢失。小心函数定义域。例如在log(x)中如果x.val接近或小于0不仅函数值无定义导数也会出问题。需要在计算前验证输入点。维度不匹配错误确保dep_indices的大小等于方程F的维度对于方程组。一个标量方程只能确定一个隐函数。确保all_vars向量的长度等于总变量数自变量因变量。数值稳定性当J_Fy的条件数很大时接近奇异求逆或求解线性系统会放大舍入误差。可以使用更稳定的分解如 SVD或者添加正则化项Tikhonov 正则化来求解(J_Fy^T * J_Fy λI) * X -J_Fy^T * J_Fx其中λ是一个小的正数。6. 在科学与工程中的实际应用场景隐函数求导不仅仅是数学练习它在许多领域是至关重要的工具。物理仿真与多体动力学在机器人学或游戏物理中物体常受到约束如关节连接、接触点。这些约束通常表示为C(q) 0其中q是广义坐标。为了计算速度dq/dt和加速度d²q/dt²需要对约束方程求导dC/dt J * dq/dt 0这直接用到隐函数定理。我们的代码可以用于计算约束雅可比矩阵J。优化问题中的等式约束在拉格朗日乘子法中处理等式约束g(x)0时需要计算约束的梯度。如果约束本身是隐式定义的就需要隐函数求导。在序列二次规划等算法中这很常见。计算机图形学中的隐式曲面隐式曲面由F(x, y, z)0定义如 Metaballs。要渲染或对其进行操作需要计算法向量而法向量正是梯度∇F (∂F/∂x, ∂F/∂y, ∂F/∂z)。如果F本身很复杂自动微分就能派上用场。经济学与均衡模型市场均衡由一组供需方程共同决定形成隐函数关系。分析政策变化外生变量对均衡价格和数量内生变量的影响即计算比较静态导数正是隐函数定理的典型应用。机器学习中的微分方程参数学习在物理信息神经网络中损失函数可能包含满足微分方程G(u, θ)0的解u。为了通过梯度下降优化参数θ需要计算du/dθ这又是一个隐函数求导问题。最后一点个人体会实现这个工具的过程让我对“导数”有了更深的理解。自动微分将求导从符号演算和有限差分近似变成了一个精确、可编程的机械过程。而隐函数定理则提供了一套框架让我们能在不知道函数显式形式的情况下仅通过其“零水平集”来操控其微分性质。将两者结合你获得的能力是对于任何能用代码描述的、满足一定光滑性的隐式关系你都可以在运行时精确地获取其局部微分信息。这种能力在构建高度动态和自适应的数值系统时是无价的。刚开始实现时我纠结于抽象的矩阵求逆公式但当我看到第一个单位圆例子正确输出-0.75时那种理论落地成代码的实在感是单纯推公式无法比拟的。建议你在理解原理后一定要亲手实现一遍哪怕从一个简单的Dual类开始遇到的每一个编译错误和数值异常都会让你对隐函数定理的条件和自动微分的细节有更刻骨铭心的认识。
1. 项目概述当方程无法显式求解时我们如何求导在工程优化、物理模拟和机器学习等领域我们常常会遇到这样的问题一个变量依赖于另一个变量但它们之间的关系并非一个简单的y f(x)公式而是被一个复杂的方程F(x, y) 0所隐含地定义。比如一个机械臂的末端位置(x, y)可能满足一个包含关节角度的复杂三角方程或者在经济学模型中均衡价格P和数量Q由市场供需方程组共同决定无法直接解出P f(Q)。这时候如果你想分析y如何随x变化即求dy/dx或者需要计算梯度进行数值优化直接求导的路就走不通了因为你连y的显式表达式都没有。隐函数定理Implicit Function Theorem就是解决这类问题的“尚方宝剑”。它不关心y具体等于什么而是直接告诉你在满足一定条件下y作为x的函数是存在的、可微的并且其导数可以通过原方程F的偏导数直接计算出来。这个项目就是要把这把理论上的“宝剑”打磨成一把可以在 C 代码中挥舞的实用“利刃”。我们将从最核心的一元隐函数求导公式dy/dx -F_x / F_y出发一直深入到多元方程组情形下的雅可比矩阵求逆运算并最终用 C 实现一个通用的、支持自动微分的隐函数偏导数计算模块。无论你是正在学习《高等数学》或《数学分析》的学生还是需要在仿真、优化算法中处理隐式约束的工程师这篇文章都将带你从理论到实践彻底掌握隐函数求导及其 C 实现的核心要领。2. 隐函数定理的核心思想与公式推导理解隐函数定理关键在于抓住其“局部存在性”和“导数可计算性”这两个核心。它不像代数求解那样追求全局的、显式的解而是保证在某个特定的点附近解是存在的、唯一的并且其微分行为完全由原函数在该点的局部性质决定。2.1 从直观例子理解定理条件让我们先看一个经典的例子单位圆方程x^2 y^2 - 1 0。我们想知道在点(0, 1)即圆的最顶端附近y能否表示为x的函数以及其导数。连续性条件C¹方程F(x, y) x^2 y^2 - 1本身及其偏导数F_x 2x,F_y 2y在(0, 1)附近都是连续且光滑的。这保证了函数行为是“良好”的没有奇点或突变。零点条件显然F(0, 1) 0点(0, 1)在圆上。非退化条件计算F_y(0, 1) 2 * 1 2 ≠ 0。这是最关键的条件。F_y ≠ 0意味着在(0, 1)点方程F(x, y)0在y方向上是“非奇异”的y的变化对F的值有显著影响。从几何上看圆在该点的切线不是竖直的。当这三个条件满足时隐函数定理断言在(0, 1)附近的一个小邻域内存在唯一的函数y φ(x)使得F(x, φ(x)) ≡ 0。并且这个隐函数φ(x)也是连续可微的。注意这个“局部”特性非常重要。对于整个圆来说y不能写成x的单一函数因为一个x对应两个y。但在(0, 1)这个点的附近比如右上半圆的一小段它是可以的。同样在(1, 0)点F_y 0定理条件不满足此时圆在该点的切线是竖直的y无法表示为x的函数但x可以表示为y的函数。2.2 核心公式的推导链式法则的巧妙应用定理不仅保证了隐函数存在还给出了计算其导数的公式。推导过程简洁而优美是链式法则的经典应用。假设由F(x, y) 0确定了隐函数y φ(x)且F和φ都可微。我们将方程F(x, φ(x)) 0的两边对x求导。这里F是一个二元函数x是直接变量y通过φ(x)也依赖于x。应用多元函数的链式法则d/dx [ F(x, φ(x)) ] (∂F/∂x) * (dx/dx) (∂F/∂y) * (dφ/dx) 0由于dx/dx 1上式简化为F_x(x, φ(x)) F_y(x, φ(x)) * φ(x) 0只要F_y ≠ 0我们就可以立即解出φ(x) - F_x(x, φ(x)) / F_y(x, φ(x))这就是一元隐函数求导的万能公式。你不需要知道φ(x)的具体形式只需要知道原函数F和该点的坐标(x, y)就能算出隐函数的导数。实操心得这个推导过程揭示了隐函数求导的本质——对方程两端同时关于自变量求全微分。对于更复杂的多元隐函数或方程组这个思想依然适用只是从标量方程升级到了矩阵方程。2.3 推广到多元隐函数与方程组实际问题往往是多维的。例如由方程组F(u, v, x, y) 0 G(u, v, x, y) 0确定了两个隐函数u u(x, y),v v(x, y)。我们想求∂u/∂x,∂v/∂x等。此时我们将所有方程对x求偏导视y为常数再次利用链式法则∂F/∂x (∂F/∂u)*(∂u/∂x) (∂F/∂v)*(∂v/∂x) (∂F/∂x)*(1) (∂F/∂y)*(0) 0 ∂G/∂x (∂G/∂u)*(∂u/∂x) (∂G/∂v)*(∂v/∂x) (∂G/∂x)*(1) (∂G/∂y)*(0) 0注意这里的∂F/∂x是指F对第一个显式变量x的偏导而∂u/∂x才是我们要求的隐函数偏导。整理后得到一个关于∂u/∂x和∂v/∂x的线性方程组[ ∂F/∂u ∂F/∂v ] [ ∂u/∂x ] [ -∂F/∂x ] [ ∂G/∂u ∂G/∂v ] [ ∂v/∂x ] [ -∂G/∂x ]左边的矩阵正是雅可比矩阵J_{(u,v)}F函数组(F, G)对变量组(u, v)的偏导数矩阵。只要这个雅可比矩阵在求解点可逆即行列式非零这是高维下的“非退化条件”我们就可以解出[ ∂u/∂x ] -1 [ -∂F/∂x ] [ ∂v/∂x ] [J_{(u,v)}F] * [ -∂G/∂x ]这就是隐函数定理的矩阵形式。它告诉我们隐函数组的导数等于原函数组对隐变量雅可比矩阵的逆乘以原函数组对显式变量的偏导数的负值。3. C实现的核心架构与设计选择将数学定理转化为代码我们需要做出几个关键的设计决策。我们的目标是构建一个灵活、高效且易于使用的ImplicitDifferentiator类。3.1 核心计算流程设计整个计算流程可以分解为以下几步这构成了我们类方法的骨架定义方程与变量用户提供方程F(x, y, ...) 0并指定哪些是自变量x哪些是因变量y即隐函数。计算偏导数计算F对所有变量的偏导数。这是最核心也最耗时的步骤。构造雅可比矩阵根据因变量和自变量的划分从偏导数集合中分别提取出J_yF对因变量的雅可比和J_xF对自变量的偏导向量/矩阵。求解线性系统计算- (J_yF)^{-1} * J_xF。对于单方程单隐函数这就是一个除法对于方程组这就是一个矩阵求逆或求解线性方程组的问题。输出结果返回隐函数关于自变量的偏导数标量、向量或矩阵。3.2 关键设计决策如何表示和计算偏导数这里有三个主流方案各有优劣方案一符号微分Symbolic Differentiation思路像 Mathematica 或 SymPy 一样解析地推导出偏导数的表达式。C实现需要构建表达式树实现符号化简和求导规则。可以使用像SymEngine这样的库。优点求导精确一次推导后可重复使用表达式对于复杂函数可能更快。缺点实现极其复杂表达式膨胀可能导致性能下降且难以处理控制流复杂的函数。结论对于专注于隐函数求导的专用工具来说过于重型不推荐。方案二数值微分Numerical Differentiation思路用有限差分近似偏导数例如∂F/∂x ≈ (F(xh, y) - F(x, y)) / h。C实现非常简单只需要能调用函数F即可。优点实现简单对函数F的形式没有任何要求即使是黑箱函数。缺点存在截断误差和舍入误差选择步长h需要技巧太小会放大舍入误差太大则截断误差大。对于条件数大的雅可比矩阵误差会被放大可能导致求解不稳定。结论适合快速原型验证或当F是黑箱函数时使用。但作为通用库的核心精度和稳定性是短板。方案三自动微分Automatic Differentiation, AD思路介于符号和数值之间。通过操作符重载在计算函数值的同时精确地计算其导数。它计算的是精确的导数值没有截断误差。C实现需要定义一个特殊的“双数”类型例如Dualdouble重载所有算术和初等函数操作。优点计算精度高机器精度效率通常优于符号微分且易于使用。缺点需要将函数F用重载了 AD 的类型重写。有轻微的计算和内存开销。结论这是本项目的最佳选择。它在精度、效率和易用性之间取得了最佳平衡。现代科学计算库如 Ceres Solver, Stan的核心都依赖于 AD。我们选择使用前向模式自动微分Forward Mode AD来实现。对于计算一个标量函数对多个变量的梯度前向模式非常直观高效。我们将实现一个简单的Dual类作为基础。3.3 类与接口设计基于以上分析我们设计以下核心类// 前向模式自动微分的双数类 templatetypename T class Dual { public: T value; // 函数值 T deriv; // 导数值 Dual(T v T(0), T d T(0)) : value(v), deriv(d) {} // 算术运算符重载 - * / 等 Dual operator(const Dual other) const { return Dual(value other.value, deriv other.deriv); } Dual operator*(const Dual other) const { return Dual(value * other.value, deriv * other.value value * other.deriv); } // ... 其他运算符和数学函数sin, cos, exp, log等 }; // 隐函数微分器主类 class ImplicitDifferentiator { public: // 构造函数可以接受一个代表方程的函数对象 templatetypename Func ImplicitDifferentiator(Func f) : func_(f) {} // 核心方法计算在点 p 处隐函数 y 关于自变量 x 的偏导数 // 输入 // - all_vars: 包含所有变量包括自变量和因变量取值的向量。 // - dep_indices: 因变量隐函数在 all_vars 中的索引。 // - indep_indices: 自变量在 all_vars 中的索引。 // 输出一个矩阵第 i 行对应第 i 个因变量第 j 列对应第 j 个自变量即 ∂y_i/∂x_j。 Eigen::MatrixXd computeDerivative(const Eigen::VectorXd all_vars, const std::vectorint dep_indices, const std::vectorint indep_indices); private: std::functiondouble(const Eigen::VectorXd) func_; // 方程 F // 内部会使用自动微分来计算雅可比矩阵 Eigen::MatrixXd computeJacobian(const Eigen::VectorXd point, const std::vectorint wrt_indices); };我们选择使用Eigen库来处理线性代数运算向量、矩阵、求逆因为它是 C 中事实上的标准线性代数库性能优异且接口友好。4. 基于自动微分的核心实现详解现在我们深入代码层面看看如何将理论公式一步步实现。4.1 实现前向模式自动微分Dual类自动微分的前向模式核心是二元数Dual Number算术。对于一个函数f(x)我们不仅计算f(x)还同时计算f(x)。规则来源于导数的基本运算法则。templatetypename T class Dual { public: T val; // 值部分 T grad; // 梯度导数部分 Dual(const T v T(0), const T g T(0)) : val(v), grad(g) {} // 基本算术运算 Dual operator(const Dual other) const { return Dual(val other.val, grad other.grad); } Dual operator-(const Dual other) const { return Dual(val - other.val, grad - other.grad); } Dual operator*(const Dual other) const { // 乘积法则: (uv) uv uv return Dual(val * other.val, grad * other.val val * other.grad); } Dual operator/(const Dual other) const { // 商法则: (u/v) (uv - uv) / v^2 T denom other.val * other.val; return Dual(val / other.val, (grad * other.val - val * other.grad) / denom); } // 与常数的运算 Dual operator(const T scalar) const { return Dual(val scalar, grad); } Dual operator-(const T scalar) const { return Dual(val - scalar, grad); } Dual operator*(const T scalar) const { return Dual(val * scalar, grad * scalar); } Dual operator/(const T scalar) const { return Dual(val / scalar, grad / scalar); } friend Dual operator(const T scalar, const Dual d) { return d scalar; } friend Dual operator*(const T scalar, const Dual d) { return d * scalar; } // ... 其他友元函数 // 初等函数重载 friend Dual sin(const Dual d) { return Dual(std::sin(d.val), std::cos(d.val) * d.grad); // sin(x)cos(x) } friend Dual cos(const Dual d) { return Dual(std::cos(d.val), -std::sin(d.val) * d.grad); // cos(x)-sin(x) } friend Dual exp(const Dual d) { T exp_val std::exp(d.val); return Dual(exp_val, exp_val * d.grad); // exp(x)exp(x) } friend Dual log(const Dual d) { return Dual(std::log(d.val), d.grad / d.val); // log(x)1/x } // ... sqrt, pow 等 };关键技巧如何计算多变量函数的梯度假设我们要计算f(x, y, z)在点(a, b, c)处关于x的偏导数。我们构造Dual变量x Dual(a, 1.0),y Dual(b, 0.0),z Dual(c, 0.0)。这里x.grad1表示我们对x求导y.grad0和z.grad0表示y, z视为常数。然后计算f(x, y, z)得到的Dual结果的grad成员就是∂f/∂x在(a,b,c)处的值。要计算关于y的偏导只需将y的grad设为 1其他设为 0。4.2 计算雅可比矩阵对于一个向量值函数F: R^n - R^m其雅可比矩阵J是一个m x n的矩阵其中J[i][j] ∂F_i / ∂x_j。利用前向模式 AD我们可以一次计算一列即关于一个变量x_j的所有偏导。Eigen::MatrixXd ImplicitDifferentiator::computeJacobian( const Eigen::VectorXd point, // 输入点长度为 n const std::vectorint wrt_indices) // 需要对其求导的变量索引 { int m 1; // 假设 F 是标量函数 (m1)。如果是向量函数需要知道其维度。 int n point.size(); int k wrt_indices.size(); // 实际求导的变量个数 Eigen::MatrixXd jacobian(m, k); std::vectorDualdouble dual_vars(n); // 首先将所有变量初始化为普通值导数为0 for (int i 0; i n; i) { dual_vars[i] Dualdouble(point(i), 0.0); } // 对每个需要求导的变量进行一次前向传播 for (int col 0; col k; col) { int var_idx wrt_indices[col]; // 将当前变量的导数分量设为1 dual_vars[var_idx].grad 1.0; // 使用 dual_vars 计算函数 F 的值现在是 Dual 类型 Dualdouble f_dual func_(dual_vars); // 需要 func_ 能处理 Dual 向量 // 雅可比矩阵的这一列就是 f_dual.grad (对于标量F) jacobian(0, col) f_dual.grad; // 重置当前变量的导数为0为下一个变量做准备 dual_vars[var_idx].grad 0.0; } return jacobian; }注意上面的代码假设func_是一个标量函数。如果F是向量函数m 1那么func_需要返回一个Eigen::VectorXDualdouble并且我们需要遍历这个向量的每个分量F_i取出其grad来填充雅可比矩阵的第i行。为了清晰我们先处理标量方程F0的情形。4.3 组装并求解隐函数导数有了计算雅可比矩阵的能力实现隐函数定理的公式就水到渠成了。Eigen::MatrixXd ImplicitDifferentiator::computeDerivative( const Eigen::VectorXd all_vars, const std::vectorint dep_indices, // 隐变量索引例如 [1] 表示 y const std::vectorint indep_indices) // 自变量索引例如 [0] 表示 x { // 1. 计算雅可比矩阵 J_Fy (F 对因变量 y 的偏导) // 这里 all_vars 的顺序是 [x, y, ...]我们需要计算在 all_vars 点处F 对 dep_indices 的偏导 Eigen::MatrixXd J_Fy computeJacobian(all_vars, dep_indices); // 对于标量FJ_Fy 是一个 1 x (dep_num) 的行向量。 // 但在多元方程组中F和y都是向量J_Fy是方阵。 // 2. 计算雅可比矩阵 J_Fx (F 对自变量 x 的偏导) Eigen::MatrixXd J_Fx computeJacobian(all_vars, indep_indices); // 对于标量FJ_Fx 是一个 1 x (indep_num) 的行向量。 // 3. 应用公式: dy/dx - (J_Fy)^{-1} * J_Fx // 对于标量F单隐函数J_Fy 是 1x1 矩阵即标量求逆就是取倒数。 // 对于标量F多隐函数不可能方程数应等于隐函数数。 // 我们假设是方程组即 F 和 y 维度相同。 int m dep_indices.size(); int n indep_indices.size(); // 检查 J_Fy 是否为方阵且可逆行列式不为零 if (J_Fy.rows() ! J_Fy.cols()) { throw std::invalid_argument(Number of equations must equal number of dependent variables for implicit function theorem.); } double det J_Fy.determinant(); if (std::fabs(det) 1e-10) { // 使用一个小的阈值 throw std::runtime_error(Jacobian w.r.t. dependent variables is singular (det ~ 0). Implicit function theorem condition violated.); } // 计算导数矩阵 Eigen::MatrixXd derivative -J_Fy.inverse() * J_Fx; // 维度: m x n return derivative; }4.4 完整示例计算一个具体隐函数的导数让我们用代码实现开头的单位圆例子。我们想求由F(x, y) x^2 y^2 - 1 0确定的隐函数y(x)在点(0.6, 0.8)处的导数dy/dx。解析解为dy/dx -x/y -0.6/0.8 -0.75。首先我们需要用Dual类型来定义函数F。#include iostream #include Eigen/Dense #include vector #include cmath // ... 此处插入 Dual 类和 ImplicitDifferentiator 类的实现 ... // 定义方程 F(x, y) x^2 y^2 - 1 templatetypename T T circle_equation(const Eigen::MatrixT, Eigen::Dynamic, 1 vars) { // vars[0] x, vars[1] y T x vars(0); T y vars(1); return x * x y * y - T(1.0); } int main() { // 1. 创建微分器传入方程函数 auto func [](const auto v) { return circle_equation(v); }; ImplicitDifferentiator diff(func); // 2. 定义点 (x0, y0) (0.6, 0.8) Eigen::VectorXd point(2); point 0.6, 0.8; // 3. 指定变量索引因变量 y 索引为 1自变量 x 索引为 0 std::vectorint dep_indices {1}; // y std::vectorint indep_indices {0}; // x // 4. 计算导数 dy/dx try { Eigen::MatrixXd dydx diff.computeDerivative(point, dep_indices, indep_indices); std::cout At point ( point(0) , point(1) )\n; std::cout dy/dx computed dydx(0,0) std::endl; std::cout dy/dx exact -point(0)/point(1) std::endl; } catch (const std::exception e) { std::cerr Error: e.what() std::endl; } return 0; }运行这段代码输出应该非常接近At point (0.6, 0.8) dy/dx computed -0.75 dy/dx exact -0.755. 高级应用、性能优化与边界情况处理一个基础的实现已经完成但要投入实际应用我们还需要考虑更多。5.1 处理多元方程组一个完整案例考虑一个更复杂的系统它定义了u(x, y)和v(x, y)F(u, v, x, y) u*x v*y - exp(x) 0 G(u, v, x, y) u*y - v*x - sin(y) 0我们想在点(x, y, u, v) (1.0, 2.0, ? , ?)附近计算∂u/∂x,∂v/∂x等。首先需要找到满足方程组的u, v可能需要用牛顿法求解假设我们已求得u≈0.8, v≈0.6。// 定义方程组 F_vec [F1, F2]^T templatetypename T Eigen::MatrixT, 2, 1 example_system(const Eigen::MatrixT, Eigen::Dynamic, 1 vars) { T u vars(0); T v vars(1); T x vars(2); T y vars(3); Eigen::MatrixT, 2, 1 result; result(0) u * x v * y - exp(x); result(1) u * y - v * x - sin(y); return result; } int main() { // 假设在 (x1.0, y2.0) 处解出 (u≈0.8, v≈0.6) Eigen::VectorXd point(4); point 0.8, 0.6, 1.0, 2.0; // 顺序: [u, v, x, y] auto func [](const auto v) { return example_system(v); }; ImplicitDifferentiator diff(func); // 因变量是 u, v (索引 0, 1)自变量是 x, y (索引 2, 3) std::vectorint dep_indices {0, 1}; std::vectorint indep_indices {2, 3}; Eigen::MatrixXd derivative diff.computeDerivative(point, dep_indices, indep_indices); // derivative 是一个 2x2 矩阵: // [ ∂u/∂x ∂u/∂y ] // [ ∂v/∂x ∂v/∂y ] std::cout Derivative matrix (dy/dx):\n derivative std::endl; return 0; }我们的computeJacobian和computeDerivative函数需要稍作修改以支持向量值函数F但核心逻辑不变分别计算J_Fy(2x2) 和J_Fx(2x2)然后计算-J_Fy^{-1} * J_Fx。5.2 性能优化策略避免矩阵求逆对于线性系统A * X B直接求逆X A^{-1} B在计算上不是最优的尤其是当A很大时。应使用矩阵分解法如 LU、QR、Cholesky。// 将 derivative -J_Fy.inverse() * J_Fx; // 替换为 Eigen::MatrixXd derivative -J_Fy.lu().solve(J_Fx); // 使用LU分解求解这通常更快、更数值稳定。稀疏性利用许多工程问题中的雅可比矩阵是稀疏的很多零元素。使用Eigen的稀疏矩阵模块SparseMatrix和对应的求解器可以极大节省内存和计算时间。自动微分模式选择对于计算一个标量函数对大量输入变量的梯度n很大前向模式需要n次计算。此时反向模式自动微分Reverse Mode AD通常更高效它只需要两次计算一次前向计算值一次反向传播梯度。但这需要记录计算图实现更复杂。可以考虑集成现有的 AD 库如CppAD或Adept。5.3 常见陷阱与排查指南“雅可比矩阵奇异”错误这是最常遇到的问题意味着det(J_Fy) ≈ 0。原因不满足隐函数定理的条件。在求解点附近隐函数可能不存在或不唯一如圆锥的顶点。排查检查你提供的点(x0, y0)是否确实满足方程F(x0, y0)0即使近似满足如果偏差大也可能在“错误”的点计算。检查方程和变量划分是否正确。你是否错误地将一个自变量指定为了因变量尝试对变量进行微小扰动或者换一个初始点重新计算。可能你刚好位于一个奇点。自动微分的精度问题确保所有运算都使用Dual类型。在定义方程F时必须使用重载过的数学函数如std::sin要替换为对Dual重载的sin。混用普通double会导致导数信息丢失。小心函数定义域。例如在log(x)中如果x.val接近或小于0不仅函数值无定义导数也会出问题。需要在计算前验证输入点。维度不匹配错误确保dep_indices的大小等于方程F的维度对于方程组。一个标量方程只能确定一个隐函数。确保all_vars向量的长度等于总变量数自变量因变量。数值稳定性当J_Fy的条件数很大时接近奇异求逆或求解线性系统会放大舍入误差。可以使用更稳定的分解如 SVD或者添加正则化项Tikhonov 正则化来求解(J_Fy^T * J_Fy λI) * X -J_Fy^T * J_Fx其中λ是一个小的正数。6. 在科学与工程中的实际应用场景隐函数求导不仅仅是数学练习它在许多领域是至关重要的工具。物理仿真与多体动力学在机器人学或游戏物理中物体常受到约束如关节连接、接触点。这些约束通常表示为C(q) 0其中q是广义坐标。为了计算速度dq/dt和加速度d²q/dt²需要对约束方程求导dC/dt J * dq/dt 0这直接用到隐函数定理。我们的代码可以用于计算约束雅可比矩阵J。优化问题中的等式约束在拉格朗日乘子法中处理等式约束g(x)0时需要计算约束的梯度。如果约束本身是隐式定义的就需要隐函数求导。在序列二次规划等算法中这很常见。计算机图形学中的隐式曲面隐式曲面由F(x, y, z)0定义如 Metaballs。要渲染或对其进行操作需要计算法向量而法向量正是梯度∇F (∂F/∂x, ∂F/∂y, ∂F/∂z)。如果F本身很复杂自动微分就能派上用场。经济学与均衡模型市场均衡由一组供需方程共同决定形成隐函数关系。分析政策变化外生变量对均衡价格和数量内生变量的影响即计算比较静态导数正是隐函数定理的典型应用。机器学习中的微分方程参数学习在物理信息神经网络中损失函数可能包含满足微分方程G(u, θ)0的解u。为了通过梯度下降优化参数θ需要计算du/dθ这又是一个隐函数求导问题。最后一点个人体会实现这个工具的过程让我对“导数”有了更深的理解。自动微分将求导从符号演算和有限差分近似变成了一个精确、可编程的机械过程。而隐函数定理则提供了一套框架让我们能在不知道函数显式形式的情况下仅通过其“零水平集”来操控其微分性质。将两者结合你获得的能力是对于任何能用代码描述的、满足一定光滑性的隐式关系你都可以在运行时精确地获取其局部微分信息。这种能力在构建高度动态和自适应的数值系统时是无价的。刚开始实现时我纠结于抽象的矩阵求逆公式但当我看到第一个单位圆例子正确输出-0.75时那种理论落地成代码的实在感是单纯推公式无法比拟的。建议你在理解原理后一定要亲手实现一遍哪怕从一个简单的Dual类开始遇到的每一个编译错误和数值异常都会让你对隐函数定理的条件和自动微分的细节有更刻骨铭心的认识。
![[C++]内存管理：串顺序存储的内存回收](http://pic.xiahunao.cn/yaotu/[C++]内存管理：串顺序存储的内存回收)
