Find the max element of a numeric vector iteratively with mask and early break

[WalterMadelim]

This is an elementary question about algorithm.
My algorithm is probably already usable but
I wonder if there is a faster algorithm that can be implemented with Julia (e.g. even not using findmax)

import Random; Random.seed!(6);
N = 10; # length of the value vector
val_vec = rand(-9:.1:9, N); # generate a test value vector
bit_vec = trues(N); # entry == 1 implies it is not fathomed yet

function is_task_accomplished_at(ind) # representing an external function
    # The body is merely an artificial one
    # for demo purpose only
    rand() > .95 && return true 
    return false
end

let # This let-block is for reference purpose only
    # It is not an algorithm because presumably `sort` is slow when N is large
    println("val_vec is $val_vec")
    println("sorted val_vec is $(reverse(sort(val_vec)))")
    println("sortperm is $(reverse(sortperm(val_vec)))")
end

# This is my algorithm
while true
    if any(bit_vec) == false # This if-block is the motivation of introducing `bit_vec`
        @info "eventually return due to all-fathoming"
        return
    end
    ind = findmax(val_vec)[2]
    println("currently visiting $ind, whose max value is $(val_vec[ind])")
    if is_task_accomplished_at(ind)
        # This early break might happen very soon at the early period of the algorithm
        # Therefore using `sort` is not economical
        @info "early break due to accomplishment"
        break
    end
    bit_vec[ind] = false
    val_vec[ind] = -Inf # mask this entry due to fathomed
end

[JADekker]

I’m not entirely sure what exactly it is that you are trying to do, but can’t you get the indices for sorting the vector once outside the loop and loop through these? This replaces repeated findmax’es with one sort, which I think should be much faster? You can then also drop bitvec (as this corresponds simply with running through the entire vector) and I think you may even not need to write the -Inf’s (unless you need them later on)

[WalterMadelim]

@JADekker

This replaces repeated findmax’es with one sort

Thanks. I’m not a computer science student, I though findmax should be way more efficient than sort? (I’m not sure, and not professional).

@JADekker

you may even not need to write the -zing’s

Sorry, I didn’t get that (What’s -zing’s).

[JADekker]

My bad, I meant “-Inf’s”, but apparently autocorrect on my phone wanted to write -zing’s instead…

[JADekker]

If you only want to find the maximum once, I think findmax is indeed more efficient. However, you are interested in running through all values from highest to lowest (until some condition is met), which means that repeated calls to findmax require you to loop through all values often, whereas sorting only requires you to do this once.

It does depend a bit on how many iterations you expect to need and how long the vectors are, it might sometimes be quicker to use findmax if you only need a few iterations. However, I think using sort (or sort!) avoids repeated findmax calls, which I guess has good potential to be more efficient.

[WalterMadelim]

@JADekker

I’m not entirely sure what exactly it is that you are trying to do,

I have N models to train, from nascent to eventually mature.

At the early stages, I can even randomly pick (instead of findmax) one index, then strictly enrich that model, thus the early break (There is still one outer loop enclosing that while in my code above). So at this time, using sort seems uneconomical. (findmax can be viewed as a good heuristic that the indexed model—which owns the max value—can be strictly enriched).

But as time goes by, the vector of models tend to become more mature, and some of them can no longer be strictly enriched, thus using my algorithm above, I may repeatedly do findmax a lot of times.

Eventually, all models are finally mature and thus that return in my code above (return can get rid of all loops, whereas break can only escape from the cloesest loop.)

---

Here is a 3rd implementation which uses neither sort nor findmax

julia> import Random; Random.seed!(6);

julia> N = 10; # length of the value vector

julia> val_vec = rand(-9:.1:9, N); # generate a test value vector

julia> bit_vec = trues(N); # entry == 1 implies it is not fathomed yet

julia> function drawmax!(bit_vec, val_vec)                                                               
           still_has_undrawn, max_i = false, 0                                                           
           if any(bit_vec) # still has an unfathomed entry                                               
               still_has_undrawn, max_v = true, -Inf                                                     
               for (i, undrawn) in enumerate(bit_vec)                                                    
                   if undrawn && max_v < val_vec[i]                                                      
                       max_v, max_i = val_vec[i], i                                                      
                   end                                                                                   
               end                                                                                       
               bit_vec[max_i] = false # fathomed this time                                               
           end                                                                                           
           return still_has_undrawn, max_i # `max_i` is valid only if `still_has_undrawn` == true        
       end;

julia> @info "for your reference" reverse(sort(val_vec)) reverse(sortperm(val_vec))
┌ Info: for your reference
│   reverse(sort(val_vec)) =
│    10-element Vector{Float64}:
│      8.5
│      6.1
│      1.1
│     -2.2
│     -2.2
│     -5.8
│     -7.1
│     -7.6
│     -8.5
│     -8.8
│   reverse(sortperm(val_vec)) =
│    10-element Vector{Int64}:
│      7
│      1
│      3
│     10
│      8
│      6
│      5
│      2
│      9
└      4

julia> drawmax!(bit_vec, val_vec)
(true, 7)

julia> drawmax!(bit_vec, val_vec)
(true, 1)

julia> drawmax!(bit_vec, val_vec)
(true, 3)

julia> drawmax!(bit_vec, val_vec)
(true, 8)

julia> drawmax!(bit_vec, val_vec)
(true, 10)

julia> drawmax!(bit_vec, val_vec)
(true, 6)

julia> drawmax!(bit_vec, val_vec)
(true, 5)

julia> drawmax!(bit_vec, val_vec)
(true, 2)

