SciPy数值优化

SciPy 数值优化：从最小二乘到非线性求解，遇到什么问题该用什么

写代码做数据处理，很多问题最终都能归结成一个核心需求：找到一组最优的参数。

• 拟合一条曲线 → 找参数使误差最小
• 控制一个系统 → 找参数使性能最好
• 解一个方程组 → 找一组值使所有等式成立

SciPy 的 optimize 模块就是干这个的。但问题是：算法太多了——leastsq、curve_fit、fmin、fmin_bfgs、fsolve……新手看了直接懵。

---

最小二乘拟合：数据点 + 模型 = 找参数

90% 的优化需求都能归到这一类：你有数据点，你有一个带参数的模型，你想让模型尽可能贴合数据。

标准写法（老派）：

from scipy.optimize import leastsq
import numpy as np

x = np.array([1, 2, 3, 4, 5])
y = 2.5 * x + 1.3 + np.random.randn(5) * 0.5

def residuals(p, x, y):
    k, b = p
    return y - (k * x + b)

p0 = [1.0, 0.0]
result = leastsq(residuals, p0, args=(x, y))
k, b = result[0]
print(f"k={k:.3f}, b={b:.3f}")

更友好的写法（推荐）：

from scipy.optimize import curve_fit

def model(x, k, b):
    return k * x + b

popt, pcov = curve_fit(model, x, y, p0=[1.0, 0.0])
k, b = popt

curve_fit 比 leastsq 好用得多——不用自己写残差函数，不用操心 args 传参方式，直接定义模型函数就行。新代码优先用 curve_fit。

拟合正弦波（复杂模型示例）：

def model(x, A, freq, phase, offset):
    return A * np.sin(2 * np.pi * freq * x + phase) + offset

popt, pcov = curve_fit(model, x, y, p0=[0.8, 4.5, 0.0, 0.0])

曲线拟合的核心原则：初始值要给对。 给错了，结果可能完全不收敛。如果不知道初始值怎么给，先画图看一眼数据的大致形态，再反推参数范围。

函数最小值：没有数据点，只有函数

有些问题没有”数据点”，只有一个目标函数——你找一组参数，让这个函数的值最小。

典型例子：Rosenbrock 函数——优化算法的”试金石”。

f(x,y)=(1−x)2+100(y−x2)2f(x,y)=(1−x)2+100(y−x2)2

它在 (1,1) 处有最小值 0，但有一条狭窄的”山谷”，算法很容易在里面打转。

不用导数（Nelder-Mead，简单问题用）：

from scipy.optimize import fmin

def rosenbrock(p):
    x, y = p
    return (1 - x)**2 + 100 * (y - x**2)**2

result = fmin(rosenbrock, [-2, 2])

用导数（BFGS，收敛快）：

from scipy.optimize import fmin_bfgs

def rosenbrock_prime(p):
    x, y = p
    return np.array([
        -2*(1-x) - 400*x*(y - x**2),
        200*(y - x**2)
    ])

result = fmin_bfgs(rosenbrock, [-2, 2], fprime=rosenbrock_prime)

关键判断：如果你的目标函数能求导，一定要给导数。 BFGS 比 Nelder-Mead 快得多，变量越多差距越大。

算法选择的决策逻辑：

问题是什么类型？
│
├── 有数据点 + 有模型 → curve_fit（最小二乘拟合）
│
├── 有目标函数（最小值问题）
│   ├── 能求导 → fmin_bfgs（变量 < 100）/ fmin_cg（变量 > 100）
│   └── 不能求导 → fmin（变量 < 20）/ fmin_powell（变量 > 20）
│
└── 要求解方程组 → fsolve

非线性方程组：找一组值让多个等式同时成立

你需要找一组 (x, y, z)，让一堆方程同时等于 0。

from scipy.optimize import fsolve

def equations(p):
    x, y, z = p
    return [
        x**2 + y**2 - 1,   # 单位圆
        x - z,             # x = z
        y - np.sin(z)      # y = sin(z)
    ]

solution = fsolve(equations, [0.5, 0.5, 0.5])
# 约等于 [0.786, 0.618, 0.786]

fsolve 的两个核心问题：

1. 初始值很关键。 给的不好，收敛到奇怪的结果或者干脆不收敛。

2. 方程组越大，越要提供雅可比矩阵。 雅可比矩阵就是”每个方程对每个变量的偏导数”。

def jacobian(p):
    x, y, z = p
    return [
        [2*x, 2*y, 0],
        [1, 0, -1],
        [0, 1, -np.cos(z)]
    ]

solution = fsolve(equations, [0.5, 0.5, 0.5], fprime=jacobian)

实际经验： 50 个未知数、每个方程只跟 6 个变量相关的稀疏系统，提供雅可比矩阵能让求解速度提高 4-5 倍。