Merge 5591bf5 into b75665d

JuliaSmoothOptimizers · Nov 18, 2017 · 737f1d1 · 737f1d1
2 parents b75665d + 5591bf5
commit 737f1d1
Show file tree

Hide file tree

Showing 4 changed files with 75 additions and 71 deletions.
diff --git a/src/solver/tron.jl b/src/solver/tron.jl
@@ -5,6 +5,8 @@ using LinearOperators, NLPModels
 
 export tron
 
+tronlogger = get_logger("optimize.tron")
+
 """`tron(nlp)`
 
 A pure Julia implementation of a trust-region solver for bound-constrained
@@ -21,7 +23,6 @@ function tron(nlp :: AbstractNLPModel;
               μ₀ :: Real=1e-2,
               μ₁ :: Real=1.0,
               σ :: Real=10.0,
-              verbose :: Bool=false,
               itmax :: Int=10_000 * nlp.meta.nvar,
               max_cgiter :: Int=nlp.meta.nvar,
               cgtol :: Real=0.1,
@@ -61,35 +62,44 @@ function tron(nlp :: AbstractNLPModel;
   tired = iter >= itmax || el_time > timemax
   unbounded = fx < fmin
   stalled = false
-  status = ""
+  status = "unknown"
 
   αC = 1.0
   tr = TRONTrustRegion(min(max(1.0, 0.1 * norm(πx)), 100))
-  if verbose
-    @printf("%4s  %9s  %7s  %7s  %7s  %s\n", "Iter", "f", "π", "Radius", "Ratio", "CG-status")
-    @printf("%4d  %9.2e  %7.1e  %7.1e\n", iter, fx, πx, get_property(tr, :radius))
-  end
+  @info(tronlogger,
+        @sprintf("%4s  %9s  %7s  %7s  %7s  %s",
+                 "Iter", "f", "π", "Radius", "Ratio", "CG-status"))
+  @info(tronlogger,
+        @sprintf("%4d  %9.2e  %7.1e  %7.1e",
+                 iter, fx, πx, get_property(tr, :radius)))
   while !(optimal || tired || stalled || unbounded)
     # Current iteration
     xc .= x
     fc = fx
     Δ = get_property(tr, :radius)
     H = hess_op!(nlp, xc, temp)
 
-    αC, s = cauchy(x, H, gx, Δ, αC, ℓ, u, μ₀=μ₀, μ₁=μ₁, σ=σ)
+    αC, s, cauchy_status = cauchy(x, H, gx, Δ, αC, ℓ, u, μ₀=μ₀, μ₁=μ₁, σ=σ)
 
+    if cauchy_status != "success"
+      @critical(tronlogger, "Cauchy step returned: $cauchy_status")
+      status = cauchy_status
+      stalled = true
+      continue
+    end
     s, Hs, cgits, cginfo = projected_newton!(x, H, gx, Δ, cgtol, s, ℓ, u, max_cgiter=max_cgiter)
     slope = dot(gx, s)
     qs = 0.5 * dot(s, Hs) + slope
     fx = f(x)
 
-    try
-      ratio!(tr, nlp, fc, fx, qs, x, s, slope)
-    catch exc
-      status = exc.msg
+    ared, pred, quad_min = aredpred(tr, nlp, fc, fx, qs, x, s, slope)
+    if pred ≥ 0
+      status = "nonnegative predicted reduction"
       stalled = true
       continue
     end
+    tr.ratio = ared / pred
+    tr.quad_min = quad_min
 
     s_norm = norm(s)
     if num_success_iters == 0
@@ -119,21 +129,17 @@ function tron(nlp :: AbstractNLPModel;
     optimal = πx <= ϵ
     unbounded = fx < fmin
 
-    verbose && @printf("%4d  %9.2e  %7.1e  %7.1e  %7.1e  %s\n", iter, fx, πx, Δ, Ared / Pred, cginfo)
+    @info(tronlogger,
+          @sprintf("%4d  %9.2e  %7.1e  %7.1e  %7.1e  %s",
+                   iter, fx, πx, Δ, tr.ratio, cginfo))
   end
 
-  status = if tired
-    iter >= itmax ? "maximum number of iterations" : "maximum elapsed time"
+  if tired
+    status = iter >= itmax ? "maximum number of iterations" : "maximum elapsed time"
   elseif optimal
-    "first-order stationary"
+    status = "first-order stationary"
   elseif unbounded
-    "objective function is unbounded from below"
-  elseif status != ""
-    status
-  elseif stalled
-    "stalled"
-  else
-    "unknown"
+    status = "objective function is unbounded from below"
   end
 
   return x, fx, πx, iter, optimal, tired, status, el_time
@@ -206,6 +212,8 @@ function cauchy{T <: Real}(x::AbstractVector{T},
   s = zeros(n)
   Hs = zeros(n)
 
+  status = "success"
+
   project_step!(s, x, g, ℓ, u, -α)
 
   # Interpolate or extrapolate
@@ -229,10 +237,9 @@ function cauchy{T <: Real}(x::AbstractVector{T},
       end
       # TODO: Correctly assess why this fails
       if α < sqrt(nextfloat(zero(α)))
-        throw("α too small (qs = $qs, slope = $slope)")
-        #stalled = true
-        #status = "α too small"
-        #break
+        stalled = true
+        search = false
+        status = "Failure to achieve sufficient decrease: α too small"
       end
     end
   else
@@ -255,7 +262,7 @@ function cauchy{T <: Real}(x::AbstractVector{T},
     α = αs
     s = project_step!(s, x, g, ℓ, u, -α)
   end
-  return α, s
+  return α, s, status
 end
 
 """`projected_newton!(x, H, g, Δ, gctol, s, max_cgiter, ℓ, u)`

