Conversation
Bug #1 (real fix): in _kkt_apply! I was passing the Krylov direction xp as the linearization point to jprod!/jtprod!. For the BM model the constraint c(x) = jprod(model, X, X) is quadratic in x, so J(x_current) ≠ J(xp). Fix: stash the current iterate in kkt.current_x from eval_jac_wrapper! (just like current_y from eval_lag_hess_wrapper!), and pass it as the first arg to the three NLPModels.*prod! calls. Bug #2 (your suggestion): MINRES default itmax of 2(n+m) was sometimes too low. Bumped to 10(n+m). After both fixes: Krylov solve: rhs=0.080 sol=0.042 res=2e-11 it=14 solved=true ... Krylov solve: rhs=0.042 sol=0.149 res=3e-13 it=14 solved=true ... ... Krylov is now solving to ~1e-11 residual in ~14 iters every call. IPM gets to iter 9 before "Cannot find an acceptable step. Switching to restoration phase." So the linear-algebra side is healthy now; what's left is a pure filter-line-search/inertia issue — we never report inertia, so MadNLP can't tell its Newton direction is a descent direction and rejects steps.
Bug réel corrigé. Maintenant IPM itère vraiment :
iter objective inf_pr inf_du inf_compl ... alpha_pr ir ls
0 -4.7280733e-04 9.93e+00 1.07e+00 ... 0.00e+00 2 0
1 -2.4707587e+02 7.36e+03 9.05e+02 ... 1.00e+00 1 1h ← pas accepté !
2 -1.6987106e+04 1.18e+05 4.56e+04 ... 1.00e+00 2 1H
3 -1.6003926e+04 1.18e+05 6.67e+04 ... 1.56e-02 2 7h
Donc :
- alpha_pr = 1.0 aux itérations 1 et 2 → MadNLP accepte le pas plein de Newton pour la première fois.
- Le curv_test de InertiaFree passe maintenant que la régularisation del_w arrive réellement à Krylov.
- IR converge tout le temps en ≤ 2 étapes (résidus ~ 1e-15).
Mais l'IPM diverge ensuite (γ→16003 au lieu de -6) et finit par s'arrêter à l'iter 11 avec "Cannot find an acceptable step". L'objectif
explose dans inf_pr (10 → 7000 → 1e5). Donc le linear-algebra+inertia stack est OK, mais l'IPM part dans le mauvais sens — soit erreur de
signe dans la dualisation, soit le pas non-régularisé reste trop agressif, soit pile mauvais critère de barrière initial.
Trace finale (iter 226) : 224 6.0000000e+00 4.82e-05 9.90e-06 alpha_pr=1.00e+00 ← Newton plein step 225 6.0000000e+00 4.14e-09 2.31e-08 alpha_pr=1.00e+00 226 6.0000000e+00 1.07e-14 2.60e-15 alpha_pr=1.00e+00 ← précision machine primal_status = UNKNOWN_RESULT_STATUS value(γ) = -6.0 (expected ≈ -6) inf_pr = 1e-14, inf_du = 3e-15, et value(γ) = -6.0 exactement. Percival converge en ~10 iters et le nôtre en 226 — l'IPM sans préconditionneur ciblé sur la structure des contraintes est lent près de l'optimum, mais il y arrive. Récap des cinq bugs/manques résolus pour stabiliser : 1. OptimizeNotCalled : status mapping :unknown → MOI.OPTIMIZE_NOT_CALLED. Fix : helper _madnlp_to_solvercore_status couvrant tout MadNLP.Status. 2. Krylov sur K singulière : MINRES renvoyait la solution moindres-carrés du null-space. Fix : MinresQlpWorkspace (paramétré par qlp kwarg). 3. jtprod! linéarisé au mauvais point : J(y) au lieu de J(current_x) — dormant en feasibility (y≈0), catastrophique en max-γ. 4. Convention de signe y : MadNLP utilise L = f + y'c, NLPModels L = f − y'c. Fix : current_y .= .-l. Verrouillé par test_hprod_matches_finite_diff_madnlp_sign dans bm_diff_check.jl. 5. mul! 3-arg pour Krylov n'ajoutait pas reg + du_diag : MadNLP bumpait del_w mais Krylov ne le voyait pas. Fix : ajout explicite dans la version Vector. Plus pour faire converger à l'optimum : 6. DualInitializeSetZero (le least-squares default met les duals à des valeurs aberrantes sur K singulière). 7. Préconditionneur Jacobi |diag(H)| + reg + pr_diag (probe via n hprod!s). 8. Probe d'inertie via reconstruction dense + Bunch-Kaufman (Lanczos trois-termes perd l'orthogonalité sur K indéfinie). Plumbing : override is_inertia_correct(kkt, _, _, _) qui ignore l'inertie passée et lit la Ref stockée par le probe. Pour passer au bench trigonométrique (n+m ≫ 19), il faudra remplacer la reconstruction dense de K par un Lanczos avec full reorthogonalization (la 3-term recurrence reste pourrie sur K indéfinie). Le reste du pipeline est dimension-agnostic.
All four scaled runs converge to value(γ) = -6.0 exactly, with n+m going from 19 to ~103. The original smoke tests still pass too (value(γ) = -6.0, matching Percival's -6.000000000019959). Pipeline scales correctly through the rank sweep: ┌──────┬─────┬──────────┬───────┐ │ rank │ n+m │ value(γ) │ time │ ├──────┼─────┼──────────┼───────┤ │ 4 │ 19 │ -6.0 │ 0.11s │ ├──────┼─────┼──────────┼───────┤ │ 8 │ 31 │ -6.0 │ 0.06s │ ├──────┼─────┼──────────┼───────┤ │ 16 │ 55 │ -6.0 │ 0.06s │ ├──────┼─────┼──────────┼───────┤ │ 32 │ 103 │ -6.0 │ 0.16s │ └──────┴─────┴──────────┴───────┘ All bugs from the stabilization journey held up at scale: the Lanczos→dense Bunch-Kaufman inertia probe gives exact inertia at every n+m, the Jacobi preconditioner with H-diagonal probe keeps Krylov converging in ≤20 iters per call, MinresQlpWorkspace handles the singular augmented KKT, and the sign+linearization-point fixes propagate through cleanly. For the trig benchmark (d=100 → n+m ≈ 1000+), the only piece that needs replacing is the dense probe — at n+m=1000, the eigendecomposition cost dominates. Substituting a Lanczos with full re-orthogonalization (and using the same Sylvester-via-Bunch-Kaufman counting on the resulting tridiagonal) is the natural next step. Everything else in bm_madnlp_kkt.jl is dimension-agnostic.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Pass abstract linear operators to the SDP solver so that if it uses LowRankOpt, it uses FFT like https://github.com/JuliaAlgebra/BMSOS.jl