Specifying the IPNewton method in Optim.jl

Tianhao Zhao

Introduction

Optim.jl is a powerful Julia package that offers a wide range of optimization algorithms. Among these, interior-point Newton (IPNewton) is the only gradient-based method capable of solving non-linearly constrained problems:

where $n$$n$ represents the number of control variables and $p$$p$ denotes the number of non-linear inequality constraints (with equality constraints satisfying $l_c=u_c$$l_c = u_c$).

Despite its flexibility, IPNewton can magically accelerate the algorithm by requiring users to provide the Jacobians and Hessians for both objective function $f(x)$$f(x)$ and constraint functions $c(x)$$c(x)$. Even though automatic differentiation is available today, user-provided information still has the best performance when it is available.

Notes: Automatic differentiation is available by option autodiff=true and implemented by ForwardDiff.jl

In this post, I clearly outline the ideas behind and math formulas for these data structures (primarily matrices) to help readers better understand the official documentation and tailor the method to their specific problems.

What to provide

For readability, we define the following alias:

• F64 = Float64

• Vec = AbstractVector

• Mat = AbstractMatrix

• V64 = Vector{Float64}

• M64 = Matrix{Float64}

• Num = Union{Real, Int, ...}

To solve the problem of Eq ($\text{1}$$\ref{eq:01}$) using IPNewton API, one need to manually specify the following functions:

Illustrative example

Let's illustrate how to specify these functions using a non-trivial example of economic modeling:

Notes: to avoid confusion, I use $C$$C$ to denote consumption but use $c$$c$ to denote inequality constraints.

Notes: Box constraints on the control variables are separately processed by the solver API, so there is no need to convert them into some inequality constraints of $c(x)$$c(x)$.

Notes: The IPNewton API allows custom but constant lower bounds and uppwer bounds for $c(x)$$c(x)$. This may help to simplify the definition of some $c(x)$$c(x)$ constraints.

This is a typical "expert" who produces consumption goods using housing wealth and labor. The expert also borrows in form of bond $a^{\prime }$$a'$ at rate and face collateral constraint that depends on $h^{\prime }$$h'$. This problem has many features that are common in many models. Meanwhile, we also make some technical assumptions regarding computer programming:

1. When solving this problem, the expected value function $\overline{v}(a^{\prime },h^{\prime },z)$$\bar{v}(a',h',z)$ has been interpolated. It is wrapped up as a callable function such that we do not need to explicitly write the expectation operation over $v(a,h,z)$$v(a,h,z)$.

2. The callable interpolant of $\overline{v}(\cdot )$$\bar{v}(\cdot)$ can smoothly handle any point in ${\mathbb{R}}^3$$\mathbb{R}^3$. If a point is not in the feasible state space, then $\overline{v}(\cdot )$$\bar{v}(\cdot)$ can still return a number (e.g. a very negative number).

Reformulating the problem to a minimization one and cancelling out $c$$c$ using the equality constraint, we re-write the problem as:

where $M$$M$ is the information set we have for this problem (including current states, parameters, computational space, and other model details); $ϵ$$\epsilon$ is a small amount to avoid exact zero; the substituted consumption $C(x)$$C(x)$ is:

and the objective function is:

Explain by part

The official documentation of Optim.jl does not explicitly explain the pipeline of IPNnewton. Thus, it is a little bit less intuitive to understand the behavior of required functions. In short, a pipeline is like:

1. User provides a pure objective function $f(x)$$f(x)$ and a pure constraint function $c(x)$$c(x)$

2. The algorithm allocates Jacobians and Hessians for $f(x)$$f(x)$ and $c(x)$$c(x)$ internally

3. Users provide "pipe" functions that receive and modify these internally created Jacobian and Hessian objects in-place without manually managing their construction, use, and destruction.

Objective function

The objective function $f(x)$$f(x)$ is set up as a scalar function call fun(x::Vec)::Num that receives a $n$$n$-vector and returns a scalar of the function value. This is standard and required even if the autodiff is turned on.

Hint: If your objective function requires extra parameters, then consider wrapping it with the information set $M$$M$ as a closure.

In the minimization API, the objective function is the negative of Eq ($\text{18}$$\ref{eq:020}$).

Jacobian of objective function

The Jacobian $\mathrm{\nabla }f(x)$$\nabla f(x)$ is a column vector:

In the illustrative example:

The function fun_grad!(g,x) should modify each element of vector g using the values in the above equation at exactly the same position. This function only updates g without returning it.

Note: One can quickly notice that: evaluating $\mathrm{\nabla }f(x)$$\nabla f(x)$ requires computing the partial derivatives of the expected value function $\overline{v}$$\bar{v}$. This can be done explicitly (e.g. if you are using Chebyshev polynomials) or using finite difference approximation (most cases), or using automatic differentiation on the callable function $\overline{v}(\cdot )$$\bar{v}(\cdot)$. Readers may have to write the difference program themselves and properly handle the case on the boundary of feasible domains. This idea also applies to all the following matrices.

