Lower bound improvement based on sharpness of the objective (#195)

* Move tightening functions into a separate file.

* Tighten bounds using sharpness.

* Add sharpness constant and exponent to the interface.

* Minor fixes.

* Tests for strong convexity and sharpness on the Normbox example.

* Use the node dual gap instead of root dual gap.

* Make sure that we collect the statistics, even if in case of timeout.

* We only need this if the key is not in the result dictionary.

* Move the BnB Callback into the appriopiate file.

* We don't need the root dual gap.

* Update src/custom_bonobo.jl

Co-authored-by: Mathieu Besançon <[email protected]>

* Add a bit of give to the sharpness constant to avoid cutting of optimal solutions.

* Added small definition of sharpness.

* Update src/callbacks.jl

* Update src/callbacks.jl

* Don't create child node if the domain oracle is infeasible.

* Set version up after bug fix.

* Optimal Experiment Design example.

* Fix syntax issues.

* Apply suggestions from code review

* Add corrected sharpness constants.

* In case the vertex is domain infeasible, use x to get the data type of f(x).

* Syntax issue fix.

* Corrections.

* Print statement in case assert fails.

* Print logs for better bug hunting.

* Another try to reproduce error.

* Test lower bound against tree.incumbent plus current fw_epsilon.

* Have the same seed as in runtests.jl.

* Rename A to Ex_mat to avoid constant redefinition.

* Clean up another constant redefinition.

* Add test wrapper.

* Debug statements.

* More debug statements.

* If heuristic return infeasible solution, ignore it.

* Minor change.

* Clean up example.

* Show solutions.

* Don't use vertex storage for the heuristics.

* Final clean up.

* Run tests with one specific seed.

* Added time limit.

* Corrected sharpness constant.

* Corrected sharpness constant and slight relaxation of test for A-Optimal.

* Update src/tightenings.jl

---------

Co-authored-by: Hendrych <[email protected]>
Co-authored-by: Mathieu Besançon <[email protected]>

examples/birkhoff.jl

@@ -12,12 +12,6 @@ import HiGHS
 # https://arxiv.org/pdf/2011.02752.pdf
 # https://www.sciencedirect.com/science/article/pii/S0024379516001257
 
-# For bug hunting:
-seed = rand(UInt64)
-@show seed
-seed = 0x3eb09305cecf69f0
-Random.seed!(seed)
-
 
 # min_{X, θ} 1/2 * || ∑_{i in [k]} θ_i X_i - Xhat ||^2
 # θ ∈ Δ_k (simplex)

examples/int_sparse_reg.jl

@@ -21,12 +21,6 @@ const MOI = MathOptInterface
 # r - controls how often we have to maximal split on a index.
 # k - is the sparsity parameter. We only want a few non zero entries.
 
-# For bug hunting:
-seed = rand(UInt64)
-@show seed
-#seed = 0xeadb922ca734998b  
-Random.seed!(seed)
-
 n = 10
 m = 30
 l = 5

examples/mps-example.jl

@@ -16,10 +16,6 @@ const MOI = MathOptInterface
 # Number of binaries    231
 # Number of constraints 198
 
-seed = rand(UInt64)
-@show seed
-Random.seed!(seed)
-
 src = MOI.FileFormats.Model(filename="22433.mps")
 MOI.read_from_file(src, joinpath(@__DIR__, "mps-examples/mps-files/22433.mps"))
 

examples/mps-examples/mip-examples.jl

@@ -21,10 +21,6 @@ import Ipopt
 # ran14x18-disj-8     https://miplib.zib.de/instance_details_ran14x18-disj-8.html
 # timtab1             https://miplib.zib.de/instance_details_timtab1.html   (Takes LONG!)
 
-seed = rand(UInt64)
-@show seed
-Random.seed!(seed)
-
 # To see debug statements
 #ENV["JULIA_DEBUG"] = "Boscia"
 

examples/nonlinear.jl

@@ -11,10 +11,6 @@ using Dates
 # using SCIP
 # const MOI = MathOptInterface
 
-n = 40
-seed = 10
-
-Random.seed!(seed)
 
 ################################################################
 # alternative implementation of LMO using MOI and SCIP
@@ -113,6 +109,6 @@ heu2 = Boscia.Heuristic(Boscia.rounding_lmo_01_heuristic, 0.8, :lmo_rounding)
 heuristics = [heu, heu2]
 # heuristics = []
 
-x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, print_iter=500, custom_heuristics=heuristics)
+x, _, _ = Boscia.solve(f, grad!, lmo, verbose=true, print_iter=500, custom_heuristics=heuristics, time_limit=300)
 
 @show x

