Merge branch 'master' into quantumrelentr

jump-dev · Aug 24, 2020 · 4826b69 · 4826b69
2 parents 2a43527 + 83f72dc
commit 4826b69
Show file tree

Hide file tree

Showing 9 changed files with 184 additions and 131 deletions.
diff --git a/examples/stabilitynumber/JuMP.jl b/examples/stabilitynumber/JuMP.jl
@@ -0,0 +1,54 @@
+#=
+Copyright 2020, Chris Coey, Lea Kapelevich and contributors
+
+strengthening of the theta function towards the stability number of a graph
+
+TODO add sparse PSD formulation
+=#
+
+using SparseArrays
+
+struct StabilityNumber{T <: Real} <: ExampleInstanceJuMP{T}
+    side::Int
+    use_doublynonnegative::Bool
+end
+
+function build(inst::StabilityNumber{T}) where {T <: Float64} # TODO generic reals
+    side = inst.side
+    sparsity = 1 - inv(side)
+    inv_graph = tril!(sprand(Bool, side, side, sparsity) + I)
+    (row_idxs, col_idxs, _) = findnz(inv_graph)
+    diags = (row_idxs .== col_idxs)
+
+    model = JuMP.Model()
+    JuMP.@variable(model, X[1:length(row_idxs)])
+    X_diag = X[diags]
+    JuMP.@objective(model, Max, 2 * sum(X) - sum(X_diag))
+    JuMP.@constraint(model, sum(X_diag) == 1)
+    X_lifted = sparse(row_idxs, col_idxs, X, side, side)
+    X_vec = JuMP.AffExpr[X_lifted[i, j] for i in 1:side for j in 1:i]
+    if inst.use_doublynonnegative
+        cone_dim = length(X_vec)
+        JuMP.@constraint(model, X_vec .* ModelUtilities.vec_to_svec!(ones(cone_dim)) in Hypatia.DoublyNonnegativeTriCone{T}(cone_dim))
+    else
+        JuMP.@constraint(model, X_vec in MOI.PositiveSemidefiniteConeTriangle(side))
+        JuMP.@constraint(model, X[.!(diags)] .>= 0)
+    end
+
+    return model
+end
+
+instances[StabilityNumber]["minimal"] = [
+    ((2, true),),
+    ((2, false),),
+    ]
+instances[StabilityNumber]["fast"] = [
+    ((20, true),),
+    ((20, false),),
+    ((50, true),),
+    ((50, false),),
+    ]
+instances[StabilityNumber]["slow"] = [
+    ((500, true),),
+    ((500, false),),
+    ]
diff --git a/src/Cones/hypoperlog.jl b/src/Cones/hypoperlog.jl
@@ -5,7 +5,7 @@ Copyright 2019, Chris Coey, Lea Kapelevich and contributors
 (u in R, v in R_+, w in R_+^d) : u <= v*sum(log.(w/v))
 
 barrier modified from "Primal-Dual Interior-Point Methods for Domain-Driven Formulations" by Karimi & Tuncel, 2019
--log(sum_i v*log(w_i/v) - u) - sum_i log(w_i) - d*log(v)
+-log(sum_i v*log(w_i/v) - u) - sum_i log(w_i) - log(v)
 =#
 
 mutable struct HypoPerLog{T <: Real} <: Cone{T}
@@ -71,7 +71,7 @@ function setup_data(cone::HypoPerLog{T}) where {T <: Real}
     return
 end
 
-get_nu(cone::HypoPerLog) = 1 + 2 * (cone.dim - 2)
+get_nu(cone::HypoPerLog) = cone.dim
 
 function set_initial_point(arr::AbstractVector, cone::HypoPerLog)
     (arr[1], arr[2], w) = get_central_ray_hypoperlog(cone.dim - 2)
@@ -117,7 +117,7 @@ function update_grad(cone::HypoPerLog)
 
     g[1] = inv(cone.vlwvu)
     cone.lvwnivlwvu = (d - cone.lwv) / cone.vlwvu
-    g[2] = cone.lvwnivlwvu - d / v
+    g[2] = cone.lvwnivlwvu - inv(v)
     gden = -1 - v / cone.vlwvu
     @. g[3:end] = gden / w
 
@@ -141,7 +141,7 @@ function update_hess(cone::HypoPerLog)
     H[1, 1] = abs2(g[1])
     H[1, 2] = lvwnivlwvu / cone.vlwvu
     @. H[1, 3:end] = -tmpw / cone.vlwvu
-    H[2, 2] = abs2(lvwnivlwvu) + d * (g[1] + inv(v)) / v
+    H[2, 2] = abs2(lvwnivlwvu) + (d * g[1] + inv(v)) / v
     hden = (-v * lvwnivlwvu - 1) / cone.vlwvu
     @. H[2, 3:end] = hden / w
     @inbounds for j in 1:d
@@ -165,21 +165,22 @@ function update_inv_hess(cone::HypoPerLog)
     H = cone.inv_hess.data
     lwv = cone.lwv
     vlwvu = cone.vlwvu
-    const1 = vlwvu / v
-    const2 = vlwvu + v
-    denom = (const1 + 2) * d / v
+    const1 = vlwvu + v
+    denom = vlwvu + (d + 1) * v
+    vw = cone.tmpw
 
-    H[1, 2] = v * lwv * (lwv - d + 1) - u * lwv + d * u
-    @. @views H[1, 3:end] = w * (d * const1 + lwv)
+    H[1, 2] = v * (v * abs2(lwv) - (u + (d - 1) * v) * lwv + d * u) * v
+    @. @views H[1, 3:end] = w * v * (v * lwv + vlwvu) # = (2 * v * lwv - u)
     @views H[1, 1] = (denom - dot(H[1, 2:end], cone.hess[1, 2:end])) / cone.hess[1, 1] # TODO complicated but doable
-    H[2, 2] = const2
-    @. @views H[2, 3:end] = w
-    @views mul!(H[3:end, 3:end], w, w')
+    H[2, 2] = v * const1 * v
+    @. vw = w * v
+    @. @views H[2, 3:end] = v * vw
+    @views mul!(H[3:end, 3:end], vw, vw')
     H ./= denom
     for j in 3:cone.dim
         H[j, j] += abs2(w[j - 2]) * vlwvu
     end
-    @views H[3:end, :] ./= const2
+    @views H[3:end, :] ./= const1
 
     cone.inv_hess_updated = true
     return cone.inv_hess
@@ -219,7 +220,7 @@ function correction(cone::HypoPerLog{T}, primal_dir::AbstractVector{T}) where {T
     corr[1] = (const9 - sumwdw * v_dir + w_dim * vd2 / v) / z / z
 
     corr[2] = ((-dzdv * const9 - w_dim * v_dir * (const5 + const4) / v + sumwdw * (v_dir * const7 - u_dir)) / z + const3) / z -
-        abs2(v_dir / v) * w_dim * (inv(z) + 2 / v)
+        abs2(v_dir / v) * (w_dim * inv(z) + 2 / v)
 
     @. w_corr = const11 .+ const12 * wdw
     w_corr .*= wdw
@@ -231,7 +232,7 @@ function correction(cone::HypoPerLog{T}, primal_dir::AbstractVector{T}) where {T
     return corr
 end
 
-# see analysis in https://github.com/lkapelevich/HypatiaBenchmarks.jl/tree/master/centralpoints
+# see analysis in https://github.com/lkapelevich/HypatiaSupplements.jl/tree/master/centralpoints
 function get_central_ray_hypoperlog(w_dim::Int)
     if w_dim <= 10
         # lookup points where x = f'(x)
@@ -240,28 +241,28 @@ function get_central_ray_hypoperlog(w_dim::Int)
     # use nonlinear fit for higher dimensions
     x = inv(w_dim)
     if w_dim <= 70
-        u = -1.974777 * x ^ 2 + 0.413520 * x - 0.706751
-        v = -1.213389 * x + 1.413551
-        w = -0.406380 * x + 1.411894
+        u = 4.657876 * x ^ 2 - 3.116192 * x + 0.000647
+        v = 0.424682 * x + 0.553392
+        w = 0.760412 * x + 1.001795
     else
-        u = 0.405290 * x - 0.707011
-        v = -1.238597 * x + 1.414216
-        w = -0.511055 * x + 1.414163
+        u = -3.011166 * x - 0.000122
+        v = 0.395308 * x + 0.553955
+        w = 0.837545 * x + 1.000024
     end
     return [u, v, w]
 end
 
 const central_rays_hypoperlog = [
-    -0.827838399  0.805102005  1.290927713;
-    -0.751337431  0.980713381  1.317894791;
-    -0.716423551  1.079796942  1.331762729;
-    -0.699644766  1.144036715 1.341797042;
-    -0.69134357  1.188706149  1.349742329;
-    -0.687251501  1.221310686  1.3562255;
-    -0.685353717  1.246016352  1.361602711;
-    -0.684641818  1.265307905  1.366119586;
-    -0.684585293  1.280747581  1.369956554;
-    -0.684893372  1.293360445  1.373249434;
+    -0.827838387  0.805102007  1.290927686
+    -0.689607388  0.724605082  1.224617936
+    -0.584372665  0.68128058  1.182421942
+    -0.503499342  0.65448622  1.153053152
+    -0.440285893  0.636444224  1.131466926
+    -0.389979809  0.623569352  1.114979519
+    -0.349255921  0.613978276  1.102013921
+    -0.315769104  0.606589839  1.091577908
+    -0.287837744  0.600745284  1.083013
+    -0.264242734  0.596019009  1.075868782
     ]
 
 # TODO add hess prod, inv hess etc functions

diff --git a/src/Cones/hypoperlogdettri.jl b/src/Cones/hypoperlogdettri.jl
@@ -5,13 +5,11 @@ Copyright 2018, Chris Coey, Lea Kapelevich and contributors
 (u in R, v in R_+, w in S_+) : u <= v*logdet(W/v)
 (see equivalent MathOptInterface LogDetConeTriangle definition)
 
-barrier (self-concordance follows from theorem 5.1.4, "Interior-Point Polynomial Algorithms in Convex Programming" by Y. Nesterov and A. Nemirovski):
-theta^2 * (-log(v*logdet(W/v) - u) - logdet(W) - (n + 1) log(v))
-we use theta = 16
+barrier analogous to hypoperlog cone
+-log(v*logdet(W/v) - u) - logdet(W) - log(v)
 
 TODO
 - describe complex case
-- try to tune theta parameter
 - investigate numerics in inverse Hessian oracles
 =#
 
@@ -25,7 +23,6 @@ mutable struct HypoPerLogdetTri{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
     point::Vector{T}
     dual_point::Vector{T}
     rt2::T
-    sc_const::T
     timer::TimerOutput
 
     feas_updated::Bool
@@ -62,8 +59,6 @@ mutable struct HypoPerLogdetTri{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
     function HypoPerLogdetTri{T, R}(
         dim::Int;
         use_dual::Bool = false,
-        # sc_const::Real = 256, # TODO reduce this
-        sc_const::Real = 25 / T(9), # NOTE not SC but works well (same as rootdet)
         use_heuristic_neighborhood::Bool = default_use_heuristic_neighborhood(),
         max_neighborhood::Real = default_max_neighborhood(),
         hess_fact_cache = hessian_cache(T),
@@ -85,7 +80,6 @@ mutable struct HypoPerLogdetTri{T <: Real, R <: RealOrComplex{T}} <: Cone{T}
             cone.is_complex = false
         end
         cone.side = side
-        cone.sc_const = sc_const
         cone.hess_fact_cache = hess_fact_cache
         return cone
     end
@@ -117,12 +111,12 @@ function setup_data(cone::HypoPerLogdetTri{T, R}) where {R <: RealOrComplex{T}}
     return
 end
 
-get_nu(cone::HypoPerLogdetTri) = 2 * cone.sc_const * (cone.side + 1)
+get_nu(cone::HypoPerLogdetTri) = cone.side + 2
 
 function set_initial_point(arr::AbstractVector{T}, cone::HypoPerLogdetTri{T, R}) where {R <: RealOrComplex{T}} where {T <: Real}
     arr .= 0
-    # NOTE if not using theta = 16, rescaling the ray yields central ray
-    (arr[1], arr[2], w) = (sqrt(cone.sc_const) / T(16)) * get_central_ray_hypoperlogdettri(cone.side)
+    # central point data are the same as for hypoperlog
+    (arr[1], arr[2], w) = get_central_ray_hypoperlog(cone.side)
     incr = (cone.is_complex ? 2 : 1)
     k = 3
     @inbounds for i in 1:cone.side
@@ -178,11 +172,11 @@ function update_grad(cone::HypoPerLogdetTri)
     cone.nLz = (cone.side - cone.ldWv) / z
     cone.ldWvuv = cone.ldWv - u / v
     cone.vzip1 = 1 + v / z
-    cone.grad[1] = cone.sc_const / z
-    cone.grad[2] = cone.sc_const * (cone.nLz - (cone.side + 1) / v)
+    cone.grad[1] = inv(z)
+    cone.grad[2] = cone.nLz - inv(v)
     gend = view(cone.grad, 3:cone.dim)
     smat_to_svec!(gend, cone.Wi, cone.rt2)
-    gend .*= -cone.vzip1 * cone.sc_const
+    gend .*= -cone.vzip1
 
     cone.grad_updated = true
     return cone.grad
@@ -201,7 +195,6 @@ function update_hess(cone::HypoPerLogdetTri)
     Hww .*= cone.vzip1
     Wivzi_vec = smat_to_svec!(cone.tmpw, Wivzi, cone.rt2)
     mul!(Hww, Wivzi_vec, Wivzi_vec', true, true)
-    @. @views cone.hess.data[3:cone.dim, :] *= cone.sc_const
 
     cone.hess_updated = true
     return cone.hess
@@ -217,15 +210,14 @@ end
 #     v = cone.point[2]
 #     @views w = cone.point[3:end]
 #     W = Hermitian(svec_to_smat!(cone.mat2, w, cone.rt2), :U)
-#
 #     update_inv_hess_prod(cone)
-#     denom = (2 * side + 1) * v + (side + 1) * z
+#     denom = z + (side + 1) * v
+#
 #     @views Hww = H[3:cone.dim, 3:cone.dim]
 #     symm_kron(Hww, W, cone.rt2)
 #     Hww .*= z
 #     mul!(Hww, w, w', abs2(v) / denom, true)
 #     Hww ./= (z + v)
-#     Hww ./= cone.sc_const
 #
 #     cone.inv_hess_updated = true
 #     return cone.inv_hess
@@ -245,11 +237,10 @@ function update_hess_prod(cone::HypoPerLogdetTri)
     h1end = view(H, 1, 3:cone.dim)
     smat_to_svec!(h1end, cone.Wivzi, cone.rt2)
     h1end ./= -z
-    H[2, 2] = abs2(cone.nLz) + (cone.side / z + inv(v) * (cone.side + 1)) / v
+    H[2, 2] = abs2(cone.nLz) + (cone.side / z + inv(v)) / v
     h2end = view(H, 2, 3:cone.dim)
     smat_to_svec!(h2end, cone.Wi, cone.rt2)
     h2end .*= ((cone.ldWv - cone.side) / cone.ldWvuv - 1) / z
-    @. @views cone.hess.data[1:2, :] *= cone.sc_const
 
     cone.hess_prod_updated = true
     return
@@ -267,15 +258,14 @@ end
 #     u = cone.point[1]
 #     v = cone.point[2]
 #     w = cone.point[3:end]
-#     const1 = (side + 1) * z / v
-#     denom = (const1 + 2 * side + 1) / v
+#     denom = z + (side + 1) * v
 #
-#     H[1, 2] = ldWv * (v * (ldWv - side + 1) - u) + side * u
-#     @. @views H[1, 3:end] = w * (const1 + ldWv)
-#     @views H[1, 1] = (denom * cone.sc_const - dot(H[1, 2:end], cone.hess[1, 2:end])) / cone.hess[1, 1] # TODO complicated but doable
-#     H[2, 2] = z + v
-#     @views H[2, 3:end] .= w
-#     @views H[1:2, :] ./= cone.sc_const * denom
+#     H[1, 2] = v * (v * abs2(ldWv) - (u + (side - 1) * v) * ldWv + side * u) * v
+#     @. @views H[1, 3:end] = w * v * (v * ldWv + z) # = (2 * v * ldWv - u)
+#     @views H[1, 1] = (denom - dot(H[1, 2:end], cone.hess[1, 2:end])) / cone.hess[1, 1] # TODO complicated but doable
+#     H[2, 2] = v * (z + v) * v
+#     @views H[2, 3:end] .= v * w * v
+#     H ./= denom
 #
 #     cone.inv_hess_prod_updated = true
 #     return
@@ -299,7 +289,7 @@ function hess_prod!(prod::AbstractVecOrMat, arr::AbstractVecOrMat, cone::HypoPer
         ldiv!(cone.fact_mat, cone.mat2)
         smat_to_svec!(view(prod, 3:cone.dim, i), cone.mat2, cone.rt2)
     end
-    @views mul!(prod[3:cone.dim, :], cone.hess[3:cone.dim, 1:2], arr[1:2, :], true, cone.vzip1 * cone.sc_const)
+    @views mul!(prod[3:cone.dim, :], cone.hess[3:cone.dim, 1:2], arr[1:2, :], true, cone.vzip1)
 
     return prod
 end
@@ -314,7 +304,7 @@ end
 #     @views w = cone.point[3:end]
 #     W = Hermitian(svec_to_smat!(cone.mat4, w, cone.rt2), :U)
 #
-#     denom = (2 * side + 1) * v + (side + 1) * z
+#     denom = z + (side + 1) * v
 #     @views mul!(prod[1:2, :], cone.inv_hess[1:2, :], arr)
 #     @inbounds for i in 1:size(arr, 2)
 #         @views arr_w = arr[3:end, i]
@@ -325,7 +315,7 @@ end
 #         mul!(cone.mat2, W, cone.mat3, z, false)
 #         smat_to_svec!(prod_w, cone.mat2, cone.rt2)
 #         prod_w .+= dot(w, arr_w) * abs2(v) .* w / denom
-#         prod_w ./= (z + v) * cone.sc_const
+#         prod_w ./= (z + v)
 #     end
 #     @views mul!(prod[3:cone.dim, :], cone.inv_hess[3:cone.dim, 1:2], arr[1:2, :], true, true)
 #
@@ -380,38 +370,11 @@ function correction(cone::HypoPerLogdetTri{T}, primal_dir::AbstractVector{T}) wh
     const2 = -2 * udz * (pi - side) / z / z
     const3 = (2 * abs2(nLz) + side / z / v) * vdz
     const5 = side / v / z
-    const4 = nLz * (2 * abs2(nLz) + 3 * const5) - (const5 + 2 * (side + 1) / abs2(v)) / v
+    const4 = nLz * (2 * abs2(nLz) + 3 * const5) - (const5 + 2 * inv(v) / v) / v
     corr[1] = 2 * abs2(udz) / z + (2 * const2 + const3) * v_dir + (uuw_scal + const7) * dot_Wi_S + vz * t4awd
     corr[2] = const4 * abs2(v_dir) + (2 * vvw_scal + uvw_scal * u_dir) * dot_Wi_S + u_dir * (const2 + 2 * const3) + const6 * t4awd
 
-    corr *= cone.sc_const / -2
+    corr /= -2
 
     return corr
 end
-
-# see analysis in https://github.com/lkapelevich/HypatiaSupplements.jl/tree/master/centralpoints
-function get_central_ray_hypoperlogdettri(Wside::Int)
-    if Wside <= 10
-        # lookup points where x = f'(x)
-        return central_rays_hypoperlogdettri[Wside, :]
-    end
-    # use nonlinear fit for higher dimensions
-    x = log10(Wside)
-    u = -0.102485 * x ^ 4 + 0.908632 * x ^ 3 - 3.029054 * x ^ 2 + 4.528779 * x - 13.901470
-    v = 0.358933 * x ^ 3 - 2.592002 * x ^ 2 + 6.354740 * x + 17.280377
-    w = 0.027883 * x ^ 3 - 0.231444 * x ^ 2 + 0.652673 * x + 21.997811
-    return [u, v, w]
-end
-
-const central_rays_hypoperlogdettri = [
-    -14.06325335  17.86151855  22.52090275
-    -13.08878205  18.91121795  22.4393585
-    -12.54888342  19.60639116  22.40621157
-    -12.22471372  20.09640151  22.39805249
-    -12.01656536  20.45698931  22.40140061
-    -11.87537532  20.73162694  22.4097267
-    -11.77522327  20.9468238  22.41993593
-    -11.70152722  21.11948341  22.4305642
-    -11.64562635  21.26079849  22.44092946
-    -11.6021318  21.37842775  22.45073131
-    ]