Fixing a parameter with Optimization.jl

Suppose I define the following problem which computes the maximum likelihood estimates (MLEs) for a regression problem using Optimization.jl:

``````using Optimization, OptimizationOptimJL, Random, Distributions
β₀, β₁, β₂, β₃, n, σ = -1.0, 1.0, 0.5, 3.0, 300, 0.1
Random.seed!(12301)
x₁ = rand(Uniform(-1, 1), n)
x₂ = rand(Normal(1.0, 0.5), n)
ε = rand(Normal(0.0, σ), n)
y = @. β₀ + β₁ * x₁ + β₂ * x₂ + β₃ * x₁ * x₂ + ε
function loglik(θ, p)
β = @views θ[1:4]
σ² = θ[5]
y, x₁, x₂ = p
ℓℓ = -n / 2 * log(2π * σ²)
for i in 1:length(y)
ℓℓ = ℓℓ - 0.5 / σ² * (y[i] - β[1] - β[2] * x₁[i] - β[3] * x₂[i] - β[4] * x₁[i] * x₂[i])^2
end
return ℓℓ
end
func = OptimizationFunction((u, p) -> -loglik(u, p), Optimization.AutoForwardDiff())
prob = OptimizationProblem(func, ones(5), (y, x₁, x₂); lb = [-Inf,-Inf,-Inf,-Inf,0.0],ub=[Inf,Inf,Inf,Inf,Inf])
``````

Is there an easy way to extend the existing `OptimizationProblem`, `prob`, to fix a certain parameter? I’d want to fix say `σ` and optimise for the other parameters, over a large range of `σ` so it would have to be reasonably efficient. (This is computing a profile likelihood.) This would be done for many parameters, and the computations that will follow it will require the MLE for the full problem, so I do need this full `prob` to start with.

I thought about just fixing `lb` and `ub` for this parameter so that `lb = ub = the fixed value of σ`, where I’d use `remake` to make this change to `prob`, but this seems like it could cause problems. Any ideas?

Do `(u, p) -> -loglik(u, p, σ²)` and enclose it, then just optimize over 4 values.

1 Like

Thank you @ChrisRackauckas. Does that not cause significant overhead if this is done many times since `OptimizationFunction` and `OptimizationProblem` have to be reconstructed everytime? Maybe it’s unavoidable anyway for this type of thing.

The value is fixed right? So the optimization should be able to run without constructing any now OptimizationX

1 Like

Sorry, by “many times” I mean doing this same thing for say `σ² ∈ 0.01:0.01:1.00`, so that a new problem is constructed for each individual `σ²` and then optimised, and same for all the next values (so 100 problems in this example).

ParameterHandling.jl is a fantastic package for managing optimization parameters. It allows you to assign bounds and fix parameters, among other things. It even automatically handles scaling parameters. I have used it alongside Optimization.jl to great effect.

3 Likes

This package looks very useful, I’ll definitely look into it. Do you have any good examples of how you’ve used it alongside `Optimization.jl`? If you don’t mind sharing; no problem otherwise.

1 Like

Sure, here’s a version of your example using ParameterHandling:

``````using ParameterHandling
using Optimization
using OptimizationOptimJL
using Random
using Distributions

β₀, β₁, β₂, β₃, n, σ = -1.0, 1.0, 0.5, 3.0, 300, 0.1
Random.seed!(12301)
x₁ = rand(Uniform(-1, 1), n)
x₂ = rand(Normal(1.0, 0.5), n)
ε = rand(Normal(0.0, σ), n)
y = @. β₀ + β₁ * x₁ + β₂ * x₂ + β₃ * x₁ * x₂ + ε

params = (σ² = bounded(1.0, 0.0, 1e9), β = ones(4)) # β doesn't have bounds, just initial value
vals, unflatten = value_flatten(params)
function loglik(θ, data)
(; σ², β) = θ
y, x₁, x₂ = data
ℓℓ = -n / 2 * log(2π * σ²)
for i in 1:length(y)
ℓℓ = ℓℓ - 0.5 / σ² * (y[i] - β[1] - β[2] * x₁[i] - β[3] * x₂[i] - β[4] * x₁[i] * x₂[i])^2
end
return ℓℓ
end
func = OptimizationFunction((θ, data) -> -loglik(unflatten(θ), data), Optimization.AutoForwardDiff())
prob = OptimizationProblem(func, vals, (y, x₁, x₂))
soln = unflatten(collect(θhat))

# Fix σ²:
params = (σ² = fixed(1.0), β = ones(4))
vals, unflatten = value_flatten(params)
func = OptimizationFunction((θ, data) -> -loglik(unflatten(θ), data), Optimization.AutoForwardDiff())
prob = OptimizationProblem(func, vals, (y, x₁, x₂))
soln = unflatten(collect(θhat))
``````

Note that I prefer to call the second argument in the objective function `data` and keep the term `params` for parameters like θ to be estimated in the statistical sense.

1 Like

Very slick, thanks so much for that. Can I just ask what is the purpose of this quoted code? I would just write `σ², β = θ` (as you do for the line `y, x₁, x₂ = data` - is there a difference?

Sure, the `(; σ², β) = θ` syntax does multiple assignment by fieldnames. This can be useful for unpacking NamedTuples like I have above.

In contrast, `σ², β = θ` only unpacks by order, and can be used even for non-named objects. This is a bit more fragile. I prefer to use names instead of relying on order when possible.

edit: I couldn’t find anything in the docs about `(; σ², β) = θ`, but it’s the first new language feature listed in the 1.7 release notes: Julia v1.7 Release Notes · The Julia Language

1 Like