%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%23%0A%23%20The%20planted-optimum%20generators%20documented%20here%20live%20in%20%60tests%2Fproblems.py%60%0A%23%20(outside%20the%20installed%20package).%20They%20are%20reproduced%20inline%20below%20so%20this%0A%23%20notebook%20is%20self-contained%20and%20runs%20anywhere%20without%20touching%20sys.path.%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%20planted-optimum%20test%20problems%0A%0A%20%20%20%20This%20notebook%20documents%20%60tests%2Fproblems.py%60%20%E2%80%94%20the%20synthetic%20problem%0A%20%20%20%20generators%20that%20power%20%60nncg%60's%20test%20suite%20*and*%20the%20numerical%20study%20of%20the%0A%20%20%20%20accompanying%20paper.%20There%20are%20four%2C%20each%20engineered%20to%20stress%20a%20different%0A%20%20%20%20part%20of%20the%20solver%3A%0A%0A%20%20%20%20%7C%20generator%20%7C%20plants%20%7C%20stresses%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20%60make_problem%60%20%7C%20%24(x%5E%5Cstar%2C%20s%5E%5Cstar)%24%20%7C%20correctness%20across%20condition%20numbers%20%24%5Ckappa%24%20%7C%0A%20%20%20%20%7C%20%60make_eq_problem%60%20%7C%20%24(x%5E%5Cstar%2C%20%5Clambda%5E%5Cstar%2C%20s%5E%5Cstar)%24%20%7C%20the%20equality-augmented%20%2F%20Schur%20path%20%7C%0A%20%20%20%20%7C%20%60make_adversarial%60%20%7C%20*nothing*%20%7C%20the%20Bland%20fallback%20(BPP%20cycles%20here)%20%7C%0A%20%20%20%20%7C%20%60make_scaled_problem%60%20%7C%20%24(x%5E%5Cstar)%24%20%7C%20Jacobi%20preconditioning%20(PCG%20vs%20CG)%20%7C%0A%0A%20%20%20%20The%20unifying%20idea%20is%20**planting%20a%20known%20optimum**.%20Rather%20than%20solve%20a%0A%20%20%20%20problem%20and%20hope%20the%20answer%20is%20right%2C%20we%20*construct*%20the%20answer%20first%20and%0A%20%20%20%20derive%20a%20problem%20that%20has%20it%20%E2%80%94%20then%20a%20solver%20is%20correct%20iff%20it%20recovers%20the%0A%20%20%20%20plant.%20It%20is%20the%20honest%20way%20to%20test%20a%20bound-constrained%20QP%20solver.%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_(np)%3A%0A%20%20%20%20%23%20Inlined%20copies%20of%20the%20four%20generators%20from%20tests%2Fproblems.py%2C%20so%20this%0A%20%20%20%20%23%20notebook%20is%20self-contained.%20Kept%20in%20step%20with%20the%20test%20suite%20%E2%80%94%20if%20you%0A%20%20%20%20%23%20change%20a%20plant%2C%20change%20it%20in%20both%20places.%0A%20%20%20%20def%20make_problem(n%2C%20kappa%2C%20support_frac%3D0.5%2C%20seed%3D0)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Random%20SPD%20problem%20with%20prescribed%20condition%20number%20and%20planted%20optimum.%22%22%22%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%20a%20%3D%20(q%20*%20eig)%20%40%20q.T%0A%20%20%20%20%20%20%20%20a%20%3D%200.5%20*%20(a%20%2B%20a.T)%0A%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%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%20off%20%3D%20perm%5Bk%3A%5D%0A%20%20%20%20%20%20%20%20s_star%5Boff%5D%20%3D%20rng.uniform(0.5%2C%201.5%2C%20size%3Dn%20-%20k)%0A%20%20%20%20%20%20%20%20b%20%3D%20a%20%40%20x_star%20-%20s_star%0A%20%20%20%20%20%20%20%20return%20a%2C%20b%2C%20x_star%2C%20s_star%0A%0A%20%20%20%20def%20make_eq_problem(n%2C%20kappa%2C%20p%2C%20support_frac%3D0.5%2C%20seed%3D0)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Equality-augmented%20planted%20problem%20for%20min%20f(x)%20s.t.%20Bx%20%3D%20c%2C%20x%20%3E%3D%200.%22%22%22%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%20a%20%3D%20(q%20*%20eig)%20%40%20q.T%0A%20%20%20%20%20%20%20%20a%20%3D%200.5%20*%20(a%20%2B%20a.T)%0A%20%20%20%20%20%20%20%20k%20%3D%20max(p%20%2B%201%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%2C%20off%20%3D%20perm%5B%3Ak%5D%2C%20perm%5Bk%3A%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%5Boff%5D%20%3D%20rng.uniform(0.5%2C%201.5%2C%20size%3Dn%20-%20k)%0A%20%20%20%20%20%20%20%20b_eq%20%3D%20rng.standard_normal((p%2C%20n))%0A%20%20%20%20%20%20%20%20lam_star%20%3D%20rng.standard_normal(p)%0A%20%20%20%20%20%20%20%20b%20%3D%20a%20%40%20x_star%20-%20b_eq.T%20%40%20lam_star%20-%20s_star%0A%20%20%20%20%20%20%20%20c_eq%20%3D%20b_eq%20%40%20x_star%0A%20%20%20%20%20%20%20%20return%20a%2C%20b%2C%20b_eq%2C%20c_eq%2C%20x_star%2C%20lam_star%2C%20s_star%0A%0A%20%20%20%20def%20make_adversarial(n%2C%20seed%3D0%2C%20noise%3D1e-2%2C%20ridge%3D1e-6)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Anti-correlated%20design%20on%20which%20the%20unguarded%20batch%20path%20cycles.%22%22%22%0A%20%20%20%20%20%20%20%20rng%20%3D%20np.random.default_rng(seed)%0A%20%20%20%20%20%20%20%20m0%20%3D%20rng.standard_normal((n%2C%20n%20%2F%2F%202))%0A%20%20%20%20%20%20%20%20m%20%3D%20np.hstack(%5Bm0%2C%20-m0%20%2B%20noise%20*%20rng.standard_normal((n%2C%20n%20%2F%2F%202))%5D)%0A%20%20%20%20%20%20%20%20a%20%3D%20m.T%20%40%20m%20%2B%20ridge%20*%20np.eye(n)%0A%20%20%20%20%20%20%20%20return%20a%2C%20m.T%20%40%20rng.standard_normal(n)%0A%0A%20%20%20%20def%20make_scaled_problem(n%2C%20kappa_core%2C%20spread%2C%20support_frac%3D0.5%2C%20seed%3D0)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Well-conditioned%20core%20under%20a%20bad%20diagonal%20scaling%2C%20with%20planted%20optimum.%22%22%22%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_core%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%20core%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%20d%20%3D%20rng.permutation(np.geomspace(1.0%2C%20spread%2C%20n))%0A%20%20%20%20%20%20%20%20a%20%3D%20core%20*%20np.sqrt(np.outer(d%2C%20d))%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%20x_star%20%3D%20np.zeros(n)%0A%20%20%20%20%20%20%20%20x_star%5Bperm%5B%3Ak%5D%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%20a%2C%20a%20%40%20x_star%20-%20s_star%2C%20x_star%0A%0A%20%20%20%20return%20make_adversarial%2C%20make_eq_problem%2C%20make_problem%2C%20make_scaled_problem%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%20The%20planting%20recipe%20(why%20it%20works)%0A%0A%20%20%20%20The%20KKT%20system%20of%20%24%5Cmin_%7Bx%5Cge0%7D%5Ctfrac12%20x%5E%5Ctop%20A%20x%20-%20b%5E%5Ctop%20x%24%20is%20the%20LCP%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x%20%5Cge%200%2C%5Cquad%20s%20%3D%20Ax%20-%20b%20%5Cge%200%2C%5Cquad%20x%20%5Codot%20s%20%3D%200%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20To%20*plant*%20a%20solution%2C%20choose%20a%20complementary%20pair%20%24(x%5E%5Cstar%2C%20s%5E%5Cstar)%24%0A%20%20%20%20directly%3A%20pick%20a%20support%20%24S%24%2C%20set%20%24x%5E%5Cstar%20%3E%200%24%20on%20%24S%24%20and%20%240%24%20off%20it%2C%20and%0A%20%20%20%20%24s%5E%5Cstar%20%3D%200%24%20on%20%24S%24%20and%20%24%3E0%24%20off%20it%20%E2%80%94%20so%20%24x%5E%5Cstar%20%5Codot%20s%5E%5Cstar%20%3D%200%24%20by%0A%20%20%20%20construction.%20Then%20define%20the%20linear%20term%20by%20**reverse-engineering**%20the%0A%20%20%20%20stationarity%20relation%20%24s%20%3D%20Ax%20-%20b%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20b%20%5C%3B%3D%5C%3B%20A%20x%5E%5Cstar%20-%20s%5E%5Cstar%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Now%20%24(x%5E%5Cstar%2C%20s%5E%5Cstar)%24%20satisfies%20the%20LCP%20exactly%2C%20and%20since%20%24A%20%5Csucc%200%24%0A%20%20%20%20makes%20the%20minimiser%20unique%2C%20%24x%5E%5Cstar%24%20*is*%20the%20answer.%20Every%20generator%0A%20%20%20%20below%20is%20a%20variation%20on%20this%20one%20move.%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%201.%20%60make_problem%60%20%E2%80%94%20the%20workhorse%0A%0A%20%20%20%20%24A%20%3D%20Q%20%5Coperatorname%7Bdiag%7D(%5Cmathrm%7Beig%7D)%5C%2C%20Q%5E%5Ctop%24%20with%20a%20Haar-random%0A%20%20%20%20orthogonal%20%24Q%24%20and%20a%20**geometric%20spectrum**%20on%20%24%5B1%2C%20%5Ckappa%5D%24%2C%20so%20the%0A%20%20%20%20condition%20number%20is%20exactly%20the%20knob%20%60kappa%60.%20A%20support%20of%20size%0A%20%20%20%20%24%5Coperatorname%7Bround%7D(%5Ctexttt%7Bsupport%5C_frac%7D%5Ccdot%20n)%24%20is%20drawn%3B%20%24x%5E%5Cstar%24%20is%0A%20%20%20%20uniform%20in%20%24%5B0.5%2C%201.5%5D%24%20there%20and%20%240%24%20elsewhere%2C%20%24s%5E%5Cstar%24%20the%20mirror%2C%20and%0A%20%20%20%20%24b%20%3D%20A%20x%5E%5Cstar%20-%20s%5E%5Cstar%24.%0A%0A%20%20%20%20Returns%20%60(a%2C%20b%2C%20x_star%2C%20s_star)%60.%20It%20drives%20the%20correctness%20sweeps%20%E2%80%94%20a%0A%20%20%20%20solver%20must%20recover%20%60x_star%60%20at%20every%20%24%5Ckappa%24.%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%20make_problem%2C%20np%2C%20solve_nnqp)%3A%0A%20%20%20%20a1%2C%20b1%2C%20x1%2C%20s1%20%3D%20make_problem(n%3D60%2C%20kappa%3D1e4%2C%20support_frac%3D0.5%2C%20seed%3D0)%0A%20%20%20%20op1%20%3D%20DenseOperator(a1)%0A%20%20%20%20r1%20%3D%20solve_nnqp(op1%2C%20b1)%0A%0A%20%20%20%20print(%22make_problem(n%3D60%2C%20kappa%3D1e4%2C%20support_frac%3D0.5)%22)%0A%20%20%20%20print(%22%20%20planted%20support%20size%20%3A%22%2C%20int((x1%20%3E%200).sum()))%0A%20%20%20%20print(%22%20%20recovered%20support%20%20%20%20%3A%22%2C%20int((r1.x%20%3E%201e-9).sum()))%0A%20%20%20%20print(%22%20%20%7C%7C%20x%20-%20x_star%20%7C%7C%20%20%20%20%20%3A%22%2C%20float(np.linalg.norm(r1.x%20-%20x1)))%0A%20%20%20%20print(%22%20%20complementarity%20plant%3A%22%2C%20float(np.max(np.abs(x1%20*%20s1)))%2C%20%22(should%20be%200)%22)%0A%20%20%20%20print(%22%20%20KKT%20violation%20%20%20%20%20%20%20%20%3A%22%2C%20kkt_violation(op1%2C%20b1%2C%20r1.x))%0A%20%20%20%20print(%22%20%20outer%20%2F%20inner%20%20%20%20%20%20%20%20%3A%22%2C%20r1.outer%2C%20%22%2F%22%2C%20r1.inner)%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%20Sweep%20the%20condition%20number%20and%20watch%20the%20solver%20stay%20exact%20while%20the%20inner%0A%20%20%20%20CG%20iterations%20climb%20at%20the%20Krylov%20rate%20%24O(%5Csqrt%7B%5Ckappa%7D)%24.%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%20make_problem%2C%20np%2C%20plt%2C%20solve_nnqp)%3A%0A%20%20%20%20kappas%20%3D%20np.geomspace(1e1%2C%201e6%2C%2012)%0A%20%20%20%20errs%2C%20inner_iters%20%3D%20%5B%5D%2C%20%5B%5D%0A%20%20%20%20for%20_kap%20in%20kappas%3A%0A%20%20%20%20%20%20%20%20_a%2C%20_b%2C%20_x%2C%20_s%20%3D%20make_problem(n%3D80%2C%20kappa%3Dfloat(_kap)%2C%20seed%3D1)%0A%20%20%20%20%20%20%20%20_op%20%3D%20DenseOperator(_a)%0A%20%20%20%20%20%20%20%20_r%20%3D%20solve_nnqp(_op%2C%20_b)%0A%20%20%20%20%20%20%20%20errs.append(kkt_violation(_op%2C%20_b%2C%20_r.x))%0A%20%20%20%20%20%20%20%20inner_iters.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.6))%0A%20%20%20%20_a1.loglog(kappas%2C%20np.maximum(errs%2C%201e-16)%2C%20%22o-%22%2C%20color%3D%22%2311D48E%22)%0A%20%20%20%20_a1.set_xlabel(r%22condition%20number%20%24%5Ckappa%24%22)%0A%20%20%20%20_a1.set_ylabel(%22KKT%20violation%22)%0A%20%20%20%20_a1.set_title(%22Correct%20at%20every%20conditioning%22)%0A%20%20%20%20_a2.semilogx(kappas%2C%20inner_iters%2C%20%22o-%22%2C%20color%3D%22%230b7d55%22)%0A%20%20%20%20_a2.set_xlabel(r%22condition%20number%20%24%5Ckappa%24%22)%0A%20%20%20%20_a2.set_ylabel(%22total%20inner%20CG%20iterations%22)%0A%20%20%20%20_a2.set_title(r%22CG%20cost%20%24%5Csim%20O(%5Csqrt%7B%5Ckappa%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%202.%20%60make_eq_problem%60%20%E2%80%94%20the%20equality-augmented%20plant%0A%0A%20%20%20%20The%20KKT%20system%20now%20carries%20a%20multiplier%3A%20with%20%24B%20%5Cin%20%5Cmathbb%7BR%7D%5E%7Bp%5Ctimes%0A%20%20%20%20n%7D%24%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x%20%5Cge%200%2C%5Cquad%20s%20%3D%20Ax%20-%20b%20-%20B%5E%5Ctop%5Clambda%20%5Cge%200%2C%5Cquad%20x%5Codot%20s%20%3D%200%2C%5Cquad%20Bx%20%3D%20c%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20generator%20plants%20%24(x%5E%5Cstar%2C%20%5Clambda%5E%5Cstar%2C%20s%5E%5Cstar)%24%3A%20a%20support%20of%20size%0A%20%20%20%20**at%20least%20%24p%2B1%24**%20(so%20%24B_F%24%20is%20generically%20full%20row%20rank%20%E2%80%94%20a%20prerequisite%0A%20%20%20%20for%20the%20Schur%20complement%20to%20be%20SPD)%2C%20%24s%5E%5Cstar%20%3D%200%24%20there%20and%20positive%20off%0A%20%20%20%20it%2C%20an%20arbitrary%20%24%5Clambda%5E%5Cstar%24%2C%20and%20then%0A%0A%20%20%20%20%24%24%0A%20%20%20%20b%20%3D%20A%20x%5E%5Cstar%20-%20B%5E%5Ctop%5Clambda%5E%5Cstar%20-%20s%5E%5Cstar%2C%20%5Cqquad%20c%20%3D%20B%20x%5E%5Cstar%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Returns%20%60(a%2C%20b%2C%20b_eq%2C%20c_eq%2C%20x_star%2C%20lam_star%2C%20s_star)%60.%20Note%20the%20solver%0A%20%20%20%20recovers%20%60x_star%60%20exactly%2C%20but%20%60lam%60%20need%20not%20equal%20%60lam_star%60%20unless%20the%0A%20%20%20%20planted%20support%20is%20the%20*unique*%20optimal%20one%20%E2%80%94%20the%20multiplier%20is%20pinned%20by%0A%20%20%20%20the%20active%20support%2C%20which%20is%20what%20actually%20matters.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(DenseOperator%2C%20make_eq_problem%2C%20np%2C%20solve_nnqp_eq)%3A%0A%20%20%20%20a_e%2C%20b_e%2C%20b_eq_e%2C%20c_e%2C%20x_e%2C%20lam_e%2C%20_s_e%20%3D%20make_eq_problem(n%3D60%2C%20kappa%3D1e3%2C%20p%3D3%2C%20seed%3D2)%0A%20%20%20%20op_e%20%3D%20DenseOperator(a_e)%0A%20%20%20%20r_e%20%3D%20solve_nnqp_eq(op_e%2C%20b_e%2C%20b_eq_e%2C%20c_e)%0A%0A%20%20%20%20print(%22make_eq_problem(n%3D60%2C%20kappa%3D1e3%2C%20p%3D3)%22)%0A%20%20%20%20print(%22%20%20%7C%7C%20x%20-%20x_star%20%7C%7C%20%20%20%3A%22%2C%20float(np.linalg.norm(r_e.x%20-%20x_e)))%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_e%20%40%20r_e.x%20-%20c_e)))%0A%20%20%20%20print(%22%20%20planted%20lam_star%20%20%20%3A%22%2C%20np.round(lam_e%2C%204))%0A%20%20%20%20print(%22%20%20recovered%20lam%20%20%20%20%20%20%3A%22%2C%20np.round(r_e.lam%2C%204))%0A%20%20%20%20print(%22%20%20support%20%3E%3D%20p%2B1%20%20%20%20%20%3A%22%2C%20int((x_e%20%3E%200).sum())%2C%20%22%3E%3D%22%2C%203%20%2B%201)%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.%20%60make_adversarial%60%20%E2%80%94%20forcing%20the%20Bland%20fallback%0A%0A%20%20%20%20This%20is%20the%20one%20generator%20that%20plants%20**no**%20optimum%3B%20it%20exists%20to%20break%0A%20%20%20%20the%20fast%20path.%20Columns%20arrive%20in%20near-**anti-parallel**%20pairs%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20M%20%3D%20%5B%5C%2CM_0%20%5Cmid%20-M_0%20%2B%20%5Ctexttt%7Bnoise%7D%5Ccdot%20E%5C%2C%5D%2C%20%5Cqquad%20A%20%3D%20M%5E%5Ctop%20M%20%2B%20%5Ctexttt%7Bridge%7D%5Ccdot%20I%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20so%20pushing%20one%20variable%20to%20its%20bound%20flips%20the%20sign%20of%20its%20partner.%20A%20pure%0A%20%20%20%20block-principal-pivoting%20exchange%20systematically%20over-shoots%20and%20**cycles**%0A%20%20%20%20%E2%80%94%20it%20revisits%20a%20working%20set%20it%20has%20already%20seen%20and%20loops%20forever.%20The%20ridge%0A%20%20%20%20keeps%20%24A%24%20a%20strictly%20positive-definite%20P-matrix.%0A%0A%20%20%20%20The%20guarded%20loop%20in%20%60nncg%60%20terminates%20on%20all%20of%20these%3A%20when%20a%20batch%20step%0A%20%20%20%20stops%20making%20progress%20and%20patience%20runs%20out%2C%20it%20takes%20a%20single%20least-index%0A%20%20%20%20Bland%20pivot%2C%20which%20cannot%20cycle.%20Returns%20%60(a%2C%20b)%60%20only%20%E2%80%94%20certify%20the%20solve%0A%20%20%20%20with%20%60kkt_violation%60.%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%20make_adversarial%2C%20solve_nnqp)%3A%0A%20%20%20%20%23%20Scan%20seeds%3B%20count%20how%20many%20need%20the%20fallback%2C%20and%20confirm%20ALL%20converge.%0A%20%20%20%20n_seeds%20%3D%2040%0A%20%20%20%20fired%2C%20all_ok%2C%20worst_kkt%20%3D%200%2C%20True%2C%200.0%0A%20%20%20%20example_fallback%20%3D%20None%0A%20%20%20%20for%20_seed%20in%20range(n_seeds)%3A%0A%20%20%20%20%20%20%20%20_a%2C%20_b%20%3D%20make_adversarial(n%3D40%2C%20seed%3D_seed)%0A%20%20%20%20%20%20%20%20_op%20%3D%20DenseOperator(_a)%0A%20%20%20%20%20%20%20%20_r%20%3D%20solve_nnqp(_op%2C%20_b)%0A%20%20%20%20%20%20%20%20_k%20%3D%20kkt_violation(_op%2C%20_b%2C%20_r.x)%0A%20%20%20%20%20%20%20%20worst_kkt%20%3D%20max(worst_kkt%2C%20_k)%0A%20%20%20%20%20%20%20%20all_ok%20%3D%20all_ok%20and%20_r.converged%20and%20_k%20%3C%201e-6%0A%20%20%20%20%20%20%20%20if%20_r.fallback%20%3E%200%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20fired%20%2B%3D%201%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20example_fallback%20is%20None%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20example_fallback%20%3D%20(_seed%2C%20_r.fallback%2C%20_r.outer)%0A%0A%20%20%20%20print(f%22scanned%20%7Bn_seeds%7D%20adversarial%20seeds%20(n%3D40)%22)%0A%20%20%20%20print(f%22%20%20seeds%20that%20triggered%20the%20Bland%20fallback%20%3A%20%7Bfired%7D%2F%7Bn_seeds%7D%22)%0A%20%20%20%20print(f%22%20%20all%20converged%20to%20KKT%20%3C%201e-6%20%20%20%20%20%20%20%20%20%20%20%20%3A%20%7Ball_ok%7D%22)%0A%20%20%20%20print(f%22%20%20worst%20KKT%20violation%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%3A%20%7Bworst_kkt%3A.2e%7D%22)%0A%20%20%20%20if%20example_fallback%3A%0A%20%20%20%20%20%20%20%20_sd%2C%20_fb%2C%20_ot%20%3D%20example_fallback%0A%20%20%20%20%20%20%20%20print(f%22%20%20e.g.%20seed%20%7B_sd%7D%3A%20%7B_fb%7D%20fallback%20pivots%20in%20%7B_ot%7D%20outer%20steps%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%20On%20a%20sizeable%20fraction%20of%20seeds%20the%20fallback%20fires%20%E2%80%94%20and%20on%20**every**%20seed%0A%20%20%20%20the%20solve%20still%20terminates%20at%20a%20certified%20optimum.%20That%20is%20the%20termination%0A%20%20%20%20guarantee%20earning%20its%20keep%3B%20%60tests%2Ftest_fallback.py%60%20locks%20this%20behaviour%0A%20%20%20%20in%20as%20a%20regression%2C%20because%20the%20fallback%20path%20is%20the%20load-bearing%20component%0A%20%20%20%20of%20the%20finite-termination%20proof.%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.%20%60make_scaled_problem%60%20%E2%80%94%20a%20case%20for%20preconditioning%0A%0A%20%20%20%20%24A%20%3D%20D%5E%7B1%2F2%7D(Q%5CLambda%20Q%5E%5Ctop)%20D%5E%7B1%2F2%7D%24%3A%20a%20**well-conditioned%20core**%20%24Q%5CLambda%0A%20%20%20%20Q%5E%5Ctop%24%20(condition%20number%20%60kappa_core%60)%20wrapped%20in%20a%20**bad%20diagonal%0A%20%20%20%20scaling**%20%24D%24%20whose%20entries%20span%20%24%5B1%2C%20%5Ctexttt%7Bspread%7D%5D%24.%20The%20full%20matrix%20is%0A%20%20%20%20badly%20conditioned%20%E2%80%94%20roughly%20%60kappa_core%20*%20spread%60%20%E2%80%94%20so%20plain%20CG%20suffers.%20But%0A%20%20%20%20the%20ill-conditioning%20is%20*entirely*%20the%20diagonal%2C%20so%20**Jacobi%0A%20%20%20%20preconditioning**%20(%24M%5E%7B-1%7D%3D%5Coperatorname%7Bdiag%7D(A)%5E%7B-1%7D%24)%20removes%20it%20and%20PCG%0A%20%20%20%20runs%20at%20the%20*core's*%20condition%20number%20regardless%20of%20the%20spread.%0A%0A%20%20%20%20%60solve_nnqp%60%20selects%20the%20inner%20solver%20with%20%60inner%3D%22cg%22%60%20(default)%20or%0A%20%20%20%20%60inner%3D%22jacobi%22%60.%20The%20sweep%20below%20contrasts%20their%20inner%20iteration%20counts%20as%20the%0A%20%20%20%20diagonal%20spread%20grows.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(DenseOperator%2C%20make_scaled_problem%2C%20np%2C%20plt%2C%20solve_nnqp)%3A%0A%20%20%20%20spreads%20%3D%20np.geomspace(1e0%2C%201e6%2C%2010)%0A%20%20%20%20cg_it%2C%20pcg_it%20%3D%20%5B%5D%2C%20%5B%5D%0A%20%20%20%20for%20_sp%20in%20spreads%3A%0A%20%20%20%20%20%20%20%20_a%2C%20_b%2C%20_x%20%3D%20make_scaled_problem(n%3D80%2C%20kappa_core%3D1e2%2C%20spread%3Dfloat(_sp)%2C%20seed%3D4)%0A%20%20%20%20%20%20%20%20_op%20%3D%20DenseOperator(_a)%0A%20%20%20%20%20%20%20%20cg_it.append(solve_nnqp(_op%2C%20_b%2C%20inner%3D%22cg%22).inner)%0A%20%20%20%20%20%20%20%20pcg_it.append(solve_nnqp(_op%2C%20_b%2C%20inner%3D%22jacobi%22).inner)%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(8%2C%204))%0A%20%20%20%20_ax.loglog(spreads%2C%20cg_it%2C%20%22o-%22%2C%20color%3D%22%230b7d55%22%2C%20label%3D'inner%3D%22cg%22')%0A%20%20%20%20_ax.loglog(spreads%2C%20pcg_it%2C%20%22s-%22%2C%20color%3D%22%2311D48E%22%2C%20label%3D'inner%3D%22jacobi%22')%0A%20%20%20%20_ax.set_xlabel(r%22diagonal%20spread%20of%20%24D%24%22)%0A%20%20%20%20_ax.set_ylabel(%22total%20inner%20iterations%22)%0A%20%20%20%20_ax.set_title(%22Jacobi%20PCG%20is%20immune%20to%20diagonal%20scaling%22)%0A%20%20%20%20_ax.legend()%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%20Plain%20CG's%20iteration%20count%20climbs%20with%20the%20spread%3B%20Jacobi-PCG%20stays%20flat%20at%0A%20%20%20%20the%20core's%20cost.%20This%20is%20the%20numerical%20evidence%20behind%20the%20%60inner%3D%22jacobi%22%60%0A%20%20%20%20option%20and%20the%20%60GramOperator%60%2F%60diag%60%20machinery.%0A%0A%20%20%20%20%23%23%23%20Recap%0A%0A%20%20%20%20These%20four%20generators%20are%20deliberately%20kept%20**outside**%20the%20installed%0A%20%20%20%20package%20(in%20%60tests%2Fproblems.py%60%2C%20next%20to%20the%20tests)%20yet%20importable%20%E2%80%94%20so%20they%0A%20%20%20%20serve%20both%20the%20CI%20test%20suite%20and%20exploratory%20notebooks%20like%20this%20one.%20Every%0A%20%20%20%20mathematical%20claim%20in%20the%20paper%20that%20%60nncg%60%20implements%20has%20a%20test%20built%20on%0A%20%20%20%20one%20of%20them%3A%0A%0A%20%20%20%20-%20%60make_problem%60%20%E2%86%92%20correctness%20across%20%24%5Ckappa%24%2C%0A%20%20%20%20-%20%60make_eq_problem%60%20%E2%86%92%20the%20equality-augmented%20Schur%20path%20for%20%24p%5Cin%5C%7B1%2C3%2C8%5C%7D%24%2C%0A%20%20%20%20-%20%60make_adversarial%60%20%E2%86%92%20the%20provably%20necessary%20Bland%20fallback%2C%0A%20%20%20%20-%20%60make_scaled_problem%60%20%E2%86%92%20the%20Jacobi-preconditioning%20win.%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
1cd957df2733b6454d61672ee02b9248