Computationally efficient way of maximizing a complicated multivariable function

user_13992974 0
Viewed 134

I've been trying to find an efficient way of maximizing the following monster function in four variables but the program is taking ages to run and I'm not even sure if the results are correct. Can anyone help me code it better in Python? Here's the function:

enter image description here

where

enter image description here

enter image description here

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

enter image description here

Y is the measured data sampled at 30 points.

Here's my code.

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)  ```
scipy math matrix optimization python
1 Answers

There is an issue in cov initialization that is why it does not converge. Also an issue on bound for damping ratio, was (0, 1) now (0.0001, 0.999) the ratio should not be 0 or 1 because if it is there will be division by zero error in R(). Code is fixed now see also the output.

Code

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.0f}s')

Output

     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 minimize is very fast

Changes:

from scipy.optimize import minimize

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

Output:

      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

Optuna after 1000 trials is right there too. value is positive here because I use maximize direction. In both scipy's DE and minimize values have to be negated.

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