CODEWITHSUNDEEP
pythonnumpyoptimizationscipyscipy-optimize-minimize
GarfieldH
GarfieldH

Reputation: 33

The minimize function of scipy.optimize doesn't give the right answer

I'm trying the solve a minimization problem using the minimize function of Scipy. The objective function is simply the ratio of two multivariate normal distributions with different mean and variance. I'm hoping to find the maximum of the function g_func, which is equivalent to find the minimum of the function g_optimization. Also, I added a constraint of x[0] = 0. Here, x is a vector with 8 elements. The objective function g_optimization is as following:

import numpy as np
from scipy.optimize import minimize
# Set up mean and variance for two MVN distributions
n_trait = 8
sigma = np.full((n_trait, n_trait),0.0005)
np.fill_diagonal(sigma,0.005)
omega = np.full((n_trait, n_trait),0.0000236)
np.fill_diagonal(omega,0.0486)
sigma_pos = np.linalg.inv(np.linalg.inv(sigma)+np.linalg.inv(omega))
mu_pos = np.array([-0.01288244,0.08732091,0.01049617,0.0860966,0.10055626,0.07952922,0.04363669,-0.0061975])
mu_pri = 0
sigma_pri = omega

#objective function
def g_func(beta,mu_sim_pos):
    g1 = ((np.linalg.det(sigma_pri))**(1/2))/((np.linalg.det(sigma_pos))**(1/2))
    g2 = (-1/2)*np.linalg.multi_dot([np.transpose(beta-mu_sim_pos),np.linalg.inv(sigma_pos),beta-mu_sim_pos])
    g3 = (1/2)*np.linalg.multi_dot([np.transpose(beta-mu_pri),np.linalg.inv(sigma_pri),beta-mu_pri])
    g = g1*np.exp(g2+g3)
    return g
def g_optimization(beta,mu_sim_pos):
    return -1*g_func(beta,mu_sim_pos)

#optimization
start_point = np.full(8,0)
cons = ({'type': 'eq',
         'fun' : lambda x: np.array([x[0]])})

anws = minimize (g_optimization, [start_point], args=(mu_pos),
                 constraints=cons, options={'maxiter': 50}, tol=0.001)
anws

The optimization stops after two iterations, and the minimum value that the function gives is 0, at the point np.array([0,10.32837891,-1.62396508,10.13790152,12.38752653,9.11615259,3.53201544,-4.22115517]). This cannot be true because even we plug in the starting point np.zeros(8) to the g_optimization function, the result given is -657.0041125829354, which is smaller than 0. So the solution provided is definitely not minimal.

g_optimization(np.zeros(8),mu_pos) #gives solution of -657.0041125829354

I'm not sure where did I go wrong.

Upvotes: 3

Views: 1771

Answers (1)

tenhjo
tenhjo

Reputation: 4537

I would try a different solver. For example L-BFGS-B works well. You can look at all options here.

anws = minimize (g_optimization, [start_point], args=(mu_pos), method='L-BFGS-B',
                 constraints=cons, options={'maxiter': 50}, tol=0.001)

print(anws)

# success: True
# message: b'CONVERGENCE: REL_REDUCTION_OF_F_<=_FACTR*EPSMCH'
# fun: -21688.00879938617
# x: array([-0.0101048, 0.09937778, 0.01543875, 0.0980401, 0.11383878, 0.09086455, 0.05164822, -0.00280081])

EDIT: L-BFGS-B can not handle general constraints h(x)=0, only bounding boxes on the variables:

Bounds on variables for L-BFGS-B, TNC, SLSQP, Powell, and trust-constr methods. There are two ways to specify the bounds: Instance of Bounds class. Sequence of (min, max) pairs for each element in x. None is used to specify no bound.

In your case you have to define 8 pairs of lower and upper limits. For x[0] you have to make a tight bound as the method can not handle x_low == x_high.

bounds = [(None, None)] * 8
bounds[0] = (0, 0.00001)
anws = minimize (g_optimization, [start_point], args=(mu_pos), method='L-BFGS-B', bounds=bounds,
                 options={'maxiter': 50}, tol=0.001)
# fun: -21467.48153792194
# x: array([0., 0.10039832, 0.01641271, 0.0990599, 0.11486735, 0.09188037, 0.05264228, -0.00183697])

Another alternative is to exclude the value x[0] from your optimisation problem:

def g_optimization(beta,mu_sim_pos):
    beta2 = np.empty(8)
    beta2[0] = 0
    beta2[1:] = beta
    return -1*g_func(beta2, mu_sim_pos)

start_point = np.zeros(7)  # exclude x[0]

anws = minimize(g_optimization, [start_point], args=(mu_pos), method='L-BFGS-B',
                 options={'maxiter': 50}, tol=0.001)
# fun: -21467.47686079844
# x: array([0.10041797, 0.01648995, 0.09908046, 0.11487707, 0.09190585, 0.05269467, -0.00174722])
#          ^ missing x[0]

Upvotes: 2

Related Questions