julia> drawmax!(bit_vec, val_vec)
(true, 9)

julia> drawmax!(bit_vec, val_vec)
(true, 4)

julia> drawmax!(bit_vec, val_vec)
(false, 0)

[stevengj]

@WalterMadelim

I have N models to train, from nascent to eventually mature.

At the early stages, I can even randomly pick (instead of findmax) one index, then strictly enrich that model, thus the early break (There is still one outer loop enclosing that while in my code above). So at this time, using sort seems uneconomical. (findmax can be viewed as a good heuristic that the indexed model—which owns the max value—can be strictly enriched).

Why not just use a data structure like a heap or a priority queue that allows you to pop off the maximum element, change/enrich it, and push it back? A vector is a very suboptimal data structure for applications where you have a collection of items and want to repeatedly find maximum elements and/or modify one element and/or insert new elements.

(There is a whole subfield of computer science that centers on designing data structures optimized for specific kinds of operations on collections.)

(Also, is the cost of findmax really significant compared to the cost of “enriching” a model, which I guess is some kind of optimization step?)

I’m still pretty confused about what you are actually trying to do here.

[WalterMadelim]

@stevengj

(Also, is the cost of findmax really significant compared to the cost of “enriching” a model, which I guess is some kind of optimization step?)

You are absolutely right. It’s not at the heart of my global algorithm. It’s a heuristic, to identify a promising block.

The (already runnable, and successful) code is here

Cutting Stock Problem-Dual Decomposition: FindMax approach

import JuMP, Gurobi; GRB_ENV = Gurobi.Env();
import LinearAlgebra.⋅ as ⋅
import Random; Random.seed!(1)

# A Heterogeneous Cutting Stock Problem
I, J, w, d, W, c = let # 📄 load data
    δ = 1.7e-12 # used in test case generation
    I = 20; # [demand side] number of types of demand
    J = 400; # [resource side] number of raw rolls
    wl, wh = rand(13:δ:35), rand(77:δ:97)
    w = rand(wl:δ:wh, I) # width of each demand type
    dl, dh = rand(7:14), rand(57:77)
    d = rand(dl:dh, I) # quantity::Int of each demand type
    Wl, Wh = rand(70:δ:100), rand(400:δ:500)
    W = rand(Wl:δ:Wh, J) # width of each raw roll
    c = 0.01rand(0.99:δ:1.99, J) .* W # used cost of raw rolls (c)
    I, J, w, d, W, c
end;

if false # 📘 The original compact Mip formulation
    Mip = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); 
    JuMP.@variable(Mip, 0 ≤ b[1:J] ≤ 1, Int);
    JuMP.@objective(Mip, Min, c ⋅ b);
    JuMP.@variable(Mip, 0 ≤ n[1:I, 1:J], Int);
    JuMP.@constraint(Mip, [j = 1:J], w ⋅ view(n, :, j) ≤ b[j]W[j]);
    JuMP.@constraint(Mip, [i = 1:I], sum(n[i, j] for j in 1:J) == d[i]); # 🟡 [π]
    JuMP.unset_integer.(b)
    JuMP.unset_integer.(n)
    JuMP.optimize!(Mip) # [Natural LP Relaxation] Optimal objective  5.232228077e+02
    JuMP.set_integer.(b)
    JuMP.set_integer.(n)
    JuMP.optimize!(Mip) # The normal solve, # ub = 543.86951 | 523.22976 = lb after 1 seconds
    # ub = 527.85975 | 523.22976 = lb  0.88%  10.7 4340s
end

# Dual Decomposition
function Min_j_fixed_part(j)
    Min_j = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV))
    JuMP.@variable(Min_j, b_j, Bin);
    JuMP.@expression(Min_j, prim_obj_tbMin, c[j]b_j)
    JuMP.@variable(Min_j, 0 ≤ n_j[1:I], Int);
    JuMP.@constraint(Min_j, w ⋅ n_j ≤ W[j]b_j)
    return Min_j
end; function penalty(j, π, n_j::Array{JuMP.VariableRef})
    model = JuMP.owner_model(n_j[1])
    return JuMP.@expression(model, -π ⋅ n_j)
end; function penalty(j, π::Array{JuMP.VariableRef}, n_j)
    model = JuMP.owner_model(π[1])
    return JuMP.@expression(model, -π ⋅ n_j)
end; function mydist(x, x_old)
    d = x - x_old
    return d > 0 ? 2d : d
end; function draw_mydistmax!(bit_vec, Θ, Θ_old)
    max_i, max_v = 0, -Inf
    for (i, undrawn) in enumerate(bit_vec)
        if undrawn
            mydist_i = mydist(Θ[i], Θ_old[i])
            if max_v < mydist_i
                max_v, max_i = mydist_i, i
            end
        end
    end
    if max_i == 0 # cannot draw one due to all-fathomed
        return false, max_i
    else
        bit_vec[max_i] = false # indicating that max_i is drawn this time
        return true, max_i # can be drawn, a valid index
    end
end;

# build subproblems (MIPs)
Min_vec = [Min_j_fixed_part(j) for j = 1:J]; # The inner subproblems
COT = 5; # cut-off tolerance ⚠️ CANNOT by too small
# 0.5^18 * 5 ≈ 1.9073486328125e-5 which is considered small for this case

# build master problem (LP)
an_UB = 543.86951 # a valid ub obtained from Gurobi's trial solve
out = JuMP.Model(() -> Gurobi.Optimizer(GRB_ENV)); # The outer dual problem
JuMP.@variable(out, π[1:I]);
JuMP.@variable(out, θ[1:J]);
JuMP.@expression(out, common, d ⋅ π);
JuMP.@expression(out, outer_obj_tbMax, common + sum(θ));
JuMP.@constraint(out, outer_obj_tbMax ≤ an_UB);
JuMP.@objective(out, Max, outer_obj_tbMax);

# start training
JuMP.set_silent.(Min_vec); JuMP.optimize!.(Min_vec);
JuMP.set_silent.(out); JuMP.optimize!.(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false)
Θ_old = JuMP.value.(θ);
t00 = time()
cnt_outˈs_cut = 0
for COT ∈ 5 * (0.5 .^ (0:18))
    t0 = time()
    re_opt_COT = false # reset to default
    while true
        JuMP.optimize!(out); JuMP.assert_is_solved_and_feasible(out; allow_local = false)
        ub = JuMP.objective_bound(out)
        Π, Θ = JuMP.value.(π), JuMP.value.(θ)
        Θˈs_bitvec = trues(J)
        while true # see if there is any block who can bring us a COT-violating cut
            still_has_undrawn, j = draw_mydistmax!(Θˈs_bitvec, Θ, Θ_old)
            if still_has_undrawn == false
                re_opt_COT = true
                break
            end
            Min_j = Min_vec[j]; prim_obj_tbMin, n_j = Min_j[:prim_obj_tbMin], Min_j[:n_j]
            JuMP.@objective(Min_j, Min, prim_obj_tbMin + penalty(j, Π, n_j))
            JuMP.optimize!(Min_j); JuMP.assert_is_solved_and_feasible(Min_j; allow_local = false)
            ObjVal_j = JuMP.objective_value(Min_j)
            if ObjVal_j + COT < Θ[j] # 🟣 the hallmark of CTPLN's progress
                # Since we implement a single-cut algorithm, we can compress Logging into one line
                @info "#cut = $cnt_outˈs_cut, ub ▶ $ub | cut found at $j with vio = $(Θ[j] - ObjVal_j)"
                Θ_old = Θ # save the RHS, as it is to be updated
                Prim_obj_tbMin, N_j = JuMP.value(prim_obj_tbMin), JuMP.value.(n_j)
                JuMP.@constraint(out, θ[j] ≤ Prim_obj_tbMin + penalty(j, π, N_j))
                cnt_outˈs_cut += 1
                break # Go back to outer model to update Π
            end
        end
        if re_opt_COT
            @info "After $(time() - t0) seconds at COT = $COT, outˈs ub is cut to $ub ◀◀◀"
            break
        end
    end
end
@info "total time elapsed is $(time() - t00)"

If you don’t have a solver, the output from my computer is listed here, it takes around 2533 seconds to converge (the Lagrangian dual bound is found).

Test Results (abridged)

