哈喽，大家好！~
牛顿法，先来一句最通俗的解释。
假设我们在一座山谷里找最低点（最小化损失），常见的方法是顺着坡度往下走：这就是梯度下降。
但有时候你想更聪明一点，不仅看坡度（一阶信息），还想看看地面是平坦还是陡峭（二阶信息，曲率）。
牛顿法就是：在当前位置把函数用一个“二阶抛物面”（二阶泰勒展开的二次近似）代替，然后直接跳到这个抛物面最低点上。
因为二次近似把局部信息用得更充分，收敛通常能快得多。
更正式一点，给定目标函数 （我们想最小化），在点 处做二阶泰勒展开：
其中 是梯度， 是 Hessian（海森矩阵，二阶导构成的矩阵）。把上式对 最小化，导数为零，得到牛顿更新：
直观上，牛顿步把梯度按曲率缩放：在曲率大的方向（函数陡峭）步子会被缩小，在曲率小的方向（比较平坦）步子相对增大，因此往往能在较少步数内收敛到极小点。
机器学习里我们常见的 Newton-Raphson、IRLS 等都属于牛顿法家族。
优点：
-
局部二阶收敛（接近极小点时收敛速度很快）。 -
步长自适应（不用像梯度下降那样反复调学习率）。
缺点 / 陷阱：
-
需要计算并反转 Hessian，计算量一般是 O(d^3)，内存 O(d^2)，对高维模型（如深度网络）不实用。 -
Hessian 可能奇异或不正定（尤其在非凸问题），需要加阻尼（Levenberg–Marquardt）或用牛顿CG 等方法解决。 -
对初始点敏感（尤其非凸问题可能跳到鞍点或局部最大）。
我们在机器学习实践中，牛顿法常用于中小规模问题（如逻辑回归、二次曲线拟合等），或者作为二阶信息近似（L-BFGS，Newton-CG，K-FAC 等）。
实战案例
下面我们给一个端到端的案例：用 PyTorch 在一个二维合成数据集上做二分类，比较牛顿法（标准牛顿 + 阻尼）和梯度下降。
后面我们会使用图形化，说明牛顿法的行为：决策边界演化、参数空间中的损失等高线与优化轨迹、损失随迭代变化、Hessian 特征值（条件数）随迭代变化等。
数据集
两类 Gaussian blob，带一点重叠，让问题不是太容易：
import torch
import torch.nn.functional as F
import numpy as np
import matplotlib.pyplot as plt
torch.set_default_dtype(torch.float64)
def make_blobs(n0=200, n1=200, seed=0):
np.random.seed(seed)
# 类别0中心 (-1, -1), 类别1中心 (1, 1)，但我们不加 bias，以保持通过原点的决策面可行
x0 = np.random.randn(n0, 2) * 0.6 + np.array([-1.0, -1.0])
x1 = np.random.randn(n1, 2) * 0.6 + np.array([1.0, 1.0])
X = np.vstack([x0, x1])
y = np.hstack([np.zeros(n0), np.ones(n1)])
return torch.from_numpy(X), torch.from_numpy(y).unsqueeze(1)
X, y = make_blobs()
N, D = X.shape
print("数据样本数:", N, "维度:", D)
我们生成二维数据以便可视化。模型将是简单的线性逻辑回归： ，这里不放 bias（为了后面在参数平面画损失更直观），因此参数维度等于 2。
定义 Logistic 损失、梯度与 Hessian
对数似然的负对数似然：
# Logistic negative log-likelihood (sum over dataset)
def nll_and_stats(w, X, y):
# w: (D,1), X: (N,D), y: (N,1)
z = X @ w # (N,1)
p = torch.sigmoid(z)
eps = 1e-12
# negative log-likelihood (sum)
nll = - (y * torch.log(p+eps) + (1-y) * torch.log(1-p+eps)).sum()
# gradient: X^T (p - y)
grad = X.t() @ (p - y) # (D,1)
# Hessian: X^T diag(p*(1-p)) X
S = (p * (1-p)).squeeze() # (N,)
# construct weighted X
XS = X * S.unsqueeze(1) # (N,D)
H = X.t() @ XS # (D,D)
return nll, grad, H, p
对于 Logistic 损失，我们可以给出闭式的 Hessian： ，这是对称正定（在非退化情况下），非常方便用于牛顿步。
优化器
牛顿法（带阻尼）与标准梯度下降。我们记录每次迭代的参数轨迹、损失、梯度范数、Hessian 的最大/最小特征值等用于分析与可视化。
def run_optim(X, y, method='newton', max_iter=30, tol=1e-8, lr=0.1, damp=1e-3):
# 初始化参数
w = torch.zeros((D,1), dtype=X.dtype)
traj = []
losses = []
grads = []
eigs_hist = []
for it in range(max_iter):
nll, grad, H, p = nll_and_stats(w, X, y)
loss = nll.item()
gnorm = torch.norm(grad).item()
# compute eigenvalues for diagnostics (numpy)
H_np = H.detach().numpy()
eigs = np.linalg.eigvalsh(H_np)
eigs_hist.append(eigs)
traj.append(w.detach().numpy().flatten().copy())
losses.append(loss)
grads.append(gnorm)
# 优化步
if method == 'gd':
w = w - lr * grad
elif method == 'newton':
# 阻尼（Levenberg-like）：H + damp * I 保证正定且数值稳定
H_damped = H + damp * torch.eye(D, dtype=X.dtype)
# solve H * delta = grad => delta = H^{-1} grad
# Newton step: w_new = w - delta
# 使用线性求解而不是矩阵逆
delta, _ = torch.solve(grad, H_damped) # torch.solve(b, A) solves Ax=b, returns (x, LU)
# 为了稳健起见，可以做步长缩放（但这里默认全步）
w = w - delta
else:
raise ValueError("Unknown method")
# 收敛判断
if gnorm < tol:
break
return {
'w': w,
'traj': np.array(traj),
'losses': np.array(losses),
'grads': np.array(grads),
'eigs': np.array(eigs_hist),
}
res_newton = run_optim(X, y, method='newton', max_iter=30, lr=0.5, damp=1e-3)
res_gd = run_optim(X, y, method='gd', max_iter=200, lr=0.05, damp=0.0)
我们在牛顿法里加了一个小的阻尼项 来避免 Hessian 奇异或条件数过差导致数值不稳定。这是 Levenberg–Marquardt 风格的一种简单处理。
直接用 torch.solve 求解线性方程比计算显式逆更稳定。
对比梯度下降（学习率 lr 需要调参），我们能看到收敛速度差别。
可视化
数据及决策边界演化、参数空间损失等高线及轨迹、损失随迭代、梯度范数随迭代、Hessian 特征值随迭代~
# 辅助：在输入空间绘制决策边界
def plot_decision_boundaries(X_np, y_np, trajs):
xs = X_np[:,0]; ys = X_np[:,1]
plt.figure(figsize=(6,6))
plt.scatter(xs[y_np.ravel()==0], ys[y_np.ravel()==0], s=40, c='orange', edgecolor='k', label='class 0', alpha=0.9)
plt.scatter(xs[y_np.ravel()==1], ys[y_np.ravel()==1], s=40, c='deepskyblue', edgecolor='k', label='class 1', alpha=0.9)
# 绘制决策边界对应多条参数点
its = [0, 1, 2, 4, 9, -1] # 迭代索引（-1 表示最后）
colors = ['purple','red','green','yellow','magenta','black']
xlim = (xs.min()-0.5, xs.max()+0.5)
ylim = (ys.min()-0.5, ys.max()+0.5)
xx = np.linspace(*xlim, 200)
yy = np.linspace(*ylim, 200)
Xg, Yg = np.meshgrid(xx, yy)
grid = np.vstack([Xg.ravel(), Yg.ravel()]).T
for idx, c in zip(its, colors):
if idx == -1:
w = trajs[-1]
label = 'final'
else:
w = trajs[idx]
label = f'it={idx}'
# decision line w0*x + w1*y = 0 => y = -(w0/w1) x
w0, w1 = w
if abs(w1) < 1e-6:
continue
y_line = -(w0/w1) * xx
plt.plot(xx, y_line, color=c, linewidth=2.0, label=f'boundary {label}')
plt.xlim(xlim); plt.ylim(ylim)
plt.legend()
plt.title('Data and decision boundary evolution')
plt.show()
plot_decision_boundaries(X.numpy(), y.numpy(), res_newton['traj'])
决策边界演化，在数据点上叠加了若干条模型在不同迭代时的决策边界（线）。
可以直观看到牛顿法如何快速把边界移动到合适位置。若用梯度下降绘制对应轨迹，会发现步子更小、需要更多次迭代才接近最终边界。
def plot_loss_contour(X, y, res_newton, res_gd):
# 网格
w0_vals = np.linspace(-4, 4, 200)
w1_vals = np.linspace(-4, 4, 200)
W0, W1 = np.meshgrid(w0_vals, w1_vals)
Z = np.zeros_like(W0)
X_np = X.numpy()
y_np = y.numpy()
for i in range(W0.shape[0]):
for j in range(W0.shape[1]):
w = np.array([W0[i,j], W1[i,j]]).reshape(2,1)
# compute NLL
z = X_np @ w
p = 1/(1+np.exp(-z))
nll = - (y_np * np.log(p+1e-12) + (1-y_np) * np.log(1-p+1e-12)).sum()
Z[i,j] = nll
plt.figure(figsize=(8,6))
cs = plt.contourf(W0, W1, Z, levels=50, cmap='RdYlBu_r')
plt.colorbar(cs)
# overlay trajectories
traj_n = res_newton['traj']
traj_g = res_gd['traj']
plt.plot(traj_n[:,0], traj_n[:,1], '-o', color='magenta', label='Newton', markersize=6)
plt.plot(traj_g[:,0], traj_g[:,1], '-o', color='lime', label='GradientDescent', markersize=4)
plt.scatter(traj_n[0,0], traj_n[0,1], color='black', marker='x', s=80, label='start')
plt.title('Loss surface (parameter space) and optimization trajectories')
plt.xlabel('w0'); plt.ylabel('w1')
plt.legend()
plt.show()
plot_loss_contour(X, y, res_newton, res_gd)
参数空间损失等高线 + 轨迹：显示在参数空间上损失函数的形状（等高线），牛顿法轨迹通常在很少几个点就跳到谷底附近（轨迹更粗、更短），而梯度下降会沿着等高线缓慢迂回（因为它只用一阶信息）。
plt.figure(figsize=(7,4))
plt.plot(res_newton['losses'], '-o', color='magenta', label='Newton', linewidth=2)
plt.plot(res_gd['losses'], '-s', color='lime', label='GD', linewidth=1.5)
plt.yscale('log')
plt.xlabel('iteration'); plt.ylabel('negative log-likelihood (log scale)')
plt.title('Loss vs iteration (Newton vs GD)')
plt.legend()
plt.show()
损失 vs 迭代，通常我们会在对数 y 轴上看到牛顿法损失在几步内就下降到非常低的值，而梯度下降需要很多步。注意如果使用不合适的阻尼或在非凸问题中，牛顿也可能振荡或失败。
# 梯度范数
plt.figure(figsize=(7,4))
plt.plot(res_newton['grads'], '-o', color='magenta', label='Newton grad norm')
plt.plot(res_gd['grads'], '--', color='lime', label='GD grad norm')
plt.yscale('log')
plt.xlabel('iteration'); plt.ylabel('||grad|| (log)')
plt.title('Gradient norm vs iteration')
plt.legend()
plt.show()
# Hessian eigenvalues：取每次迭代的最小与最大特征值
def plot_hessian_eigs(eigs_array):
# eigs_array: (iters, D)
max_eig = eigs_array.max(axis=1)
min_eig = eigs_array.min(axis=1)
cond = max_eig / (min_eig + 1e-12)
iters = np.arange(len(max_eig))
plt.figure(figsize=(8,4))
plt.plot(iters, max_eig, '-o', color='red', label='max eig')
plt.plot(iters, min_eig, '-s', color='blue', label='min eig')
plt.plot(iters, cond, '-^', color='green', label='condition number')
plt.yscale('log')
plt.xlabel('iteration'); plt.title('Hessian eigenvalues and condition number')
plt.legend()
plt.show()
plot_hessian_eigs(res_newton['eigs'])
梯度范数、Hessian 特征值：
梯度范数曲线说明收敛度：牛顿法的梯度范数通常更迅速地下降到接近 0。
Hessian 特征值显示参数点处曲率信息：最大特征值代表最陡方向的曲率，最小特征值代表最平坦方向。条件数（max/min）高说明问题沿某些方向非常扁（造成优化困难）。牛顿法等二阶方法会“自动”缩放不同方向的步长，因而能更好地处理条件数差异，但也更依赖于 Hessian 的精确性和稳定性（这就是我们加阻尼的原因）。
最后
当维度很小（d = 几十、几百），且需要快速收敛时，牛顿法或准牛顿（L-BFGS）是非常有效的。对 Logistic 回归等凸问题，Newton/IRLS 能在很少迭代内收敛。
对高维模型（深度网络），直接计算 Hessian 不现实。常用的方法是使用近似二阶信息（如 K-FAC、Hessian-vector product + CG 解线性系统），或使用 L-BFGS 等有限内存方法。
数值稳定很重要：Hessian 可能接近奇异、或者在非凸问题中出现负特征值。通用做法是对 Hessian 加阻尼项（H + λI）或做共轭梯度法（只用 Hessian-vector product）。
另外，牛顿步虽然收敛快，但每步成本高；在工程中常把牛顿作为“精炼器”：先用一阶方法快速找到大致区域，再用牛顿精细收敛。