使用scipy L-BFGS-B优化时出现Line search定位失败问题求助

使用scipy.optimize.minimize()的L-BFGS-B优化方法时，遇到如下报错：

Line search cannot locate an adequate point after MAXLS
function and gradient evaluations.
Previous x, f and g restored.
Possible causes: 1 error in function or gradient evaluation;
2 rounding error dominate computation.

主程序代码 #

import scipy
import scipy.optimize 
import numpy as np
import matplotlib.pyplot as plt
import timeit
from funcion_objetivo import funcion_objetivo
from funcion_gradiente import gradiente

Nx = 210
Nz = 68  
tEnd = 2.5 
dt = 0.004
t = np.arange(0,tEnd,dt)
frec = 3
Nt = tEnd/dt 
Nt= int(Nt)
dz = 25
dx = 25
Sx1 = 25
Sx2 = 185
Sx3 = 105
Sx4 = 65
Sx5 = 145
Sz = 3
offset_max = 4000 
a = (np.pi*frec)**2
t0 = 1
g1 = (-2*a*(t-t0)*np.exp(-a*(t-t0)**2)).T
g1 = np.reshape(g1,(np.size(g1),1))

######### Modelo original
vp_ori = np.load('modelo_original.npy')
######### Modelo inicial
vp_ite = np.load('modelo_inicial1.npy')


start=timeit.default_timer()
p = scipy.optimize.minimize(funcion_objetivo, vp_ite.flatten(), 
                            method='L-BFGS-B', jac=True, args=(vp_ori,),
                            options={'disp':True, 'maxiter':5})
stop=timeit.default_timer()
print('time:', stop-start) 

其中vp_ori和vp_ite是(210, 68)的矩阵，FWI_Grad和propagator为预定义函数。

目标函数funcion_objetivo代码 #

from propagator import propagator
import numpy as np
from funcion_gradiente import gradiente
from numba import jit

@jit
def funcion_objetivo(vp_ite, vp_ori):
    
    Sx1 = 25
    Sx2=185
    Sx3=105
    Sx4=65
    Sx5=145
    Sz=3
    dx=25
    dz=25
    dt=0.004
    frec=3
    a = (np.pi*frec)**2
    t0 = 1
    tEnd = 2.5 
    dt = 0.004
    t = np.arange(0,tEnd,dt)
    g1 = (-2*a*(t-t0)*np.exp(-a*(t-t0)**2)).T
    g1 = np.reshape(g1,(np.size(g1),1))
    
    
    Pt_obs1 = propagator(vp_ori, g1, Sx1, Sz, dx, dz, dt, 4000, frec)
    print('Punto observado 1 done')
    Pt_obs2 = propagator(vp_ori, g1, Sx2, Sz, dx, dz, dt, 4000, frec)
    print('Punto observado 2 done')
    Pt_obs3 = propagator(vp_ori, g1, Sx3, Sz, dx, dz, dt, 2125, frec)
    print('Punto observado 3 done')
    Pt_obs4 = propagator(vp_ori, g1, Sx4, Sz, dx, dz, dt, 3125, frec)
    print('Punto observado 4 done')
    Pt_obs5 = propagator(vp_ori, g1, Sx5, Sz, dx, dz, dt, 3125, frec)
    print('Punto observado 5 done')

    Pt_mod1 = propagator(vp_ite, g1, Sx1, Sz, dx, dz, dt, 4000, frec)
    print('Punto modelado 1 done')
    Pt_mod2 = propagator(vp_ite, g1, Sx2, Sz, dx, dz, dt, 4000, frec)
    print('Punto modelado 2 done')
    Pt_mod3 = propagator(vp_ite, g1, Sx3, Sz, dx, dz, dt, 2125, frec)
    print('Punto modelado 3 done')
    Pt_mod4 = propagator(vp_ite, g1, Sx4, Sz, dx, dz, dt, 3125, frec)
    print('Punto modelado 4 done')
    Pt_mod5 = propagator(vp_ite, g1, Sx5, Sz, dx, dz, dt, 3125, frec)
    print('Punto modelado 5 done')
        
    r1 = np.sum(Pt_mod1[0] - Pt_obs1[0])
    r2 = np.sum(Pt_mod2[0] - Pt_obs2[0])
    r3 = np.sum(Pt_mod3[0] - Pt_obs3[0])
    r4 = np.sum(Pt_mod4[0] - Pt_obs4[0])
    r5 = np.sum(Pt_mod5[0] - Pt_obs5[0])
    
    r = (r1**2 + r2**2 + r3**2 + r4**2 + r5**2)
    norma = 0.5*np.sqrt(r)
    print('la norma es: ', norma)
    
    grad = gradiente(vp_ite, vp_ori)
    print('el gradiente es: ', grad)
    
    return norma, grad

梯度函数gradiente代码 #

from FWI_GRAD import FWI_GRAD
import numpy as np
from numba import jit
from propagator import propagator

