nnls

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.nnls%20import%20(%0A%20%20%20%20%20%20%20%20certify_nnls_pepit_radius%2C%0A%20%20%20%20%20%20%20%20make_nnls_matrix%2C%0A%20%20%20%20%20%20%20%20pgd_step%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_nnls_pepit_radius%2C%20make_nnls_matrix%2C%20mo%2C%20np%2C%20pgd_step%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%20Nonnegative%20Least-squares%20(NNLS)%0A%0A%20%20%20%20%5BBack%20to%20the%20verification%20tutorials%5D(.%2F)%0A%0A%20%20%20%20In%20this%20tutorial%2C%20we%20show%20how%20to%20use%20our%20package%20to%20verify%20the%20performance%20of%20first-order%20methods%20solving%20the%20nonnegative%20least-squares%20(NNLS)%20problem%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%7CAz%20-%20x%5C%7C_2%5E2%20%5C%5C%20%5Ctext%7Bsubject%20to%7D%20%26%20z%20%5Cge%200%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20where%20%24A%20%5Cin%20%5Cmathbf%7BR%7D%5E%7Bm%20%5Ctimes%20n%7D%24%20with%20%24m%20%5Cge%20n%24%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%5Em%24%20is%20the%20problem%20parameter.%0A%0A%20%20%20%20%23%23%20Problem%20parameters%0A%0A%20%20%20%20In%20this%20example%2C%20we%20choose%20%24m%20%3D%205%2C%20n%20%3D%203%24%2C%20randomly%20sample%20%24A%24%2C%20and%20parameterize%20%24x%24%20by%20the%20box%20%24%5Cmathcal%7BX%7D%20%3D%20%5B10%2C%2011%5D%5Em%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(make_nnls_matrix)%3A%0A%20%20%20%20m%2C%20n%20%3D%205%2C%203%0A%20%20%20%20A%20%3D%20make_nnls_matrix(m%3Dm%2C%20n%3Dn)%0A%20%20%20%20x_lb%2C%20x_ub%20%3D%2010.0%2C%2011.0%20%20%20%20%23%20parameter%20set%0A%20%20%20%20return%20A%2C%20m%2C%20n%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%20%23%23%20Performance%20metric%0A%0A%20%20%20%20Given%20the%20sequence%20of%20iterates%20%24%5C%7Bz%5Ek%5C%7D_%7Bk%3D0%7D%5EK%24%20generated%20from%20the%20algorithm%20of%20our%20choice%2C%20we%20quantify%20its%20performance%20through%20measuring%20how%20fast%20the%20fixed-point%20residual%20%24r%5EK%20%3D%20z%5EK%20-%20z%5E%7BK-1%7D%24%20decreases%2C%20given%20the%20iteration%20budget%20%24K%3E0%24.%0A%20%20%20%20We%20impose%20the%20initial%20condition%20%24z%5E0%20%5Cin%20%5Cmathcal%7BZ%7D%20%3A%3D%20%5Bz%5E0_%5Cmathrm%7Blb%7D%2C%5C%2C%20z%5E0_%5Cmathrm%7Bub%7D%5D%24%20(%60z0_lb%60%2C%20%60z0_ub%60%20in%20the%20code).%0A%0A%20%20%20%20This%20amounts%20to%20finding%20the%20worst-case%20performance%20by%20solving%20the%20verification%20problem%20(VP)%20of%20the%20following%20form%3A%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%20%20%20%20%5Ctext%7Bmaximize%7D%20%26%20%5C%7Cr%5EK%5C%7C_%5Cinfty%20%5C%5C%0A%20%20%20%20%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%26%20%5C%7Bz%5Ek%5C%7D_%7Bk%3D0%7D%5EK%20%5Ctext%7B%20given%20by%20algorithm%7D%20%5C%5C%0A%20%20%20%20%20%20%20%20%26%20z%5E0%20%5Cin%20%5Cmathcal%7BZ%7D%20%5C%5C%0A%20%20%20%20%20%20%20%20%26%20x%20%5Cin%20%5Cmathcal%7BX%7D%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20where%20%24%5Cmathcal%7BX%7D%20%3A%3D%20%5Bx_%5Cmathrm%7Blb%7D%2C%5C%2C%20x_%5Cmathrm%7Bub%7D%5D%24%20(%60x_lb%60%2C%20%60x_ub%60%20in%20the%20code)%20is%20the%20space%20of%20problem%20parameters.%0A%20%20%20%20As%20in%20%5BExact%20Verification%20of%20First-Order%20Methods%20via%20Mixed-Integer%20Linear%20Programming%5D(https%3A%2F%2Farxiv.org%2Fabs%2F2412.11330)%2C%20we%20measure%20the%20residual%20in%20the%20infinity%20norm%2C%20which%20keeps%20the%20VP%20a%20mixed-integer%20*linear*%20program.%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%20z0_lb%2C%20z0_ub%20%3D%204.0%2C%206.0%20%20%20%20%23%20initial%20condition%0A%20%20%20%20return%20z0_lb%2C%20z0_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%20%23%23%20First-order%20method%3A%20Projected%20gradient%20descent%0A%0A%20%20%20%20We%20solve%20this%20problem%20using%20projected%20gradient%20descent%20with%20step%20size%20%24%5Ceta%20%3E%200%24%2C%20whose%20iterations%20are%20given%20by%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%5CPi_%7B%5Cmathbf%7BR%7D%5En_%2B%7D(z%5E%7Bk%7D%20-%20%5Ceta%20A%5ET(Az%5E%7Bk%7D%20-%20x))%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20where%20%24%5CPi_%7B%5Cmathbf%7BR%7D%5En_%2B%7D(v)%20%3D%20%5Cmax%5C%7Bv%2C%200%5C%7D%24%2C%20taken%20element-wise%2C%20is%20the%20projection%20onto%20the%20nonnegative%20orthant%20(the%20%60relu%60%20step).%0A%20%20%20%20We%20choose%20%24%5Ceta%20%3D%201%2F(2L)%24%2C%20where%20%24L%24%20is%20the%20maximum%20eigenvalue%20of%20%24A%5ET%20A%24.%0A%0A%20%20%20%20We%20also%20consider%20a%20momentum%20variant%20of%20projected%20gradient%20descent%20with%20Nesterov-type%20extrapolation%3A%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Brcl%7D%0A%20%20%20%20%20%20%20%20w%5E%7Bk%2B1%7D%20%26%3D%26%20%5CPi_%7B%5Cmathbf%7BR%7D%5En_%2B%7D%5Cleft(z%5E%7Bk%7D%20-%20%5Ceta%20A%5ET(Az%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%20w%5E%7Bk%2B1%7D%20%2B%20%5Cbeta_%7Bk%2B1%7D%20(w%5E%7Bk%2B1%7D%20-%20w%5E%7Bk%7D)%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20with%20%24w%5E0%20%3D%20z%5E0%24%20and%20%24%5Cbeta_k%20%3D%20(k-1)%2F(k%2B2)%24%2C%20so%20that%20%24%5Cbeta_1%20%3D%200%24.%0A%0A%20%20%20%20The%20shared%20implementation%20used%20throughout%20this%20tutorial%20is%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20def%20pgd_step(vp%2C%20A%2C%20eta%2C%20z_prev%2C%20w_prev%2C%20x_param%2C%20k%2C%20momentum%3DFalse)%3A%0A%20%20%20%20%20%20%20%20w_next%20%3D%20vp.relu(%0A%20%20%20%20%20%20%20%20%20%20%20%20z_prev%20-%20eta%20*%20A.T%20%40%20(A%20%40%20z_prev%20-%20x_param)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20if%20momentum%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20beta_k%20%3D%20(k%20-%201)%20%2F%20(k%20%2B%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20z_next%20%3D%20w_next%20%2B%20beta_k%20*%20(w_next%20-%20w_prev)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20z_next%20%3D%20w_next%0A%20%20%20%20%20%20%20%20residual%20%3D%20z_next%20-%20z_prev%0A%20%20%20%20%20%20%20%20return%20z_next%2C%20w_next%2C%20residual%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20Thus%2C%20%60momentum%3DFalse%60%20gives%20ordinary%20projected%20gradient%20descent%2C%20while%0A%20%20%20%20%60momentum%3DTrue%60%20applies%20the%20momentum%20extrapolation%20before%20the%20next%20iteration.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(A%2C%20np)%3A%0A%20%20%20%20L%20%3D%20np.max(np.linalg.eigvalsh(A.T%20%40%20A))%20%20%23%20symmetric%20routine%20%E2%86%92%20real%20eigenvalues%0A%20%20%20%20eta%20%3D%200.5%20%2F%20L%0A%20%20%20%20return%20L%2C%20eta%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%60pgd_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%20momentum%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%20We%20now%20describe%20how%20we%20can%20use%20our%20package%20%60fomtools%60%20to%20verify%20the%20worst-case%20convergence%20rate%20of%20projected%20gradient%20descent%20and%20its%20accelerated%20variant%2C%20applied%20on%20nonnegative%20least-squares%20problem.%20We%20compare%20the%20simulated%20results%20obtained%20from%20randomly%20generated%20problem%20instances%20to%20the%20obtained%20worst-case%20bound%20from%20PEPit%20and%20our%20approach%2C%20up%20to%20%24K%24%20iterations.%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%2025%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%20The%20%60Verifier%60%20class%20offers%20functionality%20of%20running%20simulations%20of%20the%20algorithm%20given%20randomly%20generated%20problem%20instances%2C%20with%20the%20option%20%60method%3D'simulator'%60.%0A%20%20%20%20This%20randomly-generated%20problem%20instance%20consists%20of%20two%20parameters%3A%20randomly%20generated%20problem%20parameter%20%24x%20%5Cin%20%5Cmathcal%7BX%7D%24%20and%20initial%20iterate%20%24z%5E0%20%5Cin%20%5Cmathcal%7BZ%7D%24.%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_(m%2C%20n%2C%20np%2C%20x_lb%2C%20x_ub%2C%20z0_lb%2C%20z0_ub)%3A%0A%20%20%20%20n_samples%20%3D%2010000%20%20%23%20the%20number%20of%20samples%20to%20simulate%0A%0A%20%20%20%20%23%20Generate%20the%20sample%20of%20problem%20parameter%20x%20and%20initial%20iterate%20z0%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%20m))%0A%20%20%20%20z0_samples%20%3D%20_rng.uniform(low%3Dz0_lb%2C%20high%3Dz0_ub%2C%20size%3D(n_samples%2C%20n))%0A%20%20%20%20return%20x_samples%2C%20z0_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%20Meanwhile%2C%20the%20simulation%20will%20be%20run%20through%20following%20function.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(A%2C%20eta%2C%20m%2C%20n%2C%20np%2C%20pgd_step)%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%20pgd_simulation(x_samples%2C%20z0_samples%2C%20K%2C%20momentum%3DFalse%2C%20ord%3Dnp.inf)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Simulate%20projected%20gradient%20descent%60%60.%0A%0A%20%20%20%20%20%20%20%20Builds%20a%20simulator-mode%20%3Aclass%3A%60Verifier%60%2C%20registers%20the%20parameter%20%60%60x%60%60%20and%0A%20%20%20%20%20%20%20%20initial%20iterate%20%60%60z0%60%60%20as%20sampled%20parameters%2C%20rolls%20out%20the%20recurrence%20with%0A%20%20%20%20%20%20%20%20%60%60pgd_step%60%60%2C%20and%20simulates%20via%20%60%60vp.solve()%60%60.%20Each%20%60%60*_samples%60%60%20array%0A%20%20%20%20%20%20%20%20has%20shape%20%60%60(n_samples%2C%20dim)%60%60%20(a%20single-row%20array%20replays%20one%20instance)%2C%0A%20%20%20%20%20%20%20%20and%20the%20returned%20%60%60residuals%60%60%20has%20shape%20%60%60(n_samples%2C%20K)%60%60.%20Summary%20statistics%0A%20%20%20%20%20%20%20%20(%60%60mean%60%60%2C%20%60%60min%60%60%2C%20%60%60max%60%60%2C%20...)%20are%20available%20on%20the%20returned%20instance%0A%20%20%20%20%20%20%20%20%60%60vp%60%60%20via%20%60%60vp%5B...%5D%60%60.%0A%0A%20%20%20%20%20%20%20%20%60%60ord%60%60%20defaults%20to%20%60%60inf%60%60%2C%20matching%20%3Afunc%3A%60pgd_verification%60.%20With%0A%20%20%20%20%20%20%20%20%60%60ord%3D2%60%60%20the%20objective%20is%20%60%60ft.sum_squares%60%60%2C%20the%20*squared*%20Euclidean%20norm%2C%0A%20%20%20%20%20%20%20%20so%20the%20returned%20residuals%20are%20in%20squared%20units.%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%3Dm%2C%20samples%3Dx_samples)%0A%20%20%20%20%20%20%20%20z0%20%3D%20vp.parameter(ndim%3Dn%2C%20samples%3Dz0_samples)%0A%0A%20%20%20%20%20%20%20%20z%2C%20w%20%3D%20z0%2C%20z0%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%20w%2C%20residual%20%3D%20pgd_step(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vp%2C%20A%2C%20eta%2C%20z%2C%20w%2C%20x_param%2C%20k%3Dk%2C%20momentum%3Dmomentum%2C%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%20maximize_residual%2C%20pgd_simulation%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%20Given%20the%20parameter%20sample%20%60x_samples%60%20and%20the%20initial%20iterate%20sample%20%60z0_samples%60%2C%20simulator%20mode%20runs%20the%20algorithm%20on%20each%20problem%20instance%20%60x_samples%5Bi%5D%60%20starting%20from%20%60z0_samples%5Bi%5D%60.%0A%20%20%20%20Running%20the%20solve%20method%20%60_vp.solve()%60%20outputs%20the%20history%20of%20residuals%20at%20objective%20(fixed-point%20residual).%20After%20calling%20%60solve()%60%2C%20summary%20statistics%20are%20also%20available%20directly%20on%20%60_vp%60%2C%20including%20%60mean%60%2C%20%60std%60%2C%20%60min%60%2C%20%60max%60%2C%20%60median%60%2C%20%60empirical_risk%60%2C%20%60risk_bound%60%2C%20%60n_samples%60%2C%20%60n_iterations%60%2C%20%60runtime%60%2C%20%60time_per_sample%60%2C%20and%20%60throughput%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(K%2C%20np%2C%20pgd_simulation%2C%20x_samples%2C%20z0_samples)%3A%0A%20%20%20%20_vp%2C%20_residuals%20%3D%20pgd_simulation(x_samples%2C%20z0_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%20%23%20Take%20average%20of%20residuals%20for%20iterate%20K%3D1%20to%20K%3DK%0A%20%20%20%20simulation_residuals%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%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%2C%20simulation_residuals%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%20run%20the%20simulation%20of%20momentum-accelerated%20projected%20gradient%20descent%20as%20well%2C%20reusing%20the%20same%20%60pgd_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%20np%2C%20pgd_simulation%2C%20x_samples%2C%20z0_samples)%3A%0A%20%20%20%20_%2C%20_residuals%20%3D%20pgd_simulation(x_samples%2C%20z0_samples%2C%20K%2C%20momentum%3DTrue%2C%20ord%3Dnp.inf)%0A%20%20%20%20simulation_residuals_momentum%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_momentum%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%20(%0A%20%20%20%20%20%20%20%20simulation_residual_quantiles_momentum%2C%0A%20%20%20%20%20%20%20%20simulation_residuals_momentum%2C%0A%20%20%20%20)%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%20guarantee%20by%20solving%20the%20corresponding%20verification%20problem%20(VP).%0A%20%20%20%20Passing%20the%20option%20%60method%3D'mip'%60%20formulates%20the%20VP%20as%20a%20mixed-integer%20linear%20program%20(MILP)%2C%20based%20on%20%5BExact%20Verification%20of%20First-Order%20Methods%20via%20Mixed-Integer%20Linear%20Programming.%5D(https%3A%2F%2Farxiv.org%2Fabs%2F2412.11330)%0A%0A%0A%20%20%20%20The%20following%20function%20rolls%20out%20the%20same%20%60pgd_step%60%20recurrence%20over%20the%0A%20%20%20%20horizon%20%24K%24%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%20A%2C%0A%20%20%20%20eta%2C%0A%20%20%20%20ft%2C%0A%20%20%20%20m%2C%0A%20%20%20%20maximize_residual%2C%0A%20%20%20%20n%2C%0A%20%20%20%20np%2C%0A%20%20%20%20pgd_step%2C%0A%20%20%20%20x_lb%2C%0A%20%20%20%20x_ub%2C%0A%20%20%20%20z0_lb%2C%0A%20%20%20%20z0_ub%2C%0A)%3A%0A%20%20%20%20def%20pgd_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%20projected-gradient%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%5Bz0_lb%2C%20z0_ub%5D%60%60%2C%0A%20%20%20%20%20%20%20%20rolls%20out%20%60%60pgd_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%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(m%2C%20lb%3Dx_lb%2C%20ub%3Dx_ub)%0A%20%20%20%20%20%20%20%20z%20%3D%20vp.iterate_trajectory(K%2B1%2C%20collapse_expr%3D'auto')%0A%20%20%20%20%20%20%20%20z%5B0%5D%20%3D%20vp.iterate(n%2C%20lb%3Dz0_lb%2C%20ub%3Dz0_ub)%0A%20%20%20%20%20%20%20%20w%20%3D%20z%5B0%5D%0A%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%20w%2C%20residual%20%3D%20pgd_step(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vp%2C%20A%2C%20eta%2C%20z%2C%20w%2C%20x_param%2C%20k%3Dk%2C%20momentum%3Dmomentum%2C%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(pgd_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%0A%20%20%20%20The%20following%20example%20formulates%20the%20verification%20problem%20as%20a%20mixed-integer%20linear%20program%20(MILP)%20and%20solves%20it%20to%20obtain%20exact%20worst-case%20fixed-point%20residual%20after%20%24K%24%20iterations%2C%20for%20both%20projected%20gradient%20descent%20and%20its%20accelerated%20variant.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(K%2C%20np%2C%20pgd_verification)%3A%0A%20%20%20%20_%2C%20all_residuals%20%3D%20pgd_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_momentum%20%3D%20pgd_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_momentum%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%20As%20we%20compare%20our%20results%20to%20the%20worst-case%20performance%20of%20projected%20gradient%20descent%2C%20we%20use%20performance%20estimation%20problem%20(PEP)%20to%20obtain%20such%20worst-case%20rate%20using%20the%20shared%20%60PEPitRunner%60.%0A%0A%20%20%20%20%23%23%23%23%20Certified%20initial-condition%20radius%0A%0A%20%20%20%20PEP%20bounds%20the%20worst%20case%20over%20every%20%24L%24-smooth%20convex%20%24f%24%20and%20convex%20indicators%20%24g%24%20whose%20initial%20iterate%20lies%20in%20a%20Euclidean%20ball%20of%20radius%20%24R%24%20about%20the%20fixed%20point%3A%20%24%5C%7Cz%5E0%20-%20z%5E%5Cstar%5C%7C_2%20%5Cle%20R%24.%0A%20%20%20%20For%20this%20PEP%20bound%20to%20be%20a%20valid%20upper%20bound%20on%20our%20VP%20with%20%24z%5E0%20%5Cin%20%5Cmathcal%7BZ%7D%24%2C%20%24x%20%5Cin%20%5Cmathcal%7BX%7D%24%2C%20the%20ball%20must%20cover%20the%20box%3A%0A%0A%20%20%20%20%24%24R%20%5C%3B%5Cge%5C%3B%20%5Csup_%7Bz%5E0%20%5Cin%20%5Cmathcal%7BZ%7D%2C%5C%2C%20x%20%5Cin%20%5Cmathcal%7BX%7D%7D%20%5C%7Cz%5E0%20-%20z%5E%5Cstar(x)%5C%7C_2%20.%24%24%0A%0A%20%20%20%20We%20find%20%24R%24%20via%20solving%20optimization%20problem%20which%20gives%20maximum%20value%20of%20%24%5C%7Cz%5E0%20-%20z%5E%5Cstar(x)%5C%7C_2%24%20for%20all%20%24z%5E0%20%5Cin%20%5Cmathcal%7BZ%7D%24%20and%20%24x%20%5Cin%20%5Cmathcal%7BX%7D%24.%0A%20%20%20%20Note%20that%20the%20NNLS%20solution%20map%20%24x%20%5Cmapsto%20z%5E%5Cstar(x)%24%20is%20piecewise%20linear%3B%20on%20each%20active%20set%20%24%5Cmathcal%7BI%7D%24%2C%20its%20Jacobian%20is%20%24(A_%7B%5Cmathcal%7BI%7D%7D%5ET%20A_%7B%5Cmathcal%7BI%7D%7D)%5E%7B-1%7D%20A_%7B%5Cmathcal%7BI%7D%7D%5ET%24%2C%20so%20%24%5Csigma_%7B%5Cmin%7D(A_%7B%5Cmathcal%7BI%7D%7D)%20%5Cge%20%5Csigma_%7B%5Cmin%7D(A)%24.%0A%20%20%20%20This%20mapping%20is%20globally%20%24%5Cgamma%24-Lipschitz%20with%20%24%5Cgamma%20%3D%201%2F%5Csigma_%7B%5Cmin%7D(A)%24.%0A%20%20%20%20Splitting%20the%20two%20box%20contributions%20with%20the%20triangle%20inequality%20gives%20the%20certified%20radius%0A%0A%20%20%20%20%24%24R%20%3D%20%5Cmax_%7Bz%5E0%20%5Cin%20%5Cmathcal%7BZ%7D%7D%20%5C%7Cz%5E0%20-%20z%5E%5Cstar(%5Cbar%7Bx%7D)%5C%7C_2%0A%20%20%20%20%20%20%20%20%5C%3B%2B%5C%3B%20%5Cgamma%20%5Ccdot%20%5Cmathbf%7Brad%7D(%5Cmathcal%7BX%7D)%2C%24%24%0A%0A%20%20%20%20where%20%24%5Cbar%7Bx%7D%24%20is%20the%20center%20of%20the%20%24x%24-box%20and%20%24%5Cmathbf%7Brad%7D(%5Cmathcal%7BX%7D)%20%3D%20%5Ctfrac12%5Csqrt%7B%5Csum_j%20(x_%5Cmathrm%7Bub%7D-x_%5Cmathrm%7Blb%7D)_j%5E2%7D%24%20is%20its%20circumradius.%0A%20%20%20%20The%20first%20term%20is%20the%20distance%20from%20%24z%5E%5Cstar(%5Cbar%7Bx%7D)%24%20to%20the%20farthest%20box%20vertex%2C%20evaluated%20coordinatewise.%0A%20%20%20%20This%20requires%20a%20single%20NNLS%20solve%20(at%20%24%5Cbar%7Bx%7D%24)%20and%20no%20sampling.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(A%2C%20K%2C%20L%2C%20certify_nnls_pepit_radius%2C%20eta%2C%20np%2C%20x_lb%2C%20x_ub%2C%20z0_lb%2C%20z0_ub)%3A%0A%20%20%20%20from%20docs.tutorials.helper_functions.pepit_runner%20import%20PEPitRunner%0A%0A%20%20%20%20_pepit_radius%20%3D%20certify_nnls_pepit_radius(%0A%20%20%20%20%20%20%20%20A%2C%0A%20%20%20%20%20%20%20%20x_lb%2C%0A%20%20%20%20%20%20%20%20x_ub%2C%0A%20%20%20%20%20%20%20%20z0_lb%2C%0A%20%20%20%20%20%20%20%20z0_ub%2C%0A%20%20%20%20)%0A%20%20%20%20%23%20PEPitRunner%20defaults%20to%20Clarabel%3A%20SCS%20fails%20on%20the%20k%3D1%20PEP%20here%20and%0A%20%20%20%20%23%20returns%20None%2C%20which%20used%20to%20show%20up%20as%20a%20missing%20K%3D1%20point%20in%20the%20figure.%0A%20%20%20%20_runner%20%3D%20PEPitRunner()%0A%0A%20%20%20%20%23%20PEPit's%20performance%20metric%20is%20%7C%7Cx%5Bk%5D%20-%20x%5Bk-1%5D%7C%7C%5E2%2C%20so%20its%20square%20root%0A%20%20%20%20%23%20bounds%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%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%22indicator%22%0A%20%20%20%20)%0A%20%20%20%20_pep_momentum%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%22indicator%22%0A%20%20%20%20)%0A%20%20%20%20all_residuals_PEP%20%3D%20np.sqrt(np.asarray(_pep%2C%20dtype%3Dfloat))%0A%20%20%20%20all_residuals_PEP_momentum%20%3D%20np.sqrt(np.asarray(_pep_momentum%2C%20dtype%3Dfloat))%0A%20%20%20%20return%20all_residuals_PEP%2C%20all_residuals_PEP_momentum%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%20Numerical%20result%3A%20comparison%20against%20baselines%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%0A%20%20%20%20worst-case%20instance.%20The%20simulation%20and%20MIP%20curves%20report%0A%20%20%20%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24.%20PEPit%20bounds%20%24%5C%7Cr%5Ek%5C%7C_2%24%2C%20which%20is%20also%20a%20bound%20on%0A%20%20%20%20%24%5C%7Cr%5Ek%5C%7C_%5Cinfty%24%20because%20%24%5C%7Cv%5C%7C_%5Cinfty%20%5Cle%20%5C%7Cv%5C%7C_2%24.%20The%20curves%20answer%0A%20%20%20%20different%20questions%2C%20so%20they%20are%20not%20interchangeable%3A%0A%0A%20%20%20%20-%20**sample%20worst-case**%20%E2%80%94%20a%20*lower*%20bound%20on%20the%20true%20worst%20case.%20The%0A%20%20%20%20%20%20samples%20are%20uniform%20over%20the%20same%20boxes%20the%20verifiers%20search%2C%20but%20uniform%0A%20%20%20%20%20%20sampling%20in%20%248%24%20dimensions%20concentrates%20away%20from%20the%20vertices%20and%20faces%0A%20%20%20%20%20%20where%20the%20maximizer%20sits%2C%20so%20this%20curve%20is%20loose%20by%20construction.%0A%20%20%20%20-%20**MIP-VP**%20%E2%80%94%20the%20*exact*%20worst%20case%20over%20the%20boxes%2C%20for%20this%20problem%0A%20%20%20%20%20%20instance.%0A%20%20%20%20-%20**PEPit**%20%E2%80%94%20the%20worst%20case%20of%20%24%5C%7Cr%5Ek%5C%7C_2%24%20over%20*every*%20%24L%24-smooth%20convex%0A%20%20%20%20%20%20%24f%24%20and%20convex%20indicator%20%24g%24%20with%20%24%5C%7Cz%5E0-z%5E%5Cstar%5C%7C_2%20%5Cle%20R%24.%20It%20bounds%20a%0A%20%20%20%20%20%20whole%20function%20class%20rather%20than%20this%20instance%2C%20and%20a%20larger%20norm%2C%20so%20it%0A%20%20%20%20%20%20should%20sit%20above%20MIP-VP%3B%20here%20%24R%24%0A%20%20%20%20%20%20comes%20from%20%60certify_nnls_pepit_radius%60%2C%20which%20covers%20the%20entire%0A%20%20%20%20%20%20initial-condition%20box%20rather%20than%20just%20the%20drawn%20samples.%0A%0A%20%20%20%20So%20the%20expected%20ordering%20at%20every%20%24K%24%20is%0A%20%20%20%20%24%5Ctext%7Bsample%20worst-case%7D%20%5Cle%20%5Ctext%7BMIP-VP%7D%20%5Cle%20%5Ctext%7BPEPit%7D%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%20plt%2C%0A%20%20%20%20simulation_residual_quantiles%2C%0A%20%20%20%20simulation_residuals%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'PEPit%20(bounds%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'MIP-VP')%0A%20%20%20%20_avg_line%2C%20%3D%20plt.plot(_ks%2C%20simulation_residuals%2C%20'--'%2C%20label%3D'sample%20average')%0A%20%20%20%20plt.fill_between(_ks%2C%20simulation_residual_quantiles%5B2%5D%2C%20simulation_residual_quantiles%5B0%5D%2C%20alpha%3D.18%2C%20label%3D'sample%2010%25-90%25%20band')%0A%20%20%20%20plt.plot(_ks%2C%20simulation_residual_quantiles%5B3%5D%2C%20'%3A'%2C%20color%3D_avg_line.get_color()%2C%20label%3D'sample%20worst-case')%0A%20%20%20%20plt.yscale('log')%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.xlabel(r'%24K%24')%0A%20%20%20%20plt.suptitle('NNLS%3A%20projected%20gradient%20descent%20(PGD)')%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%20run%20the%20simulation%20of%20momentum-accelerated%20projected%20gradient%20descent%20as%20well%2C%20reusing%20the%20same%20%60pgd_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_(%0A%20%20%20%20K%2C%0A%20%20%20%20all_residuals_PEP_momentum%2C%0A%20%20%20%20all_residuals_momentum%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20simulation_residual_quantiles_momentum%2C%0A%20%20%20%20simulation_residuals_momentum%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_momentum%2C%20label%3Dr'PEPit%20(bounds%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_momentum%2C%20label%3D'MIP-VP')%0A%20%20%20%20_avg_line%2C%20%3D%20plt.plot(_ks%2C%20simulation_residuals_momentum%2C%20'--'%2C%20label%3D'sample%20average')%0A%20%20%20%20plt.fill_between(_ks%2C%20simulation_residual_quantiles_momentum%5B2%5D%2C%20simulation_residual_quantiles_momentum%5B0%5D%2C%20alpha%3D.18%2C%20label%3D'sample%2010%25-90%25%20band')%0A%20%20%20%20plt.plot(_ks%2C%20simulation_residual_quantiles_momentum%5B3%5D%2C%20'%3A'%2C%20color%3D_avg_line.get_color()%2C%20label%3D'sample%20worst-case')%0A%20%20%20%20plt.yscale('log')%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.xlabel(r'%24K%24')%0A%20%20%20%20plt.suptitle('NNLS%3A%20momentum%20PGD')%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%20Finally%2C%20we%20overlay%20PGD%20and%20momentum%20PGD%20to%20compare%20the%20two%20methods%0A%20%20%20%20directly.%20Methods%20are%20distinguished%20by%20color%2C%20and%20the%20data%20series%20(MIP-VP%0A%20%20%20%20and%20sample%20average%20with%20its%2010%25-90%25%20band)%20by%20line%20style.%20PEPit%20is%0A%20%20%20%20left%20out%20here%3A%20it%20bounds%20a%20function%20class%20rather%20than%20this%20instance%2C%20so%20it%0A%20%20%20%20is%20the%20one%20curve%20that%20does%20not%20distinguish%20the%20two%20methods%20*on%20this%0A%20%20%20%20problem*.%20It%20is%20shown%20in%20the%20per-method%20figures%20above.%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%20all_residuals%2C%0A%20%20%20%20all_residuals_momentum%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20simulation_residual_quantiles%2C%0A%20%20%20%20simulation_residual_quantiles_momentum%2C%0A%20%20%20%20simulation_residuals%2C%0A%20%20%20%20simulation_residuals_momentum%2C%0A)%3A%0A%20%20%20%20from%20matplotlib.lines%20import%20Line2D%0A%0A%20%20%20%20plt.figure()%0A%20%20%20%20_ks%20%3D%20range(1%2C%20K%2B1)%0A%0A%20%20%20%20_methods%20%3D%20%5B%0A%20%20%20%20%20%20%20%20('PGD'%2C%20'C0'%2C%20all_residuals%2C%0A%20%20%20%20%20%20%20%20%20simulation_residuals%2C%20simulation_residual_quantiles)%2C%0A%20%20%20%20%20%20%20%20('momentum%20PGD'%2C%20'C1'%2C%20all_residuals_momentum%2C%0A%20%20%20%20%20%20%20%20%20simulation_residuals_momentum%2C%20simulation_residual_quantiles_momentum)%2C%0A%20%20%20%20%5D%0A%20%20%20%20for%20_%2C%20_color%2C%20_mip%2C%20_avg%2C%20_quant%20in%20_methods%3A%0A%20%20%20%20%20%20%20%20plt.plot(_ks%2C%20_mip%2C%20'--'%2C%20color%3D_color)%0A%20%20%20%20%20%20%20%20plt.plot(_ks%2C%20_avg%2C%20'%3A'%2C%20color%3D_color)%0A%20%20%20%20%20%20%20%20plt.fill_between(_ks%2C%20_quant%5B2%5D%2C%20_quant%5B0%5D%2C%20color%3D_color%2C%20alpha%3D.18)%0A%0A%20%20%20%20_method_handles%20%3D%20%5B%0A%20%20%20%20%20%20%20%20Line2D(%5B%5D%2C%20%5B%5D%2C%20color%3D_color%2C%20label%3D_name)%0A%20%20%20%20%20%20%20%20for%20_name%2C%20_color%2C%20*_%20in%20_methods%0A%20%20%20%20%5D%0A%20%20%20%20_style_handles%20%3D%20%5B%0A%20%20%20%20%20%20%20%20Line2D(%5B%5D%2C%20%5B%5D%2C%20color%3D'gray'%2C%20linestyle%3D'--'%2C%20label%3D'MIP-VP')%2C%0A%20%20%20%20%20%20%20%20Line2D(%5B%5D%2C%20%5B%5D%2C%20color%3D'gray'%2C%20linestyle%3D'%3A'%2C%20label%3D'sample%20average')%2C%0A%20%20%20%20%5D%0A%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%20_leg%20%3D%20plt.legend(handles%3D_method_handles%2C%20loc%3D'upper%20right')%0A%20%20%20%20plt.gca().add_artist(_leg)%0A%20%20%20%20plt.legend(handles%3D_style_handles%2C%20loc%3D'lower%20left')%0A%20%20%20%20plt.suptitle('NNLS%3A%20PGD%20vs.%20momentum%20PGD')%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