Merge pull request #612 from chriscoey/newconvconditions

use new convergence conditions
jump-dev · Oct 23, 2020 · a3b123e · a3b123e
2 parents bb9e774 + 7095c0d
commit a3b123e
Show file tree

Hide file tree

Showing 8 changed files with 90 additions and 83 deletions.
diff --git a/examples/runbenchmarks.jl b/examples/runbenchmarks.jl
@@ -37,9 +37,10 @@ hyp_solver = ("Hypatia", Hypatia.Optimizer, (
     verbose = true,
     iter_limit = 250,
     time_limit = solve_time_limit,
-    tol_abs_opt = tol,
+    tol_abs_opt = 1e-3 * tol,
     tol_rel_opt = tol,
     tol_feas = tol,
+    tol_infeas = 1e-3 * tol
     ))
 mosek_solver = ("Mosek", Mosek.Optimizer, (
     QUIET = false,

diff --git a/src/Cones/wsosinterpepinormone.jl b/src/Cones/wsosinterpepinormone.jl
@@ -334,9 +334,13 @@ function update_inv_hess_prod(cone::WSOSInterpEpiNormOne{T}) where {T}
             copyto!(diag_r, cone.hess_diag_blocks[r])
             increase_diag!(diag_r)
             r_fact = hess_diag_facts[r1] = cholesky!(Symmetric(diag_r, :U), check = false)
-            if !isposdef(r_fact) && T <: BlasReal
+            if !isposdef(r_fact)
                 copyto!(diag_r, cone.hess_diag_blocks[r])
-                hess_diag_facts[r1] = bunchkaufman!(Symmetric(diag_r, :U), true)
+                if T <: BlasReal # TODO refac
+                    hess_diag_facts[r1] = bunchkaufman!(Symmetric(diag_r, :U), true)
+                else
+                    hess_diag_facts[r1] = lu!(Symmetric(diag_r, :U))
+                end
             end
         end
 
@@ -357,7 +361,11 @@ function update_inv_hess_prod(cone::WSOSInterpEpiNormOne{T}) where {T}
         s_fact = cone.hess_schur_fact = cholesky!(Symmetric(schur_backup, :U), check = false)
         if !isposdef(s_fact)
             copyto!(schur_backup, schur)
-            cone.hess_schur_fact = bunchkaufman!(Symmetric(schur_backup, :U), true)
+                if T <: BlasReal # TODO refac
+                cone.hess_schur_fact = bunchkaufman!(Symmetric(schur_backup, :U), true)
+            else
+                cone.hess_schur_fact = lu!(Symmetric(schur_backup, :U))
+            end
         end
     end
 

diff --git a/src/Solvers/Solvers.jl b/src/Solvers/Solvers.jl
@@ -59,6 +59,7 @@ mutable struct Solver{T <: Real}
     tol_rel_opt::T
     tol_abs_opt::T
     tol_feas::T
+    tol_infeas::T
     tol_slow::T
     preprocess::Bool
     reduce::Bool
@@ -113,7 +114,6 @@ mutable struct Solver{T <: Real}
     primal_obj::T
     dual_obj::T
     gap::T
-    rel_gap::T
     x_feas::T
     y_feas::T
     z_feas::T
@@ -125,7 +125,6 @@ mutable struct Solver{T <: Real}
     prev_is_slow::Bool
     prev2_is_slow::Bool
     prev_gap::T
-    prev_rel_gap::T
     prev_x_feas::T
     prev_y_feas::T
     prev_z_feas::T
@@ -143,6 +142,8 @@ mutable struct Solver{T <: Real}
         tol_rel_opt::RealOrNothing = nothing,
         tol_abs_opt::RealOrNothing = nothing,
         tol_feas::RealOrNothing = nothing,
+        tol_infeas::RealOrNothing = nothing,
+        default_tol_power::RealOrNothing = nothing,
         default_tol_relax::RealOrNothing = nothing,
         tol_slow::Real = 1e-3,
         preprocess::Bool = true,
@@ -162,18 +163,27 @@ mutable struct Solver{T <: Real}
         end
         @assert !(init_use_indirect && preprocess) # cannot use preprocessing and indirect methods for initial point
 
-        default_tol = (T <: LinearAlgebra.BlasReal ? sqrt(eps(T)) : eps(T) ^ 0.4)
+        if isnothing(default_tol_power)
+            default_tol_power = (T <: LinearAlgebra.BlasReal ? 0.5 : 0.4)
+        end
+        default_tol_power = T(default_tol_power)
+        default_tol_loose = eps(T) ^ default_tol_power
+        default_tol_tight = eps(T) ^ (T(1.5) * default_tol_power)
         if !isnothing(default_tol_relax)
-            default_tol *= default_tol_relax
+            default_tol_loose *= T(default_tol_relax)
+            default_tol_tight *= T(default_tol_relax)
         end
         if isnothing(tol_rel_opt)
-            tol_rel_opt = default_tol
+            tol_rel_opt = default_tol_loose
         end
         if isnothing(tol_abs_opt)
-            tol_abs_opt = default_tol
+            tol_abs_opt = default_tol_tight
         end
         if isnothing(tol_feas)
-            tol_feas = default_tol
+            tol_feas = default_tol_loose
+        end
+        if isnothing(tol_infeas)
+            tol_infeas = default_tol_tight
         end
 
         solver = new{T}()
@@ -184,6 +194,7 @@ mutable struct Solver{T <: Real}
         solver.tol_rel_opt = tol_rel_opt
         solver.tol_abs_opt = tol_abs_opt
         solver.tol_feas = tol_feas
+        solver.tol_infeas = tol_infeas
         solver.tol_slow = tol_slow
         solver.preprocess = preprocess
         solver.reduce = reduce
@@ -218,7 +229,6 @@ function solve(solver::Solver{T}) where {T <: Real}
     solver.primal_obj = NaN
     solver.dual_obj = NaN
     solver.gap = NaN
-    solver.rel_gap = NaN
     solver.x_feas = NaN
     solver.y_feas = NaN
     solver.z_feas = NaN
@@ -258,17 +268,12 @@ function solve(solver::Solver{T}) where {T <: Real}
         solver.y_residual = zero(model.b)
         solver.z_residual = zero(model.h)
 
-        # TODO try other way of computing conv tols
-        # solver.x_conv_tol = inv(1 + norm(model.c, Inf))
-        # solver.y_conv_tol = inv(1 + norm(model.b, Inf))
-        # solver.z_conv_tol = inv(1 + norm(model.h, Inf))
-        solver.x_conv_tol = inv(max(one(T), norm(model.c)))
-        solver.y_conv_tol = inv(max(one(T), norm(model.b)))
-        solver.z_conv_tol = inv(max(one(T), norm(model.h)))
+        solver.x_conv_tol = inv(1 + norm(model.c, Inf))
+        solver.y_conv_tol = inv(1 + norm(model.b, Inf))
+        solver.z_conv_tol = inv(1 + norm(model.h, Inf))
         solver.prev_is_slow = false
         solver.prev2_is_slow = false
         solver.prev_gap = NaN
-        solver.prev_rel_gap = NaN
         solver.prev_x_feas = NaN
         solver.prev_y_feas = NaN
         solver.prev_z_feas = NaN
@@ -324,14 +329,15 @@ function solve(solver::Solver{T}) where {T <: Real}
 end
 
 function calc_mu(solver::Solver{T}) where {T <: Real}
-    solver.mu = (dot(solver.point.z, solver.point.s) + dot(solver.point.tau, solver.point.kap)) / (1 + solver.model.nu)
+    point = solver.point
+    solver.mu = (dot(point.z, point.s) + point.tau[1] * point.kap[1]) / (solver.model.nu + 1)
     return solver.mu
 end
 
 function calc_residual(solver::Solver{T}) where {T <: Real}
     model = solver.model
     point = solver.point
-    tau = solver.point.tau[1]
+    tau = point.tau[1]
 
     # x_residual = -A'*y - G'*z - c*tau
     x_residual = solver.x_residual
@@ -365,7 +371,6 @@ function calc_convergence_params(solver::Solver{T}) where {T <: Real}
     point = solver.point
 
     solver.prev_gap = solver.gap
-    solver.prev_rel_gap = solver.rel_gap
     solver.prev_x_feas = solver.x_feas
     solver.prev_y_feas = solver.y_feas
     solver.prev_z_feas = solver.z_feas
@@ -375,13 +380,6 @@ function calc_convergence_params(solver::Solver{T}) where {T <: Real}
     solver.primal_obj = solver.primal_obj_t / solver.point.tau[1] + model.obj_offset
     solver.dual_obj = solver.dual_obj_t / solver.point.tau[1] + model.obj_offset
     solver.gap = dot(point.z, point.s)
-    if solver.primal_obj < -eps(T)
-        solver.rel_gap = solver.gap / -solver.primal_obj
-    elseif solver.dual_obj > eps(T)
-        solver.rel_gap = solver.gap / solver.dual_obj
-    else
-        solver.rel_gap = NaN
-    end
 
     solver.x_feas = solver.x_norm_res * solver.x_conv_tol
     solver.y_feas = solver.y_norm_res * solver.y_conv_tol
@@ -393,42 +391,45 @@ end
 function check_convergence(solver::Solver{T}) where {T <: Real}
     # check convergence criteria
     # TODO nearly primal or dual infeasible or nearly optimal cases?
-    if max(solver.x_feas, solver.y_feas, solver.z_feas) <= solver.tol_feas &&
-        (solver.gap <= solver.tol_abs_opt || (!isnan(solver.rel_gap) && solver.rel_gap <= solver.tol_rel_opt))
+    tau = solver.point.tau[1]
+    primal_obj_t = solver.primal_obj_t
+    dual_obj_t = solver.dual_obj_t
+
+    is_feas = (max(solver.x_feas, solver.y_feas, solver.z_feas) <= solver.tol_feas)
+    is_abs_opt = (solver.gap <= solver.tol_abs_opt)
+    is_rel_opt = (min(solver.gap / tau, abs(primal_obj_t - dual_obj_t)) <= solver.tol_rel_opt * max(tau, min(abs(primal_obj_t), abs(dual_obj_t))))
+    if is_feas && (is_abs_opt || is_rel_opt)
         solver.verbose && println("optimal solution found; terminating")
         solver.status = Optimal
         return true
     end
-    if solver.dual_obj_t > eps(T)
-        infres_pr = solver.x_norm_res_t * solver.x_conv_tol / solver.dual_obj_t
-        if infres_pr <= solver.tol_feas
-            solver.verbose && println("primal infeasibility detected; terminating")
-            solver.status = PrimalInfeasible
-            solver.primal_obj = solver.primal_obj_t
-            solver.dual_obj = solver.dual_obj_t
-            return true
-        end
+
+    if dual_obj_t > eps(T) && solver.x_norm_res_t <= solver.tol_infeas * dual_obj_t
+        solver.verbose && println("primal infeasibility detected; terminating")
+        solver.status = PrimalInfeasible
+        solver.primal_obj = primal_obj_t
+        solver.dual_obj = dual_obj_t
+        return true
     end
-    if solver.primal_obj_t < -eps(T)
-        infres_du = -max(solver.y_norm_res_t * solver.y_conv_tol, solver.z_norm_res_t * solver.z_conv_tol) / solver.primal_obj_t
-        if infres_du <= solver.tol_feas
-            solver.verbose && println("dual infeasibility detected; terminating")
-            solver.status = DualInfeasible
-            solver.primal_obj = solver.primal_obj_t
-            solver.dual_obj = solver.dual_obj_t
-            return true
-        end
+
+    if primal_obj_t < -eps(T) && max(solver.y_norm_res_t, solver.z_norm_res_t) <= solver.tol_infeas * -primal_obj_t
+        solver.verbose && println("dual infeasibility detected; terminating")
+        solver.status = DualInfeasible
+        solver.primal_obj = primal_obj_t
+        solver.dual_obj = dual_obj_t
+        return true
     end
-    if solver.mu <= solver.tol_feas * T(1e-2) && solver.point.tau[1] <= solver.tol_feas * T(1e-2) * min(one(T), solver.point.kap[1])
+
+    # TODO experiment with ill-posedness check
+    if solver.mu <= solver.tol_infeas && tau <= solver.tol_infeas * min(one(T), solver.point.kap[1])
         solver.verbose && println("ill-posedness detected; terminating")
         solver.status = IllPosed
         return true
     end
 
     if expect_improvement(solver.stepper)
         max_improve = zero(T)
-        for (curr, prev) in ((solver.gap, solver.prev_gap), (solver.rel_gap, solver.prev_rel_gap),
-            (solver.x_feas, solver.prev_x_feas), (solver.y_feas, solver.prev_y_feas), (solver.z_feas, solver.prev_z_feas))
+        for (curr, prev) in ((solver.gap, solver.prev_gap), (solver.x_feas, solver.prev_x_feas), (solver.y_feas, solver.prev_y_feas), (solver.z_feas, solver.prev_z_feas))
             if isnan(prev) || isnan(curr)
                 continue
             end

diff --git a/src/Solvers/steppers/heurcomb.jl b/src/Solvers/steppers/heurcomb.jl
@@ -128,36 +128,36 @@ expect_improvement(::HeurCombStepper) = true
 function print_iteration_stats(stepper::HeurCombStepper{T}, solver::Solver{T}) where {T <: Real}
     if iszero(solver.num_iters)
         if iszero(solver.model.p)
-            @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s %9s %9s %9s\n",
-                "iter", "p_obj", "d_obj", "rel_gap", "abs_gap",
+            @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s %9s %9s\n",
+                "iter", "p_obj", "d_obj", "abs_gap",
                 "x_feas", "z_feas", "tau", "kap", "mu",
                 "gamma", "alpha",
                 )
-            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
-                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.rel_gap, solver.gap,
+            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
+                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.gap,
                 solver.x_feas, solver.z_feas, solver.point.tau[1], solver.point.kap[1], solver.mu
                 )
         else
