最小二乘法(基本原理+代码解析）

前言：

在工程应用中，经常会遇到目标函数是若干个函数平方和的最优化问题。

这类问题称为最小二乘问题。

本章节讲述无约束规划中的最小二乘法。

---

线性最小二乘法：

当  是线性函数时，即  ，如果令 ，，

则

原问题可以表示为： 

求解该方程就变成了寻找一个向量 x ，使得  达到最小。

如果向量  能够使得  达到最小，即对与所有  都有：

则称  为  的最小二乘解。当方程组  有解时，其解自然是一个最小二乘解；若  无解，则只能寻求其最小二乘解，即能够使得 Ax 和 b 之间差值的范数达到最小的向量。

下面给出两个引理：

---

引理1：

矩阵  ,,当且仅当  （即方阵  非奇异）时，。

引理2：

能够最小化  的向量  具有唯一性。可以通过求解方程组  得到，即 。

---

例题分析：

例题：求线性最小二乘问题：  的最优解

其中 ，

解：可以直接通过公式  判断是否有解，

,

又因为

所以：方程组 无解。

---

非线性最小二乘法（LM算法）：

基本思想是把函数线性化，用线性最小二乘问题的解去逼近非线性最小二乘问题的解。

它通过同时利用高斯牛顿方法和梯度下降方法来解决非线性最小二乘问题。其核心思想是在每次迭

代中，根据当前参数估计计算目标函数的梯度和海森矩阵，并使用这些信息来更新模型参数。

算法步骤：

步骤1：选取初始点 ，给定初始值 ，放大因子 ，允许误差 ；

步骤2：求出 ；

步骤3：求出雅克比矩阵 ；

雅克比矩阵

步骤4：如果，停止，  是近似最优解；否则，转步骤5；

步骤5：构造 LM 方向，；

步骤6：计算及，如果 ，转步骤8；否则，转步骤7；

步骤7：放大系数，令  ，返回步骤5；

步骤8：缩小参数，令 ，返回步骤3.

(注： 不能取得太小，否则不能保证  是下降方向； 更不能取得太大，会影响收敛速度；根据经验， =0.01， =10)

算法代码：

'''
最小二乘法（LM算法）
2023-11-13
'''
import numpy as np

# 定义目标函数
def fun(x, t, y):
    return x[0] * np.exp(-x[1] * t) - y

# 定义计算雅可比矩阵的函数
def jac(x, t):
    J = np.empty((len(t), 2))
    J[:, 0] = np.exp(-x[1] * t)
    J[:, 1] = -x[0] * t * np.exp(-x[1] * t)
    return J

# 定义最小二乘法函数
def least_squares(fun, x0, jac, t, y):
    max_iter = 100
    tol = 1e-6
    lam = 0.01
    x = x0.copy()

    for i in range(max_iter):
        r = fun(x, t, y)
        J = jac(x, t)
        JtJ = np.dot(J.T, J)
        g = np.dot(J.T, r)

        # 计算LM步长
        dx = np.linalg.solve(JtJ + lam * np.diag(np.diag(JtJ)), -g)

        x_new = x + dx
        r_new = fun(x_new, t, y)

        # 更新参数和LM参数
        if np.linalg.norm(r_new) < np.linalg.norm(r):
            x = x_new
            lam *= 0.1
        else:
            lam *= 10

        # 判断收敛条件
        if np.linalg.norm(dx) < tol:
            break

    return x

# 生成样本数据
t = np.linspace(0, 4, 50)
actual_params = np.array([2.0, 0.5])
y = fun(actual_params, t, 0) + 0.1 * np.random.randn(len(t))

# 定义初始参数值
x0 = np.array([1.0, 1.0])

# 使用最小二乘法拟合数据
res = least_squares(fun, x0, jac, t, y)

# 打印拟合结果
print("拟合参数：", res)

该代码所运算的结果为：

拟合参数： [1.98729463 0.49536682]

（行文中若有纰漏，希望大家指正）