牛顿法优化

牛顿法 (Newton’s Method) 在优化中的应用

1. 概念的直观理解 (The “Why” and “What”)

1.1 这是什么？我们为什么需要它？

在机器学习中，我们经常需要找到一个函数的最小值或最大值。例如，训练模型时，我们通常会定义一个损失函数 (Loss Function) 来衡量模型预测与真实值之间的差距。我们的目标就是找到一组模型参数，使得这个损失函数的值最小。这个寻找最小值（或最大值）的过程就是优化 (Optimization)。

牛顿法就是一种用来寻找函数局部最小值（或最大值）的迭代优化算法。它比我们可能更熟悉的梯度下降法（Gradient Descent）通常收敛得更快，因为它利用了函数的二阶导数信息。

1.2 现实类比和几何直觉

想象一下我们在一座山顶（或者山谷底部），我们想要找到最低点。

• 梯度下降法就像是我们闭着眼睛，感受脚下的坡度（梯度），然后沿着最陡峭的方向迈出一步。这一步可能很小，也可能有点大，但我们不知道最低点到底有多远，也不知道走到最低点之前坡度会如何变化。我们只能一步步摸索。
• 牛顿法则更像是我们拥有一个“超能力”，我们不仅知道脚下的坡度（一阶导数/梯度），我们还能感知到坡度变化的趋势，即“地面的弯曲程度”（二阶导数/Hessian矩阵）。有了这些信息，我们就可以预测最低点大致在哪里，然后直接“跳”过去。如果我们的预测足够准确，我们可能只需要几次跳跃就能到达最低点。

几何直觉：用抛物线逼近

对于一个单变量函数 $f(x)$，牛顿法的核心思想是：在当前点 $x_k$ 附近，我们用一个抛物线 (quadratic function) 来近似 $f(x)$。我们知道抛物线有一个唯一的最低点（或最高点），这个最低点可以直接通过公式计算出来。牛顿法就是把这个抛物线的最低点作为下一次迭代的猜测值 $x_{k+1}$。

想象我们在一个复杂的函数曲面上寻找最低点：

1. 在当前位置 $x_k$ 观察地面。
2. 假设当前位置附近的地面是一个抛物线形状。
3. 直接计算这个抛物线的最低点在哪里。
4. “跳”到这个抛物线的最低点，作为我们的新位置 $x_{k+1}$。
5. 重复这个过程，直到我们找到真正的最低点。

这种“用局部抛物线近似”的方法使得牛顿法能够更快地收敛到最小值，因为它利用了更多的函数信息（不仅仅是坡度，还有坡度变化的速率）。

---

2. 预备数学知识与符号解释 (Mathematical Preliminaries & Notation)

为了理解牛顿法，我们需要回顾一些基础的微积分和线性代数概念。

2.1 导数与偏导数

• 导数 (Derivative)：对于单变量函数 $f(x)$，导数 $f^{\prime }(x)$ 或 $\frac{df}{dx}$ 表示函数在某一点的瞬时变化率（斜率）。如果导数为0，通常表示函数在该点有局部最大值、最小值或鞍点。
• 二阶导数 (Second Derivative)：$f^{\prime \prime }(x)$ 或 $\frac{d^{2}f}{d{x}^2}$ 表示函数变化率的变化率，即函数的弯曲程度。
  • 如果 $f^{\prime \prime }(x)>0$，函数在该点是凹的（碗状向上，有局部最小值）。
  • 如果 $f^{\prime \prime }(x)<0$，函数在该点是凸的（碗状向下，有局部最大值）。
• 偏导数 (Partial Derivative)：对于多变量函数 $f(x_1,x_2,\dots ,x_n)$，偏导数 $\frac{\partial f}{\partial x_i}$ 表示当所有其他变量保持不变时，函数随 $x_i$ 的变化率。

2.2 梯度 (Gradient)

• 符号：$\nabla f(\mathbf{x})$
• 定义：对于多变量函数 $f(\mathbf{x})$，其中 $\mathbf{x}=[x_1,x_2,\dots ,x_n{]}^T$ 是一个向量，梯度是一个向量，其分量是 $f$ 对每个变量的偏导数。 $\nabla f(\mathbf{x})=\left[\begin{array}{c}\frac{\partial f}{\partial x_1}\\ \frac{\partial f}{\partial x_2}\\ ⋮\\ \frac{\partial f}{\partial x_n}\end{array}\right]$
• 意义：梯度向量指向函数值增长最快的方向，其大小表示该方向上的最大变化率。在优化中，我们通常沿着梯度的反方向（最速下降方向）移动来寻找最小值。

