{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Description**: An example of implementing a derivative-based nonlinear solver in Julia and hooking it into MathProgBase. MathProgBase connects solvers to modeling interfaces like JuMP and AMPL, which provide automatic computation of exact first and second-order derivatives.\n",
"\n",
"**Author**: Miles Lubin\n",
"\n",
"**License**:
This work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Newton's method in less than 50 lines\n",
"\n",
"In this notebook, we demonstrate how to implement [Newton's method](https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization), by querying derivatives through the [MathProgBase nonlinear interface](http://mathprogbasejl.readthedocs.org/en/latest/nlp.html). We then demonstrate using this new solver from both JuMP and AMPL."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using MathProgBase"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First we define the ``NewtonSolver`` object, which holds solver options (here there are none). Solver objects are used to instantiate a ``NewtonData`` object, which is the object which then solves the instance of the optimization problem. The ``NewtonData`` type is a subtype of the ``AbstractNonlinearModel`` type and stores all the instance data we need to run the algorithm."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"type NewtonSolver <: MathProgBase.AbstractMathProgSolver\n",
"end\n",
"\n",
"type NewtonData <: MathProgBase.AbstractNonlinearModel\n",
" numVar::Int\n",
" d # NLP evaluator\n",
" x::Vector{Float64}\n",
" hess_I::Vector{Int} # 1st component of Hessian sparsity pattern\n",
" hess_J::Vector{Int} # 2nd component of Hessian sparsity pattern\n",
" status::Symbol\n",
"end\n",
"\n",
"MathProgBase.NonlinearModel(solver::NewtonSolver) = NewtonData(0,nothing,Float64[],Int[],Int[],:Uninitialized);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This method is called when we're about to solve the optimization problem. It provides basic problem dimensions and lower and upper bound vectors, as well as the ``AbstractNLPEvaluator`` object, which is an oracle that can be used by the solver to query function and derivative evaluations."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"function MathProgBase.loadproblem!(m::NewtonData, numVar, numConstr, l, u, lb, ub, sense, d::MathProgBase.AbstractNLPEvaluator)\n",
" @assert numConstr == 0 # we don't handle constraints\n",
" @assert all(l .== -Inf) && all(u .== Inf) # or variable bounds\n",
" @assert sense == :Min # or maximization\n",
" \n",
" MathProgBase.initialize(d, [:Grad, :Hess]) # request gradient and hessian evaluations\n",
" \n",
" m.d = d\n",
" m.numVar = numVar\n",
" m.x = zeros(numVar)\n",
" m.hess_I, m.hess_J = MathProgBase.hesslag_structure(d)\n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"And here's the actual implementation of Newton's method. We assume:\n",
"\n",
"- The domain of the objective function is $\\mathbb{R}^n$\n",
"- The objective function is strictly convex, so the Hessian matrix is positive definite.\n",
"\n",
"We also fix the step size to 1 and do not perform a line search."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"function MathProgBase.optimize!(m::NewtonData)\n",
" \n",
" iteration = 0\n",
"\n",
" ∇f = Array(Float64,m.numVar)\n",
" hess_val = Array(Float64, length(m.hess_I)) # array to store nonzeros in Hessian\n",
" MathProgBase.eval_grad_f(m.d, ∇f, m.x) # writes gradient to ∇f vector\n",
" ∇f_norm = norm(∇f)\n",
" while ∇f_norm > 1e-5\n",
" println(\"Iteration $iteration. Gradient norm $∇f_norm, x = $(m.x)\");\n",
" \n",
" # Evaluate the Hessian of the Lagrangian.\n",
" # We don't have any constraints, so this is just the Hessian.\n",
" MathProgBase.eval_hesslag(m.d, hess_val, m.x, 1.0, Float64[])\n",
" \n",
" # The derivative evaluator provides the Hessian matrix\n",
" # in lower-triangular sparse (i,j,value) triplet format,\n",
" # so now we convert this into a Julia symmetric sparse matrix object.\n",
" H = Symmetric(sparse(m.hess_I,m.hess_J,hess_val),:L)\n",
" \n",
" m.x = m.x - H\\∇f # newton step\n",
" \n",
" MathProgBase.eval_grad_f(m.d, ∇f, m.x)\n",
" ∇f_norm = norm(∇f)\n",
" iteration += 1\n",
" end\n",
" \n",
" println(\"Iteration $iteration. Gradient norm $∇f_norm, x = $(m.x), DONE\");\n",
" \n",
" m.status = :Optimal \n",
"end;"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We implement a few methods to provide starting solutions and query solution status."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"MathProgBase.setwarmstart!(m::NewtonData,x) = (m.x = x)\n",
"MathProgBase.status(m::NewtonData) = m.status\n",
"MathProgBase.getsolution(m::NewtonData) = m.x\n",
"MathProgBase.getconstrduals(m::NewtonData) = [] # lagrange multipliers on constraints (none)\n",
"MathProgBase.getreducedcosts(m::NewtonData) = zeros(m.numVar) # lagrange multipliers on variable bounds\n",
"MathProgBase.getobjval(m::NewtonData) = MathProgBase.eval_f(m.d, m.x);"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using JuMP"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Iteration 0. Gradient norm 11.661903789690601, x = [0.0,0.0]\n",
"Iteration 1. Gradient norm 1.9860273225978185e-15, x = [4.999999999999999,2.9999999999999996], DONE\n"
]
}
],
"source": [
"jumpmodel = Model(solver=NewtonSolver())\n",
"@variable(jumpmodel, x)\n",
"@variable(jumpmodel, y)\n",
"@NLobjective(jumpmodel, Min, (x-5)^2 + (y-3)^2)\n",
"\n",
"status = solve(jumpmodel);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that we've converged to the optimal solution in one step. This is expected when applying Newton's method to a strictly convex quadratic objective function."
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"### Now let's hook our solver to AMPL!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"AMPL is a commercial algebraic modeling language. You can download a limited demo of AMPL [here](http://ampl.com/try-ampl/download-a-demo-version/) and extract it into the ``ampl-demo`` directory under your home directory to follow along."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"mod = \"\"\"\n",
"var x;\n",
"var y;\n",
"\n",
"minimize OBJ:\n",
" (x-5)^2 + (y-3)^2;\n",
"\n",
"write gquad;\n",
"\n",
"\"\"\"\n",
"fd = open(\"quad.mod\",\"w\")\n",
"write(fd, mod)\n",
"close(fd)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The proprietary `ampl` executable translates the `.mod` file into a `.nl` file which we can directly read from Julia."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
";~/ampl-demo/ampl quad.mod"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
";cat quad.nl | head -n 3"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The [AmplNLReader](https://github.com/dpo/AmplNLReader.jl) and [AMPLMathProgInterface](https://github.com/mlubin/AMPLMathProgInterface.jl) packages are not yet officially registered as Julia packages. They're not yet recommended for general use.\n",
"\n",
"But anyway, here's a demo of what they can do:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"using AmplNLReader, AMPLMathProgInterface"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here we read in the `.nl` file and call the ``AMPLMathProgInterface.loadamplproblem!`` method to load the `.nl` file into our solver which we just wrote. Then we just call ``optimize!`` and we're done."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"nlp = AmplNLReader.AmplModel(\"quad.nl\")\n",
"\n",
"m = MathProgBase.model(NewtonSolver())\n",
"AMPLMathProgInterface.loadamplproblem!(m, nlp)\n",
"MathProgBase.optimize!(m);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Again we converge to the optimal solution in one step, but this time we used AMPL to compute the derivatives of the model instead of JuMP."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 0.4.3",
"language": "julia",
"name": "julia-0.4"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "0.4.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}