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.
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).
OK, if you just want to get an impression what I’m dealing with
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.
Do they need to be abs2 instead of abs? For Convex, abs will be simpler, since it reformulates 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. 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) + ....
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
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
OK, so I managed to derive and implement the Hessian H of my problem explicitely.
Now my objective is just
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
Quadratic forms supported only for semidefinite matrices
Apparently, the tolerance for the posdef check is quite strict because
> 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)…
I found the assume_psd switch.
So at the end of the day I realized, with a problem of the form
I can write the solution explicitely
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…
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.
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.
Update:
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.
A QP is an optimization problem 