2.3 Hessian 矩阵 (Hessian Matrix)

• 符号：$\mathbf{H}(\mathbf{x})$ 或 ${\nabla }^{2}f(\mathbf{x})$
• 定义：对于多变量函数 $f(\mathbf{x})$，Hessian 矩阵是一个方阵，包含了所有的二阶偏导数。其第 $i$ 行第 $j$ 列的元素是 $f$ 对 $x_i$ 和 $x_j$ 的二阶偏导数。 $\mathbf{H}(\mathbf{x}{)}_{ij}=\frac{{\partial }^{2}f}{\partial x_i\partial x_j}$ 即： $\mathbf{H}(\mathbf{x})=\left[\begin{array}{cccc}\frac{{\partial }^{2}f}{\partial x_1^2}& \frac{{\partial }^{2}f}{\partial x_1\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_1\partial x_n}\\ \frac{{\partial }^{2}f}{\partial x_2\partial x_1}& \frac{{\partial }^{2}f}{\partial x_2^2}& \cdots & \frac{{\partial }^{2}f}{\partial x_2\partial x_n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\partial }^{2}f}{\partial x_n\partial x_1}& \frac{{\partial }^{2}f}{\partial x_n\partial x_2}& \cdots & \frac{{\partial }^{2}f}{\partial x_n^2}\end{array}\right]$
• 性质：如果函数 $f$ 的二阶偏导数是连续的，那么混合偏导数是相等的，即 $\frac{{\partial }^{2}f}{\partial x_i\partial x_j}=\frac{{\partial }^{2}f}{\partial x_j\partial x_i}$。这意味着 Hessian 矩阵通常是一个对称矩阵。
• 意义：Hessian 矩阵描述了函数曲率（弯曲程度）的性质。它告诉我们函数在某个方向上是加速增长还是减速增长，这对于判断一个点是局部最小值、最大值还是鞍点至关重要。

2.4 泰勒级数展开 (Taylor Series Expansion)

这是牛顿法的核心数学工具。它允许我们用多项式来近似一个在某点附近足够光滑的函数。

• 单变量函数 $f(x)$ 在 $x_k$ 处的二阶泰勒展开： $f(x)\approx f(x_k)+f^{\prime }(x_k)(x-x_k)+\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2$ 这里，$x$ 是我们想要估计的函数值点，$x_k$ 是我们知道函数值和导数值的点。

• 多变量函数 $f(\mathbf{x})$ 在 ${\mathbf{x}}_k$ 处的二阶泰勒展开： $f(\mathbf{x})\approx f({\mathbf{x}}_k)+\nabla f({\mathbf{x}}_k{)}^T(\mathbf{x}-{\mathbf{x}}_k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_k{)}^T\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)$ 这里，$\mathbf{x}$ 是一个向量，${}^T$ 表示向量或矩阵的转置。

  • $f({\mathbf{x}}_k)$：函数在当前点的值。
  • $\nabla f({\mathbf{x}}_k{)}^T(\mathbf{x}-{\mathbf{x}}_k)$：这是一个点积，表示一阶导数（梯度）项。
  • $\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_k{)}^T\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)$：这是一个二次型，表示二阶导数（Hessian）项。

2.5 矩阵逆 (Matrix Inverse)

• 对于一个方阵 $\mathbf{A}$，如果存在一个矩阵 ${\mathbf{A}}^{-1}$ 使得 $\mathbf{A}{\mathbf{A}}^{-1}={\mathbf{A}}^{-1}\mathbf{A}=\mathbf{I}$（其中 $\mathbf{I}$ 是单位矩阵），则称 ${\mathbf{A}}^{-1}$ 为 $\mathbf{A}$ 的逆矩阵。
• 在牛顿法中，我们需要计算 Hessian 矩阵的逆矩阵。

---

3. 详细分步推导 (Step-by-Step Derivation)

牛顿法的思想是，我们通过找到泰勒展开式的近似函数的最小值来更新我们的参数。对于一个光滑函数，其最小值点（或最大值点）处的导数（或梯度）为零。