Hessian of objective function

The Hessian ${\mathrm{\nabla }}^{2}f(x)$$\nabla^2 f(x)$ is an $n\times n$$n\times n$ matrix:

In the illustrative example, this matrix looks like:

The function fun_hess!(h::Mat,x::Vec) should modify each element of matrix h using the values in the above equation at exactly the same position. This function updates h without returning it.

Non-linear constraints

The non-linear constraints $c(x)$$c(x)$ is programmed as a vector-value function. The function con_c!(c::Vec,x::Vec)::Vec should modify the passed c vector and also return c. (Wierd behavior uhmmmm)

Jacobian of the non-linear constraints

Since $\mathrm{\nabla }c(x)$$\nabla c(x)$ is a vector-value function, its Jacobian is now a matrix looking like:

The function con_jacobian!(J::Mat,x::Vec)::Mat, however, now not only modifies the passed J matrix but returns J as well. Retuning a copy is the biggest difference of constraint-related API cp. objective-related API.

In the illustrative example, such a Jacobian looks like:

Hessian of the non-linear constraints

The IPNewton algorithm enforces the inequality constraints using log-barrier functions. It optimizes an augmented Lagrangian with barrier terms. Thus, the total Hessian includes:

which adds extra $p$$p$ Hessian of the constraints on the top of the orginal Hessian of the objective function. This formula is the key to understand the behavior of function con_h!(h::Mat,x::Vec,lam::Vec)::Nothing.

• h::Mat is the Hessian of the objective function $f(x)$$f(x)$ rather than the Hessian of the constraints $c(x)$$c(x)$. This object is internally created by the IPNewton API and we users don't need to manage its construction and destruction.

• lam::Vec is a $p$$p$-element vector of the Lagangian multipliers ${\lambda }_i$$\lambda_i$. This object is also internally created by the API and we just take it as given

The function con_h!() looks like:

function con_h!(h,x,lam)
  # add the hessian of the 1st constraint on the top of `h`; multiplied by the 1st multiplier `lam`
  h += lam[1] * hess_cons_1(x)  # hess_cons_1(x)::Matrix, n*n
  
  # add the hessian of the 2nd constraint on the top of `h`
  h += lam[2] * hess_cons_2(x)
  
  # ...
  
  return nothing
end

In the illustrative example, each constraint's Hessian wrt the control $x$$x$ are:

which are derived from the 1st and 2nd row of the Jacobian in Eq ($\text{24}$$\ref{eq:024}$) respectively.

Hint: Let's be more specific, step by step. The first row of Eq ($\text{24}$$\ref{eq:024}$) is the Jacobian of the 1st constraint $c_1(x):=C(x)-ϵ$$c_1(x):=C(x) -\epsilon$ :

$$\begin{array}{}\text{(28)}& \mathrm{\nabla }{c}_1(x)=\left[\begin{array}{ccc}1& -p& z(1-\alpha )h^{\alpha }{\ell }^{-\alpha }\end{array}\right]\end{array}$$

A typical Jacobian is a column vector so we transpose the above matrix row to:

$$\begin{array}{}\text{(29)}& \mathrm{\nabla }{c}_1(x)=\left[\begin{array}{c}1\\ -p\\ z(1-\alpha )h^{\alpha }{\ell }^{-\alpha }\end{array}\right]=\left[\begin{array}{c}\frac{\partial c_1(x)}{\partial a^{\prime }}\\ \frac{\partial c_1(x)}{\partial h^{\prime }}\\ \frac{\partial c_1(x)}{\partial \ell }\end{array}\right]\end{array}$$

Then, we take derivatives again to get the Hessian ${\mathrm{\nabla }}^2{c}_1(x)$$\nabla^2 c_1(x)$ of the 1st constraint as Eq ($\text{26}$$\ref{eq:026}$)

Appendix: Autodiff of $\overline{v}$$\bar{v}$

The evaluation of Jacobians and Hessians requires approximating the Jacobian and Hessian of the expected value function $\overline{v}(a^{\prime },h^{\prime },z)$$\bar{v}(a',h',z)$.

Suppose we have an guess of $\overline{v}(\cdot )$$\bar{v}(\cdot)$. Suppose we have interpolated the guess and wrapped it as a callable object, which means that it can be called like a function that returns the value of $\overline{v}(\cdot )$$\bar{v}(\cdot)$ at a given point $x^{\prime }$$x'$:

xp = [ap,hp,z] # a', h', z
val = ev(xp)

The ForwardDiff.jl package can perform automatic differentiation on such a callable object. For example, evaluting Jacobian and Hessian at point $[1.0,0.9,0.8]$$[1.0,0.9,0.8]$:

import ForwardDiff as fd

x0 = [1.0,0.9,0.8]

# Jacobian: 3-vector
fd.gradient(ev, x0)::Vector{Float64}

# Hessian: 3*3 matrix
fd.hessian(ev, x0)::Matrix{Float64}

Users may check the GitHub page of the package for more details.