diff --git a/src/solver/trunk.jl b/src/solver/trunk.jl
@@ -13,11 +13,6 @@
 
 export trunk
 
-"Exception type raised in case of error inside Trunk."
-type TrunkException <: Exception
-  msg  :: String
-end
-
 trunklogger = get_logger("optimize.trunk")
 
 function trunk(nlp :: AbstractNLPModel;
@@ -52,10 +47,10 @@ function trunk(nlp :: AbstractNLPModel;
   nm_iter = 0
 
   # Preallocate xt.
-  xt = Array{Float64}(n)
-  temp = Array{Float64}(n)
+  xt = Vector{Float64}(n)
+  temp = Vector{Float64}(n)
 
-  optimal = ∇fNorm2 <= ϵ
+  optimal = ∇fNorm2 ≤ ϵ
   tired = nlp.counters.neval_obj > max_f
   stalled = false
 
@@ -87,16 +82,17 @@ function trunk(nlp :: AbstractNLPModel;
     Δq = slope + 0.5 * curv
     ft = obj(nlp, xt)
 
-    try
-      ratio!(tr, nlp, f, ft, Δq, xt, s, slope)
-    catch exc # negative predicted reduction
-      status = exc.msg
+    ared, pred = aredpred(nlp, f, ft, Δq, xt, s, slope)
+    if pred ≥ 0
+      status = "nonnegative predicted reduction"
       stalled = true
       continue
     end
+    tr.ratio = ared / pred
 
     if !monotone
-      ρ_hist = ratio(nlp, fref, ft, σref + Δq, xt, s, slope)
+      ared_hist, pred_hist = aredpred(nlp, fref, ft, σref + Δq, xt, s, slope)
+      ρ_hist = ared_hist / pred_hist
       set_property!(tr, :ratio, max(get_property(tr, :ratio), ρ_hist))
     end
 
@@ -110,7 +106,11 @@ function trunk(nlp :: AbstractNLPModel;
       # slope *= get_property(tr, :radius) / sNorm
       # sNorm = get_property(tr, :radius)
 
-      slope < 0.0 || throw(TrunkException(@sprintf("not a descent direction: slope = %9.2e, ‖∇f‖ = %7.1e", slope, ∇fNorm2)))
+      if slope ≥ 0.0
+        status = @sprintf("not a descent direction: slope = %9.2e, ‖∇f‖ = %7.1e", slope, ∇fNorm2)
+        stalled = true
+        continue
+      end
       α = 1.0
       while (bk < bk_max) && (ft > f + β * α * slope)
         bk = bk + 1
@@ -123,9 +123,21 @@ function trunk(nlp :: AbstractNLPModel;
       BLAS.scal!(n, α, s, 1)
       slope *= α
       Δq = slope + 0.5 * α * α * curv
-      ratio!(tr, nlp, f, ft, Δq, xt, s, slope)
+      ared, pred = aredpred(nlp, f, ft, Δq, xt, s, slope)
+      if pred ≥ 0
+        status = "nonnegative predicted reduction"
+        stalled = true
+        continue
+      end
+      tr.ratio = ared / pred
       if !monotone
-        ρ_hist = ratio(nlp, fref, ft, σref + Δq, xt, s, slope)
+        ared_hist, pred_hist = aredpred(nlp, fref, ft, σref + Δq, xt, s, slope)
+        if pred_hist ≥ 0
+          status = "nonnegative predicted reduction"
+          stalled = true
+          continue
+        end
+        ρ_hist = ared_hist / pred_hist
         set_property!(tr, :ratio, max(get_property(tr, :ratio), ρ_hist))
       end
     end
@@ -149,7 +161,7 @@ function trunk(nlp :: AbstractNLPModel;
             σcur = 0.0
           end
 
-          if nm_iter >= nm_itmax
+          if nm_iter ≥ nm_itmax
             fref = fcur
             σref = σcur
           end
@@ -172,7 +184,7 @@ function trunk(nlp :: AbstractNLPModel;
 
     infoline = @sprintf("%4d  %9.2e  %7.1e  %7.1e  ", iter, f, ∇fNorm2, get_property(tr, :radius))
 
-    optimal = ∇fNorm2 <= ϵ
+    optimal = ∇fNorm2 ≤ ϵ
     tired = nlp.counters.neval_obj > max_f
   end
   @info(trunklogger, infoline)

diff --git a/src/trust-region/tron-trust-region.jl b/src/trust-region/tron-trust-region.jl
@@ -35,13 +35,13 @@ type TRONTrustRegion <: AbstractTrustRegion
   end
 end
 
-function ratio!(tr :: TRONTrustRegion, nlp :: AbstractNLPModel, f :: Float64,
-                f_trial :: Float64, Δm :: Float64, x_trial :: Vector{Float64},
-                step :: Vector{Float64}, slope :: Float64)
-  tr.ratio = ratio(nlp, f, f_trial, Δm, x_trial, step, slope)
+function aredpred(tr :: TRONTrustRegion, nlp :: AbstractNLPModel, f :: Float64,
+                  f_trial :: Float64, Δm :: Float64, x_trial :: Vector{Float64},
+                  step :: Vector{Float64}, slope :: Float64)
+  ared, pred = aredpred(nlp, f, f_trial, Δm, x_trial, step, slope)
   γ = f_trial - f - slope
-  tr.quad_min = γ <= 0.0 ? tr.increase_factor : max(tr.large_decrease_factor, -0.5 * slope / γ)
-  return tr
+  quad_min = γ <= 0.0 ? tr.increase_factor : max(tr.large_decrease_factor, -0.5 * slope / γ)
+  return ared, pred, quad_min
 end
 
 function update!(tr :: TRONTrustRegion, step_norm :: Float64)

diff --git a/src/trust-region/trust-region.jl b/src/trust-region/trust-region.jl
@@ -24,18 +24,17 @@ and the following function:
 """
 @compat abstract type AbstractTrustRegion end
 
-"""`ρ = ratio(nlp, f, f_trial, Δm, x_trial, step, slope)`
+"""`ared, pred = aredpred(nlp, f, f_trial, Δm, x_trial, step, slope)`
 
-Compute the actual vs. predicted reduction radio `∆f/Δm`, where
+Compute the actual and predicted reductions `∆f` and `Δm`, where
 `Δf = f_trial - f` is the actual reduction is an objective/merit/penalty function,
 `Δm = m_trial - m` is the reduction predicted by the model `m` of `f`.
 We assume that `m` is being minimized, and therefore that `Δm < 0`.
 """
-function ratio(nlp :: AbstractNLPModel, f :: Float64, f_trial :: Float64, Δm :: Float64,
-               x_trial :: Vector{Float64}, step :: Vector{Float64}, slope :: Float64)
+function aredpred(nlp :: AbstractNLPModel, f :: Float64, f_trial :: Float64, Δm :: Float64,
+                  x_trial :: Vector{Float64}, step :: Vector{Float64}, slope :: Float64)
   absf = abs(f)
   pred = Δm - max(1.0, absf) * 10.0 * ϵ
-  pred < 0 || throw(TrustRegionException(@sprintf("Nonnegative predicted reduction: pred = %8.1e", pred)))
 
   ared = f_trial - f + max(1.0, absf) * 10.0 * ϵ
   if (abs(Δm) < 1.0e+4 * ϵ) || (abs(ared) < 1.0e+4 * ϵ * absf)
@@ -44,23 +43,9 @@ function ratio(nlp :: AbstractNLPModel, f :: Float64, f_trial :: Float64, Δm ::
     slope_trial = dot(g_trial, step)
     ared = (slope_trial + slope) / 2
   end
-  return ared / pred
+  return ared, pred
 end
 
-"""`tr = ratio!(tr, nlp, f, f_trial, Δm, x_trial, step, slope)`
-
-Compute the actual vs. predicted reduction radio `∆f/Δm`, where
-`Δf = f_trial - f` is the actual reduction is an objective/merit/penalty function,
-`Δm = m_trial - m` is the reduction predicted by the model `m` of `f`.
-We assume that `m` is being minimized, and therefore that `Δm < 0`.
-`tr` is updated inplace.
-"""
-function ratio!(tr :: AbstractTrustRegion, nlp :: AbstractNLPModel, f :: Float64,
-                f_trial :: Float64, Δm :: Float64, x_trial :: Vector{Float64},
-                step :: Vector{Float64}, slope :: Float64)
-  tr.ratio = ratio(nlp, f, f_trial, Δm, x_trial, step, slope)
-  return tr
-end
 
 """`acceptable(tr)`