3.1 单变量函数 $f(x)$ 的牛顿法推导

我们的目标是找到 $x$ 使得 $f^{\prime }(x)=0$。

步骤 1：二阶泰勒展开 我们用一个抛物线来近似函数 $f(x)$ 在当前点 $x_k$ 附近的行为。这个抛物线就是 $f(x)$ 在 $x_k$ 处的二阶泰勒展开式： $f(x)\approx f(x_k)+f^{\prime }(x_k)(x-x_k)+\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2$ 令 $g(x)=f(x_k)+f^{\prime }(x_k)(x-x_k)+\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2$ 为这个近似函数。

步骤 2：找到近似函数的最小值 要找到 $g(x)$ 的最小值，我们对其求导，并令导数等于零。 我们对 $g(x)$ 关于 $x$ 求一阶导数： $\frac{d}{dx}g(x)=\frac{d}{dx}\left(f(x_k)+f^{\prime }(x_k)(x-x_k)+\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2\right)$ 解释： 逐项求导。

• $f(x_k)$ 是常数项，求导为 0。
• $f^{\prime }(x_k)(x-x_k)$ 是线性项，求导得到 $f^{\prime }(x_k)$ （因为 $f^{\prime }(x_k)$ 是一个常数）。
• $\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2$ 是二次项，求导利用链式法则。令 $u=(x-x_k)$，则 $\frac{d}{dx}{u}^2=2u\frac{du}{dx}=2(x-x_k)\cdot 1$。 所以，$\frac{d}{dx}\left(\frac{1}{2}{f}^{\prime \prime }(x_k)(x-x_k{)}^2\right)=\frac{1}{2}{f}^{\prime \prime }(x_k)\cdot 2(x-x_k)=f^{\prime \prime }(x_k)(x-x_k)$。

结合以上，我们得到 $g(x)$ 的导数： $\frac{d}{dx}g(x)=f^{\prime }(x_k)+f^{\prime \prime }(x_k)(x-x_k)$

为了找到 $g(x)$ 的最小值点（我们称之为 $x_{k+1}$），我们将这个导数设为零： $f^{\prime }(x_k)+f^{\prime \prime }(x_k)(x_{k+1}-x_k)=0$

步骤 3：解出 $x_{k+1}$ 现在我们从上面的方程中解出 $x_{k+1}$。 $f^{\prime \prime }(x_k)(x_{k+1}-x_k)=-f^{\prime }(x_k)$ 解释： 将 $f^{\prime }(x_k)$ 移到等式右侧。

假设 $f^{\prime \prime }(x_k)\ne 0$，我们可以两边同除以 $f^{\prime \prime }(x_k)$: $x_{k+1}-x_k=-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$ 解释： 两边除以 $f^{\prime \prime }(x_k)$。

最后，将 $x_k$ 移到等式右侧，得到牛顿法的迭代公式： $x_{k+1}=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$ 这就是单变量函数牛顿法的核心迭代公式。

3.2 多变量函数 $f(\mathbf{x})$ 的牛顿法推导

我们的目标是找到 $\mathbf{x}$ 使得 $\nabla f(\mathbf{x})=0$。

步骤 1：二阶泰勒展开 我们用一个二次曲面来近似函数 $f(\mathbf{x})$ 在当前点 ${\mathbf{x}}_k$ 附近的行为。这个二次曲面就是 $f(\mathbf{x})$ 在 ${\mathbf{x}}_k$ 处的二阶泰勒展开式： $f(\mathbf{x})\approx f({\mathbf{x}}_k)+\nabla f({\mathbf{x}}_k{)}^T(\mathbf{x}-{\mathbf{x}}_k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_k{)}^T\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)$ 令 $g(\mathbf{x})=f({\mathbf{x}}_k)+\nabla f({\mathbf{x}}_k{)}^T(\mathbf{x}-{\mathbf{x}}_k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_k{)}^T\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)$ 为这个近似函数。

步骤 2：找到近似函数的最小值 要找到 $g(\mathbf{x})$ 的最小值，我们对其求梯度，并令梯度等于零（因为最小值处的梯度为零向量）。

我们对 $g(\mathbf{x})$ 关于 $\mathbf{x}$ 求梯度： $\nabla g(\mathbf{x})={\nabla }_{\mathbf{x}}\left(f({\mathbf{x}}_k)+\nabla f({\mathbf{x}}_k{)}^T(\mathbf{x}-{\mathbf{x}}_k)+\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_k{)}^T\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)\right)$ 解释： 逐项求梯度。

