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)

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.

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)

if !isempty(x.affine_terms)
error("Cannot compute sqrt of QuadExpr with affine terms")
end
return sqrt(x.constant)
end
if term.variable_1 != term.variable_2
end
return sqrt(term.constant) * term.variable_1^2
end


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 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.

"""

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.
"""

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
linearZero   = all(iszero, values(Q.aff.terms))
constantZero = iszero(constant(Q))
# --------------------------------------------------------------------------

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
# 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