-            @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s %9s %9s %9s %9s\n",
-                "iter", "p_obj", "d_obj", "rel_gap", "abs_gap",
+            @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s %9s %9s %9s\n",
+                "iter", "p_obj", "d_obj", "abs_gap",
                 "x_feas", "y_feas", "z_feas", "tau", "kap", "mu",
                 "gamma", "alpha",
                 )
-            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
-                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.rel_gap, solver.gap,
+            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
+                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.gap,
                 solver.x_feas, solver.y_feas, solver.z_feas, solver.point.tau[1], solver.point.kap[1], solver.mu
                 )
         end
     else
         if iszero(solver.model.p)
-            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
-                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.rel_gap, solver.gap,
+            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
+                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.gap,
                 solver.x_feas, solver.z_feas, solver.point.tau[1], solver.point.kap[1], solver.mu,
                 stepper.prev_gamma, stepper.prev_alpha,
                 )
         else
-            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
-                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.rel_gap, solver.gap,
+            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
+                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.gap,
                 solver.x_feas, solver.y_feas, solver.z_feas, solver.point.tau[1], solver.point.kap[1], solver.mu,
                 stepper.prev_gamma, stepper.prev_alpha,
                 )

diff --git a/src/Solvers/steppers/predorcorr.jl b/src/Solvers/steppers/predorcorr.jl
@@ -98,37 +98,37 @@ expect_improvement(stepper::PredOrCentStepper) = stepper.prev_is_pred
 function print_iteration_stats(stepper::PredOrCentStepper{T}, solver::Solver{T}) where {T <: Real}
     if iszero(solver.num_iters)
         if iszero(solver.model.p)
