diff --git a/src/NLPModelsModifiers.jl b/src/NLPModelsModifiers.jl index 2a81e16..a713f7c 100644 --- a/src/NLPModelsModifiers.jl +++ b/src/NLPModelsModifiers.jl @@ -18,6 +18,7 @@ include("feasibility-form-nls.jl") include("feasibility-residual.jl") include("quasi-newton.jl") include("slack-model.jl") +include("scaled-model.jl") include("model-interaction.jl") end # module diff --git a/src/scaled-model.jl b/src/scaled-model.jl new file mode 100644 index 0000000..8e794ec --- /dev/null +++ b/src/scaled-model.jl @@ -0,0 +1,335 @@ +export ScaledModel + +struct ConservativeScaling{T} + max_gradient::T +end + +function _set_constraints_scaling!(scaling, Ji, Jj, Jx, max_gradient) + # Store norm(∇cᵢ, Inf) at index i of vector scaling + for (i, j, x) in zip(Ji, Jj, Jx) + scaling[i] = max(scaling[i], abs(x)) + end + # Compute scaling as min(1, max_gradient / norm(∇cᵢ, Inf) ) + for i in eachindex(scaling) + scaling[i] = min(1.0, max_gradient / scaling[i]) + end +end + +function _set_jacobian_scaling!(Jx, Ji, Jj, scaling) + k = 0 + for (i, j) in zip(Ji, Jj) + Jx[k += 1] = scaling[i] + end +end + +function scale_model!(scaling::ConservativeScaling{T}, nlp) where T + n, m = NLPModels.get_nvar(nlp), NLPModels.get_ncon(nlp) + nnzj = NLPModels.get_nnzj(nlp) + x0 = NLPModels.get_x0(nlp) + g = NLPModels.grad(nlp, x0) + scaling_obj = min(one(T), scaling.max_gradient / norm(g, Inf)) + scaling_cons = similar(x0, m) + scaling_jac = similar(x0, nnzj) + fill!(scaling_cons, zero(T)) + Ji, Jj = NLPModels.jac_structure(nlp) + NLPModels.jac_coord!(nlp, x0, scaling_jac) + _set_constraints_scaling!(scaling_cons, Ji, Jj, scaling_jac, scaling.max_gradient) + _set_jacobian_scaling!(scaling_jac, Ji, Jj, scaling_cons) + return (scaling_obj, scaling_cons, scaling_jac) +end + +@doc raw""" + ScaledModel + +Scale the nonlinear program +```math +\begin{aligned} + min_x \quad & f(x)\\ +\mathrm{s.t.} \quad &  c_L ≤ c(x) ≤ c_U,\\ + & ℓ ≤ x ≥ u, +\end{aligned} +``` +as +```math +\begin{aligned} + min_x \quad & σf . f(x)\\ +\mathrm{s.t.} \quad &  σc . c_L ≤ σc . c(x) ≤ σc . c_U, \\ + & ℓ ≤ x ≥ u, +\end{aligned} +``` +with ``σf`` a positive scalar defined as +``` +σf = min(1, max_gradient / norm(g0, Inf)) +``` +and ``σc`` a vector whose size is equal to the number of constraints in the model. +For ``i=1, ..., m``, +``` +σc[i] = min(1, max_gradient / norm(J0[i, :], Inf)) + +``` + +The vector ``g0 = ∇f(x0)`` and the matrix ``J0 = ∇c(x0)`` are resp. +the gradient and the Jacobian evaluated at the initial point ``x0``. +By default, the threshold parameter `max_gradient` is set to 100.0. + +The method has been originally proposed in Ipopt [1]. + +## Reference + +[1] Wächter, A., & Biegler, L. T. (2006). +On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. +Mathematical programming, 106(1), 25-57. + +""" +struct ScaledModel{T, S, M} <: NLPModels.AbstractNLPModel{T, S} + nlp::M + meta::NLPModels.NLPModelMeta{T, S} + counters::NLPModels.Counters + scaling_obj::T + scaling_cons::S # [size m] + scaling_cons_lin::S # [size nlin] + scaling_cons_nln::S # [size nnln] + scaling_jac::S # [size nnzj] + scaling_jac_lin::S # [size lin_nnzj] + scaling_jac_nln::S # [size nln_nnzj] + buffer_cons::S # [size m] +end + +function ScaledModel( + nlp::NLPModels.AbstractNLPModel{T, S}; + scaling=ConservativeScaling(T(100)), +) where {T, S} + n, m = NLPModels.get_nvar(nlp), NLPModels.get_ncon(nlp) + x0 = NLPModels.get_x0(nlp) + buffer_cons = S(undef, m) + + # Compute scaling for the problem as a whole. + scaling_obj, scaling_cons, scaling_jac = scale_model!(scaling, nlp) + + # Get scaling for linear and nonlinear constraints. + scaling_cons_lin = scaling_cons[nlp.meta.lin] + scaling_cons_nln = scaling_cons[nlp.meta.nln] + scaling_jac_lin = zeros(T, nlp.meta.lin_nnzj) + Jlin_i, Jlin_j = NLPModels.jac_lin_structure(nlp) + k = 0 + for (i, j) in zip(Jlin_i, Jlin_j) + scaling_jac_lin[k += 1] = scaling_cons_lin[i] + end + scaling_jac_nln = zeros(T, nlp.meta.nln_nnzj) + Jnln_i, Jnln_j = NLPModels.jac_nln_structure(nlp) + k = 0 + for (i, j) in zip(Jnln_i, Jnln_j) + scaling_jac_nln[k += 1] = scaling_cons_nln[i] + end + + # Copy metadata from original problem, with some modifications. + meta = NLPModels.NLPModelMeta( + nlp.meta; + y0 = NLPModels.get_y0(nlp) .* scaling_cons, + lcon = NLPModels.get_lcon(nlp) .* scaling_cons, + ucon = NLPModels.get_ucon(nlp) .* scaling_cons, + name="scaled-" * nlp.meta.name, + ) + + return ScaledModel( + nlp, + meta, + NLPModels.Counters(), + scaling_obj, + scaling_cons, + scaling_cons_lin, + scaling_cons_nln, + scaling_jac, + scaling_jac_lin, + scaling_jac_nln, + buffer_cons, + ) +end + +function NLPModels.obj(nlp::ScaledModel{T, S}, x::AbstractVector) where {T, S <: AbstractVector{T}} + @lencheck nlp.meta.nvar x + return nlp.scaling_obj * NLPModels.obj(nlp.nlp, x) +end + +function NLPModels.grad!(nlp::ScaledModel, x::AbstractVector, g::AbstractVector) + @lencheck nlp.meta.nvar x g + NLPModels.grad!(nlp.nlp, x, g) + g .*= nlp.scaling_obj + return g +end + +function NLPModels.cons!(nlp::ScaledModel, x::AbstractVector, c::AbstractVector) + @lencheck nlp.meta.nvar x + @lencheck nlp.meta.ncon c + NLPModels.cons!(nlp.nlp, x, c) + c .*= nlp.scaling_cons + return c +end + +function NLPModels.cons_lin!(nlp::ScaledModel, x::AbstractVector, c::AbstractVector) + @lencheck nlp.meta.nvar x + @lencheck nlp.meta.nlin c + NLPModels.cons_lin!(nlp.nlp, x, c) + c .*= nlp.scaling_cons_lin + return c +end + +function NLPModels.cons_nln!(nlp::ScaledModel, x::AbstractVector, c::AbstractVector) + @lencheck nlp.meta.nvar x + @lencheck nlp.meta.nnln c + NLPModels.cons_nln!(nlp.nlp, x, c) + c .*= nlp.scaling_cons_nln + return c +end + +function NLPModels.jprod!(nlp::ScaledModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector) + @lencheck nlp.meta.nvar x v + @lencheck nlp.meta.ncon Jv + NLPModels.jprod!(nlp.nlp, x, v, Jv) + Jv .*= nlp.scaling_cons + return Jv +end + +function NLPModels.jprod_lin!(nlp::ScaledModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector) + @lencheck nlp.meta.nvar x v + @lencheck nlp.meta.nlin Jv + NLPModels.jprod_lin!(nlp.nlp, x, v, Jv) + Jv .*= nlp.scaling_cons_lin + return Jv +end + +function NLPModels.jprod_nln!(nlp::ScaledModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector) + @lencheck nlp.meta.nvar x v + @lencheck nlp.meta.nnln Jv + NLPModels.jprod_lin!(nlp.nlp, x, v, Jv) + Jv .*= nlp.scaling_cons_nln + return Jv +end + +function NLPModels.jtprod!(nlp::ScaledModel, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector) + @lencheck nlp.meta.nvar x Jtv + @lencheck nlp.meta.ncon v + v_scaled = nlp.buffer_cons + v_scaled .= v .* nlp.scaling_cons + NLPModels.jtprod!(nlp.nlp, x, v_scaled, Jtv) + return Jtv +end + +function NLPModels.jtprod_lin!(nlp::ScaledModel, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector) + @lencheck nlp.meta.nvar x Jtv + @lencheck nlp.meta.nlin v + v_scaled = view(nlp.buffer_cons, 1:nlp.meta.nlin) + v_scaled .= v .* nlp.scaling_cons_lin + NLPModels.jtprod_lin!(nlp.nlp, x, v_scaled, Jtv) + return Jtv +end + +function NLPModels.jtprod_nln!(nlp::ScaledModel, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector) + @lencheck nlp.meta.nvar x Jtv + @lencheck nlp.meta.nnln v + v_scaled = view(nlp.buffer_cons, 1:nlp.meta.nnln) + v_scaled .= v .* nlp.scaling_cons_nln + NLPModels.jtprod_nln!(nlp.nlp, x, v_scaled, Jtv) + return Jtv +end + +function NLPModels.jac_structure!(nlp::ScaledModel, jrows::AbstractVector, jcols::AbstractVector) + @lencheck nlp.meta.nnzj jrows jcols + NLPModels.jac_structure!(nlp.nlp, jrows, jcols) + return jrows, jcols +end + +function NLPModels.jac_lin_structure!(nlp::ScaledModel, jrows::AbstractVector, jcols::AbstractVector) + NLPModels.jac_lin_structure!(nlp.nlp, jrows, jcols) + return jrows, jcols +end + +function NLPModels.jac_nln_structure!(nlp::ScaledModel, jrows::AbstractVector, jcols::AbstractVector) + NLPModels.jac_nln_structure!(nlp.nlp, jrows, jcols) + return jrows, jcols +end + +function NLPModels.jac_coord!(nlp::ScaledModel, x::AbstractVector, jac::AbstractVector) + NLPModels.jac_coord!(nlp.nlp, x, jac) + jac .*= nlp.scaling_jac + return jac +end + +function NLPModels.jac_lin_coord!(nlp::ScaledModel, x::AbstractVector, jac::AbstractVector) + NLPModels.jac_lin_coord!(nlp.nlp, x, jac) + jac .*= nlp.scaling_jac_lin + return jac +end + +function NLPModels.jac_nln_coord!(nlp::ScaledModel, x::AbstractVector, jac::AbstractVector) + NLPModels.jac_nln_coord!(nlp.nlp, x, jac) + jac .*= nlp.scaling_jac_nln + return jac +end + +function NLPModels.hess_structure!(nlp::ScaledModel, hrows::AbstractVector, hcols::AbstractVector) + @lencheck nlp.meta.nnzh hrows hcols + NLPModels.hess_structure!(nlp.nlp, hrows, hcols) + return hrows, hcols +end + +function NLPModels.hess_coord!( + nlp::ScaledModel, + x::AbstractVector, + vals::AbstractVector; + obj_weight::Real=one(eltype(x)), +) + @lencheck nlp.meta.nvar x + @lencheck nlp.meta.nnzh vals + σ = obj_weight * nlp.scaling_obj + NLPModels.hess_coord!(nlp.nlp, x, vals; obj_weight=σ) + return vals +end + +function NLPModels.hess_coord!( + nlp::ScaledModel, + x::AbstractVector, + y::AbstractVector, + vals::AbstractVector; + obj_weight::Real=one(eltype(x)), +) + @lencheck nlp.meta.nvar x + @lencheck nlp.meta.ncon y + @lencheck nlp.meta.nnzh vals + y_scaled = nlp.buffer_cons + y_scaled .= y .* nlp.scaling_cons + σ = obj_weight * nlp.scaling_obj + NLPModels.hess_coord!(nlp.nlp, x, y_scaled, vals; obj_weight=σ) + return vals +end + +function NLPModels.hprod!( + nlp::ScaledModel, + x::AbstractVector, + v::AbstractVector, + hv::AbstractVector; + obj_weight::Real = one(eltype(x)), +) + @lencheck nlp.meta.nvar x v hv + σ = obj_weight * nlp.scaling_obj + NLPModels.hprod!(nlp.nlp, x, v, hv; obj_weight = σ) + return hv +end + +function NLPModels.hprod!( + nlp::ScaledModel, + x::AbstractVector, + y::AbstractVector, + v::AbstractVector, + hv::AbstractVector; + obj_weight::Real = one(eltype(x)), +) + @lencheck nlp.meta.nvar x v hv + @lencheck nlp.meta.ncon y + y_scaled = nlp.buffer_cons + y_scaled .= y .* nlp.scaling_cons + σ = obj_weight * nlp.scaling_obj + NLPModels.hprod!(nlp.nlp, x, y, v, hv; obj_weight = σ) + return hv +end + diff --git a/test/nlp/scaled-model.jl b/test/nlp/scaled-model.jl new file mode 100644 index 0000000..418516c --- /dev/null +++ b/test/nlp/scaled-model.jl @@ -0,0 +1,40 @@ +@testset "ScaledModel NLP tests" begin + @testset "API" for T in [Float64, Float32], M in [NLPModelMeta, SimpleNLPMeta] + original_nlp = SimpleNLPModel(T, M) + nlp = ScaledModel(original_nlp) + σ_obj, σ_cons = nlp.scaling_obj, nlp.scaling_cons + + # Hand-code the scaled problem from the original NLP. + f(x) = σ_obj * NLPModels.obj(original_nlp, x) + ∇f(x) = σ_obj .* NLPModels.grad(original_nlp, x) + H(x) = σ_obj .* NLPModels.hess(original_nlp, x) + c(x) = σ_cons .* NLPModels.cons(original_nlp, x) + J(x) = Diagonal(σ_cons) * NLPModels.jac(original_nlp, x) + H(x, y) = NLPModels.hess(original_nlp, x, σ_cons .* y; obj_weight=σ_obj) + + n = nlp.meta.nvar + m = nlp.meta.ncon + @test nlp.meta.x0 == T[2; 2] + + x = randn(T, n) + y = randn(T, m) + v = randn(T, n) + w = randn(T, m) + Jv = zeros(T, m) + Jtw = zeros(T, n) + Hv = zeros(T, n) + Hvals = zeros(T, nlp.meta.nnzh) + + # Basic methods + @test obj(nlp, x) ≈ f(x) + @test grad(nlp, x) ≈ ∇f(x) + @test hess(nlp, x) ≈ H(x) + @test hprod(nlp, x, v) ≈ H(x) * v + @test cons(nlp, x) ≈ c(x) + @test jac(nlp, x) ≈ J(x) + @test jprod(nlp, x, v) ≈ J(x) * v + @test jtprod(nlp, x, w) ≈ J(x)' * w + @test hess(nlp, x, y) ≈ H(x, y) + @test hprod(nlp, x, y, v) ≈ H(x, y) * v + end +end diff --git a/test/nlp/simple-model.jl b/test/nlp/simple-model.jl index 56ab503..ea0caa5 100644 --- a/test/nlp/simple-model.jl +++ b/test/nlp/simple-model.jl @@ -60,6 +60,18 @@ mutable struct SimpleNLPMeta{T, S} <: AbstractNLPModelMeta{T, S} minimize::Bool islp::Bool name::String + variable_bounds_analysis::Bool + constraint_bounds_analysis::Bool + + sparse_jacobian::Bool + sparse_hessian::Bool + + grad_available::Bool + jac_available::Bool + hess_available::Bool + jprod_available::Bool + jtprod_available::Bool + hprod_available::Bool function SimpleNLPMeta{T, S}( nvar::Int; x0::S = fill!(S(undef, nvar), zero(T)), @@ -81,6 +93,16 @@ mutable struct SimpleNLPMeta{T, S} <: AbstractNLPModelMeta{T, S} minimize = true, islp = false, name = "Generic", + variable_bounds_analysis::Bool = true, + constraint_bounds_analysis::Bool = true, + sparse_jacobian::Bool = true, + sparse_hessian::Bool = true, + grad_available::Bool = true, + jac_available::Bool = (ncon > 0), + hess_available::Bool = true, + jprod_available::Bool = (ncon > 0), + jtprod_available::Bool = (ncon > 0), + hprod_available::Bool = true, ) where {T, S} if (nvar < 1) || (ncon < 0) error("Nonsensical dimensions") @@ -147,6 +169,16 @@ mutable struct SimpleNLPMeta{T, S} <: AbstractNLPModelMeta{T, S} minimize, islp, name, + variable_bounds_analysis, + constraint_bounds_analysis, + sparse_jacobian, + sparse_hessian, + grad_available, + jac_available, + hess_available, + jprod_available, + jtprod_available, + hprod_available, ) end end diff --git a/test/runtests.jl b/test/runtests.jl index 4696795..e0fd9b9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -4,6 +4,7 @@ using LinearOperators, NLPModels, NLPModelsModifiers include("nlp/simple-model.jl") include("nlp/quasi-newton.jl") include("nlp/slack-model.jl") +include("nlp/scaled-model.jl") include("nls/simple-model.jl") include("nls/feasibility-form-nls.jl")