Add chebyquad

src/ADNLPProblems/chebyquad.jl

@@ -0,0 +1,69 @@
+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::Type{T} = 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::Type{T} = 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
+

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

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