Inner product grammar is neater than sum(index) in JuMP modeling, But triggers Warning?

Specific Domains Optimization (Mathematical)
1

Continuing the discussion from JuMP slow model generation:

I figured out a novel neat way to write a particular kind of constraint, however it trigger warnings like the above post did. Can we somehow improve this?

see the following code

############## An motivating example ##############
import JuMP
I = 4 # change 4 to 20004 to see the Warning
M_vector = [rand(2, 3) for _ in 1:I];
M_mean = rand(2, 3);
model = JuMP.Model();
# The following line create a probability decision vector
JuMP.@variable(model, p[1:I] >= 0); JuMP.@constraint(model, sum(p) == 1);
# The following are meant to be 4 equivalent ways to write a Mean constraint
# 1st: usable, fast
JuMP.@constraint(model, sum(p[i] * M_vector[i] for i in 1:I) .== M_mean) # correct, but cumbersome!
# 2nd: unusable currently
import LinearAlgebra: dot as ip
JuMP.@constraint(model, ip(p, M_vector) .== M_mean) # elegantly neat, but LinearAlgebra doesn't support (regrettably)
# 3rd (my novel idea): usable, but currently slow with Warning when I is large
function my_ip(p, vector) return mapreduce(prod, +, zip(p, vector)) end
JuMP.@constraint(model, my_ip(p, M_vector) .== M_mean)
# 4th: If we want to piggyback on ip (inner product) like the 2nd method above, we could write
foreach([[M_vector[i][e] for i in eachindex(p)] for e in eachindex(M_mean)], M_mean) do l, r
    JuMP.@constraint(model, ip(p, l) == r)
end 

# if I = 20004, the 1st and 4th method are fast, the 3rd method is slow
# Can we somehow improve the 3rd method?

############## Another example ##############
I = 4
M_vector = rand(2, 3, I);
M_mean = rand(2, 3);
model = JuMP.Model();
JuMP.@variable(model, p[1:I] >= 0); JuMP.@constraint(model, sum(p) == 1);
# instead of writing 
JuMP.@constraint(model, sum(p[i] * M_vector[:, :, i] for i in 1:I) .== M_mean)
# we could write
JuMP.@constraint(model, my_ip(p, eachslice(M_vector; dims = 3)) .== M_mean)
2

An unexpected new finding is that This Warning persists, even if I opt not to use mapreduce. (I only used LinearAlgebra.dot)

The whole code have 310 lines (copy paste to REPL if you want to try). The Warning should appear within 2 minutes.

import JuMP
import Gurobi
import Random
import Distributions: Uniform as UDistri
import LinearAlgebra: dot as ip
import LinearAlgebra: norm as norm

const UPDTOL, argZ_TimeLimit = 0.01, 8.0
GRB_ENV = Gurobi.Env()
Random.seed!(3)
function GurobiDirectModel() return JuMP.direct_model(Gurobi.Optimizer(GRB_ENV)) end
function jvr6(x) return round(JuMP.value(x); digits = 6) end
function jdr6(x) return round( JuMP.dual(x); digits = 6) end
function optimise(m) return (JuMP.optimize!(m); JuMP.termination_status(m)) end
# function my_ip(x::AbstractVector{JuMP.VariableRef}, y)
#     JuMP.@expression(JuMP.owner_model(first(x)), sum(xi * yi for (xi, yi) in zip(x, y)))
# end
function opt_ass_opt(m, str)
    sta = optimise(m)
    sta == JuMP.OPTIMAL || error(str * ": " * string(sta))
end
function opt_ass_optortime(m, str)
    sta = optimise(m)
    sta == JuMP.OPTIMAL && return
    sta != JuMP.TIME_LIMIT && error(str * ": " * string(sta))
    JuMP.has_values(m) || error(str * ": no primal solution")
end
function brcs(v) return ones(T) * transpose(v) end
function gen_load_pattern(fLM)
    fL = length(fLM)
    fLP = rand(UDistri(0.7 * fLM[1], 1.4 * fLM[1]), T, 1)
    for i in 2:fL
        fLP = [fLP rand(UDistri(0.7 * fLM[i], 1.4 * fLM[i]), T)]
    end
    return fLP
end

begin # base case: such that the first stage solution is to turn on all the generators
    G = 2;
    UT = DT = 3;
    PI = [0.45, 0.375, 0.5];
    PS = [4.5, 4, 5.5];
    ZS = [0, 0,   1];
    CST = [0.63, 0.60, 0.72];
    CG = [0.8, 0.68, 0.72];
    CL_BASE = [1.6, 1.3776, 1.4442];
    NW = [2, 3];
    W = 2;
    NL = [4, 5, 6];
    L = 3;
    LM = [4.0, 3.5, 3.0];
    T = 24;
    (CL = brcs(CL_BASE); CL[end, :] *= 2);
