%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%20Active-set%20methods%20for%20non-negative%20variables%0A%0A%20%20%20%20This%20notebook%20explains%2C%20in%20detail%2C%20the%20active-set%20machinery%20that%20%60nncg%60%0A%20%20%20%20uses%20to%20solve%20the%20strictly%20convex%20**non-negative%20quadratic%20program**%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmin_%7Bx%20%5Cge%200%7D%5C%20%5Ctfrac12%5C%2C%20x%5E%5Ctop%20A%20x%20-%20b%5E%5Ctop%20x%2C%0A%20%20%20%20%5Cqquad%20A%20%3D%20A%5E%5Ctop%20%5Csucc%200%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Everything%20here%20concerns%20the%20bound-only%20problem%20solved%20by%0A%20%20%20%20%60nncg.solve_nnqp%60.%20Equality%20constraints%20get%20their%20own%20notebook%3B%20the%0A%20%20%20%20planted-optimum%20test%20problems%20get%20a%20third.%0A%0A%20%20%20%20We%20build%20the%20story%20bottom-up%3A%0A%0A%20%20%20%201.%20the%20KKT%20%2F%20complementarity%20conditions%20and%20why%20they%20*are*%20the%20solution%2C%0A%20%20%20%202.%20the%20free%20set%20%2F%20active%20set%20and%20the%20reduced%20system%2C%0A%20%20%20%203.%20the%20primal-dual%20violator%20tests%20that%20drive%20the%20iteration%2C%0A%20%20%20%204.%20block%20principal%20pivoting%20(the%20fast%20batch%20exchange)%2C%0A%20%20%20%205.%20why%20the%20batch%20path%20can%20cycle%20and%20how%20the%20**Bland%20fallback**%20guarantees%0A%20%20%20%20%20%20%20finite%20termination%2C%0A%20%20%20%206.%20a%20live%2C%20interactive%20solve%20with%20the%20free-set%20trajectory%20visualised.%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%0A%0A%20%20%20%20return%20DenseOperator%2C%20kkt_violation%2C%20np%2C%20plt%2C%20solve_nnqp%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.%20The%20optimality%20conditions%0A%0A%20%20%20%20The%20objective%20%24f(x)%20%3D%20%5Ctfrac12%20x%5E%5Ctop%20A%20x%20-%20b%5E%5Ctop%20x%24%20is%20strictly%20convex%0A%20%20%20%20(%24A%20%5Csucc%200%24)%2C%20so%20the%20feasible%20set%20%24%5C%7Bx%20%5Cge%200%5C%7D%24%20meets%20it%20at%20a%20**unique**%0A%20%20%20%20global%20minimiser.%20Introduce%20the%20*reduced%20gradient*%0A%0A%20%20%20%20%24%24%0A%20%20%20%20s%20%5C%3B%3D%5C%3B%20%5Cnabla%20f(x)%20%5C%3B%3D%5C%3B%20A%20x%20-%20b%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20A%20point%20%24x%24%20is%20optimal%20iff%20it%20satisfies%20the%20**Karush%E2%80%93Kuhn%E2%80%93Tucker%20(KKT)**%0A%20%20%20%20system%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x%20%5Cge%200%2C%20%5Cqquad%20s%20%5Cge%200%2C%20%5Cqquad%20x%20%5Codot%20s%20%3D%200%20%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24%5Codot%24%20is%20the%20elementwise%20product.%20The%20last%20equation%20is%0A%20%20%20%20**complementarity**%3A%20for%20every%20coordinate%20%24i%24%2C%20either%20%24x_i%20%3D%200%24%20(the%20bound%0A%20%20%20%20is%20*active*)%20or%20%24s_i%20%3D%200%24%20(the%20variable%20is%20*free*%20and%20interior)%2C%20never%20both%0A%20%20%20%20strictly%20positive.%20This%20is%20exactly%20a%20**linear%20complementarity%20problem**%2C%0A%20%20%20%20%24%5Cmathrm%7BLCP%7D(A%2C%20-b)%24%3A%20find%20%24x%20%5Cge%200%24%2C%20%24s%20%3D%20Ax%20-%20b%20%5Cge%200%24%20with%20%24x%5E%5Ctop%20s%20%3D%200%24.%0A%0A%20%20%20%20%60nncg%60%20ships%20this%20test%20as%20a%20certificate%20%E2%80%94%20the%20maximum%20violation%20of%20the%0A%20%20%20%20three%20conditions.%20A%20value%20of%20zero%20*proves*%20global%20optimality.%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%23%20A%20tiny%2C%20hand-checkable%20problem%3A%20a%20%3D%20diag-ish%20SPD%2C%20b%20chosen%20so%20the%20optimum%0A%20%20%20%20%23%20keeps%20some%20coordinates%20at%20the%20bound%20and%20some%20free.%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((6%2C%206)))%0A%20%20%20%20a_small%20%3D%20(_q%20*%20np.geomspace(1.0%2C%2050.0%2C%206))%20%40%20_q.T%0A%20%20%20%20a_small%20%3D%200.5%20*%20(a_small%20%2B%20a_small.T)%0A%20%20%20%20b_small%20%3D%20np.array(%5B2.0%2C%20-1.0%2C%203.0%2C%20-0.5%2C%201.0%2C%200.4%5D)%0A%0A%20%20%20%20op_small%20%3D%20DenseOperator(a_small)%0A%20%20%20%20res_small%20%3D%20solve_nnqp(op_small%2C%20b_small%2C%20track%3DTrue)%0A%0A%20%20%20%20x%20%3D%20res_small.x%0A%20%20%20%20s%20%3D%20a_small%20%40%20x%20-%20b_small%0A%20%20%20%20print(%22x%20%20%20%20%20%20%20%20%3D%22%2C%20np.round(x%2C%204))%0A%20%20%20%20print(%22s%20%3D%20Ax-b%20%20%3D%22%2C%20np.round(s%2C%204))%0A%20%20%20%20print(%22x%20.*%20s%20%20%20%20%3D%22%2C%20np.round(x%20*%20s%2C%206)%2C%20%22%20(complementarity%3A%20all%20~0)%22)%0A%20%20%20%20print(%22KKT%20violation%20%3D%22%2C%20kkt_violation(op_small%2C%20b_small%2C%20x))%0A%20%20%20%20print(%22converged%20%20%20%20%20%3D%22%2C%20res_small.converged)%0A%20%20%20%20return%20(x%2C)%0A%0A%0A%40app.cell%0Adef%20_(mo%2C%20x)%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%20Read%20the%20numbers%20above%20coordinate%20by%20coordinate%3A%20wherever%20%24x_i%20%3E%200%24%20the%0A%20%20%20%20%20%20%20%20reduced%20gradient%20%24s_i%20%5Capprox%200%24%2C%20and%20wherever%20%24s_i%20%3E%200%24%20the%20variable%20sits%0A%20%20%20%20%20%20%20%20at%20its%20bound%20%24x_i%20%3D%200%24.%20That%20is%20complementarity%20in%20action.%0A%0A%20%20%20%20%20%20%20%20-%20**Free%20set**%20%24F%20%3D%20%5C%7B%7B%5C%2C%20i%20%3A%20x_i%20%3E%200%20%5C%2C%5C%7D%7D%24%20%E2%80%94%20the%20interior%20variables.%0A%20%20%20%20%20%20%20%20%20%20Here%20%60x%20%3E%200%60%20at%20%7B%5Bint(i)%20for%20i%20in%20range(len(x))%20if%20x%5Bi%5D%20%3E%201e-9%5D%7D.%0A%20%20%20%20%20%20%20%20-%20**Active%20set**%20%24%5Cmathcal%7B%7BA%7D%7D%20%3D%20%5C%7B%7B%5C%2C%20i%20%3A%20x_i%20%3D%200%20%5C%2C%5C%7D%7D%24%20%E2%80%94%20the%20bound%0A%20%20%20%20%20%20%20%20%20%20variables%2C%20at%20%7B%5Bint(i)%20for%20i%20in%20range(len(x))%20if%20x%5Bi%5D%20%3C%3D%201e-9%5D%7D.%0A%0A%20%20%20%20%20%20%20%20On%20the%20free%20set%20the%20bound%20is%20*inactive*%20and%20the%20problem%20is%20locally%20an%0A%20%20%20%20%20%20%20%20**unconstrained**%20quadratic%20%E2%80%94%20that%20is%20the%20whole%20idea%20an%20active-set%20method%0A%20%20%20%20%20%20%20%20exploits.%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%202.%20The%20free%20set%20and%20the%20reduced%20system%0A%0A%20%20%20%20Suppose%20we%20*knew*%20the%20optimal%20free%20set%20%24F%24.%20Then%20%24x_i%20%3D%200%24%20for%20%24i%20%5Cnotin%20F%24%2C%0A%20%20%20%20and%20on%20%24F%24%20the%20variables%20are%20strictly%20interior%2C%20so%20the%20KKT%20system%20collapses%0A%20%20%20%20to%20a%20plain%20**unconstrained**%20stationarity%20condition%20on%20that%20block%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20A_%7BFF%7D%5C%2C%20x_F%20%5C%3B%3D%5C%3B%20b_F%20%2C%0A%20%20%20%20%5Cqquad%20x_%7BF%5Ec%7D%20%3D%200%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Because%20%24A%20%5Csucc%200%24%2C%20every%20principal%20submatrix%20%24A_%7BFF%7D%24%20is%20SPD%2C%20so%20this%0A%20%20%20%20reduced%20system%20has%20a%20unique%20solution%20and%20can%20be%20solved%20by%20**conjugate%0A%20%20%20%20gradients**%20%E2%80%94%20matrix-free%2C%20at%20the%20%24O(%5Csqrt%7B%5Ckappa%7D)%24%20Krylov%20rate.%20This%20is%0A%20%20%20%20the%20*inner%20solve*%3B%20%60nncg%60%20never%20forms%20or%20factorises%20%24A_%7BFF%7D%24%2C%20it%20only%20needs%0A%20%20%20%20the%20action%20%24v%20%5Cmapsto%20A_%7BFF%7D%20v%24%20(see%20%60nncg.krylov.cg%60%20and%20the%0A%20%20%20%20%60apply_free%60%20%2F%20%60restricted%60%20operator%20hooks).%0A%0A%20%20%20%20The%20catch%2C%20of%20course%2C%20is%20that%20we%20**do%20not**%20know%20%24F%24%20in%20advance.%20The%0A%20%20%20%20active-set%20loop%20is%20the%20search%20for%20it%3A%20guess%20a%20free%20set%2C%20solve%20the%20reduced%0A%20%20%20%20system%2C%20inspect%20who%20violates%20optimality%2C%20and%20correct%20the%20guess.%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%20two%20violator%20tests%0A%0A%20%20%20%20Fix%20a%20candidate%20free%20set%20%24F%24%2C%20solve%20%24A_%7BFF%7D%20x_F%20%3D%20b_F%24%2C%20set%20%24x_%7BF%5Ec%7D%3D0%24%2C%20and%0A%20%20%20%20form%20the%20reduced%20gradient%20%24s%20%3D%20Ax%20-%20b%24%20over%20*all*%20coordinates.%20Two%20things%0A%20%20%20%20can%20be%20wrong%20with%20this%20guess%3A%0A%0A%20%20%20%20-%20**Primal%20violators**%20%24D%20%3D%20%5C%7B%5C%2C%20i%20%5Cin%20F%20%3A%20x_i%20%3C%200%20%5C%2C%5C%7D%24.%20We%20*assumed*%20%24i%24%0A%20%20%20%20%20%20was%20interior%2C%20but%20the%20reduced%20solve%20pushed%20it%20negative%20%E2%80%94%20infeasible.%20Such%0A%20%20%20%20%20%20an%20%24i%24%20must%20be%20**dropped**%20to%20its%20bound.%0A%0A%20%20%20%20-%20**Dual%20violators**%20%24V%20%3D%20%5C%7B%5C%2C%20i%20%5Cnotin%20F%20%3A%20s_i%20%3C%200%20%5C%2C%5C%7D%24.%20We%20*assumed*%20%24i%24%0A%20%20%20%20%20%20was%20pinned%20at%20the%20bound%2C%20but%20its%20reduced%20gradient%20is%20negative%2C%20i.e.%20the%0A%20%20%20%20%20%20objective%20would%20*decrease*%20if%20we%20let%20%24x_i%24%20grow.%20Such%20an%20%24i%24%20must%20be%0A%20%20%20%20%20%20**added**%20to%20the%20free%20set.%0A%0A%20%20%20%20(Both%20tests%20use%20a%20tolerance%20%60tol%60%3B%20%60solve_nnqp%60%20defaults%20to%20%601e-8%60.)%0A%0A%20%20%20%20If%20**neither**%20set%20is%20non-empty%2C%20the%20current%20%24x%24%20satisfies%20%24x%20%5Cge%200%24%2C%0A%20%20%20%20%24s%20%5Cge%200%24%2C%20and%20complementarity%20by%20construction%20(%24x_%7BF%5Ec%7D%3D0%24%20and%2C%20on%20%24F%24%2C%0A%20%20%20%20%24s_F%20%3D%20A_%7BFF%7Dx_F%20-%20b_F%20%3D%200%24).%20That%20is%20the%20KKT%20system%20%E2%80%94%20we%20stop%2C%20and%20the%0A%20%20%20%20point%20is%20the%20unique%20global%20minimiser.%20Otherwise%20we%20exchange%20indices%20and%0A%20%20%20%20repeat.%20In%20%60nncg.solver._active_set_loop%60%20these%20are%20exactly%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20s%20%20%20%20%3D%20reduced_gradient(x%2C%20lam)%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%23%20s%20%3D%20Ax%20-%20b%0A%20%20%20%20prim%20%3D%20np.flatnonzero(free%20%26%20(x%20%3C%20-tol))%20%20%20%20%20%20%20%23%20D%3A%20free%20but%20negative%0A%20%20%20%20dual%20%3D%20np.flatnonzero((~free)%20%26%20(s%20%3C%20-tol))%20%20%20%20%23%20V%3A%20bound%20but%20s%20%3C%200%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_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%204.%20Block%20principal%20pivoting%20%E2%80%94%20the%20fast%20path%0A%0A%20%20%20%20The%20cheapest%20way%20to%20use%20the%20violator%20sets%20is%20a%20**batch%20exchange**%3A%20drop%0A%20%20%20%20*every*%20primal%20violator%20and%20add%20*every*%20dual%20violator%20in%20a%20single%20step%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20F%20%5C%3B%5Cleftarrow%5C%3B%20%5Cbigl(F%20%5Csetminus%20D%5Cbigr)%20%5Ccup%20V%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20This%20is%20**block%20principal%20pivoting**%20(BPP)%3A%20each%20toggle%20is%20a%20*principal%0A%20%20%20%20pivot*%20of%20the%20LCP%2C%20and%20swapping%20a%20whole%20block%20at%20once%20typically%20reaches%20the%0A%20%20%20%20optimal%20free%20set%20in%20a%20handful%20of%20outer%20steps%20%E2%80%94%20far%20fewer%20than%20the%0A%20%20%20%20one-at-a-time%20active-set%20exchanges%20of%20classical%20QP%20solvers.%20In%20the%20driver%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20free%5Bprim%5D%20%3D%20False%20%20%20%23%20drop%20all%20primal%20violators%20D%0A%20%20%20%20free%5Bdual%5D%20%3D%20True%20%20%20%20%23%20add%20%20all%20dual%20%20violators%20V%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20Empirically%20BPP%20converges%20in%20%24O(1)%24%E2%80%93%24O(%5Clog%20n)%24%20outer%20steps%20on%0A%20%20%20%20well-behaved%20problems%2C%20and%20each%20outer%20step%20is%20a%20single%20CG%20solve.%20That%20is%0A%20%20%20%20the%20regime%20%60nncg%60%20lives%20in%20almost%20all%20of%20the%20time.%0A%0A%20%20%20%20But%20%22typically%22%20is%20not%20%22always%22%20%E2%80%94%20a%20pure%20batch%20exchange%20can%20**cycle**%2C%0A%20%20%20%20revisiting%20a%20free%20set%20it%20has%20already%20seen%20and%20looping%20forever.%20Section%205%0A%20%20%20%20is%20the%20fix.%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%205.%20Cycling%20and%20the%20Bland%20fallback%20%E2%80%94%20finite%20termination%0A%0A%20%20%20%20BPP%20can%20over-shoot%3A%20on%20adversarial%20designs%2C%20adding%20a%20block%20of%20variables%0A%20%20%20%20flips%20the%20signs%20of%20their%20partners%2C%20which%20get%20dropped%20next%20step%2C%20which%20flips%0A%20%20%20%20the%20first%20block%20back%20%E2%80%94%20a%20stable%20loop.%20To%20*guarantee*%20termination%20%60nncg%60%0A%20%20%20%20guards%20the%20batch%20path%20with%20a%20**patience%20counter**%20and%20a%20**least-index%0A%20%20%20%20(Bland-type)%20fallback**%2C%20following%20Murty's%20anti-cycling%20rule%20for%20LCPs.%0A%0A%20%20%20%20The%20bookkeeping%20is%20a%20single%20scalar%20%24%5Cbar%20n%24%20%3D%20the%20best%20(smallest)%20number%20of%0A%20%20%20%20violators%20seen%20so%20far%2C%20plus%20a%20patience%20budget%20%60p_max%60%3A%0A%0A%20%20%20%20-%20If%20the%20current%20violator%20count%20%24n_%7B%5Ctext%7Bviol%7D%7D%24%20**improves**%20on%20%24%5Cbar%20n%24%2C%0A%20%20%20%20%20%20take%20the%20batch%20step%20and%20refill%20patience.%20Real%20progress%20is%20never%0A%20%20%20%20%20%20throttled.%0A%20%20%20%20-%20If%20it%20does%20**not**%20improve%20but%20patience%20remains%2C%20take%20the%20batch%20step%0A%20%20%20%20%20%20anyway%20and%20spend%20one%20unit%20of%20patience.%20(Batch%20steps%20are%20cheap%3B%20give%20them%0A%20%20%20%20%20%20a%20few%20tries.)%0A%20%20%20%20-%20If%20it%20does%20not%20improve%20**and%20patience%20is%20exhausted**%2C%20abandon%20the%20batch%0A%20%20%20%20%20%20step%20and%20take%20a%20single%20**Bland%20pivot**%3A%20toggle%20only%20the%20*least-indexed*%0A%20%20%20%20%20%20violator%2C%0A%0A%20%20%20%20%20%20%24%24%0A%20%20%20%20%20%20i%5E%5Cstar%20%3D%20%5Cmin%5Cbigl(D%20%5Ccup%20V%5Cbigr)%2C%20%5Cqquad%20F%20%5Cleftarrow%20F%20%5Ctriangle%20%5C%7Bi%5E%5Cstar%5C%7D.%0A%20%20%20%20%20%20%24%24%0A%0A%20%20%20%20The%20single%20least-index%20pivot%20is%20what%20classical%20Bland%2FMurty%20anti-cycling%0A%20%20%20%20proves%20cannot%20repeat%20a%20working%20set%2C%20so%20the%20loop%20**cannot%20cycle**%20%E2%80%94%20it%0A%20%20%20%20terminates%20in%20finitely%20many%20steps%20at%20the%20unique%20minimiser%2C%20with%20*no%0A%20%20%20%20non-degeneracy%20assumption*%20(Theorem%205.1%20of%20the%20paper).%20This%20is%20the%20exact%0A%20%20%20%20driver%20logic%3A%0A%0A%20%20%20%20%60%60%60python%0A%20%20%20%20if%20n_viol%20%3C%20n_bar%20or%20patience%20%3E%200%3A%20%20%20%20%20%20%23%20fast%20path%0A%20%20%20%20%20%20%20%20if%20n_viol%20%3C%20n_bar%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20n_bar%20%3D%20n_viol%3B%20patience%20%3D%20p_max%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20patience%20-%3D%201%0A%20%20%20%20%20%20%20%20free%5Bprim%5D%20%3D%20False%3B%20free%5Bdual%5D%20%3D%20True%0A%20%20%20%20else%3A%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%20anti-cycling%20fallback%0A%20%20%20%20%20%20%20%20i_star%20%3D%20int(viol.min())%0A%20%20%20%20%20%20%20%20free%5Bi_star%5D%20%3D%20not%20free%5Bi_star%5D%0A%20%20%20%20%60%60%60%0A%0A%20%20%20%20The%20fallback%20is%20not%20a%20theoretical%20ornament%3A%20%60tests%2Ftest_fallback.py%60%0A%20%20%20%20exercises%20it%20on%20the%20%60make_adversarial%60%20family%2C%20and%20it%20is%20the%20load-bearing%0A%20%20%20%20component%20of%20the%20termination%20guarantee.%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%206.%20A%20live%20solve%2C%20with%20the%20free-set%20trajectory%0A%0A%20%20%20%20Turn%20the%20knobs%20below.%20We%20build%20a%20random%20SPD%20problem%20with%20a%20*planted*%0A%20%20%20%20optimum%20of%20known%20support%20(see%20the%20problems%20notebook)%2C%20solve%20it%20with%0A%20%20%20%20%60solve_nnqp(...%2C%20track%3DTrue)%60%2C%20and%20plot%20which%20coordinates%20are%20free%20at%20each%0A%20%20%20%20outer%20step.%20%60track%3DTrue%60%20records%20%60Result.traj%60%20%E2%80%94%20the%20sequence%20of%20visited%0A%20%20%20%20free%20sets.%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%20n_slider%20%3D%20mo.ui.slider(6%2C%2060%2C%20value%3D24%2C%20label%3D%22dimension%20%24n%24%22)%0A%20%20%20%20kappa_slider%20%3D%20mo.ui.slider(1%2C%206%2C%20value%3D3%2C%20step%3D1%2C%20label%3D%22condition%20number%20%24%5C%5Clog_%7B10%7D%5C%5Ckappa%24%22)%0A%20%20%20%20supp_slider%20%3D%20mo.ui.slider(0.1%2C%200.9%2C%20value%3D0.5%2C%20step%3D0.05%2C%20label%3D%22support%20fraction%22)%0A%20%20%20%20seed_slider%20%3D%20mo.ui.slider(0%2C%2030%2C%20value%3D0%2C%20step%3D1%2C%20label%3D%22seed%22)%0A%20%20%20%20mo.vstack(%5Bn_slider%2C%20kappa_slider%2C%20supp_slider%2C%20seed_slider%5D)%0A%20%20%20%20return%20kappa_slider%2C%20n_slider%2C%20seed_slider%2C%20supp_slider%0A%0A%0A%40app.cell%0Adef%20_(kappa_slider%2C%20n_slider%2C%20np%2C%20seed_slider%2C%20supp_slider)%3A%0A%20%20%20%20%23%20Planted-optimum%20generator%20(same%20construction%20as%20tests%2Fproblems.make_problem).%0A%20%20%20%20def%20make_problem(n%2C%20kappa%2C%20support_frac%2C%20seed)%3A%0A%20%20%20%20%20%20%20%20rng%20%3D%20np.random.default_rng(seed)%0A%20%20%20%20%20%20%20%20eig%20%3D%20np.geomspace(1.0%2C%20kappa%2C%20n)%0A%20%20%20%20%20%20%20%20q%2C%20_%20%3D%20np.linalg.qr(rng.standard_normal((n%2C%20n)))%0A%20%20%20%20%20%20%20%20mat%20%3D%200.5%20*%20((q%20*%20eig)%20%40%20q.T%20%2B%20((q%20*%20eig)%20%40%20q.T).T)%0A%20%20%20%20%20%20%20%20k%20%3D%20max(1%2C%20round(support_frac%20*%20n))%0A%20%20%20%20%20%20%20%20perm%20%3D%20rng.permutation(n)%0A%20%20%20%20%20%20%20%20supp%20%3D%20perm%5B%3Ak%5D%0A%20%20%20%20%20%20%20%20x_star%20%3D%20np.zeros(n)%0A%20%20%20%20%20%20%20%20x_star%5Bsupp%5D%20%3D%20rng.uniform(0.5%2C%201.5%2C%20size%3Dk)%0A%20%20%20%20%20%20%20%20s_star%20%3D%20np.zeros(n)%0A%20%20%20%20%20%20%20%20s_star%5Bperm%5Bk%3A%5D%5D%20%3D%20rng.uniform(0.5%2C%201.5%2C%20size%3Dn%20-%20k)%0A%20%20%20%20%20%20%20%20return%20mat%2C%20mat%20%40%20x_star%20-%20s_star%2C%20x_star%0A%0A%20%20%20%20n%20%3D%20n_slider.value%0A%20%20%20%20kappa%20%3D%2010.0**kappa_slider.value%0A%20%20%20%20a%2C%20b%2C%20x_star%20%3D%20make_problem(n%2C%20kappa%2C%20supp_slider.value%2C%20seed_slider.value)%0A%20%20%20%20return%20a%2C%20b%2C%20kappa%2C%20n%2C%20x_star%0A%0A%0A%40app.cell%0Adef%20_(DenseOperator%2C%20a%2C%20b%2C%20kkt_violation%2C%20solve_nnqp)%3A%0A%20%20%20%20op%20%3D%20DenseOperator(a)%0A%20%20%20%20res%20%3D%20solve_nnqp(op%2C%20b%2C%20track%3DTrue)%0A%20%20%20%20kkt%20%3D%20kkt_violation(op%2C%20b%2C%20res.x)%0A%20%20%20%20return%20kkt%2C%20res%0A%0A%0A%40app.cell%0Adef%20_(kappa%2C%20kkt%2C%20mo%2C%20n%2C%20np%2C%20res%2C%20x_star)%3A%0A%20%20%20%20recovered%20%3D%20np.allclose(res.x%2C%20x_star%2C%20atol%3D1e-6)%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%20dimension%20%24n%24%20%7C%20%7Bn%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20condition%20number%20%24%5Ckappa%24%20%7C%20%7Bkappa%3A.0f%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20outer%20active-set%20steps%20%7C%20**%7Bres.outer%7D**%20%7C%0A%20%20%20%20%20%20%20%20%7C%20total%20inner%20CG%20iterations%20%7C%20%7Bres.inner%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20Bland%20fallback%20pivots%20%7C%20%7Bres.fallback%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20converged%20(KKT%20satisfied)%20%7C%20%7Bres.converged%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20KKT%20violation%20%7C%20%7Bkkt%3A.2e%7D%20%7C%0A%20%20%20%20%20%20%20%20%7C%20recovered%20planted%20optimum%20%7C%20%7Brecovered%7D%20%7C%0A%0A%20%20%20%20%20%20%20%20Note%20how%20few%20**outer**%20steps%20it%20takes%20even%20in%20high%20dimension%20%E2%80%94%20that%20is%20the%0A%20%20%20%20%20%20%20%20block-pivot%20fast%20path.%20On%20these%20random%20SPD%20problems%20the%20fallback%20count%20is%0A%20%20%20%20%20%20%20%20normally%20**0**%3B%20you%20need%20the%20adversarial%20family%20to%20force%20it%20positive.%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_(n%2C%20plt%2C%20res)%3A%0A%20%20%20%20%23%20Free-set%20trajectory%3A%20rows%20%3D%20outer%20steps%2C%20columns%20%3D%20coordinates.%0A%20%20%20%20%23%20A%20filled%20cell%20means%20that%20coordinate%20was%20in%20the%20free%20set%20at%20that%20step.%0A%20%20%20%20_traj%20%3D%20res.traj%0A%20%20%20%20_grid%20%3D%20%5B%5B1%20if%20i%20in%20set(step)%20else%200%20for%20i%20in%20range(n)%5D%20for%20step%20in%20_traj%5D%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(9%2C%200.4%20*%20len(_traj)%20%2B%201.5))%0A%20%20%20%20_ax.imshow(_grid%2C%20aspect%3D%22auto%22%2C%20cmap%3D%22Greens%22%2C%20interpolation%3D%22nearest%22)%0A%20%20%20%20_ax.set_xlabel(%22coordinate%20%24i%24%22)%0A%20%20%20%20_ax.set_ylabel(%22outer%20step%22)%0A%20%20%20%20_ax.set_yticks(range(len(_traj)))%0A%20%20%20%20_ax.set_title(%22Free-set%20trajectory%20(green%20%3D%20free%20%2F%20interior%20at%20that%20step)%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%20Each%20row%20is%20one%20outer%20iteration%3B%20a%20green%20cell%20marks%20a%20coordinate%20that%20was%0A%20%20%20%20*free*%20(interior)%20at%20that%20step.%20Watch%20the%20columns%20settle%3A%20after%20the%20first%0A%20%20%20%20couple%20of%20batch%20exchanges%20the%20free%20set%20stops%20changing%20%E2%80%94%20that%20stationary%20row%0A%20%20%20%20pattern%20is%20the%20optimal%20support%2C%20and%20the%20loop%20exits%20on%20the%20KKT%20test%20the%20next%0A%20%20%20%20time%20it%20finds%20no%20violators.%0A%0A%20%20%20%20%23%23%23%20Where%20to%20go%20next%0A%0A%20%20%20%20-%20**%5BEquality%20constraints%5D(02_equality_constraints.html)**%20%E2%80%94%20add%20a%20linear%0A%20%20%20%20%20%20system%20%24Bx%20%3D%20c%24%20(the%20simplex%20%2F%20budget%20constraint)%2C%20solved%20on%20each%20free%0A%20%20%20%20%20%20set%20through%20a%20%24p%5Ctimes%20p%24%20Schur%20complement.%0A%20%20%20%20-%20**%5BTest%20problems%5D(03_test_problems.html)**%20%E2%80%94%20the%20planted-optimum%0A%20%20%20%20%20%20generators%20used%20above%2C%20including%20the%20adversarial%20family%20that%20*forces*%20the%0A%20%20%20%20%20%20Bland%20fallback.%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
212b7838731bbf2848e6e13327765c8e