Sqrt of QuadExpr

Specific Domains Optimization (Mathematical)
1

Hi there,

Do you think sqrt of special cases of QuadExpr is useful or dangerous?
Currrently the user gets a warning when using the unsupported nonlinear operator sqrt and has to do manual reformulation (where possible).
Pro: possible reformulation / simplification can be done automatically
Con: potential modeling issues are hidden from the user, confusion when building models with parameters not yielding special case QuadExpr (sometimes works vs. sometimes errors)

The generic quadratic expression is

Q = \sum_{i,j} a_{i,j} x_i x_j + \sum_i b_i x_i + c

for which we cannot take the squareroot symbolically.
However, for constant and single quadratic variable, it’s possible:

Q_0 = 0 + 0 + c \\ \sqrt{Q_0} = \sqrt{c} \\ \\ Q_{i,i} = a_{i,i} x_i x_i + 0 + 0 \\ \sqrt{Q_{i,i}} = \sqrt{a_{i,i}} \cdot x_i

The solution of Q_{i,i} = 0 has a double solution in \mathbb{R}, so I consider this reformulation to be save (i.e. doesnt change feasible set).

Context:

I had to play around with this a little bit when using Second Order Cones and Vector norms / sums of squares in the context of chance constraint optimization.

Using PolyChaos.jl, I build sometihng like

# mop is struct needed for numerical evaluation
# Y is PCE coefficient vector with entries AffExpr
t = mean(Y, mop)
x = std(Y, mop) # = sqrt.(var(Y, mop))
@constraint(m, [t ; x] in SecondOrderCone())

var generates the QuadExpr. However, they are always in a single variable without affine parts.

I implemented a reformulation as an alternative while debugging a model. My problem is solved with first calculating \sqrt{a_{i,i}} and setting x = sq_a .* Y.

2

In the upcoming nonlinear interface, sqrt(::QuadExpr) will return a ::NonlinearExpr. We will never try to convert the quadratic expression back to affine.

If you want this, you could always write your own function my_sqrt(::QuadExpr) = ... and then use that?

I didn’t test, but perhaps something like:

my_sqrt(x) = sqrt(x)

function my_sqrt(x::QuadExpr)
    if !isempty(x.affine_terms)
        error("Cannot compute sqrt of QuadExpr with affine terms")
    elseif length(x.quadratic_terms) > 1
        error("Cannot compute sqrt of QuadExpr with multiple quadratic terms")
    end
    if length(x.quadratic_terms) == 0
        return sqrt(x.constant)
    end
    term = x.quadratic_terms[1]
    if term.variable_1 != term.variable_2
        error("Cannot compute sqrt of QuadExpr with bilinear quadratic terms")
    end
    return sqrt(term.constant) * term.variable_1^2
end
3

Ok thanks, that was my question.

For my problem a reformulation does the trick, which avoids creating intermediate QuadExpr.

As mentioned, I implemented an own version and tested it.
Your implementation looks a bit more compact than mine, but I think some lines will not work.
(term = x.quadratic_terms[1] is unreliable, term.constant should be constant(term) and
return sqrt(term.constant) * term.variable_1 not squared)

If anyone wants to use an own function for this reduction, this code is tested.

"""
    special_sqrt(Q::GenericQuadExpr) -> GenericAffExpr

Take the sqrt of a `Q` iff it represents a special case.
Special cases can be solved easily and the sqrt is at most affine.
`GenericAffExpr` is returned.

Special cases are:

- constant: \$ c \$
- single monic quadratic expression: \$ q_{i,i} x_i^2 \$

An error is thrown if `Q` is none of these special cases, or equivalently
a sum of more than one term.
"""
function special_sqrt(Q::GenericQuadExpr)

    tooGeneric =
        "sqrt of $(typeof(Q)) can only be taken for special cases with a single quadratic term in a single variable! " *
        "Given expression has more than one term / variable or is affine instead of quadratic:"
    tooGenericErrorString = "$tooGeneric $Q\n$(JuMP._OP_HINT)"

    # check each degree individually
    quadZero     = all(iszero, values(Q.terms))
    linearZero   = all(iszero, values(Q.aff.terms))
    constantZero = iszero(constant(Q))
    # --------------------------------------------------------------------------

    onlyConstant = quadZero && linearZero
    if onlyConstant
        # ! no further evaluation needed
        # FIXME qualified Base.sqrt, strip with correct export
        c = Base.sqrt(constant(Q))
        # ! fixing concrete type for variables: VariableRef
        return GenericAffExpr{typeof(c), VariableRef}(c)
    end

    ## Check if we just have q_{i,i} x_i x_i (+ 0) for one arbitrary i
    N_quads = count(!iszero, values(Q.terms))
    singleQuadTerm = N_quads <= 1 && linearZero && constantZero
    if !singleQuadTerm
        # sqrt(V + W + ...) cannot be done analytically
        error(tooGenericErrorString)
    end

    # we have a single quadratic term and need to check if it is in the same variable
    xy = only(keys(Q.terms))
    sameVariable = xy.a == xy.b
    if !sameVariable
        # sqrt(q x*y) cannot be done analytically
        error(tooGenericErrorString)
    end
    # --------------------------------------------------------------------------

    ## Expr is q x^2, we sqrt the operands to √q x
    # same variable squared; a=x_i == b=x_i
    x = xy.a
    # - D == 0: one  real solution, "double solution"
    k = Base.sqrt(a)
    return GenericAffExpr(zero(k), [x => k])
end
1 Like