Add parametric sensitivity calculations from derivatives wrt parameters by jsphchoi · Pull Request #647 · madsuite-org/MadNLP.jl

jsphchoi · 2026-09-21T22:26:02Z

Adds 2 new functions to export/public API:

sensitivity(solver, p)

where solver = MadNLPSolver(...) after solve!ing and p is from add_par with sensitivity = true.
calculates ds*/dp where s* = (x*, y*, zl*, zu*), returned as (; dx, dy, dzl, dzu)

sensitivity_result(p , p_new)

evaluates w*(p_new) by taking the step from p -> p_new, so basically just s* + ds*/dp*(p_new - p), returns a MadNLPExecutionStats struct w/ the model evaluated at the new solution

the actual derivative calculations dc/dp and d^2L/dxdp happens on the ExaModels side, since the AD happens there (also so support is general for other NLPModels solvers)

minimal example of what it looks like in practice, tested and confirmed on cpu and gpu, combined with the corresponding ExaModels PR "Add differentiate wrt parameters for parametric sensitivity calculations":

# Toy example from Pirnay 2012 sIPOPT paper
function pirnay(θ0; backend = nothing)
    c = ExaCore(; backend)
    @add_par(c, θ, θ0; sensitivity = true) # <--- New on ExaModels.jl side
    @add_var(c, x, 3; lvar = 0.0, start = 0.5)
    @add_obj(c, x[i]^2 for i = 1:3)
    @add_con(c, 6x[1] + 3x[2] + 2x[3] - θ[1])
    @add_con(c, θ[2] * x[1] + x[2] - x[3] - 1)
    return ExaModel(c), θ
end

# Solve w/ CPU
m, θ = pirnay([4.5, 1.0])
solver = MadNLPSolver(m; tol = 1e-8)
solve!(solver)
sensitivity(solver, θ)  # <--- New on MadNLP.jl side, returns (; dx/dθ, dy/dθ, dzl/dθ, dzu/dθ)
sensitivity_result(solver, θ, [4.0, 1.0]).solution  # <--- New on MadNLP.jl side, returns soln at new theta

# Solve w/ GPU
m, θ = pirnay([4.5, 1.0]; backend = CUDABackend())
solver = MadNLPSolver(m; tol = 1e-8, linear_solver = MadNLPGPU.CUDSSSolver)
solve!(solver)
sensitivity(solver, θ)  # <--- New on MadNLP.jl side, returns (; dx/dθ, dy/dθ, dzl/dθ, dzu/dθ)
sensitivity_result(solver, θ, [4.0, 1.0]).solution  # <--- New on MadNLP.jl side, returns soln at new theta

klamike · 2026-09-21T23:45:37Z

Super happy to see there is some interest in this!! You may be interested to check out (the now rather stale, unfortunately) work in https://github.com/madsuite-org/MadDiff.jl (and the forks linked in its readme). It would be great to finally land those features properly!

frapac · 2026-09-22T15:04:37Z

Very nice work! I am glad we are pursuing this direction in MadNLP.

We had a brief discussion together with @klamike . I think it would be good to provide also the reverse pass (vjp) for the user, that can be useful for some applications (e.g. machine learning).

The only part I am a bit unsure with is the user interface. We should think a little bit about what does the user expect exactly. I don't have any better idea for the moment, but I think it's worth discussing.

frapac

Overall looks good to me, I have only a few minor comments. I think this code should belongs to MadNLP, and that we should make sensitivity analysis a first-class citizen in our package suite.

I would suggest making explicit that we are computing the sensitivities using forward-mode AD, e.g. by adding _jvp in the name (I welcome any better suggestion of course). The reverse mode is slightly different, as the unreduced KKT system is not symmetric (hence the order in which we reduce the RHS matters). MadDiff figures this out.

frapac · 2026-09-22T19:00:17Z

+- `px`: variable block, `nvar × k`
+- `py`: constraint block, `ncon × k`
+"""
+function backsolve_kkt!(solver::AbstractMadNLPSolver{T}, px::AbstractMatrix, py::AbstractMatrix) where T


genuine question: why keeping px and py split?

frapac · 2026-09-22T19:01:04Z

+- `py`: constraint block, `ncon × k`
+"""
+function backsolve_kkt!(solver::AbstractMadNLPSolver{T}, px::AbstractMatrix, py::AbstractMatrix) where T
+    get_status(solver) in (SOLVE_SUCCEEDED, SOLVED_TO_ACCEPTABLE_LEVEL) ||


I would suggest reformulating as:

if get_status(solver) in (SOLVE_SUCCEEDED, SOLVED_TO_ACCEPTABLE_LEVEL) ... end

(I think it is more readable)

frapac · 2026-09-22T19:02:22Z

