最小二乘拟合的理解

本文最后更新于：2026年8月29日 晚上

最小二乘法是一种数学回归分析方法，通常用于数据拟合，在误差估计、不确定度、系统辨识及预测、预报等数据处理诸多学科领域得到广泛应用。起源于十八世纪的大航海探索时期，发展于天文领域和航海领域。法国科学家勒让德于 1805 年首先提出了最小二乘法，后来在 1822 年由高斯证明了最小二乘法的优越性，也就是高斯 - 马尔可夫定理。其中最小二乘其实是日语的翻译，二乘在日语中也就是平方的翻译，中文其实可以称为最小平方法。

参考：

最小二乘法的本质是什么？ - 知乎 (zhihu.com)

Least Squares Method: What It Means, How to Use It, With Examples (investopedia.com)

数据拟合——非线性最小二乘法 - 知乎 (zhihu.com)

非线性最小二乘 - 知乎 (zhihu.com)

非线性最小二乘 | LeijieZhang

牛顿法 高斯牛顿法 | Cheng Wei's Blog (scm_mos.gitlab.io)

非线性优化整理 -3.Levenberg-Marquardt 法 (LM 法)

老饼讲解 |【原理】LM 算法思想与介绍 (bbbdata.com)

基本定义

最小二乘法（又称最小平方法）通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合，其他一些优化问题也可通过最小化能量或最大化熵用最小二乘法来表达。

若存在一对观测量 ，其中 ，，且满足以下函数：

其中 为待定参数。

为了寻找函数 的最佳估计值，对于给定的 组（通常 ）的观测数据 ，其中 ，求解目标函数：

取最小值时的参数 。求解这类问题称为最小二乘问题，也被称为最小二乘拟合。当 为线性函数时，则该问题是线性最小二乘，否则称为非线性最小二乘。

线性最小二乘问题

当拟合函数 为线性函数时，则称求解上式目标函数的问题为线性最小二乘问题。

例如当给定 组观测值 其中 ，且为简化分析此处 维度为一，求拟合方程 中的最佳参数 和 使得目标函数 最小。

将多个观测值带入可获得多个等式方程：

求解参数为 即 ，求解参数为 2 个，则观测值至少要有 3 组，否则不存在优化问题，因为该参数已经可以唯一确定了。所以考虑当给定 3 组观测值时，则有：

写成向量形式：

则可以理解为 1 个三维向量由另外 2 个三维的基向量组成，两个基向量分别是 和 ，若在三维空间里（对应 3 组观测值），由于 和 可以取任意值，所以这两个基向量可以理解为形成了一个空间平面，而当 不在此平面内，则我们无法求得与它重合的向量，那么我们需要寻求一个最接近它的点（距离最近），而一个点与一个平面最近的点就是这个点与这个平面的垂线相交点，交点为 且有向量 与平面 垂直则有：

原理示意如下图所示：

pPRAj6f

当观测数据变多时，则基向量由原来三维变为高维向量，但其垂直关系不变。仍以三维空间为例，由上式进一步化简：

因为 为任意向量，所以：

最终可求解最佳参数 得：

其中：

$$ ^{3}, ^{3},

, =

. $$

都是可由观测值获得，由于 为非方阵，所以需要在式中求逆时求的是其广义逆，即 Moore-Penrose Matrix Inverse。其中当 的条件数比较大，矩阵呈现病态。

若拟合方程为 ，且给定 组观测数据，则：

最终通过该式可求得最佳拟合参数。

针对线性最小二乘法，建议保持与前面 高斯牛顿法（G-N） 和 列文伯格 - 马夸尔特法（L-M） 一致的章节结构：

  • 算法核心实现：实现通用线性最小二乘求解器；
  • 算法应用示例：针对具体线性模型进行参数拟合验证。

推荐标题：

（1）算法核心实现：线性最小二乘法

（2）算法应用示例：基于线性最小二乘法的一元线性回归拟合

具体内容如下。

算法核心实现：线性最小二乘法

线性最小二乘法用于求解线性模型参数估计问题，其目标是寻找最优参数：

其中： 为观测数据构成的设计矩阵， 为观测值向量， 为待估计模型参数。

对于一般线性模型：

可以写成矩阵形式：

其中：

根据最小二乘几何意义，最优解满足残差：

与矩阵空间正交：

