Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Lower bound improvement based on sharpness of the objective #195

Merged
merged 52 commits into from
Sep 20, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
52 commits
Select commit Hold shift + click to select a range
367e73f
Move tightening functions into a separate file.
Aug 30, 2024
4b55f60
Tighten bounds using sharpness.
Aug 30, 2024
7f9b4c7
Add sharpness constant and exponent to the interface.
Sep 2, 2024
5819c74
Minor fixes.
Sep 2, 2024
81000b6
Tests for strong convexity and sharpness on the Normbox example.
Sep 2, 2024
df13669
Use the node dual gap instead of root dual gap.
Sep 2, 2024
a17f993
Make sure that we collect the statistics, even if in case of timeout.
Sep 3, 2024
06a0e25
We only need this if the key is not in the result dictionary.
Sep 3, 2024
c8add28
Move the BnB Callback into the appriopiate file.
Sep 3, 2024
0fad6d2
We don't need the root dual gap.
Sep 3, 2024
48ccb45
Update src/custom_bonobo.jl
dhendryc Sep 3, 2024
5650287
Add a bit of give to the sharpness constant to avoid cutting of optim…
Sep 4, 2024
f1b676f
Merge changes from repository
Sep 4, 2024
26b52a0
Added small definition of sharpness.
Sep 6, 2024
6cd74a7
Update src/callbacks.jl
matbesancon Sep 10, 2024
7aab85b
Update src/callbacks.jl
matbesancon Sep 10, 2024
32627ee
Don't create child node if the domain oracle is infeasible.
Sep 11, 2024
414a3f9
Set version up after bug fix.
Sep 11, 2024
056e088
Merge branch 'main' into sharpness
Sep 12, 2024
d1b0441
Merge changes from the remote repository.
Sep 12, 2024
53b2c52
Merge branch 'domain-branch' into sharpness
Sep 12, 2024
d79263c
Optimal Experiment Design example.
Sep 12, 2024
802f3b8
Fix syntax issues.
Sep 12, 2024
c993562
Apply suggestions from code review
matbesancon Sep 13, 2024
a068fc8
Add corrected sharpness constants.
Sep 13, 2024
b1b5273
In case the vertex is domain infeasible, use x to get the data type o…
Sep 13, 2024
b676ef6
Syntax issue fix.
Sep 13, 2024
1f26759
Corrections.
Sep 13, 2024
98783d5
Print statement in case assert fails.
Sep 13, 2024
8dd9ea4
Merge changes from github repository
Sep 13, 2024
c55223f
Print logs for better bug hunting.
Sep 16, 2024
b5e7e24
Another try to reproduce error.
Sep 16, 2024
722cb2d
Test lower bound against tree.incumbent plus current fw_epsilon.
Sep 16, 2024
c6aaf84
Have the same seed as in runtests.jl.
Sep 16, 2024
dc65962
Rename A to Ex_mat to avoid constant redefinition.
Sep 16, 2024
4c16b12
Clean up another constant redefinition.
Sep 16, 2024
3a4da67
Add test wrapper.
Sep 16, 2024
015b4a7
Debug statements.
Sep 17, 2024
9030728
More debug statements.
Sep 17, 2024
ee85ecd
Merge branch 'domain-branch' into sharpness
Sep 17, 2024
1cc7254
If heuristic return infeasible solution, ignore it.
Sep 17, 2024
2936851
Minor change.
Sep 17, 2024
961c1a7
Clean up example.
Sep 17, 2024
e1d298b
Show solutions.
Sep 17, 2024
4e42b57
Don't use vertex storage for the heuristics.
Sep 17, 2024
9837858
Final clean up.
Sep 18, 2024
c0c822f
Run tests with one specific seed.
Sep 18, 2024
f0dedc0
Merge branch 'main' into sharpness
Sep 18, 2024
7d93f84
Added time limit.
Sep 18, 2024
c5b0019
Corrected sharpness constant.
Sep 19, 2024
069c4d8
Corrected sharpness constant and slight relaxation of test for A-Opti…
Sep 19, 2024
db20f7e
Update src/tightenings.jl
matbesancon Sep 20, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 0 additions & 6 deletions examples/birkhoff.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 0 additions & 6 deletions examples/int_sparse_reg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 0 additions & 4 deletions examples/mps-example.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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"))

Expand Down
4 changes: 0 additions & 4 deletions examples/mps-examples/mip-examples.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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"

Expand Down
6 changes: 1 addition & 5 deletions examples/nonlinear.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
293 changes: 293 additions & 0 deletions examples/oed_utils.jl
Original file line number Diff line number Diff line change
@@ -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
Loading
Loading