# How is Optim.jl's IPNewton so good compared to (exponentiated) gradient descent?

**URL:** <https://discourse.julialang.org/t/how-is-optim-jls-ipnewton-so-good-compared-to-exponentiated-gradient-descent/85685>\
**Category:** Optimization (Mathematical)\
**Created:** [August 12, 2022, 9:48pm UTC](https://discourse.julialang.org/t/how-is-optim-jls-ipnewton-so-good-compared-to-exponentiated-gradient-descent/85685 "2022-08-12T21:48:52Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [August 12, 2022, 9:48pm UTC](https://discourse.julialang.org/t/how-is-optim-jls-ipnewton-so-good-compared-to-exponentiated-gradient-descent/85685/1 "2022-08-12T21:48:52Z")

</div>

This question is about implementing an optimization algorithm in Julia and comparing it with `IPNewton` from Optim.jl. It’s kind of broad, so not sure if it fits here. Let me know if it doesn’t.

I wrote some code to minimize a function where some parameters need to be on the probability simplex, so this is constrained minimization:

```julia-auto
minimize f(p1, p2 other_stuff)
s.t. p1, p2 >= 0 and p1 + p2 == 1

```

I know that I can substitute `p2 = 1 - p1` and use something like projected gradient descent to satisfy the non-negativity constraints, but that won’t work for higher-dimensional simplexes.

I wanted to compare `IPNewton` (interior point method for constrained optimization) to basic gradient descent. Since some parameters are constrained to the probability simplex, I’m using _exponentiated_ gradient descent (EGD). I felt like using a serious interior-point optimizer for such a simple constraint was overkill, so I wanted to try something simpler, like EGD.

The problem is that EGD is much less consistent, much slower than `IPNewton` and sometimes even produces completely wrong results. The “optimum” found by EGD also results on higher loss than that of `IPNewton`.

**Questions** : is the interior point method just this good? Is my implementation of EGD wrong? Is gradient descent in general this bad for constrained optimization? Am I even allowed to use EGD here? I’ve never seen it used outside online learning.

# Code

The main functions are `minimize_optim` and `minimize_expgrad`.

```julia
import Random
using ComponentArrays
import Optim, ForwardDiff

const AV = AbstractVector{T} where T

# ========== Loss function to minimize ==========
function get_params(params::AV{<:Real})
    K = length(params) ÷ 3

    ax = Axis((p=1:K, mu=(K+1):2K, sigma=(2K+1):3K))
    ComponentVector(params, ax)
end

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

function loss(p::AV{<:Real}, mu::AV{<:Real}, sigma::AV{<:Real}, b::Real, data::AV{<:Real})
    N = length(data)

    penalty = sum(
        p[k] * p[h] * normal_pdf(mu[k], mu[h], 2b + sigma[k]^2 + sigma[h]^2)
        for k in eachindex(p), h in eachindex(p)
    )

    -2/N * sum(
        sum(
            p_k * normal_pdf(x_n, mu_k, 2b + sigma_k^2)
            for (p_k, mu_k, sigma_k) in zip(p, mu, sigma)
        )
        for x_n in data
    ) + penalty
end

function loss(params::AV{<:Real}, b::Real, data::AV{<:Real})
    par = get_params(params)
    loss(par.p, par.mu, par.sigma, b, data)
end

# ========== Minimization algorithms ==========
function constraint_probability!(c::AV{<:Real}, params::AV{<:Real})
    K = length(params) ÷ 3
    p = @view params[1:K]

    c[1] = 1 - sum(p)
end

# Minimize using interior point method from Optim.jl
function minimize_optim(params0::AV{<:Real}, b::Real, data::AV{<:Real}; tol=1e-5)
    K = length(params0) ÷ 3

    # [weights; means; stds]
    # stds can be in R because they're squared later
    lo_params = [zeros(K); fill(-Inf, K); fill(-Inf, K)]
    hi_params = [ones(K); fill(Inf, K); fill(Inf, K)]
    lo_constr = hi_constr = [0.0]

    dfc = Optim.TwiceDifferentiableConstraints(
        constraint_probability!,
        lo_params, hi_params, lo_constr, hi_constr, :forward
    )

    res = Optim.optimize(
        params->loss(params, b, data), dfc,
        params0, Optim.IPNewton(),
        Optim.Options(
            successive_f_tol=20,
            f_tol=tol, f_reltol=tol
        ),
        autodiff=:forward
    )
    Optim.minimizer(res) |> get_params
end

# Minimize using exponentiated gradient descent
function minimize_expgrad(params0::AV{<:Real}, b::Real, data::AV{<:Real}; lr=1e-2, tol=1e-5, max_iter=1_000_000)
    K = length(params0) ÷ 3
    objective = params->loss(params, b, data)

    grad_cfg = ForwardDiff.GradientConfig(objective, params0)

    params = copy(params0)
    gradient = similar(params)
    obj_values = fill(Inf, 50)
    for it in 1:max_iter
        # 1. Compute gradient
        ForwardDiff.gradient!(gradient, objective, params)

        # 2. Exponentiated gradient descent w.r.t. weights
        @. params[1:K] *= exp(-lr * gradient[1:K])
        params[1:K] ./= sum(params[1:K])

        # 3. Regular gradient descent w.r.t. other params
        @. params[K+1:end] -= lr * gradient[K+1:end]

        # 4. Track convergence
        obj_values[1:end-1] .= @view obj_values[2:end]
        obj_values[end] = objective(params)

        metric = if all(isfinite, obj_values)
            @views maximum(
                abs(old - new)
                for (old, new) in zip(obj_values[1:end-1], obj_values[2:end])
            )
        else
            1e6
        end

        if (metric < tol)
            @info "Converged" it
            break
        end
    end

    get_params(params)
end

function random_params(rng, K::Integer)
    p = rand(K)
    p ./= sum(p)

    mu = randn(K)
    sigma = rand(K)

    [p; mu; sigma]
end

# ========== Run minimization ==========
rng = Random.MersenneTwister(42)
data = [randn(rng, 200) .- 3; 0.3 * randn(rng, 400) .+ 2]
params0 = random_params(rng, 2) #[0.5, 0.5, 0,0, 1e-3, 5e-3]
b = 0.01

@info "Minimizing w/ Optim..."
par_optim = minimize_optim(params0, b, data, tol=1e-6)
@show par_optim.p
@show par_optim.mu
@show par_optim.sigma
@show loss(par_optim, b, data)

println()

@info "Minimizing w/ exponentiated gradient descent..."
par_expgrad = minimize_expgrad(params0, b, data, tol=1e-6, lr=1e-3)
@show par_expgrad.p
@show par_expgrad.mu
@show par_expgrad.sigma
@show loss(par_expgrad, b, data)

```

The correct answer is one of these:

```julia-auto
(p=[0.33333, 0.66666], mu=[-3, 2], sigma=[1, 0.3])
(p=[0.66666, 0.33333], mu=[2, -3], sigma=[0.3, 1])

```

# Results

- `IPNewton` produces correct results given any random initial parameters. Moreover, it does so _fast_, literally in a blink of an eye.
- EGD often converges to terrible points. It also takes a couple of seconds to converge. Given random initial points, it mostly converges to correct points, but not as consistently and precisely as `IPNewton`.

### Some examples of bad results by EGD

```julia-auto
julia> params0 = random_params(rng, 2)
6-element Vector{Float64}:
  0.13901488648354637
  0.8609851135164537
 -1.5498971658168026
 -1.9026076357002404
  0.27366742133875555
  0.9996700977843878

julia> minimize_expgrad(params0, b, data, tol=1e-6)
┌ Info: Converged
└ it = 19336
ComponentVector{Float64}(p = [8.126461586000774e-14, 0.9999999999999187], mu = [-1.5702887008520934, 0.5103155146607207], sigma = [0.42607377649097655, 3.0958258461104533])

```

This is terrible because `p[1]=0`, but it shouldn’t be. Given the same initial point, Optim produces results that are correct to the 3rd decimal digit, which is fantastic:

```julia-auto
julia> minimize_optim(params0, b, data, tol=1e-6)
ComponentVector{Float64}(p = [0.3330468090222882, 0.6669531909777119], mu = [-3.088148801196878, 1.9922352674701609], sigma = [0.9469738468169181, 0.3000181581983476])

```

Another one:

```julia-auto
julia> params0 = random_params(rng, 2)
6-element Vector{Float64}:
  0.34854447603600636
  0.6514555239639936
 -0.8641166349035045
  1.136448930313033
  0.45127735273921554
  0.1342420503488576

julia> minimize_expgrad(params0, b, data, tol=1e-6)
┌ Info: Converged
└ it = 12080
ComponentVector{Float64}(p = [0.3460074950904417, 0.6539925049095583], mu = [-2.520558492806493, 1.9930655538825621], sigma = [1.4743322284971923, 0.2953434943385471])

julia> minimize_optim(params0, b, data, tol=1e-6)
ComponentVector{Float64}(p = [0.3330468076989589, 0.666953192301041], mu = [-3.0881488012483316, 1.99223526745685], sigma = [0.9469738440566653, 0.30001815864792564])

```

Here, EGD overestimated `sigma[1]=1.47` and `mu[1]` is also not quite correct. `IPNewton` is spot on, as always.

* * *

**Questions repeated** : is the interior point method just this good? Is my implementation of EGD wrong? Is gradient descent in general this bad for constrained optimization? Am I even allowed to use EGD here? I’ve never seen it used outside online learning.

---

<div class="post-metadata">

**Author:** ![johnmyleswhite](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/johnmyleswhite/32/31_2.png) [@johnmyleswhite](https://discourse.julialang.org/u/johnmyleswhite)\
**Post date:** [August 13, 2022, 12:33am UTC](https://discourse.julialang.org/t/how-is-optim-jls-ipnewton-so-good-compared-to-exponentiated-gradient-descent/85685/2 "2022-08-13T00:33:46Z")

</div>

You’ve got a lot of code going on, so the comparison isn’t trivial, but it looks like you’re doing constant step size gradient descent. That algorithm, despite showing up a lot in ML, is very bad for actually optimizing functions. Were you to gradient descent code using the function in Optim, you’d have a line search, which your code seems to not do.

---

<div class="post-metadata">

**Author:** ![cvanaret](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cvanaret/32/11594_2.png) [@cvanaret](https://discourse.julialang.org/u/cvanaret)\
**Post date:** [August 14, 2022, 5:05pm UTC](https://discourse.julialang.org/t/how-is-optim-jls-ipnewton-so-good-compared-to-exponentiated-gradient-descent/85685/3 "2022-08-14T17:05:53Z")

</div>

I’m honestly not familiar with EGD. However here are a few key points:

- Gradient descent methods work with first-order derivatives. Therefore, they can be attracted to non-minimum stationary points (local max, saddle points);
- Newton-based methods use second-order derivatives, that is they capture local curvature. The local model is therefore more precise. If you’re “close enough” to the solution, you converge very fast (Newton’s method has a quadratic rate of convergence).
- Instead of just moving in the direction opposite the gradient, the descent direction in interior-point methods is obtained by solving a (strictly convex) equality-constrained Quadratic Problem. In practice, this is equivalent to solving a linear system of equations (= the first-order optimality conditions). Therefore, you can guarantee that the direction is a descent direction for your merit function.
- A robust code implements globalization strategies (global convergence = converging to a local minimum whatever the initial point), e.g. a line search + sufficient decrease condition as pointed out by @johnmyleswhite.

Why do you have a penalty term in your code? Have you tried writing explicit constraints instead?

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [August 17, 2022, 9:17pm UTC](https://discourse.julialang.org/t/how-is-optim-jls-ipnewton-so-good-compared-to-exponentiated-gradient-descent/85685/4 "2022-08-17T21:17:31Z")

</div>

I implemented a simple backtracking line search like this:

```julia
# Exponentiated gradient update
function project_exp!(x::AV{<:Real}, direction::AV{<:Real}, lr::Real)::AV{<:Real}
    K = length(x) ÷ 3

    @. x[1:K] *= exp(-lr * direction[1:K])
    x[1:K] ./= sum(x[1:K])

    @. x[K+1:end] -= lr * direction[K+1:end]

    x
end

project_exp(x::AV{<:Real}, direction::AV{<:Real}, lr::Real) =
    project_exp!(copy(x), direction, lr)

function linesearch(f::Function, x::AV{<:Real}, grad::AV{<:Real}, direction::AV{<:Real}, lr0::Real, rho::Real, c::Real)
    @assert lr0 > 0
    @assert 0 < rho < 1
    @assert 0 < c < 1

    f_x = f(x)
    local_slope::Real = grad' * direction
    a = lr0
    while f(project_exp(x, direction, a)) > f_x + a * c * local_slope
        a *= rho
    end

    a
end

```

(I’m not particularly sure if I should be doing `f(project_exp(x, direction, a))` in the loop, since the regular line search uses `x + a*direction`, but I’m trying to satisfy the constraints. Anyway, I have some paper that supposedly explains how to do this properly, so I should read it first.)

Then used the line search inside the optimization loop like this:

```julia
# 1. Compute gradient
ForwardDiff.gradient!(gradient, objective, params, grad_cfg)

# 2. Find optimal step size
lr_opt = linesearch(objective, params, gradient, gradient, 1.0, 0.9, 1e-3)

# 3. Exponentiated gradient descent w.r.t. weights
# and regular gradient descent w.r.t. other params
project_exp!(params, gradient, lr_opt)

```

I tested it with random starting points `params0`, random `data`, and it seems to work great! I’m getting consistent results, nice speed, nice everything! I’m not 100% sure whether this is a fluke or not (because it seems to be such a huge improvement!), but it seems to work much better than without the line search.

TBH, I didn’t think a line search could improve optimization this much…

---

<div class="post-metadata">

**Author:** ![ForceBru](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/forcebru/32/21389_2.png) [@ForceBru](https://discourse.julialang.org/u/ForceBru)\
**Post date:** [August 17, 2022, 9:27pm UTC](https://discourse.julialang.org/t/how-is-optim-jls-ipnewton-so-good-compared-to-exponentiated-gradient-descent/85685/5 "2022-08-17T21:27:45Z")

</div>

I ended up spending this entire week reading Boyd’s “Convex optimization” and also Nocedal, tried to implement my own barrier method, but it… didn’t really have any convergence at all, so I abandoned this idea because Optim.jl exists and I really don’t need to implement my own interior point method. Originally, I wanted my code to be self-contained and remove the dependence on Optim, but now this seems too hard to do properly (EGD is simple, but it doesn’t work as well as `IPNewton`).

The “penalty term” isn’t the penalty in the optimization sense. It’s part of the objective function that I like to call “penalty” for no particular reason.