# Help speeding up problem with Convex.jl

**URL:** <https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005>\
**Category:** Optimization (Mathematical)\
**Tags:** question, convex-optimization\
**Created:** [April 4, 2022, 1:31pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005 "2022-04-04T13:31:27Z")\
**Posts on this page:** 12\
**Page:** 2

<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:** [April 6, 2022, 9:59am UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/21 "2022-04-06T09:59:46Z")

</div>

> then it makes more sense to me to do it with _COSMO_

If you only want to solve this one problem, and you only want to use COSMO.jl, then sure, it makes sense to code in COSMO directly.

> who knows if it will be faster, or slower, or if it will work at all

Again, if we had a reproducible example, someone might be able to point out how to make Convex.jl faster, or suggest alternatives. There isn’t much to do without one.

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [April 6, 2022, 11:08am UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/22 "2022-04-06T11:08:13Z")

</div>

If the delay is really only from some kind of joint effect when all the terms from 100 lines of code are present, it seems unlikely this will ever really get solved unless you or someone else digs into the internals with the full code. However, I would guess that either portion of it counts for a lot of the slowdown, or else a bunch of pieces each contribute roughly the same order of magnitude slowdown and the 10s slowdown is the sum of all of them. In either case, there would be something small you could share, e.g. “each of these pieces takes 20 ms and I’ve got a zillion of them so overall it adds 10s”, or “this piece takes 8s when some parameter n = 10000”, or whatever.

I get it’s frustrating that it’s not just already fast, or there isn’t some quick tweak to make, but it seems like that’s the scenario here, and if you don’t want to share something so that other folks can try to debug and improve the software, it’ll probably never be fast for this (unless someone else comes along with a similar problem and figures out what bit is slow).

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [April 6, 2022, 11:24am UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/23 "2022-04-06T11:24:11Z")

</div>

OK, if you just want to get an impression what I’m dealing with

```julia

function kinetic_energy_triangles(xs, vs)
    return abs2(
        energy_gradient_triangles(xs[1], xs[2], xs[3]) ⋅ vs[1] +
        energy_gradient_triangles(xs[2], xs[3], xs[1]) ⋅ vs[2] +
        energy_gradient_triangles(xs[3], xs[1], xs[2]) ⋅ vs[3]
    )
end

_vector(edge, vertices) = vertices[edge[2]] - vertices[edge[1]]
_corner_vectors(triangle, vertices) = map(i -> vertices[i], triangle)

# first simplified implementation with constant target vertex distance
function kinetic_energy(
        edges::AbstractVector,
        triangles::AbstractVector{T},
        vertices,
        velocities,
        square_ideal_edge_lengths::AbstractVector,
        compactness_weight = 0.4
) where {T <: Triangle}
    t1 = map(zip(edges, square_ideal_edge_lengths) |> collect) do (edge, square_ideal)
        return kinetic_energy_edges(
            _vector(edge, vertices),
            _vector(edge, velocities),
            square_ideal
        )
    end
    t2 = map(triangles) do triangle
        return kinetic_energy_triangles(
            _corner_vectors(triangle, vertices),
            _corner_vectors(triangle, velocities)
        )
    end
    return sum(t1) + compactness_weight * sum(t2)
end

```

As you can see I’m map-reducing over thousands of `abs2()` terms. So I think the call-graph just gets quite large.

The quantity I’m optimizing for is the array of 3-vectors called `velocities`.

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [April 6, 2022, 12:57pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/24 "2022-04-06T12:57:19Z")

</div>