+    ifree, ufree = _free_indices(cb)
+    x0 = get_x0(nlp)
+    dx, dy, dzl, dzu = (fill!(similar(x0, n, k), zero(T)) for n in (nvar, ncon, nvar, nvar))
+    for j in 1:k


to decide together with @sshin23 : should we keep the for loop inside the function, or outside? If the user wants to assemble the full Jacobian, we have to make explicit that this is an expensive operation in the large-scale regime (we store a dense matrix with size (n+m, k))

frapac · 2026-09-22T19:02:51Z

+            unpack_z!(view(dz, :, j), cb, view(zbuf, 1:nx))
+        end
+    end
+    return (; dx, dy, dzl, dzu)


suggestion: maybe a named-tuple?

frapac · 2026-09-22T19:04:01Z

+"""
+sensitivity(solver::AbstractMadNLPSolver, Hxθ::AbstractMatrix, Jθ::AbstractMatrix) =
+    backsolve_kkt!(solver, .-Hxθ, .-Jθ)
+function sensitivity(solver::AbstractMadNLPSolver, θ...)


I would prefer keeping the argument θ explicit instead of using the splatting θ.... What is the type of θ exactly?

frapac · 2026-09-22T19:06:28Z

+- `restore_parameter`: restore `nlp[θ]` after the call
+- `recompute_residuals`: recompute `primal_feas` and `dual_feas` at the new solution
+"""
+function sensitivity_result(solver::AbstractMadNLPSolver, θ, θnew; restore_parameter = true, recompute_residuals = false)


I am not sure that this particular operation should be implemented. I would let the user calls sensitivitity directly. Let me know what do you think.

frapac · 2026-09-22T19:06:46Z

+    θ0 = copy(nlp[θ])
+    dθ = copyto!(similar(θ0), θnew) .- θ0
+    s = sensitivity(solver, θ)
+    nlp[θ] = θnew


This looks unfamiliar to me. Is it also particular to ExaModels?

frapac · 2026-09-22T19:07:27Z

+    stats.multipliers_L .+= s.dzl * dθ
+    stats.multipliers_U .+= s.dzu * dθ
+    stats.objective = NLPModels.obj(nlp, stats.solution)
+    get_ncon(nlp) > 0 && NLPModels.cons!(nlp, stats.solution, stats.constraints)


for the two following lines, I would write explicit if statement

frapac · 2026-09-22T19:07:48Z

+end
+
+_violation(v, l, u) = max(maximum(l .- v; init = zero(eltype(v))), maximum(v .- u; init = zero(eltype(v))))
+function _residuals!(stats::MadNLPExecutionStats, nlp)


Maybe add a comment to explain what this function is doing?

frapac · 2026-09-22T19:10:15Z

+(nlp::ParametricModel)(x, y) = (nlp.Hxθ, nlp.Jθ)
+(nlp::ParametricModel)(::Θ, x, y) = nlp(x, y)
+Base.getindex(nlp::ParametricModel, ::Θ) = nlp.θ
+Base.setindex!(nlp::ParametricModel, v, ::Θ) = (nlp.θ .= v; nlp)


This is exactly the kind of overloading operations that is not standard in NLPModels. This can be a solution in the short-term, but I think it's worth discussing a long-term solution that can suit anybody. Maybe time to revive the discussion in JuliaSmoothOptimizers/NLPModels.jl#557 ?

jsphchoi · 2026-09-22T23:18:37Z

Super happy to see there is some interest in this!! You may be interested to check out (the now rather stale, unfortunately) work in https://github.com/madsuite-org/MadDiff.jl (and the forks linked in its readme). It would be great to finally land those features properly!

This is exactly the kind of overloading operations that is not standard in NLPModels. This can be a solution in the short-term, but I think it's worth discussing a long-term solution that can suit anybody. Maybe time to revive the discussion in JuliaSmoothOptimizers/NLPModels.jl#557 ?

I spoke w/ sungho earlier and we agree it may be a good time to bring back ParametricNLPModels.jl, i moved over the relevant api from these PRs and also merged some of ur ideas from ur older version to a madsuite-org version @klamike

jsphchoi added 2 commits September 21, 2026 14:46

add parametric sensitivity calc from solver

ac0fd27

edit test tol

96742fe

sshin23 requested a review from frapac September 22, 2026 13:57

frapac reviewed Sep 22, 2026

View reviewed changes

PR v2 w/ ParametricNLPModels.jl

2578fd8

Conversation

jsphchoi commented Sep 21, 2026

Uh oh!

klamike commented Sep 21, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

frapac commented Sep 22, 2026

Uh oh!

frapac left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jsphchoi commented Sep 22, 2026

Uh oh!

This branch has not been deployed

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

klamike commented Sep 21, 2026 •

edited

Loading