Python中多Lorentzian曲线拟合的优化方法咨询

1. 简化拟合函数，用向量化替代冗余手写项 #

手动写10个Lorentzian项不仅代码冗余，还无法利用numpy的向量化加速。改成批量参数处理的向量化函数，既简洁又能提升计算效率：

原来的冗余写法（示例）：

def lorentzian_model(x, a1, xc1, w1, a2, xc2, w2, ..., a10, xc10, w10, bg):
    return (a1/(1 + ((x - xc1)/w1)**2) + 
            a2/(1 + ((x - xc2)/w2)**2) + 
            ... + 
            a10/(1 + ((x - xc10)/w10)**2) + bg)

优化后的向量化版本：

import numpy as np

def lorentzian_model(x, *params):
    params_arr = np.array(params)
    n_peaks = 10
    # 拆分峰参数（每组3个：振幅、中心、宽度）和背景
    peak_params = params_arr[:-1].reshape(n_peaks, 3)
    bg = params_arr[-1]

    # 提取所有峰的参数数组
    amplitudes = peak_params[:, 0]
    centers = peak_params[:, 1]
    widths = peak_params[:, 2]

    # 利用numpy广播批量计算所有峰的贡献，再求和
    peaks = amplitudes[:, np.newaxis] / (1 + ((x - centers[:, np.newaxis])/widths[:, np.newaxis])**2)
    return peaks.sum(axis=0) + bg

2. 优化初始参数，减少迭代次数 #

curve_fit的速度极大依赖初始参数的合理性，初始值偏离最优解会导致迭代次数暴增。可以通过峰值检测自动生成初始参数：

from scipy.signal import find_peaks

# x为波长数组，y为光谱强度数组
peaks, properties = find_peaks(y, height=np.percentile(y, 90), distance=20)  # 按需调整检测参数

# 生成初始参数列表
initial_params = []
for peak_idx in peaks[:10]:  # 确保取到10个峰
    initial_params.append(properties['peak_heights'][peak_idx])  # 振幅初始值
    initial_params.append(x[peak_idx])  # 中心初始值
    initial_params.append(2)  # 宽度初始值，根据数据分辨率调整
initial_params.append(np.min(y))  # 背景初始值

# 拟合时传入初始参数
popt, pcov = curve_fit(lorentzian_model, x, y, p0=initial_params)

3. 限制参数范围，缩小搜索空间 #

结合物理意义给参数设置上下界，避免无效的参数搜索，加速收敛：

# 定义每个参数的上下界
lower_bounds = []
upper_bounds = []
for _ in range(10):
    lower_bounds.extend([0, x.min(), 0.1])  # 振幅≥0，中心在光谱范围内，宽度≥0.1
upper_bounds.extend([y.max()*1.2, x.max(), 10])  # 振幅不超过峰值1.2倍，宽度上限按需调整
lower_bounds.append(y.min() - 10)  # 背景下界
upper_bounds.append(y.min() + 10)  # 背景上界

# 拟合时传入边界
popt, pcov = curve_fit(lorentzian_model, x, y, p0=initial_params, bounds=(lower_bounds, upper_bounds))

4. 更换更高效的拟合算法或库 #

• 切换scipy拟合算法：默认的Levenberg-Marquardt算法在多参数场景下，可尝试trf或dogbox算法，带边界时表现更优：

  popt, pcov = curve_fit(lorentzian_model, x, y, p0=initial_params, bounds=(lower_bounds, upper_bounds), method='trf')
  

• 改用lmfit库：专门针对曲线拟合优化，支持参数分组、边界管理，代码更简洁且速度更快：

  from lmfit import Model, Parameters
  
  # 定义单个Lorentzian函数
  def single_lorentzian(x, amp, cen, wid):
      return amp / (1 + ((x - cen)/wid)**2)
  
  # 构建多峰模型+背景
  model = Model(single_lorentzian, prefix='p1_')
  for i in range(2, 11):
      model += Model(single_lorentzian, prefix=f'p{i}_')
  model += Model(lambda x, bg: bg, name='bg')
  
  # 初始化参数
  params = Parameters()
  for i in range(1, 11):
      params.add(f'p{i}_amp', value=initial_params[3*(i-1)], min=0)
      params.add(f'p{i}_cen', value=initial_params[3*(i-1)+1], min=x.min(), max=x.max())
      params.add(f'p{i}_wid', value=initial_params[3*(i-1)+2], min=0.1, max=10)
  params.add('bg', value=initial_params[-1], min=y.min()-10, max=y.min()+10)
  
  # 执行拟合
  result = model.fit(y, params, x=x)
  popt = list(result.params.valuesdict().values())
  

5. 数据降采样减少计算量 #

如果光谱数据点过多（如上万条），可在保留峰特征的前提下降采样，减少拟合计算量：

# 降采样到2000个点（按需调整）
downsampled_x = np.linspace(x.min(), x.max(), 2000)
downsampled_y = np.interp(downsampled_x, x, y)

# 用降采样数据拟合
popt, pcov = curve_fit(lorentzian_model, downsampled_x, downsampled_y, p0=initial_params, bounds=(lower_bounds, upper_bounds))

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