In [1]:
using LinearAlgebra
function solve_trust_region_subproblem(Bk, g, delta, lam0=-1)
#A pedagogical solver for the trust region subproblem
n = size(Bk, 1)

okay, R, p, q = getchol(Bk, 0, g)
if okay
    if norm(p) < delta
        lam = 0
        return p, lam
    end
end

if lam0 > 0
    lam = lam0
else
    lam = 1
end
first = true
oldlam = lam

for i = 1:100 # maxiter should be an option
    @assert lam > 0
    okay, ~, p, q = getchol(Bk, lam, g)
    if !okay && !first
        #then lam is too small, bisect back to oldlam, which was good
        lam = lam + (oldlam - lam)/2;
        continue
    elseif first && !okay
        #keep increasing until we get somewhere
        found = false
        for j = 1:20 #maxfinditer should be an option
            lam = 2*lam
            okay, ~, p, q = getchol(Bk, lam, g)
            if okay
                found = true
                break
            end
        end
        @assert found
    end
    first = false
    if abs(norm(p) - delta) < 1e-6 #tol should be an option
        return p, lam
    end

    oldlam = lam
    lam = lam + (norm(p)/norm(q))^2 * (norm(p)-delta)/delta
end
error("too many iterations")
end

function getchol(Bk, lam, g)
n = size(Bk, 1)
try
    F = cholesky(Bk + lam*I)
    R = F.U
    p1 = R'\(-g)
    p = R\p1
    q = R'\p
    okay = true
    return okay, R, p, q                
catch e
    okay = false
    R = []
    p = []
    q = []
    return okay, R, p, q
end
end
getchol([1.0 0; 0 1.0], 1.0, [2.0,2.0])
getchol([-2.0 0; 0 1.0], 1.0, [2.0,2.0])
getchol([-2.0 0; 0 1.0], 3.0, [2.0,2.0])
Out[1]:
(true, [1.0 0.0; 0.0 2.0], [-2.0, -0.5], [-2.0, -0.25])
In [5]:
using OptimTestProblems
# These codes turn f and g and H into one function ... 
function opt_combine(x, f, g!, h!)
    g = Array{Float64,1}(undef,length(x))
    g!(g,x)
    H = Array{Float64,2}(undef,length(x),length(x))
    h!(H,x)
    return (f(x), g, H)
end
function opt_problem(p::OptimTestProblems.UnconstrainedProblems.OptimizationProblem)
    return x -> opt_combine(x, p.f, p.g!, p.h!), p
end
# this just makes it easy to use
opt_problem(s::AbstractString) = opt_problem(
    OptimTestProblems.UnconstrainedProblems.examples[s])

# Here's an example 
fgH, p = opt_problem("Rosenbrock")
f, g, H = fgH([0.0,0.0]) # show the function, gradient, and Hessian
p.f([0.0,0.0]) # get the function
Out[5]:
1.0
In [8]:
using Plots
ezcontour(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot(x, y, Z, st=:contour)
end
ezcontour(-1:0.05:1, -1:0.05:1, 
   opt_problem("Rosenbrock")[2])

MethodError: objects of type OptimTestProblems.MultivariateProblems.OptimizationProblem{Nothing, Nothing, Float64, String, Nothing} are not callable

Stacktrace:
 [1] (::var"#9#10"{OptimTestProblems.MultivariateProblems.OptimizationProblem{Nothing, Nothing, Float64, String, Nothing}})(x::Float64, y::Float64)
   @ Main ./In[8]:6
 [2] (::Base.var"#4#5"{var"#9#10"{OptimTestProblems.MultivariateProblems.OptimizationProblem{Nothing, Nothing, Float64, String, Nothing}}})(a::Tuple{Float64, Float64})
   @ Base ./generator.jl:36
 [3] iterate
   @ ./generator.jl:47 [inlined]
 [4] collect
   @ ./array.jl:787 [inlined]
 [5] map
   @ ./abstractarray.jl:3055 [inlined]
 [6] ezcontour(x::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, y::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, f::OptimTestProblems.MultivariateProblems.OptimizationProblem{Nothing, Nothing, Float64, String, Nothing})
   @ Main ./In[8]:6
 [7] top-level scope
   @ In[8]:9

Show an example with the Rosenbrock function¶

In [9]:
# show an example with the Rosenbrock function
x0 = [-1.0,0.5]
rx = -1.2:0.01:-0.8
ry = 0.2:0.01:1
fgH, p = opt_problem("Rosenbrock")
ezcontour(rx,ry,p.f)
Out[9]:
In [10]:
deltas = 0.01:0.01:1
f0,g,H = fgH(x0)
ps = map(d -> solve_trust_region_subproblem(H,g,d)[1], deltas)
ezcontour(rx,ry,x -> log10.(p.f(x)))
scatter!(map(first, ps).+x0[1], map(x->x[2], ps).+x0[2], 
        marker_z = deltas,  label="TR")
newton = -H\g
plot!([x0[1],x0[1] + newton[1]], [x0[2],x0[2]+newton[2]], label="Newton")
plot!([x0[1],x0[1] - g[1]], [x0[2],x0[2]-g[2]], label="GD")
xlims!(rx[1],rx[end])
ylims!(ry[1],ry[end])
Out[10]:
In [11]:
## Show the contours of the Trust region function
trmodel(x) = f0 .+ (x.-x0)'*g .+ 1.0/2.0*(x.-x0)'*H*(x.-x0)
ezcontour(rx,ry,x -> log10.(trmodel(x)))
Out[11]:
In [12]:
## Show the Cauchy point as Delta varies
ps = map(d -> -g/norm(g).*d.*min.(norm(g)^3.0./(d.*g'*H*g),d), deltas)
scatter!(map(first, ps).+x0[1], map(x->x[2], ps).+x0[2], 
        marker_z = deltas,  label="TR")
Out[12]:
In [13]:
(-g'*g ./ (g'*H*g)).*g
Out[13]:
2-element Vector{Float64}:
 0.1754381662197877
 0.08599910108813122
In [15]:
using Pkg; Pkg.add("Polynomials")

    Updating registry at `~/.julia/registries/General`
    Updating git-repo `https://github.com/JuliaRegistries/General.git`
   Resolving package versions...
   Installed Polynomials ─ v3.2.8
    Updating `~/.julia/environments/v1.8/Project.toml`
  [f27b6e38] + Polynomials v3.2.8
    Updating `~/.julia/environments/v1.8/Manifest.toml`
  [f27b6e38] ↑ Polynomials v3.2.4 ⇒ v3.2.8
Precompiling project...
  ✓ Polynomials
  ✓ DSP
  ✓ ToeplitzMatrices
  ✓ FastTransforms
  ✓ InfiniteLinearAlgebra
  ✓ ApproxFunBase
  ✓ ApproxFunFourier
  ✓ ApproxFunOrthogonalPolynomials
  ✓ ApproxFunSingularities
  ✓ ApproxFun
  10 dependencies successfully precompiled in 49 seconds. 408 already precompiled. 5 skipped during auto due to previous errors.
In [17]:
using Polynomials
function dogleg(delta,H,g)
    pfull = -H\g 
    pu = (-g'*g ./ (g'*H*g)).*g
    pdiff = pfull .- pu
    if norm(pfull) < delta
        return pfull
    elseif norm(pu) > delta
        return pu/norm(pu)*delta 
    else # find the intersection with the dogleg path
        # test 
        tmin = roots(ImmutablePolynomial([pu'*pu - delta*delta,2*pdiff'*pu,pdiff'*pdiff]))
        filter!(x -> x >= 0 && x <= 1, tmin)                
        @assert length(tmin) == 1
        tmin1 = tmin[1]
        return pu + (tmin1)*pdiff        
    end            
end

ps = map(d -> dogleg(d,H,g), deltas)
ezcontour(rx,ry,x -> log10.(p.f(x)))        
scatter!(map(first, ps).+x0[1], map(x->x[2], ps).+x0[2], 
        marker_z = deltas,  label="Dogleg")
Out[17]:

Try on a crazier function¶

In [7]:
x0 = [0,3.6]
rx = -1:0.01:1
ry = 2:0.01:4
fgH, p = opt_problem("Himmelblau")
ezcontour(rx,ry,p.f)
Out[7]:
-1.0 -0.5 0.0 0.5 1.0 2.0 2.5 3.0 3.5 4.0 50 60 70 80 90 100 110 120 130
In [8]:
deltas = 0.01:0.01:1
f0,g,H = fgH(x0)
ps = map(d -> solve_trust_region_subproblem(H,g,d)[1], deltas)
ezcontour(rx,ry,x -> (p.f(x)))
scatter!(map(first, ps).+x0[1], map(x->x[2], ps).+x0[2], 
        marker_z = deltas,  label="TR")
newton = -H\g
plot!([x0[1],x0[1] + newton[1]], [x0[2],x0[2]+newton[2]], label="Newton")
plot!([x0[1],x0[1] - g[1]], [x0[2],x0[2]-g[2]], label="GD")
xlims!(rx[1],rx[end])
ylims!(ry[1],ry[end])
Out[8]:
-1.0 -0.5 0.0 0.5 1.0 2.0 2.5 3.0 3.5 4.0 25 50 75 100 125 25 50 75 100 125 TR Newton GD
In [140]:
using Polynomials
function dogleg(delta,H,g)
    pfull = -H\g 
    pu = (-g'*g ./ (g'*H*g)).*g
    pdiff = pfull .- pu
    #println(norm(pfull))
    #println(norm(pu))
    if norm(pfull) < delta
        return pfull
    elseif norm(pu) > delta
        return (pu/norm(pu))*delta 
    else # find the intersection with the dogleg path
        # test 
        tmin = roots(ImmutablePolynomial([pu'*pu - delta*delta,2*pdiff'*pu,pdiff'*pdiff]))
        filter!(x -> x >= 0 && x <= 1, tmin)                
        @assert length(tmin) == 1
        tmin1 = tmin[1]
        return pu .+ (tmin1)*pdiff        
    end            
end

ps = map(d -> dogleg(d,H,g), deltas)
ezcontour(rx,ry,x -> log10.(p.f(x)))
scatter!(map(first, ps).+x0[1], map(x->x[2], ps).+x0[2], 
        marker_z = deltas,  label="Dogleg")
Out[140]:
-1.0 -0.5 0.0 0.5 1.0 2.0 2.5 3.0 3.5 4.0 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 Dogleg