Merge pull request #545 from chriscoey/arrowhessprod

WIP: improve oracles for epinorminf and entropy, misc numerical fixes

examples/lyapunovstability/JuMP.jl

@@ -37,19 +37,16 @@ example_tests(::Type{LyapunovStabilityJuMP{Float64}}, ::MinimalInstances) = [
     ((2, 2, false, false),),
     ]
 example_tests(::Type{LyapunovStabilityJuMP{Float64}}, ::FastInstances) = begin
-    options = (tol_feas = 1e-7, tol_rel_opt = 1e-6, tol_abs_opt = 1e-6)
-    relaxed_options = (tol_feas = 1e-3, tol_rel_opt = 1e-4, tol_abs_opt = 1e-4)
+    options = (tol_feas = 1e-6, tol_rel_opt = 1e-6, tol_abs_opt = 1e-6)
     return [
-    ((5, 6, true, true), nothing, relaxed_options),
+    ((5, 6, true, true), nothing, options),
     ((5, 6, true, false), nothing, options),
     ((5, 5, false, true), nothing, options),
     ((5, 5, false, false), nothing, options),
-    # ((10, 20, true, true), nothing, relaxed_options),
+    ((10, 20, true, true), nothing, options),
     ((10, 20, true, false), nothing, options),
     ((15, 15, false, true), nothing, options),
     ((15, 15, false, false), nothing, options),
-    ((25, 30, true, true), nothing, options),
-    ((25, 30, true, false), nothing, options),
     ((30, 30, false, false), nothing, options),
     ]
 end

examples/polymin/JuMP.jl

@@ -128,7 +128,7 @@ function test_extra(inst::PolyMinJuMP{T}, model::JuMP.Model) where T
     @test JuMP.termination_status(model) == MOI.OPTIMAL
     if JuMP.termination_status(model) == MOI.OPTIMAL && !isnan(inst.true_min)
         # check objective value is correct
-        tol = eps(T)^0.2
+        tol = eps(T)^0.1
         @test JuMP.objective_value(model) ≈ inst.true_min atol = tol rtol = tol
     end
 end

examples/polymin/native.jl

@@ -240,7 +240,7 @@ function test_extra(inst::PolyMinNative{T}, result::NamedTuple) where T
     @test result.status == :Optimal
     if result.status == :Optimal && !isnan(inst.true_min)
         # check objective value is correct
-        tol = eps(T)^0.2
+        tol = eps(T)^0.1
         true_min = (inst.use_primal ? -1 : 1) * inst.true_min
         @test result.primal_obj ≈ true_min atol = tol rtol = tol
     end

examples/robustgeomprog/JuMP.jl

@@ -28,10 +28,10 @@ example_tests(::Type{RobustGeomProgJuMP{Float64}}, ::FastInstances) = begin
     options = (tol_feas = 1e-6, tol_rel_opt = 1e-6, tol_abs_opt = 1e-6)
     return [
     ((5, 10), nothing, options),
-    # ((5, 10), ClassicConeOptimizer), nothing, options),
+    ((5, 10), ClassicConeOptimizer), nothing, options),
     ((10, 20), nothing, options),
     ((20, 40), nothing, options),
