最大化复杂多变量函数的计算有效方法

Posted 2023-03-12

【中文标题】最大化复杂多变量函数的计算有效方法

【英文标题】：Computationally efficient way of maximizing a complicated multivariable function

【发布时间】：2021-12-13 10:59:57

【问题描述】：

我一直试图找到一种有效的方法来最大化四个变量中的以下怪物函数，但程序运行需要很长时间，我什至不确定结果是否正确。谁能帮我用 Python 更好地编写代码？ 函数如下：

在哪里

a=[p,q,r,s]。

Y 是在 30 个点采样的测量数据。

这是我的代码。

import numpy as np
import math

Y=Y_t    #Y_t is a predefined column vector with 30 entries.
tstep=0.05               #in s
N=30

cov=np.zeros([30,30])

def R(p,q,r,t):
    om_D=p*np.sqrt(1-q**2)
    return np.pi*r*(np.exp(-q*p*abs(t)))*(np.cos(om_D*t)+(q/(np.sqrt(1-q**2)))*(np.sin(om_D*abs(t))))/(2*q*(p**3))

def I(m,p):
    if m==p:
        return 1
    else:
        return 0

def func(a):   
    a1=a[0]    #natural angular frequency   bounds=[3,20]
    a2=a[1]    #damping ratio               bounds=[0,1]
    a3=a[2]    #psd of forcing signal       bounds=[300,600]
    a4=a[3]    #variance of noise           bounds=[0,0.0001] in m
    #assuming uniform prior for a, we only have to maximise the likelihood function
    
   
    for i in range(30):
        for j in range(30):
            cov[i,j]+=R(a1,a2,a3,(j-i)*tstep)+a4*I(i,j)
    P=((2*np.pi)**(-N/2)) * ((np.linalg.det(cov))**(-0.5)) * np.exp((-0.5) *np.linalg.multi_dot([np.transpose(Y),np.linalg.inv(cov),Y])) 
    return (-1)*P[0]
a_start=[5,0.05,100,0.00001]
bnds=((5,20),(0,1),(300,600),(0,0.0001))
result=spo.differential_evolution(func,bounds=bnds)
print(result.x)  ```

 


  

【问题讨论】：

有趣的问题。但是，我建议对符号进行更多详细说明。 Γ 似乎是一个协方差矩阵。什么是 R_x[ (p-m)Δt |一种] ？条件概率？

您可能需要阅读 toeplitz 矩阵，因为这就是您的协方差，有特别有效的例程可以对它们进行反转和分解。

Γ 是协方差矩阵； |Γ|是它的行列式，R_x 是自相关函数。

我建议你让这个例子可以重现。第二件事 - 您是在寻求另一种方法（更好的优化算法）还是只想加快代码速度、提高性能？

@dankal444 如果有更好的优化算法，请分享。该程序仅需要 30 个数据点即可运行半小时。在未来的实践中，我必须增加数据点和变量的数量。所以，它不应该给我带来麻烦。

【参考方案1】：

cov 初始化存在问题，这就是它不收敛的原因。阻尼比的界限也是一个问题，(0, 1) 现在是(0.0001, 0.999)，比率不应为 0 或 1，因为如果是，则在 R() 中将除以零错误。代码已修复，请参见输出。

代码

import time

import numpy as np
from scipy.optimize import differential_evolution


Y = [[-0.00445551], [-0.01164452], [-0.02171495], [-0.03475491], [-0.00770873], [ 0.0492236 ],
[ 0.07264838], [ 0.03066707], [-0.02457141], [-0.04065968], [-0.01135125], [ 0.02677074], [ 0.06517749],
[ 0.09611112], [ 0.12300657], [ 0.0923581 ], [ 0.03982604], [-0.01473844], [-0.09024497], [-0.14304097],
[-0.17447606], [-0.16926952], [-0.12006193], [-0.00120763], [ 0.11006087], [ 0.19978283], [ 0.24388584],
[ 0.18768875], [ 0.12844553], [ 0.03099409]]    #Y_t is a predefined column vector with 30 entries.
tstep = 0.05               #in s
N = 30


def R(p,q,r,t):
    om_D = p*np.sqrt(1-q**2)
    return np.pi*r*(np.exp(-q*p*abs(t)))*(np.cos(om_D*t)+(q/(np.sqrt(1-q**2)))*(np.sin(om_D*abs(t))))/(2*q*(p**3))


def I(m,p):
    if m==p:
        return 1
    else:
        return 0


def func(a):
    cov=np.zeros([N,N])

    a1=a[0]    #natural angular frequency   bounds=[3,20]
    a2=a[1]    #damping ratio               bounds=[0,1]
    a3=a[2]    #psd of forcing signal       bounds=[300,600]
    a4=a[3]    #variance of noise           bounds=[0,0.0001] in m
    #assuming uniform prior for a, we only have to maximise the likelihood function
   
    for i in range(N):
        for j in range(N):
            cov[i,j]+=R(a1,a2,a3,(j-i)*tstep)+a4*I(i,j)
    P=((2*np.pi)**(-N/2)) * ((np.linalg.det(cov))**(-0.5)) * np.exp((-0.5) *np.linalg.multi_dot([np.transpose(Y),np.linalg.inv(cov),Y])) 
    return (-1)*P[0]


if __name__ == '__main__':
    t0 = time.perf_counter()

    a_start = [5, 0.05, 350, 0.00001]
    bnds = ((5, 20), (0.0001, 0.999), (300, 600), (0, 0.0001))

    result=differential_evolution(func, x0=a_start, bounds=bnds, maxiter=1000)
    print(result)

    print(f'elapse: time.perf_counter() - t0:0.0fs')

输出

     fun: array([-2.76736878e+11])
     jac: array([-2.91459845e+11, -4.55652161e+12,  1.27377279e+10,  3.34234132e+14])
 message: 'Optimization terminated successfully.'
    nfev: 3430
     nit: 56
 success: True
       x: array([ 20.   ,   0.999, 300.   ,   0.   ])
elapse: 55s

Scipy 最小化速度非常快

变化：

from scipy.optimize import minimize

result = minimize(func, x0=a_start, bounds=bnds, options='maxiter': 100, 'disp': True)

输出：

      fun: array([-2.76736878e+11])
 hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>
      jac: array([-2.91459845e+11, -4.55652161e+12,  1.27377279e+10,  3.34234132e+14])
  message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'
     nfev: 30
      nit: 4
     njev: 6
   status: 0
  success: True
        x: array([ 20.   ,   0.999, 300.   ,   0.   ])
elapse: 0.5s

奥图纳

Optuna 在 1000 次试验后也在那里。 这里的值是正的，因为我使用最大化方向。在 scipy 的 DE 和最小化值中都必须取反。

best param: 'a1': 20, 'a2': 0.9989999999999999, 'a3': 300, 'a4': 0.0
best value: 276736878140.3103
best trial num: 73
elapse: 22s