• $f({\mathbf{x}}_k)$ 是常数项，求梯度为 $0$。

• $\nabla f({\mathbf{x}}_k{)}^T(\mathbf{x}-{\mathbf{x}}_k)$ 是线性项（例如，如果 $\nabla f({\mathbf{x}}_k)=[a,b{]}^T$ 且 $\mathbf{x}-{\mathbf{x}}_k=[u,v{]}^T$，则该项为 $au+bv$）。对于一个向量 ${\mathbf{c}}^T\mathbf{x}$，其梯度是 $\mathbf{c}$。所以这一项的梯度是 $\nabla f({\mathbf{x}}_k)$。

  • 复习：向量微积分中的梯度规则
    • ${\nabla }_{\mathbf{x}}({\mathbf{c}}^T\mathbf{x})=\mathbf{c}$
    • ${\nabla }_{\mathbf{x}}({\mathbf{x}}^T\mathbf{Ax})=(\mathbf{A}+{\mathbf{A}}^T)\mathbf{x}$
    • 如果 $\mathbf{A}$ 是对称矩阵，则 ${\nabla }_{\mathbf{x}}({\mathbf{x}}^T\mathbf{Ax})=2\mathbf{Ax}$

• $\frac{1}{2}(\mathbf{x}-{\mathbf{x}}_k{)}^T\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)$ 是二次型。令 $\mathbf{y}=\mathbf{x}-{\mathbf{x}}_k$。则该项为 $\frac{1}{2}{\mathbf{y}}^T\mathbf{H}({\mathbf{x}}_k)\mathbf{y}$。 由于 $\mathbf{H}({\mathbf{x}}_k)$ 是对称矩阵，根据上述梯度规则，其梯度是 $2\cdot \frac{1}{2}\mathbf{H}({\mathbf{x}}_k)\mathbf{y}=\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)$。

结合以上，我们得到 $g(\mathbf{x})$ 的梯度： $\nabla g(\mathbf{x})=\nabla f({\mathbf{x}}_k)+\mathbf{H}({\mathbf{x}}_k)(\mathbf{x}-{\mathbf{x}}_k)$