-    # ((40, 80),),
+    ((40, 80),),
     ]
 end
 example_tests(::Type{RobustGeomProgJuMP{Float64}}, ::SlowInstances) = [

src/Cones/Cones.jl

@@ -71,6 +71,11 @@ inv_hess(cone::Cone) = (cone.inv_hess_updated ? cone.inv_hess : update_inv_hess(
 
 reset_data(cone::Cone) = (cone.feas_updated = cone.grad_updated = cone.hess_updated = cone.inv_hess_updated = cone.hess_fact_updated = false)
 
+function use_sqrt_oracles(cone::Cone)
+    cone.hess_fact_updated || update_hess_fact(cone; recover = true)
+    return (cone.hess_fact_cache isa DensePosDefCache)
+end
+
 use_correction(::Cone) = true
 
 update_hess_prod(cone::Cone) = nothing
@@ -169,36 +174,30 @@ use_heuristic_neighborhood(cone::Cone) = cone.use_heuristic_neighborhood
 function in_neighborhood(cone::Cone{T}, rtmu::T, max_nbhd::T) where {T <: Real}
     # norm(H^(-1/2) * (z + mu * grad))
     nbhd_tmp = cone.nbhd_tmp
+    nbhd_tmp2 = cone.nbhd_tmp2
     g = grad(cone)
     @. nbhd_tmp = cone.dual_point + rtmu * g
 
     if use_heuristic_neighborhood(cone)
         nbhd = norm(nbhd_tmp, Inf) / norm(g, Inf)
         # nbhd = maximum(abs(dj / gj) for (dj, gj) in zip(nbhd_tmp, g)) # TODO try this neighborhood
+    elseif Cones.use_sqrt_oracles(cone)
+        inv_hess_sqrt_prod!(nbhd_tmp2, nbhd_tmp, cone)
+        nbhd = norm(nbhd_tmp2)
     else
-        has_hess_fact_cache = hasfield(typeof(cone), :hess_fact_cache)
-        if has_hess_fact_cache && !update_hess_fact(cone)
+        inv_hess_prod!(nbhd_tmp2, nbhd_tmp, cone)
+        nbhd_sqr = dot(nbhd_tmp2, nbhd_tmp)
+        if nbhd_sqr < -eps(T) # TODO possibly loosen
+            # @warn("numerical failure: cone neighborhood is $nbhd_sqr")
             return false
         end
-        nbhd_tmp2 = cone.nbhd_tmp2
-        if has_hess_fact_cache && cone.hess_fact_cache isa DenseSymCache{T}
-            inv_hess_prod!(nbhd_tmp2, nbhd_tmp, cone)
-            nbhd_sqr = dot(nbhd_tmp2, nbhd_tmp)
-            if nbhd_sqr < -eps(T) # TODO possibly loosen
-                # @warn("numerical failure: cone neighborhood is $nbhd_sqr")
-                return false
-            end
-            nbhd = sqrt(abs(nbhd_sqr))
-        else
-            inv_hess_sqrt_prod!(nbhd_tmp2, nbhd_tmp, cone)
-            nbhd = norm(nbhd_tmp2)
-        end
+        nbhd = sqrt(abs(nbhd_sqr))
     end
 
     return (nbhd < rtmu * max_nbhd)
 end
 
-
+# TODO cleanup / remove if not using
 # function in_neighborhood(cone::Cone{T}, dual_point::AbstractVector{T}, rtmu::T, max_nbhd::T) where {T <: Real}
 #     # norm(H^(-1/2) * (z + mu * grad))
 #     nbhd_tmp = cone.nbhd_tmp
@@ -458,18 +457,21 @@ function sparse_upper_arrow(T::Type{<:Real}, w_dim::Int)
 end
 
 function factor_upper_arrow(uu, uw, ww, nzval)
-    minrt = abs2(eps(uu))
+    minval = sqrt(eps(uu)) # TODO tune
     nzidx = 2
     @inbounds for i in eachindex(ww)
-        wwi = sqrt(max(ww[i], minrt))
+        ww1i = ww[i]
+        ww1i < eps(uu) && return false
+        wwi = sqrt(ww1i)
         uwi = uw[i] / wwi
         uu -= abs2(uwi)
+        uu < minval && return false
         nzval[nzidx] = uwi
         nzval[nzidx + 1] = wwi
         nzidx += 2
     end
-    nzval[1] = sqrt(max(uu, minrt))
-    return nzval
+    nzval[1] = sqrt(uu)
+    return true
 end
 
 function sparse_upper_arrow_block2(T::Type{<:Real}, w_dim::Int)
@@ -491,20 +493,25 @@ function sparse_upper_arrow_block2(T::Type{<:Real}, w_dim::Int)
 end
 
 function factor_upper_arrow_block2(uu, uv, uw, vv, vw, ww, nzval)
-    minrt = abs2(eps(uu))
+    minval = sqrt(eps(uu)) # TODO tune
     nzidx = 1
     @inbounds for i in eachindex(ww)
-        wwi = sqrt(max(ww[i], minrt))
+        ww1i = ww[i]
+        ww1i < eps(uu) && return false
+        wwi = sqrt(ww1i)
         vwi = vw[i] / wwi
         uwi = uw[i] / wwi
-        vvi = sqrt(max(vv[i] - abs2(vwi), minrt))
+        vv2i = vv[i] - abs2(vwi)
+        vv2i < eps(uu) && return false
+        vvi = sqrt(vv2i)
         uvi = (uv[i] - vwi * uwi) / vvi
         uu -= abs2(uwi) + abs2(uvi)
+        uu < minval && return false
         @. nzval[nzidx .+ (1:5)] = (uvi, vvi, uwi, vwi, wwi)
         nzidx += 5
     end
-    nzval[1] = sqrt(max(uu, minrt))
-    return nzval
+    nzval[1] = sqrt(uu)
+    return true
 end
 
 end

src/Cones/epinormeucl.jl

@@ -52,6 +52,8 @@ use_heuristic_neighborhood(cone::EpiNormEucl) = false
 
 reset_data(cone::EpiNormEucl) = (cone.feas_updated = cone.grad_updated = cone.hess_updated = cone.inv_hess_updated = false)
 
+use_sqrt_oracles(cone::EpiNormEucl) = true
+
 # TODO only allocate the fields we use
 function setup_data(cone::EpiNormEucl{T}) where {T <: Real}
     reset_data(cone)