{ "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**: \"Creative
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 }