import%20marimo%0A%0A__generated_with%20%3D%20%220.25.2%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20import%20sys%0A%20%20%20%20from%20pathlib%20import%20Path%0A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%0A%20%20%20%20start%20%3D%20Path(__file__).resolve().parent%20if%20%22__file__%22%20in%20globals()%20else%20Path.cwd().resolve()%0A%20%20%20%20for%20candidate%20in%20(start%2C%20*start.parents)%3A%0A%20%20%20%20%20%20%20%20if%20(candidate%20%2F%20%22fomtools%22).is_dir()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20repo_root%20%3D%20candidate%0A%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20repo_root%20%3D%20start%0A%0A%20%20%20%20if%20str(repo_root)%20not%20in%20sys.path%3A%0A%20%20%20%20%20%20%20%20sys.path.insert(0%2C%20str(repo_root))%0A%0A%20%20%20%20from%20docs.tutorials.helper_functions.lasso%20import%20(%0A%20%20%20%20%20%20%20%20certify_lasso_pepit_radius%2C%0A%20%20%20%20%20%20%20%20ista_step%2C%0A%20%20%20%20%20%20%20%20make_lasso_data%2C%0A%20%20%20%20)%0A%0A%20%20%20%20np.set_printoptions(precision%3D5)%20%20%23%20Print%20few%20decimal%20places%0A%20%20%20%20np.set_printoptions(suppress%3DTrue)%20%20%23%20Suppress%20scientific%20notation%0A%20%20%20%20return%20certify_lasso_pepit_radius%2C%20ista_step%2C%20make_lasso_data%2C%20mo%2C%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Problem%20setup%3A%20Sparse%20coding%0A%0A%20%20%20%20%5BBack%20to%20the%20verification%20tutorials%5D(.%2F)%0A%0A%20%20%20%20In%20this%20tutorial%2C%20we%20verify%20first-order%20methods%20for%20the%20sparse%20coding%20example%20from%20Section%204.1%20of%20%5BExact%20Verification%20of%20First-Order%20Methods%20via%20Mixed-Integer%20Linear%20Programming%5D(https%3A%2F%2Farxiv.org%2Fabs%2F2412.11330).%0A%0A%20%20%20%20Consider%20the%20following%20problem%20of%20recovering%20the%20signal%20%24z_%5Cmathrm%7Btrue%7D%24%20from%20its%20noisy%20measurement%20%24x%20%3D%20D%20z_%5Cmathrm%7Btrue%7D%20%2B%20%5Cvarepsilon%20%5Cin%20%5Cmathbf%7BR%7D%5Ep%24%3A%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20(1%2F2)%20%5C%7CDz%20-%20x%5C%7C_2%5E2%20%2B%20%5Clambda%20%5C%7Cz%5C%7C_1%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20where%20%24D%20%5Cin%20%5Cmathbf%7BR%7D%5E%7Bp%20%5Ctimes%20n%7D%24%20is%20the%20dictionary%20matrix%2C%20%24%5Cvarepsilon%20%5Cin%20%5Cmathbf%7BR%7D%5Ep%24%2C%20%24%5Clambda%20%3E%200%24%20is%20a%20sparsity-inducing%20regularization%20parameter%2C%20%24z%20%5Cin%20%5Cmathbf%7BR%7D%5En%24%20is%20the%20decision%20variable%2C%20and%20%24x%20%5Cin%20%5Cmathbf%7BR%7D%5Ep%24%20is%20the%20problem%20parameter.%0A%0A%20%20%20%20%23%23%20First-order%20method%3A%20ISTA%20and%20FISTA%0A%0A%20%20%20%20This%20problem%20can%20be%20solved%20with%20ISTA%20and%20its%20accelerated%20variant%20FISTA.%0A%20%20%20%20We%20provide%20the%20definition%20of%20ISTA%20iteration%20below%3A%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%20%20%20%20z%5E%7Bk%2B1%7D%20%3D%20%5Cmathbf%7Bsoft%7D_%7B%5Clambda%20%5Ceta%7D%5Cleft(z%5E%7Bk%7D%20-%20%5Ceta%20D%5ET(Dz%5E%7Bk%7D%20-%20x)%5Cright)%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20where%20%24%5Cmathbf%7Bsoft%7D_%5Ckappa(v)_i%20%3D%20%5Cmathbf%7Bsign%7D(v_i)%5Cmax%5C%7B%7Cv_i%7C%20-%20%5Ckappa%2C%200%5C%7D%24%20is%20the%20soft-thresholding%20operator%2C%20applied%20element-wise%20(the%20%60soft_threshold%60%20step).%0A%20%20%20%20We%20choose%20%24%5Ceta%20%3D%201%2FL%24%2C%20where%20%24L%24%20is%20the%20maximum%20eigenvalue%20of%20%24D%5ET%20D%24.%0A%20%20%20%20The%20FISTA%20iteration%20can%20be%20written%20as%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Brcl%7D%0A%20%20%20%20%20%20%20%20y%5E%7Bk%2B1%7D%20%26%3D%26%20%5Cmathbf%7Bsoft%7D_%7B%5Clambda%20%5Ceta%7D%5Cleft(z%5E%7Bk%7D%20-%20%5Ceta%20D%5ET(Dz%5E%7Bk%7D%20-%20x)%5Cright)%2C%20%5C%5C%0A%20%20%20%20%20%20%20%20z%5E%7Bk%2B1%7D%20%26%3D%26%20y%5E%7Bk%2B1%7D%20%2B%20%5Cbeta_%7Bk%2B1%7D(y%5E%7Bk%2B1%7D%20-%20y%5E%7Bk%7D)%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20where%20%24y%5E0%20%3D%20z%5E0%24%2C%20%24t_0%20%3D%201%24%2C%20%24t_k%20%3D%20(1%20%2B%20%5Csqrt%7B1%20%2B%204t_%7Bk-1%7D%5E2%7D)%2F2%24%2C%20and%20%24%5Cbeta_k%20%3D%20(t_%7Bk-1%7D-1)%2Ft_k%24%2C%20so%20that%20%24%5Cbeta_1%20%3D%200%24.%0A%0A%20%20%20%20The%20performance%20metric%20is%20the%20infinity%20norm%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24%20of%20the%20fixed-point%20residual%20%24r%5Ek%20%3D%20y%5Ek%20-%20y%5E%7Bk-1%7D%24%20between%20consecutive%20proximal%20iterates%2C%20as%20in%20the%20paper%3B%20it%20keeps%20the%20verification%20problem%20a%20mixed-integer%20*linear*%20program.%0A%20%20%20%20For%20ISTA%2C%20%24y%5Ek%20%3D%20z%5Ek%24%2C%20so%20%24r%5Ek%20%3D%20z%5Ek%20-%20z%5E%7Bk-1%7D%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Problem%20parameters%0A%0A%20%20%20%20For%20the%20strongly%20convex%20example%2C%20we%20choose%20%24p%3D20%24%20and%20%24n%3D15%24%2C%20and%20for%20the%20nonstrongly%20convex%20example%2C%20we%20choose%20%24p%20%3D%2015%24%2C%20%24n%20%3D%2020%24%2C%20and%20%24%5Clambda%20%3D%2010%5E%7B-1%7D%24.%0A%20%20%20%20To%20construct%20%24D%24%2C%20we%20randomly%20sample%20from%20i.i.d.%20Gaussian%20%24%5Cmathcal%7BN%7D(0%2C%201%2Fp)%24%20the%20%2420%24%25%20nonzeros%20of%20%24D%24%20and%20normalize%20its%20columns.%0A%20%20%20%20The%20true%20signal%20%24z_%5Cmathrm%7Btrue%7D%24%20is%20constructed%20by%20sampling%20each%20entry%20from%20i.i.d.%20Gaussian%20%24%5Cmathcal%7BN%7D(0%2C%201)%24%20with%20probability%20%240.1%24%20and%20setting%20%240%24%20otherwise.%20We%20sample%20%24m%20%3D%20100%24%20of%20those%20random%20sparse%20signals%20%24z_%5Cmathrm%7Btrue%7D%24.%0A%20%20%20%20Finally%2C%20we%20let%20the%20parameter%20set%20%24%5Cmathcal%7BX%7D%20%3D%20%5Bx_%5Cmathrm%7Blb%7D%2C%20x_%5Cmathrm%7Bub%7D%5D%24%20be%20the%20tightest%20box%20containing%20all%20%24x%24's.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(make_lasso_data)%3A%0A%20%20%20%20%23%20p%2C%20n%2C%20lambd%20%3D%2015%2C%2020%2C%201e-1%0A%20%20%20%20p%2C%20n%2C%20lambd%20%3D%203%2C%204%2C%201e-1%0A%20%20%20%20%23%20num_star%20%3D%20100%0A%20%20%20%20num_star%20%3D%2020%0A%20%20%20%20D%2C%20x_lb%2C%20x_ub%2C%20L%2C%20eta%2C%20lambd%20%3D%20make_lasso_data(%0A%20%20%20%20%20%20%20%20p%3Dp%2C%20n%3Dn%2C%20lambd%3Dlambd%2C%20num_signals%3Dnum_star%2C%20seed%3D0%0A%20%20%20%20)%0A%20%20%20%20s0_lb%2C%20s0_ub%20%3D%20-1.0%2C%201.0%0A%20%20%20%20return%20D%2C%20L%2C%20eta%2C%20lambd%2C%20n%2C%20p%2C%20s0_lb%2C%20s0_ub%2C%20x_lb%2C%20x_ub%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Every%20verification%20mode%20below%20(simulator%20and%20MIP-VP)%20drives%20this%0A%20%20%20%20*same*%20%60ista_step%60%20recurrence%20and%20residual%3B%20only%20the%20%60Verifier%60%20backend%20and%0A%20%20%20%20whether%20%24(x%2C%20z%5E0)%24%20are%20sampled%20instances%20or%20bounded%20search%20variables%20change.%0A%20%20%20%20Setting%20%60momentum%3DTrue%60%20selects%20the%20FISTA%20variant.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Verification%20Problem%0A%0A%20%20%20%20In%20this%20experiment%2C%20we%20verify%20the%20performance%20of%20ISTA%20and%20FISTA%20for%20sparse%20coding%20given%20the%20iteration%20budget%20%24K%3E0%24.%20We%20compare%20simulator%20statistics%2C%20exact%20MIP%20verification%2C%20PEPit%20bounds%2C%20and%20the%20in-sample%20worst-case%20instances.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20K%20%3D%2022%20%20%23%20Maximum%20iteration%20number%0A%20%20%20%20return%20(K%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Simulator%20class%0A%0A%20%20%20%20We%20first%20simulate%20ISTA%20on%20given%20problem%20instances%2C%20which%20are%20parameterized%20by%20the%20parameter%20%24x%24.%0A%20%20%20%20Each%20problem%20instance%20is%20generated%20from%20randomly%20sampling%20%24x%24%20from%20uniform%20distribution%20over%20%24%5Cmathcal%7BX%7D%20%3D%20%5Bx_%5Cmathrm%7Blb%7D%2C%20x_%5Cmathrm%7Bub%7D%5D%24%20and%20the%20initial%20iterate%20%24z%5E0%24%20from%20uniform%20distribution%20over%20%24%5Cmathcal%7BZ%7D%20%3D%20%5Bz%5E0_%5Cmathrm%7Blb%7D%2C%20z%5E0_%5Cmathrm%7Bub%7D%5D%24%20(%60s0_lb%60%2C%20%60s0_ub%60%2C%20and%20%60s0_samples%60%20in%20the%20code).%20The%20samples%20are%20drawn%20once%20below%20and%20reused%20by%20every%20simulation%20cell%20so%20that%20all%20comparisons%20share%20the%20same%20instances.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(n%2C%20np%2C%20p%2C%20s0_lb%2C%20s0_ub%2C%20x_lb%2C%20x_ub)%3A%0A%20%20%20%20n_samples%20%3D%2010000%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(0)%0A%20%20%20%20x_samples%20%3D%20_rng.uniform(low%3Dx_lb%2C%20high%3Dx_ub%2C%20size%3D(n_samples%2C%20p))%0A%20%20%20%20s0_samples%20%3D%20_rng.uniform(low%3Ds0_lb%2C%20high%3Ds0_ub%2C%20size%3D(n_samples%2C%20n))%0A%20%20%20%20return%20s0_samples%2C%20x_samples%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20simulation%20is%20run%20through%20the%20following%20function.%20The%20output%20of%0A%20%20%20%20%60_vp.solve()%60%20is%20an%20array%20of%20size%20%60(n_samples%2C%20K)%60%20whose%20%24i%24-th%20row%20lists%0A%20%20%20%20the%20residual%20norms%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24%20for%20%24k%3D1%2C%5Cldots%2CK%24%20of%20the%20%24i%24-th%0A%20%20%20%20problem%20instance.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(D%2C%20eta%2C%20ista_step%2C%20lambd%2C%20n%2C%20np%2C%20p)%3A%0A%20%20%20%20import%20fomtools%20as%20ft%0A%0A%20%20%20%20def%20maximize_residual(vp%2C%20residual%2C%20ord)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Maximize%20%60%60%7C%7Cresidual%7C%7C_ord%60%60%2C%20squared%20for%20%60%60ord%3D2%60%60%20(%60%60ft.sum_squares%60%60).%22%22%22%0A%20%20%20%20%20%20%20%20vp.maximize(ft.sum_squares(residual)%20if%20ord%20%3D%3D%202%20else%20ft.norm(residual%2C%20ord))%0A%0A%20%20%20%20def%20ista_simulation(x_samples%2C%20s0_samples%2C%20K%2C%20momentum%3DFalse%2C%20ord%3Dnp.inf)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Simulate%20ISTA%20on%20sampled%20instances.%0A%0A%20%20%20%20%20%20%20%20Builds%20a%20simulator-mode%20%3Aclass%3A%60Verifier%60%2C%20registers%20the%20parameter%0A%20%20%20%20%20%20%20%20%60%60x%60%60%20and%20initial%20iterate%20%60%60s0%60%60%20as%20sampled%20parameters%2C%20rolls%20out%20the%0A%20%20%20%20%20%20%20%20recurrence%20with%20%60%60ista_step%60%60%2C%20and%20simulates%20via%20%60%60vp.solve()%60%60.%20Each%0A%20%20%20%20%20%20%20%20%60%60*_samples%60%60%20array%20has%20shape%20%60%60(n_samples%2C%20dim)%60%60%20(a%20single-row%20array%0A%20%20%20%20%20%20%20%20replays%20one%20instance)%2C%20and%20the%20returned%20%60%60residuals%60%60%20has%20shape%0A%20%20%20%20%20%20%20%20%60%60(n_samples%2C%20K)%60%60.%20%60%60ord%60%60%20defaults%20to%20%60%60inf%60%60%2C%20matching%0A%20%20%20%20%20%20%20%20%3Afunc%3A%60ista_verification%60%3B%20with%20%60%60ord%3D2%60%60%20the%20objective%20is%0A%20%20%20%20%20%20%20%20%60%60ft.sum_squares%60%60%2C%20the%20*squared*%20Euclidean%20norm%2C%20so%20the%20returned%0A%20%20%20%20%20%20%20%20residuals%20are%20in%20squared%20units.%20Summary%20statistics%20(%60%60mean%60%60%2C%20%60%60min%60%60%2C%20%60%60max%60%60%2C%0A%20%20%20%20%20%20%20%20...)%20are%20available%20on%20the%20returned%20instance%20%60%60vp%60%60%20via%20%60%60vp%5B...%5D%60%60.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20vp%20%3D%20ft.Verifier(method%3D'simulator'%2C%20obbt_solver%3DNone)%0A%20%20%20%20%20%20%20%20x_param%20%3D%20vp.parameter(ndim%3Dp%2C%20samples%3Dx_samples)%0A%20%20%20%20%20%20%20%20s0%20%3D%20vp.parameter(n%2C%20samples%3Ds0_samples)%0A%0A%20%20%20%20%20%20%20%20z%2C%20y%2C%20t%20%3D%20s0%2C%20s0%2C%201.0%0A%20%20%20%20%20%20%20%20for%20k%20in%20range(1%2C%20K%20%2B%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20z%2C%20y%2C%20t%2C%20residual%20%3D%20ista_step(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vp%2C%20D%2C%20eta%2C%20lambd%2C%20z%2C%20y%2C%20x_param%2C%20t%2C%20k%3Dk%2C%20momentum%3Dmomentum%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20maximize_residual(vp%2C%20residual%2C%20ord)%0A%0A%20%20%20%20%20%20%20%20residuals%20%3D%20vp.solve()%0A%20%20%20%20%20%20%20%20return%20vp%2C%20residuals%0A%0A%20%20%20%20return%20ft%2C%20ista_simulation%2C%20maximize_residual%0A%0A%0A%40app.cell%0Adef%20_(K%2C%20ista_simulation%2C%20np%2C%20s0_samples%2C%20x_samples)%3A%0A%20%20%20%20_vp%2C%20_residuals%20%3D%20ista_simulation(x_samples%2C%20s0_samples%2C%20K%2C%20momentum%3DFalse%2C%20ord%3Dnp.inf)%0A%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22Summary%20statistics%20(at%20K%3D%7BK%7D)%3A%5Cn%22%0A%20%20%20%20%20%20%20%20f%22mean%3A%20%7B_vp%5B'mean'%5D%3A.6e%7D%5Cn%22%0A%20%20%20%20%20%20%20%20f%22min%3A%20%7B_vp%5B'min'%5D%3A.6e%7D%5Cn%22%0A%20%20%20%20%20%20%20%20f%22max%3A%20%7B_vp%5B'max'%5D%3A.6e%7D%5Cn%22%0A%20%20%20%20%20%20%20%20f%22median%3A%20%7B_vp%5B'median'%5D%3A.6e%7D%5Cn%22%0A%20%20%20%20)%0A%0A%20%20%20%20simulation_residuals_ista%20%3D%20np.mean(_residuals%2C%20axis%3D0)%0A%20%20%20%20%23%20Quantile%201.0%20gives%20the%20in-sample%20worst%20case%3B%20indices%20are%20%5B90%25%2C%2050%25%2C%2010%25%2C%20max%5D.%0A%20%20%20%20simulation_residual_quantiles_ista%20%3D%20np.quantile(_residuals%2C%20%5B0.90%2C%200.50%2C%200.10%2C%201.0%5D%2C%20axis%3D0)%0A%20%20%20%20return%20simulation_residual_quantiles_ista%2C%20simulation_residuals_ista%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20further%20simulate%20FISTA%2C%20reusing%20the%20same%20%60ista_simulation%60%20helper%20and%20samples%20with%20%60momentum%3DTrue%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(K%2C%20ista_simulation%2C%20np%2C%20s0_samples%2C%20x_samples)%3A%0A%20%20%20%20_%2C%20_residuals%20%3D%20ista_simulation(x_samples%2C%20s0_samples%2C%20K%2C%20momentum%3DTrue%2C%20ord%3Dnp.inf)%0A%20%20%20%20simulation_residuals_fista%20%3D%20np.mean(_residuals%2C%20axis%3D0)%0A%20%20%20%20%23%20Quantile%201.0%20gives%20the%20in-sample%20worst%20case%3B%20indices%20are%20%5B90%25%2C%2050%25%2C%2010%25%2C%20max%5D.%0A%20%20%20%20simulation_residual_quantiles_fista%20%3D%20np.quantile(_residuals%2C%20%5B0.90%2C%200.50%2C%200.10%2C%201.0%5D%2C%20axis%3D0)%0A%20%20%20%20return%20simulation_residual_quantiles_fista%2C%20simulation_residuals_fista%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Algorithm%20verification%20through%20verification%20problem%20(VP)%0A%0A%20%20%20%20We%20now%20obtain%20worst-case%20guarantees%20by%20solving%20the%20verification%20problem%0A%20%20%20%20with%20a%20MIP%20formulation%2C%20which%20gives%20the%20exact%20worst%20case.%20The%20following%0A%20%20%20%20function%20rolls%20out%20the%20same%20%60ista_step%60%20recurrence%20over%20the%20horizon%20%24K%24%0A%20%20%20%20and%20maximizes%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24%20at%20each%20iteration.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20D%2C%0A%20%20%20%20eta%2C%0A%20%20%20%20ft%2C%0A%20%20%20%20ista_step%2C%0A%20%20%20%20lambd%2C%0A%20%20%20%20maximize_residual%2C%0A%20%20%20%20n%2C%0A%20%20%20%20np%2C%0A%20%20%20%20p%2C%0A%20%20%20%20s0_lb%2C%0A%20%20%20%20s0_ub%2C%0A%20%20%20%20x_lb%2C%0A%20%20%20%20x_ub%2C%0A)%3A%0A%20%20%20%20def%20ista_verification(K%2C%20method%2C%20momentum%3DFalse%2C%20ord%3Dnp.inf%2C%20**verifier_kwargs)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Solve%20the%20ISTA%2FFISTA%20verification%20problem%20for%20one%20backend.%0A%0A%20%20%20%20%20%20%20%20Builds%20a%20%60%60Verifier%60%60%20in%20%60%60method%60%60%20mode%2C%20registers%20the%20parameter%20box%0A%20%20%20%20%20%20%20%20%60%60%5Bx_lb%2C%20x_ub%5D%60%60%20and%20the%20initial%20iterate%20box%20%60%60%5Bs0_lb%2C%20s0_ub%5D%60%60%2C%20rolls%0A%20%20%20%20%20%20%20%20out%20%60%60ista_step%60%60%20over%20the%20horizon%20%60%60K%60%60%2C%20and%20returns%20%60%60(vp%2C%20residuals)%60%60%0A%20%20%20%20%20%20%20%20with%20one%20worst-case%20residual%20per%20iteration.%20Extra%20%60%60verifier_kwargs%60%60%0A%20%20%20%20%20%20%20%20(%60%60solver%60%60%2C%20%60%60params%60%60%2C%20...)%20are%20forwarded%20to%20the%20%60%60Verifier%60%60%0A%20%20%20%20%20%20%20%20constructor.%0A%0A%20%20%20%20%20%20%20%20%60%60ord%60%60%20defaults%20to%20%60%60inf%60%60%2C%20which%20keeps%20the%20MIP%20a%20MILP.%20A%202-norm%0A%20%20%20%20%20%20%20%20objective%20is%20squared%2C%20so%20its%20value%20is%20square-rooted%20on%20return.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20vp%20%3D%20ft.Verifier(method%3Dmethod%2C%20**verifier_kwargs)%0A%20%20%20%20%20%20%20%20x_param%20%3D%20vp.parameter(p%2C%20lb%3Dx_lb%2C%20ub%3Dx_ub)%0A%20%20%20%20%20%20%20%20z%20%3D%20vp.iterate_trajectory(K%20%2B%201%2C%20collapse_expr%3D'auto')%0A%20%20%20%20%20%20%20%20z%5B0%5D%20%3D%20vp.iterate(n%2C%20lb%3Ds0_lb%2C%20ub%3Ds0_ub)%0A%20%20%20%20%20%20%20%20y%2C%20t%20%3D%20z%5B0%5D%2C%201.0%0A%20%20%20%20%20%20%20%20residuals%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20k%20in%20range(1%2C%20K%20%2B%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20z%2C%20y%2C%20t%2C%20residual%20%3D%20ista_step(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vp%2C%20D%2C%20eta%2C%20lambd%2C%20z%2C%20y%2C%20x_param%2C%20t%2C%20k%3Dk%2C%20momentum%3Dmomentum%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20maximize_residual(vp%2C%20residual%2C%20ord)%0A%20%20%20%20%20%20%20%20%20%20%20%20val%20%3D%20vp.solve()%0A%20%20%20%20%20%20%20%20%20%20%20%20residuals.append(np.sqrt(val)%20if%20ord%20%3D%3D%202%20else%20val)%0A%20%20%20%20%20%20%20%20return%20vp%2C%20residuals%0A%0A%20%20%20%20return%20(ista_verification%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20MIP-VP%0A%20%20%20%20We%20solve%20the%20exact%20MIP%20verification%20problem%20for%20both%20ISTA%20and%20its%0A%20%20%20%20accelerated%20variant%20FISTA.%20The%20MIP%20searches%20over%20the%20measurement%20parameter%0A%20%20%20%20%24x%24%20and%20initial%20iterate%20%24z%5E0%24%2C%20and%20maximizes%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24%20at%20each%0A%20%20%20%20horizon.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(K%2C%20ista_verification%2C%20np)%3A%0A%20%20%20%20_%2C%20all_residuals%20%3D%20ista_verification(%0A%20%20%20%20%20%20%20%20K%2C%20'mip'%2C%20momentum%3DFalse%2C%20ord%3Dnp.inf%2C%20solver%3D'gurobi'%2C%0A%20%20%20%20)%0A%20%20%20%20_%2C%20all_residuals_fista%20%3D%20ista_verification(%0A%20%20%20%20%20%20%20%20K%2C%20'mip'%2C%20momentum%3DTrue%2C%20ord%3Dnp.inf%2C%20solver%3D'gurobi'%2C%0A%20%20%20%20)%0A%20%20%20%20return%20all_residuals%2C%20all_residuals_fista%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20PEP-obtained%20worst-case%20convergence%20rate%0A%0A%20%20%20%20We%20compare%20our%20ISTA%20results%20to%20the%20worst-case%20performance%20returned%20by%20the%0A%20%20%20%20performance%20estimation%20problem%20(PEP)%2C%20using%20the%20shared%20%60PEPitRunner%60.%20The%0A%20%20%20%20radius%20is%20certified%20rather%20than%20measured%20on%20the%20samples%2C%20so%20the%20curve%20is%20a%0A%20%20%20%20valid%20upper%20bound%20for%20the%20whole%20box%20the%20verification%20problem%20searches%3B%20the%0A%20%20%20%20price%20is%20that%20it%20is%20far%20looser%20than%20the%20sampled%20radius%20used%20previously%0A%20%20%20%20(%24R%20%5Capprox%2023.7%24%20instead%20of%20%243.0%24%2C%20and%20the%20PEP%20scales%20as%20%24R%5E2%24).%0A%0A%20%20%20%20The%20FISTA%20curve%20uses%20a%20PEP%20written%20to%20match%20%60ista_step%60%20rather%20than%20the%0A%20%20%20%20generic%20accelerated%20template.%20%60PEPitRunner%60%20exposes%20this%20through%0A%20%20%20%20%60variant%3D%22fista%22%60%2C%20which%20changes%20two%20things%20relative%20to%20the%20default%0A%20%20%20%20%60variant%3D%22simple%22%60%20used%20by%20the%20NNLS%20tutorial%3A%20the%20extrapolation%20weights%20come%0A%20%20%20%20from%20the%20classical%20sequence%20%24t_0%20%3D%201%24%2C%0A%20%20%20%20%24t_%7Bk%7D%20%3D%20(1%2B%5Csqrt%7B1%2B4t_%7Bk-1%7D%5E2%7D)%2F2%24%2C%20%24%5Cbeta_k%20%3D%20(t_%7Bk-1%7D-1)%2Ft_%7Bk%7D%24%20instead%0A%20%20%20%20of%20the%20hardcoded%20%24%5Cbeta_k%20%3D%20(k-1)%2F(k%2B2)%24%20(they%20differ%20precisely%20over%20the%0A%20%20%20%20small%20%24K%24%20plotted%20here%3A%20%24%5Cbeta_2%20%3D%200.282%24%20vs%20%240.250%24%2C%0A%20%20%20%20%24%5Cbeta_3%20%3D%200.434%24%20vs%20%240.400%24)%2C%20and%20the%20performance%20metric%20is%20the%20difference%0A%20%20%20%20of%20the%20*proximal*%20iterates%20%24y%5Ek%20-%20y%5E%7Bk-1%7D%24%20rather%20than%20of%20the%20extrapolated%0A%20%20%20%20ones%2C%20which%20is%20the%20residual%20%60ista_step%60%20returns.%20Since%20%24y%5E0%20%3D%20z%5E0%24%20in%0A%20%20%20%20the%20notebook%2C%20the%20same%20initial%20condition%20%24%5C%7Cz%5E0-z%5E%5Cstar%5C%7C_2%20%5Cle%20R%24%20covers%0A%20%20%20%20both%20sequences.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(%0A%20%20%20%20K%2C%0A%20%20%20%20L%2C%0A%20%20%20%20certify_lasso_pepit_radius%2C%0A%20%20%20%20eta%2C%0A%20%20%20%20lambd%2C%0A%20%20%20%20n%2C%0A%20%20%20%20np%2C%0A%20%20%20%20s0_lb%2C%0A%20%20%20%20s0_ub%2C%0A%20%20%20%20x_lb%2C%0A%20%20%20%20x_ub%2C%0A)%3A%0A%20%20%20%20from%20docs.tutorials.helper_functions.pepit_runner%20import%20PEPitRunner%0A%0A%20%20%20%20%23%20Certified%20initial-condition%20radius%3A%20covers%20the%20whole%20box%2C%20not%20just%20the%0A%20%20%20%20%23%20drawn%20samples%2C%20so%20the%20PEP%20value%20below%20is%20a%20genuine%20upper%20bound.%0A%20%20%20%20_pepit_radius%20%3D%20certify_lasso_pepit_radius(lambd%2C%20s0_lb%2C%20s0_ub%2C%20x_lb%2C%20x_ub%2C%20n)%0A%0A%20%20%20%20_runner%20%3D%20PEPitRunner()%0A%20%20%20%20_pep%2C%20_%20%3D%20_runner.proximal_gradient_residuals(%0A%20%20%20%20%20%20%20%20K%3DK%2C%20L%3Dfloat(L)%2C%20eta%3Dfloat(eta)%2C%20R%3D_pepit_radius%2C%20nonsmooth%3D%22convex%22%0A%20%20%20%20)%0A%20%20%20%20%23%20variant%3D%22fista%22%20matches%20%60%60ista_step%60%60%3A%20classical%20t_k%20weights%2C%20and%20the%0A%20%20%20%20%23%20metric%20taken%20on%20the%20proximal%20iterates%20y%5Ek%20rather%20than%20the%20extrapolated%20z%5Ek.%0A%20%20%20%20_pep_fista%2C%20_%20%3D%20_runner.accelerated_proximal_gradient_residuals(%0A%20%20%20%20%20%20%20%20K%3DK%2C%20L%3Dfloat(L)%2C%20eta%3Dfloat(eta)%2C%20R%3D_pepit_radius%2C%20nonsmooth%3D%22convex%22%2C%0A%20%20%20%20%20%20%20%20variant%3D%22fista%22%2C%0A%20%20%20%20)%0A%20%20%20%20%23%20PEPit's%20metric%20is%20%7C%7Cx%5Bk%5D%20-%20x%5Bk-1%5D%7C%7C%5E2%2C%20so%20its%20square%20root%20bounds%0A%20%20%20%20%23%20%7C%7Cr%5Ek%7C%7C_2%2C%20and%20hence%20%7C%7Cr%5Ek%7C%7C_inf%2C%20the%20norm%20the%20MIP%20maximizes.%0A%20%20%20%20all_residuals_PEP%20%3D%20np.sqrt(np.asarray(_pep%2C%20dtype%3Dfloat))%0A%20%20%20%20all_residuals_PEP_fista%20%3D%20np.sqrt(np.asarray(_pep_fista%2C%20dtype%3Dfloat))%0A%20%20%20%20return%20all_residuals_PEP%2C%20all_residuals_PEP_fista%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Numerical%20Result%0A%0A%20%20%20%20Now%20we%20plot%20the%20verification%20results%20against%20the%20PEPit%20bound%20and%20the%0A%20%20%20%20in-sample%20simulation%20statistics%2C%20including%20the%20sample%20worst-case%20instance.%0A%20%20%20%20The%20simulation%20and%20MIP%20curves%20report%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24.%20PEPit%20bounds%0A%20%20%20%20%24%5C%7Cr%5Ek%5C%7C_2%24%2C%20which%20is%20also%20a%20bound%20on%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24%20because%0A%20%20%20%20%24%5C%7Cv%5C%7C_%5Cinfty%20%5Cle%20%5C%7Cv%5C%7C_2%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20K%2C%0A%20%20%20%20all_residuals%2C%0A%20%20%20%20all_residuals_PEP%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20simulation_residual_quantiles_ista%2C%0A%20%20%20%20simulation_residuals_ista%2C%0A)%3A%0A%20%20%20%20plt.figure()%0A%20%20%20%20_ks%20%3D%20range(1%2C%20K%2B1)%0A%20%20%20%20plt.plot(_ks%2C%20all_residuals_PEP%2C%20label%3Dr'ISTA%20(PEPit%2C%20bounds%20%24%5C%7Cr%5Ek%5C%7C_2%20%5Cgeq%20%5C%7Cr%5Ek%5C%7C_%5Cinfty%24)')%0A%20%20%20%20plt.plot(_ks%2C%20all_residuals%2C%20label%3D'ISTA%20(MIP%20worst-case)')%0A%20%20%20%20_avg_line%2C%20%3D%20plt.plot(_ks%2C%20simulation_residuals_ista%2C%20'--'%2C%20label%3D'ISTA%20(in-sample%20avg)')%0A%20%20%20%20plt.fill_between(_ks%2C%20simulation_residual_quantiles_ista%5B2%5D%2C%20simulation_residual_quantiles_ista%5B0%5D%2C%20alpha%3D.18%2C%20label%3D'ISTA%20sample%2010%25-90%25%20band')%0A%20%20%20%20plt.plot(_ks%2C%20simulation_residual_quantiles_ista%5B3%5D%2C%20'%3A'%2C%20color%3D_avg_line.get_color()%2C%20label%3D'ISTA%20sample%20worst-case')%0A%20%20%20%20%23%20Bound%20the%20axis%20by%20the%20plotted%20*lines*%2C%20not%20by%20the%2010%25-90%25%20band%3A%20the%20band's%0A%20%20%20%20%23%20lower%20edge%20decays%20to%20~0%2C%20which%20would%20otherwise%20squeeze%20every%20curve%20into%20the%0A%20%20%20%20%23%20top%20of%20the%20figure.%20Nothing%20drawn%20as%20a%20line%20is%20clipped.%0A%20%20%20%20_lines%20%3D%20np.concatenate(%5B%0A%20%20%20%20%20%20%20%20all_residuals%2C%0A%20%20%20%20%20%20%20%20simulation_residuals_ista%2C%0A%20%20%20%20%20%20%20%20simulation_residual_quantiles_ista%5B3%5D%2C%0A%20%20%20%20%5D)%0A%20%20%20%20plt.ylim(bottom%3D0.5%20*%20float(np.min(_lines)))%0A%20%20%20%20plt.suptitle('sparse%20coding%20example')%0A%20%20%20%20plt.yscale('log')%0A%20%20%20%20plt.xlabel(r'%24K%24')%0A%20%20%20%20plt.ylabel(r'Worst-case%20fixed-point%20residual%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24')%0A%20%20%20%20plt.legend()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20We%20also%20plot%20the%20verification%20results%20for%20FISTA.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20K%2C%0A%20%20%20%20all_residuals_PEP_fista%2C%0A%20%20%20%20all_residuals_fista%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20simulation_residual_quantiles_fista%2C%0A%20%20%20%20simulation_residuals_fista%2C%0A)%3A%0A%20%20%20%20plt.figure()%0A%20%20%20%20_ks%20%3D%20range(1%2C%20K%2B1)%0A%20%20%20%20plt.plot(_ks%2C%20all_residuals_PEP_fista%2C%20label%3Dr'FISTA%20(PEPit%2C%20bounds%20%24%5C%7Cr%5Ek%5C%7C_2%20%5Cgeq%20%5C%7Cr%5Ek%5C%7C_%5Cinfty%24)')%0A%20%20%20%20plt.plot(_ks%2C%20all_residuals_fista%2C%20label%3D'FISTA%20(MIP%20worst-case)')%0A%20%20%20%20_avg_line%2C%20%3D%20plt.plot(_ks%2C%20simulation_residuals_fista%2C%20'--'%2C%20label%3D'FISTA%20(in-sample%20avg)')%0A%20%20%20%20plt.fill_between(_ks%2C%20simulation_residual_quantiles_fista%5B2%5D%2C%20simulation_residual_quantiles_fista%5B0%5D%2C%20alpha%3D.18%2C%20label%3D'FISTA%20sample%2010%25-90%25%20band')%0A%20%20%20%20plt.plot(_ks%2C%20simulation_residual_quantiles_fista%5B3%5D%2C%20'%3A'%2C%20color%3D_avg_line.get_color()%2C%20label%3D'FISTA%20sample%20worst-case')%0A%20%20%20%20%23%20Bound%20the%20axis%20by%20the%20plotted%20*lines*%2C%20not%20by%20the%2010%25-90%25%20band%3A%20the%20band's%0A%20%20%20%20%23%20lower%20edge%20decays%20to%20~0%2C%20which%20would%20otherwise%20squeeze%20every%20curve%20into%20the%0A%20%20%20%20%23%20top%20of%20the%20figure.%20Nothing%20drawn%20as%20a%20line%20is%20clipped.%0A%20%20%20%20_lines%20%3D%20np.concatenate(%5B%0A%20%20%20%20%20%20%20%20all_residuals_fista%2C%0A%20%20%20%20%20%20%20%20simulation_residuals_fista%2C%0A%20%20%20%20%20%20%20%20simulation_residual_quantiles_fista%5B3%5D%2C%0A%20%20%20%20%5D)%0A%20%20%20%20plt.ylim(bottom%3D0.5%20*%20float(np.min(_lines)))%0A%20%20%20%20plt.suptitle('sparse%20coding%20FISTA%20example')%0A%20%20%20%20plt.yscale('log')%0A%20%20%20%20plt.xlabel(r'%24K%24')%0A%20%20%20%20plt.ylabel(r'Worst-case%20fixed-point%20residual%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24')%0A%20%20%20%20plt.legend()%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
1cd69754327e9dd9c62dbd0174a36a97