# Sqrt of QuadExpr

**URL:** <https://discourse.julialang.org/t/sqrt-of-quadexpr/102797>\
**Category:** Optimization (Mathematical)\
**Created:** [August 14, 2023, 6:04pm UTC](https://discourse.julialang.org/t/sqrt-of-quadexpr/102797 "2023-08-14T18:04:02Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![i9e1](https://avatars.discourse-cdn.com/v4/letter/i/ac91a4/32.png) [@i9e1](https://discourse.julialang.org/u/i9e1)\
**Post date:** [August 14, 2023, 6:04pm UTC](https://discourse.julialang.org/t/sqrt-of-quadexpr/102797/1 "2023-08-14T18:04:02Z")

</div>

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

```julia
# 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`.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [August 14, 2023, 7:13pm UTC](https://discourse.julialang.org/t/sqrt-of-quadexpr/102797/2 "2023-08-14T19:13:04Z")

</div>

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:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![i9e1](https://avatars.discourse-cdn.com/v4/letter/i/ac91a4/32.png) [@i9e1](https://discourse.julialang.org/u/i9e1)\
**Post date:** [August 14, 2023, 7:45pm UTC](https://discourse.julialang.org/t/sqrt-of-quadexpr/102797/3 "2023-08-14T19:45:21Z")

</div>

> [@odow](#):
>
> We will never try to convert the quadratic expression back to affine.

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.

```julia
"""
    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

```