{
"cells": [
{
"cell_type": "code",
"execution_count": 43,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:02:58.731000-04:00",
"start_time": "2020-04-16T18:02:58.659Z"
}
},
"outputs": [],
"source": [
"using Random, LinearAlgebra, Statistics, StatsPlots"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:03:41.063000-04:00",
"start_time": "2020-04-16T18:03:41.050Z"
}
},
"outputs": [],
"source": [
"# JuMP is a front-end for optimization. The rest are backends.\n",
"using JuMP, SCS, GLPK, Mosek, MosekTools"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:03:41.618000-04:00",
"start_time": "2020-04-16T18:03:41.610Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"MathOptInterface.OptimizerWithAttributes(Mosek.Optimizer, Pair{MathOptInterface.AbstractOptimizerAttribute,Any}[MathOptInterface.RawParameter(\"QUIET\") => true])"
]
},
"execution_count": 45,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# These simplify choice of backends with the parameters I want.\n",
"opt_scs = optimizer_with_attributes(SCS.Optimizer, \"verbose\"=>0)\n",
"opt_glpk = optimizer_with_attributes(GLPK.Optimizer, \"msg_lev\" => GLPK.OFF)\n",
"opt_mosek = optimizer_with_attributes(Mosek.Optimizer, \"QUIET\" => true)"
]
},
{
"cell_type": "code",
"execution_count": 46,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:03:45.893000-04:00",
"start_time": "2020-04-16T18:03:45.848Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"6.0"
]
},
"execution_count": 46,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"onenorm(x) = norm(x,1)\n",
"onenorm([1 2 3])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Incoherence"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:04:01.346000-04:00",
"start_time": "2020-04-16T18:04:01.324Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"5×20 Array{Float64,2}:\n",
" 0.679107 0.297336 -0.688907 … 0.556741 0.441573 1.3378 \n",
" 0.828413 0.0649475 -0.762804 -0.568854 1.67689 1.21963 \n",
" -0.353007 -0.109017 0.397482 -0.47694 0.238481 -1.22187 \n",
" -0.134854 -0.51421 0.81163 0.174594 2.31379 0.256914\n",
" 0.586617 1.57433 -0.346355 -1.51761 -0.081851 0.799711"
]
},
"execution_count": 47,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d = 5\n",
"n = 20\n",
"Random.seed!(0)\n",
"A = randn(d,n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Compute Incoherence"
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:04:30.980000-04:00",
"start_time": "2020-04-16T18:04:30.957Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"20×20 Array{Float64,2}:\n",
" 2.22045e-16 0.596628 -0.860283 … 0.335478 0.930847 \n",
" 0.596628 0.0 -0.52904 -0.2255 0.439922 \n",
" -0.860283 -0.52904 -2.22045e-16 0.101946 -0.728318 \n",
" -0.159119 -0.113124 0.405184 0.173016 0.180499 \n",
" -0.473822 0.128045 0.439549 -0.263275 -0.670492 \n",
" -0.672912 -0.813001 0.620785 … 0.107993 -0.587188 \n",
" 0.684779 0.546908 -0.618257 0.139142 0.523754 \n",
" 0.425619 0.363483 -0.0511253 0.401633 0.654449 \n",
" -0.420792 -0.784395 0.165879 -0.119628 -0.337706 \n",
" 0.734215 0.413084 -0.579377 0.438772 0.58571 \n",
" -0.283173 -0.0166007 0.0309163 … -0.678922 -0.11669 \n",
" -0.342546 -0.624092 0.140511 -0.230608 -0.0875969\n",
" -0.622748 -0.467635 0.16071 -0.747335 -0.596043 \n",
" -0.55152 0.0298985 0.740323 -0.110472 -0.37182 \n",
" 0.313471 0.643136 -0.214969 0.0101717 0.101295 \n",
" -0.292367 -0.907657 0.155125 … 0.237937 -0.190913 \n",
" -0.564701 -0.783707 0.643148 0.423117 -0.558533 \n",
" -0.366999 -0.762058 0.2093 -0.0565986 -0.12792 \n",
" 0.335478 -0.2255 0.101946 -1.11022e-16 0.423089 \n",
" 0.930847 0.439922 -0.728318 0.423089 0.0 "
]
},
"execution_count": 48,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"normed_A = similar(A)\n",
"for i in 1:n\n",
" normed_A[:,i] = A[:,i] / norm(A[:,i])\n",
"end\n",
"Gram = normed_A' * normed_A - I"
]
},