如何提升SciPy最小二乘优化精度？多侧距非线性方程组求解

问题背景 #

我尝试用SciPy的scipy.optimize.least_squares求解一个与多侧距（multilateration）相关的35个非线性方程组：

• 空间中有5个坐标未知的固定参考点、7个坐标未知的固定测量点，目标是在以其中一个参考点为原点的局部坐标系下求解7个测量点的坐标。
• 测量数据为参考点到测量点与参考点之间某虚拟点的距离（纳米级高精度），测量点到该虚拟点的距离称为「死区（deadzone）」，所有测量的死区长度相同但未知。
• 核心问题可概括为：在3D空间中对7组5个球体求交点，仅已知球体半径，且额外存在未知死区约束。
• 距离关系：参考点i坐标为$(x_i; y_i; z_i)$，测量点j坐标为$(x_j; y_j; z_j)$，$D_{ij}$为测量距离，$Dz$为死区，实际距离满足$\text{参考点到测量点距离} = D_{ij} + Dz$。

我期望输入初始值精度为$1e0$，输出坐标精度达到$1e-6$，通过最小化测量距离与计算距离的差值实现优化。

当前困境 #

程序可运行并返回近似合理的解，但精度仅能达到$1e-1$，远未达标。我通过在CAD模拟的真实测量点坐标中添加噪声（如$z_7 = 28.52317743 + 1$）进行测试，尝试过调整SciPy参数、集成mpmath、改用残差总和优化等手段，但效果有限。

现有代码 #

残差计算代码 #

# Calculate the residuals
residuals = []
for j, (xj, yj, zj) in enumerate(measured_points):
    
    residual = (
        np.sqrt(xj**2 + yj**2 + zj**2) - (dA[j] + Dzone) +
        np.sqrt((xj - xB)**2 + yj**2 + zj**2) - (dB[j] + Dzone) +
        np.sqrt((xj - xC)**2 + (yj - yC)**2 + zj**2) - (dC[j] + Dzone) +
        np.sqrt((xj - xD)**2 + (yj - yD)**2 + (zj - zD)**2) - (dD[j] + Dzone) +
        np.sqrt((xj - xE)**2 + (yj - yE)**2 + (zj - zE)**2) - (dE[j] + Dzone)
    )**2
    residuals.append(residual)

return np.array(residuals)

最小二乘调用代码 #

# Perform non-linear least squares regression
result = least_squares(
    objective_function, initial_params, jac='3-point', diff_step=1e-8,
    ftol=1e-15, xtol=1e-15, gtol=1e-15, x_scale=1, method='dogbox',
    loss='cauchy', tr_solver='lsmr', max_nfev=100000000, args=(distances,), verbose=2
)
optimal_params = result.x

已知相同方程和最小二乘回归在Maple中添加单变量噪声后精度可达$1e-12$，说明目标精度可实现。

疑问 #

1. 我是否正确使用了SciPy的最小二乘函数？
2. 在Python中有哪些方法可以提升求解精度？

一、当前代码的核心问题 #

1. 残差定义错误 #

least_squares要求输入的残差是每个独立约束的差值，而非差值的平方和。当前代码将每个测量点对应的5个约束差值求和后再平方，丢失了单个约束的误差信息，导致优化器无法正确感知每个方程的残差，严重影响收敛精度。

2. 损失函数与方法选择不当 #

• loss='cauchy'会降低大残差的权重，但高精度收敛需求下应优先使用默认的loss='linear'（标准最小二乘），避免损失函数引入的近似。
• method='dogbox'适合小规模问题，对于35维的优化变量，method='trf'（Trust Region Reflective）通常表现更稳定、收敛精度更高。

3. 数值微分的精度限制 #

使用jac='3-point'虽比'2-point'准确，但数值微分本身会引入误差。提供解析雅可比矩阵是提升优化精度和速度的关键。

二、精度提升的具体措施 #

1. 修正残差计算逻辑 #

重新构造残差数组，每个独立约束对应一个残差项（最终生成35个残差，对应35个方程组）：

def objective_function(params, distances):
    # 解析参数：参考点坐标、测量点坐标、死区Dz
    idx = 0
    xB, yB, zB = params[idx:idx+3]; idx +=3
    xC, yC, zC = params[idx:idx+3]; idx +=3
    xD, yD, zD = params[idx:idx+3]; idx +=3
    xE, yE, zE = params[idx:idx+3]; idx +=3
    measured_points = params[idx:idx+7*3].reshape(7,3); idx +=21
    Dzone = params[idx]
    
    dA, dB, dC, dD, dE = distances
    residuals = []
    for j in range(7):
        xj, yj, zj = measured_points[j]
        # 每个参考点对应一个残差项
        residuals.append(np.sqrt(xj**2 + yj**2 + zj**2) - (dA[j] + Dzone))
        residuals.append(np.sqrt((xj - xB)**2 + (yj - yB)**2 + (zj - zB)**2) - (dB[j] + Dzone))
        residuals.append(np.sqrt((xj - xC)**2 + (yj - yC)**2 + (zj - zC)**2) - (dC[j] + Dzone))
        residuals.append(np.sqrt((xj - xD)**2 + (yj - yD)**2 + (zj - zD)**2) - (dD[j] + Dzone))
        residuals.append(np.sqrt((xj - xE)**2 + (yj - yE)**2 + (zj - zE)**2) - (dE[j] + Dzone))
    return np.array(residuals)

2. 调整least_squares参数 #

result = least_squares(
    objective_function, initial_params, 
    jac='3-point', diff_step=1e-10,  # 调小微分步长提升数值精度
    ftol=1e-12, xtol=1e-12, gtol=1e-12,  # 设置严格的收敛阈值
    x_scale='jac',  # 自动根据雅可比矩阵平衡变量梯度
    method='trf',  # 更适合高维问题的信任域方法
    loss='linear',  # 标准最小二乘，无近似损失
    tr_solver='exact',  # 小型问题用精确求解器，精度更高
    max_nfev=10000,  # 合理设置迭代次数，无需过大
    args=(distances,), 
    verbose=2
)

3. 实现解析雅可比矩阵 #

手动推导每个残差对优化变量的偏导数，构造雅可比矩阵。例如，对于残差项$r = \sqrt{(x_j - x_i)^2 + (y_j - y_i)^2 + (z_j - z_i)^2} - (D_{ij} + Dz)$：

• 对$x_j$的偏导数：$\frac{x_j - x_i}{\text{计算距离}}$
• 对$x_i$的偏导数：$\frac{x_i - x_j}{\text{计算距离}}$
• 对$Dz$的偏导数：$-1$

将jac参数设置为解析雅可比函数，可完全消除数值微分的误差，大幅提升收敛精度。

4. 优化初始值与变量缩放 #

• 若变量量级差异较大（如坐标为$1e4$量级，死区为$1e0$量级），手动设置x_scale为各变量的量级（如x_scale=[1e4,1e4,1e4,...1]），帮助优化器平衡梯度。
• 初始值尽量接近真实解：可使用CAD模拟的真实值作为初始值，或通过线性近似先求解粗略解作为输入。

5. 高精度数值计算 #

若上述方法仍无法满足精度，可使用numpy.float128（需平台支持）或mpmath进行高精度计算，将所有变量和计算切换为更高精度的浮点数。

三、验证与调试 #

• 优化完成后，检查result.cost（残差平方和）是否足够小，result.optimality（最优性条件度量）是否达到设定阈值。
• 将优化结果代入原方程，计算每个约束的残差值，确认其量级是否达到$1e-6$以下。

内容的提问来源于stack exchange，提问作者Tom