# Performance of optimization problem in a loop

**URL:** <https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374>\
**Category:** Optimization (Mathematical)\
**Created:** [September 10, 2020, 12:13pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374 "2020-09-10T12:13:18Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![Oliver\_Lylloff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oliver_lylloff/32/9022_2.png) [@Oliver\_Lylloff](https://discourse.julialang.org/u/Oliver_Lylloff)\
**Post date:** [September 10, 2020, 12:13pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/1 "2020-09-10T12:13:18Z")

</div>

Hello. I’m trying to improve the performance of an optimization problem. I have a Hermitian matrix and want to minimize the diagonal terms while keeping the matrix positive semidefinite.

Given a matrix S,

\text{minimize}\sum d\\ \text{subject to S}+diag(d) \geq 0 

where d is a vector variable. This is solved using Convex.jl and COSMO.jl:

```julia
using Convex, COSMO
n = 4
x = rand(ComplexF64,n); 
noise = diagm(rand(n))
S = x*x' + noise

d = Convex.Variable(n)
p = minimize(sum(d), S+Diagonal(d) in :SDP)
solve!(p, () -> COSMO.Optimizer())

Sx = S.+diagm(d.value[:]) # Solution

diag(Sx) + diag(noise) ≈ diag(S) # approx true

```

First question, is this an optimal formulation or should it be reformulated?

Second question, my real problem is larger n=80 and needs to be solved 1000s of times since my matrix S is in fact 3D. Setting up a loop:

```julia
N,N,M = size(S) # (80,80,1000)
d = Convex.Variable(N)
A = similar(S)
for i = 1:M
    p = minimize(sum(d), S[:,:,i]+Diagonal(d) in :SDP)
    solve!(p, () -> COSMO.Optimizer())
    A[:,:,i] = Hermitian(S[:,:,i].+diagm(d.value[:]))
end

```

How can I avoid to setup the problem at each iteration? Warm-starting seems like to perfect solution but I’m not sure how to do that when my constraints change at each iteration.

Thanks!

---

<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:** [September 10, 2020, 12:24pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/2 "2020-09-10T12:24:28Z")

</div>

> How can I avoid to setup the problem at each iteration? Warm-starting seems like to perfect solution but I’m not sure how to do that when my constraints change at each iteration.

Unfortunately, Convex.jl does not have a good answer to this right now.

One thing you can do is use `fix!`'d variables:

```julia
N,N,M = size(S) # (80,80,1000)
d = Convex.Variable(N)
A = similar(S)
B = Variable(size(S)[1:2])
p = minimize(sum(d), B+Diagonal(d) in :SDP)

for i = 1:M
    fix!(B, S[:, :, i])
    solve!(p, () -> COSMO.Optimizer())
    A[:,:,i] = Hermitian(S[:,:,i].+diagm(d.value[:]))
end

```

(untested, but should work). However, I’m not sure it will actually save you much time, because on every `solve!` call, Convex recalculates its [extended formulations](https://jump.dev/Convex.jl/stable/#Extended-formulations-and-the-DCP-ruleset-1). (Let me know if it does help though!). Ideally, these could be reused, but it’s not implemented ([Cache conic forms between `solve!`s? · Issue #318 · jump-dev/Convex.jl · GitHub](https://github.com/jump-dev/Convex.jl/issues/318)) (and also might need a lot of thought to implement correctly).

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [September 10, 2020, 1:07pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/3 "2020-09-10T13:07:20Z")

</div>

Could you instead parameterize the Cholesky factor directly? That way your constraint is just a box constraint on the diagonal of the factor (for uniqueness). Getting the diagonal of a matrix from its Cholesky factor requires a little bit of computation, but I wouldn’t guess that it would be a bottleneck.

EDIT: read the problem wrong, sorry. I don’t think this is a helpful suggestion.

---

<div class="post-metadata">

**Author:** ![Oliver\_Lylloff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oliver_lylloff/32/9022_2.png) [@Oliver\_Lylloff](https://discourse.julialang.org/u/Oliver_Lylloff)\
**Post date:** [September 10, 2020, 1:52pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/4 "2020-09-10T13:52:09Z")

</div>

Thanks for the clarification! Looking into the timings - i does seem to improve but need to check if the two actually finds the same result. Will update…

---

<div class="post-metadata">

**Author:** ![Oliver\_Lylloff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oliver_lylloff/32/9022_2.png) [@Oliver\_Lylloff](https://discourse.julialang.org/u/Oliver_Lylloff)\
**Post date:** [September 10, 2020, 1:53pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/5 "2020-09-10T13:53:53Z")

</div>

Interesting suggestion. Could you elaborate a bit, I’m not quite sure how to do what you propose. Thanks!

---

<div class="post-metadata">

**Author:** ![Oliver\_Lylloff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oliver_lylloff/32/9022_2.png) [@Oliver\_Lylloff](https://discourse.julialang.org/u/Oliver_Lylloff)\
**Post date:** [September 10, 2020, 2:28pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/6 "2020-09-10T14:28:43Z")

</div>

Looks like a 2X speed-up. Do get a few `Problem status INFEASIBLE` with the `fix!` not sure why but overall it works. Thanks!

---

<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:** [September 10, 2020, 2:44pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/7 "2020-09-10T14:44:53Z")

</div>

Nice! Glad to hear there was a benefit. Convex’s code _should_ treat `fix!`'d variable exactly the same as if you put the matrix there instead, so I’d suspect the INFEASIBLE status is unrelated, but there could be a bug or something else strange going on. If you can get a reproducible example where there is a difference (ideally on a small/simple problem), please file an issue so I can look into it in more detail ([https://github.com/jump-dev/Convex.jl/issues/new](https://github.com/jump-dev/Convex.jl/issues/new)).

---

<div class="post-metadata">

**Author:** ![cgeoga](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cgeoga/32/216186_2.png) [@cgeoga](https://discourse.julialang.org/u/cgeoga)\
**Post date:** [September 11, 2020, 2:14am UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/8 "2020-09-11T02:14:55Z")

</div>

Apologies! I read the problem too quickly and thought you were trying to find a matrix that minimized something. I’ve edited my above comment to say the same.

---

<div class="post-metadata">

**Author:** ![Oliver\_Lylloff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oliver_lylloff/32/9022_2.png) [@Oliver\_Lylloff](https://discourse.julialang.org/u/Oliver_Lylloff)\
**Post date:** [September 11, 2020, 6:48am UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/9 "2020-09-11T06:48:03Z")

</div>

Inspired by this post [Near positive definite (NearPD) of a symmetric matrix in Julia - #2 by lostella](https://discourse.julialang.org/t/near-positive-definite-nearpd-of-a-symmetric-matrix-in-julia/45793/2) I was wondering if ProximalAlgorithms.jl could enable warm-starting for this particular problem. @lostella if you don’t mind me asking, could you suggest strategies to solve this?

---

<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:** [September 11, 2020, 9:41am UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/10 "2020-09-11T09:41:50Z")

</div>

> [@Oliver\_Lylloff](#):
>
> enable warm-starting

Not to sidetrack this, but you can try warmstarting in Convex with `solve!(problem, solver; warmstart=true)` by the way.

---

<div class="post-metadata">

**Author:** ![Oliver\_Lylloff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oliver_lylloff/32/9022_2.png) [@Oliver\_Lylloff](https://discourse.julialang.org/u/Oliver_Lylloff)\
**Post date:** [September 11, 2020, 11:02am UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/11 "2020-09-11T11:02:46Z")

</div>

Not sidetracking at all. Should have mentioned that I did try that without any benefit (to my particular problem)

---

<div class="post-metadata">

**Author:** ![migarstka](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/migarstka/32/3828_2.png) [@migarstka](https://discourse.julialang.org/u/migarstka)\
**Post date:** [September 11, 2020, 11:13am UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/12 "2020-09-11T11:13:49Z")

</div>

Hi, I am one of the developers of COSMO.jl. It seems like you could make the model creation very fast because you are just changing one side of the PSD-constraint:  
Internally in COSMO it will look like this:

- diag(d) + X = S, X \succeq 0

and just `S` changes. This means you can recycle a lot of the steps of the algorithm, i.e. just setup the problem once and change that part of the constraint.  
However, the dominating cost of the computation is to repeatedly compute the eigendecomposition of X at each iteration. You won’t get around this so warm-starting your initial guess of d could help to decrease the number of times you have to do that.

Furthermore, your S is dense, right? If S would be sparse, COSMO could get some further speed-ups by decomposing the PSD constraint.

One issue is that at the moment each model change in Convex.jl or JuMP triggers a completely new model copy that is passed to COSMO. I haven’t gotten around to support incremental model changes yet.

If you are sure your problems are always feasible, you can lower the infeasibility check tolerance in the solver settings with `eps_prim_inf = 1e-8` and `eps_dual_inf = 1e-8`.

---

<div class="post-metadata">

**Author:** ![Oliver\_Lylloff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oliver_lylloff/32/9022_2.png) [@Oliver\_Lylloff](https://discourse.julialang.org/u/Oliver_Lylloff)\
**Post date:** [September 11, 2020, 11:47am UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/13 "2020-09-11T11:47:45Z")

</div>

Thank you for the suggestions. Yes, my matrix is dense and changing tolerances is indeed an easy way to speed of things. Will need to experiment more with this.

Seeing how COSMO treats the problem, it does look obvious that a lot can be reused. Is your suggestion not similar to the `fix!` in the Convex code above? Looking at the COSMO examples, e.g., [https://github.com/oxfordcontrol/COSMO.jl/blob/master/examples/closest\_correlation\_matrix.jl](https://github.com/oxfordcontrol/COSMO.jl/blob/master/examples/closest_correlation_matrix.jl)  
I see that the example is using JuMP, do you think that switching to that with a `@constraint` in the loop could make a significant change? Will need to test it out…

---

<div class="post-metadata">

**Author:** ![lostella](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lostella/32/356_2.png) [@lostella](https://discourse.julialang.org/u/lostella)\
**Post date:** [September 11, 2020, 7:38pm UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/14 "2020-09-11T19:38:33Z")

</div>

> [@Oliver\_Lylloff](#):
>
> I was wondering if ProximalAlgorithms.jl could enable warm-starting for this particular problem

I think all implemented methods in ProximalAlgorithms require the initial iterate to be provided, see for example [here](https://github.com/kul-forbes/ProximalAlgorithms.jl/blob/d96bc871552bb412bc239763fa71f6961ef4b23b/src/algorithms/forwardbackward.jl#L167). Whether the (approximate) solution to each of your problems is a good starting point for the next one is hard to say, of course.

If I’m not wrong the dual problem to yours should be solvable via proximal gradient or Douglas-Rachford splitting? With any of these the costly operations will definitely be the SVD computation, you might want to give some of the quasi-Newton variants of these a try (see [this](https://github.com/kul-forbes/ProximalAlgorithms.jl/blob/master/src/algorithms/panoc.jl) and [this](https://github.com/kul-forbes/ProximalAlgorithms.jl/blob/master/src/algorithms/drls.jl)).

I’m sorry the package doesn’t yet have documentation, I would really need to duplicate myself to take care of that!

_Edit:_ fixed link to quasi-Newton Douglas-Rachford.

---

<div class="post-metadata">

**Author:** ![migarstka](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/migarstka/32/3828_2.png) [@migarstka](https://discourse.julialang.org/u/migarstka)\
**Post date:** [September 15, 2020, 7:30am UTC](https://discourse.julialang.org/t/performance-of-optimization-problem-in-a-loop/46374/15 "2020-09-15T07:30:37Z")

</div>

I think it’s worth checking where you spent most of the time. If the eigenvalue decomposition of your matrix iterate takes \>80% of the time, you probably won’t be able to achieve significant improvements by assembling the problem more efficiently.

I don’t think switching to JuMP would make a difference as both JuMP and Convex pass a completely new model to COSMO once you change the model. I am also not sure if JuMP can handle complex variables.