{"id":6589,"date":"2019-09-02T23:32:44","date_gmt":"2019-09-02T22:32:44","guid":{"rendered":"http:\/\/www.walkingrandomly.com\/?p=6589"},"modified":"2019-09-02T23:32:44","modified_gmt":"2019-09-02T22:32:44","slug":"second-order-cone-programming-socp-using-the-nag-library-for-python","status":"publish","type":"post","link":"https:\/\/walkingrandomly.com\/?p=6589","title":{"rendered":"Second Order Cone Programming (SOCP) using the NAG Library for Python"},"content":{"rendered":"<p><script src=\"https:\/\/cdnjs.cloudflare.com\/ajax\/libs\/mathjax\/2.7.5\/latest.js?config=TeX-AMS_HTML\"><\/script><br \/>\n<script type=\"text\/x-mathjax-config\">\n    MathJax.Hub.Config({\n        tex2jax: {\n            inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n            displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ],\n            processEscapes: true,\n            processEnvironments: true\n        },\n        \/\/ Center justify equations in code and markdown cells. Elsewhere\n        \/\/ we use CSS to left justify single line equations in code cells.\n        displayAlign: 'center',\n        \"HTML-CSS\": {\n            styles: {'.MathJax_Display': {\"margin\": 0}},\n            linebreaks: { automatic: true }\n        }\n    });\n    <\/script><\/p>\n<h2>What is Second Order Cone Programming (SOCP)?<\/h2>\n<p>&nbsp;<\/p>\n<p><a href=\"https:\/\/en.wikipedia.org\/wiki\/Second-order_cone_programming\">Second Order Cone Programming (SOCP) problems<\/a> are a type of optimisation problem that have applications in many areas of science, finance and engineering. A summary of the type of problems that can make use of SOCP, including things as diverse as designing antenna arrays, finite impulse response (FIR) filters and structural equilibrium problems can be found in the paper <a href=\"http:\/\/www.seas.ucla.edu\/~vandenbe\/publications\/socp.pdf\">&#8216;Applications of Second Order Cone Programming&#8217; by Lobo et al<\/a>. There are also a couple of examples of using SOCP for portfolio optimisation in the <a href=\"https:\/\/github.com\/numericalalgorithmsgroup\/NAGPythonExamples\/tree\/master\/local_optimisation\/SOCP\">GitHub repository of the Numerical Algorithms Group (NAG)<\/a>.<\/p>\n<p>A large scale SOCP solver was one of the highlights of the Mark 27 release of the NAG library (<a href=\"https:\/\/www.nag.co.uk\/market\/posters\/socp.pdf\">See here for a poster about its performance<\/a>). Those who have used the NAG library for years will expect this solver to have interfaces <a href=\"https:\/\/www.nag.com\/numeric\/nl\/nagdoc_latest\/flhtml\/e04\/e04ptf.html\">in Fortran<\/a> <a href=\"https:\/\/www.nag.com\/numeric\/nl\/nagdoc_latest\/clhtml\/e04\/e04ptc.html\">and C<\/a>\u00a0and, of course, they are there.\u00a0In addition to this is the fact\u00a0that Mark 27 of the\u00a0<a href=\"https:\/\/www.nag.co.uk\/nag-library-python\">NAG Library for <em>Python<\/em><\/a>\u00a0was released at the same time as the Fortran and C interfaces which reflects the importance of Python in today&#8217;s numerical computing landscape.<\/p>\n<p>Here&#8217;s a quick demo of how the new SOCP solver works in Python. The code is based on a notebook in NAG&#8217;s <a href=\"https:\/\/github.com\/numericalalgorithmsgroup\/NAGPythonExamples\">PythonExamples GitHub repository<\/a>.<\/p>\n<p>NAG&#8217;s <a href=\"https:\/\/www.nag.co.uk\/numeric\/py\/nagdoc_latest\/naginterfaces.library.opt.html#naginterfaces.library.opt.handle_solve_socp_ipm\"><strong>handle_solve_socp_ipm<\/strong><\/a> function (also known as e04pt) is a solver from the NAG optimization modelling suite for large-scale second-order cone programming (SOCP) problems based on an interior point method (IPM).<\/p>\n<p>$$<br \/>\n\\begin{array}{ll}<br \/>\n{\\underset{x \\in \\mathbb{R}^{n}}{minimize}\\ } &amp; {c^{T}x} \\\\<br \/>\n\\text{subject to} &amp; {l_{A} \\leq Ax \\leq u_{A}\\text{,}} \\\\<br \/>\n&amp; {l_{x} \\leq x \\leq u_{x}\\text{,}} \\\\<br \/>\n&amp; {x \\in \\mathcal{K}\\text{,}} \\\\<br \/>\n\\end{array}<br \/>\n$$<\/p>\n<p>where $\\mathcal{K} = \\mathcal{K}^{n_{1}} \\times \\cdots \\times \\mathcal{K}^{n_{r}} \\times \\mathbb{R}^{n_{l}}$ is a Cartesian product of quadratic (second-order type) cones and $n_{l}$-dimensional real space, and $n = \\sum_{i = 1}^{r}n_{i} + n_{l}$ is the number of decision variables. Here $c$, $x$, $l_x$ and $u_x$ are $n$-dimensional vectors.<\/p>\n<p>$A$ is an $m$ by $n$ sparse matrix, and $l_A$ and $u_A$ and are $m$-dimensional vectors. Note that $x \\in \\mathcal{K}$ partitions subsets of variables into quadratic cones and each $\\mathcal{K}^{n_{i}}$ can be either a quadratic cone or a rotated quadratic cone. These are defined as follows:<\/p>\n<ul>\n<li>Quadratic cone:<\/li>\n<\/ul>\n<p>$$<br \/>\n\\mathcal{K}_{q}^{n_{i}} := \\left\\{ {z = \\left\\{ {z_{1},z_{2},\\ldots,z_{n_{i}}} \\right\\} \\in {\\mathbb{R}}^{n_{i}} \\quad\\quad : \\quad\\quad z_{1}^{2} \\geq \\sum\\limits_{j = 2}^{n_{i}}z_{j}^{2},\\quad\\quad\\quad z_{1} \\geq 0} \\right\\}\\text{.}<br \/>\n$$<\/p>\n<ul>\n<li>Rotated quadratic cone:<\/li>\n<\/ul>\n<p>$$<br \/>\n\\mathcal{K}_{r}^{n_{i}} := \\left\\{ {z = \\left\\{ {z_{1},z_{2},\\ldots,z_{n_{i}}} \\right\\} \\in {\\mathbb{R}}^{n_{i}}\\quad\\quad:\\quad \\quad\\quad 2z_{1}z_{2} \\geq \\sum\\limits_{j = 3}^{n_{i}}z_{j}^{2}, \\quad\\quad z_{1} \\geq 0, \\quad\\quad z_{2} \\geq 0} \\right\\}\\text{.}<br \/>\n$$<\/p>\n<p>For a full explanation of this routine, refer to <a href=\"https:\/\/www.nag.co.uk\/numeric\/nl\/nagdoc_27\/clhtml\/e04\/e04ptc.html\">e04ptc in the NAG Library Manual<\/a><\/p>\n<h2>Using the NAG SOCP Solver from Python<\/h2>\n<p>&nbsp;<\/p>\n<p>This example, derived from the documentation for the <a href=\"https:\/\/www.nag.co.uk\/numeric\/nl\/nagdoc_27\/clhtml\/e04\/e04rbc.html#example\"><strong>handle_set_group<\/strong><\/a> function solves the following SOCP problem<\/p>\n<p>minimize $${10.0x_{1} + 20.0x_{2} + x_{3}}$$<\/p>\n<pre><span class=\"kn\">from<\/span> <span class=\"nn\">naginterfaces.base<\/span> <span class=\"k\">import<\/span> <span class=\"n\">utils<\/span>\r\n<span class=\"kn\">from<\/span> <span class=\"nn\">naginterfaces.library<\/span> <span class=\"k\">import<\/span> <span class=\"n\">opt<\/span>\r\n\r\n<span class=\"c1\"># The problem size:<\/span>\r\n<span class=\"n\">n<\/span> <span class=\"o\">=<\/span> <span class=\"mi\">3<\/span>\r\n\r\n<span class=\"c1\"># Create the problem handle:<\/span>\r\n<span class=\"n\">handle<\/span> <span class=\"o\">=<\/span> <span class=\"n\">opt<\/span><span class=\"o\">.<\/span><span class=\"n\">handle_init<\/span><span class=\"p\">(<\/span><span class=\"n\">nvar<\/span><span class=\"o\">=<\/span><span class=\"n\">n<\/span><span class=\"p\">)<\/span>\r\n\r\n<span class=\"c1\"># Set objective function<\/span>\r\n<span class=\"n\">opt<\/span><span class=\"o\">.<\/span><span class=\"n\">handle_set_linobj<\/span><span class=\"p\">(<\/span><span class=\"n\">handle<\/span><span class=\"p\">,<\/span> <span class=\"n\">cvec<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span><span class=\"mf\">10.0<\/span><span class=\"p\">,<\/span> <span class=\"mf\">20.0<\/span><span class=\"p\">,<\/span> <span class=\"mf\">1.0<\/span><span class=\"p\">])<\/span>\r\n<\/pre>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\"><\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>subject to the bounds<br \/>\n$$<br \/>\n\\begin{array}{rllll}<br \/>\n{- 2.0} &amp; \\leq &amp; x_{1} &amp; \\leq &amp; 2.0 \\\\<br \/>\n{- 2.0} &amp; \\leq &amp; x_{2} &amp; \\leq &amp; 2.0 \\\\<br \/>\n\\end{array}<br \/>\n$$<\/p>\n<\/div>\n<\/div>\n<\/div>\n<pre><span class=\"c1\"># Set box constraints<\/span>\r\n<span class=\"n\">opt<\/span><span class=\"o\">.<\/span><span class=\"n\">handle_set_simplebounds<\/span><span class=\"p\">(<\/span>\r\n    <span class=\"n\">handle<\/span><span class=\"p\">,<\/span>\r\n    <span class=\"n\">bl<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span><span class=\"o\">-<\/span><span class=\"mf\">2.0<\/span><span class=\"p\">,<\/span> <span class=\"o\">-<\/span><span class=\"mf\">2.0<\/span><span class=\"p\">,<\/span> <span class=\"o\">-<\/span><span class=\"mf\">1.e20<\/span><span class=\"p\">],<\/span>\r\n    <span class=\"n\">bu<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span><span class=\"mf\">2.0<\/span><span class=\"p\">,<\/span> <span class=\"mf\">2.0<\/span><span class=\"p\">,<\/span> <span class=\"mf\">1.e20<\/span><span class=\"p\">]<\/span>\r\n<span class=\"p\">)<\/span>\r\n<\/pre>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\"><\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>the general linear constraints<\/p>\n<\/div>\n<\/div>\n<\/div>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\"><\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\\begin{array}{lcrcrcrclcl}<br \/>\n&amp; &amp; {- 0.1x_{1}} &amp; &#8211; &amp; {0.1x_{2}} &amp; + &amp; x_{3} &amp; \\leq &amp; 1.5 &amp; &amp; \\\\<br \/>\n1.0 &amp; \\leq &amp; {- 0.06x_{1}} &amp; + &amp; x_{2} &amp; + &amp; x_{3} &amp; &amp; &amp; &amp; \\\\<br \/>\n\\end{array}<\/div>\n<\/div>\n<\/div>\n<pre><span class=\"c1\"># Set linear constraints<\/span>\r\n<span class=\"n\">opt<\/span><span class=\"o\">.<\/span><span class=\"n\">handle_set_linconstr<\/span><span class=\"p\">(<\/span>\r\n    <span class=\"n\">handle<\/span><span class=\"p\">,<\/span>\r\n    <span class=\"n\">bl<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span><span class=\"o\">-<\/span><span class=\"mf\">1.e20<\/span><span class=\"p\">,<\/span> <span class=\"mf\">1.0<\/span><span class=\"p\">],<\/span>\r\n    <span class=\"n\">bu<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span><span class=\"mf\">1.5<\/span><span class=\"p\">,<\/span> <span class=\"mf\">1.e20<\/span><span class=\"p\">],<\/span>\r\n    <span class=\"n\">irowb<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span> <span class=\"mi\">1<\/span><span class=\"p\">,<\/span> <span class=\"mi\">1<\/span><span class=\"p\">,<\/span> <span class=\"mi\">2<\/span><span class=\"p\">,<\/span> <span class=\"mi\">2<\/span><span class=\"p\">,<\/span> <span class=\"mi\">2<\/span><span class=\"p\">],<\/span>\r\n    <span class=\"n\">icolb<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span> <span class=\"mi\">2<\/span><span class=\"p\">,<\/span> <span class=\"mi\">3<\/span><span class=\"p\">,<\/span> <span class=\"mi\">1<\/span><span class=\"p\">,<\/span> <span class=\"mi\">2<\/span><span class=\"p\">,<\/span> <span class=\"mi\">3<\/span><span class=\"p\">],<\/span>\r\n    <span class=\"n\">b<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span><span class=\"o\">-<\/span><span class=\"mf\">0.1<\/span><span class=\"p\">,<\/span> <span class=\"o\">-<\/span><span class=\"mf\">0.1<\/span><span class=\"p\">,<\/span> <span class=\"mf\">1.0<\/span><span class=\"p\">,<\/span> <span class=\"o\">-<\/span><span class=\"mf\">0.06<\/span><span class=\"p\">,<\/span> <span class=\"mf\">1.0<\/span><span class=\"p\">,<\/span> <span class=\"mf\">1.0<\/span><span class=\"p\">]<\/span>\r\n    <span class=\"p\">);<\/span>\r\n<\/pre>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\"><\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>and the cone constraint<\/p>\n<p>$$\\left( {x_{3},x_{1},x_{2}} \\right) \\in \\mathcal{K}_{q}^{3}\\text{.}$$<\/p>\n<\/div>\n<\/div>\n<\/div>\n<pre><span class=\"c1\"># Set cone constraint<\/span>\r\n<span class=\"n\">opt<\/span><span class=\"o\">.<\/span><span class=\"n\">handle_set_group<\/span><span class=\"p\">(<\/span>\r\n    <span class=\"n\">handle<\/span><span class=\"p\">,<\/span>\r\n    <span class=\"n\">gtype<\/span><span class=\"o\">=<\/span><span class=\"s1\">'Q'<\/span><span class=\"p\">,<\/span>\r\n    <span class=\"n\">group<\/span><span class=\"o\">=<\/span><span class=\"p\">[<\/span> <span class=\"mi\">3<\/span><span class=\"p\">,<\/span><span class=\"mi\">1<\/span><span class=\"p\">,<\/span> <span class=\"mi\">2<\/span><span class=\"p\">],<\/span>\r\n    <span class=\"n\">idgroup<\/span><span class=\"o\">=<\/span><span class=\"mi\">0<\/span>\r\n<span class=\"p\">);<\/span>\r\n<\/pre>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\"><\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>We set some algorithmic options. For more details on the options available, <a href=\"https:\/\/www.nag.com\/numeric\/nl\/nagdoc_27\/clhtml\/e04\/e04ptc.html\">refer to the routine documentation<\/a><\/p>\n<\/div>\n<\/div>\n<\/div>\n<pre><span class=\"c1\"># Set some algorithmic options.<\/span>\r\n<span class=\"k\">for<\/span> <span class=\"n\">option<\/span> <span class=\"ow\">in<\/span> <span class=\"p\">[<\/span>\r\n        <span class=\"s1\">'Print Options = NO'<\/span><span class=\"p\">,<\/span>\r\n        <span class=\"s1\">'Print Level = 1'<\/span>\r\n<span class=\"p\">]:<\/span>\r\n    <span class=\"n\">opt<\/span><span class=\"o\">.<\/span><span class=\"n\">handle_opt_set<\/span><span class=\"p\">(<\/span><span class=\"n\">handle<\/span><span class=\"p\">,<\/span> <span class=\"n\">option<\/span><span class=\"p\">)<\/span>\r\n\r\n<span class=\"c1\"># Use an explicit I\/O manager for abbreviated iteration output:<\/span>\r\n<span class=\"n\">iom<\/span> <span class=\"o\">=<\/span> <span class=\"n\">utils<\/span><span class=\"o\">.<\/span><span class=\"n\">FileObjManager<\/span><span class=\"p\">(<\/span><span class=\"n\">locus_in_output<\/span><span class=\"o\">=<\/span><span class=\"kc\">False<\/span><span class=\"p\">)<\/span>\r\n<\/pre>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\"><\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>Finally, we call the solver<\/p>\n<\/div>\n<\/div>\n<\/div>\n<pre><span class=\"c1\"># Call SOCP interior point solver<\/span>\r\n<span class=\"n\">result<\/span> <span class=\"o\">=<\/span> <span class=\"n\">opt<\/span><span class=\"o\">.<\/span><span class=\"n\">handle_solve_socp_ipm<\/span><span class=\"p\">(<\/span><span class=\"n\">handle<\/span><span class=\"p\">,<\/span> <span class=\"n\">io_manager<\/span><span class=\"o\">=<\/span><span class=\"n\">iom<\/span><span class=\"p\">)<\/span>\r\n<\/pre>\n<pre> ------------------------------------------------\r\n  E04PT, Interior point method for SOCP problems\r\n ------------------------------------------------\r\n\r\n Status: converged, an optimal solution found\r\n Final primal objective value -1.951817E+01\r\n Final dual objective value   -1.951817E+01\r\n<\/pre>\n<div class=\"cell border-box-sizing text_cell rendered\">\n<div class=\"prompt input_prompt\"><\/div>\n<div class=\"inner_cell\">\n<div class=\"text_cell_render border-box-sizing rendered_html\">\n<p>The optimal solution is<\/p>\n<\/div>\n<\/div>\n<\/div>\n<pre><span class=\"n\">result<\/span><span class=\"o\">.<\/span><span class=\"n\">x<\/span>\r\n<\/pre>\n<pre>array([-1.26819151, -0.4084294 ,  1.3323379 ])<\/pre>\n<p>and the objective function value is<\/p>\n<pre><span class=\"n\">result<\/span><span class=\"o\">.<\/span><span class=\"n\">rinfo<\/span><span class=\"p\">[<\/span><span class=\"mi\">0<\/span><span class=\"p\">]<\/span>\r\n<\/pre>\n<pre>-19.51816515094211<\/pre>\n<p>Finally, we clean up after ourselves by destroying the handle<\/p>\n<pre><span class=\"c1\"># Destroy the handle:<\/span>\r\n<span class=\"n\">opt<\/span><span class=\"o\">.<\/span><span class=\"n\">handle_free<\/span><span class=\"p\">(<\/span><span class=\"n\">handle<\/span><span class=\"p\">)<\/span>\r\n<\/pre>\n<p>As you can see, the way to use the NAG Library for Python interface follows the mathematics quite closely.<br \/>\nNAG also recently added support for the popular <a href=\"https:\/\/www.cvxpy.org\">cvxpy modelling language<\/a> that I&#8217;ll discuss another time.<\/p>\n<h2>Links<\/h2>\n<ul>\n<li><a href=\"https:\/\/github.com\/numericalalgorithmsgroup\/NAGPythonExamples\/tree\/master\/local_optimisation\/SOCP\">SOCP examples on NAG&#8217;s GitHub page<\/a><\/li>\n<li><a href=\"https:\/\/www.nag.com\/market\/posters\/socp.pdf\">Poster discussing performance of the NAG SOCP solver<\/a><\/li>\n<\/ul>\n","protected":false},"excerpt":{"rendered":"<p>What is Second Order Cone Programming (SOCP)? &nbsp; Second Order Cone Programming (SOCP) problems are a type of optimisation problem that have applications in many areas of science, finance and engineering. A summary of the type of problems that can make use of SOCP, including things as diverse as designing antenna arrays, finite impulse response [&hellip;]<\/p>\n","protected":false},"author":3,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"jetpack_post_was_ever_published":false,"footnotes":"","jetpack_publicize_message":"","jetpack_publicize_feature_enabled":true,"jetpack_social_post_already_shared":true,"jetpack_social_options":{"image_generator_settings":{"template":"highway","default_image_id":0,"font":"","enabled":false},"version":2}},"categories":[4,28,31,42],"tags":[],"class_list":["post-6589","post","type-post","status-publish","format-standard","hentry","category-math-software","category-nag-library","category-python","category-tutorials"],"jetpack_publicize_connections":[],"jetpack_featured_media_url":"","jetpack_shortlink":"https:\/\/wp.me\/p3swhs-1Ih","jetpack_sharing_enabled":true,"_links":{"self":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/6589","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/users\/3"}],"replies":[{"embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=6589"}],"version-history":[{"count":31,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/6589\/revisions"}],"predecessor-version":[{"id":6622,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=\/wp\/v2\/posts\/6589\/revisions\/6622"}],"wp:attachment":[{"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=6589"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Fcategories&post=6589"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/walkingrandomly.com\/index.php?rest_route=%2Fwp%2Fv2%2Ftags&post=6589"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}