Better integration with JuMP

The examples in this documentation generally used the simple and straightforward method of converting the trees and values to JuMP system, which depends on algebraic operators working transparently with JuMP values within function substitute.

Substitution folding problem

Despite the simplicity, this approach is sometimes sub-optimal, especially in cases when the result of the substitution is recalculated with added values. For example, in the naive case, JuMP is forced to successively build representations for all intermediate expressions with incomplete variables, until all variables are in place. In turn, this may very easily reach a quadratic computational complexity.

More generally, any representation of substitution result that "does not reduce() easily" will suffer from this problem. A different (often specialized) approach is thus needed.

Solution: Prevent successive folding

For such cases, it is recommended to replace the substitute calls with a custom function that can interpret the required Values itself, and converts them without the overhead of creating temporary values.

import ConstraintTrees as C

First, let's create a lot of variables, and a constraint that will usually trigger this problem (and a JuMP warning) if used with normal substitute:

x = :vars^C.variables(keys = Symbol.("x$i" for i = 1:1000), bounds = C.Between(0, 10))
x *= :sum^C.Constraint(sum(C.value.(values(x.vars))))

ConstraintTrees.Tree{ConstraintTrees.Constraint} with 2 elements:
  :sum  => ConstraintTrees.Constraint(ConstraintTrees.LinearValue(#= ... =#))
  :vars => ConstraintTrees.Tree{ConstraintTrees.Constraint}(#= 1000 elements =#)

Now, imagine the expressions are represented e.g. by sparse vectors of fixed size (as common in linear-algebraic systems). We can produce the vectors efficiently as follows:

import SparseArrays: sparsevec
v = x.sum.value

value_in_a_vector = sparsevec(v.idxs, v.weights, 1000)

1000-element SparseArrays.SparseVector{Float64, Int64} with 1000 stored entries:
  [1   ]  =  1.0
  [2   ]  =  1.0
  [3   ]  =  1.0
  [4   ]  =  1.0
  [5   ]  =  1.0
  [6   ]  =  1.0
  [7   ]  =  1.0
          ⋮
  [993 ]  =  1.0
  [994 ]  =  1.0
  [995 ]  =  1.0
  [996 ]  =  1.0
  [997 ]  =  1.0
  [998 ]  =  1.0
  [999 ]  =  1.0
  [1000]  =  1.0

This usually requires only a single memory allocation, and runs in time linear with the number of variables in the value. As an obvious downside, you need to implement this functionality for all kinds of Values you encounter.

Solution for JuMP

LinearValues can be translated to JuMP's AffExprs:

using JuMP, GLPK

function substitute_jump(val::C.LinearValue, vars)
    e = AffExpr() # unfortunately @expression(model, 0) is not type stable and gives an Int
    for (i, w) in zip(val.idxs, val.weights)
        if i == 0
            add_to_expression!(e, w)
        else
            add_to_expression!(e, w, vars[i])
        end
    end
    e
end

model = Model(GLPK.Optimizer)
@variable(model, V[1:1000])
jump_value = substitute_jump(x.sum.value, V)

\[ V_{1} + V_{2} + V_{3} + V_{4} + V_{5} + V_{6} + V_{7} + V_{8} + V_{9} + V_{10} + V_{11} + V_{12} + V_{13} + V_{14} + V_{15} + V_{16} + V_{17} + V_{18} + V_{19} + V_{20} + V_{21} + V_{22} + V_{23} + V_{24} + V_{25} + V_{26} + V_{27} + V_{28} + V_{29} + V_{30} + [[\ldots\text{940 terms omitted}\ldots]] + V_{971} + V_{972} + V_{973} + V_{974} + V_{975} + V_{976} + V_{977} + V_{978} + V_{979} + V_{980} + V_{981} + V_{982} + V_{983} + V_{984} + V_{985} + V_{986} + V_{987} + V_{988} + V_{989} + V_{990} + V_{991} + V_{992} + V_{993} + V_{994} + V_{995} + V_{996} + V_{997} + V_{998} + V_{999} + V_{1000} \]

This function can be re-used in functions like optimized_vars as shown in other examples in the documentation.

For QuadraticValues, the same approach extends only with a minor modification:

function substitute_jump(val::C.QuadraticValue, vars)
    e = QuadExpr() # unfortunately @expression(model, 0) is not type stable and gives an Int
    for ((i, j), w) in zip(val.idxs, val.weights)
        if i == 0 && j == 0
            add_to_expression!(e, w)
        elseif i == 0 # the symmetric case is prohibited
            add_to_expression!(e, w, vars[j])
        else
            add_to_expression!(e, w, vars[i], vars[j])
        end
    end
    e
end

qvalue = 123 + (x.vars.x1.value + x.vars.x2.value) * (x.vars.x3.value - 321)
jump_qvalue = substitute_jump(qvalue, V)

\[ V_{1}\times V_{3} + V_{2}\times V_{3} - 321 V_{1} - 321 V_{2} + 123 \]

This page was generated using Literate.jl.