@jit
def gradiente(vp_ite, vp_ori):
    Nx = 210
    Nz = 68  
    dz = 25
    Sx1 = 25
    Sx2=185
    Sx3=105
    Sx4=65
    Sx5=145
    Sz=3
    dx=25
    dz=25
    dt=0.004
    frec=3
    a = (np.pi*frec)**2
    t0 = 1
    tEnd = 2.5 
    dt = 0.004
    t = np.arange(0,tEnd,dt)
    g1 = (-2*a*(t-t0)*np.exp(-a*(t-t0)**2)).T
    g1 = np.reshape(g1,(np.size(g1),1))
    dt = 0.004
    frec = 3
    beta = 1
    Nt = tEnd/dt 
    Nt= int(Nt)
    
    Pt_obs1 = propagator(vp_ori, g1, Sx1, Sz, dx, dz, dt, 4000, frec)
    print('Punto observado 1 done')
    Pt_obs2 = propagator(vp_ori, g1, Sx2, Sz, dx, dz, dt, 4000, frec)
    print('Punto observado 2 done')
    Pt_obs3 = propagator(vp_ori, g1, Sx3, Sz, dx, dz, dt, 2125, frec)
    print('Punto observado 3 done')
    Pt_obs4 = propagator(vp_ori, g1, Sx4, Sz, dx, dz, dt, 3125, frec)
    print('Punto observado 4 done')
    Pt_obs5 = propagator(vp_ori, g1, Sx5, Sz, dx, dz, dt, 3125, frec)
    print('Punto observado 5 done')
        
    f1,grad1 = FWI_GRAD(vp_ite, Nx, Nz, Nt, g1, Sx1, Sz, dx, dz, dt, 4000, frec, Pt_obs1[0], 1) 
    f2,grad2 = FWI_GRAD(vp_ite, Nx, Nz, Nt, g1, Sx2, Sz, dx, dz, dt, 4000, frec, Pt_obs2[0], 2) 
    f3,grad3 = FWI_GRAD(vp_ite, Nx, Nz, Nt, g1, Sx3, Sz, dx, dz, dt, 2125, frec, Pt_obs3[0], 3)
    f4,grad4 = FWI_GRAD(vp_ite, Nx, Nz, Nt, g1, Sx4, Sz, dx, dz, dt, 3125, frec, Pt_obs4[0], 4) 
    f5,grad5 = FWI_GRAD(vp_ite, Nx, Nz, Nt, g1, Sx5, Sz, dx, dz, dt, 3125, frec, Pt_obs5[0], 5) 

    grad = beta*grad1 + beta*grad2 + beta*grad3 + beta*grad4 + beta*grad5
    
    return grad

已尝试让目标函数仅返回norma并设置jac=gradiente，但问题仍未解决。请问该错误的可能原因是什么？

1. 目标函数与梯度的一致性问题 #

L-BFGS-B对目标函数和梯度的一致性要求极高，若梯度计算存在偏差，会直接导致线搜索失败：

• 用scipy.optimize.check_grad对比解析梯度和数值梯度的差异，若误差超过1e-4，说明梯度计算存在错误，需检查FWI_GRAD的实现逻辑。
• 当前目标函数的残差计算逻辑是0.5*sqrt(sum((sum(residual))²))，这种形式会降低目标函数的光滑性，建议改成FWI标准的平方和形式：0.5 * (np.sum((Pt_mod1[0]-Pt_obs1[0])**2) + np.sum((Pt_mod2[0]-Pt_obs2[0])**2) + ... + np.sum((Pt_mod5[0]-Pt_obs5[0])**2))，光滑的目标函数更利于线搜索收敛。

2. 重复计算与数值稳定性问题 #

• 目标函数和梯度函数中多次重复计算Pt_obs，不仅浪费计算资源，还可能引入微小数值误差。建议将Pt_obs的计算提前到主程序中，作为参数传入目标函数和梯度函数，避免重复计算。
• 检查propagator和FWI_GRAD函数是否存在数值不稳定情况，比如除以极小值、指数溢出等，可加入数值截断或精度控制逻辑。

3. 线搜索参数设置问题 #

L-BFGS-B的默认线搜索参数可能不匹配当前问题：

• 增大options中的maxls参数（默认20），给线搜索更多尝试次数，比如设置options={'disp':True, 'maxiter':5, 'maxls':50}。
• 适当调整ftol参数，放宽线搜索的终止条件，比如设置ftol=1e-3。

4. 初始模型与变量尺度问题 #

• 若初始模型vp_ite与真实模型vp_ori差异过大，会导致目标函数值过高、梯度方向不稳定，建议先使用更接近真实模型的初始值做测试。
• 变量尺度不一致会影响L-BFGS-B的收敛效果，可对vp_ite做归一化处理（比如除以均值），优化完成后再还原，提升线搜索的稳定性。

5. Numba编译的潜在问题 #

• 目标函数和梯度函数使用了@jit装饰器，需检查Numba编译是否存在类型推断错误或未覆盖的代码路径。可尝试去掉@jit装饰器，用纯Python代码运行测试，排除编译导致的数值异常。

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