Problem set 1: Optimization solvers results based on initial guess

I couldn't raise an issue in the PS1 exercise repository, so I am creating it here instead. Feel free to (re)move the post if necessary.

My issue/question is regarding the solvers in PS1: As the utility is dependent on the initial guess how should we define the initial guess so it is robust to a large variance in the denominator?

To illustrate consider the following:

from scipy import optimize

def utility_function(x,alpha):
    u = 1
    for x_now,alpha_now in zip(x,alpha):
        u *= np.max(x_now,0)**alpha_now
    return u

def expenditures(x,p):
    E = 0
    for x_now,p_now in zip(x,p):
        E += p_now*x_now
    return E

def print_solution(x,alpha,I,p):

    # a. x values
    text = 'x = ['
    for x_now in x:
        text += f'{x_now:.2f} '
    text += f']\n'

    # b. utility
    u = utility_function(x,alpha)    
    text += f'utility = {u:.3f}\n'

    # c. expenditure vs. income
    E =  expenditures(x,p)
    text += f'E = {E:.2f} <= I = {I:.2f}\n'

    # d. expenditure shares
    e = p*x/I
    text += 'expenditure shares = ['
    for e_now in e:
        text += f'{e_now:.2f} '
    text += f']'        

    print(text)

alpha = np.ones(5)/5
p = np.array([1,2,3,4,5])
I = 10

Solution using a constrained optimizer:

# a. contraint function (negative if violated)
constraints = ({'type': 'ineq', 'fun': lambda x:  I-expenditures(x,p)})
bounds = [(0,I/p_now) for p_now in p]

# b. call optimizer
initial_guess = (I/p)/100# some guess
res = optimize.minimize(
    lambda x: -utility_function(x,alpha),initial_guess,
    method='SLSQP',bounds=bounds,constraints=constraints)

# c. print result
print_solution(res.x,alpha,I,p)

Changing initial_guess = (I/p)/100 to initial_guess = (I/p)/1 changes the utility from 0.768 to 0.000.

Conversely, for the solution using an unconstrained optimizer:

# a. define objective function
def unconstrained_objective(x,alpha,I,p):

    penalty = 0
    E = expenditures(x,p)
    if E >= I:
        ratio = I/E
        x *= ratio # now p*x = I
        penalty = 1000*(E-I)**2

    u = utility_function(x,alpha)
    return -u + penalty

# b. call optimizer
initial_guess = (I/p)/10000
res = optimize.minimize(
    unconstrained_objective,initial_guess,
    method='Nelder-Mead',args=(alpha,I,p))

# c. print result
print_solution(res.x,alpha,I,p)

Changing the code from initial_guess = (I/p)/1 to initial_guess = (I/p)/10000 changes the utility from 0.768 to 0.000

Looking instead at the loops we consider the code snippet using raw loops:

N = 15 # number of points in each dimension
fac = np.linspace(0,1,N) # vector betweein 0 and 1
x_max = I/p # maximum x so E = I

u_best = -np.inf
x_best = np.empty(5)
for x1 in fac:
    for x2 in fac:
        for x3 in fac:
            for x4 in fac:
                for x5 in fac:
                    x = np.array([x1,x2,x3,x4,x5])*x_max
                    E = expenditures(x,p)
                    if E <= I:
                        u_now = utility_function(x,alpha)
                        if u_now > u_best:
                            x_best = x
                            u_best = u_now

print_solution(x_best,alpha,I,p)

If we change the original value N = 15 to N = 6 the utility changes from 0.758 to 0.768 (same as the the solvers above), however changing it to N = 5 yields a utility of 0.000. How do we determine the "optimal" or correct number of points in each dimension for our linespacing?

NumEconCopenhagen / lectures-2019

Problem set 1: Optimization solvers results based on initial guess #5