How can I get sco.minimize to give me a solution that isn't the initial guesses?

Question

I have a piece of code that worked well when I optimized advertising budget with 2 variables (channels) but when I added aditional channels, it stopped optimizing with no error messages.

import numpy as np
import scipy.optimize as sco

# setup variables
media_budget = 100000 # total media budget
media_labels = ['launchvideoviews', 'conversion', 'traffic', 'videoviews', 'reach'] # channel names
media_coefs = [0.3524764781, 5.606903166, -0.1761937775, 5.678596017, 10.50445914] # 
# model coefficients
media_drs = [-1.15, 2.09, 6.7, -0.201, 1.21] # diminishing returns
const = -243.1018144

# the function for our model
def model_function(x, media_coefs, media_drs, const):
    # transform variables and multiply them by coefficients to get contributions
    channel_1_contrib = media_coefs[0] * x[0]**media_drs[0]
    channel_2_contrib = media_coefs[1] * x[1]**media_drs[1]
    channel_3_contrib = media_coefs[2] * x[2]**media_drs[2]
    channel_4_contrib = media_coefs[3] * x[3]**media_drs[3]
    channel_5_contrib = media_coefs[4] * x[4]**media_drs[4]

    # sum contributions and add constant
    y = channel_1_contrib + channel_2_contrib + channel_3_contrib + channel_4_contrib + channel_5_contrib + const 

    # return negative conversions for the minimize function to work
    return -y 

# set up guesses, constraints and bounds
num_media_vars = len(media_labels)
guesses = num_media_vars*[media_budget/num_media_vars,] # starting guesses: divide budget evenly

args = (media_coefs, media_drs, const) # pass non-optimized values into model_function

con_1 = {'type': 'eq', 'fun': lambda x: np.sum(x) - media_budget} # so we can't go over budget
constraints = (con_1)

bound = (0, media_budget) # spend for a channel can't be negative or higher than budget
bounds = tuple(bound for x in range(5))

# run the SciPy Optimizer
solution = sco.minimize(model_function, x0=guesses, args=args, method='SLSQP', constraints=constraints, bounds=bounds)

# print out the solution
print(f"Spend: ${round(float(media_budget),2)}\n")
print(f"Optimized CPA: ${round(media_budget/(-1 * solution.fun),2)}")
print("Allocation:")
for i in range(len(media_labels)):
    print(f"-{media_labels[i]}: ${round(solution.x[i],2)} ({round(solution.x[i]/media_budget*100,2)}%)")

And the result is

Spend: $100000.0

Optimized CPA: $-0.0
Allocation:
  -launchvideoviews: $20000.0 (20.0%)
  -conversion: $20000.0 (20.0%)
  -traffic: $20000.0 (20.0%)
  -videoviews: $20000.0 (20.0%)
  -reach: $20000.0 (20.0%)

Which is the same as the initial guesses argument.

Thank you very much!

Update: Following @joni comment, I passed the gradient function explicitly, but still no result. I don't know how to change the constrains to test @chthonicdaemon comment yet.


import numpy as np
import scipy.optimize as sco

# setup variables
media_budget = 100000 # total media budget
media_labels = ['launchvideoviews', 'conversion', 'traffic', 'videoviews', 'reach'] # channel names
media_coefs = [0.3524764781, 5.606903166, -0.1761937775, 5.678596017, 10.50445914] # 
# model coefficients
media_drs = [-1.15, 2.09, 6.7, -0.201, 1.21] # diminishing returns
const = -243.1018144

# the function for our model
def model_function(x, media_coefs, media_drs, const):
    # transform variables and multiply them by coefficients to get contributions
    channel_1_contrib = media_coefs[0] * x[0]**media_drs[0]
    channel_2_contrib = media_coefs[1] * x[1]**media_drs[1]
    channel_3_contrib = media_coefs[2] * x[2]**media_drs[2]
    channel_4_contrib = media_coefs[3] * x[3]**media_drs[3]
    channel_5_contrib = media_coefs[4] * x[4]**media_drs[4]

    # sum contributions and add constant (objetive function)
    y = channel_1_contrib + channel_2_contrib + channel_3_contrib + channel_4_contrib + channel_5_contrib + const 

    # return negative conversions for the minimize function to work
    return -y

# partial derivative of the objective function
def fun_der(x, media_coefs, media_drs, const):
   d_chan1 = 1
   d_chan2 = 1
   d_chan3 = 1
   d_chan4 = 1
   d_chan5 = 1
   
   return np.array([d_chan1, d_chan2, d_chan3, d_chan4, d_chan5])
    
# set up guesses, constraints and bounds
num_media_vars = len(media_labels)
guesses = num_media_vars*[media_budget/num_media_vars,] # starting guesses: divide budget evenly

args = (media_coefs, media_drs, const) # pass non-optimized values into model_function

con_1 = {'type': 'eq', 'fun': lambda x: np.sum(x) - media_budget} # so we can't go over budget
constraints = (con_1)

bound = (0, media_budget) # spend for a channel can't be negative or higher than budget
bounds = tuple(bound for x in range(5))

