Problems with gradient definition in Python

[SBolo]

Hi, I'm trying to use NLOpt to minimize a function of the kind

T = sum_k sum_n [ R_kn - sum_i w_i R_kni ]^2

with respect to parameters wi, under the constraints 0 <= w_i < 1 and sum_i w_i =1. Numpy arrays R_kn and R_kni are known since the beginning of the minimization and are just accessed from the function to be minimized. This is the code I'm using to do so:

import numpy as np
import nlopt

def exp_data( K, N, mu, sigma ):

    R_exp = np.empty( (K, N) )
    
    for k in range(K):
        for n in range(N):
            R_exp[k,n] = np.random.normal( mu, sigma )
            
    return R_exp

def md_data( R_exp, B ):
    
    K = R_exp.shape[0]
    N = R_exp.shape[1]
    R_md = np.empty( (K, N, B) )
    
    for k in range(K):
        for n in range(N):

            mu_noise = np.random.uniform(-10,10)
            sigma_noise = np.random.uniform(1,10)
            av = R_exp[k,n] + np.random.normal(mu_noise, sigma_noise, size = B)
            rands = np.random.uniform( size = B )
            s = np.mean(rands)
            rands = rands / s * av
            
            R_md[k,n,:] = rands 
            
    return R_md

def residuals_nlopt( w, grad, R_exp, R_md ):
    
    R_md_av = np.average( R_md, weights = w, axis = 2 )
    T = np.sum( (R_exp - R_md_av) **2 )
    
    if grad.size > 0:
        
        K = R_exp.shape[0]
        N = R_exp.shape[1]
        
        for i in range(len(grad)):
            grad[i] = 0
            for k in range(K):
                for n in range(N):
                    grad[i] -= 2. * ( R_exp[k,n] - R_md_av[k,n] ) * R_md[k,n,i]
            
    return T

def normalization( w, grad ):
    
    if grad.size > 0:
        grad[:] = np.arange(0, len(w), 1)
        
    return np.sum(w) - 1

def call_nlopt( w0, R_exp, R_md ):
    
    opt = nlopt.opt( nlopt.LD_SLSQP, B )
    
    lb = [0] * B
    ub = [1] * B
    
    opt.set_lower_bounds(lb)
    opt.set_upper_bounds(ub)

    opt.set_min_objective( lambda x,grad: residuals_nlopt(x,grad,R_exp,R_md) )
    opt.add_equality_constraint( normalization, 1e-10 )
    opt.set_xtol_rel( 1e-4 )
    x = opt.optimize(w0)
    minf = opt.last_optimum_value()
    
    print( "minimum value = ", minf )
    print( "result code = ", opt.last_optimize_result() )
    
    return x

K = 1   
N = 47  
B = 37
mu = 80.
sigma = 5.

R_exp = exp_data( K, N, mu, sigma )
R_md = md_data( R_exp, B )                 
w0 = np.full( B, 1./B )                       
T = residuals_nlopt( w0, np.array([]), R_exp, R_md )

print("Initial residuals: ", T)
x = call_nlopt( w0, R_exp, R_md )

exp_data and md_data are used to generate synthetic data, and it's not really important how they do that. In this code, the minimization doesn't work: the value of T after the minimization is identical to the initial one, even though the minimization returns code 4. I then tried COBYLA, which seems to work better and with which I can get the same result I obtain by employing the scipy implementation of SLSQP. Since the scipy version is working, I believe the numerical formulation of my gradient should be correct, but because COBYLA is working and the NLOpt version of SLSQP is not, I guess there might be something wrong with the implementation which I cannot understand. I read the documentation and I could't find the issue (I'm using the [:] operator to assign new gradient values as suggested). Can you please help understand what I'm doing wrong? Moreover, do you have any suggestions on what could be the best algorithm to solve this kind of minimization problem?

[stevengj]

You might try

f = lambda x,grad: residuals_nlopt(x,grad,R_exp,R_md)

and check that if you pass in an x and grad array to f(x, grad) then upon output grad is set correctly.

The most common bug is simply that your gradient is incorrect — it is easy to get this wrong. Usually it is a good idea to validate each component of the grad output against finite-difference approximations.