end
MZ = gen_load_pattern(LM)
MY, yM, NY = let
    YABSMAX = [1.6, 2.3]
    Pr = UDistri.(0.2 * rand(W), YABSMAX / 1.1)
    Ybaseline = [rand(Pr[w]) for t in 1:T, w in 1:W]
    YMIN, YMAX = 0.9 * Ybaseline, 1.1 * Ybaseline # used to construct vertices
    all0() = Bool.(zero(YMIN))
    vertexList = [all0() for _ in eachindex(YMIN)]
    for i in eachindex(YMIN)
        vertexList[i][i] = true
    end
    pushfirst!(vertexList, all0())
    for i in 1:1+length(YMIN)
        push!(vertexList, .!(vertexList[i]))
    end
    ref = vertexList[1]
    yM = [(ref[i,j] ? YMAX[i,j] : YMIN[i,j]) for i in eachindex(eachrow(YMIN)), j in eachindex(eachcol(YMIN))]
    yM = reshape(yM, (T, W, 1))
    popfirst!(vertexList)
    for ref in vertexList
        yM = cat(yM, [(ref[i,j] ? YMAX[i,j] : YMIN[i,j]) for i in eachindex(eachrow(YMIN)), j in eachindex(eachcol(YMIN))]; dims = 3)
    end # âś… after this line, yM is already gened
    NY = size(yM, 3)
    PrVec = rand(UDistri(.4, .6), NY)
    PrVec = PrVec / sum(PrVec)
    MY = sum(yM[:, :, i] * PrVec[i] for i in 1:NY)
    @assert all(MY .< YMAX) && all(MY .> YMIN)
    function is_strict_in(MY, yM)
        NY = size(yM, 3)
        ø = GurobiDirectModel()
        JuMP.set_silent(ø)
        JuMP.@variable(ø, c[1:NY] >= 1e-5)
        JuMP.@constraint(ø, sum(c) == 1)
        JuMP.@constraint(ø, sum(yM[:, :, i] * c[i] for i in 1:NY) .== MY)
        status = optimise(ø)
        status == JuMP.OPTIMAL && return true
        return false
    end
    @assert is_strict_in(MY, yM)
    MY, yM, size(yM)[3]
end;
Z_list = [MZ]; # must involve MZ to ensure feasibility
function recruit(Z_list, Z_new)
    for Z in Random.shuffle(Z_list)
        norm(Z .- Z_new, Inf) < 5e-5 && return false
    end
    push!(Z_list, Z_new); return true
end

trh = GurobiDirectModel(); JuMP.set_objective_sense(trh, JuMP.MIN_SENSE); JuMP.@variable(trh, oLh); JuMP.set_silent(trh);
JuMP.@variable(trh, trh_u[t = 1:T, g = 1:G+1], Bin);
JuMP.@variable(trh, trh_v[t = 1:T, g = 1:G+1], Bin);
JuMP.@variable(trh, trh_x[t = 1:T, g = 1:G+1], Bin);
let # 1st stage feasible region
    JuMP.@constraint(trh, trh_x .- [transpose(ZS); trh_x[1:T-1, :]] .== trh_u .- trh_v)
    JuMP.@constraint(trh, [g = 1:G+1, t = 1:T-UT+1], sum(trh_x[i, g] for i in t:t+UT-1) >= UT * trh_u[t, g])
    JuMP.@constraint(trh, [g = 1:G+1, t = T-UT+1:T], sum(trh_x[i, g] - trh_u[t, g] for i in t:T) >= 0)
    JuMP.@constraint(trh, [g = 1:G+1, t = 1:T-DT+1], sum(1 - trh_x[i, g] for i in t:t+DT-1) >= DT * trh_v[t, g])
    JuMP.@constraint(trh, [g = 1:G+1, t = T-DT+1:T], sum(1 - trh_x[i, g] - trh_v[t, g] for i in t:T) >= 0)
end;
u, v, x = let
    JuMP.set_objective_function(trh, ip(brcs(CST), trh_u)) # to derive a feas (u, v, x)
    opt_ass_opt(trh, "trh")
    u = jvr6.(trh_u);
    v = jvr6.(trh_v);
    x = jvr6.(trh_x);
    u, v, x
end

