How many cuts used

[zhenweilin]

Thank you for sharing the code for the proximal bundle method. It demonstrates strong performance.

I have a question that I would like to explore further.

Upon reviewing the IPOPT logs by

set_optimizer_attribute(model, "print_level", 5)

I observed that the maximum number of inequality constraints remained at dimension of variable + 2, even when I set the parameter ncuts_per_iter to very large, like

BM.set_parameter(pm.params, "ncuts_per_iter", 1000)

I would like to know whether there is an upper limit on the number of cuts that can be added.

In fact, it seems that parameter ncuts_per_iter does not work.

[kibaekkim]

Can you try to increase max_age parameter value? If that does not work for your purpose, I can take a look into the code.

[zhenweilin]

Thank you very much for you answer. The parameter max_age does not work.

I use the following scripts:

using Pkg
Pkg.activate("min_max_quadratic")
using JuMP, Ipopt
using LinearAlgebra
using Random  # Add this import
using Ipopt

using SparseArrays
using MAT

# include("./BundleMethod.jl/src/BundleMethod.jl")
using BundleMethod

const BM = BundleMethod

function generate_data(n::Int, num_A::Int, lower_eigenvalue::Float64, upper_eigenvalue::Float64)
    Random.seed!(num_A * n) 
    # generate num_A random matrices A_i in R^{n x n} with symmetric positive definite
    A = Vector{Matrix{Float64}}(undef, num_A)
    B = Vector{Vector{Float64}}(undef, num_A)
    C = Vector{Float64}(undef, num_A)
    max_eig = Vector{Float64}(undef, num_A)
    min_eig = Vector{Float64}(undef, num_A)
    for i in 1:num_A
        A_i = randn(n, n)
        A_i = A_i * A_i'
        A_i .= A_i ./ norm(A_i, Inf)
        ## generate a random eigenvalue between lower_eigenvalue and upper_eigenvalue with constant difference
        eig_val = [lower_eigenvalue + (i - 1) * (upper_eigenvalue - lower_eigenvalue) / (n - 1) for i = 1:n]
        eig_vec = eigen(A_i)
        A_i = eig_vec.vectors * diagm(eig_val) * eig_vec.vectors'
        A[i] = A_i
        if n < 5
            print("A_i: ", A_i, "\n")
        end
        B[i] = vec(randn(n, 1))
        if n < 5
            print("B[i]: ", B[i], "\n")
        end
    end
    C = randn(num_A)
    return A, B, C
end


# Randomly generate problem data
Random.seed!(1)
n = 30
m = 100
A, B, C = generate_data(n, m, 0.1, 1.0)
N = 1


function evaluate_f(x::Vector{Float64})
    fvals = zeros(N)
    grads = Dict{Int,SparseVector}()
    fObj = -Inf;
    index = 0;
    for i = 1:m
        obj = 0.5 * x' * A[i] * x + B[i]' * x + C[i]
        if obj > fObj
            fObj = obj
            index = i   
        end
    end
    if index == 0
        error("No valid index found, obj: ", fObj, " x: ", x)
    else
        grad = A[index] * x + B[index]
    end
    fvals[1] = fObj
    grads[1] = grad
    return fvals, grads
end

# This initializes the proximal bundle method with required arguments.
pm = BM.ProximalMethod(n, N, evaluate_f)

# Set optimization solver to the internal JuMP.Model
model = BM.get_jump_model(pm)
set_optimizer(model, Ipopt.Optimizer)
set_optimizer_attribute(model, "print_level", 5)

# We overwrite the function to have column bounds.
function BM.add_variables!(method::BM.ProximalMethod)
    bundle = BM.get_model(method)
    model = BM.get_model(bundle)
    @variable(model, -1000000 <= x[i=1:bundle.n] <= 1000000)
    @variable(model, θ[j=1:bundle.N])
end

# This builds the bundle model.
BM.build_bundle_model!(pm)
pm.params.ϵ_s = 1e-6
BM.set_parameter(pm.params, "ncuts_per_iter", 20)
BM.set_parameter(pm.params, "max_age", 20)
# This runs the bundle method.
BM.run!(pm)

