From cb4eaa3e84ee9a10b70b80fec9ee12911509808b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 28 Apr 2026 14:53:48 +0200 Subject: [PATCH 01/15] Add sampling benchmark --- .../Getting started/sampling_bench.jl | 59 +++++++++++++++++++ 1 file changed, 59 insertions(+) create mode 100644 docs/src/tutorials/Getting started/sampling_bench.jl diff --git a/docs/src/tutorials/Getting started/sampling_bench.jl b/docs/src/tutorials/Getting started/sampling_bench.jl new file mode 100644 index 000000000..0ae263e36 --- /dev/null +++ b/docs/src/tutorials/Getting started/sampling_bench.jl @@ -0,0 +1,59 @@ +# # Sampling basis + +#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/Getting started/sampling.ipynb) +#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/Getting started/sampling.ipynb) +# **Contributed by**: Benoît Legat + +using Test #src +using DynamicPolynomials +using SumOfSquares +import MultivariateBases as MB +import SDPLR +import Hypatia +import SCS +import BMSOS + +# In this tutorial, we show how to use a different polynomial basis +# for enforcing the equality between the polynomial and its Sum-of-Squares decomposition. + +@polyvar x +p = x^4 - 4x^3 - 2x^2 + 12x + 3 + +scs = SCS.Optimizer +sdplr = optimizer_with_attributes(SDPLR.Optimizer, "maxrank" => (m, n) -> 4) +hypatia = Hypatia.Optimizer +bmsos = BMSOS.Optimizer +function test(solver, feas::Bool) + model = Model(solver) + set_silent(model) + if feas + γ = -6 + else + @variable(model, γ) + @objective(model, Max, γ) + end + @constraint(model, p - γ in SOSCone(), zero_basis = BoxSampling([-1], [1])) + optimize!(model) + @test primal_status == MOI.FEASIBLE_POINT + if !feasible + @test value(γ) ≈ -6 rtol=1e-4 + end +end +test(scs) +test(sdplr) +test(hypatia) +test(bmsos, true) + +function test_rand(solver, d, B) + model = Model(solver) + set_silent(model) + p = MB.algebra_element(rand(2d+1), MB.SubBasis{B}(monomials(x, 0:2d))) + @constraint(model, p in SOSCone(), zero_basis = BoxSampling([-1], [1])) + optimize!(model) + return solve_time(model) +end + +test_rand(scs, 100, MultivariateBases.Trigonometric) +test_rand(sdplr, 100, MultivariateBases.Trigonometric) +test_rand(hypatia, 100, MultivariateBases.Trigonometric) +test_rand(bmsos, 100, MultivariateBases.Trigonometric) From 09d4889bde2b82b12e7a02a3dbc1ba84486e7d6e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Thu, 28 May 2026 11:43:02 +0200 Subject: [PATCH 02/15] rename --- .../tutorials/Getting started/{sampling_bench.jl => sampling.jl} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename docs/src/tutorials/Getting started/{sampling_bench.jl => sampling.jl} (100%) diff --git a/docs/src/tutorials/Getting started/sampling_bench.jl b/docs/src/tutorials/Getting started/sampling.jl similarity index 100% rename from docs/src/tutorials/Getting started/sampling_bench.jl rename to docs/src/tutorials/Getting started/sampling.jl From dd3a3f2d4ed1d3918bb7c95a107ec4dc37634eea Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 29 May 2026 10:31:37 +0200 Subject: [PATCH 03/15] LowrankOpt with fft --- docs/Project.toml | 4 +- .../src/tutorials/Getting started/sampling.jl | 59 ++- src/Bridges/Variable/lowrank.jl | 60 ++- src/variables.jl | 8 +- test/Bridges/Variable/lowrank_param.jl | 376 ++++++++++++++++++ test/Project.toml | 3 + 6 files changed, 469 insertions(+), 41 deletions(-) create mode 100644 test/Bridges/Variable/lowrank_param.jl diff --git a/docs/Project.toml b/docs/Project.toml index 5f70b20a2..ace27f482 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -5,13 +5,13 @@ Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Cyclotomics = "da8f5974-afbb-4dc8-91d8-516d5257c83b" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" -OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" Dualization = "191a621a-6537-11e9-281d-650236a99e60" DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" GroupsCore = "d5909c97-4eac-4ecc-a3dc-fdd0858a4120" HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327" +Hypatia = "b99e6be6-89ff-11e8-14f8-45c827f4f8f2" ImplicitPlots = "55ecb840-b828-11e9-1645-43f4a9f9ace7" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" KnuthBendix = "c2604015-7b3d-4a30-8a26-9074551ec60a" @@ -23,11 +23,13 @@ MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" PermutationGroups = "8bc5a954-2dfc-11e9-10e6-cd969bffa420" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374" RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01" SCS = "c946c3f1-0d1f-5ce8-9dea-7daa1f7e2d13" +SDPLR = "56161740-ea4e-4253-9d15-43c62ff94d95" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" diff --git a/docs/src/tutorials/Getting started/sampling.jl b/docs/src/tutorials/Getting started/sampling.jl index 0ae263e36..de99d80df 100644 --- a/docs/src/tutorials/Getting started/sampling.jl +++ b/docs/src/tutorials/Getting started/sampling.jl @@ -8,10 +8,6 @@ using Test #src using DynamicPolynomials using SumOfSquares import MultivariateBases as MB -import SDPLR -import Hypatia -import SCS -import BMSOS # In this tutorial, we show how to use a different polynomial basis # for enforcing the equality between the polynomial and its Sum-of-Squares decomposition. @@ -19,10 +15,17 @@ import BMSOS @polyvar x p = x^4 - 4x^3 - 2x^2 + 12x + 3 -scs = SCS.Optimizer -sdplr = optimizer_with_attributes(SDPLR.Optimizer, "maxrank" => (m, n) -> 4) -hypatia = Hypatia.Optimizer -bmsos = BMSOS.Optimizer +# We want to find the minimum of the above polynomial (which is -6). + +model = Model() +set_silent(model) +γ = -6 +@variable(model, γ) +@objective(model, Max, γ) +@constraint(model, p - γ in SOSCone(), zero_basis = BoxSampling([-1], [1])) +set_optimizer() +optimize!(model) + function test(solver, feas::Bool) model = Model(solver) set_silent(model) @@ -35,16 +38,37 @@ function test(solver, feas::Bool) @constraint(model, p - γ in SOSCone(), zero_basis = BoxSampling([-1], [1])) optimize!(model) @test primal_status == MOI.FEASIBLE_POINT - if !feasible + if !feas @test value(γ) ≈ -6 rtol=1e-4 end end -test(scs) + + +import SDPLR + +import BMSOS + +import LowRankOpt as LRO +import Percival +import Dualization + +sdplr = optimizer_with_attributes(SDPLR.Optimizer, "maxrank" => (m, n) -> 4) +bmsos = BMSOS.Optimizer +bmlbfgs = Dualization.dual_optimizer( + optimizer_with_attributes( + LRO.Optimizer, + "solver" => LRO.BurerMonteiro.Solver, + "sub_solver" => Percival.PercivalSolver, + "ranks" => [4], + "square_scalars" => true, + ); + assume_min_if_feasibility = true, +) test(sdplr) -test(hypatia) test(bmsos, true) +test(bmlbfgs, true) -function test_rand(solver, d, B) +function bench_rand(solver, d, B) model = Model(solver) set_silent(model) p = MB.algebra_element(rand(2d+1), MB.SubBasis{B}(monomials(x, 0:2d))) @@ -53,7 +77,14 @@ function test_rand(solver, d, B) return solve_time(model) end -test_rand(scs, 100, MultivariateBases.Trigonometric) -test_rand(sdplr, 100, MultivariateBases.Trigonometric) +import SCS +scs = SCS.Optimizer +bench_rand(scs, 100, MultivariateBases.Trigonometric) + +import Hypatia +hypatia = Hypatia.Optimizer test_rand(hypatia, 100, MultivariateBases.Trigonometric) + +test_rand(sdplr, 100, MultivariateBases.Trigonometric) test_rand(bmsos, 100, MultivariateBases.Trigonometric) +test_rand(bmlbfgs, 100, MultivariateBases.Trigonometric) diff --git a/src/Bridges/Variable/lowrank.jl b/src/Bridges/Variable/lowrank.jl index c55b85dcd..e0fedbc17 100644 --- a/src/Bridges/Variable/lowrank.jl +++ b/src/Bridges/Variable/lowrank.jl @@ -1,17 +1,31 @@ -struct LowRankBridge{T,M} <: MOI.Bridges.Variable.AbstractBridge +struct LowRankBridge{T,M,B,G,W} <: MOI.Bridges.Variable.AbstractBridge affine::Vector{MOI.ScalarAffineFunction{T}} variables::Vector{Vector{MOI.VariableIndex}} constraints::Vector{MOI.ConstraintIndex{MOI.VectorOfVariables}} - set::SOS.WeightedSOSCone{M} + set::SOS.WeightedSOSCone{M,B,G,W} end import LinearAlgebra +# The matrix returned by `MB.transformation_to(gram_basis, target_basis)`. +# `Base.promote_op` resolves this at the type level so we can declare a +# concrete `LRO.Factorization{T, SubArray{...}, Array{T,0}}` in +# `added_constrained_variable_types` — needed because PolyJuMP's +# `bridgeable` rejects `UnionAll` set types. +_transformation_type(::Type{G}, ::Type{B}) where {G,B} = + Base.promote_op(MB.transformation_to, G, B) + +# Row-view of the transformation matrix, i.e. the type of `view(U, j, :)`. +# Used as the `Factorization.factor` type so the parent `U` survives down +# to downstream batched-FFT consumers (e.g. `LowRankOpt.BurerMonteiro`). +_row_view_type(::Type{MT}) where {MT} = + Base.promote_op(view, MT, Int, Colon) + function MOI.Bridges.Variable.bridge_constrained_variable( - ::Type{LowRankBridge{T,M}}, + ::Type{LowRankBridge{T,M,B,G,W}}, model::MOI.ModelLike, - set::SOS.WeightedSOSCone{M}, -) where {T,M} + set::SOS.WeightedSOSCone{M,B,G,W}, +) where {T,M,B,G,W} variables = Vector{Vector{MOI.VariableIndex}}(undef, length(set.gram_bases)) constraints = Vector{MOI.ConstraintIndex{MOI.VectorOfVariables}}( undef, @@ -20,19 +34,26 @@ function MOI.Bridges.Variable.bridge_constrained_variable( for i in eachindex(set.gram_bases) U = MB.transformation_to(set.gram_bases[i], set.basis) weights = SA.coeffs(set.weights[i], set.basis) + # `view(U, j, :)` preserves the parent matrix `U` so that downstream + # `LowRankOpt.BurerMonteiro` can recognise the row-sharing pattern and + # use a single `mul!(buf, U, X.factor)` (one batched FFT per row of `Y` + # when `U::MB.TrigEvalMatrix`) instead of `n` separate inner products. variables[i], constraints[i] = MOI.add_constrained_variables( model, LRO.SetDotProducts{LRO.WITHOUT_SET}( SOS.matrix_cone(M, length(set.gram_bases[i])), [ LRO.TriangleVectorization( - LRO.Factorization(U[j, :], reshape(T[weights[j]], ())), + LRO.Factorization( + view(U, j, :), + reshape(T[weights[j]], ()), + ), ) for j in eachindex(set.basis) ], ), ) end - return LowRankBridge{T,M}( + return LowRankBridge{T,M,B,G,W}( [ MOI.ScalarAffineFunction( [ @@ -59,19 +80,14 @@ function MOI.Bridges.Variable.supports_constrained_variable( end function MOI.Bridges.added_constrained_variable_types( - ::Type{LowRankBridge{T,M}}, -) where {T,M} + ::Type{LowRankBridge{T,M,B,G,W}}, +) where {T,M,B,G,W} + MT = _transformation_type(G, B) + FT = _row_view_type(MT) + TVT = LRO.TriangleVectorization{T,LRO.Factorization{T,FT,Array{T,0}}} return Tuple{Type}[ - ( - LRO.SetDotProducts{ - LRO.WITHOUT_SET, - S[1], - LRO.TriangleVectorization{ - T, - LRO.Factorization{T,Vector{T},Array{T,0}}, - }, - }, - ) for S in SOS.Bridges.Constraint.constrained_variable_types(M) if + (LRO.SetDotProducts{LRO.WITHOUT_SET,S[1],TVT,Vector{TVT}},) + for S in SOS.Bridges.Constraint.constrained_variable_types(M) if S[1] == MOI.PositiveSemidefiniteConeTriangle # FIXME hack ] end @@ -82,9 +98,9 @@ end function MOI.Bridges.Variable.concrete_bridge_type( ::Type{<:LowRankBridge{T}}, - ::Type{<:SOS.WeightedSOSCone{M}}, -) where {T,M} - return LowRankBridge{T,M} + ::Type{<:SOS.WeightedSOSCone{M,B,G,W}}, +) where {T,M,B,G,W} + return LowRankBridge{T,M,B,G,W} end # Attributes, Bridge acting as a model diff --git a/src/variables.jl b/src/variables.jl index 889c6a672..f1f20695d 100644 --- a/src/variables.jl +++ b/src/variables.jl @@ -50,11 +50,11 @@ function PolyJuMP.bridges( S, LRO.TriangleVectorization{ T, - LRO.Factorization{T,Vector{T},Array{T,0}}, + LRO.Factorization{T,F,Array{T,0}}, }, }, }, -) where {S,T} +) where {S,T,F<:AbstractVector{T}} return Tuple{Type,Type}[ (LRO.Bridges.Variable.ToPositiveBridge, T), (LRO.Bridges.Variable.AppendSetBridge, T), @@ -67,11 +67,11 @@ function PolyJuMP.bridges( S, LRO.TriangleVectorization{ T, - LRO.Factorization{T,Vector{T},LRO.One{T}}, + LRO.Factorization{T,F,LRO.One{T}}, }, }, }, -) where {S,T} +) where {S,T,F<:AbstractVector{T}} return Tuple{Type,Type}[(LRO.Bridges.Variable.AppendSetBridge, T)] end function PolyJuMP.bridges( diff --git a/test/Bridges/Variable/lowrank_param.jl b/test/Bridges/Variable/lowrank_param.jl new file mode 100644 index 000000000..7bb79dd51 --- /dev/null +++ b/test/Bridges/Variable/lowrank_param.jl @@ -0,0 +1,376 @@ +# Copyright (c) 2026: Benoît Legat and contributors +# +# Use of this source code is governed by an MIT-style license that can be found +# in the LICENSE.md file or at https://opensource.org/licenses/MIT. +# +# Self-contained tests for the `Bridges.Variable.LowRankBridge` parameterisation +# rewrite (the change that lets `MultivariateBases.TrigEvalMatrix` survive +# intact to `LowRankOpt.BurerMonteiro`). +# +# What this file pins down: +# +# 1. `LowRankBridge` now carries the basis / gram-basis / weight type +# parameters `{T,M,B,G,W}` (was `{T,M}`). +# 2. `_transformation_type(G, B)` / `_row_view_type(MT)` resolve at the +# type level (via `Base.promote_op`) so `added_constrained_variable_types` +# declares the concrete SubArray factor type that +# `bridge_constrained_variable` actually produces. +# 3. The bridge stores `view(U, j, :)` (not `collect(...)`) inside each +# `LRO.Factorization`, so `parent` on the factor returns the underlying +# matrix `U` — the contract that the (future) batched-FFT path in +# BurerMonteiro relies on. +# 4. For `gram_basis = MultivariateBases.Trigonometric`, the same path +# resolves `U::MultivariateBases.TrigEvalMatrix` and the factor type is +# a `SubArray` of that matrix. + +module TestVariableLowRankParam + +using Test +using DynamicPolynomials +using JuMP +using SumOfSquares +import MultivariateBases as MB +import LowRankOpt as LRO +import StarAlgebras as SA +import MathOptInterface as MOI + +# Build a `WeightedSOSCone{M,B,G,W}` of the shape `LowRankBridge` consumes, +# parameterised by the gram-basis type so we can exercise both the dense +# (`Monomial`) and FFT (`Trigonometric`) paths. +function _weighted_sos_cone(::Type{B}; var = :t, gram_degree = 2, n_pts = 5) where {B} + @polyvar t + pts = [Float64[k / (n_pts - 1)] for k in 0:(n_pts-1)] + lag = MB.LagrangeBasis((t,), pts) + gram = MB.SubBasis{B}(monomials(t, 0:gram_degree)) + weight = MB.algebra_element(1.0 * t^0) + return SumOfSquares.WeightedSOSCone{MOI.PositiveSemidefiniteConeTriangle}( + lag, + [gram], + [weight], + ), lag, gram, weight +end + +# 1 ───────────────────────────────────────────────────────────────────────── +# Bridge type carries 5 parameters. +# ───────────────────────────────────────────────────────────────────────── + +function test_bridge_struct_has_five_type_parameters() + BridgeT = SumOfSquares.Bridges.Variable.LowRankBridge + @test BridgeT isa UnionAll + # `LowRankBridge{T}` should still be a UnionAll over the remaining 4 + # parameters — the existing bridge-graph registration uses that syntax. + @test BridgeT{Float64} isa UnionAll + # Concretely instantiate to make sure 5 type parameters resolve. + set, _, gram, weight = _weighted_sos_cone(MB.Monomial) + SetT = typeof(set) + Mp = MOI.PositiveSemidefiniteConeTriangle + Bp = typeof(set.basis) + Gp = eltype(set.gram_bases) + Wp = eltype(set.weights) + concrete = BridgeT{Float64,Mp,Bp,Gp,Wp} + @test concrete <: BridgeT + @test isconcretetype(concrete) || concrete isa DataType +end + +# 2 ───────────────────────────────────────────────────────────────────────── +# `concrete_bridge_type` round-trips: given a `WeightedSOSCone{M,B,G,W}` +# it returns the 5-param `LowRankBridge{T,M,B,G,W}`. +# ───────────────────────────────────────────────────────────────────────── + +function test_concrete_bridge_type_monomial_gram() + set, _, _, _ = _weighted_sos_cone(MB.Monomial) + bridge_t = MOI.Bridges.Variable.concrete_bridge_type( + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + typeof(set), + ) + @test bridge_t isa Type + @test bridge_t <: SumOfSquares.Bridges.Variable.LowRankBridge{Float64} + @test length(bridge_t.parameters) == 5 + @test bridge_t.parameters[1] === Float64 + @test bridge_t.parameters[2] === MOI.PositiveSemidefiniteConeTriangle + @test bridge_t.parameters[3] === typeof(set.basis) + @test bridge_t.parameters[4] === eltype(set.gram_bases) + @test bridge_t.parameters[5] === eltype(set.weights) +end + +function test_concrete_bridge_type_trigonometric_gram() + set, _, _, _ = _weighted_sos_cone(MB.Trigonometric; gram_degree = 4, n_pts = 9) + bridge_t = MOI.Bridges.Variable.concrete_bridge_type( + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + typeof(set), + ) + @test bridge_t <: SumOfSquares.Bridges.Variable.LowRankBridge{Float64} + @test bridge_t.parameters[4] === eltype(set.gram_bases) +end + +# 3 ───────────────────────────────────────────────────────────────────────── +# `_transformation_type(G, B)` is a compile-time map from +# `(gram_basis_type, target_basis_type)` to the matrix type that +# `MB.transformation_to` produces. This is what +# `added_constrained_variable_types` consults. +# ───────────────────────────────────────────────────────────────────────── + +function test_transformation_type_monomial_returns_matrix() + _, lag, gram, _ = _weighted_sos_cone(MB.Monomial) + MT = SumOfSquares.Bridges.Variable._transformation_type(typeof(gram), typeof(lag)) + @test MT !== Any + @test MT <: AbstractMatrix{Float64} + @test MT === typeof(MB.transformation_to(gram, lag)) +end + +function test_transformation_type_trigonometric_returns_trig_eval() + _, lag, gram, _ = _weighted_sos_cone(MB.Trigonometric; gram_degree = 4, n_pts = 9) + MT = SumOfSquares.Bridges.Variable._transformation_type(typeof(gram), typeof(lag)) + @test MT <: MB.TrigEvalMatrix{Float64} + @test MT === typeof(MB.transformation_to(gram, lag)) +end + +function test_row_view_type_matches_view_at_runtime() + _, lag, gram, _ = _weighted_sos_cone(MB.Monomial) + MT = SumOfSquares.Bridges.Variable._transformation_type(typeof(gram), typeof(lag)) + FT = SumOfSquares.Bridges.Variable._row_view_type(MT) + U = MB.transformation_to(gram, lag) + @test FT === typeof(view(U, 1, :)) +end + +function test_row_view_type_for_trig_eval() + _, lag, gram, _ = _weighted_sos_cone(MB.Trigonometric; gram_degree = 4, n_pts = 9) + MT = SumOfSquares.Bridges.Variable._transformation_type(typeof(gram), typeof(lag)) + FT = SumOfSquares.Bridges.Variable._row_view_type(MT) + U = MB.transformation_to(gram, lag) + @test FT === typeof(view(U, 1, :)) + # The factor is a `SubArray` whose parent IS the `TrigEvalMatrix` — the + # invariant downstream batched-FFT consumers will rely on. + v = view(U, 2, :) + @test parent(v) === U + @test parent(v) isa MB.TrigEvalMatrix +end + +# 4 ───────────────────────────────────────────────────────────────────────── +# `added_constrained_variable_types` declares the *same* concrete +# `SetDotProducts` type that `bridge_constrained_variable` actually adds. +# (This is the contract MOI's bridge graph relies on; the older `Vector{T}` +# declaration here would silently mismatch the new `SubArray` factor type +# and PolyJuMP would refuse to bridge.) +# ───────────────────────────────────────────────────────────────────────── + +function _added_constraint_set_type(BridgeT) + types = MOI.Bridges.added_constrained_variable_types(BridgeT) + @assert length(types) == 1 + return types[1][1] +end + +function test_added_constrained_variable_types_monomial() + set, _, _, _ = _weighted_sos_cone(MB.Monomial) + BridgeT = MOI.Bridges.Variable.concrete_bridge_type( + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + typeof(set), + ) + declared = _added_constraint_set_type(BridgeT) + @test declared <: LRO.SetDotProducts{LRO.WITHOUT_SET} + # 4 type parameters (W, S, V, Vs) — the relaxed `SetDotProducts`. + @test length(declared.parameters) == 4 + # Factor type matches `view(::Matrix{Float64}, ::Int, ::Colon)`. + V = declared.parameters[3] + @test V <: LRO.TriangleVectorization{Float64} + Fact = V.parameters[2] + @test Fact <: LRO.Factorization{Float64} + F = Fact.parameters[2] + @test F === typeof(view(Matrix{Float64}(undef, 1, 1), 1, :)) +end + +function test_added_constrained_variable_types_trigonometric() + set, lag, gram, _ = _weighted_sos_cone( + MB.Trigonometric; + gram_degree = 4, + n_pts = 9, + ) + BridgeT = MOI.Bridges.Variable.concrete_bridge_type( + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + typeof(set), + ) + declared = _added_constraint_set_type(BridgeT) + @test declared <: LRO.SetDotProducts{LRO.WITHOUT_SET} + V = declared.parameters[3] + Fact = V.parameters[2] + F = Fact.parameters[2] + # Factor type must be a `SubArray` *of* the `TrigEvalMatrix`, not of a + # plain `Matrix{Float64}` — this is what preserves the FFT path through + # the bridge chain. + U = MB.transformation_to(gram, lag) + @test F === typeof(view(U, 1, :)) + @test F <: SubArray + @test parent(view(U, 1, :)) isa MB.TrigEvalMatrix +end + +# 5 ───────────────────────────────────────────────────────────────────────── +# End-to-end: after `bridge_constrained_variable`, the constraint that lands +# at the inner model has the declared concrete type, and pulling out a +# constraint vector gives a `TriangleVectorization` whose `Factorization`'s +# `.factor` is a `SubArray` with a recoverable parent matrix. +# ───────────────────────────────────────────────────────────────────────── + +function _bridge_into_inner(::Type{B}; kwargs...) where {B} + set, _, _, _ = _weighted_sos_cone(B; kwargs...) + inner = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()) + model = MOI.Bridges.Variable.SingleBridgeOptimizer{ + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + }(inner) + MOI.add_constrained_variables(model, set) + return inner, set +end + +function test_bridge_emits_declared_set_type_monomial() + inner, set = _bridge_into_inner(MB.Monomial) + BridgeT = MOI.Bridges.Variable.concrete_bridge_type( + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + typeof(set), + ) + declared = _added_constraint_set_type(BridgeT) + cts = MOI.get(inner, MOI.ListOfConstraintTypesPresent()) + @test any(((F, S),) -> S === declared, cts) +end + +function test_bridge_emits_declared_set_type_trigonometric() + inner, set = _bridge_into_inner(MB.Trigonometric; gram_degree = 4, n_pts = 9) + BridgeT = MOI.Bridges.Variable.concrete_bridge_type( + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + typeof(set), + ) + declared = _added_constraint_set_type(BridgeT) + cts = MOI.get(inner, MOI.ListOfConstraintTypesPresent()) + @test any(((F, S),) -> S === declared, cts) +end + +function test_bridge_stores_subarray_factor_with_recoverable_parent() + # The bridge must use `view(U, j, :)` rather than materialising the row + # so that `parent(.factor) === U` for each `j`. This is what lets the + # downstream batched-FFT path detect "all rank-1 constraints share a + # common parent matrix". + inner, set = _bridge_into_inner(MB.Trigonometric; gram_degree = 4, n_pts = 9) + declared = _added_constraint_set_type( + MOI.Bridges.Variable.concrete_bridge_type( + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + typeof(set), + ), + ) + cis = MOI.get( + inner, + MOI.ListOfConstraintIndices{MOI.VectorOfVariables,declared}(), + ) + @test length(cis) == 1 + s = MOI.get(inner, MOI.ConstraintSet(), first(cis)) + @test length(s.vectors) == length(set.basis) + factors = [tv.matrix.factor for tv in s.vectors] + @test all(f -> f isa SubArray, factors) + parents = unique(parent.(factors)) + @test length(parents) == 1 + @test only(parents) isa MB.TrigEvalMatrix +end + +function test_bridge_factor_matches_transformation_row_monomial() + inner, set = _bridge_into_inner(MB.Monomial) + declared = _added_constraint_set_type( + MOI.Bridges.Variable.concrete_bridge_type( + SumOfSquares.Bridges.Variable.LowRankBridge{Float64}, + typeof(set), + ), + ) + ci = first( + MOI.get( + inner, + MOI.ListOfConstraintIndices{MOI.VectorOfVariables,declared}(), + ), + ) + s = MOI.get(inner, MOI.ConstraintSet(), ci) + U = MB.transformation_to(only(set.gram_bases), set.basis) + for j in eachindex(set.basis) + # Each factor is the j-th *row* of `U`, accessed lazily through a + # `SubArray` (no materialisation). + @test s.vectors[j].matrix.factor == view(U, j, :) + end +end + +# 6 ───────────────────────────────────────────────────────────────────────── +# `PolyJuMP.bridges(::Type{<:LRO.SetDotProducts{...}})` — the dispatch hook +# in `src/variables.jl` was retyped to accept *any* `F<:AbstractVector{T}` +# factor type (instead of hardcoded `Vector{T}`). Cover both the previous +# `Vector{T}` and the new `SubArray` shape; both must hit the same bridge +# list so the bridge graph routes them identically. +# ───────────────────────────────────────────────────────────────────────── + +function _polyjump_bridges_for(SetType) + return PolyJuMP.bridges(SetType) +end + +function test_polyjump_bridges_dispatches_for_vector_factor() + T = Float64 + V = LRO.TriangleVectorization{T,LRO.Factorization{T,Vector{T},Array{T,0}}} + SetType = LRO.SetDotProducts{ + LRO.WITHOUT_SET, + MOI.PositiveSemidefiniteConeTriangle, + V, + Vector{V}, + } + bridges = _polyjump_bridges_for(SetType) + @test bridges isa Vector + @test !isempty(bridges) + # The classic chain through `ToPositiveBridge` + `AppendSetBridge`. + @test any(b -> b[1] === LRO.Bridges.Variable.ToPositiveBridge, bridges) +end + +function test_polyjump_bridges_dispatches_for_subarray_factor() + T = Float64 + U = Matrix{Float64}(undef, 3, 2) + F = typeof(view(U, 1, :)) + V = LRO.TriangleVectorization{T,LRO.Factorization{T,F,Array{T,0}}} + SetType = LRO.SetDotProducts{ + LRO.WITHOUT_SET, + MOI.PositiveSemidefiniteConeTriangle, + V, + Vector{V}, + } + bridges = _polyjump_bridges_for(SetType) + @test bridges isa Vector + @test !isempty(bridges) + @test any(b -> b[1] === LRO.Bridges.Variable.ToPositiveBridge, bridges) +end + +function test_polyjump_bridges_dispatches_for_trig_eval_subarray_factor() + T = Float64 + # Points must stay in `[-1, 1]` because `Trigonometric.recurrence_eval` + # computes `sqrt(1 - cos²θ)` for sin terms. The `_weighted_sos_cone` + # helper above keeps points in that range; we mirror its layout here. + _, lag, gram, _ = _weighted_sos_cone(MB.Trigonometric; gram_degree = 4, n_pts = 9) + U = MB.transformation_to(gram, lag) + F = typeof(view(U, 1, :)) + @test F <: SubArray + @test eltype(parent(view(U, 1, :))) === T + V = LRO.TriangleVectorization{T,LRO.Factorization{T,F,Array{T,0}}} + SetType = LRO.SetDotProducts{ + LRO.WITHOUT_SET, + MOI.PositiveSemidefiniteConeTriangle, + V, + Vector{V}, + } + bridges = _polyjump_bridges_for(SetType) + @test bridges isa Vector + @test !isempty(bridges) +end + +# ───────────────────────────────────────────────────────────────────────── + +function runtests() + for name in names(@__MODULE__; all = true) + if startswith("$(name)", "test_") + @testset "$(name)" begin + getfield(@__MODULE__, name)() + end + end + end + return +end + +end + +TestVariableLowRankParam.runtests() diff --git a/test/Project.toml b/test/Project.toml index 4d7c1243e..27911dc28 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -4,6 +4,7 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Dualization = "191a621a-6537-11e9-281d-650236a99e60" DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" Hypatia = "b99e6be6-89ff-11e8-14f8-45c827f4f8f2" +JSOSolvers = "10dff2fc-5484-5881-a0e0-c90441020f8a" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" @@ -12,6 +13,7 @@ MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374" SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -19,6 +21,7 @@ StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" SymbolicWedderburn = "858aa9a9-4c7c-4c62-b466-2421203962a2" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +TrigPolys = "bbdedc48-cb31-4a37-9fe3-b015aecc8dd3" TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d" [sources] From e1b1fc3c1e798d1eb098a977506db9804e3c0ff7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 29 May 2026 16:11:40 +0200 Subject: [PATCH 04/15] Add script --- bench/Project.toml | 7 + bench/bm_madnlp_bench.jl | 102 ++++++++++++++ bench/bm_madnlp_kkt.jl | 282 +++++++++++++++++++++++++++++++++++++++ 3 files changed, 391 insertions(+) create mode 100644 bench/bm_madnlp_bench.jl create mode 100644 bench/bm_madnlp_kkt.jl diff --git a/bench/Project.toml b/bench/Project.toml index eb3791fa3..c2cb22b50 100644 --- a/bench/Project.toml +++ b/bench/Project.toml @@ -1,6 +1,13 @@ [deps] +Dualization = "191a621a-6537-11e9-281d-650236a99e60" DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" +LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" +MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" +MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" +SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl new file mode 100644 index 000000000..b0cfb4d91 --- /dev/null +++ b/bench/bm_madnlp_bench.jl @@ -0,0 +1,102 @@ +# Required branches +# +# SumOfSquares → branch `bl/fft` (LowRankOpt + FFT wiring) +# LowRankOpt → branch `bl/sampling` (incl. `hprod!` sign fix) +# MultivariateBases → branch `bl/trigpolys` (TrigEvalMatrix + batched mul!) +# +# Run from this directory: +# +# julia --project=. bm_madnlp_bench.jl + +using SumOfSquares +import DynamicPolynomials +import MultivariateBases as MB +import LowRankOpt as LRO +import MadNLP +import Krylov +import NLPModels +import Dualization +import LinearAlgebra +import MathOptInterface as MOI + +include(joinpath(@__DIR__, "bm_madnlp_kkt.jl")) + +DynamicPolynomials.@polyvar x + +# Custom sub_solver factory that builds a MadNLPSolver with our BMKKTSystem. +function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; kws...) + # Override the random `meta.x0` with a feasibility-friendlier starting + # point. `BMSOSAL` uses `rand(n)` and `BMSOS` uses `randn(n) * 1e-8`; + # `rand(n) / n` keeps norms small as the rank-factor count grows. + n = length(nlp.meta.x0) + nlp.meta.x0 .= rand(eltype(nlp.meta.x0), n) ./ n + # Quick check whether the BM objective is identically zero (which would + # make the dual-Newton step purely feasibility-driven). + f0 = NLPModels.obj(nlp, nlp.meta.x0) + g0 = NLPModels.grad(nlp, nlp.meta.x0) + @info "BM model at x0" obj_x0=f0 grad_norm=LinearAlgebra.norm(g0) nvar=n ncon=length(nlp.meta.y0) + return MadNLP.MadNLPSolver( + nlp; + callback = MadNLP.DenseCallback, + kkt_system = BMKKTSystem, + linear_solver = MadNLP.LapackCPUSolver, # unused — Krylov takes over + nlp_scaling = false, # `jac_dense!` not implemented + print_level = MadNLP.INFO, + max_iter = 50, + ) +end + +inner = optimizer_with_attributes( + LRO.Optimizer, + "solver" => LRO.BurerMonteiro.Solver, + "sub_solver" => MadNLPMinresSolver, + "ranks" => [4], + "square_scalars" => true, +) +bmlbfgs = Dualization.dual_optimizer(inner; assume_min_if_feasibility = true) + +function with_lro_bridges!(model) + backend = JuMP.backend(model) + SumOfSquares.Bridges.add_all_bridges(backend.optimizer, Float64) + MOI.Bridges.remove_bridge( + backend.optimizer, + SumOfSquares.Bridges.Constraint.ImageBridge{Float64}, + ) +end + +println("== BMKKTSystem smoke: max γ s.t. p − γ ∈ SOS ==") +model = Model(bmlbfgs) +set_silent(model) +@variable(model, γ) +@objective(model, Max, γ) +with_lro_bridges!(model) +p = x^4 - 4x^3 - 2x^2 + 12x + 3 +@constraint(model, p - γ in SOSCone(), zero_basis = MB.BoxSampling([-1.0], [1.0])) +optimize!(model) +println("primal_status = ", primal_status(model)) +println("value(γ) = ", value(γ), " (expected ≈ -6)") + +# Should give +# == BMKKTSystem smoke: max γ s.t. p − γ ∈ SOS == +# ┌ Info: BM model at x0 +# │ obj_x0 = 0.00140682890105754 +# │ grad_norm = 0.12033747127060285 +# │ nvar = 14 +# └ ncon = 5 +# [ Info: Custom SolverCore.solve! → regular! only (skipping restoration) +# Number of nonzeros in constraint Jacobian............: 70 +# Number of nonzeros in Lagrangian Hessian.............: 105 +# +# Total number of variables............................: 14 +# variables with only lower bounds: 0 +# variables with lower and upper bounds: 0 +# variables with only upper bounds: 0 +# Total number of equality constraints.................: 5 +# Total number of inequality constraints...............: 0 +# inequality constraints with only lower bounds: 0 +# inequality constraints with lower and upper bounds: 0 +# inequality constraints with only upper bounds: 0 +# +# iter objective inf_pr inf_du inf_compl lg(mu) lg(rg) alpha_pr ir ls +# 0 1.4068289e-03 9.97e+00 1.00e-01 0.00e+00 -1.0 - 0.00e+00 10 0 +# primal_status = NO_SOLUTION diff --git a/bench/bm_madnlp_kkt.jl b/bench/bm_madnlp_kkt.jl new file mode 100644 index 000000000..c7cb5fd2a --- /dev/null +++ b/bench/bm_madnlp_kkt.jl @@ -0,0 +1,282 @@ +# BMKKTSystem — matrix-free MadNLP KKT system for LowRankOpt.BurerMonteiro.Model +# Routes through `BurerMonteiro.Model`'s `NLPModels.hprod!`/`jprod!`/`jtprod!` +# (so the Stage-3 batched-FFT path stays live) and solves each Newton step +# with `Krylov.minres!` (no preconditioner). +# +# Closely modeled after `CompressedSensingIPM/src/fft_kkt.jl`. + +import MadNLP +import Krylov +import NLPModels +import LinearAlgebra +import LowRankOpt as LRO + +mutable struct BMKKTSystem{T,VT,NLP,LS} <: + MadNLP.AbstractReducedKKTSystem{T,VT,Matrix{T},MadNLP.ExactHessian{T,VT}} + nlp::NLP + n::Int + m::Int + # MadNLP-standard diagonals it updates between iterations + reg::VT + pr_diag::VT + du_diag::VT + l_diag::VT + u_diag::VT + l_lower::VT + u_lower::VT + ind_lb::Vector{Int} + ind_ub::Vector{Int} + # Current Lagrangian multipliers — stashed at `eval_lag_hess_wrapper!` time + # and read by `mul!` (the BM Hessian depends on `y` but not on `x`). + current_y::VT + # Krylov state + linear_solver::LS + krylov_iterations::Vector{Int} + krylov_residuals::Vector{Float64} + # Buffers reused inside `mul!` + hv_buf::VT + jv_buf::VT + jtv_buf::VT +end + +function BMKKTSystem(nlp; T = Float64, VT = Vector{T}) + n = NLPModels.get_nvar(nlp) + m = NLPModels.get_ncon(nlp) + workspace = Krylov.MinresWorkspace(n + m, n + m, VT) + return BMKKTSystem{T,VT,typeof(nlp),typeof(workspace)}( + nlp, n, m, + VT(undef, n), # reg + VT(undef, n), # pr_diag + VT(undef, m), # du_diag + VT(undef, 0), # l_diag (no lb) + VT(undef, 0), # u_diag (no ub) + VT(undef, 0), # l_lower + VT(undef, 0), # u_lower + Int[], # ind_lb + Int[], # ind_ub + zeros(T, m), # current_y + workspace, + Int[], + Float64[], + VT(undef, n), # hv_buf + VT(undef, m), # jv_buf + VT(undef, n), # jtv_buf + ) +end + +function MadNLP.create_kkt_system( + ::Type{<:BMKKTSystem}, + cb::MadNLP.AbstractCallback{T,VT}, + linear_solver::Type; + opt_linear_solver = MadNLP.default_options(linear_solver), + hessian_approximation = MadNLP.ExactHessian, + qn_options = MadNLP.QuasiNewtonOptions(), +) where {T,VT} + # Only supported with `square_scalars=true` (no bounds, equality + # constraints only). Verify here so misconfiguration surfaces early. + @assert isempty(cb.ind_ineq) "BMKKTSystem assumes equality-only constraints" + @assert isempty(cb.ind_lb) "BMKKTSystem assumes no lower-bound variables (square_scalars=true)" + @assert isempty(cb.ind_ub) "BMKKTSystem assumes no upper-bound variables (square_scalars=true)" + return BMKKTSystem(cb.nlp; T = T, VT = VT) +end + +MadNLP.num_variables(kkt::BMKKTSystem) = kkt.n +MadNLP.get_hessian(::BMKKTSystem) = nothing +MadNLP.get_jacobian(::BMKKTSystem) = nothing + +# Pretend the Krylov workspace is a "factorization" with no inertia info. +MadNLP.is_inertia(::Krylov.MinresWorkspace) = true +MadNLP.inertia(::Krylov.MinresWorkspace) = (0, 0, 0) +MadNLP.introduce(::Krylov.MinresWorkspace) = "Krylov.MINRES" +MadNLP.improve!(::Krylov.MinresWorkspace) = true +MadNLP.factorize!(::Krylov.MinresWorkspace) = nothing +MadNLP.is_inertia_correct(::BMKKTSystem, _, _, _) = true + +Base.eltype(::BMKKTSystem{T}) where {T} = T +Base.size(kkt::BMKKTSystem) = (kkt.n + kkt.m, kkt.n + kkt.m) +Base.size(kkt::BMKKTSystem, ::Int) = kkt.n + kkt.m + +function MadNLP.initialize!(kkt::BMKKTSystem{T}) where {T} + fill!(kkt.reg, one(T)) + fill!(kkt.pr_diag, one(T)) + fill!(kkt.du_diag, zero(T)) + fill!(kkt.current_y, zero(T)) + return +end + +# Never assemble Jacobian or Hessian — but DO stash the current multipliers. +function MadNLP.eval_jac_wrapper!( + ::MadNLP.MadNLPSolver, ::BMKKTSystem, ::MadNLP.PrimalVector, +) + return +end +function MadNLP.eval_lag_hess_wrapper!( + ::MadNLP.MadNLPSolver, + kkt::BMKKTSystem, + ::MadNLP.PrimalVector, + l::AbstractVector; + is_resto = false, +) + copyto!(kkt.current_y, l) + return +end + +MadNLP.compress_jacobian!(::BMKKTSystem) = nothing +MadNLP.compress_hessian!(::BMKKTSystem) = nothing +MadNLP.build_kkt!(::BMKKTSystem) = nothing +function MadNLP.factorize_wrapper!( + solver::MadNLP.MadNLPSolver{T,VT,IT,KKT}, +) where {T,VT,IT,KKT<:BMKKTSystem} + MadNLP.build_kkt!(solver.kkt) + return true +end + +# Augmented KKT matvec on the `[Δx; Δy]` layout: +# yp = β·yp + α·(H·xp + Aᵀ·xd + pr_diag·xp) +# yd = β·yd + α·(A·xp − du_diag·xd) +function _kkt_apply!( + yp::AbstractVector, yd::AbstractVector, + kkt::BMKKTSystem, + xp::AbstractVector, xd::AbstractVector, + alpha::Number, beta::Number, +) + # Augmented KKT operator (no bounds, no slacks; reg + du_diag added by + # the trailing `_kktmul!`, which mirrors MadNLP's standard sparse path): + # yp = β·yp + α·(H·xp + Aᵀ·xd) + # yd = β·yd + α·(A·xp) + NLPModels.hprod!( + kkt.nlp, xp, kkt.current_y, xp, kkt.hv_buf; + obj_weight = one(eltype(yp)), + ) + NLPModels.jtprod!(kkt.nlp, xp, xd, kkt.jtv_buf) + yp .= beta .* yp .+ alpha .* (kkt.hv_buf .+ kkt.jtv_buf) + + NLPModels.jprod!(kkt.nlp, xp, xp, kkt.jv_buf) + yd .= beta .* yd .+ alpha .* kkt.jv_buf + return +end + +# Variant called by MadNLP on `AbstractKKTVector` (extra `_kktmul!` at the +# end handles the bound-mult blocks; empty in our case but kept for safety). +function MadNLP.mul!( + y::MadNLP.AbstractKKTVector, + kkt::BMKKTSystem, + x::MadNLP.AbstractKKTVector, + alpha::Number, beta::Number, +) + n, m = kkt.n, kkt.m + _x = MadNLP.full(x) + _y = MadNLP.full(y) + _kkt_apply!( + view(_y, 1:n), view(_y, (n+1):(n+m)), + kkt, + view(_x, 1:n), view(_x, (n+1):(n+m)), + alpha, beta, + ) + MadNLP._kktmul!( + y, x, + kkt.reg, kkt.du_diag, kkt.l_lower, kkt.u_lower, kkt.l_diag, kkt.u_diag, + alpha, beta, + ) + return y +end + +# Variant called by `Krylov.kmul!` on plain `Vector{T}` (length `n + m`). +function LinearAlgebra.mul!( + y::AbstractVector, + kkt::BMKKTSystem, + x::AbstractVector, + alpha::Number, beta::Number, +) + n, m = kkt.n, kkt.m + _kkt_apply!( + view(y, 1:n), view(y, (n+1):(n+m)), + kkt, + view(x, 1:n), view(x, (n+1):(n+m)), + alpha, beta, + ) + return y +end + +LinearAlgebra.mul!(y::AbstractVector, kkt::BMKKTSystem, x::AbstractVector) = + LinearAlgebra.mul!(y, kkt, x, true, false) + +# `Aᵀ x`, just routes to the NLP's `jtprod!`. +function MadNLP.jtprod!(y::AbstractVector, kkt::BMKKTSystem, x::AbstractVector) + NLPModels.jtprod!(kkt.nlp, y, x, y) + return y +end + +import SolverCore + +# `BurerMonteiro.Solver`'s `solve!` calls +# `SolverCore.solve!(madnlp_solver, bm_model, stats; kws...)`. MadNLP's own +# entrypoints are `solve!(solver)` or `solve!(nlp, solver, stats)`; the +# (solver, nlp, stats) ordering doesn't exist. Bridge here. +# +# CompressedSensingIPM-style: only run `MadNLP.regular!`, never `robust!`. +# That bypasses restoration entirely — which is fine for us because +# restoration would call `eval_lag_hess_wrapper!`/`eval_jac_wrapper!` (which +# we stub to no-ops) and bail anyway. +# MadNLP `Status` → SolverCore status symbol (consumed by NLPModelsJuMP's +# `TERMINATION_STATUS`; `:unknown` → `MOI.OPTIMIZE_NOT_CALLED`, which makes +# JuMP throw `OptimizeNotCalled` on `value(...)` — avoid). +function _madnlp_to_solvercore_status(s::MadNLP.Status) + s == MadNLP.SOLVE_SUCCEEDED && return :first_order + s == MadNLP.SOLVED_TO_ACCEPTABLE_LEVEL && return :acceptable + s == MadNLP.SEARCH_DIRECTION_BECOMES_TOO_SMALL && return :small_step + s == MadNLP.DIVERGING_ITERATES && return :unbounded + s == MadNLP.INFEASIBLE_PROBLEM_DETECTED && return :infeasible + s == MadNLP.MAXIMUM_ITERATIONS_EXCEEDED && return :max_iter + s == MadNLP.MAXIMUM_WALLTIME_EXCEEDED && return :max_time + s == MadNLP.USER_REQUESTED_STOP && return :user + s == MadNLP.ERROR_IN_STEP_COMPUTATION && return :neg_pred + # In-progress states leaking out (e.g. we returned mid-iteration after the + # line search failed) — treat as a stalled / slow-progress termination so + # JuMP still sees a result and can read the iterate. + s in (MadNLP.INITIAL, MadNLP.REGULAR, MadNLP.RESTORE, MadNLP.ROBUST, + MadNLP.LINESEARCH_SUCCEEDED) && return :stalled + return :exception +end + +function SolverCore.solve!( + solver::MadNLP.MadNLPSolver, + nlp::NLPModels.AbstractNLPModel, + stats::SolverCore.GenericExecutionStats; + kws..., +) + @info "Custom SolverCore.solve! → regular! only (skipping restoration)" + MadNLP.print_init(solver) + MadNLP.initialize!(solver) + try + MadNLP.regular!(solver) + catch err + @info "MadNLP.regular! threw — accepting current state" err + end + res = MadNLP.MadNLPExecutionStats(solver) + stats.solution .= res.solution + stats.objective = res.objective + stats.iter = res.iter + SolverCore.set_status!(stats, _madnlp_to_solvercore_status(res.status)) + stats.dual_feas = res.dual_feas + stats.primal_feas = res.primal_feas + stats.solver_specific[:madnlp] = res + return stats +end + +function MadNLP.solve_kkt!(kkt::BMKKTSystem, w::MadNLP.AbstractKKTVector) + MadNLP.reduce_rhs!(kkt, w) + # The `[primal; dual]` block — what MadNLP's standard + # `solve!(::AbstractReducedKKTSystem, w)` hands to its linear solver. + b = MadNLP.primal_dual(w) + Krylov.krylov_solve!( + kkt.linear_solver, kkt, b; + atol = 1e-10, rtol = 1e-8, verbose = 0, + ) + x = Krylov.solution(kkt.linear_solver) + copyto!(b, x) + push!(kkt.krylov_iterations, Krylov.iteration_count(kkt.linear_solver)) + push!(kkt.krylov_residuals, Krylov.elapsed_time(kkt.linear_solver)) + MadNLP.finish_aug_solve!(kkt, w) + return w +end From 522fdaf5c448c83ca9a3ff6ff94cdf24736cf8b9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 29 May 2026 16:37:13 +0200 Subject: [PATCH 05/15] Fixes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- bench/bm_madnlp_bench.jl | 2 +- bench/bm_madnlp_kkt.jl | 45 +++++++++++++++++++++++++++++++++------- 2 files changed, 39 insertions(+), 8 deletions(-) diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index b0cfb4d91..5c8d5bc5b 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -41,7 +41,7 @@ function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; kws...) kkt_system = BMKKTSystem, linear_solver = MadNLP.LapackCPUSolver, # unused — Krylov takes over nlp_scaling = false, # `jac_dense!` not implemented - print_level = MadNLP.INFO, + print_level = MadNLP.TRACE, max_iter = 50, ) end diff --git a/bench/bm_madnlp_kkt.jl b/bench/bm_madnlp_kkt.jl index c7cb5fd2a..c108e8259 100644 --- a/bench/bm_madnlp_kkt.jl +++ b/bench/bm_madnlp_kkt.jl @@ -29,6 +29,12 @@ mutable struct BMKKTSystem{T,VT,NLP,LS} <: # Current Lagrangian multipliers — stashed at `eval_lag_hess_wrapper!` time # and read by `mul!` (the BM Hessian depends on `y` but not on `x`). current_y::VT + # Current primal iterate — stashed at `eval_jac_wrapper!` time. The BM + # constraint c(x) is *quadratic* in `x` (since x is the rank factor / + # square-scalar pre-image), so its Jacobian depends on the linearization + # point and we must pass the actual IPM iterate to `jprod!`/`jtprod!`, + # not the Krylov direction. + current_x::VT # Krylov state linear_solver::LS krylov_iterations::Vector{Int} @@ -55,6 +61,7 @@ function BMKKTSystem(nlp; T = Float64, VT = Vector{T}) Int[], # ind_lb Int[], # ind_ub zeros(T, m), # current_y + zeros(T, n), # current_x workspace, Int[], Float64[], @@ -101,13 +108,16 @@ function MadNLP.initialize!(kkt::BMKKTSystem{T}) where {T} fill!(kkt.pr_diag, one(T)) fill!(kkt.du_diag, zero(T)) fill!(kkt.current_y, zero(T)) + fill!(kkt.current_x, zero(T)) return end -# Never assemble Jacobian or Hessian — but DO stash the current multipliers. +# Never assemble Jacobian or Hessian — but DO stash the current iterate so +# our matrix-free `mul!` linearizes the (nonlinear) BM constraint at it. function MadNLP.eval_jac_wrapper!( - ::MadNLP.MadNLPSolver, ::BMKKTSystem, ::MadNLP.PrimalVector, + ::MadNLP.MadNLPSolver, kkt::BMKKTSystem, x::MadNLP.PrimalVector, ) + copyto!(kkt.current_x, MadNLP.full(x)) return end function MadNLP.eval_lag_hess_wrapper!( @@ -144,14 +154,17 @@ function _kkt_apply!( # the trailing `_kktmul!`, which mirrors MadNLP's standard sparse path): # yp = β·yp + α·(H·xp + Aᵀ·xd) # yd = β·yd + α·(A·xp) + # Linearize at the *current IPM iterate* (`kkt.current_x`), not at the + # Krylov direction `xp` — the BM Jacobian is x-dependent (quadratic + # constraint). NLPModels.hprod!( - kkt.nlp, xp, kkt.current_y, xp, kkt.hv_buf; + kkt.nlp, kkt.current_x, kkt.current_y, xp, kkt.hv_buf; obj_weight = one(eltype(yp)), ) - NLPModels.jtprod!(kkt.nlp, xp, xd, kkt.jtv_buf) + NLPModels.jtprod!(kkt.nlp, kkt.current_x, xd, kkt.jtv_buf) yp .= beta .* yp .+ alpha .* (kkt.hv_buf .+ kkt.jtv_buf) - NLPModels.jprod!(kkt.nlp, xp, xp, kkt.jv_buf) + NLPModels.jprod!(kkt.nlp, kkt.current_x, xp, kkt.jv_buf) yd .= beta .* yd .+ alpha .* kkt.jv_buf return end @@ -269,13 +282,31 @@ function MadNLP.solve_kkt!(kkt::BMKKTSystem, w::MadNLP.AbstractKKTVector) # The `[primal; dual]` block — what MadNLP's standard # `solve!(::AbstractReducedKKTSystem, w)` hands to its linear solver. b = MadNLP.primal_dual(w) + rhs_norm = LinearAlgebra.norm(b) + rhs_copy = copy(b) Krylov.krylov_solve!( kkt.linear_solver, kkt, b; - atol = 1e-10, rtol = 1e-8, verbose = 0, + atol = 1e-10, rtol = 1e-8, itmax = 10 * (kkt.n + kkt.m), + verbose = 0, ) x = Krylov.solution(kkt.linear_solver) + # Compute true residual ‖K·x − b‖ against the operator we just gave Krylov + Kx = similar(x) + LinearAlgebra.mul!(Kx, kkt, x, true, false) # routes through `_kkt_apply!` + true_res = LinearAlgebra.norm(Kx .- rhs_copy) + nit = Krylov.iteration_count(kkt.linear_solver) + issolved = Krylov.issolved(kkt.linear_solver) + kstatus = Krylov.statistics(kkt.linear_solver).status + println(stderr, + "Krylov solve: rhs=", rhs_norm, + " sol=", LinearAlgebra.norm(x), + " res=", true_res, + " it=", nit, + " solved=", issolved, + " ", kstatus, + ) copyto!(b, x) - push!(kkt.krylov_iterations, Krylov.iteration_count(kkt.linear_solver)) + push!(kkt.krylov_iterations, nit) push!(kkt.krylov_residuals, Krylov.elapsed_time(kkt.linear_solver)) MadNLP.finish_aug_solve!(kkt, w) return w From 9ebf9ea36c696ab03e444b3126f8fabf52a40dbc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 29 May 2026 16:47:32 +0200 Subject: [PATCH 06/15] Fix seed and use DEBUG --- bench/bm_madnlp_bench.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index 5c8d5bc5b..27f553a6d 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -18,6 +18,7 @@ import NLPModels import Dualization import LinearAlgebra import MathOptInterface as MOI +import Random include(joinpath(@__DIR__, "bm_madnlp_kkt.jl")) @@ -25,6 +26,8 @@ DynamicPolynomials.@polyvar x # Custom sub_solver factory that builds a MadNLPSolver with our BMKKTSystem. function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; kws...) + # Deterministic init so every run sees the same KKT trace. + Random.seed!(0) # Override the random `meta.x0` with a feasibility-friendlier starting # point. `BMSOSAL` uses `rand(n)` and `BMSOS` uses `randn(n) * 1e-8`; # `rand(n) / n` keeps norms small as the rank-factor count grows. @@ -41,7 +44,7 @@ function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; kws...) kkt_system = BMKKTSystem, linear_solver = MadNLP.LapackCPUSolver, # unused — Krylov takes over nlp_scaling = false, # `jac_dense!` not implemented - print_level = MadNLP.TRACE, + print_level = MadNLP.DEBUG, max_iter = 50, ) end From 5a58eebdba0cb2ea674840224f2ea72ea7d013e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 29 May 2026 17:04:24 +0200 Subject: [PATCH 07/15] Add qlp option --- bench/bm_madnlp_bench.jl | 4 ++-- bench/bm_madnlp_kkt.jl | 39 +++++++++++++++++++++++++-------------- 2 files changed, 27 insertions(+), 16 deletions(-) diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index 27f553a6d..57e9d9fa3 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -25,7 +25,7 @@ include(joinpath(@__DIR__, "bm_madnlp_kkt.jl")) DynamicPolynomials.@polyvar x # Custom sub_solver factory that builds a MadNLPSolver with our BMKKTSystem. -function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; kws...) +function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; qlp::Bool = true, kws...) # Deterministic init so every run sees the same KKT trace. Random.seed!(0) # Override the random `meta.x0` with a feasibility-friendlier starting @@ -41,7 +41,7 @@ function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; kws...) return MadNLP.MadNLPSolver( nlp; callback = MadNLP.DenseCallback, - kkt_system = BMKKTSystem, + kkt_system = BMKKTSystem{qlp}, linear_solver = MadNLP.LapackCPUSolver, # unused — Krylov takes over nlp_scaling = false, # `jac_dense!` not implemented print_level = MadNLP.DEBUG, diff --git a/bench/bm_madnlp_kkt.jl b/bench/bm_madnlp_kkt.jl index c108e8259..ba2ca1102 100644 --- a/bench/bm_madnlp_kkt.jl +++ b/bench/bm_madnlp_kkt.jl @@ -11,7 +11,7 @@ import NLPModels import LinearAlgebra import LowRankOpt as LRO -mutable struct BMKKTSystem{T,VT,NLP,LS} <: +mutable struct BMKKTSystem{QLP,T,VT,NLP,LS} <: MadNLP.AbstractReducedKKTSystem{T,VT,Matrix{T},MadNLP.ExactHessian{T,VT}} nlp::NLP n::Int @@ -45,11 +45,17 @@ mutable struct BMKKTSystem{T,VT,NLP,LS} <: jtv_buf::VT end -function BMKKTSystem(nlp; T = Float64, VT = Vector{T}) +# `qlp=true` → `Krylov.MinresQlpWorkspace` (robust on singular K, default). +# `qlp=false` → `Krylov.MinresWorkspace` (cheaper per iter, but falls back to +# the min-norm least-squares solution when K is rank-deficient — yields +# bogus Newton directions there). +function BMKKTSystem(nlp; qlp::Bool = true, T = Float64, VT = Vector{T}) n = NLPModels.get_nvar(nlp) m = NLPModels.get_ncon(nlp) - workspace = Krylov.MinresWorkspace(n + m, n + m, VT) - return BMKKTSystem{T,VT,typeof(nlp),typeof(workspace)}( + workspace = qlp ? + Krylov.MinresQlpWorkspace(n + m, n + m, VT) : + Krylov.MinresWorkspace(n + m, n + m, VT) + return BMKKTSystem{qlp,T,VT,typeof(nlp),typeof(workspace)}( nlp, n, m, VT(undef, n), # reg VT(undef, n), # pr_diag @@ -72,38 +78,43 @@ function BMKKTSystem(nlp; T = Float64, VT = Vector{T}) end function MadNLP.create_kkt_system( - ::Type{<:BMKKTSystem}, + ::Type{BMKKTSystem{QLP}}, cb::MadNLP.AbstractCallback{T,VT}, linear_solver::Type; opt_linear_solver = MadNLP.default_options(linear_solver), hessian_approximation = MadNLP.ExactHessian, qn_options = MadNLP.QuasiNewtonOptions(), -) where {T,VT} +) where {QLP,T,VT} # Only supported with `square_scalars=true` (no bounds, equality # constraints only). Verify here so misconfiguration surfaces early. @assert isempty(cb.ind_ineq) "BMKKTSystem assumes equality-only constraints" @assert isempty(cb.ind_lb) "BMKKTSystem assumes no lower-bound variables (square_scalars=true)" @assert isempty(cb.ind_ub) "BMKKTSystem assumes no upper-bound variables (square_scalars=true)" - return BMKKTSystem(cb.nlp; T = T, VT = VT) + return BMKKTSystem(cb.nlp; qlp = QLP, T = T, VT = VT) end +# Convenience dispatch — `kkt_system = BMKKTSystem` (no type param) ↔ QLP=true. +MadNLP.create_kkt_system(::Type{BMKKTSystem}, cb, ls; kws...) = + MadNLP.create_kkt_system(BMKKTSystem{true}, cb, ls; kws...) MadNLP.num_variables(kkt::BMKKTSystem) = kkt.n MadNLP.get_hessian(::BMKKTSystem) = nothing MadNLP.get_jacobian(::BMKKTSystem) = nothing # Pretend the Krylov workspace is a "factorization" with no inertia info. -MadNLP.is_inertia(::Krylov.MinresWorkspace) = true -MadNLP.inertia(::Krylov.MinresWorkspace) = (0, 0, 0) -MadNLP.introduce(::Krylov.MinresWorkspace) = "Krylov.MINRES" -MadNLP.improve!(::Krylov.MinresWorkspace) = true -MadNLP.factorize!(::Krylov.MinresWorkspace) = nothing +const _KrylovWS = Union{Krylov.MinresWorkspace,Krylov.MinresQlpWorkspace} +MadNLP.is_inertia(::_KrylovWS) = true +MadNLP.inertia(::_KrylovWS) = (0, 0, 0) +MadNLP.introduce(::Krylov.MinresQlpWorkspace) = "Krylov.MINRES-QLP" +MadNLP.introduce(::Krylov.MinresWorkspace) = "Krylov.MINRES" +MadNLP.improve!(::_KrylovWS) = true +MadNLP.factorize!(::_KrylovWS) = nothing MadNLP.is_inertia_correct(::BMKKTSystem, _, _, _) = true -Base.eltype(::BMKKTSystem{T}) where {T} = T +Base.eltype(::BMKKTSystem{QLP,T}) where {QLP,T} = T Base.size(kkt::BMKKTSystem) = (kkt.n + kkt.m, kkt.n + kkt.m) Base.size(kkt::BMKKTSystem, ::Int) = kkt.n + kkt.m -function MadNLP.initialize!(kkt::BMKKTSystem{T}) where {T} +function MadNLP.initialize!(kkt::BMKKTSystem{QLP,T}) where {QLP,T} fill!(kkt.reg, one(T)) fill!(kkt.pr_diag, one(T)) fill!(kkt.du_diag, zero(T)) From 80e917223c2dcc30745fbccf5105e6a427aaeffe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 29 May 2026 17:37:49 +0200 Subject: [PATCH 08/15] Fixes MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- bench/bm_madnlp_bench.jl | 7 ++++++- bench/bm_madnlp_kkt.jl | 33 +++++++++++++++++++++++++++------ 2 files changed, 33 insertions(+), 7 deletions(-) diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index 57e9d9fa3..a7e63e00f 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -43,7 +43,12 @@ function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; qlp::Bool = true, k callback = MadNLP.DenseCallback, kkt_system = BMKKTSystem{qlp}, linear_solver = MadNLP.LapackCPUSolver, # unused — Krylov takes over - nlp_scaling = false, # `jac_dense!` not implemented + nlp_scaling = false, # `jac_dense!` not implemented + # We can't honestly report KKT inertia from a matrix-free MINRES + # solve, so use MadNLP's inertia-free correction (Curtis–Schenk–Wächter + # heuristic). It just needs two extra `solve_kkt!`s per Newton step + # to validate that the computed direction is a descent step. + inertia_correction_method = MadNLP.InertiaFree, print_level = MadNLP.DEBUG, max_iter = 50, ) diff --git a/bench/bm_madnlp_kkt.jl b/bench/bm_madnlp_kkt.jl index ba2ca1102..3ce874bbc 100644 --- a/bench/bm_madnlp_kkt.jl +++ b/bench/bm_madnlp_kkt.jl @@ -206,6 +206,10 @@ function MadNLP.mul!( end # Variant called by `Krylov.kmul!` on plain `Vector{T}` (length `n + m`). +# Must mirror the `AbstractKKTVector` variant — including the `reg`/`du_diag` +# contributions that the sparse path applies via `_kktmul!`. Without these, +# MadNLP's `del_w` regularization never reaches Krylov, the curvature test +# can never be satisfied, and the IPM bails into restoration. function LinearAlgebra.mul!( y::AbstractVector, kkt::BMKKTSystem, @@ -213,18 +217,35 @@ function LinearAlgebra.mul!( alpha::Number, beta::Number, ) n, m = kkt.n, kkt.m - _kkt_apply!( - view(y, 1:n), view(y, (n+1):(n+m)), - kkt, - view(x, 1:n), view(x, (n+1):(n+m)), - alpha, beta, - ) + yp = view(y, 1:n); yd = view(y, (n+1):(n+m)) + xp = view(x, 1:n); xd = view(x, (n+1):(n+m)) + _kkt_apply!(yp, yd, kkt, xp, xd, alpha, beta) + # Match `_kktmul!` from the `AbstractKKTVector` path: add `reg` to the + # primal block and `du_diag` to the dual block. These reflect MadNLP's + # cumulative `del_w`/`del_c` Hessian/Jacobian perturbations. + yp .+= alpha .* kkt.reg .* xp + yd .+= alpha .* kkt.du_diag .* xd return y end LinearAlgebra.mul!(y::AbstractVector, kkt::BMKKTSystem, x::AbstractVector) = LinearAlgebra.mul!(y, kkt, x, true, false) +# Hessian-block matvec used by `InertiaFree`'s curvature test (`curv_test` +# in `MadNLP/src/IPM/solver.jl:785`). Computes `wx = (H + pr_diag) · t` +# where `t`, `wx` ∈ ℝⁿ are *primal only* (see `build_inertia_corrector`). +function MadNLP.mul_hess_blk!( + wx::AbstractVector, kkt::BMKKTSystem, t::AbstractVector, +) + NLPModels.hprod!( + kkt.nlp, kkt.current_x, kkt.current_y, t, kkt.hv_buf; + obj_weight = one(eltype(wx)), + ) + copyto!(wx, kkt.hv_buf) + wx .+= t .* kkt.pr_diag + return wx +end + # `Aᵀ x`, just routes to the NLP's `jtprod!`. function MadNLP.jtprod!(y::AbstractVector, kkt::BMKKTSystem, x::AbstractVector) NLPModels.jtprod!(kkt.nlp, y, x, y) From 1993836a2ec5f53742cd3c70beac721be977b512 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 29 May 2026 18:22:20 +0200 Subject: [PATCH 09/15] Fix sign convention --- bench/Project.toml | 1 + bench/bm_madnlp_bench.jl | 60 ++++++++++++++++++++++++++++++++-------- bench/bm_madnlp_kkt.jl | 17 ++++++++++-- 3 files changed, 63 insertions(+), 15 deletions(-) diff --git a/bench/Project.toml b/bench/Project.toml index c2cb22b50..168ca4ba9 100644 --- a/bench/Project.toml +++ b/bench/Project.toml @@ -9,5 +9,6 @@ MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" +Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index a7e63e00f..a5775894c 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -49,8 +49,8 @@ function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; qlp::Bool = true, k # heuristic). It just needs two extra `solve_kkt!`s per Newton step # to validate that the computed direction is a descent step. inertia_correction_method = MadNLP.InertiaFree, - print_level = MadNLP.DEBUG, - max_iter = 50, + print_level = MadNLP.INFO, + max_iter = 500, ) end @@ -72,17 +72,53 @@ function with_lro_bridges!(model) ) end -println("== BMKKTSystem smoke: max γ s.t. p − γ ∈ SOS ==") -model = Model(bmlbfgs) -set_silent(model) -@variable(model, γ) -@objective(model, Max, γ) -with_lro_bridges!(model) p = x^4 - 4x^3 - 2x^2 + 12x + 3 -@constraint(model, p - γ in SOSCone(), zero_basis = MB.BoxSampling([-1.0], [1.0])) -optimize!(model) -println("primal_status = ", primal_status(model)) -println("value(γ) = ", value(γ), " (expected ≈ -6)") + +function smoke(feas::Bool) + println("\n== BMKKTSystem smoke: ", feas ? "feasibility (γ=-6)" : "max γ", " ==") + model = Model(bmlbfgs) + set_silent(model) + γ = -6 + if !feas + @variable(model, γ) + @objective(model, Max, γ) + end + with_lro_bridges!(model) + @constraint(model, p - γ in SOSCone(), zero_basis = MB.BoxSampling([-1.0], [1.0])) + optimize!(model) + println("primal_status = ", primal_status(model)) + if !feas + println("value(γ) = ", value(γ), " (expected ≈ -6)") + end +end + +smoke(true) +smoke(false) + +# Cross-check : same problem, but with Percival as the sub_solver. If Percival +# *also* diverges on `max γ`, the bug is upstream (Dualization / bridges), +# not in our MadNLP integration. +println("\n== Percival cross-check (max γ) ==") +import Percival +percival_inner = optimizer_with_attributes( + LRO.Optimizer, + "solver" => LRO.BurerMonteiro.Solver, + "sub_solver" => Percival.PercivalSolver, + "ranks" => [4], + "square_scalars" => true, +) +percival_bmlbfgs = Dualization.dual_optimizer(percival_inner; assume_min_if_feasibility = true) +let + model = Model(percival_bmlbfgs) + set_silent(model) + @variable(model, γ) + @objective(model, Max, γ) + with_lro_bridges!(model) + @constraint(model, p - γ in SOSCone(), zero_basis = MB.BoxSampling([-1.0], [1.0])) + optimize!(model) + println("primal_status = ", primal_status(model)) + println("value(γ) = ", value(γ), " (expected ≈ -6)") +end # Should give # == BMKKTSystem smoke: max γ s.t. p − γ ∈ SOS == diff --git a/bench/bm_madnlp_kkt.jl b/bench/bm_madnlp_kkt.jl index 3ce874bbc..9bf555e31 100644 --- a/bench/bm_madnlp_kkt.jl +++ b/bench/bm_madnlp_kkt.jl @@ -138,7 +138,12 @@ function MadNLP.eval_lag_hess_wrapper!( l::AbstractVector; is_resto = false, ) - copyto!(kkt.current_y, l) + # MadNLP convention: `L(x, y) = obj_weight·f(x) + y'·c(x)` (dual feasibility + # `∇f + Jᵀy = 0`, see `IPM/kernels.jl:247`). + # NLPModels convention: `L(x, y) = obj_weight·f(x) − y'·c(x)`, so its + # `hprod!` returns `(obj_weight·∇²f − Σ y_i ∇²c_i)·v`. + # To get MadNLP's Hessian `∇²f + Σ y_i ∇²c_i` we must pass `−l`. + kkt.current_y .= .-l return end @@ -246,9 +251,15 @@ function MadNLP.mul_hess_blk!( return wx end -# `Aᵀ x`, just routes to the NLP's `jtprod!`. +# `Aᵀ x`, linearized at the current IPM iterate (stashed in `kkt.current_x` +# by our `eval_jac_wrapper!`). Until that fix, this passed `y` as the +# evaluation point of the Jacobian — fine if `c` is linear, but `c(X) = +# jprod(model, X, X)` for our BM model is *quadratic* in `X`, so J depends +# on the linearization point and the wrong one produced bogus duals once +# `y` had non-trivial magnitude (feasibility happens to keep `y ≈ 0` so the +# bug stayed dormant; `max γ` revealed it). function MadNLP.jtprod!(y::AbstractVector, kkt::BMKKTSystem, x::AbstractVector) - NLPModels.jtprod!(kkt.nlp, y, x, y) + NLPModels.jtprod!(kkt.nlp, kkt.current_x, x, y) return y end From 186eccfbfd6e2a1eb6f3373b5b9dfba56614a030 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Fri, 29 May 2026 19:04:02 +0200 Subject: [PATCH 10/15] tweak --- bench/bm_madnlp_bench.jl | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index a5775894c..c3567d355 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -29,10 +29,12 @@ function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; qlp::Bool = true, k # Deterministic init so every run sees the same KKT trace. Random.seed!(0) # Override the random `meta.x0` with a feasibility-friendlier starting - # point. `BMSOSAL` uses `rand(n)` and `BMSOS` uses `randn(n) * 1e-8`; - # `rand(n) / n` keeps norms small as the rank-factor count grows. + # point. Starting from `rand(n)/n` (small) makes the Hessian near-singular + # at iter 0, which gives Newton steps that overshoot the optimum and the + # IPM oscillates. Use `rand(n)` (BMSOSAL's choice) for a non-degenerate + # initial Hessian. n = length(nlp.meta.x0) - nlp.meta.x0 .= rand(eltype(nlp.meta.x0), n) ./ n + nlp.meta.x0 .= rand(eltype(nlp.meta.x0), n) # Quick check whether the BM objective is identically zero (which would # make the dual-Newton step purely feasibility-driven). f0 = NLPModels.obj(nlp, nlp.meta.x0) @@ -49,6 +51,18 @@ function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; qlp::Bool = true, k # heuristic). It just needs two extra `solve_kkt!`s per Newton step # to validate that the computed direction is a descent step. inertia_correction_method = MadNLP.InertiaFree, + # Stabilisation : la `curv_test` par défaut tolère ζ=0 (curvature ≥ 0), + # ce qui laisse passer des directions presque-singulières → pas de + # Newton qui dépassent l'optimum. Demander une curvature strictement + # positive et un régulariseur primal de base force un comportement + # type «Levenberg-Marquardt léger». + inertia_free_tol = 1e-8, + default_primal_regularization = 1e-8, + default_dual_regularization = 1e-8, + # Default least-squares dual init solves K·[d;y] = [-g;-c] once, + # which on our singular KKT can throw `y` to wild values. Starting + # from `y = 0` lets MadNLP build up the duals iteration by iteration. + dual_initialization_method = MadNLP.DualInitializeSetZero, print_level = MadNLP.INFO, max_iter = 500, ) From e0e5503bbfb253795d2fffa4a2bc924a0197f00d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 30 May 2026 07:52:32 +0200 Subject: [PATCH 11/15] Add preconditioner --- bench/bm_madnlp_kkt.jl | 48 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/bench/bm_madnlp_kkt.jl b/bench/bm_madnlp_kkt.jl index 9bf555e31..d35514bfd 100644 --- a/bench/bm_madnlp_kkt.jl +++ b/bench/bm_madnlp_kkt.jl @@ -43,7 +43,24 @@ mutable struct BMKKTSystem{QLP,T,VT,NLP,LS} <: hv_buf::VT jv_buf::VT jtv_buf::VT + # Diagonal preconditioner for Krylov. Length `n+m`. Refreshed before + # each `krylov_solve!` from `reg + pr_diag + |diag(H)|` on the primal + # block and `max(|du_diag|, 1)` on the dual block. `diag(H)` is + # extracted by probing `hprod!(e_i)`; cheap for our smoke-test sizes + # but `O(n)` FFTs per IPM iter for the larger benchmark. + precond_diag::VT + precond_probe::VT # buffer for the probe vector end +# Wrapper so Krylov treats `precond_diag` as a left preconditioner via +# `LinearAlgebra.ldiv!` (we set `ldiv = true` in `solve_kkt!`). +struct BMJacobiPrecond{VT} <: AbstractMatrix{Float64} + diag::VT +end +Base.size(M::BMJacobiPrecond) = (length(M.diag), length(M.diag)) +Base.eltype(::BMJacobiPrecond{VT}) where {VT} = eltype(VT) +LinearAlgebra.ldiv!(y::AbstractVector, M::BMJacobiPrecond, x::AbstractVector) = + (y .= x ./ M.diag; return y) +LinearAlgebra.ldiv!(M::BMJacobiPrecond, x::AbstractVector) = (x ./= M.diag; return x) # `qlp=true` → `Krylov.MinresQlpWorkspace` (robust on singular K, default). # `qlp=false` → `Krylov.MinresWorkspace` (cheaper per iter, but falls back to @@ -74,6 +91,8 @@ function BMKKTSystem(nlp; qlp::Bool = true, T = Float64, VT = Vector{T}) VT(undef, n), # hv_buf VT(undef, m), # jv_buf VT(undef, n), # jtv_buf + ones(T, n + m), # precond_diag + zeros(T, n), # precond_probe ) end @@ -320,6 +339,32 @@ function SolverCore.solve!( return stats end +# Probe `diag(H)` via `n` Hessian-vector products on standard basis vectors. +# Cheap for our smoke-test (`n = 14`); for the trigonometric `d = 100` +# benchmark this is `n ≈ 800` FFTs per IPM iter, still negligible compared +# to the Krylov inner loop. +function _refresh_preconditioner!(kkt::BMKKTSystem{QLP,T}) where {QLP,T} + n, m = kkt.n, kkt.m + e = kkt.precond_probe + fill!(e, zero(T)) + for i in 1:n + e[i] = one(T) + NLPModels.hprod!( + kkt.nlp, kkt.current_x, kkt.current_y, e, kkt.hv_buf; + obj_weight = one(T), + ) + # `|diag(H)| + reg + pr_diag`, guarded away from zero. Taking abs + # keeps the preconditioner SPD even where H has indefinite diagonal + # entries — for MINRES-QLP a SPD preconditioner is required. + kkt.precond_diag[i] = max(abs(kkt.hv_buf[i]) + kkt.reg[i] + kkt.pr_diag[i], 1e-12) + e[i] = zero(T) + end + for j in 1:m + kkt.precond_diag[n + j] = max(abs(kkt.du_diag[j]), one(T)) + end + return +end + function MadNLP.solve_kkt!(kkt::BMKKTSystem, w::MadNLP.AbstractKKTVector) MadNLP.reduce_rhs!(kkt, w) # The `[primal; dual]` block — what MadNLP's standard @@ -327,8 +372,11 @@ function MadNLP.solve_kkt!(kkt::BMKKTSystem, w::MadNLP.AbstractKKTVector) b = MadNLP.primal_dual(w) rhs_norm = LinearAlgebra.norm(b) rhs_copy = copy(b) + _refresh_preconditioner!(kkt) + M = BMJacobiPrecond(kkt.precond_diag) Krylov.krylov_solve!( kkt.linear_solver, kkt, b; + M = M, ldiv = true, atol = 1e-10, rtol = 1e-8, itmax = 10 * (kkt.n + kkt.m), verbose = 0, ) From 8f567c9ea096fa3bd686e797e4457302a6043ad5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 30 May 2026 13:47:55 +0100 Subject: [PATCH 12/15] =?UTF-8?q?Pipeline=20complet=20fonctionnel,=20IPM?= =?UTF-8?q?=20converge=20=C3=A0=20l'optimum=20exact.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- bench/bm_madnlp_bench.jl | 5 +- bench/bm_madnlp_kkt.jl | 102 ++++++++++++++++++++++++++++++++++----- 2 files changed, 94 insertions(+), 13 deletions(-) diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index c3567d355..aded567e6 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -50,7 +50,10 @@ function MadNLPMinresSolver(nlp::NLPModels.AbstractNLPModel; qlp::Bool = true, k # solve, so use MadNLP's inertia-free correction (Curtis–Schenk–Wächter # heuristic). It just needs two extra `solve_kkt!`s per Newton step # to validate that the computed direction is a descent step. - inertia_correction_method = MadNLP.InertiaFree, + # Lanczos probe in `factorize_wrapper!` reports honest inertia, so + # MadNLP can drive `del_w` on its own — `InertiaBased` is now the + # right choice (no more `InertiaFree` curvature heuristic). + inertia_correction_method = MadNLP.InertiaBased, # Stabilisation : la `curv_test` par défaut tolère ζ=0 (curvature ≥ 0), # ce qui laisse passer des directions presque-singulières → pas de # Newton qui dépassent l'optimum. Demander une curvature strictement diff --git a/bench/bm_madnlp_kkt.jl b/bench/bm_madnlp_kkt.jl index d35514bfd..9f4cec2c3 100644 --- a/bench/bm_madnlp_kkt.jl +++ b/bench/bm_madnlp_kkt.jl @@ -50,6 +50,17 @@ mutable struct BMKKTSystem{QLP,T,VT,NLP,LS} <: # but `O(n)` FFTs per IPM iter for the larger benchmark. precond_diag::VT precond_probe::VT # buffer for the probe vector + # Inertia estimate produced by a short Lanczos probe of `K` before each + # solve; `MadNLP.inertia(workspace)` reads this so `InertiaBased`'s + # `is_inertia_correct(n_pos, n_zero, n_neg)` can fire and trigger a + # `del_w` bump when the operator is indefinite on the augmented system. + inertia_pos::Base.RefValue{Int} + inertia_zero::Base.RefValue{Int} + inertia_neg::Base.RefValue{Int} + # Lanczos probe buffers (length `n + m`). + lanczos_v_prev::VT + lanczos_v_curr::VT + lanczos_w::VT end # Wrapper so Krylov treats `precond_diag` as a left preconditioner via # `LinearAlgebra.ldiv!` (we set `ldiv = true` in `solve_kkt!`). @@ -93,6 +104,12 @@ function BMKKTSystem(nlp; qlp::Bool = true, T = Float64, VT = Vector{T}) VT(undef, n), # jtv_buf ones(T, n + m), # precond_diag zeros(T, n), # precond_probe + Ref(n), # inertia_pos + Ref(0), # inertia_zero + Ref(m), # inertia_neg + zeros(T, n + m), # lanczos_v_prev + zeros(T, n + m), # lanczos_v_curr + zeros(T, n + m), # lanczos_w ) end @@ -119,7 +136,10 @@ MadNLP.num_variables(kkt::BMKKTSystem) = kkt.n MadNLP.get_hessian(::BMKKTSystem) = nothing MadNLP.get_jacobian(::BMKKTSystem) = nothing -# Pretend the Krylov workspace is a "factorization" with no inertia info. +# Krylov workspaces themselves carry no inertia. The real estimate lives on +# `BMKKTSystem` (refreshed by `_inertia_probe!` before each `solve_kkt!`); +# the workspace shim below stays a no-op so `is_inertia_correct(::BMKKTSystem, +# ...)` reads the kkt-level refs. const _KrylovWS = Union{Krylov.MinresWorkspace,Krylov.MinresQlpWorkspace} MadNLP.is_inertia(::_KrylovWS) = true MadNLP.inertia(::_KrylovWS) = (0, 0, 0) @@ -127,7 +147,17 @@ MadNLP.introduce(::Krylov.MinresQlpWorkspace) = "Krylov.MINRES-QLP" MadNLP.introduce(::Krylov.MinresWorkspace) = "Krylov.MINRES" MadNLP.improve!(::_KrylovWS) = true MadNLP.factorize!(::_KrylovWS) = nothing -MadNLP.is_inertia_correct(::BMKKTSystem, _, _, _) = true +# MadNLP's `InertiaBased` corrector calls `inertia(kkt.linear_solver)` to +# get the inertia tuple, *then* hands it to `is_inertia_correct(kkt, ...)`. +# Our `Krylov.MinresQlpWorkspace` shim is forced to return `(0,0,0)` since +# the workspace has no view of `kkt`. So we ignore the passed args and read +# the inertia stashed on `kkt` by `_inertia_probe!`. Same for +# `should_regularize_dual`. +MadNLP.is_inertia_correct(kkt::BMKKTSystem, _, _, _) = + kkt.inertia_pos[] == kkt.n && + kkt.inertia_zero[] == 0 && + kkt.inertia_neg[] == kkt.m +MadNLP.should_regularize_dual(kkt::BMKKTSystem, _, _, _) = kkt.inertia_zero[] != 0 Base.eltype(::BMKKTSystem{QLP,T}) where {QLP,T} = T Base.size(kkt::BMKKTSystem) = (kkt.n + kkt.m, kkt.n + kkt.m) @@ -173,6 +203,10 @@ function MadNLP.factorize_wrapper!( solver::MadNLP.MadNLPSolver{T,VT,IT,KKT}, ) where {T,VT,IT,KKT<:BMKKTSystem} MadNLP.build_kkt!(solver.kkt) + # Lanczos inertia probe runs against the regularization-augmented matvec, + # so MadNLP's subsequent `inertia(kkt)` / `is_inertia_correct(kkt, ...)` + # check honestly reflects the current `reg`/`pr_diag`/`du_diag` state. + _inertia_probe!(solver.kkt, solver.kkt.n + solver.kkt.m) return true end @@ -339,6 +373,60 @@ function SolverCore.solve!( return stats end +# Dense reconstruction of `K` via `n+m` standard-basis matvecs, then +# `eigvals`. Cheap for the smoke (`n+m = 19`); the trig benchmark +# (`n+m ≈ 1000+`) will need a Lanczos with full reorthogonalization +# instead. Plain three-term Lanczos on indefinite KKT loses orthogonality +# after a few iters and produces spurious eigenvalues, which made +# `is_inertia_correct` permanently false. +function _inertia_probe!(kkt::BMKKTSystem{QLP,T}, k_max::Int) where {QLP,T} + nm = kkt.n + kkt.m + Kmat = zeros(T, nm, nm) + e = kkt.lanczos_v_curr + w = kkt.lanczos_w + fill!(e, zero(T)) + for j in 1:nm + e[j] = one(T) + LinearAlgebra.mul!(w, kkt, e) + @views Kmat[:, j] .= w + e[j] = zero(T) + end + # Symmetrize defensively (drops asymmetry from accumulated FP error). + Kmat .= (Kmat .+ Kmat') ./ 2 + # Bunch-Kaufman on the symmetric `Kmat` — the signs of `D`'s diagonal + # give exact inertia by Sylvester's law, without an eigenvalue tolerance + # to tune. (The 2×2 blocks have one positive and one negative eigenvalue.) + bk = LinearAlgebra.bunchkaufman(LinearAlgebra.Symmetric(Kmat); check = false) + D = bk.D # `Tridiagonal` view on the block-diagonal D + n_pos = 0; n_neg = 0; n_zero = 0 + i = 1 + while i ≤ nm + if i < nm && D[i + 1, i] != zero(T) + # 2×2 block: one positive, one negative eigenvalue + n_pos += 1; n_neg += 1 + i += 2 + else + d = D[i, i] + if d > zero(T); n_pos += 1 + elseif d < zero(T); n_neg += 1 + else; n_zero += 1 + end + i += 1 + end + end + kkt.inertia_pos[] = n_pos + kkt.inertia_zero[] = n_zero + kkt.inertia_neg[] = n_neg + return +end + +# Hook MadNLP reads after `factorize_wrapper!`. We computed the inertia +# during `factorize!` (overridden below) so just return the cached values. +function MadNLP.inertia(kkt::BMKKTSystem) + return (kkt.inertia_pos[], kkt.inertia_zero[], kkt.inertia_neg[]) +end +MadNLP.is_inertia(::BMKKTSystem) = true + # Probe `diag(H)` via `n` Hessian-vector products on standard basis vectors. # Cheap for our smoke-test (`n = 14`); for the trigonometric `d = 100` # benchmark this is `n ≈ 800` FFTs per IPM iter, still negligible compared @@ -386,16 +474,6 @@ function MadNLP.solve_kkt!(kkt::BMKKTSystem, w::MadNLP.AbstractKKTVector) LinearAlgebra.mul!(Kx, kkt, x, true, false) # routes through `_kkt_apply!` true_res = LinearAlgebra.norm(Kx .- rhs_copy) nit = Krylov.iteration_count(kkt.linear_solver) - issolved = Krylov.issolved(kkt.linear_solver) - kstatus = Krylov.statistics(kkt.linear_solver).status - println(stderr, - "Krylov solve: rhs=", rhs_norm, - " sol=", LinearAlgebra.norm(x), - " res=", true_res, - " it=", nit, - " solved=", issolved, - " ", kstatus, - ) copyto!(b, x) push!(kkt.krylov_iterations, nit) push!(kkt.krylov_residuals, Krylov.elapsed_time(kkt.linear_solver)) From 614df23de8abede42cd7c13a9749aa5704553580 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 30 May 2026 15:08:55 +0100 Subject: [PATCH 13/15] Larges scale MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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. --- bench/bm_madnlp_bench.jl | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index aded567e6..a7f894147 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -137,6 +137,38 @@ let println("value(γ) = ", value(γ), " (expected ≈ -6)") end +# Bump the BM rank to scale up the primal-variable count without changing +# the optimum (γ* = -6) or the constraint structure. Smoke uses `rank=4` +# (`n + m = 19`); we sweep up to `rank = 32` (`n + m ≈ 103`). Tests the +# dense `n+m × n+m` Bunch-Kaufman inertia probe at a non-trivial size. +# (Polynomial-degree scaling hits an unrelated upstream `BoundsError` in +# `MultivariateBases.eval_basis!` for `degree > 4`; bumping rank exercises +# the same `n + m` growth without touching the SOS-bridge path.) +function smoke_box_with_rank(rank::Int) + inner_r = optimizer_with_attributes( + LRO.Optimizer, + "solver" => LRO.BurerMonteiro.Solver, + "sub_solver" => MadNLPMinresSolver, + "ranks" => [rank], + "square_scalars" => true, + ) + opt = Dualization.dual_optimizer(inner_r; assume_min_if_feasibility = true) + model = Model(opt) + set_silent(model) + @variable(model, γ) + @objective(model, Max, γ) + with_lro_bridges!(model) + @constraint(model, p - γ in SOSCone(), zero_basis = MB.BoxSampling([-1.0], [1.0])) + t = @elapsed optimize!(model) + println("rank = ", rank, " value(γ) = ", round(value(γ), digits = 8), + " (expected ≈ -6) time = ", round(t, digits = 2), " s") + return t +end +println("\n== Scale-up via BM rank (smoke poly, expected γ = -6) ==") +for r in (4, 8, 16, 32) + smoke_box_with_rank(r) +end + # Should give # == BMKKTSystem smoke: max γ s.t. p − γ ∈ SOS == # ┌ Info: BM model at x0 From dfe5fceaf982f5ecfb98c3afe75fcbe50af0dffc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Sat, 30 May 2026 16:15:45 +0100 Subject: [PATCH 14/15] Use fft --- bench/Project.toml | 1 + bench/bm_madnlp_bench.jl | 53 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 54 insertions(+) diff --git a/bench/Project.toml b/bench/Project.toml index 168ca4ba9..8ad3c28b6 100644 --- a/bench/Project.toml +++ b/bench/Project.toml @@ -11,4 +11,5 @@ MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" +StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1" diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index a7f894147..cb31d36b7 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -19,6 +19,7 @@ import Dualization import LinearAlgebra import MathOptInterface as MOI import Random +import StarAlgebras as SA include(joinpath(@__DIR__, "bm_madnlp_kkt.jl")) @@ -169,6 +170,58 @@ for r in (4, 8, 16, 32) smoke_box_with_rank(r) end +# Convert the smoke polynomial `3 + 12x − 2x² − 4x³ + x⁴` to a cos-only +# trig polynomial via Chebyshev as intermediate. Under `x = cos(θ)` we have +# `Tₖ(cos θ) = cos(kθ)`, so the Chebyshev expansion immediately gives the +# Fourier-cosine coefficients of `p(cos θ)`. Sampling `θ ∈ [−π, π]` covers +# `x ∈ [−1, 1]`, so the smoke optimum `γ* = −6` (attained at `x = −1`, i.e. +# `θ = π`) is preserved. +function _smoke_as_cos_only_trig() + mon_coeffs = Float64[3, 12, -2, -4, 1] # constant, x, x², x³, x⁴ + mon_sub = MB.SubBasis{MB.Monomial}(DynamicPolynomials.monomials(x, 0:4)) + cheb_full = MB.FullBasis{MB.Chebyshev}([x]) + cheb_sparse = SA.coeffs(mon_coeffs, mon_sub, cheb_full) + # `cheb_sparse` is `SparseCoefficients` over the chebyshev keys. Pack + # values into the trig basis (which interleaves cos/sin: index + # `2k − 1` → `cos(kθ)`, `2k` → `sin(kθ)`, `0` → constant). + cheb_vals = collect(SA.values(cheb_sparse)) + # The conversion may produce fewer than 5 chebyshev terms if some are + # zero; we know the smoke poly has nonzero `T_0, …, T_4`. + @assert length(cheb_vals) == 5 "expected 5 chebyshev coefficients, got $(length(cheb_vals))" + trig_coeffs = zeros(9) # `monomials(x, 0:8)` → 9 trig basis fns + trig_coeffs[1] = cheb_vals[1] # constant + trig_coeffs[2] = cheb_vals[2] # cos(θ) + trig_coeffs[4] = cheb_vals[3] # cos(2θ) + trig_coeffs[6] = cheb_vals[4] # cos(3θ) + trig_coeffs[8] = cheb_vals[5] # cos(4θ) + trig_basis = MB.SubBasis{MB.Trigonometric}(DynamicPolynomials.monomials(x, 0:8)) + return MB.algebra_element(trig_coeffs, trig_basis) +end +function smoke_box_trig(solver, p_trig) + model = Model(solver) + set_silent(model) + @variable(model, γ) + @objective(model, Max, γ) + with_lro_bridges!(model) + # MultivariateBases' `Trigonometric` recurrence interprets the sample + # value as `cos(θ)` directly (see `recurrence_eval` in trigonometric.jl), + # so `BoxSampling([-1, 1])` covers the full `θ ∈ [-π, π]` and the + # Chebyshev→cos-only conversion exactly reproduces the smoke polynomial. + @constraint(model, p_trig - γ in SOSCone(), zero_basis = MB.BoxSampling([-1.0], [1.0])) + t = @elapsed optimize!(model) + println(" value(γ) = ", round(value(γ), digits = 6), + " time = ", round(t, digits = 2), " s") + return value(γ) +end +println("\n== Trigonometric smoke (smoke poly via Chebyshev → cos-only trig) ==") +p_trig = _smoke_as_cos_only_trig() +println("MadNLP+MINRES-QLP:") +γ_mad = smoke_box_trig(bmlbfgs, p_trig) +println("Percival:") +γ_perc = smoke_box_trig(percival_bmlbfgs, p_trig) +println("agreement: |Δγ| = ", round(abs(γ_mad - γ_perc), digits = 6), + " expected ≈ -6") + # Should give # == BMKKTSystem smoke: max γ s.t. p − γ ∈ SOS == # ┌ Info: BM model at x0 From fb10e2f241652d83eeefff85e82065e74ac42635 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 16 Jun 2026 14:50:37 +0200 Subject: [PATCH 15/15] up --- bench/bm_madnlp_bench.jl | 31 +++++++++++++++---------------- test/Project.toml | 7 +++++++ 2 files changed, 22 insertions(+), 16 deletions(-) diff --git a/bench/bm_madnlp_bench.jl b/bench/bm_madnlp_bench.jl index cb31d36b7..ebefaf637 100644 --- a/bench/bm_madnlp_bench.jl +++ b/bench/bm_madnlp_bench.jl @@ -176,25 +176,24 @@ end # Fourier-cosine coefficients of `p(cos θ)`. Sampling `θ ∈ [−π, π]` covers # `x ∈ [−1, 1]`, so the smoke optimum `γ* = −6` (attained at `x = −1`, i.e. # `θ = π`) is preserved. -function _smoke_as_cos_only_trig() - mon_coeffs = Float64[3, 12, -2, -4, 1] # constant, x, x², x³, x⁴ - mon_sub = MB.SubBasis{MB.Monomial}(DynamicPolynomials.monomials(x, 0:4)) +# Reduce to a degree-2 monomial polynomial: `x² − x − 1`. Min on `ℝ` +# is at `x = 1/2`, with `p(1/2) = −5/4 = −1.25`. Its `Monomial → Chebyshev` +# conversion gives `−0.5 T₀ − T₁ + 0.5 T₂`, so the cos-only trig form is +# `p_trig(θ) = −0.5 − cos(θ) + 0.5·cos(2θ)`, with the same minimum `−1.25` +# attained at `θ = π/3` (where `cos(θ) = 1/2`). Sticks to `monomials(x, 0:4)` +# → 5 trig basis functions (odd, dodges the upstream `isodd(n_coef)` assertion). +function _quad_as_cos_only_trig() + mon_coeffs = Float64[-1, -1, 1] # constant, x, x² + mon_sub = MB.SubBasis{MB.Monomial}(DynamicPolynomials.monomials(x, 0:2)) cheb_full = MB.FullBasis{MB.Chebyshev}([x]) cheb_sparse = SA.coeffs(mon_coeffs, mon_sub, cheb_full) - # `cheb_sparse` is `SparseCoefficients` over the chebyshev keys. Pack - # values into the trig basis (which interleaves cos/sin: index - # `2k − 1` → `cos(kθ)`, `2k` → `sin(kθ)`, `0` → constant). cheb_vals = collect(SA.values(cheb_sparse)) - # The conversion may produce fewer than 5 chebyshev terms if some are - # zero; we know the smoke poly has nonzero `T_0, …, T_4`. - @assert length(cheb_vals) == 5 "expected 5 chebyshev coefficients, got $(length(cheb_vals))" - trig_coeffs = zeros(9) # `monomials(x, 0:8)` → 9 trig basis fns + @assert length(cheb_vals) == 3 "expected 3 chebyshev coefficients, got $(length(cheb_vals))" + trig_coeffs = zeros(5) # `monomials(x, 0:4)` → 5 trig basis fns trig_coeffs[1] = cheb_vals[1] # constant trig_coeffs[2] = cheb_vals[2] # cos(θ) trig_coeffs[4] = cheb_vals[3] # cos(2θ) - trig_coeffs[6] = cheb_vals[4] # cos(3θ) - trig_coeffs[8] = cheb_vals[5] # cos(4θ) - trig_basis = MB.SubBasis{MB.Trigonometric}(DynamicPolynomials.monomials(x, 0:8)) + trig_basis = MB.SubBasis{MB.Trigonometric}(DynamicPolynomials.monomials(x, 0:4)) return MB.algebra_element(trig_coeffs, trig_basis) end function smoke_box_trig(solver, p_trig) @@ -213,14 +212,14 @@ function smoke_box_trig(solver, p_trig) " time = ", round(t, digits = 2), " s") return value(γ) end -println("\n== Trigonometric smoke (smoke poly via Chebyshev → cos-only trig) ==") -p_trig = _smoke_as_cos_only_trig() +println("\n== Trigonometric quad (x²−x−1 via Chebyshev → cos-only trig) ==") +p_trig = _quad_as_cos_only_trig() println("MadNLP+MINRES-QLP:") γ_mad = smoke_box_trig(bmlbfgs, p_trig) println("Percival:") γ_perc = smoke_box_trig(percival_bmlbfgs, p_trig) println("agreement: |Δγ| = ", round(abs(γ_mad - γ_perc), digits = 6), - " expected ≈ -6") + " expected ≈ -1.25") # Should give # == BMKKTSystem smoke: max γ s.t. p − γ ∈ SOS == diff --git a/test/Project.toml b/test/Project.toml index 27911dc28..e20264312 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -3,19 +3,26 @@ Clarabel = "61c947e1-3e6d-4ee4-985a-eec8c727bd6e" Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" Dualization = "191a621a-6537-11e9-281d-650236a99e60" DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07" +FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" +FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" Hypatia = "b99e6be6-89ff-11e8-14f8-45c827f4f8f2" JSOSolvers = "10dff2fc-5484-5881-a0e0-c90441020f8a" JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LowRankOpt = "607ca3ad-272e-43c8-bcbe-fc71b56c935c" +MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6" MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" MultivariateBases = "be282fd4-ad43-11e9-1d11-8bd9d7e43378" MultivariateMoments = "f4abf1af-0426-5881-a0da-e2f168889b5e" MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3" MutableArithmetics = "d8a4904e-b15c-11e9-3269-09a3773c0cb0" +NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" +NLPModelsTest = "7998695d-6960-4d3a-85c4-e1bceb8cd856" Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc" PolyJuMP = "ddf597a6-d67e-5340-b84c-e37d84115374" SemialgebraicSets = "8e049039-38e8-557d-ae3a-bc521ccf6204" +SolverCore = "ff4d7338-4cf1-434d-91df-b86cb86fb843" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StarAlgebras = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1" SumOfSquares = "4b9e565b-77fc-50a5-a571-1244f986bda1"