examples/oed_utils.jl

@@ -0,0 +1,293 @@
+"""
+    build_data
+
+seed - for the Random functions.
+m    - number of experiments.
+fusion - boolean deiciding whether we build the fusion or standard problem.
+corr - boolean deciding whether we build the independent or correlated data.   
+"""
+function build_data(m)
+    # set up
+    n = Int(floor(m/10))
+
+    B = rand(m,n)
+    B = B'*B
+    @assert isposdef(B)
+    D = MvNormal(randn(n),B)
+
+    A = rand(D, m)'
+    @assert rank(A) == n 
+
+    N = floor(1.5*n)
+    u = floor(N/3)
+    ub = rand(1.0:u, m)
+
+    return A, n, N, ub 
+end
+
+"""
+Build Probability Simplex BLMO for Boscia
+"""
+function build_blmo(m, N, ub)
+    simplex_lmo = Boscia.ProbabilitySimplexSimpleBLMO(N)
+    blmo = Boscia.ManagedBoundedLMO(simplex_lmo, fill(0.0, m), ub, collect(1:m), m)
+    return blmo
+end
+
+"""
+Build function and the gradient for the A-criterion. 
+There is the option to build a safe version of the objective and gradient. If x does not
+satisfy the domain oracle, infinity is returned. 
+If one builds the unsafe version, a FrankWolfe line search must be chosen that can take 
+a domain oracle as an input like Secant or MonotonicGenericStepSize. 
+"""
+function build_a_criterion(A; μ=1e-4, build_safe=false)
+    m, n = size(A) 
+    a = m
+    domain_oracle = build_domain_oracle(A, n)
+
+    function f_a(x)
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        U = cholesky(X)
+        X_inv = U \ I
+        return LinearAlgebra.tr(X_inv)/a 
+    end
+
+    function grad_a!(storage, x)
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X*X)
+        F = cholesky(X)
+        for i in 1:length(x)
+            storage[i] = LinearAlgebra.tr(- (F \ A[i,:]) * transpose(A[i,:]))/a
+        end
+        return storage 
+    end 
+
+    function f_a_safe(x)
+        if !domain_oracle(x)
+            return Inf
+        end
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        X_inv = LinearAlgebra.inv(X)
+        return LinearAlgebra.tr(X_inv)/a 
+    end
+
+    function grad_a_safe!(storage, x)
+        if !domain_oracle(x)
+            return fill(Inf, length(x))        
+        end
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X*X)
+        F = cholesky(X)
+        for i in 1:length(x)
+            storage[i] = LinearAlgebra.tr(- (F \ A[i,:]) * transpose(A[i,:]))/a
+        end
+        return storage
+    end
+
+    if build_safe
+        return f_a_safe, grad_a_safe!
+    end
+
+    return f_a, grad_a!
+end
+
+"""
+Build function and gradient for the D-criterion.
+There is the option to build a safe version of the objective and gradient. If x does not
+satisfy the domain oracle, infinity is returned. 
+If one builds the unsafe version, a FrankWolfe line search must be chosen that can take 
+a domain oracle as an input like Secant or MonotonicGenericStepSize. 
+"""
+function build_d_criterion(A; μ =0.0, build_safe=false)
+    m, n = size(A)
+    a = 1#m
+    domain_oracle = build_domain_oracle(A, n)
+
+    function f_d(x)
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        return float(-log(det(X))/a)
+    end
+
+    function grad_d!(storage, x)
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        F = cholesky(X) 
+        for i in 1:length(x)        
+            storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) 
+        end
+        return storage
+    end
+
+    function f_d_safe(x)
+        if !domain_oracle(x)
+            return Inf
+        end
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        return float(-log(det(X))/a)
+    end
+
+    function grad_d_safe!(storage, x)
+        if !domain_oracle(x)
+            return fill(Inf, length(x))
+        end
+        X = transpose(A)*diagm(x)*A + Matrix(μ *I, n, n)
+        X = Symmetric(X)
+        F = cholesky(X) 
+        for i in 1:length(x)        
+            storage[i] = 1/a * LinearAlgebra.tr(-(F \ A[i,:] )*transpose(A[i,:])) 
+        end
+        # https://stackoverflow.com/questions/46417005/exclude-elements-of-array-based-on-index-julia
+        return storage
+    end
+
+    if build_safe
+        return f_d_safe, grad_d_safe!
+    end
+
+    return f_d, grad_d!
+end
+
+"""
+Find n linearly independent rows of A to build the starting point.
+"""
+function linearly_independent_rows(A)
+    S = []
+    m, n = size(A)
+    for i in 1:m
+        S_i= vcat(S, i)
+        if rank(A[S_i,:])==length(S_i)
+            S=S_i
+        end
+        if length(S) == n # we only n linearly independent points
+            return S
+        end
+    end 
+    return S # then x= zeros(m) and x[S] = 1
+end
+
+function add_to_min(x, u)
+    perm = sortperm(x)
+    j = findfirst(x->x != 0, x[perm])
+
+    for i in j:length(x)
+        if x[perm[i]] < u[perm[i]]
+            x[perm[i]] += 1
+            break
+        else
+            continue
+        end
+    end
+    return x
+end
+
+function remove_from_max(x)
+    perm = sortperm(x, rev = true)
+    j = findlast(x->x != 0, x[perm])
+
+    for i in 1:j
+        if x[perm[i]] > 1
+            x[perm[i]] -= 1
+            break
+        else
+            continue
+        end
+    end
+    return x
+end
+
+"""
+Build start point used in FrankWolfe and Boscia for the A-Optimal and D-Optimal Design Problem.
+The functions are self concordant and so not every point in the feasible region
+is in the domain of f and grad!.
+""" 
+function build_start_point(A, N, ub)
+    # Get n linearly independent rows of A
+    m, n = size(A)
+    S = linearly_independent_rows(A)
+    @assert length(S) == n
+
+    x = zeros(m)
+    E = []
+    V = Vector{Float64}[]
+
+    while !isempty(setdiff(S, findall(x-> !(iszero(x)),x)))
+        v = zeros(m)
+        while sum(v) < N
+            idx = isempty(setdiff(S, findall(x-> !(iszero(x)),v))) ? rand(setdiff(collect(1:m), S)) : rand(setdiff(S, findall(x-> !(iszero(x)),v)))
+            if !isapprox(v[idx], 0.0)
+                @debug "Index $(idx) already picked"
+                continue
+            end
+            v[idx] = min(ub[idx], N - sum(v))
+            push!(E, idx)
+        end
+        push!(V,v)
+        x = sum(V .* 1/length(V)) 
+    end
+    unique!(V)
+    a = length(V)
+    x = sum(V .* 1/a)
+    active_set= FrankWolfe.ActiveSet(fill(1/a, a), V, x)
+
+    return x, active_set, S
+end
+
+"""
+Create first incumbent for Boscia and custom BB in a greedy fashion.
+"""
+function greedy_incumbent(A, N, ub)
+    # Get n linearly independent rows of A
+    m, n = size(A)
+    S = linearly_independent_rows(A)
+    @assert length(S) == n
+
+    # set entries to their upper bound
+    x = zeros(m)
+    x[S] .= ub[S]
+
+    if isapprox(sum(x), N; atol=1e-4, rtol=1e-2)
+        return x
+    elseif sum(x) > N
+        while sum(x) > N
+            remove_from_max(x)
+        end
+    elseif sum(x) < N
+        S1 = S
+        while sum(x) < N
+            jdx = rand(setdiff(collect(1:m), S1))
+            x[jdx] = min(N-sum(x), ub[jdx])
+            push!(S1,jdx)
+            sort!(S1)
+        end
+    end
+    @assert isapprox(sum(x), N; atol=1e-4, rtol=1e-2)
+    @assert sum(ub - x .>= 0) == m 
+    return x
+end
+
+"""
+Check if given point is in the domain of f, i.e. X = transpose(A) * diagm(x) * A 
+positive definite.
+
+(a) Check the rank of A restricted to the rows activated by x.
+(b) Check if the resulting information matrix A' * diagm(x) * A is psd.
+
+(b) is a bit faster for smaller dimensions (< 100). For larger (> 200) (a) is faster.
+"""
+function build_domain_oracle(A, n)
+    return function domain_oracle(x)
+        S = findall(x-> !iszero(x),x)
+        return length(S) >= n && rank(A[S,:]) == n 
+    end
+end
+
+function build_domain_oracle2(A)
+    return function domain_oracle2(x)
+        return isposdef(Symmetric(A' * diagm(x) * A))
+    end 
+end