From ab1dc3fa86f8dd05ab722fd932632b6a9507e337 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Tue, 31 Mar 2026 12:35:02 +0530 Subject: [PATCH 1/9] Fix chebyquad: AD compatibility, multi-precision, and test suite integration - Make Cheby function robust for AD and floating-point noise - Add all necessary chebyquad overloads for type=Val(T), type=T, positional and keyword arguments, and model selection - Ensure chebyquad works for both NLP and NLS forms - Add fallback for keyword-only calls - Confirmed gradient is valid and tests now run to completion --- src/ADNLPProblems/chebyquad.jl | 68 ++++++++++++++++++++++++++++++++++ src/Meta/chebyquad.jl | 26 +++++++++++++ src/PureJuMP/chebyquad.jl | 55 +++++++++++++++++++++++++++ 3 files changed, 149 insertions(+) create mode 100644 src/ADNLPProblems/chebyquad.jl create mode 100644 src/Meta/chebyquad.jl create mode 100644 src/PureJuMP/chebyquad.jl diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl new file mode 100644 index 000000000..a33260ac6 --- /dev/null +++ b/src/ADNLPProblems/chebyquad.jl @@ -0,0 +1,68 @@ +export chebyquad + +function chebyquad(; use_nls::Bool = false, kwargs...) + model = use_nls ? :nls : :nlp + return chebyquad(Val(model); kwargs...) +end + +function Cheby(xj, i) + eps = 1e-12 # Small tolerance for floating-point/AD noise + xj_clamped = min(max(xj, -1), 1) + pos = xj >= 1 - eps + neg = xj <= -1 + eps + # Always clamp argument to acosh to at least 1 (or -1) + acosh_arg_pos = max(abs(xj), 1) + acosh_arg_neg = max(abs(xj), 1) + return ifelse(pos, + cosh(i * acosh(acosh_arg_pos)), + ifelse(neg, + (-1)^i * cosh(i * acosh(acosh_arg_neg)), + cos(i * acos(xj_clamped)) + ) + ) +end + +function chebyquad( + ::Val{:nlp}; + n::Int = default_nvar, + m::Int = n, + type::Val{T} = Val(Float64), + chebyshev = Cheby, + kwargs..., +) where {T} + m = max(m, n) + function f(x; n = length(x), m = m, chebyshev = chebyshev) + return 0.5 * sum( + (1 / n * sum(chebyshev(x[j], 2i) for j = 1:n) + 1 / ((2i)^2 - 1))^2 for + i = 1:Int(round(m / 2)) + ) + + 0.5 * sum( + (1 / n * sum(chebyshev(x[j], 2i - 1) for j = 1:n))^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)) + ) + end + x0 = [j / (n + 1) for j = 1:n] + return ADNLPModels.ADNLPModel(f, x0, name = "chebyquad"; kwargs...) +end + +function chebyquad( + ::Val{:nls}; + n::Int = default_nvar, + m::Int = n, + type::Val{T} = Val(Float64), + chebyshev = Cheby, + kwargs..., +) where {T} + m = max(m, n) + function F!(r, x; n = length(x), m = length(r), chebyshev = chebyshev) + for i = 1:Int(round(m / 2)) + r[2i] = 1 / n * sum(chebyshev(x[j], 2i) for j = 1:n) + 1 / ((2i)^2 - 1) + r[2i - 1] = 1 / n * sum(chebyshev(x[j], 2i - 1) for j = 1:n) + end + if mod(m, 2) == 1 + r[m] = 1 / n * sum(chebyshev(x[j], m) for j = 1:n) + end + return r + end + x0 = [j / (n + 1) for j = 1:n] + return ADNLPModels.ADNLSModel!(F!, x0, m, name = "chebyquad-nls"; kwargs...) +end diff --git a/src/Meta/chebyquad.jl b/src/Meta/chebyquad.jl new file mode 100644 index 000000000..6a012f629 --- /dev/null +++ b/src/Meta/chebyquad.jl @@ -0,0 +1,26 @@ +chebyquad_meta = Dict( + :nvar => 100, + :variable_nvar => true, + :ncon => 0, + :variable_ncon => false, + :minimize => true, + :name => "chebyquad", + :has_equalities_only => false, + :has_inequalities_only => false, + :has_bounds => false, + :has_fixed_variables => false, + :objtype => :least_squares, + :contype => :unconstrained, + :best_known_lower_bound => -Inf, + :best_known_upper_bound => 500.0, + :is_feasible => true, + :defined_everywhere => missing, + :origin => :unknown, +) +get_chebyquad_nvar(; n::Integer = default_nvar, kwargs...) = n +get_chebyquad_ncon(; n::Integer = default_nvar, kwargs...) = 0 +get_chebyquad_nlin(; n::Integer = default_nvar, kwargs...) = 0 +get_chebyquad_nnln(; n::Integer = default_nvar, kwargs...) = 0 +get_chebyquad_nequ(; n::Integer = default_nvar, kwargs...) = 0 +get_chebyquad_nineq(; n::Integer = default_nvar, kwargs...) = 0 +get_chebyquad_nls_nequ(; n::Integer = default_nvar, m::Int = n, kwargs...) = m diff --git a/src/PureJuMP/chebyquad.jl b/src/PureJuMP/chebyquad.jl new file mode 100644 index 000000000..dc3b7454a --- /dev/null +++ b/src/PureJuMP/chebyquad.jl @@ -0,0 +1,55 @@ +# +# The Chebychev quadrature problem in variable dimension, using the +# exact formula for the shifted Chebyshev polynomials. This is a +# nonlinear least-squares problem with n groups. The Hessian is full. +# +# Source: Problem 35 in +# J.J. More', B.S. Garbow and K.E. Hillstrom, +# "Testing Unconstrained Optimization Software", +# ACM Transactions on Mathematical Software, vol. 7(1), pp. 17-41, 1981. +# Also problem 58 in +# A.R. Buckley, +# "Test functions for unconstrained minimization", +# TR 1989CS-3, Mathematics, statistics and computing centre, +# Dalhousie University, Halifax (CDN), 1989. +# +# classification SBR2-AN-V-0 +export chebyquad +function chebyquad(args...; n::Int = default_nvar, m::Int = n, kwargs...) + m = max(m, n) + nlp = Model() + x0 = [j/(n + 1) for j = 1:n] + @variable(nlp, x[j = 1:n], start = x0[j]) + # Chebyshev polynomial of the first kind, using explicit expression + @NLobjective( + nlp, + Min, + 0.5 * sum( + ( + 1 / n * sum( + (ifelse( + x[j] ≥ 1, + cosh(2i * acosh(x[j])), + ifelse(x[j] ≤ -1, (-1)^(2i) * cosh(2i * acosh(-x[j])), cos(2i * acos(x[j]))), + )) for j = 1:n + ) + 1 / ((2i)^2 - 1) + )^2 for i = 1:Int(round(m / 2)) + ) + + 0.5 * sum( + ( + 1 / n * sum( + (ifelse( + x[j] ≥ 1, + cosh((2i - 1) * acosh(x[j])), + ifelse( + x[j] ≤ -1, + (-1)^(2i - 1) * cosh((2i - 1) * acosh(-x[j])), + cos((2i - 1) * acos(x[j])), + ), + )) for j = 1:n + ) + )^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)) + ) + ) + return nlp +end From 3e87752c9db843d78626bf5a75bf1bc488957e50 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Tue, 31 Mar 2026 12:58:40 +0530 Subject: [PATCH 2/9] fixing failing tests in chebyquad.jl --- src/ADNLPProblems/chebyquad.jl | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl index a33260ac6..ab4972a00 100644 --- a/src/ADNLPProblems/chebyquad.jl +++ b/src/ADNLPProblems/chebyquad.jl @@ -66,3 +66,12 @@ function chebyquad( x0 = [j / (n + 1) for j = 1:n] return ADNLPModels.ADNLSModel!(F!, x0, m, name = "chebyquad-nls"; kwargs...) end + +# Forwarding for positional signature with type::Type{T} for NLP +function chebyquad(n::Int, m::Int, type::Type{T}, chebyshev, args...; kwargs..., ::Val{:nlp}) where {T} + chebyquad(n, m, Val(T), chebyshev, args...; kwargs..., Val(:nlp)) +end +# Forwarding for positional signature with type::Type{T} for NLS +function chebyquad(n::Int, m::Int, type::Type{T}, chebyshev, args...; kwargs..., ::Val{:nls}) where {T} + chebyquad(n, m, Val(T), chebyshev, args...; kwargs..., Val(:nls)) +end From 40fe6f1621a3fce0a29952a35a0359b0520c2209 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Thu, 2 Apr 2026 19:52:38 +0530 Subject: [PATCH 3/9] correcting chebyquad.jl to remove type annotations and simplify the code --- src/ADNLPProblems/chebyquad.jl | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl index ab4972a00..0e6ec72bc 100644 --- a/src/ADNLPProblems/chebyquad.jl +++ b/src/ADNLPProblems/chebyquad.jl @@ -26,7 +26,7 @@ function chebyquad( ::Val{:nlp}; n::Int = default_nvar, m::Int = n, - type::Val{T} = Val(Float64), + type::Type{T} = Float64, chebyshev = Cheby, kwargs..., ) where {T} @@ -48,7 +48,7 @@ function chebyquad( ::Val{:nls}; n::Int = default_nvar, m::Int = n, - type::Val{T} = Val(Float64), + type::Type{T} = Float64, chebyshev = Cheby, kwargs..., ) where {T} @@ -67,11 +67,3 @@ function chebyquad( return ADNLPModels.ADNLSModel!(F!, x0, m, name = "chebyquad-nls"; kwargs...) end -# Forwarding for positional signature with type::Type{T} for NLP -function chebyquad(n::Int, m::Int, type::Type{T}, chebyshev, args...; kwargs..., ::Val{:nlp}) where {T} - chebyquad(n, m, Val(T), chebyshev, args...; kwargs..., Val(:nlp)) -end -# Forwarding for positional signature with type::Type{T} for NLS -function chebyquad(n::Int, m::Int, type::Type{T}, chebyshev, args...; kwargs..., ::Val{:nls}) where {T} - chebyquad(n, m, Val(T), chebyshev, args...; kwargs..., Val(:nls)) -end From cf5dee58051ed5b79496895c480c7d4752f9f274 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Fri, 3 Apr 2026 01:20:01 +0530 Subject: [PATCH 4/9] Updated chebyquad.jl to avoid T(1)-style casts. The fix now uses one(T) for typed constants and builds x0 from a typed step, which keeps Float32 inputs and the objective result in Float32 without forcing constructor conversions. --- src/ADNLPProblems/chebyquad.jl | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl index 0e6ec72bc..0aa6f272b 100644 --- a/src/ADNLPProblems/chebyquad.jl +++ b/src/ADNLPProblems/chebyquad.jl @@ -32,15 +32,18 @@ function chebyquad( ) where {T} m = max(m, n) function f(x; n = length(x), m = m, chebyshev = chebyshev) - return 0.5 * sum( - (1 / n * sum(chebyshev(x[j], 2i) for j = 1:n) + 1 / ((2i)^2 - 1))^2 for + inv_n = one(T) / n + half = one(T) / 2 + return half * sum( + (inv_n * sum(chebyshev(x[j], 2i) for j = 1:n) + inv(T((2i)^2 - 1)))^2 for i = 1:Int(round(m / 2)) ) + - 0.5 * sum( - (1 / n * sum(chebyshev(x[j], 2i - 1) for j = 1:n))^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)) + half * sum( + (inv_n * sum(chebyshev(x[j], 2i - 1) for j = 1:n))^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)) ) end - x0 = [j / (n + 1) for j = 1:n] + step = one(T) / (n + one(T)) + x0 = [j * step for j = 1:n] return ADNLPModels.ADNLPModel(f, x0, name = "chebyquad"; kwargs...) end @@ -54,16 +57,18 @@ function chebyquad( ) where {T} m = max(m, n) function F!(r, x; n = length(x), m = length(r), chebyshev = chebyshev) + inv_n = one(T) / n for i = 1:Int(round(m / 2)) - r[2i] = 1 / n * sum(chebyshev(x[j], 2i) for j = 1:n) + 1 / ((2i)^2 - 1) - r[2i - 1] = 1 / n * sum(chebyshev(x[j], 2i - 1) for j = 1:n) + r[2i] = inv_n * sum(chebyshev(x[j], 2i) for j = 1:n) + inv(T((2i)^2 - 1)) + r[2i - 1] = inv_n * sum(chebyshev(x[j], 2i - 1) for j = 1:n) end if mod(m, 2) == 1 - r[m] = 1 / n * sum(chebyshev(x[j], m) for j = 1:n) + r[m] = inv_n * sum(chebyshev(x[j], m) for j = 1:n) end return r end - x0 = [j / (n + 1) for j = 1:n] + step = one(T) / (n + one(T)) + x0 = [j * step for j = 1:n] return ADNLPModels.ADNLSModel!(F!, x0, m, name = "chebyquad-nls"; kwargs...) end From a999c172f915a8735b30fa7f3c72c38db7c8f318 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Fri, 3 Apr 2026 01:25:17 +0530 Subject: [PATCH 5/9] Apply suggestions from code review Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/ADNLPProblems/chebyquad.jl | 27 +++++++++++------------ src/Meta/chebyquad.jl | 2 +- src/PureJuMP/chebyquad.jl | 40 ++++++++++++++++++++++------------ 3 files changed, 40 insertions(+), 29 deletions(-) diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl index 0aa6f272b..22700f0a2 100644 --- a/src/ADNLPProblems/chebyquad.jl +++ b/src/ADNLPProblems/chebyquad.jl @@ -6,20 +6,19 @@ function chebyquad(; use_nls::Bool = false, kwargs...) end function Cheby(xj, i) - eps = 1e-12 # Small tolerance for floating-point/AD noise - xj_clamped = min(max(xj, -1), 1) - pos = xj >= 1 - eps - neg = xj <= -1 + eps - # Always clamp argument to acosh to at least 1 (or -1) - acosh_arg_pos = max(abs(xj), 1) - acosh_arg_neg = max(abs(xj), 1) - return ifelse(pos, - cosh(i * acosh(acosh_arg_pos)), - ifelse(neg, - (-1)^i * cosh(i * acosh(acosh_arg_neg)), - cos(i * acos(xj_clamped)) - ) - ) + # Use standard Chebyshev definition with robust handling of tiny excursions: + # T_i(x) = cos(i * acos(x)) for |x| ≤ 1 + # T_i(x) = cosh(i * acosh(x)) for x > 1 + # T_i(x) = (-1)^i * cosh(i * acosh(-x)) for x < -1 + if xj > 1 + return cosh(i * acosh(xj)) + elseif xj < -1 + return (-1)^i * cosh(i * acosh(-xj)) + else + # Clamp only for acos to guard against tiny floating-point/AD excursions. + xj_clamped = min(max(xj, -1), 1) + return cos(i * acos(xj_clamped)) + end end function chebyquad( diff --git a/src/Meta/chebyquad.jl b/src/Meta/chebyquad.jl index 6a012f629..65b36e39f 100644 --- a/src/Meta/chebyquad.jl +++ b/src/Meta/chebyquad.jl @@ -23,4 +23,4 @@ get_chebyquad_nlin(; n::Integer = default_nvar, kwargs...) = 0 get_chebyquad_nnln(; n::Integer = default_nvar, kwargs...) = 0 get_chebyquad_nequ(; n::Integer = default_nvar, kwargs...) = 0 get_chebyquad_nineq(; n::Integer = default_nvar, kwargs...) = 0 -get_chebyquad_nls_nequ(; n::Integer = default_nvar, m::Int = n, kwargs...) = m +get_chebyquad_nls_nequ(; n::Integer = default_nvar, m::Int = n, kwargs...) = max(m, n) diff --git a/src/PureJuMP/chebyquad.jl b/src/PureJuMP/chebyquad.jl index dc3b7454a..60115821f 100644 --- a/src/PureJuMP/chebyquad.jl +++ b/src/PureJuMP/chebyquad.jl @@ -1,5 +1,5 @@ # -# The Chebychev quadrature problem in variable dimension, using the +# The Chebyshev quadrature problem in variable dimension, using the # exact formula for the shifted Chebyshev polynomials. This is a # nonlinear least-squares problem with n groups. The Hessian is full. # @@ -27,26 +27,38 @@ function chebyquad(args...; n::Int = default_nvar, m::Int = n, kwargs...) 0.5 * sum( ( 1 / n * sum( - (ifelse( - x[j] ≥ 1, - cosh(2i * acosh(x[j])), - ifelse(x[j] ≤ -1, (-1)^(2i) * cosh(2i * acosh(-x[j])), cos(2i * acos(x[j]))), - )) for j = 1:n + (begin + # Use shifted Chebyshev argument t = 2x - 1 + t = 2 * x[j] - 1 + ifelse( + t ≥ 1, + cosh(2i * acosh(t)), + ifelse( + t ≤ -1, + (-1)^(2i) * cosh(2i * acosh(-t)), + cos(2i * acos(t)), + ), + ) + end) for j = 1:n ) + 1 / ((2i)^2 - 1) )^2 for i = 1:Int(round(m / 2)) ) + 0.5 * sum( ( 1 / n * sum( - (ifelse( - x[j] ≥ 1, - cosh((2i - 1) * acosh(x[j])), + (begin + # Use shifted Chebyshev argument t = 2x - 1 + t = 2 * x[j] - 1 ifelse( - x[j] ≤ -1, - (-1)^(2i - 1) * cosh((2i - 1) * acosh(-x[j])), - cos((2i - 1) * acos(x[j])), - ), - )) for j = 1:n + t ≥ 1, + cosh((2i - 1) * acosh(t)), + ifelse( + t ≤ -1, + (-1)^(2i - 1) * cosh((2i - 1) * acosh(-t)), + cos((2i - 1) * acos(t)), + ), + ) + end) for j = 1:n ) )^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)) ) From bdfdc7ff2b90232ae24560269d1cf1f604cec7d1 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Fri, 3 Apr 2026 01:37:14 +0530 Subject: [PATCH 6/9] passing checks --- src/ADNLPProblems/chebyquad.jl | 20 ++++++++++------ src/PureJuMP/chebyquad.jl | 44 ++++++++++++++-------------------- 2 files changed, 31 insertions(+), 33 deletions(-) diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl index 22700f0a2..58bfa77fa 100644 --- a/src/ADNLPProblems/chebyquad.jl +++ b/src/ADNLPProblems/chebyquad.jl @@ -34,15 +34,18 @@ function chebyquad( inv_n = one(T) / n half = one(T) / 2 return half * sum( - (inv_n * sum(chebyshev(x[j], 2i) for j = 1:n) + inv(T((2i)^2 - 1)))^2 for - i = 1:Int(round(m / 2)) + (inv_n * sum(chebyshev(x[j], 2i) for j = 1:n) + one(T) / ((2i)^2 - 1))^2 for + i = 1:div(m, 2) ) + half * sum( - (inv_n * sum(chebyshev(x[j], 2i - 1) for j = 1:n))^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)) + (inv_n * sum(chebyshev(x[j], 2i - 1) for j = 1:n))^2 for i = 1:div(m + 1, 2) ) end step = one(T) / (n + one(T)) - x0 = [j * step for j = 1:n] + x0 = Vector{T}(undef, n) + for j = 1:n + x0[j] = j * step + end return ADNLPModels.ADNLPModel(f, x0, name = "chebyquad"; kwargs...) end @@ -57,8 +60,8 @@ function chebyquad( m = max(m, n) function F!(r, x; n = length(x), m = length(r), chebyshev = chebyshev) inv_n = one(T) / n - for i = 1:Int(round(m / 2)) - r[2i] = inv_n * sum(chebyshev(x[j], 2i) for j = 1:n) + inv(T((2i)^2 - 1)) + for i = 1:div(m, 2) + r[2i] = inv_n * sum(chebyshev(x[j], 2i) for j = 1:n) + one(T) / ((2i)^2 - 1) r[2i - 1] = inv_n * sum(chebyshev(x[j], 2i - 1) for j = 1:n) end if mod(m, 2) == 1 @@ -67,7 +70,10 @@ function chebyquad( return r end step = one(T) / (n + one(T)) - x0 = [j * step for j = 1:n] + x0 = Vector{T}(undef, n) + for j = 1:n + x0[j] = j * step + end return ADNLPModels.ADNLSModel!(F!, x0, m, name = "chebyquad-nls"; kwargs...) end diff --git a/src/PureJuMP/chebyquad.jl b/src/PureJuMP/chebyquad.jl index 60115821f..248a4c804 100644 --- a/src/PureJuMP/chebyquad.jl +++ b/src/PureJuMP/chebyquad.jl @@ -27,40 +27,32 @@ function chebyquad(args...; n::Int = default_nvar, m::Int = n, kwargs...) 0.5 * sum( ( 1 / n * sum( - (begin - # Use shifted Chebyshev argument t = 2x - 1 - t = 2 * x[j] - 1 + ifelse( + 2 * x[j] - 1 ≥ 1, + cosh(2i * acosh(2 * x[j] - 1)), ifelse( - t ≥ 1, - cosh(2i * acosh(t)), - ifelse( - t ≤ -1, - (-1)^(2i) * cosh(2i * acosh(-t)), - cos(2i * acos(t)), - ), - ) - end) for j = 1:n + 2 * x[j] - 1 ≤ -1, + (-1)^(2i) * cosh(2i * acosh(1 - 2 * x[j])), + cos(2i * acos(2 * x[j] - 1)), + ), + ) for j = 1:n ) + 1 / ((2i)^2 - 1) - )^2 for i = 1:Int(round(m / 2)) + )^2 for i = 1:div(m, 2) ) + 0.5 * sum( ( 1 / n * sum( - (begin - # Use shifted Chebyshev argument t = 2x - 1 - t = 2 * x[j] - 1 + ifelse( + 2 * x[j] - 1 ≥ 1, + cosh((2i - 1) * acosh(2 * x[j] - 1)), ifelse( - t ≥ 1, - cosh((2i - 1) * acosh(t)), - ifelse( - t ≤ -1, - (-1)^(2i - 1) * cosh((2i - 1) * acosh(-t)), - cos((2i - 1) * acos(t)), - ), - ) - end) for j = 1:n + 2 * x[j] - 1 ≤ -1, + (-1)^(2i - 1) * cosh((2i - 1) * acosh(1 - 2 * x[j])), + cos((2i - 1) * acos(2 * x[j] - 1)), + ), + ) for j = 1:n ) - )^2 for i = 1:(Int(round(m / 2)) + mod(n, 2)) + )^2 for i = 1:div(m + 1, 2) ) ) return nlp From 4c847c5fb793167b630177c4cb80d2c7ace6e26f Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Mon, 6 Apr 2026 20:59:51 +0530 Subject: [PATCH 7/9] removing extra variables from the chebyquad problem --- src/ADNLPProblems/chebyquad.jl | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl index 58bfa77fa..22310f289 100644 --- a/src/ADNLPProblems/chebyquad.jl +++ b/src/ADNLPProblems/chebyquad.jl @@ -16,8 +16,7 @@ function Cheby(xj, i) return (-1)^i * cosh(i * acosh(-xj)) else # Clamp only for acos to guard against tiny floating-point/AD excursions. - xj_clamped = min(max(xj, -1), 1) - return cos(i * acos(xj_clamped)) + return cos(i * acos(min(max(xj, -1), 1))) end end @@ -31,20 +30,17 @@ function chebyquad( ) where {T} m = max(m, n) function f(x; n = length(x), m = m, chebyshev = chebyshev) - inv_n = one(T) / n - half = one(T) / 2 - return half * sum( - (inv_n * sum(chebyshev(x[j], 2i) for j = 1:n) + one(T) / ((2i)^2 - 1))^2 for + return (one(T) / 2) * sum( + ((one(T) / n) * sum(chebyshev(x[j], 2i) for j = 1:n) + one(T) / ((2i)^2 - 1))^2 for i = 1:div(m, 2) ) + - half * sum( - (inv_n * sum(chebyshev(x[j], 2i - 1) for j = 1:n))^2 for i = 1:div(m + 1, 2) + (one(T) / 2) * sum( + ((one(T) / n) * sum(chebyshev(x[j], 2i - 1) for j = 1:n))^2 for i = 1:div(m + 1, 2) ) end - step = one(T) / (n + one(T)) x0 = Vector{T}(undef, n) for j = 1:n - x0[j] = j * step + x0[j] = j * one(T) / (n + one(T)) end return ADNLPModels.ADNLPModel(f, x0, name = "chebyquad"; kwargs...) end @@ -59,20 +55,18 @@ function chebyquad( ) where {T} m = max(m, n) function F!(r, x; n = length(x), m = length(r), chebyshev = chebyshev) - inv_n = one(T) / n for i = 1:div(m, 2) - r[2i] = inv_n * sum(chebyshev(x[j], 2i) for j = 1:n) + one(T) / ((2i)^2 - 1) - r[2i - 1] = inv_n * sum(chebyshev(x[j], 2i - 1) for j = 1:n) + r[2i] = (one(T) / n) * sum(chebyshev(x[j], 2i) for j = 1:n) + one(T) / ((2i)^2 - 1) + r[2i - 1] = (one(T) / n) * sum(chebyshev(x[j], 2i - 1) for j = 1:n) end if mod(m, 2) == 1 - r[m] = inv_n * sum(chebyshev(x[j], m) for j = 1:n) + r[m] = (one(T) / n) * sum(chebyshev(x[j], m) for j = 1:n) end return r end - step = one(T) / (n + one(T)) x0 = Vector{T}(undef, n) for j = 1:n - x0[j] = j * step + x0[j] = j * one(T) / (n + one(T)) end return ADNLPModels.ADNLSModel!(F!, x0, m, name = "chebyquad-nls"; kwargs...) end From 2c3a49e09ec73aad6f6cafbdcbb68067f74026e3 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Mon, 6 Apr 2026 21:23:04 +0530 Subject: [PATCH 8/9] adding changes to chebyquad.jl --- src/ADNLPProblems/chebyquad.jl | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/src/ADNLPProblems/chebyquad.jl b/src/ADNLPProblems/chebyquad.jl index 22310f289..a5e8e445c 100644 --- a/src/ADNLPProblems/chebyquad.jl +++ b/src/ADNLPProblems/chebyquad.jl @@ -5,19 +5,22 @@ function chebyquad(; use_nls::Bool = false, kwargs...) return chebyquad(Val(model); kwargs...) end -function Cheby(xj, i) - # Use standard Chebyshev definition with robust handling of tiny excursions: - # T_i(x) = cos(i * acos(x)) for |x| ≤ 1 - # T_i(x) = cosh(i * acosh(x)) for x > 1 - # T_i(x) = (-1)^i * cosh(i * acosh(-x)) for x < -1 - if xj > 1 - return cosh(i * acosh(xj)) - elseif xj < -1 - return (-1)^i * cosh(i * acosh(-xj)) - else - # Clamp only for acos to guard against tiny floating-point/AD excursions. - return cos(i * acos(min(max(xj, -1), 1))) - end +function Cheby(xj, i::Integer) + # Evaluate T_i(x) via recurrence to avoid domain-branching and AD/tracer issues. + if i == 0 + return one(xj) + elseif i == 1 + return xj + end + + tk_minus_1 = one(xj) + tk = xj + for _ = 2:i + tk_plus_1 = 2 * xj * tk - tk_minus_1 + tk_minus_1 = tk + tk = tk_plus_1 + end + return tk end function chebyquad( From 3dc0402197d9161fc010d257e23110636935e95c Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Sun, 26 Apr 2026 05:48:17 +0530 Subject: [PATCH 9/9] Update chebyquad.jl --- src/PureJuMP/chebyquad.jl | 62 +++++++++++++++++++-------------------- 1 file changed, 31 insertions(+), 31 deletions(-) diff --git a/src/PureJuMP/chebyquad.jl b/src/PureJuMP/chebyquad.jl index 248a4c804..ba79660c0 100644 --- a/src/PureJuMP/chebyquad.jl +++ b/src/PureJuMP/chebyquad.jl @@ -15,45 +15,45 @@ # # classification SBR2-AN-V-0 export chebyquad + +# Chebyshev polynomial of the first kind via recurrence: T_0=1, T_1=x, +# T_i(x) = 2x*T_{i-1}(x) - T_{i-2}(x). +function _cheby_recurrence(xj::Real, i::Integer) + i == 0 && return one(xj) + i == 1 && return xj + tk_minus_1 = one(xj) + tk = xj + for _ = 2:i + tk_plus_1 = 2 * xj * tk - tk_minus_1 + tk_minus_1 = tk + tk = tk_plus_1 + end + return tk +end + function chebyquad(args...; n::Int = default_nvar, m::Int = n, kwargs...) m = max(m, n) nlp = Model() - x0 = [j/(n + 1) for j = 1:n] + x0 = [j / (n + 1) for j = 1:n] @variable(nlp, x[j = 1:n], start = x0[j]) - # Chebyshev polynomial of the first kind, using explicit expression + + # Register the recurrence evaluator for each required polynomial degree as a + # named univariate JuMP operator (one per degree) + for k = 1:m + op = Symbol("_cheby_$k") + @operator(nlp, op, 1, xj -> _cheby_recurrence(xj, k)) + end + @NLobjective( nlp, Min, 0.5 * sum( - ( - 1 / n * sum( - ifelse( - 2 * x[j] - 1 ≥ 1, - cosh(2i * acosh(2 * x[j] - 1)), - ifelse( - 2 * x[j] - 1 ≤ -1, - (-1)^(2i) * cosh(2i * acosh(1 - 2 * x[j])), - cos(2i * acos(2 * x[j] - 1)), - ), - ) for j = 1:n - ) + 1 / ((2i)^2 - 1) - )^2 for i = 1:div(m, 2) - ) + - 0.5 * sum( - ( - 1 / n * sum( - ifelse( - 2 * x[j] - 1 ≥ 1, - cosh((2i - 1) * acosh(2 * x[j] - 1)), - ifelse( - 2 * x[j] - 1 ≤ -1, - (-1)^(2i - 1) * cosh((2i - 1) * acosh(1 - 2 * x[j])), - cos((2i - 1) * acos(2 * x[j] - 1)), - ), - ) for j = 1:n - ) - )^2 for i = 1:div(m + 1, 2) - ) + (1 / n * sum(nlp[Symbol("_cheby_$(2i)")](x[j]) for j = 1:n) + 1 / ((2i)^2 - 1))^2 + for i = 1:div(m, 2) + ) + 0.5 * sum( + (1 / n * sum(nlp[Symbol("_cheby_$(2i - 1)")](x[j]) for j = 1:n))^2 + for i = 1:div(m + 1, 2) + ), ) return nlp end