Do they need to be `abs2` instead of `abs`? For Convex, `abs` will be simpler, since it [reformulates](https://github.com/jump-dev/Convex.jl/blob/28bc9f8628b0227a7e19350339647287a1a8937f/src/atoms/lp_cone/abs.jl#L39-L58) `abs` (for real x) as |x| = min\_{t \geq 0, -x \leq t, x \leq t} t. Then for `abs2` it needs to compose that with the extended formulation for squaring, which requires [second-order cone constraints](https://github.com/jump-dev/Convex.jl/blob/e3ca4249495b836e437043f37f34d86163297ef2/src/atoms/second_order_cone/qol_elementwise.jl#L42). This adds to both problem formulation time and solve time.

Also, Convex interpets `abs` of a vector as component-wise, so you can reduce formulation time if you can do things like `sum(abs(v))` instead of `abs(x1) + abs(x2) + ...`.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [April 6, 2022, 1:18pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/25 "2022-04-06T13:18:17Z")

</div>

Since this is an established method I cannot simply change the objective function (I’m pretty sure with the constraints `abs` would give a different result than `abs2`)

I dont think I can take advantage of the vector method either—I already tried this with `sumsquares` and got

```julia
MethodError: no method matching sumsquares(::Vector{Convex.QolElemAtom})
Closest candidates are:
  sumsquares(::Convex.AbstractExpr) at /home/arbeit/.julia/packages/Convex/FQF1R/src/atoms/second_order_cone/qol_elementwise.jl:65

```

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [April 7, 2022, 12:42pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/26 "2022-04-07T12:42:02Z")

</div>

OK, so I managed to derive and implement the Hessian `H` of my problem explicitely.

Now my objective is just

```julia
quadform(x, H)

```

(I could plug this matrix directly into COSMO but I wanted to give Convex another chance and provide everyone who helped with a bit of closure.)

Unfortunately, I get

```julia
Quadratic forms supported only for semidefinite matrices

```

Apparently, the tolerance for the posdef check is quite strict because

```julia
> minimum(eigvals(H))
-5.60169996052515e-13

```

For a 2400x2400 matrix I think this is a reasonable numerical error—so maybe the tolerance should be relaxed a bit for `quadform` (or scaled with the size of the problem)…

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [April 7, 2022, 1:35pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/27 "2022-04-07T13:35:07Z")

</div>

I found the `assume_psd` switch.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [April 7, 2022, 6:50pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/28 "2022-04-07T18:50:46Z")

</div>

So at the end of the day I realized, with a problem of the form

\text{minimize}\quad xPx\\ \text{with constraint}\quad Cx=b

I can write the solution explicitely

x=P^{-1}C^T(CP^{-1}C^T)^{-1}b

and solve with BLAS (and so I need neither Convex, nor COSMO, nor JuMP).

The operation above takes about 1.5 seconds.

Needless to say: Please let me know if I’m shooting myself in the foot here…

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [April 7, 2022, 7:08pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/29 "2022-04-07T19:08:03Z")

</div>

Looks different to your original formulation, but fine if everything is invertible (especially C?). But I believe then it isn’t a minimization problem anymore. Nope, that looks like a QP.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [April 7, 2022, 7:41pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/30 "2022-04-07T19:41:17Z")

</div>

Yeah, I essentially had to re-write the whole procedure. I knew the problem had this structure but I hadn’t explicitely derived the components of P when I first raised my question here. It would have been counter-productive to phrase the problem this way without knowing the value of P.

P is positive semi-definite by construction. C is not a square matrix and at no point needs to be inverted. CP^{-1}C^T is square and invertible.

---

<div class="post-metadata">

**Author:** ![josuagrw](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/josuagrw/32/1015_2.png) [@josuagrw](https://discourse.julialang.org/u/josuagrw)\
**Post date:** [April 25, 2022, 2:43pm UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/31 "2022-04-25T14:43:21Z")

</div>

Update:

> [@josuagrw](#):
>
> I can write the solution explicitely
> 
> x=P^{-1}C^T(CP^{-1}C^T)^{-1}b
> 
> and solve with BLAS (and so I need neither Convex, nor COSMO, nor JuMP).

So in the end this is a nice expression but BLAS does not provide a good result. I plugged P, C and b into _COSMO.jl_ directly and now I get a suberb solution in ca 0.2 seconds (in total maybe 1 second when I include calculating those quantities) down from 15-20 seconds before.

I am quite happy with that result, thank you to everyone who responded and shared their knowledge.

---

<div class="post-metadata">

**Author:** ![cvanaret](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cvanaret/32/11594_2.png) [@cvanaret](https://discourse.julialang.org/u/cvanaret)\
**Post date:** [April 26, 2022, 8:11am UTC](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005/32 "2022-04-26T08:11:52Z")

</div>

> [@goerch](#):
>
> But I believe then it isn’t a minimization problem anymore. Nope, that looks like a QP.

A QP is an optimization problem 😉

[Previous page](https://discourse.julialang.org/t/help-speeding-up-problem-with-convex-jl/79005.md?page=1)