代码优化：如何使用scipy minimize优化代码？

Question

我正在尝试优化我为计算方程组的最小二乘法而编写的代码，returns 我是未知数的最佳值：a1,a2,a3,z1,z2（pottemp和 zlevels 都知道）。

我写了下面的代码，它很有效，但它不是很有效，你可以看到由于 for 循环的数量以及如果我增加 hstep 和 [=17 的数量=]，我的数组变得非常大，从 leastsquared.

计算最小二乘法需要很多时间

def leastsquared(zlevels,pottemp,z1,z2,a1,a2,a3):
    dtot=0.0
    for zi in range(len(zlevels)): #zvalue of obs points 
        if zlevels[zi]<=z1:
            F=np.tan(a1)*zlevels[zi]+pottemp[0]
        elif zlevels[zi]<=z2:
            F=(np.tan(a1)*z1+pottemp[0]) + np.tan(a2)*(zlevels[zi]-z1)
        else: #zlevels[zi]<=H
            F=((np.tan(a1)*z1+pottemp[0]) + np.tan(a2)*(z2-z1) ) + np.tan(a3)*(zlevels[zi]-z2)
        d=(pottemp[zi]-F)**2.0
        dtot+=d
    dtot=dtot/len(zlevels)
    #print("dtot is" + str(dtot))
    return dtot

def optm(hstep,astep,zlevels,pottemp):

    sqdist=np.inf
    for z1 in np.linspace(0,zlevels[-1],hstep):
        for z2 in np.linspace(0,zlevels[-1],hstep):
            for a1 in np.linspace(0,0.1,astep): # max angle
                for a2 in np.linspace(0,0.01,astep):
                    for a3 in np.linspace(0,0.01,astep):
                        sqdist_new=leastsquared(zlevels,pottemp,z1,z2,a1,a2,a3)
                        if sqdist_new<sqdist:
                            sqdist=sqdist_new
                            optimalsetting=[z1,z2,a1,a2,a3]

    return optimalsetting

代码本身非常简单但效率低下。我试图使用 scipy.optimize.minimize 来实现此代码，但我无法运行它。有谁知道如何优化它？也许，scipy.optimize.minimize 是最简单的方法，但我不确定如何调整我的代码以适应它。

Answer 1

我不太确定实际变量及其值。但下面是我尝试过的东西：

import numpy as np
from scipy.optimize import minimize

def leastsquared(zlevels,pottemp,z1,z2,a1,a2,a3):
    dtot=0.0
    for zi in range(len(zlevels)): #zvalue of obs points 
        if zlevels[zi]<=z1:
            F = np.tan(a1)*zlevels[zi]+pottemp[0]
        elif zlevels[zi]<=z2:
            F = (np.tan(a1)*z1+pottemp[0]) + np.tan(a2)*(zlevels[zi]-z1)
        else: #zlevels[zi]<=H
            F=((np.tan(a1)*z1+pottemp[0]) + np.tan(a2)*(z2-z1) ) + np.tan(a3)*(zlevels[zi]-z2)
        d=(pottemp[zi]-F)**2.0
        dtot+=d
    dtot=dtot/len(zlevels)
    #print("dtot is" + str(dtot))
    return dtot

def optm(hstep,astep,zlevels,pottemp):
    steps = np.linspace(0,0.01,astep)
    x0 = np.array([steps, steps, steps]) #a1, a2, a3

    def objective(x):
        return leastsquared(zlevels,pottemp, z1, z2, x[0], x[1], x[2])

    sqdist=np.inf
    for z1 in np.linspace(0,zlevels[-1],hstep):
        for z2 in np.linspace(0,zlevels[-1],hstep):
            sol = minimize(objective, x0, method='SLSQP', options={'disp':True})
            optimalsetting=[z1,z2,sol.x[0],sol.x[1],sol.x[2]]
    return optimalsetting


# test
# some random initialization
astep = 2
steps = 3


optm(5, 2, [1, 2], [1, 2])

您可能想在此处更新 objective 函数和变量类型。希望这能让你入门。

Answer 2

肯定有更有效的方法来处理这个问题，因为分段线性问题得到了很好的研究，但下面是根据要求用 scipy.optimize 完成的。 scipy 也无法处理动态约束，例如 z1<z2 的要求，因此如果不强制执行，我会交换求解器中变量的顺序。通过删除这些代码行并让数据说明一切，优化可能 more/less 稳定。

import numpy as np
import scipy.optimize

def f(zlevels, t0, a1, a2, a3, z1, z2, H):
    """
    Function to evaluate f, given a set of variables
    """
    # Check that no samples are out of bounds
    if any(zlevels < 0) or any(zlevels > H):
        quit("Values of zlevels out of bounds")

    pottemp = np.zeros(zlevels.shape)
    mask1 = zlevels <= z1
    mask2 = (z1 < zlevels) & (zlevels <= z2)
    mask3 = z2 < zlevels

    z1_arr1 = np.ones(sum(mask2))*z1
    z1_arr2 = np.ones(sum(mask3))*z1
    z2_arr = np.ones(sum(mask3))*z2

    pottemp[mask1] = zlevels[mask1] * a1 + t0
    pottemp[mask2] = z1_arr1 * a1 + t0 + a2 * (zlevels[mask2] - z1)
    pottemp[mask3] = z1_arr2 * a1 + t0 + a2 * (z2_arr - z1) + a3 * (zlevels[mask3] - z2)

    return pottemp

def obj_fun(x, zlevels, pottemp, t0, H):
    """
    Least-squares objective function for the problem
    """
    a1, a2, a3, z1, z2 = x
    # Swap order if z1 is larger than z2
    if z1 > z2:
        z1, z2 = z2, z1

    pottemp_pred = f(zlevels, t0, a1, a2, a3, z1, z2, H)
    return sum((pottemp_pred-pottemp)**2) 

def optimize(zlevels, pottemp, t0, H):
    """
    Optimize a1, a2, a3, z1, z2
    """

    starting_guess = [0,0,0,0.5,2.5]
    res = scipy.optimize.minimize(obj_fun, starting_guess, args=(zlevels, pottemp, t0, H), \
            bounds=[(None,None),(None,None),(None,None),(0,H),(0,H)])
    x = res.x
    # Swap order if z1 is larger than z2
    if x[-2] > x[-1]:
        x[-2], x[-1] = x[-1], x[-2]
    return x

# Make the true model for testing
t0 = 0.5
a1 = -1
a2 = +1
a3 = -2
z1 = 1
z2 = 2
H = 3

# Generate some variables
n = 1000
zlevels = np.random.random(n)*3
# Get values and add some random noise
pottemp = f(zlevels, t0, a1, a2, a3, z1, z2, H) + np.random.normal(0,0.1, size=n)

print("Correct values:", a1, a2, a3, z1, z2)
a1, a2, a3, z1, z2 = optimize(zlevels, pottemp, t0, H)
print("Fitted values:", a1, a2, a3, z1, z2)

代码优化：如何使用scipy minimize优化代码？

Code optimization: How to optimize a code using scipy minimize?

python

python-3.x

scipy-optimize