# run the SciPy Optimizer
solution = sco.minimize(model_function, x0=guesses, args=args, method='SLSQP', constraints=constraints, bounds=bounds, jac=fun_der)

# print out the solution
print(f"Spend: ${round(float(media_budget),2)}\n")
print(f"Optimized CPA: ${round(media_budget/(-1 * solution.fun),2)}")
print("Allocation:")
for i in range(len(media_labels)):
    print(f"-{media_labels[i]}: ${round(solution.x[i],2)} ({round(solution.x[i]/media_budget*100,2)}%)")

If you look at the `solution` variable, you will see a message which explains the termination. In this case it says the constraints are incompatible, but I think it's struggling to find a descent direction which doesn't violate constraints. Might be worth explicitly only optimising the fractional spend for N-1 channels and calculating the last one as 1-sum(the others). — chthonicdaemon, Dec 04 '21 at 04:58
Is the constraint `np.sum(x) == media_budget` really intended? If not, I'd try setting the constraint `np.sum(x) <= media_budget` and using the `trust-constr` method. The latter is more numerically robust in my experience than SLSQP. In addition, it's always a good idea to pass exact derivatives, i.e. at least the objective gradient and hessian. — joni, Dec 04 '21 at 08:45
@chthonicdaemon pardon my ignorance, but how might I explicitly optimize that fractional spend. — toku_mo, Dec 05 '21 at 01:10
@joni Yes, ``np.(sum) == media_budget`` is really intended as the model requires that we consume the whole budget. — toku_mo, Dec 05 '21 at 01:19
Why does this 2 variable example of the same code work, but this 5 var doesn't? https://colab.research.google.com/drive/1IBNRLux1MN6kBu9lNinyGjEPY9IDUBne — toku_mo, Dec 05 '21 at 05:24

score 2 · Accepted Answer · answered Dec 05 '21 at 11:43

The reason you are not able to solve this exact problem turns out to be all about the specific coefficients you have. For the problem as it is specified, the optimum appears to be near allocations where some spends are zero. However, at spends near zero, due to the negative coefficients in media_drs, the objective function rapidly becomes infinite. I believe this is what is causing the issues you are experiencing. I can get a solution with success = True by manipulating the 6.7 to be 0.7 in the coefficients and setting lower bound that is larger than 0 to stop the objective function from exploding. So this isn't so much of a programming issue as a problem formulation issue.

I cannot imagine it would be true that you would see more payoff when you reduce the budget on a particular item, so all the negative powers in media_dirs seem off to me.

I will also post here some improvements I made while debugging this issue. Notice that I'm using numpy arrays more to make some of the functions easier to read. Also notice how I have calculated a correct jacobian:

import numpy as np
import scipy.optimize as sco

# setup variables
media_budget = 100000  # total media budget
media_labels = ['launchvideoviews', 'conversion', 'traffic', 'videoviews', 'reach'] # channel names
media_coefs = np.array([0.3524764781, 5.606903166, -0.1761937775, 5.678596017, 10.50445914]) # 
# model coefficients
media_drs = np.array([-1.15, 2.09, 1.7, -0.201, 1.21]) # diminishing returns
const = -243.1018144

# the function for our model
def model_function(x, media_coefs, media_drs, const):
    # transform variables and multiply them by coefficients to get contributions
    channel_contrib = media_coefs * x**media_drs

    # sum contributions and add constant
    y = channel_contrib.sum() + const 

    # return negative conversions for the minimize function to work
    return -y 

def model_function_jac(x, media_coefs, media_drs, const):
    dy_dx = media_coefs * media_drs * x**(media_drs-1)
    return -dy_dx

# set up guesses, constraints and bounds
num_media_vars = len(media_labels)
guesses = num_media_vars*[media_budget/num_media_vars,] # starting guesses: divide budget evenly

args = (media_coefs, media_drs, const) # pass non-optimized values into model_function

con_1 = {'type': 'ineq', 'fun': lambda x: media_budget - sum(x)} # so we can't go over budget
constraints = (con_1,)

bound = (10, media_budget) # spend for a channel can't be negative or higher than budget
bounds = tuple(bound for x in range(5))

# run the SciPy Optimizer
solution = sco.minimize(
    model_function, x0=guesses, args=args, 
    method='SLSQP', 
    jac=model_function_jac,
    constraints=constraints, 
    bounds=bounds
)

# print out the solution
print(solution)
print(f"Spend: ${round(float(media_budget),2)}\n")
print(f"Optimized CPA: ${round(media_budget/(-1 * solution.fun),2)}")
print("Allocation:")
for i in range(len(media_labels)):
    print(f"-{media_labels[i]}: ${round(solution.x[i],2)} ({round(solution.x[i]/media_budget*100,2)}%)")

This solution at least "works" in the sense that it reports a successful solve and returns an answer different from the initial guess.

Thanks for your answer, I started to dig a bit more and it became apparent that `media_drs` need to be between [0, 1] and any media with a negative `media_coefs` will always result in 0 investment. So it seem this is pretty much a problem based on the input given to the model. — toku_mo, Dec 06 '21 at 18:33

How can I get sco.minimize to give me a solution that isn't the initial guesses?

1 Answers1