top-k-sum

Routines for Euclidean projection of an $n$-dimensional vector $x^0$ onto the top- $\!\!k$-sum constraint math \mathcal{B}_{(k)}^r \coloneqq \Bigl\{x\in\mathbb{R}^n : \sum_{i=1}^k x_{[i]} \leq r\Bigr\}, where $k\in\{1,\ldots,n\}$ and $r\in\mathbb{R}$ are parameters and where $x_{[i]}$ denotes the $i^{\text{th}}$ largest element of $x$.

The repository contains two directories related to the paper https://arxiv.org/abs/2310.07224: - src/: implementations of the algorithms. - run/: scripts for running the experiments. Experiments can be run from the run/ folder in terminal as $ julia run_tks_experiments.jl by specifying a proper DATAPATH variable.

Routines in src/ include three finite-termination algorithms: - ESGS: an early-stopping grid-search approach callable by project_topksum_esgs!; - PLCP: a parametric-LCP approach callable by project_topksum_plcp!; - GRID: a full grid-search approach callable by project_topksum_grid!;

and two inexact, QP solvers: - GRBS: a QP barrier method using Gurobi based on the sorted formulation, and callable by project_topksum_grbs; - GRBU: a QP barrier method using Gurobi based on the unsorted formulation, and callable by project_topksum_grbu. Note: there is not a native quickselect procedure in Julia, so GRBU assumes that the largest $k$ elements of the input vector are sorted.

An example calling sequence for ESGS and the Gurobi-based methods is provided below.

``` using Random Random.seed!(1234567)

using Gurobi, SparseArrays, SparseMatricesCSR # for Gurobi methods

n = 10; k = 4; r = 0.1234567; x0 = randn(n); pi = sortperm(x0, rev=true); # full sort pipartial = sortperm(x0, alg=PartialQuickSort(1:k), rev=true); # partial sort x0sort = x0[pi]; # full sort x0psort = x0[pipartial]; # partial sort

active = sum(x0sort[1:k]) > r; # binary: is the constraint active xbarsort = similar(x0sort); xbarsort .= x0sort; prepop = true; # binary: is the solution vector pre-populated with elements from sorted input vector outesgs = projecttopksumesgs!(xbarsort, x0sort, r, k, active, prepop); # solve with ESGS method outgrbs = projecttopksumgrbs(x0sort, r, k); # solve with GRBS method outgrbu = projecttopksumgrbu(x0psort, r, k); # solve with GRBU (partial sort) method which results in julia> # display sorted solution vcat(["sorted input" "esgs" "grbs" "grbu"], hcat(x0sort, xbarsort, outgrbs[:x], out_grbu[:xsort])) 11×4 Matrix{Any}: "sorted input" "esgs" "grbs" "grbu" 2.0882 0.762312 0.762312 0.762312 1.01192 -0.212952 -0.212952 -0.212952 1.00237 -0.212952 -0.212952 -0.212952 0.577227 -0.212952 -0.212952 -0.212952 0.377373 -0.212952 -0.212952 -0.212952 -0.0814697 -0.212952 -0.212952 -0.212952 -0.187464 -0.212952 -0.212952 -0.212952 -0.47593 -0.47593 -0.47593 -0.47593 -0.754928 -0.754928 -0.754928 -0.754928 -1.41154 -1.41154 -1.41154 -1.41154

julia> # verify feasibility (sum(xbarsort[1:k]), sum(outgrbs[:x][1:k]), sum(outgrbu[:xsort][1:k]), r) (0.12345669999999945, 0.12345669999491526, 0.12345669999960365, 0.1234567)

julia> # display objective value objvals = (0.5sum((xbarsort .- x0sort).^2), 0.5sum((outgrbs[:x] .- x0sort).^2), 0.5sum((outgrbu[:xsort] .- x0sort).^2)); julia> minobj = minimum(objvals); julia> objvals .- minobj (0.0, 2.2381181352670865e-9, 2.375699637013895e-10)

julia> # display index-pairs (k0,k1) (outesgs[2], getk0k1(outgrbs[:x], k), getk0k1(outgrbu[:xsort], k)) ((1, 7), (3, 4), (3, 4)) ``` Note that the Gurobi-based methods return incorrect index-pairs, even in an $n=10$ dimensional problem. In general, it is not obvious how to recover appropriate $(k0,k_1)$ information from an inexact solution.