Based the log, it seems that the number of cuts has been increased to about 120. (from Total number of inequality constraints...............: 125)

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:     3875
Number of nonzeros in Lagrangian Hessian.............:       30

Total number of variables............................:       31
                     variables with only lower bounds:        0
                variables with lower and upper bounds:       30
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
**Total number of inequality constraints...............:      125**
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:      125

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
   0  0.0000000e+00 0.00e+00 9.98e-01  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
   1  1.0978513e-01 0.00e+00 1.41e-01  -1.0 1.10e-01    -  8.76e-01 1.00e+00f  1
   2  8.4120323e-01 0.00e+00 5.07e-01  -1.0 8.05e-01    -  3.80e-01 1.00e+00f  1
   3  1.8612334e+00 0.00e+00 3.39e-01  -1.0 1.28e+00    -  4.84e-01 1.00e+00f  1
   4  3.4509310e+00 0.00e+00 1.11e-02  -1.0 2.37e+00    -  8.26e-01 1.00e+00f  1
   5  7.8859439e+00 0.00e+00 5.35e-02  -1.0 9.32e+00    -  4.01e-01 1.00e+00f  1
   6  1.0514940e+01 0.00e+00 8.46e-05  -1.0 5.86e+00    -  9.90e-01 1.00e+00f  1
   7  1.1918913e+01 0.00e+00 2.26e-05  -1.0 2.92e+00    -  9.92e-01 1.00e+00f  1
   8  2.5952988e+00 0.00e+00 4.88e-05  -1.7 1.38e+01    -  1.00e+00 9.86e-01f  1
   9  6.8014895e-01 0.00e+00 4.10e-04  -2.5 5.95e+00    -  1.00e+00 8.97e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  10  1.0272658e-01 0.00e+00 3.86e-04  -3.8 2.74e+00    -  1.00e+00 9.21e-01f  1
  11  2.5676232e-02 0.00e+00 1.50e-09  -3.8 7.95e-01    -  1.00e+00 1.00e+00f  1
  12  3.9501620e-03 0.00e+00 2.09e-04  -5.7 2.45e-01    -  1.00e+00 9.05e-01f  1
  13  9.5018164e-04 0.00e+00 1.29e-03  -5.7 1.32e-01    -  1.00e+00 8.56e-01f  1
  14  2.3497957e-04 0.00e+00 1.84e-11  -5.7 7.38e-02    -  1.00e+00 1.00e+00f  1
  15  3.9883509e-05 0.00e+00 2.15e-03  -8.6 2.89e-02    -  1.00e+00 7.72e-01f  1
  16  2.2546846e-05 0.00e+00 2.12e-02  -8.6 1.32e-02    -  9.65e-01 3.71e-01f  1
  17  5.6460187e-06 0.00e+00 1.78e-02  -8.6 9.00e-03    -  8.90e-01 6.36e-01f  1
  18  8.4526850e-08 0.00e+00 1.14e-02  -8.6 3.33e-03    -  8.67e-01 6.85e-01f  1
  19 -1.4256899e-06 0.00e+00 4.01e-03  -8.6 1.00e-03    -  1.00e+00 7.38e-01f  1
iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
  20 -1.6546975e-06 0.00e+00 2.36e-03  -8.6 1.89e-04    -  1.00e+00 6.98e-01f  1
  21 -1.6962260e-06 0.00e+00 1.31e-04  -8.6 3.09e-05    -  1.00e+00 9.49e-01f  1
  22 -1.6920905e-06 0.00e+00 2.51e-14  -8.6 4.94e-06    -  1.00e+00 1.00e+00f  1

Number of Iterations....: 22

                                   (scaled)                 (unscaled)