进一步得到正规方程：

因此：

在实际计算过程中，不直接计算：

而采用矩阵分解方法求解，例如 QR 分解或奇异值分解（SVD），避免矩阵病态导致的数值误差。

Python 核心实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import numpy as np


def linear_least_squares(A, y):
    """
    线性最小二乘法求解器
    求解：
        min ||Aw-y||^2
    参数
    ----------
    A:
        设计矩阵
        shape=(m,n)
    y:
        观测向量
        shape=(m,)
    返回
    ----------
    w:
        最优参数
    """
    A = np.asarray(A, dtype=float)
    y = np.asarray(y, dtype=float)
    # 方法 1：
    # 正规方程
    #
    # (A^T A)w=A^T y
    ATA = A.T @ A
    ATy = A.T @ y
    try:
        w = np.linalg.solve(ATA, ATy)
    except np.linalg.LinAlgError:
        # 矩阵病态时采用伪逆
        w = np.linalg.pinv(A) @ y
    return w

上述代码实现了线性最小二乘求解流程：

  1. 根据观测数据构造设计矩阵：

  2. 构造正规方程：

  3. 求解参数：

当：

接近奇异时，通过 Moore-Penrose 广义逆：

保证算法仍然可以获得稳定解。

算法应用示例：一元线性回归拟合

为了验证线性最小二乘算法的有效性，采用一元线性模型：

其中： 表示直线斜率， 表示截距。

将模型改写为：

因此：

给定多组带噪声观测数据：

通过线性最小二乘法求解最佳拟合直线。

Python 应用示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
import numpy as np
import matplotlib.pyplot as plt

# =========================
# 1. 构造实验数据
# =========================

np.random.seed(0)

# 真实模型：
# y=3x+2

x = np.linspace(0, 10, 30)
y_true = 3 * x + 2

# 添加随机噪声

y_obs = y_true + np.random.normal(0, 2, size=x.shape)

# =========================
# 2. 构造设计矩阵 A
# =========================

# A=
#
# [x1 1]
# [x2 1]
# ...

A = np.column_stack((x, np.ones_like(x)))

# =========================
# 3. 调用线性最小二乘求解
# =========================

w = linear_least_squares(A, y_obs)
a, b = w
print("拟合参数：")
print("a=", a)
print("b=", b)

# =========================
# 4. 计算拟合结果
# =========================

x_fit = np.linspace(0, 10, 100)
y_fit = a * x_fit + b

# =========================
# 5. 绘制结果
# =========================

plt.scatter(x, y_obs, label="observation")
plt.plot(x_fit, y_fit, label="least squares")
plt.plot(x_fit, 3 * x_fit + 2, linestyle="--", label="true model")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Linear Least Squares Fitting")
plt.legend()
plt.show()

实验结果分析

程序首先根据观测数据构造设计矩阵：

然后通过正规方程：

计算模型参数。

最终得到：

即：

  • 直线斜率；
  • 直线截距。

该实验展示了线性最小二乘法在线性回归、标定参数估计、传感器数据拟合以及系统辨识中的基本应用。

非线性最小二乘问题

当拟合函数 为非线性函数时，则称求解上式目标函数的问题为非线性最小二乘问题。

例如当给定 组观测值 其中 ，求拟合方程 中的最佳参数 ，令残差为 同时可得到 个参数方程， 同样是非线性函数：

所以为求最佳拟合，则其残差平方和最小，即 的二范数最小，求解的最优化问题可转变为求解目标函数 最小值，同样 是非线性函数矩阵：

求解非线性优化的方法可以主要分为线性搜索法（Line Search）和信赖域法（Trust Region Method）。

其中：高斯牛顿法属于线性搜索法，列文伯格—马夸尔特法属于信赖域法。

线性搜索核心思想是先确定搜索方向，然后确定步长进行迭代。

高斯牛顿法（G-N 法）

根据上式求得 对 的雅可比矩阵：

目标函数的梯度为：

在任意参数组合 处将目标函数线性化：

将其带入到梯度式中 ，并令式子等于零：

令 ，这里的 并非是 Hessian 矩阵，但可以理解为牛顿法里的 Hessian 矩阵的近似：

至此可以求得迭代步长 ：

其中：