fx2Z = GurobiDirectModel(); JuMP.set_objective_sense(fx2Z, JuMP.MAX_SENSE); JuMP.set_silent(fx2Z);
f_c_x, f_c_Y, f_c_Z = let
    JuMP.@variable(fx2Z, Dbal[t = 1:T]);
    JuMP.@variable(fx2Z, Dvp[t = 1:T, w = 1:W] >= 0);
    JuMP.@variable(fx2Z, Dzt[t = 1:T, l = 1:L] >= 0);
    JuMP.@variable(fx2Z, Dvr[t = 1:T, g = 1:G+1] >= 0);
    JuMP.@variable(fx2Z, Dpi[t = 1:T, g = 1:G+1] >= 0);
    JuMP.@variable(fx2Z, Dps[t = 1:T, g = 1:G+1] >= 0);
    JuMP.@constraint(fx2Z, [t = 1:T, g = 1:G+1], Dps[t, g] - Dpi[t, g] - Dvr[t, g] + CG[g] == 0);
    JuMP.@constraint(fx2Z, [t = 1:T, g = 1:G+1], Dbal[t] - CG[g] + Dvr[t, g] >= 0);
    JuMP.@constraint(fx2Z, [t = 1:T, w = 1:W], Dvp[t, w] + Dbal[t] >= 0);
    JuMP.@constraint(fx2Z, [t = 1:T, l = 1:L], Dzt[t, l] - CL[t, l] - Dbal[t] >= 0);
    JuMP.@expression(fx2Z, f_c_x, brcs(PI) .* Dpi .- brcs(PS) .* Dps); # in x2
    JuMP.@expression(fx2Z, f_c_Y, -Dvp); # in x2
    JuMP.@expression(fx2Z, f_c_Z, CL .- Dzt)
    f_c_x, f_c_Y, f_c_Z
end;
function set_f_objective(x, Y, Z)
    JuMP.set_objective_function(fx2Z,
          ip(f_c_x, x)
        + ip(f_c_Y, Y)
        + ip(f_c_Z, Z)
    )
end
function get_f_value(x, Y, Z)
    set_f_objective(x, Y, Z)
    opt_ass_opt(fx2Z, "fx2Z")
    JuMP.objective_value(fx2Z)
end
function get_f_cut(x, Y, Z)
    set_f_objective(x, Y, Z)
    opt_ass_opt(fx2Z, "fx2Z")
    cn = ip(jvr6.(f_c_Z), Z) + ip(jvr6.(f_c_Y), Y)
    px = jvr6.(f_c_x)
    return cn, px
end
function add_cut_for_oLh(cn, px)
    JuMP.@constraint(trh, oLh >= cn + ip(px, trh_x))
end
function get_cut_oLh_coeff(x, Q, P)
    cn, px = 0., zero(x)
    for j in 1:Q[1]
        for i in 1:P[1]
            p = Q[2][j] * P[2][i]
            this_cn, this_px = get_f_cut(x, Q[3][:, :, j], P[3][i])
            cn += p * this_cn
            px += p * this_px
        end
    end
    return cn, px
end
function add_cut_for_oLh(x, Q, P, oLh_check)
    cn, px = get_cut_oLh_coeff(x, Q, P)
    if oLh_check < cn + ip(px, x) - UPDTOL
        add_cut_for_oLh(cn, px)
        return true
    end
    return false
end

argZ = GurobiDirectModel(); JuMP.set_objective_sense(argZ, JuMP.MAX_SENSE); JuMP.set_attribute(argZ, "TimeLimit", argZ_TimeLimit); JuMP.set_silent(argZ);
JuMP.@variable(argZ, argZ_Z[t = 1:T, l = 1:L] == MZ[t, l]); # TBD
argZ_c_x, argZ_c_Y, argZ_c_b2, argZ_cn = let
    JuMP.@variable(argZ, Dbal[t = 1:T]);
    JuMP.@variable(argZ, Dvp[t = 1:T, w = 1:W] >= 0);
    JuMP.@variable(argZ, Dzt[t = 1:T, l = 1:L] >= 0);
    JuMP.@variable(argZ, Dvr[t = 1:T, g = 1:G+1] >= 0);
    JuMP.@variable(argZ, Dpi[t = 1:T, g = 1:G+1] >= 0);
    JuMP.@variable(argZ, Dps[t = 1:T, g = 1:G+1] >= 0);
    JuMP.@constraint(argZ, [t = 1:T, g = 1:G+1], Dps[t, g] - Dpi[t, g] - Dvr[t, g] + CG[g] == 0);
    JuMP.@constraint(argZ, [t = 1:T, g = 1:G+1], Dbal[t] - CG[g] + Dvr[t, g] >= 0);
    JuMP.@constraint(argZ, [t = 1:T, w = 1:W], Dvp[t, w] + Dbal[t] >= 0);
    JuMP.@constraint(argZ, [t = 1:T, l = 1:L], Dzt[t, l] - CL[t, l] - Dbal[t] >= 0);
    JuMP.@expression(argZ, argZ_c_x, brcs(PI) .* Dpi .- brcs(PS) .* Dps); # in x2
    JuMP.@expression(argZ, argZ_c_Y, -Dvp); # in x2
    JuMP.@expression(argZ, argZ_c_b2, -argZ_Z); # is beta2
    JuMP.@expression(argZ, argZ_cn, ip(CL .- Dzt, argZ_Z)); # the rest
    argZ_c_x, argZ_c_Y, argZ_c_b2, argZ_cn