Objective...............:  -1.6920905291666035e-06   -1.6920905291666035e-06
Dual infeasibility......:   2.5063284780912909e-14    2.5063284780912909e-14
Constraint violation....:   0.0000000000000000e+00    0.0000000000000000e+00
Variable bound violation:   0.0000000000000000e+00    0.0000000000000000e+00
Complementarity.........:   5.6698030980209089e-09    5.6698030980209089e-09
Overall NLP error.......:   5.6698030980209089e-09    5.6698030980209089e-09


Number of objective function evaluations             = 23
Number of objective gradient evaluations             = 23
Number of equality constraint evaluations            = 0
Number of inequality constraint evaluations          = 23
Number of equality constraint Jacobian evaluations   = 0
Number of inequality constraint Jacobian evaluations = 1
Number of Lagrangian Hessian evaluations             = 1
Total seconds in IPOPT                               = 0.010

I hope this is a limited memory algorithm, so it can handle large scale problem.

[kibaekkim]

Was the number of cuts increased? If you set max_age very large, the algorithm would not remove any cuts generated. Do you have logs from the bundle method?

[zhenweilin]

The log is shown like

Iter  246: ncols    31, nrows   123, fx0 +1.065055e+00, fy +1.065055e+00, m +1.065053e+00, v -2.451063e-06, u 1.000000e-01, i +1, master time    1.6s, eval time    0.3s, time    4.1s
Iter  247: ncols    31, nrows   124, fx0 +1.065055e+00, fy +1.065055e+00, m +1.065053e+00, v -2.420010e-06, u 1.000000e-01, i -1, master time    1.6s, eval time    0.3s, time    4.1s
Iter  248: ncols    31, nrows   125, fx0 +1.065055e+00, fy +1.065055e+00, m +1.065053e+00, v -1.786314e-06, u 1.000000e-01, i +1, master time    1.6s, eval time    0.3s, time    4.1s
TERMINATION: Optimal: v = -1.7863135307410635e-6

You mean nrows is the number of cuts?

My goal is to set the number of cuts at level 20.

[kibaekkim]

You can directly control the number of cuts to maintain by re-defining the following function (for the proximal bundle):

BundleMethod.jl/src/Proximal.jl

Lines 260 to 281 in e87ad0a

	function purge_bundles!(method::ProximalMethod)
	bundle = get_model(method)
	model = get_model(bundle)
	ncuts_removed = 0
	ncols = JuMP.num_variables(model)
	ncuts = length(method.cuts)
	for (ref, cut) in method.cuts
	if ncuts - ncuts_removed <= ncols
	break
	end
	if cut["age"] >= method.params.max_age
	JuMP.delete(model, ref)
	delete!(method.cuts, ref)
	ncuts_removed += 1
	end
	end
	if ncuts_removed > 0
	if (method.params.print_output)
	@printf("Removed %d inactive cuts.\n", ncuts_removed)
	end
	end
	end

You can do so similarly for the trust-region bundle.

[zhenweilin]

Thank you very much.

function purge_bundles!(method::ProximalMethod)
    bundle = get_model(method)
    model = get_model(bundle)
    ncuts_removed = 0
    ncols = JuMP.num_variables(model)
    ncuts = length(method.cuts)
    sorted_cuts = sort(collect(method.cuts), by = x -> x.second["age"])
    counter = 0
    for (ref, cut) in sorted_cuts
        counter += 1
        if cut["age"] >= method.params.max_age || counter > method.params.ncuts_per_iter
            JuMP.delete(model, ref)
            delete!(method.cuts, ref)
            ncuts_removed += 1
        end
    end
    if ncuts_removed > 0
        if (method.params.print_output)
            @printf("Removed %d inactive cuts.\n", ncuts_removed)
        end
    end
end

The code provided above is used to control the number of cuts. It is designed to meet the intended objective of removing the oldest cuts while maintaining a fixed total number of cuts.

[kibaekkim]

if cut["age"] >= method.params.max_age || counter > method.params.ncuts_per_iter

method.params.ncuts_per_iter is the number of cuts to be added every iteration. Instead, you may need additional parameter for your purpose.

[zhenweilin]

Thank you for the reminder. I will introduce a new parameter to regulate the number of cuts.