
{
"cell_type": "code",
"execution_count": 49,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:04:40.161000-04:00",
"start_time": "2020-04-16T18:04:40.151Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.9679382785673595"
]
},
"execution_count": 49,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"μ = maximum(abs.(Gram))"
]
},
{
"cell_type": "code",
"execution_count": 50,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:05:06.511000-04:00",
"start_time": "2020-04-16T18:05:06.486Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"incoherence (generic function with 1 method)"
]
},
"execution_count": 50,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function incoherence(A)\n",
" normed_A = similar(A)\n",
" for i in 1:size(A,2)\n",
" normed_A[:,i] = A[:,i] / norm(A[:,i])\n",
" end\n",
" Gram = normed_A' * normed_A - I \n",
" return maximum(abs.(Gram))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:05:31.522000-04:00",
"start_time": "2020-04-16T18:05:31.476Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.4155190340602961"
]
},
"execution_count": 57,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = randn(5*25,5*100)\n",
"incoherence(A)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Basis Pursuit by JuMP"
]
},
{
"cell_type": "code",
"execution_count": 58,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:06:28.734000-04:00",
"start_time": "2020-04-16T18:06:28.711Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"10-element Array{Float64,1}:\n",
" -0.11248681351342672\n",
" -0.5662892741372261 \n",
" -0.14948376404232092\n",
" -0.7541284390976164 \n",
" -0.35188939912748723\n",
" 1.7002809122995806 \n",
" -1.9143586614424437 \n",
" 0.310830838565587 \n",
" 0.3351183436255904 \n",
" 0.5234120006014658 "
]
},
"execution_count": 58,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d = 10\n",
"n = 20\n",
"A = randn(d,n)\n",
"x0 = zeros(n)\n",
"x0[1] = 1.0\n",
"b = A*x0"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:06:40.249000-04:00",
"start_time": "2020-04-16T18:06:40.231Z"
}
},
"outputs": [
{
"data": {
"text/latex": [
"$$ \\begin{alignat*}{1}\\text{feasibility}\\\\\n",
"\\text{Subject to} \\quad\\end{alignat*}\n",
" $$"
],
"text/plain": [
"A JuMP Model\n",
"Feasibility problem with:\n",
"Variables: 0\n",
"Model mode: AUTOMATIC\n",
"CachingOptimizer state: EMPTY_OPTIMIZER\n",
"Solver name: SCS"
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"model = Model(SCS.Optimizer)"
]
},
{
"cell_type": "code",
"execution_count": 60,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:06:47.505000-04:00",
"start_time": "2020-04-16T18:06:47.328Z"
},
"collapsed": true
},
"outputs": [
{
"data": {
"text/plain": [
"20-element Array{VariableRef,1}:\n",
" x[1] \n",
" x[2] \n",
" x[3] \n",
" x[4] \n",
" x[5] \n",
" x[6] \n",
" x[7] \n",
" x[8] \n",
" x[9] \n",
" x[10]\n",
" x[11]\n",
" x[12]\n",
" x[13]\n",
" x[14]\n",
" x[15]\n",
" x[16]\n",
" x[17]\n",
" x[18]\n",
" x[19]\n",
" x[20]"
]
},
"execution_count": 60,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@variable(model, x[1:n])"
]
},
{
"cell_type": "code",
"execution_count": 61,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:06:59.706000-04:00",
"start_time": "2020-04-16T18:06:59.690Z"
},
"collapsed": true
},
"outputs": [
{
"data": {
"text/plain": [
"10-element Array{ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.EqualTo{Float64}},ScalarShape},1}:\n",
" -0.11248681351342672 x[1] - 1.2317225437811226 x[2] + 0.60532364592491 x[3] + 0.6124028737446184 x[4] + 0.6528565097491319 x[5] + 1.6979580898265894 x[6] - 0.6862486265217821 x[7] + 0.15289903450089515 x[8] - 0.004772539509995808 x[9] + 0.9716086943526653 x[10] - 0.9042366324633847 x[11] - 0.6269118217601528 x[12] + 0.8970043567906548 x[13] - 0.4667104302119095 x[14] + 0.3817729672243532 x[15] - 0.19991603062192095 x[16] + 0.6758885311998514 x[17] + 0.44157726600978825 x[18] + 0.9334671338170232 x[19] - 0.523082339678034 x[20] = -0.11248681351342672 \n",
" -0.5662892741372261 x[1] + 0.5674902277653299 x[2] + 1.1163665397172606 x[3] + 0.4567725810041199 x[4] + 1.357259294102547 x[5] + 1.1523743380963307 x[6] - 0.11599185003107795 x[7] + 3.1476344298475682 x[8] + 1.1042089249028642 x[9] - 0.16524890818251917 x[10] - 0.9086133913056114 x[11] - 0.22708133243023448 x[12] - 1.119918619955327 x[13] - 0.6035039380546332 x[14] + 0.26272948667516594 x[15] - 1.0513935582164728 x[16] + 0.44570087758409865 x[17] + 0.5956957857056084 x[18] + 0.5821841889434244 x[19] - 0.19889138315438443 x[20] = -0.5662892741372261 \n",
" -0.14948376404232092 x[1] + 1.5517163112180166 x[2] - 1.7256775636583113 x[3] - 1.4001529409440712 x[4] - 0.40816543905259284 x[5] + 1.2256635223911487 x[6] - 0.9095984379305072 x[7] + 0.8158670916316041 x[8] - 0.33253786352787884 x[9] - 1.6199664190518903 x[10] + 2.2825810094694656 x[11] + 0.337284707387779 x[12] + 0.7416188228601653 x[13] - 0.6498398772782128 x[14] + 1.675344752615471 x[15] - 0.9999991635264687 x[16] + 0.2458102783782662 x[17] + 0.850114386589317 x[18] + 0.8493338819701169 x[19] - 1.095166371145091 x[20] = -0.14948376404232092 \n",
" -0.7541284390976164 x[1] - 1.4941525933748385 x[2] + 1.0832040669018208 x[3] - 0.15912656304085246 x[4] + 1.751498033558951 x[5] + 0.4864199247342999 x[6] + 1.0239173561441013 x[7] + 0.06718640764610144 x[8] + 0.8900716414239433 x[9] + 0.5749389317056772 x[10] - 0.8638524329128551 x[11] + 1.792164501718729 x[12] + 0.3343406337434719 x[13] + 2.235498627467402 x[14] - 1.7776175854238803 x[15] - 1.0258963869354887 x[16] - 0.4928523771749276 x[17] + 0.39524107102528344 x[18] - 1.0638415363597291 x[19] + 0.4428400393578505 x[20] = -0.7541284390976164 \n",
" -0.35188939912748723 x[1] - 0.5564516765297487 x[2] + 0.8609077087689985 x[3] - 0.4459580353344765 x[4] - 1.3787248167701522 x[5] - 0.524354973448415 x[6] - 0.30809513724737636 x[7] + 0.4870258380831458 x[8] - 1.1679921737733099 x[9] - 0.25442927020910516 x[10] + 0.5666652944229207 x[11] - 1.0221727823880569 x[12] - 1.0936761381380715 x[13] + 0.4953375902470542 x[14] + 0.7614509551498848 x[15] + 0.34888583675755624 x[16] + 0.012279935365854578 x[17] + 1.227447616778815 x[18] - 1.0170144870456475 x[19] + 0.6173099053655823 x[20] = -0.35188939912748723 \n",
" 1.7002809122995806 x[1] - 1.3080362019937761 x[2] + 0.6698914942502485 x[3] - 1.1474334649183995 x[4] + 0.11050526255250609 x[5] + 0.8995290099838263 x[6] - 0.6634478105530445 x[7] - 0.05235412305469047 x[8] - 0.02988514100020276 x[9] + 1.1521383853656486 x[10] - 0.06063301729704743 x[11] + 0.8365335929703048 x[12] - 0.31223936865583823 x[13] + 0.02602500549162202 x[14] - 0.035716247200114765 x[15] + 0.5853579238128782 x[16] - 0.6052857888172622 x[17] - 0.6293520497536604 x[18] - 0.4989225029374304 x[19] + 0.0007411896423990069 x[20] = 1.7002809122995806\n",
" -1.9143586614424437 x[1] + 0.2332720816757861 x[2] - 0.9108041941774613 x[3] + 0.6718497718092494 x[4] - 0.3864706275744775 x[5] + 0.9430354905454666 x[6] - 1.9040919552655386 x[7] + 1.9572258747799127 x[8] - 0.18713576853924896 x[9] - 1.6269066173163877 x[10] + 0.28418540958548016 x[11] - 0.04699744228290746 x[12] + 1.2886659385255383 x[13] + 0.3513225094666955 x[14] + 0.5512684747303959 x[15] - 0.08891960612658953 x[16] - 1.8808964982943022 x[17] - 0.2376802282037522 x[18] - 0.6731284617927262 x[19] - 1.2158663083813768 x[20] = -1.9143586614424437 \n",
" 0.310830838565587 x[1] + 0.8888270096032808 x[2] + 0.7577760165487661 x[3] - 0.15518885359100631 x[4] - 0.7684844714001808 x[5] + 0.7497799172643206 x[6] - 0.35062790918196557 x[7] - 2.6337340518702805 x[8] - 1.4851587993581619 x[9] - 0.43480108607650014 x[10] - 0.35573584983593276 x[11] + 0.4466884954755695 x[12] + 0.06879200480357561 x[13] + 0.34490304362831997 x[14] + 0.22451760846655378 x[15] - 1.1339988801179888 x[16] - 0.08390543983509796 x[17] - 0.3661366125564423 x[18] - 1.125539882856653 x[19] + 0.549895815884775 x[20] = 0.310830838565587 \n",
" 0.3351183436255904 x[1] - 0.20224029527764964 x[2] + 0.37287750283043997 x[3] + 0.9197793389347283 x[4] - 2.3555443079883993 x[5] + 1.2171416584700034 x[6] + 1.063852102607513 x[7] - 0.6863502981597671 x[8] + 1.3033516262799636 x[9] + 0.30786257474128864 x[10] + 0.15560180496178777 x[11] - 1.2368309277212384 x[12] - 0.9645565854600385 x[13] + 0.2816474936831833 x[14] + 0.6622404850086145 x[15] - 1.3613390912274006 x[16] - 0.323996338035808 x[17] + 0.5074502141260505 x[18] - 0.1329498331538597 x[19] - 0.29052657932964593 x[20] = 0.3351183436255904 \n",
" 0.5234120006014658 x[1] + 0.15127595075684833 x[2] + 1.9741941124149873 x[3] - 2.3005881197896825 x[4] + 0.6298493295838168 x[5] - 2.3097456337123443 x[6] + 1.927549279609831 x[7] + 0.4347452916044357 x[8] - 0.05169266614265093 x[9] - 1.179174687395125 x[10] - 2.0405758671455123 x[11] - 0.027821906466193414 x[12] - 0.8486853210219902 x[13] + 0.13773693582241403 x[14] - 0.8753206403726591 x[15] + 0.43341447986052756 x[16] + 0.13045667892447826 x[17] + 0.9359966432611737 x[18] + 0.649770818201045 x[19] + 1.2941694965906838 x[20] = 0.5234120006014658 "
]
},
"execution_count": 61,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@constraint(model, A*x .== b)"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:07:34.018000-04:00",
"start_time": "2020-04-16T18:07:33.749Z"
}
},
"outputs": [
{
"data": {
"text/latex": [
"$$ y_{1} + y_{2} + y_{3} + y_{4} + y_{5} + y_{6} + y_{7} + y_{8} + y_{9} + y_{10} + y_{11} + y_{12} + y_{13} + y_{14} + y_{15} + y_{16} + y_{17} + y_{18} + y_{19} + y_{20} $$"
],
"text/plain": [
"y[1] + y[2] + y[3] + y[4] + y[5] + y[6] + y[7] + y[8] + y[9] + y[10] + y[11] + y[12] + y[13] + y[14] + y[15] + y[16] + y[17] + y[18] + y[19] + y[20]"
]
},
"execution_count": 62,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"@variable(model, y[1:n])\n",
"@constraint(model, x .<= y)\n",
"@constraint(model, -x .<= y)\n",
"@objective(model, Min, sum(y))"
]
},
{
"cell_type": "code",
"execution_count": 63,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:07:40.945000-04:00",
"start_time": "2020-04-16T18:07:40.898Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"----------------------------------------------------------------------------\n",
"\tSCS v2.1.1 - Splitting Conic Solver\n",
"\t(c) Brendan O'Donoghue, Stanford University, 2012\n",
"----------------------------------------------------------------------------\n",
"Lin-sys: sparse-indirect, nnz in A = 280, CG tol ~ 1/iter^(2.00)\n",
"eps = 1.00e-05, alpha = 1.50, max_iters = 5000, normalize = 1, scale = 1.00\n",
"acceleration_lookback = 10, rho_x = 1.00e-03\n",
"Variables n = 40, constraints m = 50\n",
"Cones:\tprimal zero / dual free vars: 10\n",
"\tlinear vars: 40\n",
"Setup time: 1.03e-04s\n",
"----------------------------------------------------------------------------\n",
" Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)\n",
"----------------------------------------------------------------------------\n",
" 0| 3.74e+19 6.44e+19 1.00e+00 -2.10e+20 1.43e+20 9.17e+19 7.25e-05 \n",
" 40| 1.62e-07 8.12e-07 4.43e-07 1.00e+00 1.00e+00 2.06e-16 1.32e-03 \n",
"----------------------------------------------------------------------------\n",
"Status: Solved\n",
"Timing: Solve time: 1.32e-03s\n",
"\tLin-sys: avg # CG iterations: 2.76, avg solve time: 1.54e-05s\n",
"\tCones: avg projection time: 1.45e-07s\n",
"\tAcceleration: avg step time: 1.32e-05s\n",
"----------------------------------------------------------------------------\n",
"Error metrics:\n",
"dist(s, K) = 4.2368e-17, dist(y, K*) = 0.0000e+00, s'y/|s||y| = 5.8604e-18\n",
"primal res: |Ax + s - b|_2 / (1 + |b|_2) = 1.6239e-07\n",
"dual res: |A'y + c|_2 / (1 + |c|_2) = 8.1246e-07\n",
"rel gap: |c'x + b'y| / (1 + |c'x| + |b'y|) = 4.4289e-07\n",
"----------------------------------------------------------------------------\n",
"c'x = 1.0000, -b'y = 1.0000\n",
"============================================================================\n"
]
}
],
"source": [
"optimize!(model)"
]
},
{
"cell_type": "code",
"execution_count": 64,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:07:52.700000-04:00",
"start_time": "2020-04-16T18:07:52.688Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"20-element Array{Float64,1}:\n",
" 1.0000002097854823 \n",
" -3.7194676307110125e-8\n",
" 1.3057975606922636e-8\n",
" -6.652681661529527e-9 \n",
" 2.0993775947052414e-8\n",
" -6.769918215855081e-9 \n",
" -2.06070404672574e-8 \n",
" -2.1598410376842593e-8\n",
" 6.107184897131848e-9 \n",
" -3.68665891009158e-8 \n",
" 3.4011270251743704e-8\n",
" 9.123559643923159e-10\n",
" -2.31421917475553e-8 \n",
" -5.048228113906143e-8 \n",
" 8.591719552765204e-9 \n",
" -1.3121841985150527e-8\n",
" -5.4338447222664766e-8\n",
" 1.504538540165436e-8 \n",
" 2.625446251349709e-8 \n",
" -5.771370350149791e-8 "
]
},
"execution_count": 64,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xs = value.(x)"
]
},
{
"cell_type": "code",
"execution_count": 65,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:08:01.883000-04:00",
"start_time": "2020-04-16T18:08:01.846Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"6.632473941796915e-7"
]
},
"execution_count": 65,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"onenorm(xs - x0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's make a function out of that."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T13:01:28.612000-04:00",
"start_time": "2020-04-16T16:59:51.893Z"
}
},
"outputs": [],
"source": [
"function BP(A,b; optimizer = opt_scs)\n",
" \n",
" d, n = size(A)\n",
" @assert length(b) == d\n",
" \n",
" model = Model(optimizer) \n",
" @variable(model, x[1:n])\n",
" @constraint(model, A*x .== b)\n",
" \n",
" @variable(model, y[1:n])\n",
" @constraint(model, x .<= y)\n",
" @constraint(model, -x .<= y)\n",
" @objective(model, Min, sum(y))\n",
" \n",
" optimize!(model)\n",
" \n",
" xs = value.(x)\n",
" \n",
" return xs\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 66,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:08:34.961000-04:00",
"start_time": "2020-04-16T18:08:34.933Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"20-element Array{Float64,1}:\n",
" 1.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0\n",
" 0.0"
]
},
"execution_count": 66,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xs = BP(A,b, optimizer=opt_glpk)"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:08:52.088000-04:00",
"start_time": "2020-04-16T18:08:52.080Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.0"
]
},
"execution_count": 67,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"onenorm(xs-x0)"
]
},
{
"cell_type": "code",
"execution_count": 68,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:09:07.609000-04:00",
"start_time": "2020-04-16T18:09:07.587Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3.5449407029009996e-17"
]
},
"execution_count": 68,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xs = BP(A,b, optimizer=opt_mosek)\n",
"onenorm(xs-x0)"
]
},
{
"cell_type": "code",
"execution_count": 69,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:09:31.643000-04:00",
"start_time": "2020-04-16T18:09:31.566Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3-element Array{Int64,1}:\n",
" 1\n",
" 3\n",
" 6"
]
},
"execution_count": 69,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"randset(n,k) = randperm(n)[1:k]\n",
"randset(10,3)"
]
},
{
"cell_type": "code",
"execution_count": 70,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:09:55.675000-04:00",
"start_time": "2020-04-16T18:09:55.663Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"200-element Array{Float64,1}:\n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" ⋮ \n",
" 0.0 \n",
" -1.0444736865400661\n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 \n",
" 0.0 "
]
},
"execution_count": 70,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"d = 50\n",
"n = 200\n",
"k = 7\n",
"S = randset(n,k)\n",
"x0 = zeros(n)\n",
"x0[S] = randn(k)\n",
"x0"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:10:04.832000-04:00",
"start_time": "2020-04-16T18:10:04.814Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"200-element SparseVector{Float64,Int64} with 7 stored entries:\n",
" [29 ] = -0.285493\n",
" [38 ] = -0.599437\n",
" [87 ] = -1.27762\n",
" [107] = -0.580424\n",
" [123] = 0.915584\n",
" [143] = 0.321277\n",
" [190] = -1.04447"
]
},
"execution_count": 71,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sparse(x0)"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:10:25.560000-04:00",
"start_time": "2020-04-16T18:10:25.324Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"1.3954734856354438e-5"
]
},
"execution_count": 72,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A = randn(d,n)\n",
"b = A*x0\n",
"xs = BP(A,b)\n",
"onenorm(x0-xs)"
]
},
{
"cell_type": "code",
"execution_count": 73,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:10:52.764000-04:00",
"start_time": "2020-04-16T18:10:52.705Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"([1.4302075533475465 0.7520118815077361 … -1.4701788233252002 -0.35555099557525566; 0.7812316922601905 0.19914232040424418 … -0.9954739771581482 0.26420375974891974; … ; 0.7734909444638944 -0.6865218343782568 … -0.19714239408002335 -0.3838898158942502; 0.8647792696114751 1.0643908284555994 … 1.4109178925704224 0.8426275793521971], [1.7650602445585468, 4.6057093169563865, -1.7296368096215045, -0.19432879001198017, 0.9264485353718146, -0.1862857479190762, 1.0405629461544077, -1.1371974562444191, -0.3326791218206496, -6.230899766426076 … -2.734486046085591, -2.61266368229103, -3.489783791060211, 5.010377527588336, -0.731604765215425, -2.674764514776432, -1.9619302473929523, 5.069184704505968, -0.9854364907820127, 1.0245333091398552], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.1945546502452986, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])"
]
},
"execution_count": 73,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function make_problem(d,n,k) \n",
" S = randset(n,k)\n",
" x0 = zeros(n)\n",
" x0[S] = randn(k)\n",
" A = randn(d,n)\n",
" b = A*x0; \n",
" return A,b,x0\n",
"end\n",
"A, b, x0 = make_problem(d, n, k)"
]
},
{
"cell_type": "code",
"execution_count": 74,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:11:13.429000-04:00",
"start_time": "2020-04-16T18:11:08.033Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"0.3553601360338345"
]
},
"execution_count": 74,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A, b, x0 = make_problem(50, 200, 15)\n",
"xs = BP(A,b)\n",
"onenorm(x0-xs)"
]
},
{
"cell_type": "code",
"execution_count": 75,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:11:23.553000-04:00",
"start_time": "2020-04-16T18:11:23.486Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"9.25872226678312"
]
},
"execution_count": 75,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"onenorm(x0)"
]
},
{
"cell_type": "code",
"execution_count": 76,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:11:44.551000-04:00",
"start_time": "2020-04-16T18:11:38.541Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 5.881336 seconds (61.02 k allocations: 9.716 MiB)\n",
" 0.048013 seconds (30.03 k allocations: 3.071 MiB)\n",
" 0.066250 seconds (40.00 k allocations: 4.095 MiB)\n"
]
}
],
"source": [
"@time xs = BP(A,b,optimizer=opt_scs)\n",
"@time xs = BP(A,b,optimizer=opt_glpk)\n",
"@time xs = BP(A,b,optimizer=opt_mosek)\n",
";"
]
},
{
"cell_type": "code",
"execution_count": 77,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:12:08.826000-04:00",
"start_time": "2020-04-16T18:12:06.509Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.576325 seconds (79.33 k allocations: 16.935 MiB, 1.69% gc time)\n",
" 0.660002 seconds (104.77 k allocations: 19.672 MiB)\n"
]
},
{
"data": {
"text/plain": [
"(42.337713825009715, 42.33771382500943)"
]
},
"execution_count": 77,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A, b, x0 = make_problem(200, 500, 100)\n",
"@time x1 = BP(A,b,optimizer=opt_glpk)\n",
"@time x2 = BP(A,b,optimizer=opt_mosek)\n",
"onenorm(x1-x0), onenorm(x2-x0)"
]
},
{
"cell_type": "code",
"execution_count": 78,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:12:36.289000-04:00",
"start_time": "2020-04-16T18:12:34.246Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 1.098129 seconds (79.33 k allocations: 16.935 MiB)\n",
" 0.914337 seconds (104.77 k allocations: 19.700 MiB)\n"
]
},
{
"data": {
"text/plain": [
"(6.133142236098867e-13, 1.1976518734336954e-13)"
]
},
"execution_count": 78,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A, b, x0 = make_problem(200, 500, 40)\n",
"@time x1 = BP(A,b,optimizer=opt_glpk)\n",
"@time x2 = BP(A,b,optimizer=opt_mosek)\n",
"onenorm(x1-x0), onenorm(x2-x0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# BPDN"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T13:01:49.168000-04:00",
"start_time": "2020-04-16T16:59:51.924Z"
}
},
"outputs": [],
"source": [
"function make_problem(d,n,k,ϵ) \n",
" S = randset(n,k)\n",
" x0 = zeros(n)\n",
" x0[S] = randn(k)\n",
" \n",
" A = randn(d,n)\n",
" for i in 1:n\n",
" A[:,i] = A[:,i] / norm(A[:,i])\n",
" end\n",
" \n",
" e = randn(d)\n",
" e = ϵ * e / norm(e)\n",
" b = A*x0 + e\n",
" \n",
" return A, b, x0\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 79,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:13:07.969000-04:00",
"start_time": "2020-04-16T18:13:07.890Z"
},
"scrolled": false
},
"outputs": [
{
"data": {
"text/plain": [
"([-0.31455647227518035 -0.008090775210207369 … -0.17042816394751953 -0.19542786176733629; -0.05731604790835892 0.1868328526942932 … -0.06236958591823509 -0.34149074468279744; … ; 0.20562128249380562 -0.4121349455984317 … 0.14460064093060726 0.1688491475418608; 0.38116480151380894 0.4977991652864089 … 0.09710406080012117 0.11245830487696122], [-0.9056498409818027, -1.400990935048949, 0.18839297043983338, -0.05102212802650892, -0.8472918886171541, -0.03799898960726512, 0.08314299597279266, -0.4346907602594959, -0.34503047689841976, 0.6672663786293269], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.754102336582079, 0.0, 0.0])"
]
},
"execution_count": 79,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A, b, x0 = make_problem(10,50,2,0.1)"
]
},
{
"cell_type": "code",
"execution_count": 80,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:13:27.716000-04:00",
"start_time": "2020-04-16T18:13:27.609Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"BPDN (generic function with 1 method)"
]
},
"execution_count": 80,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function BPDN(A,b,ϵ; optimizer = opt_scs)\n",
" \n",
" d, n = size(A)\n",
" @assert length(b) == d\n",
" \n",
" model = Model(optimizer) \n",
" @variable(model, x[1:n])\n",
" @constraint(model, [ϵ; A*x - b] in SecondOrderCone()) # norm(Ax - b) ≤ ϵ\n",
" \n",
" @variable(model, y[1:n])\n",
" @constraint(model, x .<= y)\n",
" @constraint(model, -x .<= y)\n",
" @objective(model, Min, sum(y))\n",
" \n",
" optimize!(model)\n",
" \n",
" xs = value.(x)\n",
" \n",
" return xs\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 81,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:13:32.278000-04:00",
"start_time": "2020-04-16T18:13:31.826Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(0.114430905146159, 0.2304691248833006)"
]
},
"execution_count": 81,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xs = BPDN(A,b,0.1)\n",
"norm(xs-x0), onenorm(xs-x0)"
]
},
{
"cell_type": "code",
"execution_count": 82,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:14:16.620000-04:00",
"start_time": "2020-04-16T18:14:16.616Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(2.1189445683491557, 2.025754594527242)"
]
},
"execution_count": 82,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"norm(x0), norm(xs)"
]
},
{
"cell_type": "code",
"execution_count": 83,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:14:28.828000-04:00",
"start_time": "2020-04-16T18:14:28.776Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(0.01502772364275557, 1.2763608041825145)"
]
},
"execution_count": 83,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A, b, x0 = make_problem(10,50,2,0.01)\n",
"xs = BPDN(A,b,0.01)\n",
"norm(xs-x0), norm(x0)"
]
},
{
"cell_type": "code",
"execution_count": 84,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:14:33.746000-04:00",
"start_time": "2020-04-16T18:14:33.543Z"
}
},
"outputs": [
{
"ename": "ErrorException",
"evalue": "Constraints of type MathOptInterface.VectorAffineFunction{Float64}-in-MathOptInterface.SecondOrderCone are not supported by the solver and there are no bridges that can reformulate it into supported constraints.",
"output_type": "error",
"traceback": [
"Constraints of type MathOptInterface.VectorAffineFunction{Float64}-in-MathOptInterface.SecondOrderCone are not supported by the solver and there are no bridges that can reformulate it into supported constraints.",
"",
"Stacktrace:",
" [1] error(::String) at ./error.jl:33",
" [2] moi_add_constraint(::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.AbstractOptimizer,MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, ::MathOptInterface.VectorAffineFunction{Float64}, ::MathOptInterface.SecondOrderCone) at /Users/spielman/.julia/packages/JuMP/MnJQc/src/constraints.jl:389",
" [3] add_constraint(::Model, ::VectorConstraint{GenericAffExpr{Float64,VariableRef},MathOptInterface.SecondOrderCone,VectorShape}, ::String) at /Users/spielman/.julia/packages/JuMP/MnJQc/src/constraints.jl:403",
" [4] macro expansion at /Users/spielman/.julia/packages/JuMP/MnJQc/src/macros.jl:382 [inlined]",
" [5] #BPDN#27(::MathOptInterface.OptimizerWithAttributes, ::typeof(BPDN), ::Array{Float64,2}, ::Array{Float64,1}, ::Float64) at ./In[80]:8",
" [6] (::var\"#kw##BPDN\")(::NamedTuple{(:optimizer,),Tuple{MathOptInterface.OptimizerWithAttributes}}, ::typeof(BPDN), ::Array{Float64,2}, ::Array{Float64,1}, ::Float64) at ./none:0",
" [7] top-level scope at In[84]:2"
]
}
],
"source": [
"A, b, x0 = make_problem(100,500,20,0.1)\n",
"xs = BPDN(A,b,0.01, optimizer=opt_glpk)\n",
"norm(xs-x0), norm(x0)"
]
},
{
"cell_type": "code",
"execution_count": 85,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:14:46.072000-04:00",
"start_time": "2020-04-16T18:14:45.795Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"(0.34580576413952496, 3.77540385217923)"
]
},
"execution_count": 85,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A, b, x0 = make_problem(100,500,20,0.1)\n",
"xs = BPDN(A,b,0.01, optimizer=opt_mosek)\n",
"norm(xs-x0), norm(x0)"
]
},
{
"cell_type": "code",
"execution_count": 86,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:15:01.500000-04:00",
"start_time": "2020-04-16T18:15:01.404Z"
}
},
"outputs": [
{
"data": {
"text/plain": [
"plot_compare (generic function with 1 method)"
]
},
"execution_count": 86,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function plot_compare(x0, xs, k)\n",
" @assert any(x0 .== 0) # make sure got order of x0 vs xs correct\n",
" p = sortperm(abs.(xs) + abs.(x0),rev=true)\n",
" groupedbar(hcat(x0[p[1:k]] , xs[p[1:k]]), label=hcat(\"x0\", \"xs\"))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 87,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:15:07.061000-04:00",
"start_time": "2020-04-16T18:15:06.702Z"
}
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 87,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot_compare(x0, xs, 30)"
]
},
{
"cell_type": "code",
"execution_count": 88,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:15:47.953000-04:00",
"start_time": "2020-04-16T18:15:47.678Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(norm(xs - x0), norm(x0)) = (0.8268405421419495, 1.9021607993177583)\n"
]
},
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 88,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"A, b, x0 = make_problem(10, 30, 6)\n",
"xs = BP(A,b)\n",
"@show norm(xs-x0), norm(x0)\n",
"plot_compare(x0, xs, 10)"
]
},
{
"cell_type": "code",
"execution_count": 90,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T14:16:37.559000-04:00",
"start_time": "2020-04-16T18:16:36.758Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" 0.449484 seconds (87.56 k allocations: 17.339 MiB)\n",
" 0.276066 seconds (87.81 k allocations: 28.340 MiB, 8.63% gc time)\n"
]
}
],
"source": [
"d = 200\n",
"n = 400\n",
"k = 50\n",
"A, b = make_problem(d, n, k)\n",
"@time xs = BP(A,b, optimizer=opt_mosek)\n",
"\n",
"A, b, x0 = make_problem(d, n, k, 0.1)\n",
"@time xs = BPDN(A,b, 0.1, optimizer=opt_mosek)\n",
";\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2020-04-16T13:02:33.864000-04:00",
"start_time": "2020-04-16T16:59:51.954Z"
}
},
"outputs": [],
"source": [
"norm(xs-x0), norm(x0)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"@webio": {
"lastCommId": null,
"lastKernelId": null
},
"kernelspec": {
"display_name": "Julia 1.3.1",
"language": "julia",
"name": "julia-1.3"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.3.1"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}