%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.11%22%0A%23%20dependencies%20%3D%20%5B%0A%23%20%20%20%20%20%22marimo%3D%3D0.23.13%22%2C%0A%23%20%20%20%20%20%22numpy%3E%3D2.0.0%22%2C%0A%23%20%20%20%20%20%22matplotlib%22%2C%0A%23%20%20%20%20%20%22cvx-linalg%3E%3D0.9.6%22%2C%0A%23%20%20%20%20%20%22nncg%22%2C%0A%23%20%5D%0A%23%0A%23%20%5Btool.uv.sources%5D%0A%23%20nncg%20%3D%20%7B%20path%20%3D%20%22..%2F..%2F..%22%2C%20editable%20%3D%20true%20%7D%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.23.13%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20The%20solver%2C%20with%20and%20without%20equality%20constraints%0A%0A%20%20%20%20%60nncg%60%20exposes%20two%20entry%20points%20that%20share%20one%20active-set%20driver%3A%0A%0A%20%20%20%20-%20%60solve_nnqp(a%2C%20b)%60%20%E2%80%94%20the%20**bound-only**%20program%0A%20%20%20%20%20%20%24%5Cmin_%7Bx%20%5Cge%200%7D%20%5Ctfrac12%20x%5E%5Ctop%20A%20x%20-%20b%5E%5Ctop%20x%24.%0A%20%20%20%20-%20%60solve_nnqp_eq(a%2C%20b%2C%20b_eq%2C%20c_eq)%60%20%E2%80%94%20the%20same%20objective%20**plus**%20a%20linear%0A%20%20%20%20%20%20equality%20system%20%24Bx%20%3D%20c%24.%0A%0A%20%20%20%20The%20active-set%20notebook%20covered%20the%20bound-only%20case%20in%20detail.%20Here%20we%0A%20%20%20%20focus%20on%20what%20the%20equality%20constraints%20add%3A%20a%20**multiplier**%20%24%5Clambda%24%2C%20a%0A%20%20%20%20**saddle%20system**%20on%20each%20free%20set%2C%20and%20a%20%24p%20%5Ctimes%20p%24%20**Schur%20complement**%0A%20%20%20%20that%20eliminates%20%24%5Clambda%24%20using%20the%20same%20matrix-free%20CG%20inner%20solver.%20The%0A%20%20%20%20important%20structural%20fact%20is%20that%20*the%20outer%20loop%20does%20not%20change%20at%20all*%20%E2%80%94%0A%20%20%20%20only%20the%20per-free-set%20subproblem%20does.%0A%0A%20%20%20%20(The%20code%20uses%20the%20library's%20own%20argument%20names%3A%20%60a%60%20for%20the%20SPD%20matrix%0A%20%20%20%20%24A%24%2C%20%60b_eq%60%20for%20the%20equality%20matrix%20%24B%24%2C%20%60c_eq%60%20for%20its%20right-hand%20side.)%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%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20cvx.linalg%20import%20DenseOperator%0A%0A%20%20%20%20from%20nncg%20import%20kkt_violation%2C%20solve_nnqp%2C%20solve_nnqp_eq%0A%0A%20%20%20%20return%20DenseOperator%2C%20kkt_violation%2C%20np%2C%20plt%2C%20solve_nnqp%2C%20solve_nnqp_eq%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%201.%20Without%20equality%20constraints%20(recap)%0A%0A%20%20%20%20For%20the%20bound-only%20problem%20the%20KKT%20system%20is%20the%20LCP%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x%20%5Cge%200%2C%20%5Cqquad%20s%20%3D%20Ax%20-%20b%20%5Cge%200%2C%20%5Cqquad%20x%20%5Codot%20s%20%3D%200%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20On%20a%20fixed%20free%20set%20%24F%24%20the%20subproblem%20is%20a%20single%20reduced%20SPD%20solve%0A%20%20%20%20%24A_%7BFF%7D%20x_F%20%3D%20b_F%24%2C%20done%20matrix-free%20by%20CG.%20That%20is%20the%20whole%20per-step%20cost.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(DenseOperator%2C%20kkt_violation%2C%20np%2C%20solve_nnqp)%3A%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(0)%0A%20%20%20%20_q%2C%20_%20%3D%20np.linalg.qr(_rng.standard_normal((80%2C%2080)))%0A%20%20%20%20a%20%3D%20(_q%20*%20np.geomspace(1.0%2C%201e3%2C%2080))%20%40%20_q.T%0A%20%20%20%20a%20%3D%200.5%20*%20(a%20%2B%20a.T)%0A%20%20%20%20b%20%3D%20_rng.standard_normal(80)%0A%0A%20%20%20%20op%20%3D%20DenseOperator(a)%0A%20%20%20%20res%20%3D%20solve_nnqp(op%2C%20b)%0A%20%20%20%20print(%22bound-only%20%20solve_nnqp%22)%0A%20%20%20%20print(%22%20%20outer%20steps%20%20%20%20%20%20%20%20%3A%22%2C%20res.outer)%0A%20%20%20%20print(%22%20%20inner%20CG%20iterations%3A%22%2C%20res.inner)%0A%20%20%20%20print(%22%20%20converged%20%20%20%20%20%20%20%20%20%20%3A%22%2C%20res.converged)%0A%20%20%20%20print(%22%20%20KKT%20violation%20%20%20%20%20%20%3A%22%2C%20kkt_violation(op%2C%20b%2C%20res.x))%0A%20%20%20%20print(%22%20%20%23%20active%20(x%3D%3D0)%20%20%20%20%3A%22%2C%20int((res.x%20%3C%3D%201e-9).sum())%2C%20%22of%22%2C%20len(b))%0A%20%20%20%20print(%22%20%20sum(x)%20%20%20%20%20%20%20%20%20%20%20%20%20%3A%22%2C%20round(float(res.x.sum())%2C%204)%2C%20%22%20(unconstrained)%22)%0A%20%20%20%20return%20a%2C%20b%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%202.%20Adding%20equality%20constraints%0A%0A%20%20%20%20Now%20constrain%20the%20same%20problem%20with%20%24B%20x%20%3D%20c%24%2C%20where%20%24B%20%5Cin%20%5Cmathbb%7BR%7D%5E%7Bp%0A%20%20%20%20%5Ctimes%20n%7D%24%20has%20full%20row%20rank%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmin_%7Bx%20%5Cge%200%7D%5C%20%5Ctfrac12%20x%5E%5Ctop%20A%20x%20-%20b%5E%5Ctop%20x%0A%20%20%20%20%5Cquad%20%5Ctext%7Bs.t.%7D%5Cquad%20B%20x%20%3D%20c%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Attach%20a%20multiplier%20%24%5Clambda%20%5Cin%20%5Cmathbb%7BR%7D%5Ep%24%20to%20the%20equalities.%20The%0A%20%20%20%20Lagrangian%20is%20%24%5Ctfrac12%20x%5E%5Ctop%20A%20x%20-%20b%5E%5Ctop%20x%20-%20%5Clambda%5E%5Ctop%20(Bx%20-%20c)%24%2C%20and%0A%20%20%20%20the%20**constrained%20reduced%20gradient**%20becomes%0A%0A%20%20%20%20%24%24%0A%20%20%20%20s%20%5C%3B%3D%5C%3B%20A%20x%20-%20b%20-%20B%5E%5Ctop%20%5Clambda%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20KKT%20%2F%20complementarity%20conditions%20read%20exactly%20as%20before%20*but%20in%20this%0A%20%20%20%20shifted%20%24s%24*%2C%20together%20with%20primal%20feasibility%20of%20the%20equalities%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x%20%5Cge%200%2C%5Cquad%20s%20%5Cge%200%2C%5Cquad%20x%20%5Codot%20s%20%3D%200%2C%5Cqquad%20B%20x%20%3D%20c%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20multiplier%20%24%5Clambda%24%20is%20**free**%20(no%20sign%20restriction)%20%E2%80%94%20it%20is%0A%20%20%20%20determined%2C%20not%20searched%20over.%20The%20active-set%20loop%20still%20toggles%20only%20the%0A%20%20%20%20*bound*%20constraints%3B%20the%20equalities%20are%20enforced%20exactly%20at%20every%20step.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%203.%20The%20saddle%20system%20on%20a%20free%20set%0A%0A%20%20%20%20Fix%20a%20free%20set%20%24F%24.%20With%20%24x_%7BF%5Ec%7D%3D0%24%2C%20stationarity%20on%20%24F%24%20plus%20the%0A%20%20%20%20equalities%20give%20the%20**saddle%20(KKT)%20system**%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Bbmatrix%7D%20A_%7BFF%7D%20%26%20B_F%5E%5Ctop%20%5C%5C%20B_F%20%26%200%20%5Cend%7Bbmatrix%7D%0A%20%20%20%20%5Cbegin%7Bbmatrix%7D%20x_F%20%5C%5C%20-%5Clambda%20%5Cend%7Bbmatrix%7D%0A%20%20%20%20%3D%0A%20%20%20%20%5Cbegin%7Bbmatrix%7D%20b_F%20%5C%5C%20c%20%5Cend%7Bbmatrix%7D%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24B_F%24%20is%20%24B%24%20restricted%20to%20the%20free%20columns.%20This%20is%20indefinite%2C%20so%0A%20%20%20%20we%20do%20**not**%20hand%20it%20to%20CG%20directly.%20Instead%20we%20eliminate%20%24%5Clambda%24%20with%20a%0A%20%20%20%20**Schur%20complement**.%20From%20the%20first%20block%2C%0A%20%20%20%20%24x_F%20%3D%20A_%7BFF%7D%5E%7B-1%7D(b_F%20%2B%20B_F%5E%5Ctop%20%5Clambda)%24.%20Substituting%20into%20%24B_F%20x_F%20%3D%20c%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cunderbrace%7BB_F%20A_%7BFF%7D%5E%7B-1%7D%20B_F%5E%5Ctop%7D_%7BS%20%5C%2C%5Cin%5C%2C%20%5Cmathbb%7BR%7D%5E%7Bp%5Ctimes%20p%7D%7D%5C%2C%0A%20%20%20%20%5Clambda%20%5C%3B%3D%5C%3B%20c%20-%20B_F%20A_%7BFF%7D%5E%7B-1%7D%20b_F%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20%24S%24%20is%20%24p%20%5Ctimes%20p%24%20and%20**SPD**%20(because%20%24A_%7BFF%7D%5Csucc%200%24%20and%20%24B_F%24%20has%20full%0A%20%20%20%20row%20rank)%2C%20so%20a%20tiny%20Cholesky%20solve%20fixes%20%24%5Clambda%24%3B%20then%20%24x_F%24%20follows%20in%0A%20%20%20%20closed%20form.%20Crucially%2C%20forming%20%24S%24%20and%20the%20right-hand%20side%20needs%20only%0A%0A%20%20%20%20%24%24%0A%20%20%20%20v_0%20%3D%20A_%7BFF%7D%5E%7B-1%7D%20b_F%2C%20%5Cqquad%20v_1%5E%7B(j)%7D%20%3D%20A_%7BFF%7D%5E%7B-1%7D%20(B_F%5E%5Ctop%20e_j)%2C%5C%20j%3D1..p%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20i.e.%20**%24p%2B1%24%20SPD%20solves%20that%20all%20share%20the%20operator%20%24A_%7BFF%7D%24**%20%E2%80%94%20each%20one%20a%0A%20%20%20%20matrix-free%20CG%20call.%20So%20an%20equality-augmented%20step%20costs%20%24p%2B1%24%20CG%20solves%0A%20%20%20%20instead%20of%20one%2C%20and%20the%20matrix%20is%20still%20never%20formed.%20This%20is%20exactly%0A%20%20%20%20%60sub_solve%60%20inside%20%60solve_nnqp_eq%60%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20v0%2C%20k0%20%3D%20cg(matvec_f%2C%20b%5Bidx%5D%2C%20tol%3Dcg_tol%2C%20x0%3Dx0)%20%20%20%20%20%23%20A_F%5E%7B-1%7D%20b_F%0A%20%20%20%20v1%20%3D%20np.zeros((idx.size%2C%20p))%0A%20%20%20%20for%20j%20in%20range(p)%3A%0A%20%20%20%20%20%20%20%20v1%5B%3A%2C%20j%5D%2C%20kj%20%3D%20cg(matvec_f%2C%20b_f%5Bj%5D%2C%20tol%3Dcg_tol)%20%20%23%20A_F%5E%7B-1%7D%20B_F%5ET%20e_j%0A%20%20%20%20schur%20%3D%20b_f%20%40%20v1%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20S%20%3D%20B_F%20A_F%5E%7B-1%7D%20B_F%5ET%0A%20%20%20%20lam%20%20%20%3D%20cholesky_solve(schur%2C%20c_eq%20-%20b_f%20%40%20v0)%20%20%20%20%20%20%20%23%20S%20lam%20%3D%20c%20-%20B_F%20v0%0A%20%20%20%20xf%20%20%20%20%3D%20v0%20%2B%20v1%20%40%20lam%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20x_F%0A%20%20%20%20%60%60%60%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(DenseOperator%2C%20a%2C%20b%2C%20np%2C%20solve_nnqp_eq)%3A%0A%20%20%20%20%23%20Add%20p%20equality%20constraints%20to%20the%20SAME%20(a%2C%20b)%20as%20the%20bound-only%20solve.%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(7)%0A%20%20%20%20p%20%3D%203%0A%20%20%20%20b_eq%20%3D%20_rng.standard_normal((p%2C%20len(b)))%0A%20%20%20%20%23%20Pick%20c%20that%20is%20definitely%20reachable%20by%20some%20x%20%3E%3D%200%20(use%20a%20random%20x0%20%3E%3D%200).%0A%20%20%20%20_x0%20%3D%20np.abs(_rng.standard_normal(len(b)))%0A%20%20%20%20c%20%3D%20b_eq%20%40%20_x0%0A%0A%20%20%20%20op_eq%20%3D%20DenseOperator(a)%0A%20%20%20%20res_eq%20%3D%20solve_nnqp_eq(op_eq%2C%20b%2C%20b_eq%2C%20c)%0A%0A%20%20%20%20print(%22equality-augmented%20%20solve_nnqp_eq%20%20%20(p%20%3D%22%2C%20p%2C%20%22)%22)%0A%20%20%20%20print(%22%20%20outer%20steps%20%20%20%20%20%20%20%20%3A%22%2C%20res_eq.outer)%0A%20%20%20%20print(%22%20%20inner%20CG%20iterations%3A%22%2C%20res_eq.inner%2C%20%22%20(~%22%2C%20p%20%2B%201%2C%20%22x%20the%20bound-only%20work%20per%20step)%22)%0A%20%20%20%20print(%22%20%20converged%20%20%20%20%20%20%20%20%20%20%3A%22%2C%20res_eq.converged)%0A%20%20%20%20print(%22%20%20multipliers%20%20lam%20%20%20%3A%22%2C%20np.round(res_eq.lam%2C%204))%0A%20%20%20%20print(%22%20%20%7C%7C%20B%20x%20-%20c%20%7C%7C%20%20%20%20%20%20%3A%22%2C%20float(np.linalg.norm(b_eq%20%40%20res_eq.x%20-%20c)))%0A%20%20%20%20print(%22%20%20min(x)%20%20%20%20%20%20%20%20%20%20%20%20%20%3A%22%2C%20float(res_eq.x.min())%2C%20%22%20(feasibility%20x%20%3E%3D%200)%22)%0A%20%20%20%20return%20b_eq%2C%20c%2C%20p%2C%20res_eq%0A%0A%0A%40app.cell%0Adef%20_(b_eq%2C%20c%2C%20mo%2C%20np%2C%20p%2C%20res_eq)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20rf%22%22%22%0A%20%20%20%20%20%20%20%20The%20solve%20returns%20a%20multiplier%20vector%20%60lam%60%20of%20shape%20%60(%7Bp%7D%2C)%60%20and%20the%0A%20%20%20%20%20%20%20%20equality%20residual%20%24%5ClVert%20Bx%20-%20c%5CrVert%24%20is%0A%20%20%20%20%20%20%20%20%24%7Bnp.linalg.norm(b_eq%20%40%20res_eq.x%20-%20c)%3A.2e%7D%24%20%E2%80%94%20the%20constraint%20is%20satisfied%20to%0A%20%20%20%20%20%20%20%20solver%20tolerance.%20The%20dual%20test%20in%20the%20outer%20loop%20now%20uses%20the%20**shifted**%0A%20%20%20%20%20%20%20%20reduced%20gradient%20%24s%20%3D%20Ax%20-%20b%20-%20B%5E%5Ctop%20%5Clambda%24%3B%20%60reduced_gradient%60%20in%0A%20%20%20%20%20%20%20%20%60solve_nnqp_eq%60%20adds%20precisely%20that%20%24-B%5E%5Ctop%5Clambda%24%20correction.%20Everything%0A%20%20%20%20%20%20%20%20else%20%E2%80%94%20the%20primal%2Fdual%20violator%20tests%2C%20the%20batch%20exchange%2C%20the%20Bland%0A%20%20%20%20%20%20%20%20fallback%20%E2%80%94%20is%20the%20*identical*%20driver%20%60_active_set_loop%60%20used%20by%20the%0A%20%20%20%20%20%20%20%20bound-only%20solver.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%204.%20The%20%24p%20%3D%201%24%20case%3A%20a%20budget%20%2F%20simplex%20constraint%0A%0A%20%20%20%20The%20single%20most%20common%20equality%20is%20the%20**normalisation**%20%24%5Cmathbf%7B1%7D%5E%5Ctop%20x%0A%20%20%20%20%3D%20%5Cbeta%24%20(with%20%24B%20%3D%20%5Cmathbf%7B1%7D%5E%5Ctop%24%2C%20%24c%20%3D%20%5Cbeta%24).%20Combined%20with%20%24x%20%5Cge%200%24%0A%20%20%20%20this%20is%20the%20**scaled%20simplex**%20%E2%80%94%20the%20feasible%20set%20of%20a%20long-only%20portfolio%0A%20%20%20%20with%20a%20fixed%20budget%2C%20of%20a%20mixture%20weight%20vector%2C%20of%20any%20convex%20combination.%0A%20%20%20%20Here%20%24p%20%3D%201%24%2C%20so%20the%20Schur%20complement%20is%20a%20**scalar**%20and%20%24B_F%24%20trivially%0A%20%20%20%20has%20full%20row%20rank%20on%20any%20non-empty%20free%20set.%20This%20is%20the%20workhorse%20case%20and%0A%20%20%20%20needs%20no%20special%20handling%20%E2%80%94%20it%20is%20just%20%60solve_nnqp_eq%60%20with%20a%20one-row%20%24B%24.%0A%0A%20%20%20%20Slide%20the%20budget%20%24%5Cbeta%24%20below%20and%20watch%20the%20solution%20renormalise%20while%0A%20%20%20%20staying%20non-negative.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20beta_slider%20%3D%20mo.ui.slider(0.2%2C%205.0%2C%20value%3D1.0%2C%20step%3D0.1%2C%20label%3Dr%22budget%20%24%5Cbeta%20%3D%20%5Cmathbf%7B1%7D%5E%5Ctop%20x%24%22)%0A%20%20%20%20beta_slider%0A%20%20%20%20return%20(beta_slider%2C)%0A%0A%0A%40app.cell%0Adef%20_(DenseOperator%2C%20a%2C%20b%2C%20beta_slider%2C%20np%2C%20solve_nnqp_eq)%3A%0A%20%20%20%20beta%20%3D%20beta_slider.value%0A%20%20%20%20b_eq1%20%3D%20np.ones((1%2C%20len(b)))%0A%20%20%20%20c1%20%3D%20np.array(%5Bbeta%5D)%0A%20%20%20%20res_b%20%3D%20solve_nnqp_eq(DenseOperator(a)%2C%20b%2C%20b_eq1%2C%20c1)%0A%20%20%20%20return%20beta%2C%20res_b%0A%0A%0A%40app.cell%0Adef%20_(beta%2C%20mo%2C%20res_b)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20rf%22%22%22%0A%20%20%20%20%20%20%20%20%7C%20quantity%20%7C%20value%20%7C%0A%20%20%20%20%20%20%20%20%7C---%7C---%7C%0A%20%20%20%20%20%20%20%20%7C%20budget%20%24%5Cbeta%24%20%7C%20%7Bbeta%3A.2f%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20%24%5Csum_i%20x_i%24%20(should%20equal%20%24%5Cbeta%24)%20%7C%20**%7Bres_b.x.sum()%3A.4f%7D**%20%7C%0A%20%20%20%20%20%20%20%20%7C%20%24%5Cmin_i%20x_i%24%20(feasibility)%20%7C%20%7Bres_b.x.min()%3A.2e%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20%23%20active%20bounds%20(%24x_i%20%3D%200%24)%20%7C%20%7Bint((res_b.x%20%3C%3D%201e-9).sum())%7D%20of%20%7Blen(res_b.x)%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20scalar%20multiplier%20%24%5Clambda%24%20%7C%20%7Bfloat(res_b.lam%5B0%5D)%3A.4f%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20outer%20steps%20%7C%20%7Bres_b.outer%7D%20%7C%0A%0A%20%20%20%20%20%20%20%20The%20mass%20always%20sums%20to%20%24%5Cbeta%24%20and%20never%20goes%20negative%20%E2%80%94%20the%20constraint%0A%20%20%20%20%20%20%20%20surface%20is%20the%20simplex%20%24%5C%7B%7Bx%20%5Cge%200%2C%5C%20%5Cmathbf1%5E%5Ctop%20x%20%3D%20%5Cbeta%5C%7D%7D%24%2C%20and%20the%0A%20%20%20%20%20%20%20%20support%20(the%20non-zero%20coordinates)%20shifts%20as%20%24%5Cbeta%24%20changes.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(plt%2C%20res_b)%3A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(9%2C%203))%0A%20%20%20%20_ax.bar(range(len(res_b.x))%2C%20res_b.x%2C%20color%3D%22%2311D48E%22)%0A%20%20%20%20_ax.axhline(0.0%2C%20color%3D%22k%22%2C%20lw%3D0.8)%0A%20%20%20%20_ax.set_xlabel(%22coordinate%20%24i%24%22)%0A%20%20%20%20_ax.set_ylabel(%22%24x_i%24%22)%0A%20%20%20%20_ax.set_title(rf%22Optimal%20weights%20on%20the%20simplex%20(%24%5Csum%20x_i%20%3D%20%7Bres_b.x.sum()%3A.2f%7D%24)%22)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%205.%20Cost%20comparison%3A%20bound-only%20vs.%20%24p%24%20constraints%0A%0A%20%20%20%20Because%20each%20equality-augmented%20outer%20step%20solves%20%24p%2B1%24%20CG%20systems%20that%0A%20%20%20%20share%20%24A_%7BFF%7D%24%2C%20the%20**inner**%20work%20scales%20roughly%20like%20%24(p%2B1)%5Ctimes%24%20the%0A%20%20%20%20bound-only%20work%20per%20outer%20step%2C%20while%20the%20**outer**%20step%20count%20stays%20in%20the%0A%20%20%20%20same%20regime.%20The%20panel%20below%20sweeps%20%24p%24%20on%20one%20fixed%20%24(A%2C%20b)%24%20and%20a%20random%0A%20%20%20%20full-rank%20%24B%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(DenseOperator%2C%20np%2C%20plt%2C%20solve_nnqp%2C%20solve_nnqp_eq)%3A%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(3)%0A%20%20%20%20_n%20%3D%20120%0A%20%20%20%20_q%2C%20_%20%3D%20np.linalg.qr(_rng.standard_normal((_n%2C%20_n)))%0A%20%20%20%20_core%20%3D%20(_q%20*%20np.geomspace(1.0%2C%201e3%2C%20_n))%20%40%20_q.T%0A%20%20%20%20a2%20%3D%200.5%20*%20(_core%20%2B%20_core.T)%0A%20%20%20%20b2%20%3D%20_rng.standard_normal(_n)%0A%20%20%20%20op2%20%3D%20DenseOperator(a2)%0A%0A%20%20%20%20ps%20%3D%20%5B0%2C%201%2C%202%2C%204%2C%208%5D%0A%20%20%20%20outers%2C%20inners%20%3D%20%5B%5D%2C%20%5B%5D%0A%20%20%20%20for%20_p%20in%20ps%3A%0A%20%20%20%20%20%20%20%20if%20_p%20%3D%3D%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_r%20%3D%20solve_nnqp(op2%2C%20b2)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_b_eq%20%3D%20_rng.standard_normal((_p%2C%20_n))%0A%20%20%20%20%20%20%20%20%20%20%20%20_x0%20%3D%20np.abs(_rng.standard_normal(_n))%0A%20%20%20%20%20%20%20%20%20%20%20%20_r%20%3D%20solve_nnqp_eq(op2%2C%20b2%2C%20_b_eq%2C%20_b_eq%20%40%20_x0)%0A%20%20%20%20%20%20%20%20outers.append(_r.outer)%0A%20%20%20%20%20%20%20%20inners.append(_r.inner)%0A%0A%20%20%20%20_fig%2C%20(_a1%2C%20_a2)%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(10%2C%203.5))%0A%20%20%20%20_a1.plot(ps%2C%20outers%2C%20%22o-%22%2C%20color%3D%22%2311D48E%22)%0A%20%20%20%20_a1.set_xlabel(%22number%20of%20equality%20constraints%20%24p%24%22)%0A%20%20%20%20_a1.set_ylabel(%22outer%20active-set%20steps%22)%0A%20%20%20%20_a1.set_title(%22Outer%20steps%20stay%20flat%22)%0A%20%20%20%20_a2.plot(ps%2C%20inners%2C%20%22o-%22%2C%20color%3D%22%230b7d55%22)%0A%20%20%20%20_a2.set_xlabel(%22number%20of%20equality%20constraints%20%24p%24%22)%0A%20%20%20%20_a2.set_ylabel(%22total%20inner%20CG%20iterations%22)%0A%20%20%20%20_a2.set_title(r%22Inner%20work%20grows%20%24%5Csim%20(p%7B%2B%7D1)%5Ctimes%24%22)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Summary%0A%0A%20%20%20%20%7C%20%7C%20%60solve_nnqp%60%20%7C%20%60solve_nnqp_eq%60%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20feasible%20set%20%7C%20%24x%20%5Cge%200%24%20%7C%20%24x%20%5Cge%200%2C%5C%20Bx%20%3D%20c%24%20%7C%0A%20%20%20%20%7C%20reduced%20gradient%20%7C%20%24s%20%3D%20Ax%20-%20b%24%20%7C%20%24s%20%3D%20Ax%20-%20b%20-%20B%5E%5Ctop%5Clambda%24%20%7C%0A%20%20%20%20%7C%20per-free-set%20subproblem%20%7C%20one%20SPD%20solve%20%24A_%7BFF%7Dx_F%20%3D%20b_F%24%20%7C%20saddle%20system%20via%20%24p%5Ctimes%20p%24%20Schur%20complement%20%7C%0A%20%20%20%20%7C%20inner%20CG%20solves%20per%20outer%20step%20%7C%201%20%7C%20%24p%20%2B%201%24%20(shared%20operator%20%24A_%7BFF%7D%24)%20%7C%0A%20%20%20%20%7C%20outer%20loop%20%7C%20%60_active_set_loop%60%20%7C%20**same**%20%60_active_set_loop%60%20%7C%0A%20%20%20%20%7C%20extra%20output%20%7C%20%E2%80%94%20%7C%20multipliers%20%60Result.lam%60%20%7C%0A%0A%20%20%20%20The%20equality%20machinery%20is%20entirely%20contained%20in%20the%20per-free-set%0A%20%20%20%20subproblem.%20The%20termination%20guarantee%20%E2%80%94%20batch%20pivots%20guarded%20by%20a%0A%20%20%20%20least-index%20Bland%20fallback%20%E2%80%94%20is%20inherited%20unchanged.%0A%0A%20%20%20%20Next%3A%20**%5BTest%20problems%5D(03_test_problems.html)**%20documents%20the%0A%20%20%20%20planted-optimum%20generators%20(including%20the%20equality-augmented%0A%20%20%20%20%60make_eq_problem%60)%20that%20these%20demos%20and%20the%20paper's%20numerical%20study%20are%0A%20%20%20%20built%20on.%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%20import%20marimo%20as%20mo%0A%0A%20%20%20%20return%20(mo%2C)%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
a33c2869549fb703c368d54b7b1fd437