import%20marimo%0A%0A__generated_with%20%3D%20%220.23.16%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%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%20Choosing%20a%20Ridge%20Penalty%20from%20Validation%20Data%0A%0A%20%20%20%20Regularization%20is%20often%20selected%20by%20trying%20a%20grid%20of%20values.%20Here%20we%20pose%0A%20%20%20%20the%20same%20task%20as%20a%20bilevel%20problem%3A%20the%20**upper%20problem**%20chooses%20the%20ridge%0A%20%20%20%20penalty%20from%20validation%20performance%2C%20while%20the%20**lower%20problem**%20fits%20the%0A%20%20%20%20regression%20coefficients%20for%20that%20penalty%20using%20the%20training%20data.%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%20from%20pathlib%20import%20Path%0A%0A%20%20%20%20import%20cvxpy%20as%20cp%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%20from%20blvpy%20import%20BilevelProblem%2C%20LowerProblem%0A%0A%20%20%20%20plt.style.use(Path(__file__).resolve().parents%5B1%5D%20%2F%20%22_shared%22%20%2F%20%22zhlatex.mplstyle%22)%0A%20%20%20%20return%20BilevelProblem%2C%20LowerProblem%2C%20Path%2C%20cp%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%23%20Bilevel%20formulation%0A%0A%20%20%20%20Let%20%24(X_%7B%5Cmathrm%7Btr%7D%7D%2Cy_%7B%5Cmathrm%7Btr%7D%7D)%24%20and%0A%20%20%20%20%24(X_%7B%5Cmathrm%7Bval%7D%7D%2Cy_%7B%5Cmathrm%7Bval%7D%7D)%24%20denote%20the%20training%20and%20validation%0A%20%20%20%20samples%2C%20with%20%24m_%7B%5Cmathrm%7Btr%7D%7D%24%20and%20%24m_%7B%5Cmathrm%7Bval%7D%7D%24%20observations%2C%0A%20%20%20%20respectively.%20The%20bilevel%20problem%20is%0A%0A%20%20%20%20%5C%5B%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Cmathop%7B%5Cmathrm%7Bminimize%7D%7D_%7B%5Clambda%2Cw%7D%20%26%0A%20%20%20%20%20%20(1%2Fm_%7B%5Cmathrm%7Bval%7D%7D)%0A%20%20%20%20%20%20%5ClVert%20X_%7B%5Cmathrm%7Bval%7D%7Dw-y_%7B%5Cmathrm%7Bval%7D%7D%5CrVert_2%5E2%5C%5C%0A%20%20%20%20%5Cmathop%7B%5Cmathrm%7Bsubject%5C%20to%7D%7D%20%26%0A%20%20%20%20%20%2010%5E%7B-4%7D%5Cleq%5Clambda%5Cleq10%2C%5C%5C%0A%20%20%20%20%20%20%26w%5Cin%20S(%5Clambda)%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%5C%5D%0A%0A%20%20%20%20where%0A%0A%20%20%20%20%5C%5B%0A%20%20%20%20S(%5Clambda)%3D%5Cmathop%7B%5Cmathrm%7Bargmin%7D%7D_w%0A%20%20%20%20%20%20(1%2Fm_%7B%5Cmathrm%7Btr%7D%7D)%0A%20%20%20%20%20%20%5ClVert%20X_%7B%5Cmathrm%7Btr%7D%7Dw-y_%7B%5Cmathrm%7Btr%7D%7D%5CrVert_2%5E2%0A%20%20%20%20%20%20%2B%5Clambda%5ClVert%20w%5CrVert_2%5E2.%0A%20%20%20%20%5C%5D%0A%0A%20%20%20%20The%20lower%20problem%20is%20strongly%20convex%20because%20%24%5Clambda%3E0%24%2C%20so%20its%20response%0A%20%20%20%20is%20unique.%20BLVPY%20converts%20it%20to%20its%20exact%20SOCP%20representation%20and%20applies%0A%20%20%20%20%24%5Cepsilon%24-gap%20continuation%20to%20the%20optimistic%20single-level%20reformulation.%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%20Deterministic%20regression%20data%0A%0A%20%20%20%20We%20generate%208%20features%2C%2024%20noisy%20training%20observations%2C%20and%2080%20validation%0A%20%20%20%20observations.%20The%20validation%20set%20is%20larger%20and%20less%20noisy%2C%20so%20it%20provides%0A%20%20%20%20a%20useful%20signal%20for%20choosing%20a%20nonzero%20regularization%20weight.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20rng%20%3D%20np.random.default_rng(1907)%0A%20%20%20%20n_features%20%3D%208%0A%20%20%20%20m_tr%20%3D%2024%0A%20%20%20%20m_val%20%3D%2080%0A%0A%20%20%20%20w_true%20%3D%20np.array(%5B1.6%2C%20-1.2%2C%200.8%2C%200.0%2C%20-0.5%2C%200.35%2C%200.0%2C%200.6%5D)%0A%20%20%20%20X_tr%20%3D%20rng.normal(size%3D(m_tr%2C%20n_features))%0A%20%20%20%20X_val%20%3D%20rng.normal(size%3D(m_val%2C%20n_features))%0A%20%20%20%20feature_scale%20%3D%20np.std(X_tr%2C%20axis%3D0)%0A%20%20%20%20X_tr%20%3D%20X_tr%20%2F%20feature_scale%0A%20%20%20%20X_val%20%3D%20X_val%20%2F%20feature_scale%0A%20%20%20%20y_tr%20%3D%20X_tr%20%40%20w_true%20%2B%20rng.normal(scale%3D1.1%2C%20size%3Dm_tr)%0A%20%20%20%20y_val%20%3D%20X_val%20%40%20w_true%20%2B%20rng.normal(scale%3D0.25%2C%20size%3Dm_val)%0A%20%20%20%20return%20X_tr%2C%20X_val%2C%20m_tr%2C%20m_val%2C%20n_features%2C%20y_tr%2C%20y_val%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%20Specify%20and%20solve%20the%20bilevel%20model%0A%0A%20%20%20%20Listing%20%60lbd%60%20in%20%60LowerProblem(parameters%3D%5B...%5D)%60%20makes%20it%20fixed%0A%20%20%20%20data%20for%20training%20while%20leaving%20it%20as%20a%20decision%20of%20the%20upper%20problem.%0A%20%20%20%20Note%20that%20the%20native%20bounds%20for%20%60lbd%60%0A%20%20%20%20are%20part%20of%20the%20mathematical%20model.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20BilevelProblem%2C%0A%20%20%20%20LowerProblem%2C%0A%20%20%20%20X_tr%2C%0A%20%20%20%20X_val%2C%0A%20%20%20%20cp%2C%0A%20%20%20%20m_tr%2C%0A%20%20%20%20m_val%2C%0A%20%20%20%20n_features%2C%0A%20%20%20%20y_tr%2C%0A%20%20%20%20y_val%2C%0A)%3A%0A%20%20%20%20lbd%20%3D%20cp.Variable(%0A%20%20%20%20%20%20%20%20nonneg%3DTrue%2C%0A%20%20%20%20%20%20%20%20bounds%3D%5B1e-4%2C%2010.0%5D%2C%0A%20%20%20%20%20%20%20%20name%3D%22lbd%22%2C%0A%20%20%20%20)%0A%20%20%20%20w%20%3D%20cp.Variable(n_features%2C%20name%3D%22w%22)%0A%0A%20%20%20%20training_loss%20%3D%20cp.sum_squares(X_tr%20%40%20w%20-%20y_tr)%20%2F%20m_tr%0A%20%20%20%20lower_problem%20%3D%20LowerProblem(%0A%20%20%20%20%20%20%20%20cp.Minimize(training_loss%20%2B%20lbd%20*%20cp.sum_squares(w))%2C%0A%20%20%20%20%20%20%20%20parameters%3D%5Blbd%5D%2C%0A%20%20%20%20)%0A%20%20%20%20validation_loss%20%3D%20cp.sum_squares(X_val%20%40%20w%20-%20y_val)%20%2F%20m_val%0A%20%20%20%20problem%20%3D%20BilevelProblem(cp.Minimize(validation_loss)%2C%20lower_problem)%0A%20%20%20%20return%20lbd%2C%20problem%0A%0A%0A%40app.cell%0Adef%20_(problem)%3A%0A%20%20%20%20epsilon_target%20%3D%201e-5%0A%20%20%20%20result%20%3D%20problem.solve(%0A%20%20%20%20%20%20%20%20epsilon_initial%3D1e-2%2C%0A%20%20%20%20%20%20%20%20epsilon_target%3Depsilon_target%2C%0A%20%20%20%20%20%20%20%20verbose%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20diagnostics%20%3D%20problem.gap_diagnostics(result)%0A%20%20%20%20return%20diagnostics%2C%20result%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%20Direct-CVXPY%20comparison%20grid%0A%0A%20%20%20%20For%20context%2C%20we%20independently%20solve%20the%20fixed-penalty%20training%20problem%20on%0A%20%20%20%20a%20logarithmic%20grid%20with%20Clarabel.%20Note%20that%20this%20grid%20is%20only%20a%20visualization%0A%20%20%20%20and%20a%20baseline%2C%20i.e.%2C%20it%20is%20not%20visible%20to%20BLVPY%20when%20choosing%20the%20penalty.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(X_tr%2C%20X_val%2C%20cp%2C%20m_tr%2C%20m_val%2C%20n_features%2C%20np%2C%20y_tr%2C%20y_val)%3A%0A%20%20%20%20lbd_grid%20%3D%20np.unique(np.concatenate((np.geomspace(1e-4%2C%2010.0%2C%2070)%2C%20np.array(%5B1.0%5D))))%0A%20%20%20%20lbd_parameter%20%3D%20cp.Parameter(nonneg%3DTrue%2C%20name%3D%22lbd_grid%22)%0A%20%20%20%20w_grid%20%3D%20cp.Variable(n_features%2C%20name%3D%22w_grid%22)%0A%20%20%20%20grid_lower_objective%20%3D%20cp.sum_squares(X_tr%20%40%20w_grid%20-%20y_tr)%20%2F%20m_tr%20%2B%20lbd_parameter%20*%20cp.sum_squares(w_grid)%0A%20%20%20%20grid_problem%20%3D%20cp.Problem(cp.Minimize(grid_lower_objective))%0A%0A%20%20%20%20grid_training_mse%20%3D%20%5B%5D%0A%20%20%20%20grid_validation_mse%20%3D%20%5B%5D%0A%20%20%20%20for%20_lbd%20in%20lbd_grid%3A%0A%20%20%20%20%20%20%20%20lbd_parameter.value%20%3D%20float(_lbd)%0A%20%20%20%20%20%20%20%20grid_problem.solve(solver%3Dcp.CLARABEL%2C%20warm_start%3DTrue)%0A%20%20%20%20%20%20%20%20_w_value%20%3D%20np.asarray(w_grid.value%2C%20dtype%3Dfloat)%0A%20%20%20%20%20%20%20%20grid_training_mse.append(float(np.sum(np.square(X_tr%20%40%20_w_value%20-%20y_tr))%20%2F%20m_tr))%0A%20%20%20%20%20%20%20%20grid_validation_mse.append(float(np.sum(np.square(X_val%20%40%20_w_value%20-%20y_val))%20%2F%20m_val))%0A%0A%20%20%20%20grid_training_mse%20%3D%20np.asarray(grid_training_mse)%0A%20%20%20%20grid_validation_mse%20%3D%20np.asarray(grid_validation_mse)%0A%20%20%20%20baseline_index%20%3D%20int(np.flatnonzero(np.isclose(lbd_grid%2C%201.0))%5B0%5D)%0A%20%20%20%20baseline_validation_mse%20%3D%20float(grid_validation_mse%5Bbaseline_index%5D)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20baseline_validation_mse%2C%0A%20%20%20%20%20%20%20%20grid_training_mse%2C%0A%20%20%20%20%20%20%20%20grid_validation_mse%2C%0A%20%20%20%20%20%20%20%20lbd_grid%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(baseline_validation_mse%2C%20diagnostics%2C%20lbd%2C%20mo%2C%20result)%3A%0A%20%20%20%20assert%20result.succeeded%2C%20result.message%0A%20%20%20%20assert%20float(result.objective)%20%3C%20baseline_validation_mse%20-%200.1%2C%20(%0A%20%20%20%20%20%20%20%20%22The%20selected%20ridge%20weight%20did%20not%20improve%20validation%20error.%22%0A%20%20%20%20)%0A%0A%20%20%20%20mo.md(rf%22%22%22%0A%20%20%20%20%23%23%20Result%0A%0A%20%20%20%20%7C%20quantity%20%7C%20value%20%7C%0A%20%20%20%20%7C%20---%20%7C%20---%3A%20%7C%0A%20%20%20%20%7C%20status%20%7C%20%60%7Bresult.status%7D%60%20%7C%0A%20%20%20%20%7C%20selected%20ridge%20weight%20%7C%20%7Bfloat(lbd.value)%3A.5f%7D%20%7C%0A%20%20%20%20%7C%20validation%20MSE%20at%20%24%5Clambda%3D1%24%20%7C%20%7Bbaseline_validation_mse%3A.6f%7D%20%7C%0A%20%20%20%20%7C%20BLVPY%20validation%20MSE%20%7C%20%7Bfloat(result.objective)%3A.6f%7D%20%7C%0A%20%20%20%20%7C%20final%20epsilon%20%7C%20%7Bresult.final_epsilon%3A.3e%7D%20%7C%0A%20%20%20%20%7C%20maximum%20lifted%20residual%20%7C%20%7Bresult.residuals.max_violation%3A.3e%7D%20%7C%0A%20%20%20%20%7C%20complementarity%20%7C%20%7Bresult.complementarity%3A.3e%7D%20%7C%0A%20%20%20%20%7C%20signed%20source%20gap%20%7C%20%7Bdiagnostics.source_gap%3A.3e%7D%20%7C%0A%0A%20%20%20%20The%20selected%20penalty%20substantially%20improves%20validation%20error%20over%20the%0A%20%20%20%20conventional%20%24%5Clambda%3D1%24%20baseline.%20Its%20nonzero%20value%20also%20illustrates%20the%0A%20%20%20%20bias-variance%20tradeoff%3A%20the%20training%20curve%20worsens%20as%20regularization%20grows%2C%0A%20%20%20%20while%20the%20validation%20curve%20initially%20improves.%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%20Regularization%20path%0A%0A%20%20%20%20The%20curves%20below%20come%20from%20the%20independent%20fixed-penalty%20solves.%20The%20marker%0A%20%20%20%20shows%20the%20continuous%20penalty%20selected%20by%20BLVPY.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20Path%2C%0A%20%20%20%20grid_training_mse%2C%0A%20%20%20%20grid_validation_mse%2C%0A%20%20%20%20lbd%2C%0A%20%20%20%20lbd_grid%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20result%2C%0A)%3A%0A%20%20%20%20figure_dir%20%3D%20Path(__file__).resolve().parent%20%2F%20%22figures%22%0A%20%20%20%20figure_dir.mkdir(parents%3DTrue%2C%20exist_ok%3DTrue)%0A%0A%20%20%20%20fig%2C%20axis%20%3D%20plt.subplots(figsize%3D(6.5%2C%204.5))%0A%20%20%20%20axis.semilogx(lbd_grid%2C%20grid_training_mse%2C%20label%3D%22Training%22%2C%20color%3D%22C0%22)%0A%20%20%20%20axis.semilogx(lbd_grid%2C%20grid_validation_mse%2C%20label%3D%22Validation%22%2C%20color%3D%22C3%22)%0A%20%20%20%20axis.scatter(%0A%20%20%20%20%20%20%20%20%5Bfloat(lbd.value)%5D%2C%0A%20%20%20%20%20%20%20%20%5Bfloat(result.objective)%5D%2C%0A%20%20%20%20%20%20%20%20color%3D%22black%22%2C%0A%20%20%20%20%20%20%20%20marker%3D%22X%22%2C%0A%20%20%20%20%20%20%20%20s%3D80%2C%0A%20%20%20%20%20%20%20%20zorder%3D5%2C%0A%20%20%20%20)%0A%20%20%20%20axis.set(%0A%20%20%20%20%20%20%20%20xlabel%3Dr%22%24%5Clambda%24%22%2C%0A%20%20%20%20%20%20%20%20ylabel%3Dr%22%24(1%2Fm)%7B%5C%7CXw-y%5C%7C%7D_2%5E2%24%22%2C%0A%20%20%20%20%20%20%20%20xlim%3D(1e-4%2C%2010)%2C%0A%20%20%20%20)%0A%20%20%20%20axis.grid(alpha%3D0.25)%0A%20%20%20%20axis.legend(frameon%3DFalse%2C%20fontsize%3D13)%0A%20%20%20%20fig.tight_layout()%0A%20%20%20%20figure_path%20%3D%20figure_dir%20%2F%20%22ridge_hyperparameter.pdf%22%0A%20%20%20%20fig.savefig(figure_path%2C%20bbox_inches%3D%22tight%22)%0A%20%20%20%20plt.show()%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
41287fbcf8621ee1818d22f66de38152