Optimization Pushups

The spirit of this simple tutorial consists in learning how to write simple solution algorithms. For each algorithm, test that it works, using simple test functions whose solution is known.

Write a function fixed_point(f::Function, x0::Float64) which computes the fixed point of f starting from initial point x0.

Write a function bisection(f::Function, a::Float64, b::Float64) which computes a zero of function f within (a,b) using a bisection method.

Write a function golden(f::Function, a::Float64, b::Float64) which computes a zero of function f within (a,b) using a golden ratio method.

Write a function zero_newton(f::Function, x0::Float64) which computes the zero of function f starting from initial point x0.

Add an option zero_newton(f::Function, x0::Float64, backtracking=true) which computes the zero of function f starting from initial point x0 using backtracking in each iteration.

Write a function min_gd(f::Function, x0::Float64) which computes the minimum of function f using gradient descent. Assume f returns a scalar and a gradient.

Write a function min_nr(f::Function, x0::Float64) which computes the minimum of function f using Newton-Raphson method. Assume f returns a scalar, a gradient, and a hessian.

Write a method zero_newton(f::Function, x0::Vector{Float64}) which computes the zero of a vector valued function f starting from initial point x0.

f(z) = let 
    x,y = z
    [x^3 + y, x-y] 
end
x0 = [0.1, 0.3]

2-element Vector{Float64}:
 0.1
 0.3

f(x0)

2-element Vector{Float64}:
  0.301
 -0.19999999999999998

using ForwardDiff
ForwardDiff.jacobian(f, x0)

2×2 Matrix{Float64}:
 0.03   1.0
 1.0   -1.0

function zero_newton(f, x0; maxit=10, verbose=false, τ_η=1e-10, τ_ϵ=1e-10)

    # local x1
    x1 = x0
    
    for it ∈ 1:maxit

        f0 = f(x0)
        ϵ = maximum(abs.(f0))
        if ϵ < τ_ϵ
            verbose ? println("Iteration $(it): ϵ=$(ϵ) ") : nothing
            return x0
        end


        J0 = ForwardDiff.jacobian(f, x0)

        x1 = x0 - J0 \ f0

        η = maximum(abs.(x1-x0))
        if η<τ_η
            return x1
        end

        # cond ? do : otherwise

        verbose ? println("Iteration $(it): ϵ=$(ϵ) : η=$(η)") : nothing


        x0 = x1

    end

    return x1

end

zero_newton (generic function with 1 method)

zero_newton(f, x0; verbose=true)

Iteration 1: ϵ=0.301 : η=0.2980582524271845
Iteration 2: ϵ=0.0019417548939487858 : η=0.00194173293071458
Iteration 3: ϵ=1.4642100930805857e-8 : η=1.4642100930805848e-8
Iteration 4: ϵ=6.617444900424222e-24

2-element Vector{Float64}:
 6.617444900424222e-24
 6.617444900424222e-24

Add an method zero_newton(f::Function, x0::Vector{Float64}, backtracking=true) which computes the zero of function f starting from initial point x0 using backtracking in each iteration.

function zero_newton(f, x0; maxit=10, verbose=false, τ_η=1e-10, τ_ϵ=1e-10, backtrack=true)

    # local x1
    x1 = x0
    
    for it ∈ 1:maxit

        f0 = f(x0)
        ϵ = maximum(abs.(f0))
        if ϵ < τ_ϵ
            verbose ? println("Iteration $(it): ϵ=$(ϵ) ") : nothing
            return x0
        end


        J0 = ForwardDiff.jacobian(f, x0)
        
        # Newton steap
        Δ = - J0 \ f0
        
        if !backtrack
        
            x1 = x0 + Δ

        else
            for i = 0:10
                x1 = x0 + 2.0^(-i)*Δ

                # test function for new guess
                f1 = f(x1)
                ϵ1 = maximum(abs.(f1))
                
                if ϵ1<ϵ
                    # new guess is good
                    break
                end

            end
        end

        η = maximum(abs.(x1-x0))
        if η<τ_η
            return x1
        end

        # cond ? do : otherwise

        verbose ? println("Iteration $(it): ϵ=$(ϵ) : η=$(η)") : nothing


        x0 = x1

    end

    return x1

end

zero_newton (generic function with 1 method)

zero_newton(f, x0, verbose=true, backtrack=true)

Iteration 1: ϵ=0.301 : η=0.2980582524271845
Iteration 2: ϵ=0.0019417548939487858 : η=0.00194173293071458
Iteration 3: ϵ=1.4642100930805857e-8 : η=1.4642100930805848e-8
Iteration 4: ϵ=6.617444900424222e-24

2-element Vector{Float64}:
 6.617444900424222e-24
 6.617444900424222e-24

Add a method zero_newton(f::Function, x0::Vector{Float64}, backtracking=true, lb=Vector{Float64}) which computes the zero of function f starting from initial point x0 taking complementarity constraint into account x>=lb using the Fischer-Burmeister method.

g(z) = let
    x,y=z
    [ y-x, x^2+y^2-1 ]
end


lb = [0.0, -10000]

2-element Vector{Float64}:
      0.0
 -10000.0

φ_FB(a,b) = a+b-sqrt(a^2+b^2)  # min

φ_FB (generic function with 1 method)

function ncp(f, x0, lb; args...)
    
    # the dot calls φ_FB component-wise on its vector arguments
    # fun(u) = φ_FB.(f(u), u - lb)

    fun = u-> φ_FB.(f(u), u - lb)

    zero_newton(fun, x0; args...)

end

ncp (generic function with 2 methods)

x0

2-element Vector{Float64}:
 0.1
 0.3

ncp(g, x0, lb; maxit=100, verbose=true)

Iteration 1: ϵ=0.900040498785529 : η=0.8449227125187968
Iteration 2: ϵ=0.34501875850946817 : η=0.17623432346721998
Iteration 3: ϵ=0.04599075111946149 : η=0.022412443670669857
Iteration 4: ϵ=0.0005761770207755035 : η=0.0002880043594630788
Iteration 5: ϵ=8.413735486101359e-8 : η=4.2068675654149956e-8
Iteration 6: ϵ=0.0

2-element Vector{Float64}:
 -5.460886251568741e-17
  1.0000000000002816

ncp(g, [0.6, 0.6], lb; maxit=100, verbose=true)

Iteration 1: ϵ=0.2800039197663864 : η=0.11666503345865797
Iteration 2: ϵ=0.027217503327847226 : η=0.009494525381928298
Iteration 3: ϵ=0.0001802549850253854 : η=6.37240192483679e-5
Iteration 4: ϵ=8.119968697428703e-9 : η=2.8708424526513454e-9
Iteration 5: ϵ=0.0

2-element Vector{Float64}:
 0.7071067811866388
 0.7071067811866388

Other Formats