Column generation requiring a lot of memory?

Specific Domains Optimization (Mathematical)
1

Hi

I’m solving a (rather large) LP through two approaches. One involves solving the full LP. The second consists in doing column generation, gradually adding variables through @variable and set_normalized_coefficient.

I’m quite surprised by the resulting memory usage (measured through @allocated). For instance:

  • a problem for which the full LP has 22k constraints and 1.3m variables requires 2.3GB of memory.
  • with column generation, the final LP still has 22k constraints, but only 0.2m variables. Yet, the latter algorithm consumes 11GB of memory.

Is it an expected output? I tried an alternative approach, where I get the backend model and use instead MathOptInterface.modify, but this did not lead to improvements.

Best,

Michaël.

2

Hi @mposs, note that @allocated measures total allocations, not the peak requirement. Your code is probably creating and destroying lots of intermediate objects.

A better question is: is the column generation faster to solve?

If not, could you post a reproducible example of your code? There might be some easy improvements to make.

3

It is not. So I was hoping at least the memory usage would be interesting, but not even that.

The code is quite large, but maybe providing only the column generation loop might lead to interesting insights (as this is clearly the critical piece)? As you might see below, the pricing involves an enumeration, as the number of variables is polynomially bounded.

        while var_added
            niterations += 1
            var_added = false
            memory_separation = @allocated begin
                for s in 1:4, r in R[s]
                    σ_val[s,r] = dual(σ[s,r])
                    μ_val[s,r] = dual.(μ[s,r])
                end
                for s in 1:4, r in R[s], e in E_setdiff
                    if μ_val[s,r][e[1]]*sqrtΣ[e[2]] - μ_val[s,r][e[2]]*sqrtΣ[e[1]] > γ_mod⁺[e[1],e[2]]*σ_val[s,r] + ϵ
                        β⁺[s,r][e] = @variable(model, lower_bound = 0.0, base_name = "β⁺[$s,$r][$e]")
                        set_normalized_coefficient(μ[s,r][e[1]], β⁺[s,r][e], sqrtΣ[e[2]])
                        set_normalized_coefficient(μ[s,r][e[2]], β⁺[s,r][e], -sqrtΣ[e[1]])
                        set_normalized_coefficient(σ[s,r], β⁺[s,r][e], -γ_mod⁺[e[1],e[2]])
                        var_added = true
                        nvars += 1
                    end
                    if μ_val[s,r][e[2]]*sqrtΣ[e[1]] - μ_val[s,r][e[1]]*sqrtΣ[e[2]] > γ_mod⁻[e[1],e[2]]*σ_val[s,r] + ϵ
                        β⁻[s,r][e] = @variable(model, lower_bound = 0.0, base_name = "β⁻[$s,$r][$e]")
                        set_normalized_coefficient(μ[s,r][e[2]], β⁻[s,r][e], sqrtΣ[e[1]])
                        set_normalized_coefficient(μ[s,r][e[1]], β⁻[s,r][e], -sqrtΣ[e[2]])
                        set_normalized_coefficient(σ[s,r], β⁻[s,r][e], -γ_mod⁻[e[1],e[2]])
                        var_added = true
                        nvars += 1
                    end
                end
            end
            if var_added 
                memory_optimization = @allocated optimize!(model)
            end
        end
4

So I was hoping at least the memory usage would be interesting, but not even that.

@allocated is misleading. It does not measure peak memory usage. It measures the total allocations in Julia. It doesn’t measure the memory usage of the solver, for example.

Have you tried profiling where the time is spent? How much time is spent solving the master problem? How much time is spent computing new columns? How much time is spent adding new columns to the master?

5

Apparently, roughly half of the time is in the separation and half of the time is in re-optimizing the master. I have an additional question by the way. Here is the previous code adapted to profiling with TimerOutputs:

            @timeit to "separation" begin
                for s in 1:4, r in R[s]
                    @timeit to "getting duals" σ_val[s,r] = dual(σ[s,r])
                    @timeit to "getting duals" μ_val[s,r] = dual.(μ[s,r])
                end
                for s in 1:4, r in R[s], e in E_setdiff
                    if μ_val[s,r][e[1]]*sqrtΣ[e[2]] - μ_val[s,r][e[2]]*sqrtΣ[e[1]] > γ_mod⁺[e[1],e[2]]*σ_val[s,r] + ϵ
                        @timeit to "creating var" β⁺[s,r][e] = @variable(model, lower_bound = 0.0, base_name = "β⁺[$s,$r][$e]")
                        @timeit to "set coef" set_normalized_coefficient(μ[s,r][e[1]], β⁺[s,r][e], sqrtΣ[e[2]])
                        @timeit to "set coef" set_normalized_coefficient(μ[s,r][e[2]], β⁺[s,r][e], -sqrtΣ[e[1]])
                        @timeit to "set coef" set_normalized_coefficient(σ[s,r], β⁺[s,r][e], -γ_mod⁺[e[1],e[2]])
                        var_added = true
                        nvars += 1
                    end
                    if μ_val[s,r][e[2]]*sqrtΣ[e[1]] - μ_val[s,r][e[1]]*sqrtΣ[e[2]] > γ_mod⁻[e[1],e[2]]*σ_val[s,r] + ϵ
                        @timeit to "creating var" β⁻[s,r][e] = @variable(model, lower_bound = 0.0, base_name = "β⁻[$s,$r][$e]")
                        @timeit to "set coef" set_normalized_coefficient(μ[s,r][e[2]], β⁻[s,r][e], sqrtΣ[e[1]])
                        @timeit to "set coef" set_normalized_coefficient(μ[s,r][e[1]], β⁻[s,r][e], -sqrtΣ[e[2]])
                        @timeit to "set coef" set_normalized_coefficient(σ[s,r], β⁻[s,r][e], -γ_mod⁻[e[1],e[2]])
                        var_added = true
                        nvars += 1
                    end
                end
            end
            if var_added 
                @timeit to "re-optimization" optimize!(model)
            end

and here come the results:

─────────────────────────────────────────────────────────────────────────────────
                                        Time                    Allocations      
                               ───────────────────────   ────────────────────────
       Tot / % measured:            8.75s /  98.0%            774MiB /  96.4%    

Section                ncalls     time    %tot     avg     alloc    %tot      avg
─────────────────────────────────────────────────────────────────────────────────
re-optimization             6    4.79s   55.9%   799ms      864B    0.0%     144B
separation                  7    3.21s   37.4%   458ms    471MiB   63.2%  67.3MiB
  set coef              90.2k   93.9ms    1.1%  1.04μs   23.1MiB    3.1%     269B
  creating var          30.1k   60.5ms    0.7%  2.01μs   47.0MiB    6.3%  1.60KiB
  getting duals         3.32k   37.4ms    0.4%  11.3μs   19.7MiB    2.6%  6.07KiB
first optimization          1    425ms    5.0%   425ms   54.9MiB    7.4%  54.9MiB
creating first model        1    150ms    1.7%   150ms    220MiB   29.5%   220MiB
─────────────────────────────────────────────────────────────────────────────────

Apparently, most of the time is spent testing the conditions. I first thought this was due to using Dictionnaries for σ_val and μ_val, but turning them to multi-dimensional arrays does not seem to lead to any improvement.