From 0999f87265adec19bdf42ee182bd62029299bbc7 Mon Sep 17 00:00:00 2001
From: Theo Guyard <theo.guyard@insa-rennes.fr>
Date: Thu, 30 Sep 2021 21:54:14 +0200
Subject: [PATCH] Upload code. TODO : complete README
---
.gitignore | 13 +
LICENSE | 21 ++
Project.toml | 5 +
README.md | 72 ++++-
examples/solve_bnb.jl | 39 +++
examples/solve_mip.jl | 27 ++
exp/ICASSP/performances.jl | 171 ++++++++++++
src/BnbScreening.jl | 25 ++
src/activeset.jl | 520 +++++++++++++++++++++++++++++++++++++
src/bnb.jl | 191 ++++++++++++++
src/bvls.jl | 132 ++++++++++
src/data.jl | 87 +++++++
src/displays.jl | 53 ++++
src/mips.jl | 117 +++++++++
src/screen.jl | 52 ++++
src/types.jl | 295 +++++++++++++++++++++
tests/test.jl | 31 +++
17 files changed, 1850 insertions(+), 1 deletion(-)
create mode 100644 .gitignore
create mode 100644 LICENSE
create mode 100644 Project.toml
create mode 100644 examples/solve_bnb.jl
create mode 100644 examples/solve_mip.jl
create mode 100644 exp/ICASSP/performances.jl
create mode 100644 src/BnbScreening.jl
create mode 100644 src/activeset.jl
create mode 100644 src/bnb.jl
create mode 100644 src/bvls.jl
create mode 100644 src/data.jl
create mode 100644 src/displays.jl
create mode 100644 src/mips.jl
create mode 100644 src/screen.jl
create mode 100644 src/types.jl
create mode 100644 tests/test.jl
diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..06fd450
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,13 @@
+.vscode
+.DS_Store
+
+exp/ICASSP/*.txt
+
+# Julia's .gitignore
+Manifest.toml
+*.jl.cov
+*.jl.*.cov
+*.jl.mem
+/docs/build/
+/docs/Manifest.toml
+
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000..0b6a1b6
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 Guyard Theo
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
diff --git a/Project.toml b/Project.toml
new file mode 100644
index 0000000..fcd9fed
--- /dev/null
+++ b/Project.toml
@@ -0,0 +1,5 @@
+[deps]
+CPLEX = "a076750e-1247-5638-91d2-ce28b192dca0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
diff --git a/README.md b/README.md
index 93681b2..8423837 100644
--- a/README.md
+++ b/README.md
@@ -1,2 +1,72 @@
-# BnB screening
+# Branch-screening tests for the L0-penalized least-square problem
+This repository contains numerical procedures linked to the paper
+
+**> ADD CITATION**
+
+If you encounter a bug or something unexpected, please let me know by raising an issue on the project page or by contacting me by email at `theo.guyard@insa-rennes.fr`.
+
+## Requirements
+
+This repository works with Julia v1.5+. Please refer to the [Julia download page](https://julialang.org/downloads/) for install instructions.
+
+Dependencies :
+* [JuMP.jl](https://github.com/jump-dev/JuMP.jl)
+* [CPLEX.jl](https://github.com/jump-dev/CPLEX.jl)
+* [Distributions.jl](https://github.com/JuliaStats/Distributions.jl)
+* [StatsBase](https://github.com/JuliaStats/StatsBase.jl)
+
+Our code works with the CPLEX solver through [CPLEX.jl](https://github.com/jump-dev/CPLEX.jl). You cannot use CPLEX without having purchased and installed a copy of CPLEX Optimization Studio from [IBM](https://www.ibm.com). However, CPLEX is available for free to [academics and students](https://community.ibm.com/community/user/datascience/blogs/xavier-nodet1/2020/07/09/cplex-free-for-students). We recommend to follow carefully the install instructions described at [CPLEX.jl](https://github.com/jump-dev/CPLEX.jl).
+
+## Install instructions
+
+1. Clone the repository
+```bash
+git clone https://gitlab.insa-rennes.fr/Theo.Guyard/bnb-screening.git
+```
+
+2. Enter in the project folder
+```bash
+cd bnb-screening
+```
+
+3. Ensure that Julia v1.5+ is installed
+```bash
+julia -v
+> julia version 1.5.3 (or higher)
+```
+
+4. Install dependencies using [Julia's Pkg](https://docs.julialang.org/en/v1/stdlib/Pkg/) in REPL mode
+```julia
+(@v1.5) pkg> activate .
+(@v1.5) pkg> add JuMP
+(@v1.5) pkg> add CPLEX
+(@v1.5) pkg> add Distributions
+(@v1.5) pkg> add StatsBase
+```
+
+## ICASSP experiments
+
+Experiments presented in the paper corresponding to this repository are located in the folder `bnb-screening/exp/ICASSP/`. The available experiment are
+* `performances.jl` corresponding to table 1
+
+Run experiments directly from the project root folder :
+```bash
+julia --project=. -t 1 exp/ICASSP/performances.jl
+```
+
+The `--project=.` argument installs package dependencies if needed. The flag `-t 1` restricts computations to only 1 CPU core in order to avoid bias due to parallelization capabilities. Logs are saved in the same folder as `performances.jl`. You can edit the experimental setup directly in the script, which is well commented.
+
+## Running demo examples
+
+The `example` folder contains portions of code explaining how to use our package.
+
+## Licence
+
+This software is distributed under the [MIT Licence](https://mit-license.org).
+
+## Cite this work
+
+If you use this package for your own work, please consider citing it as :
+
+**ADD CITATION**
diff --git a/examples/solve_bnb.jl b/examples/solve_bnb.jl
new file mode 100644
index 0000000..d08ade7
--- /dev/null
+++ b/examples/solve_bnb.jl
@@ -0,0 +1,39 @@
+include("../src/BnbScreening.jl")
+
+# ---------------------------------------------------------------------------- #
+# Example : Solving a L0-penalized least-square problem using a BnB algorithm #
+# ---------------------------------------------------------------------------- #
+
+# Generate problem data (gaussian setup or toeplitz setup)
+k, m, n, snr = 5, 500, 300, 10*log10(10)
+A, y, λ, M = data_toeplitz(k,m,n,snr)
+#k, m, n, snr = 5, 500, 1000, 10*log10(10)
+#A, y, λ, M = data_gaussian(k,m,n,snr)
+
+# BnB parameters
+bnbparams = BnbParams(
+ lb_method=ACTIVESET,
+ # Relaxation solving method
+ # - ACTIVESET : the relaxation is solved using an Active-Set algorithm
+ # - LBMIP : the relaxation is solved using CPLEX convex solver
+ ub_method=BVLS,
+ # Heuristic method
+ # - BVLS : the heuristic is solved using the Bounded-Variable least-square algorithm
+ # - LBMIP : the heuristic is solved using CPLEX convex solver
+ maxtime=1000, # Maximum solution time allowed
+ tol=1e-8, # Integer tolerance, ie. |x|<=tol <=> x=0
+ maxiter=1000, # Maximum number of iterations allowed in the Active-Set method
+ screening=false, # Toogle branch-screening rules
+ verbosity=true, # Toogle verbosity
+ showevery=1, # If verbosity=true, number of nodes between two displays
+)
+
+# Solve the problem using the BnB algorithm
+result_bnb = solve_bnb(A,y,λ,M,bnbparams=BnbParams(verbosity=false))
+
+# Display BnB results and statistics
+println("BnB")
+println(" Status : $(result_bnb.termination_status)")
+println(" Objective : $(result_bnb.objective_value)")
+println(" Nodes : $(result_bnb.node_count)")
+println(" Timer : $(result_bnb.solve_time)")
\ No newline at end of file
diff --git a/examples/solve_mip.jl b/examples/solve_mip.jl
new file mode 100644
index 0000000..4630f50
--- /dev/null
+++ b/examples/solve_mip.jl
@@ -0,0 +1,27 @@
+include("../src/BnbScreening.jl")
+
+# ---------------------------------------------------------------------------- #
+# Example : Solving a L0-penalized least-square problem using a MIP solver #
+# ---------------------------------------------------------------------------- #
+
+# Generate problem data (gaussian setup or toeplitz setup)
+k, m, n, snr = 5, 500, 300, 10*log10(10)
+A, y, λ, M = data_toeplitz(k,m,n,snr)
+#k, m, n, snr = 5, 500, 1000, 10*log10(10)
+#A, y, λ, M = data_gaussian(k,m,n,snr)
+
+# MIP solver parameters
+mipparams = MipParams(
+ maxtime = 1000, # Maximum solution time in seconds
+ silent = false, # Toogle solver verbosity
+)
+
+# Solve the problem using a MIP solver
+result_mip = solve_mip(A,y,λ,M)
+
+# Display solver results and statistics
+println("MIP")
+println(" Status : $(result_mip.termination_status)")
+println(" Objective : $(result_mip.objective_value)")
+println(" Nodes : $(result_mip.node_count)")
+println(" Timer : $(result_mip.solve_time)")
\ No newline at end of file
diff --git a/exp/ICASSP/performances.jl b/exp/ICASSP/performances.jl
new file mode 100644
index 0000000..a668323
--- /dev/null
+++ b/exp/ICASSP/performances.jl
@@ -0,0 +1,171 @@
+# ICASSP experient - Performances
+
+include("../../src/BnbScreening.jl")
+using StatsBase
+
+################################ Configuration #################################
+
+# This is the only part of code you may have to edit
+
+setup = "toeplitz" # Setup to use : "gaussian" or "toeplitz"
+k = 7 # Sparsity level
+m, n = 500, 300 # Dictionnary dimension
+snr = 10. * log10(10.) # SNR of the observation noise
+maxtime = 1000 # Maximum solution time (in second) allowed
+maxiter = 1000 # Maximum Active-Set iterations allowed per relaxation resolution
+tol = 1.e-8 # Integer tolerance (ie, |x|<tol <=> x=0)
+repeats = 100 # Number of repeats in the experiment
+
+
+################################## Experiment ##################################
+
+function precompilation(
+ setup::String,
+ k::Int,
+ m::Int,
+ n::Int,
+ snr::Float64,
+ maxtime::Int,
+ maxiter::Int,
+ tol::Float64,
+ )
+ println("Precomiling methods ...")
+ try
+ if setup == "gaussian"
+ A, y, λ, M = data_gaussian(k,m,n,snr)
+ elseif setup == "toeplitz"
+ A, y, λ, M = data_toeplitz(k,m,n,snr)
+ else
+ error("Parameter `setup` must be either \"gaussian\" or \"toeplitz\".")
+ end
+ result_mip = solve_mip(A,y,λ,M,mipparams=MipParams(maxtime=maxtime))
+ result_bnb = solve_bnb(A,y,λ,M,bnbparams=BnbParams(maxtime=maxtime,maxiter=maxiter,tol=tol,verbosity=false))
+ result_bnbs = solve_bnb(A,y,λ,M,bnbparams=BnbParams(maxtime=maxtime,maxiter=maxiter,tol=tol,verbosity=false,screening=true))
+ catch e
+ e == InterruptException() && rethrow(e)
+ println("Precompilation failed with the following error")
+ throw(e)
+ end
+ return nothing
+end
+
+logfile_name = string(
+ @__DIR__,
+ "/",
+ Dates.format(Dates.now(),"yyyy-mm-dd-HH:MM:SS"),
+ ".txt"
+ )
+logfile = open(logfile_name, "w")
+
+println("--------------------------------")
+println("ICASSP experiment - Performances\n")
+println("Configuration")
+println(" setup : $(setup)")
+println(" k : $(k)")
+println(" m, n : $(m), $(n)")
+println(" snr : $(snr)")
+println(" maxtime : $(maxtime)")
+println(" maxiter : $(maxiter)")
+println(" tol : $(tol)")
+println(" repeats : $(repeats)")
+println()
+
+write(logfile,"ICASSP experiment - Performances\n\n")
+write(logfile,"Configuration\n")
+write(logfile," setup : $(setup)\n")
+write(logfile," k : $(k)\n")
+write(logfile," m, n : $(m), $(n)\n")
+write(logfile," snr : $(snr)\n")
+write(logfile," maxtime : $(maxtime)\n")
+write(logfile," maxiter : $(maxiter)\n")
+write(logfile," tol : $(tol)\n")
+write(logfile," repeats : $(repeats)\n")
+
+fail_mip = Vector{Union{Bool,Missing}}(missing,repeats)
+fail_bnb = Vector{Union{Bool,Missing}}(missing,repeats)
+fail_bnbs = Vector{Union{Bool,Missing}}(missing,repeats)
+node_mip = Vector{Union{Int64,Missing}}(missing,repeats)
+node_bnb = Vector{Union{Int64,Missing}}(missing,repeats)
+node_bnbs = Vector{Union{Int64,Missing}}(missing,repeats)
+time_mip = Vector{Union{Float64,Missing}}(missing,repeats)
+time_bnb = Vector{Union{Float64,Missing}}(missing,repeats)
+time_bnbs = Vector{Union{Float64,Missing}}(missing,repeats)
+
+precompilation(setup,k,m,n,snr,maxtime,maxiter,tol)
+
+for i in 1:repeats
+ println("Repeat $i/$repeats")
+ try
+ # Generate data
+ if setup == "gaussian"
+ A, y, λ, M = data_gaussian(k,m,n,snr)
+ elseif setup == "toeplitz"
+ A, y, λ, M = data_toeplitz(k,m,n,snr)
+ end
+
+ # Run the tree methods
+ result_mip = solve_mip(A,y,λ,M,mipparams=MipParams(maxtime=maxtime))
+ result_bnb = solve_bnb(A,y,λ,M,bnbparams=BnbParams(maxtime=maxtime,maxiter=maxiter,tol=tol,verbosity=false))
+ result_bnbs = solve_bnb(A,y,λ,M,bnbparams=BnbParams(maxtime=maxtime,maxiter=maxiter,tol=tol,verbosity=false,screening=true))
+
+ # Gather results
+ if result_mip.termination_status == MOI.OPTIMAL
+ fail_mip[i] = false
+ node_mip[i] = result_mip.node_count
+ time_mip[i] = result_mip.solve_time
+ else
+ fail_mip[i] = true
+ end
+ if result_bnb.termination_status == MOI.OPTIMAL
+ fail_bnb[i] = false
+ node_bnb[i] = result_bnb.node_count
+ time_bnb[i] = result_bnb.solve_time
+ else
+ fail_bnb[i] = true
+ end
+ if result_bnbs.termination_status == MOI.OPTIMAL
+ fail_bnbs[i] = false
+ node_bnbs[i] = result_bnbs.node_count
+ time_bnbs[i] = result_bnbs.solve_time
+ else
+ fail_bnbs[i] = true
+ end
+ catch e
+ e == InterruptException() && rethrow(e)
+ println("warning : Repeat $i failed with the following error")
+ println(e)
+ end
+end
+
+
+################################### Results ####################################
+
+println()
+println("Fails")
+println(" MIP : $(100 * mean(skipmissing(fail_mip)))%")
+println(" BnB : $(100 * mean(skipmissing(fail_bnb)))%")
+println(" SBnB : $(100 * mean(skipmissing(fail_bnbs)))%")
+println("Node count")
+println(" MIP : $(mean(skipmissing(node_mip)))")
+println(" BnB : $(mean(skipmissing(node_bnb)))")
+println(" SBnB : $(mean(skipmissing(node_bnbs)))")
+println("Solve time")
+println(" MIP : $(round(mean(skipmissing(time_mip)),digits=3)) seconds")
+println(" BnB : $(round(mean(skipmissing(time_bnb)),digits=3)) seconds")
+println(" SBnB : $(round(mean(skipmissing(time_bnbs)),digits=3)) seconds")
+println("\n--------------------------------")
+
+write(logfile,"\n")
+write(logfile,"Fails\n")
+write(logfile," MIP : $(100 * mean(skipmissing(fail_mip)))%\n")
+write(logfile," BnB : $(100 * mean(skipmissing(fail_bnb)))%\n")
+write(logfile," SBnB : $(100 * mean(skipmissing(fail_bnbs)))%\n")
+write(logfile,"Node count\n")
+write(logfile," MIP : $(mean(skipmissing(node_mip)))\n")
+write(logfile," BnB : $(mean(skipmissing(node_bnb)))\n")
+write(logfile," SBnB : $(mean(skipmissing(node_bnbs)))\n")
+write(logfile,"Solve time\n")
+write(logfile," MIP : $(round(mean(skipmissing(time_mip)),digits=3)) seconds\n")
+write(logfile," BnB : $(round(mean(skipmissing(time_bnb)),digits=3)) seconds\n")
+write(logfile," SBnB : $(round(mean(skipmissing(time_bnbs)),digits=3)) seconds\n")
+close(logfile)
\ No newline at end of file
diff --git a/src/BnbScreening.jl b/src/BnbScreening.jl
new file mode 100644
index 0000000..6c47a7b
--- /dev/null
+++ b/src/BnbScreening.jl
@@ -0,0 +1,25 @@
+# Standard libraries
+using Dates
+using Printf
+using Random
+
+# External libraries
+using CPLEX
+using Distributions
+using JuMP
+using LinearAlgebra
+
+version() = "v1.0"
+author() = "Theo Guyard"
+contact() = "theo.guyard@insa-rennes.fr"
+license() = "MIT License (https://mit-license.org)"
+
+# Module files
+include("types.jl")
+include("activeset.jl")
+include("bnb.jl")
+include("bvls.jl")
+include("data.jl")
+include("displays.jl")
+include("mips.jl")
+include("screen.jl")
\ No newline at end of file
diff --git a/src/activeset.jl b/src/activeset.jl
new file mode 100644
index 0000000..5365243
--- /dev/null
+++ b/src/activeset.jl
@@ -0,0 +1,520 @@
+tolsign(x::Float64, tol::Float64) = sign(x) * (abs(x) > tol)
+
+function symbi_app(Ainv::Matrix{Float64}, a::Vector{Float64})
+ Ainv_updated = Matrix{Float64}(undef,length(a),length(a))
+ @views begin
+ d = a[1:end-1]
+ g = a[end]
+ dw = BLAS.symv('U',Ainv,-d)
+ gw = 1. / (g + dw' * d)
+ gw_dw = gw .* dw
+ BLAS.syr!('U',gw,dw,Ainv)
+ Ainv_updated[1:end-1,1:end-1] = Ainv
+ Ainv_updated[1:end-1,end] = gw_dw
+ Ainv_updated[end,end] = gw
+ end
+ return Ainv_updated
+end
+
+function symbi_rem(Ainv::Matrix{Float64}, i::Int)
+ n = size(Ainv)[1]
+ Ainv_updated = Matrix{Float64}(undef,n-1,n-1)
+ @views begin
+ q = -1. / Ainv[i,i]
+ c1 = 1:i-1
+ r1 = Ainv[c1,i]
+ c2 = i+1:n
+ r2 = Ainv[i,c2]
+ c2red = i:n-1
+ Ainv_updated[c1,c1] = BLAS.syr!('U',q,r1,Ainv[c1,c1])
+ Ainv_updated[c1,c2red] = BLAS.ger!(q,r1,r2,Ainv[c1,c2])
+ Ainv_updated[c2red,c2red] = BLAS.syr!('U',q,r2,Ainv[c2,c2])
+ end
+ return Ainv_updated
+end
+
+function symbi_ro!(Ainv::Matrix{Float64}, r::Vector{Float64}, g::Float64)
+ Ainv_r = BLAS.symv('U',Ainv,r)
+ c = -g / (1. + g * (r' * Ainv_r))
+ BLAS.syr!('U',c,Ainv_r,Ainv)
+ return nothing
+end
+
+function swap_Sbarin_to_Sbar0!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ i::Int
+ )
+ iSbarin = findfirst(j->j==i,node.support.Sbarin)
+ node.GSbarin_inv = symbi_rem(node.GSbarin_inv,iSbarin)
+ push!(node.support.Sbar0,i)
+ deleteat!(node.support.Sbarin,iSbarin)
+ deleteat!(asviews.xSbarin,iSbarin)
+ deleteat!(asviews.sSbarin,iSbarin)
+ deleteat!(asviews.xnewSbarin,iSbarin)
+ asviews.ASbarin = A[:,node.support.Sbarin]
+ return nothing
+end
+
+function swap_Sbarin_to_Sbarbox!(node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ i::Int
+ )
+ iSbarin = findfirst(j->j==i,node.support.Sbarin)
+ node.GSbarin_inv = symbi_rem(node.GSbarin_inv,iSbarin)
+ push!(node.support.Sbarbox,i)
+ push!(asviews.xSbarbox,popat!(asviews.xSbarin,iSbarin))
+ push!(asviews.sSbarbox,popat!(asviews.sSbarin,iSbarin))
+ deleteat!(asviews.xnewSbarin,iSbarin)
+ deleteat!(node.support.Sbarin,iSbarin)
+ asviews.ASbarin = A[:,node.support.Sbarin]
+ asviews.ASbarbox = A[:,node.support.Sbarbox]
+ return nothing
+end
+
+function swap_S1in_to_S1box!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ B::Matrix{Float64},
+ i::Int
+ )
+ iS1in = findfirst(j->j==i,node.support.S1in)
+ push!(node.support.S1box,i)
+ push!(asviews.xS1box,popat!(asviews.xS1in,iS1in))
+ push!(asviews.sS1box,popat!(asviews.sS1in,iS1in))
+ deleteat!(asviews.xnewS1in,iS1in)
+ deleteat!(node.support.S1in,iS1in)
+ asviews.AS1in = A[:,node.support.S1in]
+ asviews.AS1box = A[:,node.support.S1box]
+
+ @views begin
+ a = A[:,i]
+ d = vcat(node.GS1in_inv[1:(iS1in-1),iS1in], node.GS1in_inv[iS1in,(iS1in+1):end])
+ e = node.GS1in_inv[iS1in,iS1in]
+ Ad = asviews.AS1in * d
+ r1 = (Ad./e) + a
+ BLAS.syr!('U',e,r1,B)
+ node.GS1in_inv = symbi_rem(node.GS1in_inv,iS1in)
+ symbi_ro!(node.GSbarin_inv, asviews.ASbarin'*r1, e)
+ end
+
+ return nothing
+end
+
+function swap_Sbar0_to_Sbarin!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ B::Matrix{Float64},
+ i::Int
+ )
+ @views begin
+ a = A[:,i]
+ Ba = BLAS.symv('U',B,a)
+ ABa = asviews.ASbarin' * Ba
+ aBa = a' * Ba
+ v = vcat(ABa,aBa)
+ node.GSbarin_inv = symbi_app(node.GSbarin_inv,v)
+ end
+ push!(node.support.Sbarin,i)
+ filter!(j->j!=i,node.support.Sbar0)
+ push!(asviews.xSbarin, 0.)
+ push!(asviews.sSbarin, 0.)
+ push!(asviews.xnewSbarin, 0.)
+ asviews.ASbarin = A[:,node.support.Sbarin]
+ return nothing
+end
+
+function swap_Sbarbox_to_Sbarin!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ B::Matrix{Float64},
+ i::Int
+ )
+ iSbarbox = findfirst(j->j==i,node.support.Sbarbox)
+ @views begin
+ a = A[:,i]
+ Ba = BLAS.symv('U',B,a)
+ ABa = asviews.ASbarin' * Ba
+ aBa = a' * Ba
+ v = vcat(ABa,aBa)
+ node.GSbarin_inv = symbi_app(node.GSbarin_inv,v)
+ end
+ push!(node.support.Sbarin, popat!(node.support.Sbarbox,iSbarbox))
+ push!(asviews.xSbarin, popat!(asviews.xSbarbox,iSbarbox))
+ push!(asviews.sSbarin, popat!(asviews.sSbarbox,iSbarbox))
+ push!(asviews.xnewSbarin, asviews.xSbarin[end])
+ asviews.ASbarin = A[:,node.support.Sbarin]
+ asviews.ASbarbox = A[:,node.support.Sbarbox]
+ return nothing
+end
+
+function swap_S1box_to_Sbarin!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ B::Matrix{Float64},
+ i::Int
+ )
+ iS1box = findfirst(j->j==i,node.support.S1box)
+ @views begin
+ a = A[:,i]
+ Ba = BLAS.symv('U',B,a)
+ ABa = asviews.ASbarin' * Ba
+ aBa = a' * Ba
+ v = vcat(ABa,aBa)
+ node.GSbarin_inv = symbi_app(node.GSbarin_inv,v)
+ end
+ push!(node.support.Sbarin, popat!(node.support.S1box,iS1box))
+ push!(asviews.xSbarin, popat!(asviews.xS1box,iS1box))
+ push!(asviews.sSbarin, popat!(asviews.sS1box,iS1box))
+ push!(asviews.xnewSbarin, asviews.xSbarin[end])
+ asviews.ASbarin = A[:,node.support.Sbarin]
+ asviews.AS1box = A[:,node.support.S1box]
+ return nothing
+end
+
+function swap_Sbar_to_S0!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ i::Int
+ )
+ iSbarin = findfirst(j->j==i,node.support.Sbarin)
+ if !isa(iSbarin,Nothing)
+ node.GSbarin_inv = symbi_rem(node.GSbarin_inv,iSbarin)
+ deleteat!(node.support.Sbarin,iSbarin)
+ deleteat!(asviews.xSbarin,iSbarin)
+ deleteat!(asviews.sSbarin,iSbarin)
+ deleteat!(asviews.xnewSbarin,iSbarin)
+ asviews.ASbarin = A[:,node.support.Sbarin]
+ else
+ iSbarbox = findfirst(j->j==i,node.support.Sbarbox)
+ if !isa(iSbarbox,Nothing)
+ deleteat!(node.support.Sbarbox,iSbarbox)
+ deleteat!(asviews.xSbarbox,iSbarbox)
+ deleteat!(asviews.sSbarbox,iSbarbox)
+ asviews.ASbarbox = A[:,node.support.Sbarbox]
+ else
+ filter!(j->j!=i,node.support.Sbar0)
+ end
+ end
+ return nothing
+end
+
+function swap_Sbar_to_S1!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ B::Matrix{Float64},
+ i::Int
+ )
+ iSbarbox = findfirst(j->j==i,node.support.Sbarbox)
+ if !isa(iSbarbox,Nothing)
+ push!(node.support.S1box, popat!(node.support.Sbarbox,iSbarbox))
+ push!(asviews.xS1box, popat!(asviews.xSbarbox,iSbarbox))
+ push!(asviews.sS1box, popat!(asviews.sSbarbox,iSbarbox))
+ asviews.ASbarbox = A[:,node.support.Sbarbox]
+ asviews.AS1box = A[:,node.support.S1box]
+ else
+ @views begin
+ asviews.AS1in = A[:,node.support.S1in]
+ a = A[:,i]
+ Ata = asviews.AS1in' * a
+ ata = a' * a
+ d = - BLAS.symv('U',node.GS1in_inv,Ata)
+ g = -1. / (ata + d'*Ata)
+ r1 = asviews.AS1in * d + a
+ v = vcat(Ata,ata)
+ node.GS1in_inv = symbi_app(node.GS1in_inv,v)
+ BLAS.syr!('U',g,r1,B)
+ end
+ iSbarin = findfirst(j->j==i,node.support.Sbarin)
+ if !isa(iSbarin,Nothing)
+ node.GSbarin_inv = symbi_rem(node.GSbarin_inv,iSbarin)
+ push!(node.support.S1in, popat!(node.support.Sbarin,iSbarin))
+ push!(asviews.xS1in, popat!(asviews.xSbarin,iSbarin))
+ push!(asviews.sS1in, popat!(asviews.sSbarin,iSbarin))
+ push!(asviews.xnewS1in, popat!(asviews.xnewSbarin,iSbarin))
+ asviews.ASbarin = A[:,node.support.Sbarin]
+ asviews.AS1in = A[:,node.support.S1in]
+ symbi_ro!(node.GSbarin_inv, asviews.ASbarin'*r1,g)
+ else
+ iSbar0 = findfirst(j->j==i,node.support.Sbar0)
+ push!(node.support.S1in, popat!(node.support.Sbar0,iSbar0))
+ push!(asviews.xS1in, 0.)
+ push!(asviews.sS1in, 0.)
+ push!(asviews.xnewS1in, 0.)
+ asviews.AS1in = A[:,node.support.S1in]
+ symbi_ro!(node.GSbarin_inv, asviews.ASbarin'*r1,g)
+ end
+ end
+ return nothing
+end
+
+function solve_kkt_system!(
+ node::Node,
+ asviews::ActivesetViews,
+ B::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64
+ )
+ Axbox = asviews.ASbarbox * asviews.xSbarbox + asviews.AS1box * asviews.xS1box
+ r = y - Axbox
+ rhsSbarin = asviews.ASbarin' * BLAS.symv('U',B,r) - (λ/M) .* asviews.sSbarin
+ asviews.xnewSbarin = BLAS.symv('U',node.GSbarin_inv,rhsSbarin)
+ rhsS1in = asviews.AS1in' * (r - asviews.ASbarin * asviews.xnewSbarin)
+ asviews.xnewS1in = BLAS.symv('U',node.GS1in_inv,rhsS1in)
+ return Axbox
+end
+
+function discrete_line_search!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ Axbox::Vector{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ bnbparams::BnbParams,
+ )
+
+ # TODO : optimize
+
+ xdiffSbarin = asviews.xnewSbarin - asviews.xSbarin
+ xdiffS1in = asviews.xnewS1in - asviews.xS1in
+ Ax = asviews.ASbarin * asviews.xSbarin + asviews.AS1in * asviews.xS1in + Axbox
+ Axdiff = asviews.ASbarin * xdiffSbarin + asviews.AS1in * xdiffS1in
+ u = y - Ax
+ λ_M = λ / M
+ λ_S1 = λ * length(node.branch.S1)
+
+ a = Axdiff' * Axdiff
+ b_base = u' * Axdiff
+
+ if a == 0.
+ Atu = A' * u
+ lb = 0.5*u'*u + λ_M*asviews.xSbarin'*asviews.sSbarin + λ_S1
+ return u, Atu, lb
+ end
+
+ tflip = []
+ tmax = 1.
+ for i in 1:length(node.support.Sbarin)
+ if xdiffSbarin[i] < 0.
+ tmax = min(tmax, -(M+asviews.xSbarin[i])/xdiffSbarin[i])
+ else
+ tmax = min(tmax, (M-asviews.xSbarin[i])/xdiffSbarin[i])
+ end
+ push!(tflip, -asviews.xSbarin[i]/xdiffSbarin[i])
+ end
+ for i in 1:length(node.support.S1in)
+ if xdiffS1in[i] < 0.
+ tmax = min(tmax, -(M+asviews.xS1in[i])/xdiffS1in[i])
+ else
+ tmax = min(tmax, (M-asviews.xS1in[i])/xdiffS1in[i])
+ end
+ push!(tflip, -asviews.xS1in[i]/xdiffS1in[i])
+ end
+
+ filter!(ti -> 0. < ti < tmax, sort!(unique!(tflip)))
+ pushfirst!(tflip,0.)
+ push!(tflip,tmax)
+
+ topt = 0.
+ uopt = u
+ vopt = Inf
+
+ for i in 1:(length(tflip)-1)
+
+ tl = tflip[i]
+ tu = tflip[i+1]
+ tmed = (tl+tu)/2.
+
+ thetaSbarin = sign.(asviews.xSbarin + tmed .* xdiffSbarin)
+
+ b = b_base - λ_M * thetaSbarin' * xdiffSbarin
+
+ tinf = max(tl,min(tu,b/a))
+ uinf = u - tinf .* Axdiff
+ vinf = 0.5 * norm(uinf,2)^2 + λ_M * norm(asviews.xSbarin + tinf .* xdiffSbarin,1)
+
+ if vinf < vopt
+ topt = tinf
+ uopt = uinf
+ vopt = vinf
+ end
+ end
+
+ asviews.xSbarin += topt .* xdiffSbarin
+ asviews.xS1in += topt .* xdiffS1in
+ asviews.sSbarin = tolsign.(asviews.xSbarin,bnbparams.tol)
+ asviews.sS1in = tolsign.(asviews.xS1in,bnbparams.tol)
+ Atuopt = A' * uopt
+ lb = vopt + (λ_S1 + λ * sum(node.support.Sbarbox))
+
+ return uopt, Atuopt, lb
+end
+
+function detect_blocking_behaviour(swapped::Vector{Int}, imax::Int)
+ return (length(swapped)==1) & (imax in swapped)
+end
+
+function active_optimality(
+ asviews::ActivesetViews,
+ M::Float64,
+ bnbparams::BnbParams
+ )
+ cond = (
+ all(asviews.sSbarin .== tolsign.(asviews.xnewSbarin,bnbparams.tol)) &
+ all(asviews.sS1in .== tolsign.(asviews.xnewS1in,bnbparams.tol)) &
+ (norm(asviews.xnewSbarin,Inf) - M <= bnbparams.tol) &
+ (norm(asviews.xnewS1in,Inf) - M <= bnbparams.tol)
+ )
+ return cond
+end
+
+function inactive_optimality(
+ node::Node,
+ asviews::ActivesetViews,
+ Atu::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ bnbparams::BnbParams
+ )
+ vmax = 0.
+ imax = nothing
+ smax = ""
+ λ_M = λ / M
+ λ_M_m_tol = λ_M - bnbparams.tol
+ λ_M_p_tol = λ_M + bnbparams.tol
+ m_tol = -bnbparams.tol
+ for i in node.support.Sbar0
+ v = abs(Atu[i]) - λ_M_p_tol
+ if v > vmax
+ vmax = v
+ imax = i
+ smax = "Sbar0"
+ end
+ end
+ for (iSbarbox,i) in enumerate(node.support.Sbarbox)
+ v = λ_M_m_tol - Atu[i] * asviews.sSbarbox[iSbarbox]
+ if v > vmax
+ vmax = v
+ imax = i
+ smax = "Sbarbox"
+ end
+ end
+ for (iS1box,i) in enumerate(node.support.S1box)
+ v = m_tol - Atu[i] * asviews.sS1box[iS1box]
+ if v > vmax
+ vmax = v
+ imax = i
+ smax = "S1box"
+ end
+ end
+ return imax, smax
+end
+
+function new_tentative_sign!(
+ node::Node,
+ asviews::ActivesetViews,
+ imax::Int,
+ smax::String,
+ A::Matrix{Float64},
+ B::Matrix{Float64}
+ )
+ if smax == "Sbar0"
+ swap_Sbar0_to_Sbarin!(node,asviews,A,B,imax)
+ elseif smax == "Sbarbox"
+ swap_Sbarbox_to_Sbarin!(node,asviews,A,B,imax)
+ elseif smax == "S1box"
+ swap_S1box_to_Sbarin!(node,asviews,A,B,imax)
+ end
+ return nothing
+end
+
+function update_sets!(
+ node::Node,
+ asviews::ActivesetViews,
+ A::Matrix{Float64},
+ B::Matrix{Float64},
+ M::Float64,
+ bnbparams::BnbParams
+ )
+
+ swapped = Vector{Int}()
+ for (iS1in,i) in enumerate(node.support.S1in)
+ absxi = abs(asviews.xS1in[iS1in])
+ if abs(M - absxi) <= bnbparams.tol
+ swap_S1in_to_S1box!(node,asviews,A,B,i)
+ push!(swapped,i)
+ end
+ end
+ for (iSbarin,i) in enumerate(node.support.Sbarin)
+ absxi = abs(asviews.xSbarin[iSbarin])
+ if absxi <= bnbparams.tol
+ swap_Sbarin_to_Sbar0!(node,asviews,A,i)
+ push!(swapped,i)
+ elseif abs(M - absxi) <= bnbparams.tol
+ swap_Sbarin_to_Sbarbox!(node,asviews,A,i)
+ push!(swapped,i)
+ end
+ end
+ return swapped
+end
+
+function solve_relax_activeset!(
+ tree::Tree,
+ node::Node,
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ bnbparams::BnbParams,
+ )
+
+ asviews = ActivesetViews(node,A,bnbparams)
+ B = I(size(A)[1]) - BLAS.symm('R','U',node.GS1in_inv,asviews.AS1in) * asviews.AS1in'
+ node.branch_on[2] == 0 && swap_Sbar_to_S0!(node,asviews,A,node.branch_on[1])
+ node.branch_on[2] == 1 && swap_Sbar_to_S1!(node,asviews,A,B,node.branch_on[1])
+
+ it = 0
+ lb = Inf
+
+ while true
+ it += 1
+ it >= bnbparams.maxiter && error("activetset : iteration limit")
+
+ Axbox = solve_kkt_system!(node,asviews,B,y,λ,M)
+ u, Atu, lb = discrete_line_search!(node,asviews,A,Axbox,y,λ,M,bnbparams)
+ swapped = update_sets!(node,asviews,A,B,M,bnbparams)
+
+ if bnbparams.screening
+ prune = screen!(tree,node,asviews,u,Atu,B,A,y,λ,M)
+ prune && break
+ end
+
+ if active_optimality(asviews,M,bnbparams)
+ imax, smax = inactive_optimality(node,asviews,Atu,λ,M,bnbparams)
+ isa(imax,Nothing) && break
+ new_tentative_sign!(node,asviews,imax,smax,A,B)
+ detect_blocking_behaviour(swapped,imax) && break
+ end
+ end
+
+ node.lb = lb
+ node.x[node.support.Sbar0] .= 0.
+ node.x[node.support.Sbarin] = asviews.xSbarin
+ node.x[node.support.Sbarbox] = asviews.xSbarbox
+ node.x[node.support.S1in] = asviews.xS1in
+ node.x[node.support.S1box] = asviews.xS1box
+
+ return nothing
+end
+
diff --git a/src/bnb.jl b/src/bnb.jl
new file mode 100644
index 0000000..ec8819f
--- /dev/null
+++ b/src/bnb.jl
@@ -0,0 +1,191 @@
+function best_lb(tree::Tree)
+ isempty(tree.active) && return 0.
+ return maximum([node.lb for node in tree.active])
+end
+
+function branch!(tree::Tree, node::Node)
+ prune(tree,node) && return nothing
+ j = branching_index(node)
+ isa(j,Nothing) && return nothing
+ node_j0 = Node(node,j,0)
+ node_j1 = Node(node,j,1)
+ push!(tree.queue,node_j0)
+ push!(tree.queue,node_j1)
+end
+
+function branching_index(node::Node)
+ isempty(node.branch.Sbar) && return nothing
+ return node.branch.Sbar[argmax(abs.(node.x[node.branch.Sbar]))]
+end
+
+function depth(node::Node)
+ depth = 0
+ while !isa(node,Nothing)
+ depth += 1
+ node = node.parent
+ end
+ return depth
+end
+
+function elapsed_time(tree::Tree)
+ return (Dates.now() - tree.start_time).value / 1000.
+end
+
+function fixto!(node::Node, j::Int, jval::Int)
+ jpos = findfirst(i->i==j,node.branch.Sbar)
+ isa(jpos,Nothing) && error("Branching index is already fixed")
+ if jval == 0
+ push!(node.branch.S0,popat!(node.branch.Sbar,jpos))
+ node.x[j] = 0.
+ elseif jval == 1
+ push!(node.branch.S1,popat!(node.branch.Sbar,jpos))
+ end
+ return nothing
+end
+
+function global_gap(tree::Tree)
+ return (tree.ub - best_lb(tree)) / tree.ub
+end
+
+function has_perfect_relaxation(node::Node, M::Float64, bnbparams::BnbParams)
+ return !(any(bnbparams.tol .< abs.(node.x[node.branch.Sbar]) .< M - bnbparams.tol))
+end
+
+function lower_bound!(
+ tree::Tree,
+ node::Node,
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ bnbparams::BnbParams
+ )
+ if bnbparams.lb_method == ACTIVESET
+ solve_relax_activeset!(tree,node,A,y,λ,M,bnbparams)
+ elseif bnbparams.lb_method == LBMIP
+ solve_relax_mip!(node,A,y,λ,M)
+ else
+ error("Paramter `bnbparams.lb_method` not understood.")
+ end
+ return nothing
+end
+
+function next_node!(tree::Tree)
+ return pop!(tree.queue)
+end
+
+function processed!(tree::Tree, node::Node)
+ if !prune(tree,node)
+ push!(tree.active,node)
+ end
+ tree.nexpl += 1
+ return nothing
+end
+
+function prune(tree::Tree, node::Node)
+ return (tree.ub < node.lb)
+end
+
+function terminated!(tree::Tree, bnbparams::BnbParams)
+ if isempty(tree.queue)
+ tree.termination_status = MOI.OPTIMAL
+ elseif elapsed_time(tree) > bnbparams.maxtime
+ tree.termination_status = MOI.TIME_LIMIT
+ end
+ return (tree.termination_status != MOI.OPTIMIZE_NOT_CALLED)
+end
+
+function update_upper_bound!(tree::Tree, ub::Float64, inc::Vector{Float64})
+ if ub < tree.ub
+ tree.ub = ub
+ tree.inc = inc
+ return true
+ end
+ return false
+end
+
+function update_active_nodes!(tree::Tree)
+ filter!(node->!prune(tree,node),tree.active)
+ return nothing
+end
+
+function upper_bound!(
+ tree::Tree,
+ node::Node,
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ bnbparams::BnbParams
+ )
+
+ perfect = false
+
+ if isempty(node.branch.S1)
+ ub = 0.5*y'*y
+ inc = zeros(Float64,size(A)[2])
+ elseif has_perfect_relaxation(node,M,bnbparams)
+ ub = node.lb
+ inc = node.x
+ perfect = true
+ else
+ if bnbparams.ub_method == BVLS
+ ub, inc = solve_heuristic_bvls(node,A,y,λ,M,bnbparams)
+ elseif bnbparams.ub_method == UBMIP
+ ub, inc = solve_heuristic_mip(node,A,y,λ,M)
+ else
+ error("Paramter `bnbparams.ub_method` not understood.")
+ end
+ end
+ updated = update_upper_bound!(tree,ub,inc)
+ updated && update_active_nodes!(tree)
+ return perfect
+end
+
+# ============================================================================ #
+
+"""
+ solve_bnb(
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64;
+ bnbparams::BnbParams=BnbParams(),
+ )
+
+Solve the L0-penalized least-square problem
+
+ min_x (1/2) norm(y-Ax,2)^2 + λ norm(x,0)
+ st. norm(x,Inf) <= M
+
+using a Branch-and-Bound algorithm. This algorithm can use branch-screening
+tests to enhance its efficiency. Returns a [`@BnbResults`](@ref) struct
+containing the optimization results and statistics. See [`@BnbParams`](@ref) for
+informations about additional method parameters.
+"""
+function solve_bnb(
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64;
+ bnbparams::BnbParams=BnbParams(),
+ )
+
+ validate_data(A,y,λ,M)
+ tree = Tree(A,y)
+ bnbparams.verbosity && display_head()
+
+ while !terminated!(tree,bnbparams)
+ node = next_node!(tree)
+ lower_bound!(tree,node,A,y,λ,M,bnbparams)
+ perfect = upper_bound!(tree,node,A,y,λ,M,bnbparams)
+ processed!(tree,node)
+ perfect || branch!(tree,node)
+ bnbparams.verbosity && display_node(tree,node)
+ end
+
+ results = BnbResults(tree)
+ bnbparams.verbosity && display_tail(results)
+
+ return results
+end
\ No newline at end of file
diff --git a/src/bvls.jl b/src/bvls.jl
new file mode 100644
index 0000000..b7a3aaa
--- /dev/null
+++ b/src/bvls.jl
@@ -0,0 +1,132 @@
+function compute_kkt_optimality(g::Vector{Float64}, on_bound::Vector{Int})
+ length(g) == 0 && return 0.
+ @views begin
+ g_kkt = g .* on_bound
+ free_set = (on_bound .== 0)
+ g_kkt[free_set] = abs.(g[free_set])
+ end
+ return maximum(g_kkt)
+end
+
+function solve_heuristic_bvls(
+ node::Node,
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ bnbparams::BnbParams,
+ )
+
+ AS1 = A[:,node.branch.S1]
+ n = size(AS1)[2]
+ lb = fill(-M,(n,))
+ ub = fill(M,(n,))
+
+ x = zeros(Float64,n)
+ on_bound = zeros(Int,n)
+ free_set = (on_bound .== 0)
+ active_set = .!free_set
+ free_set = findall(free_set)
+
+ r = AS1 * x - y
+ cost = r' * r
+ g = AS1' * r
+ cost_change = Inf
+
+ # ----- Initialization loop ----- #
+
+ while !isempty(free_set)
+
+ AS1_free = AS1[:, free_set]
+ y_free = y - AS1 * (x .* active_set)
+ z = AS1_free \ y_free
+ lbv = (z .< lb[free_set])
+ ubv = (z .> ub[free_set])
+ v = lbv .| ubv
+
+ if any(lbv)
+ ind = free_set[lbv]
+ x[ind] = lb[ind]
+ active_set[ind] .= true
+ on_bound[ind] .= -1
+ end
+ if any(ubv)
+ ind = free_set[ubv]
+ x[ind] = ub[ind]
+ active_set[ind] .= true
+ on_bound[ind] .= 1
+ end
+
+ ind = free_set[.!v]
+ x[ind] = z[.!v]
+ r = AS1 * x - y
+ g = AS1' * r
+ any(v) || break
+ free_set = free_set[.!v]
+ end
+
+ # ----- Main loop ----- #
+
+ optimality = compute_kkt_optimality(g,on_bound)
+
+ while true
+
+ (optimality < bnbparams.tol) && break
+ move_to_free = argmax(g .* on_bound)
+ on_bound[move_to_free] = 0
+
+ while true
+
+ free_set = (on_bound .== 0)
+ active_set = .!free_set
+ free_set = findall(free_set)
+ x_free = x[free_set]
+ lb_free = lb[free_set]
+ ub_free = ub[free_set]
+ AS1_free = AS1[:,free_set]
+ y_free = y - AS1 * (x .* active_set)
+ z = AS1_free \ y_free
+ lbv = findall(z .< lb_free)
+ ubv = findall(z .> ub_free)
+ v = vcat(lbv,ubv)
+
+ if !isempty(v)
+
+ alphas = vcat(
+ lb_free[lbv] - x_free[lbv],
+ ub_free[ubv] - x_free[ubv]
+ ) ./ (z[v] - x_free[v])
+ i = argmin(alphas)
+ i_free = v[i]
+ alpha = alphas[i]
+ x_free *= 1. - alpha
+ x_free += alpha * z
+ x[free_set] = x_free
+
+ if i < length(lbv)
+ on_bound[free_set[i_free]] = -1
+ else
+ on_bound[free_set[i_free]] = 1
+ end
+ else
+ x_free = z
+ x[free_set] = x_free
+ break
+ end
+ end
+
+ r = AS1 * x - y
+ cost_new = 0.5 * (r' * r)
+ cost_change = cost - cost_new
+
+ (cost_change < bnbparams.tol * cost) && break
+
+ cost = cost_new
+ g = AS1' * r
+ optimality = compute_kkt_optimality(g,on_bound)
+ end
+
+ obj = 0.5 * norm(r,2)^2 + λ * n
+
+ return obj, x
+end
\ No newline at end of file
diff --git a/src/data.jl b/src/data.jl
new file mode 100644
index 0000000..d1d1f55
--- /dev/null
+++ b/src/data.jl
@@ -0,0 +1,87 @@
+function validate_data(
+ A::Array{Float64,2},
+ y::Array{Float64,1},
+ λ::Float64,
+ M::Float64
+ )
+ size(A)[1] == size(y)[1] || throw(ArgumentError("Dimension missmatch"))
+ λ > 0 || throw(ArgumentError("λ must be positive"))
+ M > 0 || throw(ArgumentError("M must be positive"))
+end
+
+"""
+ data_gaussian(
+ k::Int=5,
+ m::Int=500,
+ n::Int=1000,
+ snr::Float64=10.
+ )
+
+Generate a Gaussian setup.
+"""
+function data_gaussian(
+ k::Int=5,
+ m::Int=500,
+ n::Int=1000,
+ snr::Float64=10.
+ )
+
+ x = zeros(n)
+ a = rand(Normal(0,1),k)
+ a += sign.(a)
+ x[randperm(n)[1:k]] = a
+
+ A = rand(Normal(0,1),m,n)
+ for i in 1:n
+ A[:,i] /= norm(A[:,i],2)
+ end
+
+ σ = norm(A*x,2) / sqrt(m * 10^(snr/10))
+ ϵ = rand(Normal(0,σ),m)
+ y = A * x .+ ϵ
+ λ = 2 * σ * log(n/k - 1)
+ M = 1.5 * norm(A'*y,Inf)
+
+ return A, y, λ, M
+end
+
+"""
+ data_gaussian(
+ k::Int=5,
+ m::Int=500,
+ n::Int=1000,
+ snr::Float64=10.
+ )
+
+Generate a Toeplitz setup.
+"""
+function data_toeplitz(
+ k::Int=5,
+ m::Int=500,
+ n::Int=300,
+ snr::Float64=10.
+ )
+
+ x = zeros(n)
+ a = rand(Normal(0,1),k)
+ a += sign.(a)
+ x[randperm(n)[1:k]] = a
+
+ tmax = n/2
+ t = collect(1:m)
+ A = zeros(m,n)
+ offset = (n-m)/2
+ for i in 1:n
+ A[:,i] = @. sinc((offset + (t - i)) * (tmax/m))
+ A[:,i] /= norm(A[:,i],2)
+ end
+ A += rand(m,n)./10000
+
+ σ = norm(A*x,2) / sqrt(m * 10^(snr/10))
+ ϵ = rand(Normal(0,σ),m)
+ y = A * x .+ ϵ
+ λ = 2 * σ * log(n/k - 1)
+ M = 1.5 * norm(A'*y,Inf)
+
+ return A, y, λ, M
+end
\ No newline at end of file
diff --git a/src/displays.jl b/src/displays.jl
new file mode 100644
index 0000000..391dde8
--- /dev/null
+++ b/src/displays.jl
@@ -0,0 +1,53 @@
+WIDTH = 67
+
+function display_head()
+ println(repeat("*",WIDTH))
+ println("B&B-screening solver for the L0-penalized least-square problem")
+ println("Version : $(version())")
+ println("Author : $(author())")
+ println("Contact : $(contact())")
+ println("License : $(license())")
+ println()
+ println(repeat("-",WIDTH))
+ println("| Node | Bounds | Branch |")
+ println("| nexpl depth queue| lb ub gap| S0 S1 Sbar|")
+ println(repeat("-",WIDTH))
+ return nothing
+end
+
+function display_node(tree::Tree, node::Node)
+
+ @printf(
+ "|%7d%7d%7d|%7.3f%7.3f%7.0e|%7d%7d%7d|",
+ tree.nexpl,
+ depth(node),
+ length(tree.queue),
+ best_lb(tree),
+ tree.ub,
+ global_gap(tree),
+ length(node.branch.S0),
+ length(node.branch.S1),
+ length(node.branch.Sbar),
+ )
+ println()
+ return nothing
+end
+
+function diplay_results(results::BnbResults)
+ println("B&B results")
+ println(" Status : $(results.termination_status)")
+ println(" Solve time : $(results.solve_time)")
+ println(" Nodes : $(results.node_count)")
+ println(" Objective : $(results.objective_value)")
+ return nothing
+end
+
+function display_tail(results::BnbResults)
+ println(repeat("-",WIDTH))
+ println()
+ diplay_results(results)
+ println()
+ println(repeat("*",WIDTH))
+ return nothing
+end
+
diff --git a/src/mips.jl b/src/mips.jl
new file mode 100644
index 0000000..7199214
--- /dev/null
+++ b/src/mips.jl
@@ -0,0 +1,117 @@
+function solve_relax_mip!(
+ node::Node,
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ )
+
+ relax = Model(CPLEX.Optimizer)
+ set_silent(relax)
+
+ gram = A'*A
+ Aty = A'*y
+
+ n = size(A)[2]
+ @variable(relax, x[1:n])
+ @variable(relax, xabs[1:n])
+ @objective(relax, Min,
+ 0.5*(y'*y) - Aty'*x + (0.5*x)'*(gram*x) +
+ (λ/M)*sum(xabs[node.branch.Sbar]) +
+ λ*length(node.branch.S1)
+ )
+ @constraint(relax, box, -M .<= x .<= M)
+ @constraint(relax, abs_pos, xabs .>= x)
+ @constraint(relax, abs_neg, xabs .>= -x)
+ @constraint(relax, setzero, x[node.branch.S0] .== 0)
+ @constraint(relax, abs_setzero, xabs[node.branch.S0] .== 0)
+
+ for i in 1:n
+ set_start_value(x[i], node.x[i])
+ set_start_value(xabs[i], abs(node.x[i]))
+ end
+
+ optimize!(relax)
+
+ node.lb = objective_value(relax)
+ node.x = value.(x)
+end
+
+function solve_heuristic_mip(
+ node::Node,
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ )
+
+ heur = Model(CPLEX.Optimizer)
+ set_silent(heur)
+
+ nS1 = length(node.branch.S1)
+ AS1 = A[:,node.branch.S1]
+ xS1 = node.x[node.branch.S1]
+
+ @variable(heur, x[1:nS1])
+ @objective(heur, Min, 0.5*(y-AS1*x)'*(y-AS1*x)+ λ*length(node.branch.S1))
+ @constraint(heur, box, -M .<= x .<= M)
+
+ for i in 1:nS1
+ set_start_value(x[i], xS1[i])
+ end
+
+ optimize!(heur)
+
+ inc = zeros(Float64,size(A)[2])
+ inc[node.branch.S1] = value.(x)
+
+ return objective_value(heur), inc
+end
+
+# ============================================================================ #
+
+"""
+ solve_mip(
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64;
+ mipparams::MipParams=MipParams(),
+ )
+
+Solve the L0-penalized least-square problem
+
+ min_x (1/2) norm(y-Ax,2)^2 + λ norm(x,0)
+ st. norm(x,Inf) <= M
+
+using a MIP solver (CPLEX). Returns a [`@MipResults`](@ref) struct containing
+the optimization results and statistics. See [`@MipParams`](@ref) for
+informations about additional method parameters.
+"""
+function solve_mip(
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64;
+ mipparams::MipParams=MipParams(),
+ )
+
+ problem = Model(CPLEX.Optimizer)
+ set_optimizer_attribute(problem,"CPX_PARAM_TILIM",mipparams.maxtime)
+ mipparams.silent && set_silent(problem)
+
+ gram = A'*A
+ Aty = A'*y
+
+ n = size(A)[2]
+ @variable(problem, x[1:n])
+ @variable(problem, z[1:n], Bin)
+ @objective(problem, Min, 0.5*(y'*y) - Aty'*x + (0.5*x)'*(gram*x) + λ*sum(z))
+ @constraint(problem, box_lower, -M.*z .<= x)
+ @constraint(problem, box_upper, x .<= M.*z)
+
+ optimize!(problem)
+ result = MipResults(problem)
+
+ return result
+end
diff --git a/src/screen.jl b/src/screen.jl
new file mode 100644
index 0000000..569c17d
--- /dev/null
+++ b/src/screen.jl
@@ -0,0 +1,52 @@
+function screen!(
+ tree::Tree,
+ node::Node,
+ asviews::ActivesetViews,
+ u::Vector{Float64},
+ Atu::Vector{Float64},
+ B::Matrix{Float64},
+ A::Matrix{Float64},
+ y::Vector{Float64},
+ λ::Float64,
+ M::Float64,
+ )
+
+ @views begin
+ pivotSbar = @. M * (abs(Atu[node.branch.Sbar]) - λ/M)
+ pivotSbar0 = max.(pivotSbar,0)
+ pivotSbar1 = max.(-pivotSbar,0)
+ pivotS1 = @. M * (abs(Atu[node.branch.S1]) - λ/M)
+ dobj = 0.5 * (y'y - norm(y-u,2)^2) - sum(pivotSbar0) - sum(pivotS1)
+ end
+
+ Sbar_to_S0 = Vector{Int}()
+ Sbar_to_S1 = Vector{Int}()
+ for (iSbar,i) in enumerate(node.branch.Sbar)
+ if pivotSbar1[iSbar] > tree.ub - dobj
+ push!(Sbar_to_S0,i)
+ tree.scr0 += 1
+ end
+ if pivotSbar0[iSbar] > tree.ub - dobj
+ push!(Sbar_to_S1,i)
+ tree.scr1 += 1
+ end
+ end
+
+ if !isempty(intersect(Sbar_to_S0,Sbar_to_S1))
+ node.lb = Inf
+ node.x = rand(size(A)[2])
+ return true
+ end
+
+ for i in Sbar_to_S0
+ fixto!(node,i,0)
+ swap_Sbar_to_S0!(node,asviews,A,i)
+ end
+
+ for i in Sbar_to_S1
+ fixto!(node,i,1)
+ swap_Sbar_to_S1!(node,asviews,A,B,i)
+ end
+
+ return false
+end
\ No newline at end of file
diff --git a/src/types.jl b/src/types.jl
new file mode 100644
index 0000000..16e48d4
--- /dev/null
+++ b/src/types.jl
@@ -0,0 +1,295 @@
+# ========================================== #
+# Active-Set support configuration and views #
+# ========================================== #
+
+mutable struct Support
+ Sbar0::Vector{Int} # Indices verifying 0 = xi with i in Sbar
+ Sbarin::Vector{Int} # Indices verifying 0 < xi < M with i in Sbar
+ Sbarbox::Vector{Int} # Indices verifying xi = M with i in Sbar
+ S1in::Vector{Int} # Indices verifying xi < M with i in S1
+ S1box::Vector{Int} # Indices verifying xi = M with i in S1
+end
+
+mutable struct ActivesetViews
+ # Current iterate
+ xSbarin::Vector{Float64}
+ xSbarbox::Vector{Float64}
+ xS1in::Vector{Float64}
+ xS1box::Vector{Float64}
+
+ # Sign of x
+ sSbarin::Vector{Float64}
+ sSbarbox::Vector{Float64}
+ sS1in::Vector{Float64}
+ sS1box::Vector{Float64}
+
+ # New iterate
+ xnewSbarin::Vector{Float64}
+ xnewS1in::Vector{Float64}
+
+ # Submatrices
+ ASbarin::Matrix{Float64}
+ ASbarbox::Matrix{Float64}
+ AS1in::Matrix{Float64}
+ AS1box::Matrix{Float64}
+end
+
+
+# =========================== #
+# Branch-and-Bound components #
+# =========================== #
+
+mutable struct Branch
+ S0::Vector{Int} # Indices fixed to zero
+ S1::Vector{Int} # Indices fixed to non-zero
+ Sbar::Vector{Int} # Indices not fixed
+end
+
+mutable struct Node
+ parent::Union{Nothing,Node} # Parent node
+ branch::Branch # Current branch (fixed variables)
+ branch_on::Tuple{Int,Int} # Branching index above node
+ x::Vector{Float64} # Relaxation solution
+ lb::Float64 # Relaxation objective value
+
+ # Active set data cashed in each node
+ support::Support
+ GS1in_inv::Matrix{Float64}
+ GSbarin_inv::Matrix{Float64}
+end
+
+mutable struct Tree
+ termination_status::MOI.TerminationStatusCode # Optimization status
+ start_time::DateTime # Start time
+ nexpl::Int # Number of nodes explored
+ queue::Vector{Node} # Nodes to explore
+ active::Vector{Node} # Active nodes (with node.lb < tree.ub)
+ ub::Float64 # Best upper bound
+ inc::Vector{Float64} # Incumbent solution
+ scr0::Int # Indices fixed to 0 with screening tests
+ scr1::Int # Indices fixed to 1 with screening tests
+end
+
+
+# ======================================= #
+# Branch-and-Bound parameters and outputs #
+# ======================================= #
+
+# Available solving methods for the relaxed problem (lower bounding step)
+@enum LbMethod begin
+ LBMIP # MIP solver
+ ACTIVESET # Active-Set algorihtm
+end
+
+# Available heuristic methods (upper bounding step)
+@enum UbMethod begin
+ UBMIP # MIP solver
+ BVLS # Bounded Variable Least-Square algorithm
+end
+
+struct BnbParams
+ lb_method::LbMethod # Method used to solve node relaxations
+ ub_method::UbMethod # Heuristic to compute node uppper bounds
+ maxtime::Int # Maximum solving time in seconds
+ tol::Float64 # Integer tolerance
+ maxiter::Int # Max. iter. of the relaxation solution algorithm
+ screening::Bool # Toogle screening
+ verbosity::Bool # Whether to enable iterations displays
+ showevery::Int # Node display interval
+end
+
+struct BnbResults
+ termination_status::MOI.TerminationStatusCode # Optimization status
+ solve_time::Float64 # Solution time in seconds
+ node_count::Int # Number of nodes explored
+ objective_value::Float64 # Objective value (if optimal)
+ x::Vector{Float64} # Optimizer (if optimal)
+ scr0::Int # Number of screning tests passed (xi = 0)
+ scr1::Int # Number of screning tests passed (xi != 0)
+end
+
+
+# ================================= #
+# MIP solver parameters and outputs #
+# ================================= #
+
+struct MipParams
+ maxtime::Int # Maximum solution time in seconds
+ silent::Bool # Set solver silent
+end
+
+struct MipResults
+ termination_status::MOI.TerminationStatusCode # Optimization status
+ solve_time::Float64 # Solution time in seconds
+ node_count::Int # Number of nodes explored
+ objective_value::Float64 # Objective value (if optimal)
+ x::Vector{Float64} # Optimizer (if optimal)
+end
+
+
+# ============ #
+# Constructors #
+# ============ #
+
+function Support(n::Int)
+ return Support(
+ collect(1:n),
+ Vector{Int}(),
+ Vector{Int}(),
+ Vector{Int}(),
+ Vector{Int}()
+ )
+end
+
+function ActivesetViews(node::Node, A::Matrix{Float64}, bnbparams::BnbParams)
+ xSbarin = node.x[node.support.Sbarin]
+ xSbarbox = node.x[node.support.Sbarbox]
+ xS1in = node.x[node.support.S1in]
+ xS1box = node.x[node.support.S1box]
+ sSbarin = tolsign.(xSbarin, bnbparams.tol)
+ sSbarbox = tolsign.(xSbarbox, bnbparams.tol)
+ sS1in = tolsign.(xS1in, bnbparams.tol)
+ sS1box = tolsign.(xS1box, bnbparams.tol)
+ xnewSbarin = copy(xSbarin)
+ xnewS1in = copy(xS1in)
+ ASbarin = A[:,node.support.Sbarin]
+ ASbarbox = A[:,node.support.Sbarbox]
+ AS1in = A[:,node.support.S1in]
+ AS1box = A[:,node.support.S1box]
+ return ActivesetViews(
+ xSbarin,
+ xSbarbox,
+ xS1in,
+ xS1box,
+ sSbarin,
+ sSbarbox,
+ sS1in,
+ sS1box,
+ xnewSbarin,
+ xnewS1in,
+ ASbarin,
+ ASbarbox,
+ AS1in,
+ AS1box,
+ )
+end
+
+function Branch(n::Int)
+ return Branch(Vector{Int}(), Vector{Int}(), collect(1:n))
+end
+
+function Node(A::Matrix{Float64})
+ m, n = size(A)
+ return Node(
+ nothing,
+ Branch(n),
+ (-1,-1),
+ zeros(Float64,n),
+ 0.,
+ Support(n),
+ zeros(Float64,0,0),
+ zeros(Float64,0,0),
+ )
+end
+
+function Node(parent::Node, j::Int, jval::Int)
+ node = Node(
+ parent,
+ copy(parent.branch),
+ (j,jval),
+ copy(parent.x),
+ parent.lb,
+ copy(parent.support),
+ copy(parent.GS1in_inv),
+ copy(parent.GSbarin_inv),
+ )
+ fixto!(node,j,jval)
+ return node
+end
+
+function Tree(A::Array{Float64,2}, y::Vector{Float64})
+ return Tree(
+ MOI.OPTIMIZE_NOT_CALLED,
+ Dates.now(),
+ 0,
+ Vector{Node}([Node(A)]),
+ Vector{Node}(),
+ 0.5 * norm(y, 2)^2,
+ zeros(Float64,size(A)[2]),
+ 0,
+ 0,
+ )
+end
+
+function BnbParams(;
+ lb_method=ACTIVESET,
+ ub_method=BVLS,
+ maxtime=1000,
+ tol=1e-8,
+ maxiter=1000,
+ screening=false,
+ verbosity=true,
+ showevery=1,
+ )
+ maxtime > 0 || error("Parameter `maxtime` must be positive")
+ tol > 0 || error("Parameter `tol` must be positive")
+ maxiter > 0 || error("Parameter `maxiter` must be positive")
+ showevery > 0 || error("Parameter `showevery` must be positive")
+ return BnbParams(
+ lb_method,
+ ub_method,
+ maxtime,
+ tol,
+ maxiter,
+ screening,
+ verbosity,
+ showevery
+ )
+end
+
+function BnbResults(tree::Tree)
+ return BnbResults(
+ tree.termination_status,
+ elapsed_time(tree),
+ tree.nexpl,
+ tree.ub,
+ tree.inc,
+ tree.scr0,
+ tree.scr1,
+ )
+end
+
+function MipParams(;maxtime=1000,silent=true)
+ maxtime > 0 || error("Parameter `maxtime` must be positive")
+ return MipParams(maxtime, silent)
+end
+
+function MipResults(problem::JuMP.Model)
+ return MipResults(
+ termination_status(problem),
+ solve_time(problem),
+ convert(Int,node_count(problem)),
+ objective_value(problem),
+ value.(problem[:x])
+ )
+end
+
+
+# ====================== #
+# Base module extensions #
+# ====================== #
+
+
+function Base.copy(support::Support)
+ return Support(
+ copy(support.Sbar0),
+ copy(support.Sbarin),
+ copy(support.Sbarbox),
+ copy(support.S1in),
+ copy(support.S1box)
+ )
+end
+
+function Base.copy(branch::Branch)
+ return Branch(copy(branch.S0), copy(branch.S1), copy(branch.Sbar))
+end
diff --git a/tests/test.jl b/tests/test.jl
new file mode 100644
index 0000000..b204ae1
--- /dev/null
+++ b/tests/test.jl
@@ -0,0 +1,31 @@
+include("../src/BnbScreening.jl")
+
+k, m, n, snr = 5, 500, 300, 10*log10(10)
+A, y, λ, M = data_toeplitz(k,m,n,snr)
+#k, m, n, snr = 5, 500, 1000, 10*log10(10)
+#A, y, λ, M = data_gaussian(k,m,n,snr)
+
+result_mip = solve_mip(A,y,λ,M)
+println("MIP")
+println(" Status : $(result_mip.termination_status)")
+println(" Objective : $(result_mip.objective_value)")
+println(" Nodes : $(result_mip.node_count)")
+println(" Timer : $(result_mip.solve_time)")
+
+result_bnb = solve_bnb(A,y,λ,M,bnbparams=BnbParams(verbosity=false))
+println("BnB")
+println(" Status : $(result_bnb.termination_status)")
+println(" Objective : $(result_bnb.objective_value)")
+println(" Nodes : $(result_bnb.node_count)")
+println(" Timer : $(result_bnb.solve_time)")
+println("Error : $(abs(result_mip.objective_value - result_bnb.objective_value))")
+
+result_bnb_scr = solve_bnb(A,y,λ,M,bnbparams=BnbParams(verbosity=false,screening=true))
+println("BnB + scr")
+println(" Status : $(result_bnb_scr.termination_status)")
+println(" Objective : $(result_bnb_scr.objective_value)")
+println(" Nodes : $(result_bnb_scr.node_count)")
+println(" Timer : $(result_bnb_scr.solve_time)")
+println(" Scr-0 : $(result_bnb_scr.scr0)")
+println(" Scr-1 : $(result_bnb_scr.scr1)")
+println("Error : $(abs(result_mip.objective_value - result_bnb_scr.objective_value))")
\ No newline at end of file
--
GitLab