end;
function set_argZ_objective(x, Y, b2) # (x2, b2)
    JuMP.set_objective_function(argZ, argZ_cn
        + ip(argZ_c_x, x)
        + ip(argZ_c_Y, Y)
        + ip(argZ_c_b2, b2)
    )
end
tr2 = GurobiDirectModel(); JuMP.@variable(tr2, oΛ2); JuMP.set_silent(tr2); # ✅
JuMP.@variable(tr2, tr2_b2[eachindex(eachrow(MZ)), eachindex(eachcol(MZ))]);
JuMP.@variable(tr2, tr2_x[t = 1:T, g = 1:G+1] == x[t, g]);
JuMP.@variable(tr2, tr2_Y[t = 1:T, w = 1:W] == MY[t, w]);
function add_cut_for_oΛ2(cn, px, pY, pb)
    JuMP.@constraint(tr2, oΛ2 >= cn + ip(px, tr2_x) + ip(pY, tr2_Y) + ip(pb, tr2_b2)) 
end
let
    set_argZ_objective(x, (local Y = MY), (local b2 = zero(MZ)));
    opt_ass_optortime(argZ, "argZ")
    cn, px, pY = JuMP.value(argZ_cn), jvr6.(argZ_c_x), jvr6.(argZ_c_Y)
    add_cut_for_oΛ2(cn, px, pY, (local pb2 = -JuMP.fix_value.(argZ_Z)))
    JuMP.@objective(tr2, Min, ip(MZ, tr2_b2) + oΛ2); # ❄️ stable obj
    try
        opt_ass_opt(tr2, "tr2")
    finally # finish establishing Z's support
        JuMP.unfix.(argZ_Z)
        (JuMP.set_lower_bound.(argZ_Z, 0.9 * MZ); JuMP.set_upper_bound.(argZ_Z, 1.1 * MZ))
        JuMP.@variable(argZ, argZ_a[eachindex(eachrow(MZ)), eachindex(eachcol(MZ))])
        JuMP.@constraint(argZ, argZ_a .>= argZ_Z .- MZ)
        JuMP.@constraint(argZ, argZ_a .>= MZ .- argZ_Z)
        JuMP.@constraint(argZ, [t = 1:T], sum(argZ_a[t, :] ./ (0.1 * MZ[t, :])) <= L/3)
    end
end
function psi_via_19(x, i)
    JuMP.fix.(tr2_x, x)
    JuMP.fix.(tr2_Y, yM[:, :, i])
    opt_ass_opt(tr2, "tr2")
    JuMP.objective_value(tr2)
end
function get_proper_Q(x)
    JuMP.set_objective_function(fr9, ip(fr9_p, [psi_via_19(x, i) for i in 1:NY]));
    opt_ass_opt(fr9, "fr9")
    prob_Y = jvr6.(fr9_p)
    itmp = prob_Y .> 1e-6
    proper_Q_p = prob_Y[itmp]
    proper_Q_n = length(proper_Q_p)
    proper_Q_Y = yM[:, :, itmp]
    proper_Q_n, proper_Q_p, proper_Q_Y
end
function get_proper_P(x, Y, Z_list)
    fr10 = GurobiDirectModel(); JuMP.set_silent(fr10); JuMP.set_objective_sense(fr10, JuMP.MAX_SENSE); # JuMP.set_attribute(fr10, "Presolve", 0);
    JuMP.@variable(fr10, p[eachindex(Z_list)] >= 0); JuMP.@constraint(fr10, sum(p) == 1);
    # JuMP.@constraint(fr10, my_ip(p, Z_list) .== MZ) # this style triggers Warning
    JuMP.@constraint(fr10, sum(p[i] * Z_list[i] for i in eachindex(p)) .== MZ) # this equivalent style also triggers
    JuMP.set_objective_function(fr10, ip(p, [get_f_value(x, Y, Z) for Z in Z_list]));
    opt_ass_opt(fr10, "fr10")
    prob_Z = jvr6.(p)
    itmp = prob_Z .> 1e-6
    proper_P_p = prob_Z[itmp]
    proper_P_n = length(proper_P_p)
    proper_P_Z = Z_list[itmp]
    proper_P_n, proper_P_p, proper_P_Z
end
fr9 = GurobiDirectModel(); JuMP.set_objective_sense(fr9, JuMP.MAX_SENSE); JuMP.set_silent(fr9);
JuMP.@variable(fr9, fr9_p[1:NY] >= 0); JuMP.@constraint(fr9, sum(fr9_p) == 1);
# JuMP.@constraint(fr9, my_ip(fr9_p, [yM[:, :, i] for i in 1:size(yM, 3)]) .== MY);
JuMP.@constraint(fr9, sum(fr9_p[i] * yM[:, :, i] for i in 1:NY) .== MY);
proper_Q = get_proper_Q(x) # âś…
JuMP.fix.(tr2_x, x);
i = 1; # for one Y
Y = proper_Q[3][:, :, i];
JuMP.fix.(tr2_Y, Y);
opt_ass_opt(tr2, "tr2");
b2 = jvr6.(tr2_b2); # 🥑
set_argZ_objective(x, Y, b2);
opt_ass_optortime(argZ, "argZ");
Z_new = jvr6.(argZ_Z);
recruit(Z_list, Z_new);
proper_P = get_proper_P(x, Y, Z_list) # âś…
cn, px = get_cut_oLh_coeff(x, proper_Q, proper_P); add_cut_for_oLh(cn, px)
JuMP.set_objective_function(trh, ip(brcs(CST), trh_u) + oLh) # ❄️ stable obj
opt_ass_opt(trh, "trh")
@info "trh's initialization succeed, main loop begins"

updVec = trues(2)
while true
    opt_ass_opt(trh, "trh")
    x = jvr6.(trh_x); oLh_check = JuMP.value(oLh); lb = JuMP.objective_bound(trh);
    JuMP.fix.(tr2_x, x) # 🖌️
    proper_Q = get_proper_Q(x) # âś…
    updVec .= false
    for j in 1:proper_Q[1]
        Y = proper_Q[3][:, :, j]
        JuMP.fix.(tr2_Y, Y) # 🖌️
        opt_ass_opt(tr2, "tr2")
        b2 = jvr6.(tr2_b2); oΛ2_check = JuMP.value(oΛ2) # 🥑✂️
        set_argZ_objective(x, Y, b2); opt_ass_optortime(argZ, "argZ")
        Z_new = jvr6.(argZ_Z)
        cn, px, pY, pb = JuMP.value(argZ_cn), jvr6.(argZ_c_x), jvr6.(argZ_c_Y), jvr6.(argZ_c_b2)
        if oΛ2_check < cn + ip(px, x) + ip(pY, Y) + ip(pb, b2) - UPDTOL
            add_cut_for_oΛ2(cn, px, pY, pb)
            updVec[2] = true
        end   
        recruit(Z_list, Z_new)
        proper_P = get_proper_P(x, Y, Z_list) # âś…
        add_cut_for_oLh(x, proper_Q, proper_P, oLh_check) && (updVec[1] = true)
    end
    @info "$updVec, lb = $lb"
    all(updVec .== false) && break
end

The result should be like

