Interior-point Newton seems to be non-deterministic

[ForceBru]

Code

import Pkg
Pkg.activate(temp=true, io=devnull)
Pkg.add(name="Optim", version="1.6.0", io=devnull)
Pkg.status()

const AV = AbstractVector{T} where T
using DelimitedFiles, Optim

normal_pdf(x::Real, μ::Real, var::Real) =
    exp(-(x-μ)^2 / (2var)) / sqrt(2π * var)

function D(θ::AV{<:Real}, b::Real, r::AV{<:Real})
	@assert length(θ) % 3 == 0

	N = length(r)
	K = length(θ) ÷ 3
	α, μ, σ = θ[1:K], θ[K+1:2K], θ[2K+1:end]
	var = σ .^2

	sum(
		α[k] * (
			sum(
				α[h] * normal_pdf(μ[k], μ[h], 2b + var[k] + var[h])
				for h ∈ 1:K
			)
			- 2*sum(
				1/N * normal_pdf(μ[k], r[j], 2b + var[k])
				for j ∈ eachindex(r)
			)
		)
		for k ∈ 1:K
	)
end

data = open("sample_data.csv", "r") do io
	DelimitedFiles.readdlm(io, ',')
end
data = data[:, 1]

result = let
	θ0 = Float64[1/3,1/3,1/3, 0,0,0, 1,1,1]
	obj(θ::AV{<:Real}) = D(θ, 1.0, data)

	con_c!(c, θ::AV{<:Real}) = (c[1] = sum(θ[1:(length(θ)÷3)]) - 1; c)
	con_jac!(J, θ::AV{<:Real}) = (J[1, 1:(length(θ)÷3)] .= 1; J)
	con_h!(H, θ, λ) = H

	lx = Float64[0,0,0, -Inf,-Inf,-Inf, 0,  0,  0]
	ux = Float64[1,1,1,  Inf, Inf, Inf, Inf,Inf,Inf]
	lc = [0.0]; uc = [0.0]
	dfc = TwiceDifferentiableConstraints(
		con_c!, con_jac!, con_h!,
		lx, ux, lc, uc
	)
	optimize(
		obj, dfc, θ0, IPNewton(),
		Optim.Options(store_trace=true, show_trace=true, show_every=1, iterations=10),
		autodiff=:forward
	)
end

@show Optim.minimizer(result)

• The objective function is obj(θ::AV{<:Real}) = D(θ, 1.0, data)
• One affine constraint: θ[1] + θ[2] + θ[3] - 1 == 0
• Inequality constraints:
  • θ[1]>=0 && θ[2]>=0 && θ[3]>=0
  • θ[1]<=1 && θ[2]<=1 && θ[3]<=1
  • θ[7]>=0 && θ[8]>=0 && θ[9]>=0
• Starting point: θ0 = Float64[1/3,1/3,1/3, 0,0,0, 1,1,1]

100 data points: sample_data.csv

I'm not particularly sure about the specification of the Hessian in con_h!(H, θ, λ) = H. AFAIK, this function should add the Hessian of the constraints to H. Since the constraints are linear, I should be adding zero, which is equivalent to doing nothing and simply returning the original Hessian, right?

Problem

If I run this code multiple times (with the same starting point and the same settings!), I randomly get NaNs everywhere:

~/t/Optim $ julia-1.7 report.jl                                                                                                               (base) 
      Status `/private/var/folders/ys/3h0gnqns4b98zb66_vl_m35m0000gn/T/jl_uNRF2Q/Project.toml`
  [429524aa] Optim v1.6.0