高斯牛顿法计算过程步骤：

  1. 给定初值 ，根据上式求解迭代步长 ；
  2. 更新迭代 ；
  3. 判断 ，若成立则收敛，完成迭代，取得最佳参数 ；
  4. 若不收敛，则更新的参数 返回 1，2 步骤继续迭代。

相比于牛顿法，高斯牛顿法不用计算 Hessian 矩阵，直接用 近似，所以节省了计算量。但是高斯牛顿法要求矩阵 是可逆且正定的，而实际计算中是半正定的，所以 会出现奇异或病态的情况，此时增量的稳定性就会变差，导致迭代发散。另一方面，增量较大时，近似替代函数式 (8) 就会产生较大的误差，也会导致迭代发散。这是高斯牛顿法的缺陷。

理解。针对你的章节结构，“算法核心实现”应给出一个通用 G-N 求解器代码，不依赖具体模型；“算法应用示例”应选择一个具体非线性最小二乘问题，并调用前面的求解器完成参数估计。

因此建议调整为：

算法核心实现：高斯牛顿法（G-N 法）

本部分实现一个通用的高斯牛顿优化算法，用于求解非线性最小二乘问题：

其中： - ：残差函数； - ：残差函数关于参数的雅可比矩阵。

根据高斯牛顿迭代公式：

求得参数增量：

然后更新：

实际计算过程中避免显式求逆矩阵，而采用线性方程求解方法，提高数值稳定性。

Python 核心实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
import numpy as np


def gauss_newton(residual_func, jacobian_func, initial_w, max_iter=100, tol=1e-8):
    """
    高斯牛顿算法
    参数
    ----------
    residual_func:
        残差函数 g(w)
    jacobian_func:
        雅可比矩阵 J(w)
    initial_w:
        参数初值
    max_iter:
        最大迭代次数
    tol:
        收敛阈值
    返回
    ----------
    w:
        最优参数
    history:
        迭代过程
    """

    w = np.array(initial_w, dtype=float)
    history = []
    for k in range(max_iter):
        # 计算残差
        g = residual_func(w)
        # 计算雅可比矩阵
        J = jacobian_func(w)
        # 构造高斯牛顿方程
        H = J.T @ J
        b = -J.T @ g
        try:
            # 求解：
            # (J^T J) Δw = -J^T g
            delta_w = np.linalg.solve(H, b)
        except np.linalg.LinAlgError:
            # 矩阵奇异时采用最小二乘求解
            delta_w = np.linalg.lstsq(H, b, rcond=None)[0]
        # 更新参数
        w_new = w + delta_w
        # 保存迭代信息
        cost = 0.5 * np.sum(g**2)
        history.append(
            {"iteration": k, "parameter": w.copy(), "cost": cost, "step": delta_w}
        )
        # 判断收敛
        if np.linalg.norm(delta_w) < tol:
            w = w_new
            break
        w = w_new
    return w, history

上述代码实现了完整的 G-N 算法流程：

  1. 输入当前参数；
  2. 计算残差向量；
  3. 计算雅可比矩阵；
  4. 构造：
  5. 求解参数增量；
  6. 更新参数；
  7. 判断收敛。

算法应用示例：指数曲线拟合

为了验证高斯牛顿法的有效性，选择指数模型进行参数估计。

假设观测数据满足：

其中：

  • ：幅值参数；
  • ：指数变化率；
  • ：偏置项。

给定观测数据：

定义残差：

对应雅可比矩阵：

利用前面的高斯牛顿优化器，对模型参数进行迭代求解。

Python 应用示例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
import numpy as np
import matplotlib.pyplot as plt

# ==========================
# 1. 生成实验数据
# ==========================

np.random.seed(0)
x = np.linspace(0, 3, 50)

# 真实参数

true_w = np.array([2.5, -1.2, 0.5])
y = true_w[0] * np.exp(true_w[1] * x) + true_w[2]

# 加入随机噪声
y_noise = y + np.random.normal(0, 0.05, size=x.shape)

# ==========================
# 2. 定义残差函数
# ==========================

def residual(w):
    w0, w1, w2 = w
    return w0 * np.exp(w1 * x) + w2 - y_noise

# ==========================
# 3. 定义雅可比矩阵
# ==========================

def jacobian(w):
    w0, w1, w2 = w
    e = np.exp(w1 * x)
    J = np.zeros((len(x), 3))
    # 对三个参数求偏导
    J[:, 0] = e
    J[:, 1] = w0 * x * e
    J[:, 2] = 1
    return J

