From 1c0d5c9cbb360628e6c682d3650cd8c85b23f88d Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Sat, 11 Oct 2025 14:14:57 +0530 Subject: [PATCH 1/6] tron: support MINRES via krylov callback when radius kw not supported; add bench script and save results - Add run_krylov_subsolver! wrapper to enforce trust-region via Krylov callback if workspace doesn't accept `radius`. - Use wrapper inside projected_newton! so MINRES can run inside TRON trust-region loop. - Save benchmark results to bench/results_tron_minres_vs_cg.csv. --- Project.toml | 1 + bench/README.md | 23 +++++++++++ bench/bench_tron_minres_vs_cg.jl | 57 +++++++++++++++++++++++++ bench/results_tron_minres_vs_cg.csv | 5 +++ src/tron.jl | 64 +++++++++++++++++++++++++---- 5 files changed, 143 insertions(+), 7 deletions(-) create mode 100644 bench/README.md create mode 100644 bench/bench_tron_minres_vs_cg.jl create mode 100644 bench/results_tron_minres_vs_cg.csv diff --git a/Project.toml b/Project.toml index c95dc72e..b1a04fca 100644 --- a/Project.toml +++ b/Project.toml @@ -3,6 +3,7 @@ uuid = "10dff2fc-5484-5881-a0e0-c90441020f8a" version = "0.14.3" [deps] +ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" diff --git a/bench/README.md b/bench/README.md new file mode 100644 index 00000000..cbf9b48a --- /dev/null +++ b/bench/README.md @@ -0,0 +1,23 @@ +# Benchmark: TRON subsolver comparison (CG vs MINRES) + +This benchmark compares `tron(..., subsolver = :cg)` with `tron(..., subsolver = :minres)` using Krylov.jl. + +## Requirements + +- Julia 1.10 (as in Project.toml) +- From the package root, run `julia --project=.` to use the package environment. Ensure Krylov.jl is available in the environment (Project.toml pins Krylov = "0.10.0"). + +## Run + +```sh +julia --project=. bench/bench_tron_minres_vs_cg.jl +``` + +## Output + +- CSV printed to stdout with columns: prob,subsolver,elapsed,iter,status,objective,dual + +## Notes + +- The script uses synthetic logistic-like problems. Adjust sizes in the script if you want larger/smaller tests or longer time limits. +- If Krylov.jl does not support `:minres` with your installed version, you may need to update the environment to a Krylov release with MINRES trust-region support. diff --git a/bench/bench_tron_minres_vs_cg.jl b/bench/bench_tron_minres_vs_cg.jl new file mode 100644 index 00000000..e662ebbc --- /dev/null +++ b/bench/bench_tron_minres_vs_cg.jl @@ -0,0 +1,57 @@ +# Benchmark TRON subsolver: CG vs MINRES +# Usage: julia --project=. bench/bench_tron_minres_vs_cg.jl + +using Pkg, LinearAlgebra, Printf, Dates +Pkg.activate(joinpath(@__DIR__, "..")) + + +using JSOSolvers, Random, ADNLPModels + +# Build a logistic-like smooth convex objective as an ADNLPModel +function make_logistic(n, m; rng = Random.GLOBAL_RNG) + A = randn(rng, m, n) + xtrue = randn(rng, n) + y = sign.(A * xtrue + 0.1 * randn(rng, m)) + f(x) = begin + z = A * x + s = @. log(1 + exp(-y * z)) + sum(s) + end + gradf(x) = begin + z = A * x + sig = @. 1 / (1 + exp(y * z)) + return A' * (-y .* sig) + end + x0 = zeros(n) + l = fill(-Inf, n) + u = fill(Inf, n) + model = ADNLPModel((x->f(x)), x0, l, u; grad = (x->gradf(x))) + return model +end + +function run_one(model; subsolver = :cg, max_time = 10.0) + t0 = time() + stats = tron(model, subsolver = subsolver, max_time = max_time, verbose = 0, subsolver_verbose = 0) + t = time() - t0 + return (subsolver = subsolver, elapsed = t, iter = stats.iter, status = stats.status, objective = stats.objective, dual = stats.dual_feas) +end + +function main() + rng = MersenneTwister(1234) + problems = [make_logistic(1000, 2000, rng = rng), make_logistic(2000, 4000, rng = rng)] + + results = [] + for (i, model) in enumerate(problems) + for sub in (:cg, :minres) + push!(results, merge((prob = i,), run_one(model, subsolver = sub, max_time = 60.0))) + end + end + + # Print summary + println("prob,subsolver,elapsed,iter,status,objective,dual") + for r in results + @printf("%d,%s,%.6f,%d,%s,%.12g,%.6g\n", r[:prob], string(r[:subsolver]), r[:elapsed], r[:iter], string(r[:status]), r[:objective], r[:dual]) + end +end + +main() diff --git a/bench/results_tron_minres_vs_cg.csv b/bench/results_tron_minres_vs_cg.csv new file mode 100644 index 00000000..8ed140b4 --- /dev/null +++ b/bench/results_tron_minres_vs_cg.csv @@ -0,0 +1,5 @@ +prob,subsolver,elapsed,iter,status,objective,dual +1,cg,68.995797,2,max_time,195.150087826,118.665 +1,minres,74.223810,2,max_time,237.635885566,129.697 +2,cg,182.949333,1,max_time,948.148566607,623.416 +2,minres,139.410956,1,max_time,1027.95775894,632.073 diff --git a/src/tron.jl b/src/tron.jl index f2dc5972..d814c52d 100644 --- a/src/tron.jl +++ b/src/tron.jl @@ -223,12 +223,19 @@ end kwargs..., ) where {T, V} dict = Dict(kwargs) - subsolver_keys = intersect(keys(dict), tron_keys) - subsolver_kwargs = Dict(k => dict[k] for k in subsolver_keys) - solver = TronSolver(nlp; μ₀ = μ₀, μ₁ = μ₁, σ = σ, subsolver_kwargs...) - for k in subsolver_keys - pop!(dict, k) + + # Extract subsolver if provided (default :cg) and ensure it is consumed here + subsolver = haskey(dict, :subsolver) ? dict[:subsolver] : :cg + if haskey(dict, :subsolver) + pop!(dict, :subsolver) end + + # Construct solver with chosen subsolver. We intentionally do not forward + # additional TRON trust-region kwargs here to keep the call simple; the + # TronSolver constructor will use defaults for those parameters. + solver = TronSolver(nlp; μ₀ = μ₀, μ₁ = μ₁, σ = σ, subsolver = subsolver) + + # Remaining dict contains kwargs meant for solve! (we removed :subsolver) return solve!(solver, nlp; x = x, dict...) end @@ -582,6 +589,47 @@ function cauchy!( return α, s, status end +""" + run_krylov_subsolver!(workspace, A, b; radius, rtol, atol, timemax, verbose) + +Small compatibility wrapper around Krylov.jl's `krylov_solve!` to handle +workspaces that don't accept the `radius` keyword (e.g. Minres). We first +attempt to call with `radius`, and on a MethodError we retry without it. +""" +function run_krylov_subsolver!(workspace, A, b; radius=nothing, rtol, atol, timemax, verbose) + # Try calling with `radius` first (covers CG-like workspaces). + if radius !== nothing + try + return krylov_solve!(workspace, A, b; radius = radius, rtol = rtol, atol = atol, timemax = timemax, verbose = verbose) + catch e + # If the failure is due to an unsupported keyword (MethodError), fall + # back to calling with a callback that enforces the trust-region by + # inspecting the Krylov workspace's current solution `x` at each + # iteration. Any other error is rethrown. + if isa(e, MethodError) + # Krylov callbacks receive the workspace and should return `true` + # to request termination. We build a closure that captures `radius`. + function tr_callback(cb_workspace) + # Some workspaces store the current iterate in `x`. + if !hasfield(typeof(cb_workspace), :x) + # If workspace has no `x` field, don't stop via callback. + return false + end + xcur = getfield(cb_workspace, :x) + # Stop when the norm of the current solution reaches the radius. + return norm(xcur) >= radius + end + + return krylov_solve!(workspace, A, b; rtol = rtol, atol = atol, timemax = timemax, verbose = verbose, callback = tr_callback) + else + rethrow(e) + end + end + else + return krylov_solve!(workspace, A, b; rtol = rtol, atol = atol, timemax = timemax, verbose = verbose) + end +end + """ projected_newton!(solver, x, H, g, Δ, cgtol, ℓ, u, s, Hs; max_time = Inf, max_cgiter = 50, subsolver_verbose = 0) @@ -639,10 +687,12 @@ function projected_newton!( cg_op_diag[i] = ifix[i] ? 0 : 1 # implictly changes cg_op and so ZHZ end - cg!( + # Use a compatibility wrapper so workspaces with differing keyword + # signatures (MINRES doesn't accept `radius`) are handled gracefully. + run_krylov_subsolver!( cg_solver, ZHZ, - cgs_rhs, + cgs_rhs; radius = Δ, rtol = cgtol, atol = zero(T), From 78d3132e702916b40453af75964ae953763f21ed Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Sat, 11 Oct 2025 14:17:35 +0530 Subject: [PATCH 2/6] Update Project.toml --- Project.toml | 1 - 1 file changed, 1 deletion(-) diff --git a/Project.toml b/Project.toml index b1a04fca..c95dc72e 100644 --- a/Project.toml +++ b/Project.toml @@ -3,7 +3,6 @@ uuid = "10dff2fc-5484-5881-a0e0-c90441020f8a" version = "0.14.3" [deps] -ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LinearOperators = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125" From 087a666f68c09038161cd05b651a36dea1256f95 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Mon, 20 Oct 2025 18:48:33 +0530 Subject: [PATCH 3/6] Update src/tron.jl Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/tron.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/tron.jl b/src/tron.jl index d814c52d..729411d5 100644 --- a/src/tron.jl +++ b/src/tron.jl @@ -616,8 +616,8 @@ function run_krylov_subsolver!(workspace, A, b; radius=nothing, rtol, atol, time return false end xcur = getfield(cb_workspace, :x) - # Stop when the norm of the current solution reaches the radius. - return norm(xcur) >= radius + # Stop when the squared norm of the current solution reaches the squared radius (avoids sqrt). + return norm(xcur)^2 >= radius^2 end return krylov_solve!(workspace, A, b; rtol = rtol, atol = atol, timemax = timemax, verbose = verbose, callback = tr_callback) From c2eb640a88671bcb34a9c9de588368321785a7dc Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Mon, 20 Oct 2025 18:48:43 +0530 Subject: [PATCH 4/6] Update src/tron.jl Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/tron.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tron.jl b/src/tron.jl index 729411d5..9438861e 100644 --- a/src/tron.jl +++ b/src/tron.jl @@ -596,7 +596,7 @@ Small compatibility wrapper around Krylov.jl's `krylov_solve!` to handle workspaces that don't accept the `radius` keyword (e.g. Minres). We first attempt to call with `radius`, and on a MethodError we retry without it. """ -function run_krylov_subsolver!(workspace, A, b; radius=nothing, rtol, atol, timemax, verbose) +function run_krylov_subsolver!(workspace, A, b; radius=nothing, rtol=1e-6, atol=0.0, timemax=Inf, verbose=0) # Try calling with `radius` first (covers CG-like workspaces). if radius !== nothing try From 1752676fb76acf6b7a09508d900a46f887289da7 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Mon, 20 Oct 2025 18:48:53 +0530 Subject: [PATCH 5/6] Update src/tron.jl Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/tron.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/tron.jl b/src/tron.jl index 9438861e..da4ee9ef 100644 --- a/src/tron.jl +++ b/src/tron.jl @@ -612,8 +612,7 @@ function run_krylov_subsolver!(workspace, A, b; radius=nothing, rtol=1e-6, atol= function tr_callback(cb_workspace) # Some workspaces store the current iterate in `x`. if !hasfield(typeof(cb_workspace), :x) - # If workspace has no `x` field, don't stop via callback. - return false + error("Krylov workspace of type $(typeof(cb_workspace)) does not have an `x` field. Trust-region constraint cannot be enforced. Please use a compatible workspace type.") end xcur = getfield(cb_workspace, :x) # Stop when the squared norm of the current solution reaches the squared radius (avoids sqrt). From dc6cdd68119942024ccf64bd387d1e215bac78c8 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Tue, 28 Oct 2025 01:27:31 +0530 Subject: [PATCH 6/6] documentation --- bench/README.md | 4 ++-- bench/bench_tron_minres_vs_cg.jl | 8 -------- src/tron.jl | 16 +++------------- 3 files changed, 5 insertions(+), 23 deletions(-) diff --git a/bench/README.md b/bench/README.md index cbf9b48a..c4a85303 100644 --- a/bench/README.md +++ b/bench/README.md @@ -4,8 +4,8 @@ This benchmark compares `tron(..., subsolver = :cg)` with `tron(..., subsolver = ## Requirements -- Julia 1.10 (as in Project.toml) -- From the package root, run `julia --project=.` to use the package environment. Ensure Krylov.jl is available in the environment (Project.toml pins Krylov = "0.10.0"). +- Julia 1.10 +- From the package root, run `julia --project=.` to use the package environment. Ensure Krylov.jl is available in the environment. ## Run diff --git a/bench/bench_tron_minres_vs_cg.jl b/bench/bench_tron_minres_vs_cg.jl index e662ebbc..f376191a 100644 --- a/bench/bench_tron_minres_vs_cg.jl +++ b/bench/bench_tron_minres_vs_cg.jl @@ -1,13 +1,6 @@ -# Benchmark TRON subsolver: CG vs MINRES -# Usage: julia --project=. bench/bench_tron_minres_vs_cg.jl - using Pkg, LinearAlgebra, Printf, Dates -Pkg.activate(joinpath(@__DIR__, "..")) - - using JSOSolvers, Random, ADNLPModels -# Build a logistic-like smooth convex objective as an ADNLPModel function make_logistic(n, m; rng = Random.GLOBAL_RNG) A = randn(rng, m, n) xtrue = randn(rng, n) @@ -47,7 +40,6 @@ function main() end end - # Print summary println("prob,subsolver,elapsed,iter,status,objective,dual") for r in results @printf("%d,%s,%.6f,%d,%s,%.12g,%.6g\n", r[:prob], string(r[:subsolver]), r[:elapsed], r[:iter], string(r[:status]), r[:objective], r[:dual]) diff --git a/src/tron.jl b/src/tron.jl index da4ee9ef..bd8c892f 100644 --- a/src/tron.jl +++ b/src/tron.jl @@ -231,11 +231,9 @@ end end # Construct solver with chosen subsolver. We intentionally do not forward - # additional TRON trust-region kwargs here to keep the call simple; the - # TronSolver constructor will use defaults for those parameters. + # additional TRON trust-region kwargs here to keep the call simple. solver = TronSolver(nlp; μ₀ = μ₀, μ₁ = μ₁, σ = σ, subsolver = subsolver) - # Remaining dict contains kwargs meant for solve! (we removed :subsolver) return solve!(solver, nlp; x = x, dict...) end @@ -596,23 +594,15 @@ Small compatibility wrapper around Krylov.jl's `krylov_solve!` to handle workspaces that don't accept the `radius` keyword (e.g. Minres). We first attempt to call with `radius`, and on a MethodError we retry without it. """ -function run_krylov_subsolver!(workspace, A, b; radius=nothing, rtol=1e-6, atol=0.0, timemax=Inf, verbose=0) - # Try calling with `radius` first (covers CG-like workspaces). +function run_krylov_subsolver!(workspace, A, b; radius=nothing, rtol, atol, timemax, verbose) if radius !== nothing try return krylov_solve!(workspace, A, b; radius = radius, rtol = rtol, atol = atol, timemax = timemax, verbose = verbose) catch e - # If the failure is due to an unsupported keyword (MethodError), fall - # back to calling with a callback that enforces the trust-region by - # inspecting the Krylov workspace's current solution `x` at each - # iteration. Any other error is rethrown. if isa(e, MethodError) - # Krylov callbacks receive the workspace and should return `true` - # to request termination. We build a closure that captures `radius`. function tr_callback(cb_workspace) - # Some workspaces store the current iterate in `x`. if !hasfield(typeof(cb_workspace), :x) - error("Krylov workspace of type $(typeof(cb_workspace)) does not have an `x` field. Trust-region constraint cannot be enforced. Please use a compatible workspace type.") + return false end xcur = getfield(cb_workspace, :x) # Stop when the squared norm of the current solution reaches the squared radius (avoids sqrt).