Iter     Lagrangian value Function value   Gradient norm    |==constr.|      μ
     0   -1.966311e-01    -1.966871e-01    1.096659e-03     0.000000e+00     1.24e-05
 * time: 1.642772858940843e9
     1   -1.967483e-01    -1.967538e-01    5.965825e-04     4.191989e-21     1.24e-06
 * time: 1.642772860832457e9
     2   -1.967599e-01    -1.967656e-01    1.449232e-03     0.000000e+00     1.24e-06
 * time: 1.642772860833223e9
     3   -1.967671e-01    -1.967714e-01    8.505952e-05     3.350116e-23     9.33e-07
 * time: 1.642772860833715e9
     4   -1.967796e-01    -1.967802e-01    4.338795e-04     1.761310e-22     1.30e-07
 * time: 1.642772860834255e9
     5   -1.967911e-01    -1.967917e-01    4.314652e-04     1.548437e-22     1.25e-07
 * time: 1.642772860834772e9
     6   -1.967975e-01    -1.967984e-01    6.682571e-05     2.934913e-22     1.63e-07
 * time: 1.642772860835274e9
     7   -1.968029e-01    -1.968034e-01    1.454685e-04     6.996569e-23     9.13e-08
 * time: 1.64277286083576e9
     8   -1.968043e-01    -1.968048e-01    1.681820e-04     3.166961e-22     8.18e-08
 * time: 1.642772860836246e9
     9   -1.968060e-01    -1.968065e-01    4.130831e-04     2.883045e-22     7.89e-08
 * time: 1.642772860836716e9
    10   -1.968087e-01    -1.968091e-01    6.875015e-04     1.149389e-22     7.09e-08
 * time: 1.642772860837178e9
Optim.minimizer(result) = [0.12046220033401919, 0.21628965837338163, 0.6632481412925991, 0.3006553011305413, 0.36322704596103234, -0.15810174712000133, 1.3014704098568006, 0.561227621176915, 1.0813300576532012]
~/t/Optim $ julia-1.7 report.jl                                                                                                               (base) 
      Status `/private/var/folders/ys/3h0gnqns4b98zb66_vl_m35m0000gn/T/jl_JcPxe4/Project.toml`
  [429524aa] Optim v1.6.0
Iter     Lagrangian value Function value   Gradient norm    |==constr.|      μ
     0   NaN              -1.966871e-01    NaN              NaN              1.00e+00
 * time: 1.642772890437607e9
     1   NaN              NaN              NaN              NaN              1.00e-01
 * time: 1.642772892260335e9
     2   NaN              NaN              NaN              NaN              NaN   
 * time: 1.64277289226151e9
     3   NaN              NaN              NaN              NaN              NaN   
 * time: 1.642772892262521e9
     4   NaN              NaN              NaN              NaN              NaN   
 * time: 1.642772892263575e9
     5   NaN              NaN              NaN              NaN              NaN   
 * time: 1.642772892264575e9
     6   NaN              NaN              NaN              NaN              NaN   
 * time: 1.642772892265608e9
     7   NaN              NaN              NaN              NaN              NaN   
 * time: 1.642772892266602e9
     8   NaN              NaN              NaN              NaN              NaN   
 * time: 1.642772892267635e9
     9   NaN              NaN              NaN              NaN              NaN   
 * time: 1.642772892268625e9
    10   NaN              NaN              NaN              NaN              NaN   
 * time: 1.642772892269617e9
Optim.minimizer(result) = [NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN]
~/t/Optim $

Looks like NaN first appears in Lagrangian value, and then everything becomes NaN.

Is this expected behavior? Am I maybe specifying constraints incorrectly?

[gwater]

It looks like you are following this Tutorial — is that correct? https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/#generic-nonlinear-constraints

[gwater]

I suspect the following: IPNewton allocates some memory without zero-ing it. (eg here

Optim.jl/src/multivariate/solvers/constrained/ipnewton/ipnewton.jl

Line 111 in cdebfca

constr_c = Array{T}(undef, mc)

)That memory can contain NaN, infinities, etc.

Depending on the way you define your methods con_c!, con_jac! etc that undefined memory can "slip through" without ever getting properly set up and then pollutes the rest of the data (because NaN propagate, for example in sums).

[ForceBru]

It looks like you are following this Tutorial — is that correct? https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/#generic-nonlinear-constraints

Yes, that's right.

IPNewton allocates some memory without zero-ing it

Looks like this is the issue. I initialized the Jacobian to zeros, and now everything works fine:

con_jac!(J, θ::AV{<:Real}) = (J .= 0; J[1, 1:(length(θ)÷3)] .= 1; J)

---

I think the documentation could explain more clearly how users are expected to fill the Jacobian and the Hessian:

• What's the value of the incoming Jacobian? Say, in con_jac!(J, θ::AV{<:Real}) = ..., what's the initial value of J? Turns out, it's actually random (IMHO, this should be either explicitly stated in the docs or changed to zero-initialization), so the user must fill the entire matrix with valid values - otherwise they can get random results or failures as I did in this issue.
• Filling out the Hessian manually is extremely error-prone:
  • It's confusing that I have to update the Hessian and manually multiply my second derivatives by some lambdas. I understand that lambdas come from the Lagrangian, but IMO, the optimization library probably should have another simpler interface where the user doesn't need to care about lambdas.
  • It's unclear why Optim can't calculate the Hessian of constraints itself. It can easily differentiate the objective function, but constraints, which are no different - just more functions! - need special treatment for no apparent reason. Why not have the user provide multiple functions like constr1(theta), constr2(theta) etc, and then Optim could differentiate them and build the Hessian - all by itself?

[gwater]

After looking into it more deeply: These allocations all happen in initial_state() which means they only happen once. I think it will not affect performance too much to replace them with zeros(…).

[pkofod]

I suspect the following: IPNewton allocates some memory without zero-ing it. (eg here

Yes, this is bad. I have changed undef initialization everywhere else. That's an oversight.

Depending on the way you define your methods con_c!, con_jac! etc that undefined memory can "slip through" without ever getting properly set up and then pollutes the rest of the data (because NaN propagate, for example in sums).

Yes, true. In other places (ie NLSolversBase) I initialize H, J caches to NaN for this to happen with guarantee. That way, it's more obvious if you forgot to update all elements: you'll get an error right away. The undef-way makes this random which is bad.

Looks like this is the issue. I initialized the Jacobian to zeros, and now everything works fine:

This is what we expect users to do, but maybe it isn't clear enough. We re-use the caches between iterations, so it's 100% on the user to make sure all elements are set correctly. This has not caused any issues or a significant number of bug reports in the past (that doesn't mean we don't have to change of course).

• What's the value of the incoming Jacobian?

You cannot assume that you can re-use any element of any cache. The idea is that the user has to set everything if the caches are dense, even "structural zeros" because these are not tracked in a dense array.

• It's confusing that I have to update the Hessian and manually multiply my second derivatives by some lambdas. I understand that lambdas come from the Lagrangian, but IMO, the optimization library probably should have another simpler interface where the user doesn't need to care about lambdas.

I accept the critique. We could consider changing it.

It's unclear why Optim can't calculate the Hessian of constraints itself.

It can, although I realize now that I merged that PR thinking that autodiff and chunk were kwargs, but they are positional..

θ0 = Float64[1/3,1/3,1/3, 0,0,0, 1,1,1]
	obj(θ::AV{<:Real}) = D(θ, 1.0, data)

	function con_c!(c, θ::AV{<:Real})
		c .= 0
		c[1] = sum(θ[1:(length(θ)÷3)]) - 1
		c
	end

	lx = Float64[0,0,0, -Inf,-Inf,-Inf, 0,  0,  0]
	ux = Float64[1,1,1,  Inf, Inf, Inf, Inf,Inf,Inf]
	lc = [0.0]; uc = [0.0]
	dfc_fin = TwiceDifferentiableConstraints(
		con_c!, lx, ux, lc, uc
	);
	dfc_ad = TwiceDifferentiableConstraints(
		con_c!, lx, ux, lc, uc, :forwarddiff
	);

	optimize(
		obj, dfc_fin, θ0, IPNewton(),
		Optim.Options(store_trace=true, show_trace=true, show_every=1, iterations=10),
		autodiff=:forward
	)
		optimize(
		obj, dfc_ad, θ0, IPNewton(),
		Optim.Options(store_trace=true, show_trace=true, show_every=1, iterations=10),
		autodiff=:forward
	)

[ForceBru]

It's unclear why Optim can't calculate the Hessian of constraints itself.

It can

Wow, that's a game changer! I've been calculating Hessians by hand for way too long, and now I'm free! (Right?)

---

However, the documentation about autodiff in Optim never mentions that it's possible to use autodiff for constraints. The IPNewton tutorial doesn't mention this either, even though autodiff seems to be most helpful when one needs to compute higher-order derivatives, like the Hessian.

I think it would be helpful if the documentation explained that it's possible to use autodiff with constraints.

[pkofod]

Yes, certainly, it the documentation should be improved. It think it would also be useful to add a constraints_autodiff=... at the optimize level.