# ==========================
# 4. 设置初值
# ==========================

w0 = np.array([1.0, -0.5, 0.0])

# ==========================
# 5. 调用高斯牛顿算法
# ==========================

w_opt, history = gauss_newton(residual, jacobian, w0, max_iter=100)
print("真实参数：")
print(true_w)
print("\n估计参数：")
print(w_opt)
print("\n迭代次数：", len(history))

# ==========================
# 6. 绘制拟合结果
# ==========================

x_fit = np.linspace(0, 3, 200)
y_fit = w_opt[0] * np.exp(w_opt[1] * x_fit) + w_opt[2]
plt.scatter(x, y_noise, label="data")
plt.plot(x_fit, y_fit, label="Gauss Newton")
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.title("Nonlinear Curve Fitting Using Gauss-Newton")
plt.show()

实验结果分析

运行程序后，高斯牛顿法会从初始参数：

开始迭代，通过不断计算：逐渐逼近真实模型参数。

最终得到：

并使残差： 达到最小。

该实例展示了高斯牛顿法在非线性曲线拟合、参数估计以及最小二乘优化问题中的应用过程。同时也说明，当模型初值较合理且雅可比矩阵条件良好时，高斯牛顿法具有较快的收敛速度；而当矩阵 接近奇异时，则需要结合阻尼策略（如 Levenberg-Marquardt 方法）提高算法稳定性。

Trust Region

信赖域法核心思想则是先确定最大距离，再确定方向。具体为每次迭代都会给出一个半径为 的信赖域使得模型在该信赖域内足够精确。在该信赖域内求解问题得到一个步长 h，接着通过某一评价标准判断该步长的好坏来决定是否要接受该步长。如果接受则更新当前位置，否则则当前位置不变。而新信赖域的半径大小则依照判断预测的好坏，来增大或减小。

列文伯格 - 马夸尔特法（L-M 法）

Line Search 依赖线性化近似有较高的拟合度，但是有时候线性近似效果较差，导致迭代不稳定；Trust Region 就是解决了这种问题。高斯牛顿法中采用的近似二阶泰勒展开只在该点附近有较好的近似结果，对 添加一个信赖域区域，就变为 Trust Region 方法。

LM 法迭代公式：

当 较小时，接近于高斯牛顿法：

当 较大时，接近于最速下降法：

对于信赖域 的定义，通过近似模型与实际函数之间的差距来更新这个范围；如果这个差距小，则增大信赖域，否则减小信赖域。通过如下式子定义差距 ：

pPWCXJx

来判断泰勒近似是否够好， 的分子是实际模型下降的值，分母是近似模型下降的值。如果 接近于 1，则认为近似是好的。如果 太小，说明实际减小的值远少于近似减少的值，这认为近似比较差。反之，如果 比较大，则说明实际下降的比预计的更大，我们可以放大近似范围。

列文伯格 - 马夸尔特法计算过程步骤：

  1. 给定初值 ，以及初始化信赖域的半径 ，根据该式求解迭代步长 ；
  2. 根据该式求解近似模型的差距 ：
    • 若 ，则扩大信赖域半径
    • 若 ，则缩小信赖域半径
    • 若 ，信赖域半径不变，进行更新迭代 ；
  3. 判断 ，若成立则收敛，完成迭代，取得最佳参数 ；
  4. 若不收敛，则更新的参数 返回 1，2 步骤继续迭代。
  • LM 法是一种非常高效的迭代算法，能适用于大多数的例子
  • 该算法的许多变量决定了迭代的效率和结果的收敛性，而这些参数的理想值往往依目标函数不同而不同
  • 因为存在矩阵求逆这一运算，所以对于大型数据求解需要用一些特殊技巧

算法核心实现：L-M 非线性最小二乘

下面给出一个通用的 L-M 非线性最小二乘核心例程。目标是求解：

其中 residual_func(w) 返回残差向量 ，jacobian_func(w) 返回雅可比矩阵 。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
from __future__ import annotations

from typing import Callable, Optional
import numpy as np
from numpy.typing import NDArray

Array = NDArray[np.float64]

