IPNewton hangs depending on inputs

[juliohm]

Below is a MWE with Optim.jl v1.10 (latest stable release) and Julia v1.11:

using Optim
using Random

gramian(xs, ys; σ=1) = [exp(-euclidsq(x, y) / 2σ^2) for x in xs, y in ys]
euclidsq(x, y) = sum((x[i] - y[i])^2 for i in eachindex(x))

# ------------
# main script
# ------------

rng = MersenneTwister(1234)

σ = 2.0
b = 10

x_nu, x_de = 5randn(rng, 100), randn(rng, 100)

x_ba = x_nu[1:b]
K_nu = gramian(x_nu, x_ba, σ=σ)
K_de = gramian(x_de, x_ba, σ=σ)

# number of numerator and denominator samples
n_nu, n_de = size(K_nu, 1), size(K_de, 1)

# constants for objective function
K = K_nu

# constants for equality constraints
A = sum(K_de, dims=1)
lc = uc = [n_de]

# constants for inequality constraints
T = eltype(K_de)
lx = fill(zero(T), b)
ux = fill(Inf, b)

# objective
f(α) = -sum(log, K * α)
function ∇f!(g, α)
  p = K * α
  for l in 1:b
    g[l] = -sum(K[j, l] / p[j] for j in 1:n_nu)
  end
end
function ∇²f!(h, α)
  p = K * α
  for k in 1:b, l in 1:b
    h[k, l] = sum(view(K, :, k) .* view(K, :, l) ./ p)
  end
end

# equality constraint
c!(c, α) = c .= A * α
J!(J, α) = J .= A
H!(H, α, λ) = H .+= 0

# initial guess
αₒ = fill(n_de / sum(A), b)

# optimization problem
objective = TwiceDifferentiable(f, ∇f!, ∇²f!, αₒ)
constraints = TwiceDifferentiableConstraints(c!, J!, H!, lx, ux, lc, uc)
initguess = αₒ

# solve problem with interior-point primal-dual Newton
solution = optimize(objective, constraints, initguess, IPNewton())

If you change the seed to 123, then IPNewton returns as expected.

[juliohm]

@pkofod any idea on what may be happening here? I believe this issue was introduced by recent releases.