[ Info: Bool[1, 1], lb = -208.90723332481141
┌ Warning: The addition operator has been used on JuMP expressions a large number of times. This warning is safe to ignore but may indicate that model generation is slower than necessary. For performance reasons, you should not add expressions in a loop. Instead of x += y, use add_to_expression!(x,y) to modify x in place. If y is a single variable, you may also use add_to_expression!(x, coef, y) for x += coef*y.
â”” @ JuMP K:\judepot1112\packages\JuMP\i68GU\src\operators.jl:282
[ Info: Bool[1, 1], lb = 0.6
...
3

The warning happens when you use operators that do not use the efficient in-place mutable arithmetics.

We have dot specialized for dot(::Vector, ::Vector), but not dot(::Vector, ::Vector{<:Matrix}).

One suggestion would be

function ip(model, x, y)
    @expression(model, sum(xi * yi for (xi, yi) in zip(x, y)))
end

JuMP.@constraint(model, ip(model, p, M_vector) .== M_mean)

or

function ip(x::AbstractVector{JuMP.VariableRef}, y)
    model = owner_model(first(x))
    @expression(model, sum(xi * yi for (xi, yi) in zip(x, y)))
end

JuMP.@constraint(model, ip(p, M_vector) .== M_mean)
1 Like
4

@blegat do we want to cover this case? I don’t know if it’s very common.

5

Another option that we could consider fixing is:

julia> using JuMP

julia> model = Model()
A JuMP Model
├ solver: none
├ objective_sense: FEASIBILITY_SENSE
├ num_variables: 0
├ num_constraints: 0
â”” Names registered in the model: none

julia> @variable(model, x[1:3])
3-element Vector{VariableRef}:
 x[1]
 x[2]
 x[3]

julia> y = [rand(2, 2) for _ in 1:3]
3-element Vector{Matrix{Float64}}:
 [0.5701080035016627 0.6022404547340912; 0.5605957743951222 0.6061453626876424]
 [0.362357957737371 0.5970481709139359; 0.03292266575902236 0.34169043415484046]
 [0.8061135686585252 0.764917425247723; 0.7067819448426711 0.12129596496862716]

julia> @which x' * y
*(u::Union{LinearAlgebra.Adjoint{T, var"#s997"}, LinearAlgebra.Transpose{T, var"#s997"}} where {T, var"#s997"<:(AbstractVector)}, v::AbstractVector)
     @ LinearAlgebra ~/.julia/juliaup/julia-1.10.7+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/adjtrans.jl:464

julia> x' * y
ERROR: MethodError: no method matching zero(::Type{Matrix{AffExpr}})

Closest candidates are:
  zero(::Type{Union{}}, Any...)
   @ Base number.jl:310
  zero(::Type{Pkg.Resolve.FieldValue})
   @ Pkg ~/.julia/juliaup/julia-1.10.7+0.x64.apple.darwin14/share/julia/stdlib/v1.10/Pkg/src/Resolve/fieldvalues.jl:38
  zero(::Type{Dates.Date})
   @ Dates ~/.julia/juliaup/julia-1.10.7+0.x64.apple.darwin14/share/julia/stdlib/v1.10/Dates/src/types.jl:439
  ...

Stacktrace:
 [1] neutral_element(::typeof(+), T::Type)
   @ MutableArithmetics ~/.julia/packages/MutableArithmetics/BLlgj/src/reduce.jl:25
 [2] fused_map_reduce(::typeof(MutableArithmetics.add_mul), ::LinearAlgebra.Adjoint{VariableRef, Vector{…}}, ::Vector{Matrix{…}})
   @ MutableArithmetics ~/.julia/packages/MutableArithmetics/BLlgj/src/reduce.jl:44
 [3] _dot_nonrecursive(lhs::LinearAlgebra.Adjoint{VariableRef, Vector{VariableRef}}, rhs::Vector{Matrix{Float64}})
   @ MutableArithmetics ~/.julia/packages/MutableArithmetics/BLlgj/src/dispatch.jl:49
 [4] *(u::LinearAlgebra.Adjoint{VariableRef, Vector{VariableRef}}, v::Vector{Matrix{Float64}})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.7+0.x64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/adjtrans.jl:464
 [5] top-level scope
   @ REPL[8]:1
Some type information was truncated. Use `show(err)` to see complete types.

julia> z = rand(3)
3-element Vector{Float64}:
 0.11627605066264912
 0.3058601513730612
 0.44928823383281546

julia> z' * y
2Ă—2 Matrix{Float64}:
 0.539298  0.596308
 0.392802  0.229487
6

When I was a noob I used to write inner product with x' * y, but currently I find that ip(x, y) is more decent, where ip is imported from LinearAlgebra.dot.

I’ve updated my 2nd post, yet I cannot find out where that Warning is from.
You can tell me once you’ve found.

7

Here’s the code for the warning:

There is a model-level variable that tracks how many times two AffExpr are added together, and the warning is thrown once this reaches 20,000. So you won’t find the warning on small examples.

I’m guessing your issue is:

JuMP.set_objective_function(fx2Z,
          ip(f_c_x, x)
        + ip(f_c_Y, Y)
        + ip(f_c_Z, Z)
    )

Just do:

JuMP.@objective(fx2Z, Max,
          ip(f_c_x, x)
        + ip(f_c_Y, Y)
        + ip(f_c_Z, Z)
    )
1 Like
8

Although your finding seems to be workable. But isn’t the former functional API more fundamental, intuitive and concise?
I’m going to guess you want to fix this problem in JuMP, will you? :smiling_face_with_tear:
Is there any motivation for me to resort to the macro API?

9

I’m going to guess you want to fix this problem in JuMP, will you?

I think you might mean “you won’t fix this problem”, in which case, no, we won’t be changing the current behavior.

For the motivation, see Performance tips · JuMP

Note that if you followed the text in the warning, it would suggest:

expr = ip(f_c_x, x)
add_to_expression!(expr, ip(f_c_Y, Y))
add_to_expression!(expr, ip(f_c_Z, Z))
JuMP.set_objective_function(fx2Z, expr)

The key is that as you have written the original function, JuMP doesn’t know that it can modify the results of ip(f_c_x, x) in-place. Therefore, it needs to allocate memory for each temporary expression.

I didn’t look through your code in detail, but this probably happens in a bunch of other places as well.

10

Yes, I take your point. @expression should be used for performance.
So is it true that LinearAlgebra.dot is now well supported with performance guarantees by the JuMP macros? like
@expression, @objective, @constraint, @build_constraint

11

Perhaps this is a good idea

import JuMP
import Gurobi
GRB_ENV = Gurobi.Env()
macro set_objective_function(model, f) # use macro to retain performance
    esc(:( JuMP.@objective($model, JuMP.objective_sense($model), $f) ))
end
function GurobiDirectModel(str)
    m = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
    occursin("m", str) && JuMP.@objective(m, Min, JuMP.objective_function(m))
    occursin("M", str) && JuMP.@objective(m, Max, JuMP.objective_function(m))
    occursin("s", str) && JuMP.set_silent(m)
    m
end
model = GurobiDirectModel("ms") # OBJ_SENSE setting is one-off
for i in 1:3 # but we may want to change the objective function over time
    @set_objective_function(model, i)
end
12

For some/most cases. But there might be ones that we don’t have a good fallback for. @blegat can chime in here.

But maybe we should just remove the warning. I wonder if it does more harm than good. Are you happy with the performance of dot? If so, ignore the warning.

You can do

new_obj = @expression(model, expression)
set_objective_function(model, new_obj)

But really. Just use @objective(model, sense, expression).

If I had one comment about your code, is that it looks too “clever.” I would write your last example as:

import JuMP
import Gurobi
GRB_ENV = Gurobi.Env()
function gurobi_direct_model(;sense::JuMP.OptimizationSense, silent::Bool)
    model = JuMP.direct_model(Gurobi.Optimizer(GRB_ENV))
    JuMP.set_objective_sense(model, sense)
    if silent
        JuMP.set_silent(model)
    end
    return model
end
model = gurobi_direct_model(; sense = JuMP.MIN_SENSE, silent = true)
for i in 1:3
    obj = JuMP.@expression(model, i)
    JuMP.set_objective_function(model, obj)
end

At this point, we’re essentially not adding new functionality or macros to JuMP if there is a simple alternative (even if the alternative is a few extra lines of code). The bar is very high, and we won’t be adding a macro to set the objective function without modifying the objective sense.

1 Like
13

Definitely. In the scope of (Mi)LP, the fundamental language is vector addition and scalar multiplication. Inner Product is one of the essences in this context. I’ve been writing julia code with JuMP for a long time, so I am fairly experienced.

Benefits include: 1. a lot more concise than sum(p[i] * x[i] for i in 1:I) (Do I need to add a more line to specify I, otherwise sum(p[i] * x[i] for i in eachindex(p))) which is way too lengthy.
2. By simply writing ip(p, x), we could also ensure behind the scenes that the size(p) == size(x), which is definitely wanted.

In my opinion, building an optimization model with JuMP (or whatever someelse modeling language) entails a lot of effort to assure correctness. Therefore it is necessary to maintain conciseness.

For example, If I could write

@constraint(model, x .== y)

I would never want to write

@constraint(model, [i = 1:I, j = 1:J], x[i, j] == y[i, j])

The latter introduce unnecessary indices i, j and Int64’s I, J.
The same reasonapplies for sum(p[i] * x[i] for i in 1:I).

But it is not recommended by the document :thinking:
Does it mean that I could also write

c1 = @build_constraint(x >= y + 1)
JuMP.add_constraint(model, c1)

If this usage is correct, I would suggest that it can be written into the docstring :smile:

14

Can you point to where this is?

We’re not going to be documenting @build_constraint in more detail. It has a very particular use-case: it is mainly intended for use in lazy constraint and user cut callbacks. There is no good reason to use it in favour of @constraint.

15 
help?> JuMP.set_objective_function
  set_objective_function(model::GenericModel, func::MOI.AbstractFunction)
  set_objective_function(model::GenericModel, func::AbstractJuMPScalar)
  set_objective_function(model::GenericModel, func::Real)
  set_objective_function(model::GenericModel, func::Vector{<:AbstractJuMPScalar})

  Sets the objective function of the model to the given function.

  See set_objective_sense to set the objective sense.

  These are low-level functions; the recommended way to set the objective is with the @objective macro.

See the last line of this docstring

16

I mean, yes? The recommended way is to use @objective. But you don’t have to follow our recommendations if you know what you’re doing.

17

Admittedly, given its versatility, usage of @expression should be promoted.

18

Nothing serious though, I find a bit of inconsistency between @expression and @build_constraint.
For a theoretical standpoint, the former builds an expression, which is easier.
The latter builds a constraint, which is more complex than the former.
Yet @expression demands that the model should be involved as its 1st argument, while @build_constraint doesn’t require this, thus is more straightforward to use.

19

We’re aware of the inconsistency, but we’re not making breaking changes to JuMP any more. Most people should never know or use @build_constraint.

1 Like
20

I am clearer about what you want to express.
x' * y is more versatile than the standard LinearAlgebra.dot.

About the nice solution you’ve offered in your first post, it can deal with the Vector{Matrix} case successfully. However, it could not handle the “another example” I suggested in my first post (And it indicates that an expression that exists outside the macro could fail, were it built inside the macro). To put it more straightforward,

julia> import JuMP

julia> model = JuMP.Model(); JuMP.@variable(model, x[1:4]);

julia> y = rand(2, 3, 4) # âś… This structure boasts facilitating browsing, compared to [rand(2, 3) for _ in 1:4]
2Ă—3Ă—4 Array{Float64, 3}:
[:, :, 1] =
 0.060666  0.274775  0.00366431
 0.467394  0.12283   0.864536

[:, :, 2] =
 0.697338  0.0978811  0.0464617
 0.720297  0.737535   0.466657

[:, :, 3] =
 0.795264  0.170381  0.087852
 0.824574  0.702593  0.638999

[:, :, 4] =
 0.547624  0.398487  0.170125
 0.438258  0.716708  0.764069

julia> yslices = eachslice(y; dims = 3); # exactly the 4 matrices above

julia> # the following expression could be built on its own, but not within `@expression`

julia> success_expr = sum(l * r for (l, r) in zip(x, yslices)); # âś… doesn't have ERROR

julia> fail_expr = JuMP.@expression(                                                                                                      
           JuMP.owner_model(first(x)),                                                                                                    
           sum(l * r for (l, r) in zip(x, yslices))                                                                                       
       )
ERROR: MethodError: no method matching similar_array_type(::Type{SubArray{Float64, 2, Array{…}, Tuple{…}, true}}, ::Type{JuMP.AffExpr})
The function `similar_array_type` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  similar_array_type(::Type{SparseArrays.SparseVector{Tv, Ti}}, ::Type{T}) where {T, Tv, Ti}
   @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\implementations\SparseArrays.jl:51
  similar_array_type(::Type{Array{T, N}}, ::Type{S}) where {S, T, N}
   @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\implementations\LinearAlgebra.jl:213
  similar_array_type(::Type{SparseArrays.SparseMatrixCSC{Tv, Ti}}, ::Type{T}) where {T, Tv, Ti}
   @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\implementations\SparseArrays.jl:58
  ...

Stacktrace:
  [1] promote_operation(op::typeof(*), ::Type{JuMP.VariableRef}, A::Type{SubArray{Float64, 2, Array{…}, Tuple{…}, true}})
    @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\implementations\LinearAlgebra.jl:231
  [2] promote_operation_fallback(op::typeof(MutableArithmetics.add_mul), ::Type{Matrix{…}}, ::Type{JuMP.VariableRef}, ::Type{SubArray{…}})
    @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\interface.jl:82
  [3] promote_operation(::typeof(MutableArithmetics.add_mul), ::Type, ::Type, ::Type)
    @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\interface.jl:113
  [4] mutability(::Type, ::Function, ::Type, ::Type, ::Type)
    @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\interface.jl:273
  [5] mutability(::Matrix{…}, ::Function, ::Matrix{…}, ::JuMP.VariableRef, ::SubArray{…})
    @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\interface.jl:281
  [6] operate!!(::typeof(MutableArithmetics.add_mul), ::Matrix{…}, ::JuMP.VariableRef, ::SubArray{…})
    @ MutableArithmetics K:\judepot1112\packages\MutableArithmetics\BLlgj\src\rewrite.jl:112
  [7] macro expansion
    @ K:\judepot1112\packages\MutableArithmetics\BLlgj\src\rewrite_generic.jl:307 [inlined]
  [8] macro expansion
    @ K:\judepot1112\packages\MutableArithmetics\BLlgj\src\rewrite.jl:371 [inlined]
  [9] macro expansion
    @ K:\judepot1112\packages\JuMP\i68GU\src\macros.jl:264 [inlined]
 [10] macro expansion
    @ K:\judepot1112\packages\JuMP\i68GU\src\macros\@expression.jl:86 [inlined]
 [11] top-level scope
    @ REPL[14]:1
Some type information was truncated. Use `show(err)` to see complete types.