Algorithms Homework 4: Interior Point Methods
Contents
Algorithms Homework 4: Interior Point Methods#
Tips and Tricks#
Background#
Problem Formulation#
Consider the following nonlinear program:
\[\begin{split}
\begin{align}
\min_{x} \quad & f(x) \\
\mathrm{s.t.} \quad & c(x) = 0 \\
& x_i \geq 0, \quad i \in \mathcal{I}
\end{align}
\end{split}\]
where \(x \in \mathbb{R}^{n}\), \(f(x): \mathbb{R}^{n} \rightarrow \mathbb{R}\), \(c(x): \mathbb{R}^{n} \rightarrow \mathbb{R}^{m}\) and \(|\mathcal{I} |\) = \(r\) (i.e., there are \(r\) variables with a lower bound and \(r \leq n\)). This is an extension of (6.48) in Biegler (2010).
This has the corresponding log-barrier approximation:
\[\begin{split}
\begin{align}
\min_{x} \quad & \phi_{\mu_l}(x) := f(x) - \mu_l \sum_{i \in \mathcal{I}} \log(x_i) \\
\mathrm{s.t.} \quad & c(x) = 0
\end{align}
\end{split}\]
which is an extension of (6.49) in Biegler (2010).
Let the \(n \times r\) matrix \(G\) encode which variables are bounded. If variable \(i\) corresponds to the \(j\)th bound, then \(G_{i,j} = 1\) and otherwise \(G_{i,j} = 0\). \(G\) is assembled as follows:
Notice that \(G\) is the gradient of \(x \geq 0\). \(G^T G = I\) but the converse does not hold unless \(r = n\).
Reformulation Example#
Start with:
\[\begin{split}
\begin{align}
\min_{x} \quad & x_1^2 + x_2 \\
\mathrm{s.t.} \quad & x_1 + x_2 = 1 \\
& x_1 + 1 \geq 0
\end{align}
\end{split}\]
Add slack variable \(x_3\) and convert the inequality constraint to an equality constraint and bound:
\[\begin{split}
\begin{align}
\min_{x} \quad & x_1^2 + x_2 \\
\mathrm{s.t.} \quad & x_1 + x_2 = 1 \\
& x_1 + 1 - x_3 = 0 \\
& x_3 \geq 0
\end{align}
\end{split}\]
Now assemble \(G\):
\[\begin{split}
G = \begin{bmatrix}
0 \\
0 \\
1
\end{bmatrix}
\end{split}\]
Primal Dual Optimality Conditions#
\[
\newcommand{\dims}[3]{\underbrace{#1}_{#2 \times #3}}
\newcommand{\I}[0]{\mathcal{I}}
\newcommand{\dfx}[0]{\dims{\nabla f(x^k)}{n}{1}}
\newcommand{\cx}[0]{\dims{c(x^k)}{m}{1}}
\newcommand{\dcx}[0]{\dims{\nabla c(x^k)}{n}{m}}
\newcommand{\dcxT}[0]{\dims{\nabla_x c(x^k)^T}{m}{n}}
\newcommand{\W}[0]{\dims{W^k}{n}{n}}
\newcommand{\G}[0]{\dims{G}{n}{r}}
\newcommand{\GT}[0]{\dims{G^T}{r}{n}}
\newcommand{\vk}[0]{\dims{v^k}{m}{1}}
\newcommand{\Xinv}[0]{\dims{(\hat{X}^k)^{-1}}{r}{r}}
\newcommand{\X}[0]{\dims{\hat{X}^k}{r}{r}}
\newcommand{\e}[0]{\dims{e}{r}{1}}
\newcommand{\uk}[0]{\dims{u^k}{r}{1}}
\newcommand{\dx}[0]{\dims{d_x^k}{n}{1}}
\newcommand{\dv}[0]{\dims{d_v^k}{m}{1}}
\newcommand{\du}[0]{\dims{d_u^k}{r}{1}}
\newcommand{\Uk}[0]{\dims{U^k}{r}{r}}
\newcommand{\Sk}[0]{\dims{\Sigma^k}{n}{n}}
\newcommand{\dphi}[0]{\dims{\nabla \phi_{\mu_l}}{n}{1}}
\]
Next we extend (6.51) in Biegler (2010):
\[\begin{split}
\begin{gather}
\dfx ~+~ \dcx \vk ~-~ \G \uk = \dims{0}{n}{1} \\
\X \uk = \mu \e \\
\cx = \dims{0}{m}{1}
\end{gather}
\end{split}\]
We now extend (6.56) in Biegler (2010):
\[\begin{split}
\begin{bmatrix}
\W & \dcx & \G \\
\dcxT & \dims{0}{m}{m} & \dims{0}{m}{r} \\
\Uk \GT & \dims{0}{r}{m} & \X
\end{bmatrix}
\begin{bmatrix}
\dx \\
\dv \\
\du
\end{bmatrix} = -
\begin{bmatrix}
\dfx ~+~ \dcx \vk ~-~ \G \uk \\
\cx \\
\X \uk ~-~ \mu_l \e
\end{bmatrix}
\end{split}\]
where \(U^k = \mathrm{diag}\{u^k\}\) and \(W^k = \nabla_{xx} L(x^k,v^k)\). Notice that \(W^k\) does NOT include a contribution from the barrier term:
\[
\begin{equation}
\W = \dims{\nabla^2 f(x^k)}{n}{n} + \sum_{j=1}^{m} \left( \dims{\nabla^2 c_{j}(x^k)}{n}{n} \underbrace{v^k_{j}}_{\mathrm{scalar}} \right)
\end{equation}
\]
We can verify that \(W^k\) does not include \(\hat{X}\) by showing that the KKT system above is a Newton step to solve the nonlinear system for the primal dual conditions.
The Newton step can be simplified, similar to (6.57) and (6.58) in Biegler (2010):
\[\begin{split}
\begin{equation}
\begin{bmatrix}
\W + \Sk & \dcx \\
\dcxT & 0
\end{bmatrix}
\begin{bmatrix}
\dx \\
\dv
\end{bmatrix} = -
\begin{bmatrix}
\dphi ~+~ \dcx \vk \\
\cx
\end{bmatrix}
\end{equation}
\end{split}\]
and
\[
\begin{equation} \label{eq:du}
\du = \mu_l \Xinv \e ~-~ \uk ~-~ \GT \Sk \dx
\end{equation}
\]
where $$
(2)#\[\begin{equation} \label{eq:Sk}
\Sk = \G \Xinv \Uk \GT
\end{equation}\]
\[and \]
(3)#\[\begin{equation}
\dphi = \dfx - \mu_l \G \Xinv \dims{e}{r}{1}
\end{equation}\]
\[Notice that the equation for $\du$ simplifies by substituting $G^T G = I$:\]
\du = \mu_l \Xinv \e ~-~ \uk ~-~ \Xinv \Uk \GT \dx
$$
Finally, inertia correction can be applied to simplified KKT step similar to (6.59) in Biegler (2010).
Basic Interior Point Method for Inequality and Equality Constraint NLPs#
Implement a basic interior point method for inequality and equality constrained nonlinear programs. See pg. 154 – 155 in Biegler (2010). You may skip the line search, i.e., always take a full step.
Pseudocode#
Write detailed pseudocode on paper or a whiteboard. Scan/take a photo and turn in.
Python Implementation#
Implement in Python. Hints: Reuse code from Algorithm 5.2 example.
### Load Python libraries
import numpy as np
from scipy import linalg
### Define helper functions
## Check is element of array is NaN
def check_nan(A):
return np.sum(np.isnan(A))
## Calculate gradient with central finite difference
def my_grad_approx(x,f,eps1,verbose=False):
'''
Calculate gradient of function f using central difference formula
Inputs:
x - point for which to evaluate gradient
f - function to consider
eps1 - perturbation size
Outputs:
grad - gradient (vector)
'''
n = len(x)
grad = np.zeros(n)
if(verbose):
print("***** my_grad_approx at x = ",x,"*****")
for i in range(0,n):
# Create vector of zeros except eps in position i
e = np.zeros(n)
e[i] = eps1
# Finite difference formula
my_f_plus = f(x + e)
my_f_minus = f(x - e)
# Diagnostics
if(verbose):
print("e[",i,"] = ",e)
print("f(x + e[",i,"]) = ",my_f_plus)
print("f(x - e[",i,"]) = ",my_f_minus)
grad[i] = (my_f_plus - my_f_minus)/(2*eps1)
if(verbose):
print("***** Done. ***** \n")
return grad
## Calculate gradient with central finite difference
def my_jac_approx(x,h,eps1,verbose=False):
'''
Calculate Jacobian of function h(x) using central difference formula
Inputs:
x - point for which to evaluate gradient
h - vector-valued function to consider. h(x): R^n --> R^m
eps1 - perturbation size
Outputs:
A - Jacobian (n x m matrix)
'''
# Check h(x) at x
h_x0 = h(x)
# Extract dimensions
n = len(x)
m = len(h_x0)
# Initialize Jacobian matrix
A = np.zeros((n,m))
# Calculate Jacobian by row
for i in range(0,n):
# Create vector of zeros except eps in position i
e = np.zeros(n)
e[i] = eps1
# Finite difference formula
my_h_plus = h(x + e)
my_h_minus = h(x - e)
# Diagnostics
if(verbose):
print("e[",i,"] = ",e)
print("h(x + e[",i,"]) = ",my_h_plus)
print("h(x - e[",i,"]) = ",my_h_minus)
A[i,:] = (my_h_plus - my_h_minus)/(2*eps1)
if(verbose):
print("***** Done. ***** \n")
return A
## Calculate Hessian using central finite difference
def my_hes_approx(x,grad,eps2):
'''
Calculate gradient of function my_f using central difference formula and my_grad
Inputs:
x - point for which to evaluate gradient
grad - function to calculate the gradient
eps2 - perturbation size (for Hessian NOT gradient approximation)
Outputs:
H - Hessian (matrix)
'''
n = len(x)
H = np.zeros([n,n])
for i in range(0,n):
# Create vector of zeros except eps in position i
e = np.zeros(n)
e[i] = eps2
# Evaluate gradient twice
grad_plus = grad(x + e)
grad_minus = grad(x - e)
# Notice we are building the Hessian by column (or row)
H[:,i] = (grad_plus - grad_minus)/(2*eps2)
return H
## Linear algebra calculation
def xxT(u):
'''
Calculates u*u.T to circumvent limitation with SciPy
Arguments:
u - numpy 1D array
Returns:
u*u.T
Assume u is a nx1 vector.
Recall: NumPy does not distinguish between row or column vectors
u.dot(u) returns a scalar. This functon returns an nxn matrix.
'''
n = len(u)
A = np.zeros([n,n])
for i in range(0,n):
for j in range(0,n):
A[i,j] = u[i]*u[j]
return A
## Analyze Hessian
def analyze_hes(B):
print(B,"\n")
l = linalg.eigvals(B)
print("Eigenvalues: ",l,"\n")
## Assemble KKT matrix (equality constrained only)
def assemble_check_KKT(W,Sk,A,deltaA,deltaW,verbose):
# Add your solution here
return KKT,inertia_correct,pos_ev,neg_ev,zero_ev
def barrier_subproblem(x0,v0,u0,calc_f,calc_c,var_bounds,mu,eps,max_iter=100,verbose=False):
'''
Basic Full Space Newton Method for Solving Barrier Subproblem Constrained NLP
Input:
x0 - starting point (vector)
calc_f - function to calculate objective (returns scalar)
calc_c - function to calculate constraints (returns vector)
var_bounds - list of indicies for variables with lower bound
mu - barrier penalty
eps - tolerance for termination
Histories (stored for debugging):
x - history of steps (primal variables)
v - history of steps (duals for constraints)
u - history of steps (duals for bounds)
f - history of objective evaluations
c - history of constraint evaluations
df - history of objective gradients
dL - history of Lagrange function gradients
A - history of constraint Jacobians
W - history of Lagrange Hessians
S - history of sigma matrix
Outputs:
x - final value for primal variable
v - final value for constraint duals
u - final value for bound duals
Notes:
1. For simplicity, central finite difference is used
for all gradient calculations.
'''
### Specifics for Algorithm 5.2
# Tuning parameters
delta_bar_W_min = 1E-20
delta_bar_W_0 = 1E-4
delta_bar_W_max = 1E40
delta_bar_A = 1E-8
kappa_u = 8
kappa_l = 1/3
# Declare iteration histories as empty lists
x = []
v = []
u = []
f = []
L = []
c = []
df = []
dL = []
A = []
W = []
# Add your solution here
def interior_point(x0,calc_f,calc_c,var_bounds,max_iter=20):
# Add your solution here
return x, v, u, mu, E
Test Problems#
Problem 1: Convex#
f = lambda x: x[0] + 2*x[1]
c = lambda x: (x[0] + x[1] - 1)*np.ones(1)
# Indices of variables with lower bound of zero
vb = [1]
x0 = np.ones(2)
u0 = np.ones(1)
v0 = np.ones(1)
x_, v_, u_, E_ = barrier_subproblem(x0,v0,u0,f,c,vb,1E-1,1E-10,verbose=False)
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 3.0000e+00 1.0000e+00 2.00e+00 9.06e-01 2.00e+00 3.56e-04 0.00e+00 0.00e+00
1 1.0996e+00 1.3978e-10 3.20e-04 5.03e-04 1.78e-04 3.56e-04 0.00e+00 0.00e+00
2 1.1000e+00 3.9524e-14 1.26e-07 1.70e-07 3.94e-08 5.91e-08 0.00e+00 0.00e+00
3 1.1000e+00 0.0000e+00 6.87e-11 1.48e-11 1.11e-10 1.04e-10 0.00e+00 0.00e+00
Problem 2: Convex#
f = lambda x: x[0] + 2*x[1]
c = lambda x: (x[0] + x[1] - 1)*np.ones(1)
# Indices of variables with lower bound of zero
vb = [1]
x0 = np.ones(2)
#u0 = np.ones(1)
#v0 = np.ones(1)
# x_, v_, u_, E_ = barrier_subproblem(x0,v0,u0,f,c,vb,1E-1,1E-10,verbose=False)
x, v, u, mu, E = interior_point(x0,f,c,vb)
*** Barrier Subproblem 0 ***
mu = 10.0
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 3.0000e+00 1.0000e+00 9.00e+00 1.35e+01 2.00e+00 8.45e-03 0.00e+00 0.00e+00
*** Barrier Subproblem 1 ***
mu = 2.0
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.1008e+01 1.3978e-10 7.92e+00 1.14e+01 4.22e-03 8.45e-03 0.00e+00 0.00e+00
*** Barrier Subproblem 2 ***
mu = 0.4
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 2.9318e+00 0.0000e+00 1.53e+00 2.17e+00 1.70e-04 3.40e-04 0.00e+00 0.00e+00
*** Barrier Subproblem 3 ***
mu = 0.08000000000000002
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.3993e+00 2.5520e-10 3.19e-01 4.51e-01 1.70e-04 3.40e-04 0.00e+00 0.00e+00
*** Barrier Subproblem 4 ***
mu = 0.016000000000000004
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0801e+00 5.3163e-11 6.41e-02 9.07e-02 1.07e-05 1.78e-05 0.00e+00 0.00e+00
*** Barrier Subproblem 5 ***
mu = 0.0020238577025077633
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0160e+00 0.0000e+00 1.40e-02 1.98e-02 1.15e-05 1.86e-05 0.00e+00 0.00e+00
*** Barrier Subproblem 6 ***
mu = 9.104790579399288e-05
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0020e+00 7.7571e-13 1.93e-03 2.73e-03 7.76e-07 8.83e-07 0.00e+00 0.00e+00
1 1.0001e+00 0.0000e+00 1.07e-07 2.43e-09 0.00e+00 1.07e-07 0.00e+00 0.00e+00
*** Barrier Subproblem 7 ***
mu = 8.687702517211205e-07
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0001e+00 0.0000e+00 9.02e-05 1.28e-04 1.11e-10 5.12e-09 0.00e+00 0.00e+00
1 1.0000e+00 9.9920e-15 4.95e-09 6.52e-13 3.88e-17 4.95e-09 0.00e+00 0.00e+00
import matplotlib.pyplot as plt
# number of iterations
N = len(x)
# iteration
iters = range(0,N)
plt.figure
for i in range(0,len(x0)):
plt.plot(iters,[x[j][i] for j in range(0,N)],label="x"+str(i))
plt.xlabel("Iteration")
plt.ylabel("Primal Variables")
plt.legend()
plt.show()
plt.figure
for i in range(0,len(u[0])):
plt.plot(iters,[u[j][i] for j in range(0,N)],label="u"+str(i))
for i in range(0,len(v[0])):
plt.plot(iters,[v[j][i] for j in range(0,N)],label="v"+str(i))
plt.xlabel("Iteration")
plt.ylabel("Dual Variables")
plt.legend()
plt.show()
plt.figure
plt.semilogy(iters,mu)
plt.xlabel("Iteration")
plt.ylabel("Barrier Penalty")
plt.show()
plt.figure
plt.semilogy(range(1,N),E)
plt.xlabel("Iteration")
plt.ylabel("Error")
plt.show()
Problem 3: Nonconvex#
f = lambda x: x[0] + 2*x[1] + x[2]**2
def c(x):
rhs = np.zeros(2)
rhs[0] = x[0] + x[1] - 1
rhs[1] = x[2] - x[1] - 3
return rhs
# Indices of variables with lower bound of zero
vb = [1,2]
x0 = np.ones(3)
u0 = np.ones(2)
v0 = np.ones(2)
## Test barrier subproblem
# xtest, vtest, utest, Etest = barrier_subproblem(x0,v0,u0,f,c,vb,1E-1,1E-10,max_iter=10,verbose=True)
##
x, v, u, mu, E = interior_point(x0,f,c,vb)
*** Barrier Subproblem 0 ***
mu = 10.0
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 4.0000e+00 3.1623e+00 9.00e+00 6.17e+00 8.25e+00 8.06e+00 0.00e+00 0.00e+00
*** Barrier Subproblem 1 ***
mu = 2.0
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 4.0013e+01 4.4202e-10 2.80e+01 3.79e+00 1.52e+00 3.22e+00 0.00e+00 0.00e+00
1 1.6354e+01 1.3323e-15 6.22e+00 9.26e-01 2.49e-01 1.35e+00 0.00e+00 0.00e+00
*** Barrier Subproblem 2 ***
mu = 0.4
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.2029e+01 5.9350e-11 2.31e+00 3.94e-01 1.93e-01 6.76e-01 0.00e+00 0.00e+00
*** Barrier Subproblem 3 ***
mu = 0.08000000000000002
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0363e+01 5.0465e-11 4.67e-01 7.00e-02 6.98e-02 1.66e-01 0.00e+00 0.00e+00
*** Barrier Subproblem 4 ***
mu = 0.016000000000000004
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0077e+01 8.9693e-12 7.01e-02 1.52e-02 5.71e-03 2.39e-02 0.00e+00 0.00e+00
*** Barrier Subproblem 5 ***
mu = 0.0020238577025077633
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0016e+01 4.8606e-13 1.42e-02 3.45e-03 7.38e-04 4.77e-03 0.00e+00 0.00e+00
*** Barrier Subproblem 6 ***
mu = 9.104790579399288e-05
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0002e+01 4.5590e-13 1.94e-03 4.78e-04 9.46e-05 6.54e-04 0.00e+00 0.00e+00
1 1.0000e+01 6.1284e-14 2.94e-07 6.48e-09 4.36e-07 8.52e-07 0.00e+00 0.00e+00
*** Barrier Subproblem 7 ***
mu = 8.687702517211205e-07
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0000e+01 4.4409e-16 9.02e-05 2.23e-05 4.29e-06 3.04e-05 0.00e+00 0.00e+00
1 1.0000e+01 2.6645e-15 1.67e-08 1.36e-11 1.69e-08 2.21e-08 0.00e+00 0.00e+00
*** Barrier Subproblem 8 ***
mu = 1e-09
Iter. f(x) ||c(x)|| E ||dx|| ||dv|| ||du|| delta_A delta_W
0 1.0000e+01 1.1102e-16 8.68e-07 2.15e-07 4.23e-08 2.92e-07 0.00e+00 0.00e+00
1 1.0000e+01 0.0000e+00 1.36e-09 1.33e-15 3.33e-10 1.60e-09 0.00e+00 0.00e+00
import matplotlib.pyplot as plt
# number of iterations
N = len(x)
# iteration
iters = range(0,N)
plt.figure
for i in range(0,len(x0)):
plt.plot(iters,[x[j][i] for j in range(0,N)],label="x"+str(i))
plt.xlabel("Iteration")
plt.ylabel("Primal Variables")
plt.legend()
plt.show()
plt.figure
for i in range(0,len(u[0])):
plt.plot(iters,[u[j][i] for j in range(0,N)],label="u"+str(i))
for i in range(0,len(v[0])):
plt.plot(iters,[v[j][i] for j in range(0,N)],label="v"+str(i))
plt.xlabel("Iteration")
plt.ylabel("Dual Variables")
plt.legend()
plt.show()
plt.figure
plt.semilogy(iters,mu)
plt.xlabel("Iteration")
plt.ylabel("Barrier Penalty")
plt.show()
plt.figure
plt.semilogy(range(1,N),E)
plt.xlabel("Iteration")
plt.ylabel("Error")
plt.show()
