# Add nested optimization problem in JuMP

**URL:** <https://discourse.julialang.org/t/add-nested-optimization-problem-in-jump/103378>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [August 30, 2023, 4:47pm UTC](https://discourse.julialang.org/t/add-nested-optimization-problem-in-jump/103378 "2023-08-30T16:47:17Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![leima0521](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leima0521/32/51274_2.png) [@leima0521](https://discourse.julialang.org/u/leima0521)\
**Post date:** [August 30, 2023, 4:47pm UTC](https://discourse.julialang.org/t/add-nested-optimization-problem-in-jump/103378/1 "2023-08-30T16:47:17Z")

</div>

I want to include an inner loop optimization (a fixed point iteration algorithm) inside my current optimization problem. I cannot figure out how to make it work… Any help/guidance/advice is truly appreciated!

The current optimization code is `ipopt_jump` (thanks to the help of @odow in [this post](https://discourse.julialang.org/t/large-memory-usage-in-jump/101491/2)). It tries to find `theta[1:K, 1:L]` and `delta[1:H]` to maximize the objective function. The dimension of `theta` is relatively small, `K=11, L=5`, but the dimension of delta `H` is very big, `H > 2000`. Optimizing over \>2000 parameters takes forever to run.

There is a trick that the value of `delta` can be calculated for any `theta` using a fixed point algorithm. Then I only need to maximize over `K*L = 55` parameters instead of \>2055 parameters.

`delta` can be expressed as a function of `theta`: \delta' = \delta + \ln s\_h - \ln \hat{s}\_h(\delta,\theta) . This will keep iterating until \delta' = \delta. \ln s\_h is observed from the data, and \ln \hat{s}\_h is \hat s\_h(\delta,\theta) = \frac{1}{N} \sum\_{i=1}^N \frac{\exp(\delta\_h + X\theta Z)}{\sum\_h \exp(\delta\_h + X\theta Z)} .

I want to add an inner loop to solve for \delta(\theta) given \theta, then I only need to optimize over \theta. The outer loop is to maximize the objective function in `ipopt_jump`. I tried to incorporate the fixed point algorithm in `ipopt_fp_jump`, but it doesn’t work because of type incompatibility…

A data sample can be found [here](https://www.dropbox.com/scl/fi/hujdy9o22pc6207q2ls9d/example.jld?rlkey=kbk85ugjp59lk9ifmz30wkkhb&dl=0).

```julia
function ipopt_jump(data, iter=10000, tol=1e-3)
    X, Z, choice = data
    N, L = size(Z)
    H, K = size(X)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    set_attributes(model, "tol" => tol, "max_iter" => iter)
    @variable(model, theta[1:K, 1:L])
    @variable(model, delta[1:H])
    @expression(model, v, delta .+ (X * theta * Z'))
    @NLexpression(model, evsum[i = 1:N], sum(exp(v[h, i]) for h = 1:H) + 1)
    choice_h = sum(choice; dims = 1)
    @NLobjective(
        model,
        Max,
        sum(choice[h, i] * v[h, i] for h = 1:H, i = 1:N) - 
        sum(choice_h[i] * log(evsum[i]) for i = 1:N)
    )
    optimize!(model)
    return objective_value(model), value.(theta), value.(delta)
end

function ipopt_fp_jump(data, iter=10000, tol=1e-3)
    X, Z, choice, log_obs_s_h = data
    N, L = size(Z)
    H, K = size(X)
    model = Model(Ipopt.Optimizer)
    set_silent(model)
    set_attributes(model, "tol" => tol, "max_iter" => iter)

    @variable(model, theta[1:K, 1:L])

    # does not work..
    function calc_pred_s_h(θ, δ)
        @expression(model, v, δ .+ (X * θ * Z'))
        @NLexpression(model, evsum[i = 1:N], sum(exp(v[h, i]) for h = 1:H) + 1)
        @NLexpression(model, s_ih[h = 1:H, i = 1:N], exp(v[h, i]) / evsum[i])
        @NLexpression(model, s_h[h = 1:H], sum(s_ih[h, i] for i = 1:N)/N)  
        @NLexpression(model, ln_s_h[h = 1:H], log(s_h[h]))
        return ln_s_h
    end

    # same function as above: does not work..
    # function calc_s_h(θ, δ)
    # ev = exp.(δ .+ (X * θ * Z'))
    # evsum = sum(ev, dims = 1)
    # s_ih = ev ./ evsum
    # s_h = mean(s_ih, dims = 2)
    # return log.(s_h)
    # end
    
    function fixed_point_algorithm(θ)
        old_delta = rand(H)
        normdiff = Inf

        while normdiff > tol
            new_delta = old_delta .+ log_obs_s_h .- calc_pred_s_h(θ, old_delta)
            normdiff = norm(new_delta - old_delta)
            old_delta = new_delta
        end
        
        return new_delta
    end

    delta = fixed_point_algorithm(theta)

    @expression(model, v, delta .+ (X * theta * Z'))
    @NLexpression(model, evsum[i = 1:N], sum(exp(v[h, i]) for h = 1:H) + 1)
    choice_h = sum(choice; dims = 1)

    @NLobjective(
        model,
        Max,
        sum(choice[h, i] * v[h, i] for h = 1:H, i = 1:N) - 
        sum(choice_h[i] * log(evsum[i]) for i = 1:N)
    )
    optimize!(model)
    return objective_value(model), value.(theta), value.(delta)
end

```

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/224680_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [August 30, 2023, 7:22pm UTC](https://discourse.julialang.org/t/add-nested-optimization-problem-in-jump/103378/2 "2023-08-30T19:22:55Z")

</div>

Why does it currently fail?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 30, 2023, 9:21pm UTC](https://discourse.julialang.org/t/add-nested-optimization-problem-in-jump/103378/3 "2023-08-30T21:21:06Z")

</div>

You cannot do what you’re trying to do in `ipopt_fp_jump` because you’re passing in `theta::Vector{VariableRef}`, which tries to compute the fixed point based on the JuMP variable object, _not_ on the _value_ of the variable in a solution.

The usual work-around would be to use a user-defined function, [Nonlinear Modeling · JuMP](https://jump.dev/JuMP.jl/stable/manual/nlp/#User-defined-Functions), but these require scalar outputs.

Nothing immediately comes to mind for a good way to solve this problem with JuMP. Others might have a suggestion.

I’ll note that this problem is unconstrained, so you might have better luck trying something like Optim.jl. if you write things as a pure Julia function given `theta` (including the fixed point calculation), it might just work.

---

<div class="post-metadata">

**Author:** ![gdalle](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gdalle/32/224680_2.png) [@gdalle](https://discourse.julialang.org/u/gdalle)\
**Post date:** [August 31, 2023, 7:19am UTC](https://discourse.julialang.org/t/add-nested-optimization-problem-in-jump/103378/4 "2023-08-31T07:19:03Z")

</div>

Shameless plug: my package [ImplicitDifferentiation.jl](https://github.com/gdalle/ImplicitDifferentiation.jl) is very well-suited to computing [gradients of fixed points](https://gdalle.github.io/ImplicitDifferentiation.jl/stable/examples/1_basic/#Fixed-point). Combining it with a first order algorithm from Optim.jl or Optimization.jl for the outer loop might work

---

<div class="post-metadata">

**Author:** ![leima0521](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/leima0521/32/51274_2.png) [@leima0521](https://discourse.julialang.org/u/leima0521)\
**Post date:** [September 1, 2023, 2:09pm UTC](https://discourse.julialang.org/t/add-nested-optimization-problem-in-jump/103378/5 "2023-09-01T14:09:55Z")

</div>

Hi gdalle. It fails because ln\_s\_h is a general expression that cannot be directly used inside the fixed point function.