[ Info: #cut = 0, ub ▶ 543.86951 | cut found at 1 with vio = 11.870371549786752
    [ Info: #cut = 1, ub ▶ 543.86951 | cut found at 1 with vio = 5.762397337158175
    [ Info: #cut = 18, ub ▶ 543.8695099999999 | cut found at 1 with vio = 15.191737369851793
    [ Info: #cut = 19, ub ▶ 543.86951 | cut found at 2 with vio = 10.98054942552802
    [ Info: #cut = 20, ub ▶ 543.86951 | cut found at 2 with vio = 12.617016471574242
    [ Info: After 2.697000026702881 seconds at COT = 5.0, outˈs ub is cut to 543.86951 ◀◀◀
    [ Info: #cut = 21, ub ▶ 543.86951 | cut found at 24 with vio = 2.6140881157299267
    [ Info: #cut = 212, ub ▶ 543.8695100000001 | cut found at 193 with vio = 153.31869868564294
    [ Info: #cut = 213, ub ▶ 543.8695100000001 | cut found at 194 with vio = 153.31869868564294
    [ Info: After 1.932999849319458 seconds at COT = 2.5, outˈs ub is cut to 543.86951 ◀◀◀
    [ Info: #cut = 214, ub ▶ 543.86951 | cut found at 83 with vio = 1.493928269336275
    [ Info: #cut = 215, ub ▶ 543.8695100000001 | cut found at 196 with vio = 123.24650581449272
    [ Info: #cut = 422, ub ▶ 543.8695099999995 | cut found at 25 with vio = 3.4328229931116083
    [ Info: #cut = 423, ub ▶ 543.8695099999995 | cut found at 9 with vio = 2.0647643124021444
    [ Info: After 2.759000062942505 seconds at COT = 1.25, outˈs ub is cut to 543.8695099999998 ◀◀◀
    [ Info: #cut = 424, ub ▶ 543.8695099999998 | cut found at 5 with vio = 0.6533606411371813
    [ Info: #cut = 463, ub ▶ 543.8695099999998 | cut found at 390 with vio = 0.7524202592835127
    [ Info: After 26.055999994277954 seconds at COT = 0.625, outˈs ub is cut to 543.869509999999 ◀◀◀
    [ Info: #cut = 464, ub ▶ 543.869509999999 | cut found at 78 with vio = 0.5883608670741545
    [ Info: #cut = 704, ub ▶ 527.1729260387428 | cut found at 385 with vio = 0.5893970978720195
    [ Info: #cut = 705, ub ▶ 526.5968815073952 | cut found at 398 with vio = 0.38940554750680656
    [ Info: #cut = 706, ub ▶ 526.2160343539398 | cut found at 399 with vio = 0.49137861616685163
    [ Info: After 248.89100003242493 seconds at COT = 0.3125, outˈs ub is cut to 525.7264968923939 ◀◀◀
    [ Info: #cut = 707, ub ▶ 525.7264968923939 | cut found at 261 with vio = 0.19378271594975816
    [ Info: #cut = 708, ub ▶ 525.5327141764443 | cut found at 324 with vio = 0.1587827252084647
    [ Info: #cut = 709, ub ▶ 525.374449508871 | cut found at 326 with vio = 0.19334235832094182
    [ Info: #cut = 710, ub ▶ 525.1855748735687 | cut found at 342 with vio = 0.16975234578766685
    [ Info: #cut = 711, ub ▶ 525.0202081432539 | cut found at 378 with vio = 0.2680457571026508
    [ Info: #cut = 712, ub ▶ 524.75764558499 | cut found at 388 with vio = 0.1999731071202997
    [ Info: After 60.257999897003174 seconds at COT = 0.15625, outˈs ub is cut to 524.585856285976 ◀◀◀
    [ Info: #cut = 713, ub ▶ 524.585856285976 | cut found at 231 with vio = 0.09411569109741103
    [ Info: #cut = 714, ub ▶ 524.493510595999 | cut found at 99 with vio = 0.08641177116800397
    [ Info: #cut = 715, ub ▶ 524.4106933166834 | cut found at 213 with vio = 0.12359425878824548
    [ Info: #cut = 716, ub ▶ 524.2871950596319 | cut found at 389 with vio = 0.1551771802770019
    [ Info: After 35.04699993133545 seconds at COT = 0.078125, outˈs ub is cut to 524.1380433375095 ◀◀◀
    [ Info: #cut = 717, ub ▶ 524.1380433375095 | cut found at 121 with vio = 0.07289130946642994
    [ Info: #cut = 718, ub ▶ 524.0664450340456 | cut found at 88 with vio = 0.06109768608013777
    [ Info: #cut = 719, ub ▶ 524.0602434646266 | cut found at 207 with vio = 0.04746792342495898
    [ Info: #cut = 720, ub ▶ 524.0147440349274 | cut found at 131 with vio = 0.04158498831953028
    [ Info: After 25.60800004005432 seconds at COT = 0.0390625, outˈs ub is cut to 523.9865789904918 ◀◀◀
    [ Info: #cut = 721, ub ▶ 523.9865789904918 | cut found at 146 with vio = 0.022078321366601372
    [ Info: #cut = 732, ub ▶ 523.812691201531 | cut found at 41 with vio = 0.019877049312119144
    [ Info: #cut = 733, ub ▶ 523.7952254196871 | cut found at 35 with vio = 0.021413011989790842
    [ Info: After 37.80999994277954 seconds at COT = 0.01953125, outˈs ub is cut to 523.7740964625181 ◀◀◀
    [ Info: #cut = 734, ub ▶ 523.7740964625181 | cut found at 5 with vio = 0.018884466457590898
    [ Info: #cut = 735, ub ▶ 523.76302393272 | cut found at 83 with vio = 0.01023027256361475
    [ Info: #cut = 765, ub ▶ 523.5744781843947 | cut found at 225 with vio = 0.015091169589341202
    [ Info: #cut = 766, ub ▶ 523.5600904160037 | cut found at 349 with vio = 0.014028395570364283
    [ Info: After 136.50100016593933 seconds at COT = 0.009765625, outˈs ub is cut to 523.5461076106526 ◀◀◀
    [ Info: #cut = 767, ub ▶ 523.5461076106526 | cut found at 235 with vio = 0.005194271290245578
    [ Info: #cut = 817, ub ▶ 523.3834419524652 | cut found at 25 with vio = 0.005503791238143374
    [ Info: #cut = 823, ub ▶ 523.3569835347876 | cut found at 72 with vio = 0.005375088765035119
    [ Info: #cut = 824, ub ▶ 523.3526246827648 | cut found at 291 with vio = 0.005578701596306579
    [ Info: After 170.40400004386902 seconds at COT = 0.0048828125, outˈs ub is cut to 523.3481736276733 ◀◀◀
    [ Info: #cut = 825, ub ▶ 523.3481736276733 | cut found at 8 with vio = 0.0029746368053011896
    [ Info: #cut = 826, ub ▶ 523.3468288450199 | cut found at 147 with vio = 0.0037508826174513787
    [ Info: #cut = 852, ub ▶ 523.3090039057765 | cut found at 213 with vio = 0.0026518133953270517
    [ Info: #cut = 853, ub ▶ 523.306819105528 | cut found at 175 with vio = 0.0034458242405134287
    [ Info: After 87.07299995422363 seconds at COT = 0.00244140625, outˈs ub is cut to 523.3034285414248 ◀◀◀
    [ Info: #cut = 854, ub ▶ 523.3034285414248 | cut found at 75 with vio = 0.0019928843004987373
    [ Info: #cut = 855, ub ▶ 523.302636813051 | cut found at 55 with vio = 0.0014690346490231487
    [ Info: #cut = 856, ub ▶ 523.3016203057981 | cut found at 25 with vio = 0.0014558820220373914
    [ Info: #cut = 897, ub ▶ 523.2701752027644 | cut found at 261 with vio = 0.0025875415109232747
    [ Info: #cut = 898, ub ▶ 523.2679630122562 | cut found at 156 with vio = 0.0013831000528373716
    [ Info: After 202.1050000190735 seconds at COT = 0.001220703125, outˈs ub is cut to 523.2670407086156 ◀◀◀
    [ Info: #cut = 899, ub ▶ 523.2670407086156 | cut found at 165 with vio = 0.0008338154695605171
    [ Info: #cut = 900, ub ▶ 523.2666638754026 | cut found at 283 with vio = 0.0006771581259470416
    [ Info: #cut = 930, ub ▶ 523.2544906104197 | cut found at 231 with vio = 0.001032443670632599
    [ Info: #cut = 931, ub ▶ 523.2538459783019 | cut found at 307 with vio = 0.0007227225626918465
    [ Info: After 254.76199984550476 seconds at COT = 0.0006103515625, outˈs ub is cut to 523.2536779260275 ◀◀◀
    [ Info: #cut = 932, ub ▶ 523.2536779260275 | cut found at 388 with vio = 0.0005013221796758005
    [ Info: #cut = 933, ub ▶ 523.2534109472358 | cut found at 390 with vio = 0.000313563626513913
    [ Info: #cut = 986, ub ▶ 523.2446530442485 | cut found at 82 with vio = 0.00035620807272623844
    [ Info: After 299.9390001296997 seconds at COT = 0.00030517578125, outˈs ub is cut to 523.2443384690868 ◀◀◀
    [ Info: #cut = 987, ub ▶ 523.2443384690868 | cut found at 385 with vio = 0.00024342442282865306
    [ Info: #cut = 988, ub ▶ 523.2441904876744 | cut found at 195 with vio = 0.00020671622353529706
    [ Info: #cut = 1018, ub ▶ 523.2418745652013 | cut found at 207 with vio = 0.00022796573633820927
    [ Info: #cut = 1019, ub ▶ 523.241678218329 | cut found at 305 with vio = 0.00018377982955475325
    [ Info: After 204.11999988555908 seconds at COT = 0.000152587890625, outˈs ub is cut to 523.2415159393454 ◀◀◀
    [ Info: #cut = 1020, ub ▶ 523.2415159393454 | cut found at 223 with vio = 0.00013569924776690723
    [ Info: #cut = 1021, ub ▶ 523.2414833821751 | cut found at 399 with vio = 0.00011242151491436259
    [ Info: #cut = 1022, ub ▶ 523.2414351290494 | cut found at 398 with vio = 0.00013316328034196356
    [ Info: #cut = 1023, ub ▶ 523.2414097351362 | cut found at 4 with vio = 0.00012974848384272875
    [ Info: #cut = 1024, ub ▶ 523.2413199660729 | cut found at 9 with vio = 9.582991882495229e-5
    [ Info: #cut = 1025, ub ▶ 523.2412766134983 | cut found at 29 with vio = 8.126213381220992e-5
    [ Info: #cut = 1026, ub ▶ 523.2412603390333 | cut found at 159 with vio = 8.763934125166628e-5
    [ Info: #cut = 1027, ub ▶ 523.2412265620188 | cut found at 192 with vio = 7.995814703898318e-5
    [ Info: #cut = 1028, ub ▶ 523.2412218015062 | cut found at 391 with vio = 9.027611943301928e-5
    [ Info: #cut = 1029, ub ▶ 523.2412032184443 | cut found at 95 with vio = 0.0001203012964747785
    [ Info: #cut = 1030, ub ▶ 523.2411352613677 | cut found at 14 with vio = 8.762553698860032e-5
    [ Info: #cut = 1031, ub ▶ 523.2410932661438 | cut found at 114 with vio = 0.0001241730609136127
    [ Info: #cut = 1032, ub ▶ 523.2410144793715 | cut found at 63 with vio = 0.00010697843825580033
    [ Info: #cut = 1033, ub ▶ 523.2409686003512 | cut found at 194 with vio = 0.00011662782590429743
    [ Info: #cut = 1034, ub ▶ 523.2409033361577 | cut found at 132 with vio = 0.00011963021611627322
    [ Info: #cut = 1035, ub ▶ 523.2408801103109 | cut found at 100 with vio = 0.00010669589253209377
    [ Info: #cut = 1036, ub ▶ 523.240821679687 | cut found at 213 with vio = 8.917986898093755e-5
    [ Info: #cut = 1037, ub ▶ 523.2408216796871 | cut found at 367 with vio = 0.00010616863921097952
    [ Info: #cut = 1038, ub ▶ 523.2408042892193 | cut found at 252 with vio = 7.917423658176936e-5
    [ Info: #cut = 1039, ub ▶ 523.2407807952125 | cut found at 2 with vio = 8.636630156988456e-5
    [ Info: #cut = 1040, ub ▶ 523.2407689690912 | cut found at 132 with vio = 0.00017625938635246197
    [ Info: #cut = 1041, ub ▶ 523.240757521706 | cut found at 78 with vio = 0.00010360397596143045
    [ Info: #cut = 1042, ub ▶ 523.2407142695083 | cut found at 192 with vio = 8.564255785081798e-5
    [ Info: #cut = 1043, ub ▶ 523.2406975906948 | cut found at 283 with vio = 0.00014691020591839354
    [ Info: #cut = 1044, ub ▶ 523.2406492457351 | cut found at 344 with vio = 8.739689572045961e-5
    [ Info: #cut = 1045, ub ▶ 523.2405764631116 | cut found at 41 with vio = 0.00010154481794649728
    [ Info: #cut = 1046, ub ▶ 523.240558025559 | cut found at 15 with vio = 9.966533388383603e-5
    [ Info: #cut = 1047, ub ▶ 523.2405026893395 | cut found at 76 with vio = 8.609617984878248e-5
    [ Info: #cut = 1048, ub ▶ 523.2404591827366 | cut found at 193 with vio = 8.361445968729786e-5
    [ Info: #cut = 1049, ub ▶ 523.240396195796 | cut found at 233 with vio = 8.704243709656279e-5
    [ Info: #cut = 1050, ub ▶ 523.2403729280155 | cut found at 129 with vio = 0.00010199275911948469
    [ Info: #cut = 1051, ub ▶ 523.2403107431636 | cut found at 380 with vio = 9.303311616848386e-5
    [ Info: #cut = 1052, ub ▶ 523.2402931795766 | cut found at 229 with vio = 9.908977232431226e-5
    [ Info: #cut = 1053, ub ▶ 523.2402287695182 | cut found at 142 with vio = 0.0001383911439821195
    [ Info: #cut = 1054, ub ▶ 523.240130093988 | cut found at 186 with vio = 9.563875620077766e-5
    [ Info: #cut = 1055, ub ▶ 523.240102241559 | cut found at 12 with vio = 9.907024139776954e-5
    [ Info: #cut = 1056, ub ▶ 523.2400230223338 | cut found at 112 with vio = 0.000119105138781328
    [ Info: #cut = 1057, ub ▶ 523.2399393903746 | cut found at 263 with vio = 0.00010478953650772116
    [ Info: #cut = 1058, ub ▶ 523.2398727028756 | cut found at 74 with vio = 9.078906427578692e-5
    [ Info: #cut = 1059, ub ▶ 523.2398344548566 | cut found at 53 with vio = 7.845580608045744e-5
    [ Info: #cut = 1060, ub ▶ 523.239814074643 | cut found at 281 with vio = 8.303729833203377e-5
    [ Info: #cut = 1061, ub ▶ 523.2397696787589 | cut found at 382 with vio = 9.85939465454333e-5
    [ Info: #cut = 1062, ub ▶ 523.2397352616018 | cut found at 18 with vio = 8.692458729053776e-5
    [ Info: #cut = 1063, ub ▶ 523.2396679050113 | cut found at 152 with vio = 0.00012598548905173867
    [ Info: #cut = 1064, ub ▶ 523.2395497097251 | cut found at 83 with vio = 8.50125318718753e-5
    [ Info: #cut = 1065, ub ▶ 523.2394880251247 | cut found at 131 with vio = 9.588267572413489e-5
    [ Info: #cut = 1066, ub ▶ 523.2394688947174 | cut found at 390 with vio = 7.989506432271032e-5
    [ Info: #cut = 1067, ub ▶ 523.2394107756855 | cut found at 53 with vio = 7.933976781326901e-5
    [ Info: #cut = 1068, ub ▶ 523.2393785466109 | cut found at 173 with vio = 7.82225073409637e-5
    [ Info: #cut = 1069, ub ▶ 523.2393548357724 | cut found at 83 with vio = 8.260504420753989e-5
    [ Info: #cut = 1070, ub ▶ 523.2393256770007 | cut found at 307 with vio = 0.00010639385603256057
    [ Info: After 264.7590000629425 seconds at COT = 7.62939453125e-5, outˈs ub is cut to 523.2392432822675 ◀◀◀
    [ Info: #cut = 1071, ub ▶ 523.2392432822675 | cut found at 29 with vio = 4.426142106028763e-5
    [ Info: #cut = 1072, ub ▶ 523.239233319602 | cut found at 198 with vio = 4.53107651852136e-5
    [ Info: #cut = 1073, ub ▶ 523.2392184199151 | cut found at 74 with vio = 4.1951644445925584e-5
    [ Info: #cut = 1074, ub ▶ 523.239212243405 | cut found at 131 with vio = 9.430536082088992e-5
    [ Info: #cut = 1075, ub ▶ 523.2391576818933 | cut found at 83 with vio = 7.640760908489419e-5
    [ Info: #cut = 1076, ub ▶ 523.2391372809327 | cut found at 192 with vio = 5.305300963331927e-5
    [ Info: #cut = 1077, ub ▶ 523.2391194848177 | cut found at 41 with vio = 6.379105111642414e-5
    [ Info: #cut = 1078, ub ▶ 523.2390939659858 | cut found at 182 with vio = 8.6231927776792e-5
    [ Info: #cut = 1079, ub ▶ 523.2390218416822 | cut found at 252 with vio = 5.058667591806287e-5
    [ Info: #cut = 1080, ub ▶ 523.2389939147854 | cut found at 25 with vio = 4.436827978215874e-5
    [ Info: #cut = 1081, ub ▶ 523.2389879851349 | cut found at 398 with vio = 4.2191710506456204e-5
    [ Info: #cut = 1082, ub ▶ 523.2389820173807 | cut found at 235 with vio = 4.6073566891347006e-5
    [ Info: #cut = 1083, ub ▶ 523.2389731828639 | cut found at 169 with vio = 5.327052265446941e-5
    [ Info: #cut = 1084, ub ▶ 523.2389701130397 | cut found at 388 with vio = 5.1747868878793124e-5
    [ Info: #cut = 1085, ub ▶ 523.2389616445747 | cut found at 279 with vio = 7.293151089649008e-5
    [ Info: #cut = 1086, ub ▶ 523.238896073237 | cut found at 109 with vio = 4.5594883815780474e-5
    [ Info: #cut = 1087, ub ▶ 523.2388812047641 | cut found at 323 with vio = 3.8357874375893886e-5
    [ Info: #cut = 1088, ub ▶ 523.2388718345184 | cut found at 288 with vio = 5.528659208298148e-5
    [ Info: #cut = 1089, ub ▶ 523.238863606057 | cut found at 55 with vio = 5.005559772963686e-5
    [ Info: #cut = 1090, ub ▶ 523.2388573751626 | cut found at 337 with vio = 6.045039733126867e-5
    [ Info: #cut = 1091, ub ▶ 523.238818751392 | cut found at 382 with vio = 4.069047129129366e-5
    [ Info: #cut = 1092, ub ▶ 523.2387970594327 | cut found at 73 with vio = 6.344569554483925e-5
    [ Info: #cut = 1093, ub ▶ 523.2387551588735 | cut found at 344 with vio = 5.843257174886762e-5
    [ Info: #cut = 1094, ub ▶ 523.2387273709537 | cut found at 92 with vio = 5.0427894647708804e-5
    [ Info: #cut = 1095, ub ▶ 523.23871852673 | cut found at 90 with vio = 5.6684968510301026e-5
    [ Info: #cut = 1096, ub ▶ 523.2386778703632 | cut found at 122 with vio = 5.258248175632474e-5
    [ Info: #cut = 1097, ub ▶ 523.2386465318253 | cut found at 291 with vio = 4.139749354037381e-5
    [ Info: #cut = 1098, ub ▶ 523.238611871967 | cut found at 137 with vio = 5.802865082293884e-5
    [ Info: #cut = 1099, ub ▶ 523.2385899229574 | cut found at 177 with vio = 6.307857307485953e-5
    [ Info: #cut = 1100, ub ▶ 523.2385353178263 | cut found at 4 with vio = 4.419450095838329e-5
    [ Info: #cut = 1101, ub ▶ 523.2385161072827 | cut found at 194 with vio = 4.225805279500783e-5
    [ Info: #cut = 1102, ub ▶ 523.2384895583265 | cut found at 88 with vio = 4.369615911092861e-5
    [ Info: #cut = 1103, ub ▶ 523.2384573857343 | cut found at 389 with vio = 3.876502129707138e-5
    [ Info: #cut = 1104, ub ▶ 523.238449356921 | cut found at 69 with vio = 3.8338808874049946e-5
    [ Info: #cut = 1105, ub ▶ 523.2384436232437 | cut found at 361 with vio = 5.712111977951295e-5
    [ Info: #cut = 1106, ub ▶ 523.2384099591344 | cut found at 122 with vio = 4.989615350092791e-5
    [ Info: #cut = 1107, ub ▶ 523.2383956346908 | cut found at 29 with vio = 5.328234003032506e-5
    [ Info: #cut = 1108, ub ▶ 523.2383883901622 | cut found at 54 with vio = 5.065972250828388e-5
    [ Info: #cut = 1109, ub ▶ 523.2383776129595 | cut found at 155 with vio = 5.171202056286628e-5
    [ Info: #cut = 1110, ub ▶ 523.2383625979144 | cut found at 219 with vio = 3.954330518785021e-5
    [ Info: #cut = 1111, ub ▶ 523.2383474329455 | cut found at 155 with vio = 5.041599881105974e-5
    [ Info: #cut = 1112, ub ▶ 523.2383264961516 | cut found at 193 with vio = 4.017583991489104e-5
    [ Info: #cut = 1113, ub ▶ 523.2383051489453 | cut found at 135 with vio = 3.953838588488523e-5
    [ Info: #cut = 1114, ub ▶ 523.2382968678381 | cut found at 195 with vio = 5.085172113061809e-5
    [ Info: After 263.210000038147 seconds at COT = 3.814697265625e-5, outˈs ub is cut to 523.238256412327 ◀◀◀
    [ Info: #cut = 1115, ub ▶ 523.238256412327 | cut found at 169 with vio = 2.7982802723203193e-5
    [ Info: #cut = 1116, ub ▶ 523.238256412327 | cut found at 155 with vio = 2.68989434304423e-5
    [ Info: #cut = 1117, ub ▶ 523.2382545227116 | cut found at 344 with vio = 2.85705482968579e-5
    [ Info: #cut = 1118, ub ▶ 523.2382485252559 | cut found at 342 with vio = 3.624098613330684e-5
    [ Info: #cut = 1119, ub ▶ 523.2382436354844 | cut found at 380 with vio = 2.460664920808653e-5
    [ Info: #cut = 1120, ub ▶ 523.238239498117 | cut found at 109 with vio = 2.124706960193734e-5
    [ Info: #cut = 1121, ub ▶ 523.2382365177845 | cut found at 12 with vio = 2.9296279871804387e-5
    [ Info: #cut = 1122, ub ▶ 523.2382328553996 | cut found at 288 with vio = 3.517697504573647e-5
    [ Info: #cut = 1123, ub ▶ 523.2382190574542 | cut found at 190 with vio = 2.09685199785703e-5
    [ Info: #cut = 1124, ub ▶ 523.2382148597029 | cut found at 35 with vio = 2.003825775909185e-5
    [ Info: #cut = 1125, ub ▶ 523.2382132242893 | cut found at 370 with vio = 3.436112319843421e-5
    [ Info: #cut = 1126, ub ▶ 523.2382037690924 | cut found at 399 with vio = 3.6549704637955927e-5
    [ Info: #cut = 1127, ub ▶ 523.2381793546394 | cut found at 152 with vio = 2.84965889995048e-5
    [ Info: #cut = 1128, ub ▶ 523.2381773783794 | cut found at 36 with vio = 2.3843521569921755e-5
    [ Info: #cut = 1129, ub ▶ 523.238170913081 | cut found at 137 with vio = 4.7259751448480714e-5
    [ Info: #cut = 1130, ub ▶ 523.2381567225704 | cut found at 360 with vio = 4.461406164724124e-5
    [ Info: #cut = 1131, ub ▶ 523.2381389348593 | cut found at 69 with vio = 2.6382134950542202e-5
    [ Info: #cut = 1132, ub ▶ 523.2381295611558 | cut found at 114 with vio = 3.853129831121738e-5
    [ Info: #cut = 1133, ub ▶ 523.2381173011207 | cut found at 137 with vio = 2.3879520188341274e-5
    [ Info: #cut = 1134, ub ▶ 523.2381101697115 | cut found at 95 with vio = 2.1124604836608718e-5
    [ Info: #cut = 1135, ub ▶ 523.2381054438142 | cut found at 143 with vio = 3.04018486537716e-5
    [ Info: #cut = 1136, ub ▶ 523.2381044716744 | cut found at 100 with vio = 3.0434692670833208e-5
    [ Info: #cut = 1137, ub ▶ 523.2381022966075 | cut found at 201 with vio = 3.2399399676696916e-5
    [ Info: #cut = 1138, ub ▶ 523.2381009933702 | cut found at 130 with vio = 2.2427237003186384e-5
    [ Info: #cut = 1139, ub ▶ 523.2380967384635 | cut found at 137 with vio = 2.9902523812186388e-5
    [ Info: #cut = 1140, ub ▶ 523.2380876004186 | cut found at 361 with vio = 2.8593955334832444e-5
    [ Info: #cut = 1141, ub ▶ 523.2380848079059 | cut found at 117 with vio = 2.0629731041221788e-5
    [ Info: #cut = 1142, ub ▶ 523.2380794741351 | cut found at 28 with vio = 2.5587404309357353e-5
    [ Info: #cut = 1143, ub ▶ 523.2380663752955 | cut found at 22 with vio = 2.298969215247526e-5
    [ Info: #cut = 1144, ub ▶ 523.2380593962182 | cut found at 159 with vio = 2.4487319511212746e-5
    [ Info: #cut = 1145, ub ▶ 523.2380477392574 | cut found at 257 with vio = 2.978624057448087e-5
    [ Info: #cut = 1146, ub ▶ 523.2380441795431 | cut found at 328 with vio = 2.5511451843485133e-5
    [ Info: #cut = 1147, ub ▶ 523.2380410851437 | cut found at 147 with vio = 1.9839962012047074e-5
    [ Info: #cut = 1148, ub ▶ 523.2380395970837 | cut found at 201 with vio = 3.0420364535599376e-5
    [ Info: #cut = 1149, ub ▶ 523.238037371932 | cut found at 243 with vio = 1.9465537836271807e-5
    [ Info: #cut = 1150, ub ▶ 523.2380284092757 | cut found at 375 with vio = 2.412867487899817e-5
    [ Info: After 208.92499995231628 seconds at COT = 1.9073486328125e-5, outˈs ub is cut to 523.2380048199409 ◀◀◀

    julia> @info "total time elapsed is $(time() - t00)"
    [ Info: total time elapsed is 2533.3559999465942

As a comparison, the random shuffle heuristic spends 2843 seconds to finish, which is inferior. Another point of inferior is that it adds more cuts.
The (abridged) result for the shuffle heuristic is

[ Info: #cut = 1399, ub ▶ 523.2382136929144 | cut found at 380 with vio = 1.98959518964148e-5
[ Info: #cut = 1400, ub ▶ 523.2382003885377 | cut found at 78 with vio = 2.1956024736979884e-5
[ Info: After 109.31599998474121 seconds at COT = 1.9073486328125e-5, outˈs ub is cut to 523.2381860879603 ◀◀◀

julia> @info "total time elapsed is $(time() - t00)"
[ Info: total time elapsed is 2842.6630001068115

I spot an error in my previous findmax heuristic. I’ll re-run it.
It’s interesting that it happens to be not-bad.
I indeed intended to write

function mydist(x, x_old)
    d = x - x_old
    return d > 0 ? 2d : -d
end

But it seems that the code I miswrote

function mydist(x, x_old)
    d = x - x_old
    return d > 0 ? 2d : d
end

happens to have better performance. This actually reduces to x - x_old.
I wonder if there is any better heuristics, in cases where we don’t have parallel computing resources.

[stevengj]

@WalterMadelim

You are absolutely right. It’s not at the heart of my global algorithm. It’s a heuristic, to identify a promising block.

@WalterMadelim

I wonder if there is any better heuristics, in cases where we don’t have parallel computing resources.

So, as I understand it, your question has nothing do with the speed of findmax or sort — you have some kind of combinatorial/MIP optimization problem, and you are looking for heuristic algorithms to identify promising candidates.

This seems to have created some confusion in the dialogue above — if you had written down your optimization problem in the first place, people might have been able to point you to some heuristic optimizers to experiment with.

Good luck — there’s a whole literature on this sort of thing, but it’s mostly empirical and the answers are problem-specific as I understand it.

[WalterMadelim]

I’ve already written a (somewhat) satisfying version of algorithm here.

It’s a normal single-threaded program (I don’t have a multicore server or GPU).
Typically other people may claim they do “parallel computing”. (But there is actually a lot of small issues to be noticed in this context.)

Actually if parallel computing is equipped, this is a lesser issue. Since all subproblem blocks can be executed in parallel, and I can pick one according to the violation level fast and return to the master problem. So it would be fine.

I think my single-threaded code is also fine.
I’m going to try some other decomposition algorithms and make comparisons later.
Although that problem might also admit of a Benders decomposition solution method, a block-decomposition spirit will not be embodied. Therefore my investigation is finished.