-            @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s %9s %5s %9s\n",
-                "iter", "p_obj", "d_obj", "rel_gap", "abs_gap",
+            @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s %5s %9s\n",
+                "iter", "p_obj", "d_obj", "abs_gap",
                 "x_feas", "z_feas", "tau", "kap", "mu",
                 "step", "alpha",
                 )
-            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
-                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.rel_gap, solver.gap,
+            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
+                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.gap,
                 solver.x_feas, solver.z_feas, solver.point.tau[1], solver.point.kap[1], solver.mu
                 )
         else
-            @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s %9s %9s %5s %9s\n",
-                "iter", "p_obj", "d_obj", "rel_gap", "abs_gap",
+            @printf("\n%5s %12s %12s %9s %9s %9s %9s %9s %9s %9s %5s %9s\n",
+                "iter", "p_obj", "d_obj", "abs_gap",
                 "x_feas", "y_feas", "z_feas", "tau", "kap", "mu",
                 "step", "alpha",
                 )
-            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
-                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.rel_gap, solver.gap,
+            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e\n",
+                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.gap,
                 solver.x_feas, solver.y_feas, solver.z_feas, solver.point.tau[1], solver.point.kap[1], solver.mu
                 )
         end
     else
         step = (stepper.prev_is_pred ? "pred" : "cent")
         if iszero(solver.model.p)
-            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %5s %9.2e\n",
-                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.rel_gap, solver.gap,
+            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %5s %9.2e\n",
+                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.gap,
                 solver.x_feas, solver.z_feas, solver.point.tau[1], solver.point.kap[1], solver.mu,
                 step, stepper.prev_alpha,
                 )
         else
-            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %5s %9.2e %9.2e\n",
-                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.rel_gap, solver.gap,
+            @printf("%5d %12.4e %12.4e %9.2e %9.2e %9.2e %9.2e %9.2e %9.2e %5s %9.2e %9.2e\n",
+                solver.num_iters, solver.primal_obj, solver.dual_obj, solver.gap,
                 solver.x_feas, solver.y_feas, solver.z_feas, solver.point.tau[1], solver.point.kap[1], solver.mu,
                 step, stepper.prev_alpha,
                 )

diff --git a/test/runcblibtests.jl b/test/runcblibtests.jl
@@ -26,15 +26,12 @@ instance_sets = [
     ]
 
 # options to solvers
-tol = 1e-7
 default_options = (
     verbose = false,
     # verbose = true,
     iter_limit = 150,
     time_limit = 120,
-    tol_rel_opt = tol,
-    tol_abs_opt = tol,
-    tol_feas = tol,
+    default_tol_relax = 100,
     # system_solver = Solvers.NaiveDenseSystemSolver{Float64}(),
     system_solver = Solvers.SymIndefSparseSystemSolver{Float64}(),
     # system_solver = Solvers.QRCholDenseSystemSolver{Float64}(),

diff --git a/test/runmoitests.jl b/test/runmoitests.jl
@@ -27,7 +27,7 @@ end
 default_options = (
     # verbose = true,
     verbose = false,
-    default_tol_relax = 3,
+    default_tol_relax = 2,
     )
 
 @testset "MOI.Test tests" begin

diff --git a/test/runnativetests.jl b/test/runnativetests.jl
@@ -84,7 +84,7 @@ end
 @testset "indirect solvers tests" begin
     println("\nstarting indirect solvers tests")
     for inst_name in inst_indirect, T in diff_reals
-        options = (; default_options..., init_use_indirect = true, preprocess = false, reduce = false, system_solver = Solvers.SymIndefIndirectSystemSolver{T}(), tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4)
+        options = (; default_options..., init_use_indirect = true, preprocess = false, reduce = false, system_solver = Solvers.SymIndefIndirectSystemSolver{T}(), tol_feas = 1e-4, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4, tol_infeas = 1e-6)
         test_instance_solver(inst_name, T, options)
     end
 end