import numpy as np
from scipy.optimize import minimize
from numpy.linalg import norm, solve
import logging
np.random.seed(1257)
# Setup our utility function
def u(x, α):
return np.sum(α * np.log(x))
# Make some test parameters
p = np.random.rand(10)
p[0] = 1 # random prices -- normalize
α = np.random.rand(10)
α /= np.sum(α) # ensure that α sum to 1Week 6 — Numerical Methods II — Optimisation in Economics
Programming and Numerical Methods for Economics (ECNM10115) · The University of Edinburgh
Learning goals. Applying the optimisers from Week 5 to economic problems — consumer demand and beyond: setting up the objective, choosing an algorithm, and interpreting the solution.
How to work through this notebook: run every cell in order (
Shift+Enter). When you reach a result, pause and predict it before running — that habit is what turns reading into learning. Experiment: change parameters, break things, re-run.
Consumer Demand
Consider a consumer with an income of \(I\), who is choosing between \(n\) different goods. We encode the prices of these goods as a vector \(p \in \mathbb R^n\). We assume that they have Cobb-Douglass preferences:
# \[ u(x) = \sum_{i=1}^n \alpha_i \log(x_i) \] where \(\alpha\) is a vector of parameters satisfying \(\sum \alpha_i = 1\). This looks like a constrained optimization problem: # \[\begin{equation*} \begin{aligned} \max_{x \in \R^n}\;\; & \sum_{i=1}^n \alpha_i \log(x_i) \\ \text{s.t. } \;\; & p \cdot x \leq I \\ & x_i \geq 0 \end{aligned} \end{equation*}\]
This problem is nice because it has known solutions in closed form, so we can always go back and check our work! If we take the Lagrangian for this problem, we get: # \[ \mathcal L = \sum_{i=1}^n \alpha_i \log(x_i) - \lambda (p \cdot x - I) \] and a set of first order conditions # \[\begin{align*} 0 = {\partial \mathcal L \over \partial x_i} &= {\alpha_i \over x_i} - \lambda p_i \\ 0 = {\partial \mathcal L \over \partial \lambda} &= p \cdot x - I \end{align*}\]
We can show that these imply: # \[ p_i x_i = {\alpha_i \over \lambda} \Rightarrow I = \sum_{i=1}^n p_i x_i = {1 \over \lambda} \sum_{i=1}^n \alpha_i = {1 \over \lambda} \] So we obtain that # \[ x_i = {\alpha_i \over p_i} I \] The consumers spend a constant share of their income \(\alpha_i\) on good \(i\). This should all look very familiar.
The Plan
We are going to solve this problem in three different ways: 1. Sequential Quadratic Programming 2. Penalty Function 3. Augmented Lagrangian
Sequential Quadratic Programming
Remember that SQP is just Newton’s method on the Lagrangian. Let’s recall our update rule for Newton’s method: if \(f: \mathbb R^n \to \mathbb R\), and we want to maximize \(f\), we start with \(x_0\) and update according to # \[ D^2f(x_k) (x_{k+1} - x_{k}) = -Df(x_k) \]
def gradient(f, x, eps=1e-8):
grad = np.zeros_like(x)
fx = f(x)
for i in range(len(x)):
x_eps = x.copy()
x_eps[i] += eps
grad[i] = (f(x_eps) - fx) / eps
return grad
def hessian(f, x, eps=1e-5):
n = len(x)
hess = np.zeros((n, n))
fx = f(x)
for i in range(n):
x_i_eps = x.copy()
x_i_eps[i] += eps
f_i_eps = f(x_i_eps)
for j in range(i, n):
x_ij_eps = x_i_eps.copy()
x_ij_eps[j] += eps
f_ij_eps = f(x_ij_eps)
x_j_eps = x.copy()
x_j_eps[j] += eps
f_j_eps = f(x_j_eps)
hess[i, j] = (f_ij_eps - f_i_eps - f_j_eps + fx) / (eps**2)
hess[j, i] = hess[i, j]
return hess
def newton(f, x0, tol=1e-8, itermax=100, trace=False):
x = x0.copy()
err = np.inf
iter = 0
while not (err < tol) and iter < itermax:
fx = f(x)
Df = gradient(f, x)
D2f = hessian(f, x)
try:
sk = solve(D2f, -Df)
except np.linalg.LinAlgError:
# If Hessian is singular, break
break
x_new = x + sk
err = norm(x_new - x)
x[:] = x_new
iter += 1
if trace:
logging.info(f"Trace Information iter={iter} x={x} fx={f(x)}")
return {'fx': f(x), 'x': x, 'iter': iter}
# Domain transformation helpers
def as_positive_real(x):
# transform from R to positive real via exp
return np.exp(x)
def inverse_as_positive_real(y):
# inverse transform from positive real to R via log
return np.log(y)
def transform(t, z):
# t is a tuple of (transform_func, inverse_func, length)
transform_func, _, length = t
return transform_func(z)
def inverse(t, x):
_, inverse_func, _ = t
return inverse_func(x)
def as_array_as_positive_real(length):
return (as_positive_real, inverse_as_positive_real, length)# Setting up the Lagrangian
def L(x, λ, α, p):
return u(x, α) - λ * (np.dot(x, p) - 1)
# Let's solve the problem:
def h(z):
t = as_array_as_positive_real(11)
y = transform(t, z)
return L(y[:-1], y[-1], α, p)
ret = newton(h, np.ones(11))
ret{'fx': np.float64(-1.29486120494734),
'x': array([-1.92028858e+00, -2.39513058e+00, 6.30204317e-01, -1.78572355e+00,
-2.48017678e+00, -1.35117821e+00, -2.57257838e+00, -8.59301791e-01,
-3.37912533e+00, -8.50913638e-02, 1.89060233e-07]),
'iter': 100}# ## We still need to get our solution back out:
# # Apply the domain transformation again to get
ys = transform(as_array_as_positive_real(11), ret['x'])
xs = ys[:-1]
# Compare to the true solution
# assert np.allclose(xs, α / p)
np.column_stack((xs, α / p))array([[0.14656466, 0.14656467],
[0.09116077, 0.09116077],
[1.87799425, 1.87799412],
[0.16767569, 0.16767569],
[0.08372842, 0.08372838],
[0.258935 , 0.25893501],
[0.07633846, 0.07633843],
[0.42345764, 0.42345757],
[0.03407725, 0.03407728],
[0.91842837, 0.91842848]])Penalty Method
Now, we’re going to try an approach with a penalty function. We will keep using our Newton’s method implementation on the inside.
Recall that the penalty method approach solves a sequence of problems: \[\begin{equation} \max_{x \in \mathbb R^n} f(x) - P_k \sum_{i=1}^m g_i(x)^2 \end{equation}\]
We need a function that takes in the objective, our constraints, and a penalty value \(P_k\). It should return the solved value.
def solve_penalty(P, x0):
t = as_array_as_positive_real(len(x0))
def objective(y):
x = transform(t, y)
return u(x, α) - P * (np.dot(p, x) - 1)**2
y0 = inverse(t, x0)
ret = newton(objective, y0)
ret['x'] = transform(t, ret['x'])
return ret
ret = solve_penalty(1e20, np.ones(10))
np.column_stack((ret['x'], α / p))
ret = solve_penalty(1e12, np.ones(10))
logging.info(f"We took a lot of iterations to get here iter={ret['iter']}")
np.column_stack((ret['x'], α / p))array([[0.16856219, 0.14656467],
[0.16880646, 0.09116077],
[0.19802859, 1.87799412],
[0.16861268, 0.16767569],
[0.17013083, 0.08372838],
[0.1696667 , 0.25893501],
[0.17049524, 0.07633843],
[0.22192187, 0.42345757],
[0.1688856 , 0.03407728],
[0.17814211, 0.91842848]])First, let’s observe that if you start with a very large value of \(P\), the poor conditioning of the problem shows up. We get completely incorrect answers!!
And if we happen to start with a \(P\) that’s large, but not too large, it still doesn’t save us.
We get the right answer, but look at how many iterations it takes!
Let’s try instead to solve the problem with a sequence of penalty values
def penalty_method(x0, tol=1e-12, γ=10, trace=False):
iter = 0
num_evals = 0
P = 1
x = x0.copy()
while True:
ret = solve_penalty(P, x)
num_evals += ret['iter']
x = ret['x']
gx = np.dot(p, x) - 1
if gx < tol:
break
x[:] = ret['x']
P *= γ
iter += 1
if trace:
logging.info(f"Trace iter={iter} P={P} penalty_violation={gx}")
return {'fx': u(x, α), 'x': x, 'iter': iter, 'num_evals': num_evals}
penalty_ret = penalty_method(np.ones(10), trace=True)
np.column_stack((penalty_ret['x'], α / p))array([[0.10952115, 0.14656467],
[0.10635963, 0.09116077],
[1.62785909, 1.87799412],
[0.11548382, 0.16767569],
[0.15034913, 0.08372838],
[0.16950587, 0.25893501],
[0.16153673, 0.07633843],
[2.61919519, 0.42345757],
[0.09716603, 0.03407728],
[0.60538526, 0.91842848]])Let’s check that we got the right answer!
Augmented Lagrangian Approach
Let’s finally try to implement this using the augmented lagrangian method. Recall that now we solve a sequence of problems: \[\begin{equation} \max_{x \in \mathbb \R^n} f(x) - {P_k \over 2} \sum_{i=1}^m g_i(x)^2 - \sum_{i=1}^m \lambda_i^k g_i(x) \end{equation}\] and where we update \(\lambda\) according to \[ \lambda_i^{k+1} = \lambda_i^k + P_k g_i(x_k) \]
# First, we need to solve the inner augmented lagrangian problem
# we'll use newton's method again
def augmented_lagrangian_inner(f, g, P, λ, x0):
def objective(x):
gx = g(x)
return f(x) - 0.5 * P * np.dot(gx, gx) - np.dot(λ, gx)
ret = newton(objective, x0)
return ret
def augmented_lagrangian(f, g, x0, tol=1e-8, γ=10, trace=False):
# setup
gx = g(x0)
λ = np.zeros_like(gx)
P = 1.0
x = x0.copy()
num_evals = 0
iter = 0
while True:
# Solve the inner problem
ret = augmented_lagrangian_inner(f, g, P, λ, x)
x_prime = ret.x
num_evals += ret.iter
iter += 1
# Update λ
gx = g(x_prime)
λ_prime = gx * P
# Check if we're violating the constraints
err = np.dot(g(x), g(x))
if err < tol:
break
else: # otherwise keep going
P *= γ
x[:] = x_prime
λ[:] = λ_prime
if trace:
print(f"Trace iter={iter} P={P} err={err} penalty_violation={gx}")
return {"fx": f(x), "x": x, "λ": λ, "num_evals": num_evals, "iter": iter}# Let's apply now
def obj(y):
t = as_array_as_positive_real(10)
y = transform(t, y)
return u(y, α)
def constraints(y):
t = as_array_as_positive_real(10)
y = transform(t, y)
return np.dot(y, p) - 1
ret = augmented_lagrangian(obj, constraints, np.ones(10), trace=True)
print(ret)--------------------------------------------------------------------------- AttributeError Traceback (most recent call last) Cell In[37], line 12 9 y = transform(t, y) 10 return np.dot(y, p) - 1 ---> 12 ret = augmented_lagrangian(obj, constraints, np.ones(10), trace=True) 13 print(ret) Cell In[35], line 24, in augmented_lagrangian(f, g, x0, tol, γ, trace) 21 while True: 22 # Solve the inner problem 23 ret = augmented_lagrangian_inner(f, g, P, λ, x) ---> 24 x_prime = ret.x 25 num_evals += ret.iter 26 iter += 1 AttributeError: 'dict' object has no attribute 'x'
# Compare values
t = as(Array, as_positive_real, 10)
x = transform(t, ret.x)
hcat(x, α./p)10×2 Matrix{Float64}:
0.0771455 0.0771405
0.306856 0.306836
0.743693 0.743645
0.792526 0.792475
11.5688 11.568
0.323116 0.323095
0.858697 0.858642
0.0722062 0.0722015
0.0659291 0.0659248
0.137085 0.137076—## Before the lab- Try it: change the preference parameters in the consumer problem and check the demand response — does it move the way theory says it should?- Work through PS5 with your group.- Reference if you need the maths: Math Review.Next week: Numerical Methods III — function approximation.