Add chebyquad

src/ADNLPProblems/chebyquad.jl

@@ -0,0 +1,79 @@
+export chebyquad
+
+function chebyquad(; use_nls::Bool = false, kwargs...)
+  model = use_nls ? :nls : :nlp
+  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.
+        xj_clamped = min(max(xj, -1), 1)
+        return cos(i * acos(xj_clamped))
+    end
+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)
+    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
+      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)
+    )
+  end
+  step = one(T) / (n + one(T))
+  x0 = Vector{T}(undef, n)
+  for j = 1:n
+    x0[j] = j * step
+  end
+  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)
+    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)
+    end
+    if mod(m, 2) == 1
+      r[m] = inv_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
+  end
+  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...) = max(m, n)

src/PureJuMP/chebyquad.jl

@@ -0,0 +1,59 @@
+#
+#   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.
+#
+#   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(
+            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)
+    )
+  )
+  return nlp
+end