def numerical_jacobian(
    residual_func: Callable[[Array], Array],
    w: Array,
    step: float = 1e-6,
) -> Array:
    """
    使用中心差分计算残差函数 g(w) 的雅可比矩阵。
    J[i, j] = ∂g_i / ∂w_j
    """
    w = np.asarray(w, dtype=float)
    residual = np.asarray(residual_func(w), dtype=float).reshape(-1)
    jacobian = np.empty((residual.size, w.size), dtype=float)
    for j in range(w.size):
        # 根据参数尺度自适应设置差分步长
        h = step * max(1.0, abs(w[j]))
        w_forward = w.copy()
        w_backward = w.copy()
        w_forward[j] += h
        w_backward[j] -= h
        residual_forward = np.asarray(residual_func(w_forward), dtype=float).reshape(-1)
        residual_backward = np.asarray(residual_func(w_backward), dtype=float).reshape(
            -1
        )
        jacobian[:, j] = (residual_forward - residual_backward) / (2.0 * h)
    return jacobian

def levenberg_marquardt(
    residual_func: Callable[[Array], Array],
    initial_w: Array,
    jacobian_func: Optional[Callable[Array](Array) = None,
    initial_mu: float = 1e-3,
    max_iterations: int = 100,
    residual_tolerance: float = 1e-10,
    gradient_tolerance: float = 1e-8,
    step_tolerance: float = 1e-10,
    rho_expand: float = 0.75,
    rho_shrink: float = 0.25,
    verbose: bool = False,
) -> tuple[Array, dict]:
    """
    使用 Levenberg-Marquardt 法求解非线性最小二乘问题：
        min 0.5 * ||g(w)||^2
    参数
    ----
    residual_func:
        残差函数 g(w)，返回形状为 (m,) 的残差向量。
    initial_w:
        参数初始值，形状为 (n,)。
    jacobian_func:
        雅可比函数 J(w)，返回形状为 (m, n) 的矩阵。
        未提供时使用中心差分计算。
    initial_mu:
        初始阻尼参数 μ。
    返回
    ----
    w:
        优化后的参数。
    info:
        迭代信息。
    """
    w = np.asarray(initial_w, dtype=float).reshape(-1).copy()
    mu = float(initial_mu)

    if mu <= 0:
        raise ValueError("initial_mu 必须大于 0。")
    accepted_steps = 0
    converged = False
    stop_reason = "达到最大迭代次数"
    for iteration in range(max_iterations):
        residual = np.asarray(residual_func(w), dtype=float).reshape(-1)
        if not np.all(np.isfinite(residual)):
            raise ValueError("残差函数返回了 NaN 或无穷大。")
        cost = 0.5 * residual @ residual
        if jacobian_func is None:
            jacobian = numerical_jacobian(residual_func, w)
        else:
            jacobian = np.asarray(jacobian_func(w), dtype=float)

        if jacobian.shape != (residual.size, w.size):
            raise ValueError(
                f"雅可比矩阵形状应为 {(residual.size, w.size)}，"
                f"实际为 {jacobian.shape}。"
            )
        # 一阶梯度：∇F(w) = J^T g
        gradient = jacobian.T @ residual

        if np.linalg.norm(gradient, ord=np.inf) < gradient_tolerance:
            converged = True
            stop_reason = "梯度满足收敛条件"
            break
        # LM 方程：
        # (J^T J + μI) Δw = -J^T g
        normal_matrix = jacobian.T @ jacobian
        lm_matrix = normal_matrix + mu * np.eye(w.size)
        try:
            # 不显式计算矩阵逆，数值稳定性更好
            delta_w = np.linalg.solve(lm_matrix, -gradient)
        except np.linalg.LinAlgError:
            # 矩阵病态时退化为最小二乘求解
            delta_w = np.linalg.lstsq(lm_matrix, -gradient, rcond=None)[0]
        if np.linalg.norm(delta_w) <= (
            step_tolerance * (np.linalg.norm(w) + step_tolerance)
        ):
            converged = True
            stop_reason = "参数步长满足收敛条件"
            break
        candidate_w = w + delta_w
        candidate_residual = np.asarray(
            residual_func(candidate_w), dtype=float
        ).reshape(-1)
        candidate_cost = 0.5 * (candidate_residual @ candidate_residual)
        # 实际下降量：
        # F(w) - F(w + Δw)
        actual_reduction = cost - candidate_cost
        # 线性化模型的残差：
        # g(w) + J(w)Δw
        model_residual = residual + jacobian @ delta_w
        model_cost = 0.5 * (model_residual @ model_residual)
        # 预测下降量：
        # F(w) - M(Δw)
        predicted_reduction = cost - model_cost
        if predicted_reduction > 0:
            rho = actual_reduction / predicted_reduction
        else:
            rho = -np.inf
        # 仅接受能够实际降低目标函数的迭代步
        accepted = rho > 0 and actual_reduction > 0
        if accepted:
            old_residual = residual
            w = candidate_w
            accepted_steps += 1
            residual_change = np.linalg.norm(candidate_residual - old_residual)
            # rho 较大：模型预测准确，可减小阻尼，
            # 使算法更接近高斯 - 牛顿法，相当于扩大信赖域
            if rho > rho_expand:
                mu = max(mu * 0.5, 1e-15)
            # rho 较小：虽然接受，但模型预测较差，增大阻尼
            elif rho < rho_shrink:
                mu = min(mu * 2.0, 1e15)

            if candidate_cost < residual_tolerance:
                converged = True
                stop_reason = "残差目标值满足收敛条件"
            elif residual_change < residual_tolerance:
                converged = True
                stop_reason = "相邻两次残差变化满足收敛条件"
        else:
            # 拒绝该步，并增大阻尼，使下一步更保守
            mu = min(mu * 2.0, 1e15)
        if verbose:
            print(
                f"iteration={iteration:3d}, "
                f"cost={cost:.6e}, "
                f"mu={mu:.3e}, "
                f"rho={rho:.3e}, "
                f"step_norm={np.linalg.norm(delta_w):.3e}, "
                f"accepted={accepted}"
            )
        if converged:
            break
    final_residual = np.asarray(residual_func(w), dtype=float).reshape(-1)
    info = {
        "converged": converged,
        "stop_reason": stop_reason,
        "iterations": iteration + 1,
        "accepted_steps": accepted_steps,
        "final_cost": 0.5 * final_residual @ final_residual,
        "final_residual_norm": np.linalg.norm(final_residual),
        "final_mu": mu,
    }
    return w, info

算法应用示例：指数曲线拟合

拟合模型：

残差定义为：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
import numpy as np

x_data = np.linspace(0.0, 2.0, 30)

true_parameters = np.array([2.5, -1.3, 0.4])
rng = np.random.default_rng(42)

y_data = (
    true_parameters[0] * np.exp(true_parameters[1] * x_data)
    + true_parameters[2]
    + rng.normal(0.0, 0.03, size=x_data.size)
)

def residual_func(w: np.ndarray) -> np.ndarray:
    a, b, c = w
    predicted_y = a * np.exp(b * x_data) + c
    return predicted_y - y_data


def jacobian_func(w: np.ndarray) -> np.ndarray:
    a, b, _ = w
    exp_bx = np.exp(b * x_data)

    # 每一行对应一个残差，每一列对应一个待求参数
    return np.column_stack(
        (
            exp_bx,  # ∂g/∂a
            a * x_data * exp_bx,  # ∂g/∂b
            np.ones_like(x_data),  # ∂g/∂c
        )
    )
    
initial_w = np.array([1.0, -0.5, 0.0])

optimal_w, info = levenberg_marquardt(
    residual_func=residual_func,
    initial_w=initial_w,
    jacobian_func=jacobian_func,
    initial_mu=1e-2,
    max_iterations=200,
    verbose=True,
)
print("\n最优参数：", optimal_w)
print("真实参数：", true_parameters)
print("迭代信息：", info)

需要注意， 通常称为阻尼参数，它与信赖域半径的变化方向相反： - 减小 ：更接近高斯—牛顿法，步长通常变大，相当于扩大信赖域。 - 增大 ：更接近梯度下降法，步长通常变小，相当于缩小信赖域。

因此，当 较大、近似模型可靠时，代码中采用 mu *= 0.5；当 较小或迭代步被拒绝时，采用 mu *= 2.0。这与直接把 定义为“信赖域半径”的表述有所不同。

参考 牛顿法 高斯牛顿法 | Cheng Wei's Blog (scm_mos.gitlab.io) 非线性最小二乘 | LeijieZhang (leijiezhang001.github.io)