为了找到 $g(\mathbf{x})$ 的最小值点（我们称之为 ${\mathbf{x}}_{k+1}$），我们将这个梯度设为零向量： $\nabla f({\mathbf{x}}_k)+\mathbf{H}({\mathbf{x}}_k)({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)=0$

步骤 3：解出 ${\mathbf{x}}_{k+1}$ 现在我们从上面的方程中解出 ${\mathbf{x}}_{k+1}$。 $\mathbf{H}({\mathbf{x}}_k)({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)=-\nabla f({\mathbf{x}}_k)$ 解释： 将 $\nabla f({\mathbf{x}}_k)$ 移到等式右侧。

假设 Hessian 矩阵 $\mathbf{H}({\mathbf{x}}_k)$ 是可逆的，我们可以两边同乘以 $\mathbf{H}({\mathbf{x}}_k{)}^{-1}$: $\mathbf{H}({\mathbf{x}}_k{)}^{-1}\mathbf{H}({\mathbf{x}}_k)({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)=-\mathbf{H}({\mathbf{x}}_k{)}^{-1}\nabla f({\mathbf{x}}_k)$ 解释： 两边左乘 $\mathbf{H}({\mathbf{x}}_k{)}^{-1}$。

由于 $\mathbf{H}({\mathbf{x}}_k{)}^{-1}\mathbf{H}({\mathbf{x}}_k)=\mathbf{I}$（单位矩阵），我们得到： $\mathbf{I}({\mathbf{x}}_{k+1}-{\mathbf{x}}_k)=-\mathbf{H}({\mathbf{x}}_k{)}^{-1}\nabla f({\mathbf{x}}_k)$ ${\mathbf{x}}_{k+1}-{\mathbf{x}}_k=-\mathbf{H}({\mathbf{x}}_k{)}^{-1}\nabla f({\mathbf{x}}_k)$

最后，将 ${\mathbf{x}}_k$ 移到等式右侧，得到多变量牛顿法的迭代公式： ${\mathbf{x}}_{k+1}={\mathbf{x}}_k-\mathbf{H}({\mathbf{x}}_k{)}^{-1}\nabla f({\mathbf{x}}_k)$ 这就是多变量函数牛顿法的核心迭代公式。这里的 $\mathbf{H}({\mathbf{x}}_k{)}^{-1}\nabla f({\mathbf{x}}_k)$ 常常被称为牛顿方向 (Newton Direction)。

总结： 无论是单变量还是多变量，牛顿法的迭代公式都是： 新的点 = 当前点 - (二阶导数的逆) * (一阶导数)

---

4. 实例 (Concrete Numerical Example)

让我们通过一个简单的例子来演示牛顿法如何工作。

单变量例子： 假设我们想找到函数 $f(x)=x^2-4x+7$ 的最小值。 这个函数是一个抛物线，其最小值点可以通过求导并令导数等于 0 找到：$f^{\prime }(x)=2x-4=0⟹x=2$。 我们用牛顿法来验证。

首先，我们需要计算一阶导数和二阶导数：

• $f^{\prime }(x)=2x-4$
• $f^{\prime \prime }(x)=2$

现在，我们选择一个初始点 $x_0$。让我们选择 $x_0=0$。

第一次迭代 ($k=0$)：

• 计算 $f^{\prime }(x_0)=f^{\prime }(0)=2(0)-4=-4$
• 计算 $f^{\prime \prime }(x_0)=f^{\prime \prime }(0)=2$
• 应用牛顿法更新公式： $x_1=x_0-\frac{f^{\prime }(x_0)}{f^{\prime \prime }(x_0)}$ $x_1=0-\frac{-4}{2}$ $x_1=0-(-2)$ $x_1=2$ 看！我们仅用一次迭代就直接找到了最小值点 $x=2$。这是因为我们选择的函数是一个二次函数，它的二阶泰勒展开式就是它本身，所以近似是完全准确的。

多变量例子： 假设我们想找到函数 $f(x_1,x_2)=x_1^2+2x_2^2-4x_1-8x_2+10$ 的最小值。

首先，我们需要计算梯度向量 $\nabla f(\mathbf{x})$ 和 Hessian 矩阵 $\mathbf{H}(\mathbf{x})$。

• 计算梯度： $\frac{\partial f}{\partial x_1}=2x_1-4$ $\frac{\partial f}{\partial x_2}=4x_2-8$ 所以，梯度向量是： $\nabla f(\mathbf{x})=\left[\begin{array}{c}2x_1-4\\ 4x_2-8\end{array}\right]$

• 计算 Hessian 矩阵： $\frac{{\partial }^{2}f}{\partial x_1^2}=2$ $\frac{{\partial }^{2}f}{\partial x_1\partial x_2}=0$ $\frac{{\partial }^{2}f}{\partial x_2\partial x_1}=0$ $\frac{{\partial }^{2}f}{\partial x_2^2}=4$ 所以，Hessian 矩阵是： $\mathbf{H}(\mathbf{x})=\left[\begin{array}{cc}2& 0\\ 0& 4\end{array}\right]$ 注意到 Hessian 矩阵在这个例子中是一个常数矩阵，不依赖于 $\mathbf{x}$。

现在，我们选择一个初始点 ${\mathbf{x}}_0=\left[\begin{array}{c}0\\ 0\end{array}\right]$。

第一次迭代 ($k=0$)：

• 在 ${\mathbf{x}}_0=\left[\begin{array}{c}0\\ 0\end{array}\right]$ 处计算梯度： $\nabla f({\mathbf{x}}_0)=\left[\begin{array}{c}2(0)-4\\ 4(0)-8\end{array}\right]=\left[\begin{array}{c}-4\\ -8\end{array}\right]$

• 计算 Hessian 矩阵的逆： $\mathbf{H}({\mathbf{x}}_0{)}^{-1}={\left[\begin{array}{cc}2& 0\\ 0& 4\end{array}\right]}^{-1}$ 对于对角矩阵，逆矩阵是其对角元素取倒数： $\mathbf{H}({\mathbf{x}}_0{)}^{-1}=\left[\begin{array}{cc}1/2& 0\\ 0& 1/4\end{array}\right]=\left[\begin{array}{cc}0.5& 0\\ 0& 0.25\end{array}\right]$

• 应用牛顿法更新公式： ${\mathbf{x}}_1={\mathbf{x}}_0-\mathbf{H}({\mathbf{x}}_0{)}^{-1}\nabla f({\mathbf{x}}_0)$ ${\mathbf{x}}_1=\left[\begin{array}{c}0\\ 0\end{array}\right]-\left[\begin{array}{cc}0.5& 0\\ 0& 0.25\end{array}\right]\left[\begin{array}{c}-4\\ -8\end{array}\right]$ ${\mathbf{x}}_1=\left[\begin{array}{c}0\\ 0\end{array}\right]-\left[\begin{array}{c}(0.5)(-4)+(0)(-8)\\ (0)(-4)+(0.25)(-8)\end{array}\right]$ ${\mathbf{x}}_1=\left[\begin{array}{c}0\\ 0\end{array}\right]-\left[\begin{array}{c}-2\\ -2\end{array}\right]$ ${\mathbf{x}}_1=\left[\begin{array}{c}2\\ 2\end{array}\right]$ 再次，我们仅用一次迭代就直接找到了最小值点 $\mathbf{x}=\left[\begin{array}{c}2\\ 2\end{array}\right]$。这同样是因为我们选择的函数是一个二次函数，其二阶泰勒展开式就是它本身。 对于更复杂的非二次函数，牛顿法通常需要多次迭代才能收敛，但其收敛速度通常比梯度下降法快得多。

---

5. NumPy 实现 (Code Implementation)

我们将使用 Python 和 NumPy 来实现牛顿法。

import numpy as np
 
# --- 1. 单变量函数牛顿法示例 ---
 
def f_univariate(x):
    """单变量目标函数：f(x) = x^2 - 4x + 7"""
    return x**2 - 4*x + 7
 
def df_univariate(x):
    """单变量目标函数的一阶导数：f'(x) = 2x - 4"""
    return 2*x - 4
 
def ddf_univariate(x):
    """单变量目标函数的二阶导数：f''(x) = 2"""
    return 2
 
def newton_method_univariate(initial_x, tol=1e-6, max_iter=100):
    """
    单变量牛顿法实现
    Args:
        initial_x (float): 初始猜测值 x_0
        tol (float): 收敛容差
        max_iter (int): 最大迭代次数
    Returns:
        float: 找到的最小值点 x
    """
    x_k = initial_x
    print(f"--- 单变量牛顿法开始 (初始点 x_0 = {x_k}) ---")
    for i in range(max_iter):
        grad = df_univariate(x_k)   # 对应数学公式中的 f'(x_k)
        hess = ddf_univariate(x_k)  # 对应数学公式中的 f''(x_k)
 
        if abs(hess) < 1e-10: # 避免除以零，或者Hessian为零导致无法更新
            print(f"迭代 {i+1}: 二阶导数过小，可能无法收敛或遇到鞍点。当前 x = {x_k}")
            break
        
        # 牛顿法更新公式：x_{k+1} = x_k - f'(x_k) / f''(x_k)
        delta_x = grad / hess 
        x_new = x_k - delta_x
        
        print(f"迭代 {i+1}: x_k={x_k:.6f}, f'(x_k)={grad:.6f}, f''(x_k)={hess:.6f}, delta_x={delta_x:.6f}, x_{{k+1}}={x_new:.6f}")
        
        if abs(x_new - x_k) < tol:
            print(f"--- 单变量牛顿法收敛于 {x_new:.6f}，迭代次数：{i+1} ---")
            return x_new
        
        x_k = x_new
    
    print(f"--- 单变量牛顿法未收敛 (达到最大迭代次数)，最终结果：{x_k:.6f} ---")
    return x_k
 
# 运行单变量牛顿法
# result_univariate = newton_method_univariate(initial_x=0.0) 
# result_univariate = newton_method_univariate(initial_x=5.0)
 
# --- 2. 多变量函数牛顿法示例 ---
 
def f_multivariate(x_vec):
    """多变量目标函数：f(x1, x2) = x1^2 + 2x2^2 - 4x1 - 8x2 + 10"""
    x1, x2 = x_vec[0], x_vec[1]
    return x1**2 + 2*x2**2 - 4*x1 - 8*x2 + 10
 
def grad_multivariate(x_vec):
    """多变量目标函数的梯度向量：nabla f(x) = [2x1 - 4, 4x2 - 8]^T"""
    x1, x2 = x_vec[0], x_vec[1]
    return np.array([2*x1 - 4, 4*x2 - 8])
 
def hessian_multivariate(x_vec):
    """多变量目标函数的Hessian矩阵：H(x) = [[2, 0], [0, 4]]"""
    # 这个例子中Hessian是常数，不依赖x_vec
    return np.array([[2, 0], [0, 4]])
 
def newton_method_multivariate(initial_x_vec, tol=1e-6, max_iter=100):
    """
    多变量牛顿法实现
    Args:
        initial_x_vec (np.array): 初始猜测向量 x_0
        tol (float): 收敛容差
        max_iter (int): 最大迭代次数
    Returns:
        np.array: 找到的最小值点 x 向量
    """
    x_k = np.array(initial_x_vec, dtype=float)
    print(f"\n--- 多变量牛顿法开始 (初始点 x_0 = {x_k}) ---")
    for i in range(max_iter):
        grad = grad_multivariate(x_k)    # 对应数学公式中的 nabla f(x_k)
        hess = hessian_multivariate(x_k) # 对应数学公式中的 H(x_k)
 
        # 检查Hessian是否可逆 (正定或负定，用于最小化通常要求正定)
        # 实际应用中，通常会用线性方程组求解 H * delta_x = -grad，而不是直接求逆
        # 因为求逆的计算量大且数值稳定性差。这里为了清晰展示公式，我们直接求逆。
        try:
            hess_inv = np.linalg.inv(hess) # 对应数学公式中的 H(x_k)^-1
        except np.linalg.LinAlgError:
            print(f"迭代 {i+1}: Hessian矩阵不可逆。当前 x = {x_k}")
            break
        
        # 牛顿法更新公式：x_{k+1} = x_k - H(x_k)^-1 * nabla f(x_k)
        # 这里的 @ 运算符是矩阵乘法
        delta_x = hess_inv @ grad # 对应数学公式中的 H(x_k)^-1 * nabla f(x_k)
        x_new = x_k - delta_x
        
        print(f"迭代 {i+1}: x_k={x_k}, grad={grad}, hess={hess.tolist()}, delta_x={delta_x}, x_{{k+1}}={x_new}")
        
        if np.linalg.norm(x_new - x_k) < tol: # 使用L2范数判断收敛
            print(f"--- 多变量牛顿法收敛于 {x_new}，迭代次数：{i+1} ---")
            return x_new
        
        x_k = x_new
    
    print(f"--- 多变量牛顿法未收敛 (达到最大迭代次数)，最终结果：{x_k} ---")
    return x_k
 
# 运行多变量牛顿法
result_multivariate = newton_method_multivariate(initial_x_vec=[0.0, 0.0])

代码与数学公式的关联：

• df_univariate(x) 和 grad_multivariate(x_vec) 对应数学公式中的 $f^{\prime }(x_k)$ 和 $\nabla f({\mathbf{x}}_k)$。
• ddf_univariate(x) 和 hessian_multivariate(x_vec) 对应数学公式中的 $f^{\prime \prime }(x_k)$ 和 $\mathbf{H}({\mathbf{x}}_k)$。
• np.linalg.inv(hess) 对应数学公式中的 $\mathbf{H}({\mathbf{x}}_k{)}^{-1}$。
• x_new = x_k - delta_x 对应核心的迭代公式 $x_{k+1}=x_k-\frac{f^{\prime }(x_k)}{f^{\prime \prime }(x_k)}$ 和 ${\mathbf{x}}_{k+1}={\mathbf{x}}_k-\mathbf{H}({\mathbf{x}}_k{)}^{-1}\nabla f({\mathbf{x}}_k)$。
• @ 运算符是 NumPy 中进行矩阵乘法的操作，等价于数学公式中的矩阵与向量的乘积。

请注意，在实际应用中，计算Hessian矩阵的逆（np.linalg.inv(hess)）可能非常昂贵且不稳定，特别是对于高维问题。更常用的做法是解一个线性方程组 $\mathbf{H}({\mathbf{x}}_k){\mathbf{p}}_k=-\nabla f({\mathbf{x}}_k)$ 来找到更新方向 ${\mathbf{p}}_k$，然后令 ${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\mathbf{p}}_k$。这个过程可以通过 np.linalg.solve(hess, -grad) 来完成，它比直接求逆更有效率和数值稳定性。

6. 需要注意的点