# The Optimization Problem with Nested Loops in Julia

**URL:** <https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833>\
**Category:** Performance\
**Tags:** question, performance, loops, physics, sparsearrays\
**Created:** [October 11, 2023, 8:01am UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833 "2023-10-11T08:01:27Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Umut\_Can\_Turhan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/umut_can_turhan/32/49221_2.png) [@Umut\_Can\_Turhan](https://discourse.julialang.org/u/Umut_Can_Turhan)\
**Post date:** [October 11, 2023, 8:01am UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/1 "2023-10-11T08:01:27Z")

</div>

Hi everyone, I have a problem related to the optimization problem with nested loops. I use _[QuantumOptics.jl (qojulia.org)](https://qojulia.org/)_ module to calculate for some physical situation. First of all, let me share my code and explain it.

```julia
function Hubbard_Interaction(P, Pt, cut_mb_basis, cut_off, U)
    
    # P1 and P1t are just a matrix, don't focus on it :)
    P1 = P.data
    P1t = Pt.data

    #Preety fast calculation with einsum. No problem here
    @einsum coefficient[k,l,m,n] := P1[k,i] * P1[l,i] * P1t[i,m] * P1t[i,n]

    # Sparse operator for fast modifying
    Vint_mb_cut = SparseOperator(cut_mb_basis)
    
    # The problem starts here :/    
    for k in 1:cut_off
        for l in 1:cut_off
            for m in 1:cut_off
                for n in 1:cut_off

                    # These four operators are just matrices but note that they depend on the loop indices!
                    a1t = create(cut_mb_basis, k)
                    a2t = create(cut_mb_basis, l)
                    a2 = destroy(cut_mb_basis, m)      
                    a1 = destroy(cut_mb_basis, n)

                    #Matrix multiplication pretty fast, no problem here      
                    Vint_mb_cut += U/2*coefficient[k,l,m,n]*a1t*a2t*a2*a1
                end
            end
        end
    end
    
    return Vint_mb_cut
end

```

There is a function that calculates the interaction between two particles. It’s unnecessary for non-physicist people. The main point is that there is a nested loop in my code. When cutt\_off is a big number, the code works quite slowly ☹ The big question is how can I optimize this code for big cutt\_off (or just a number from starting 1) values? For example, cut\_off = 50 😕

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [October 11, 2023, 8:20am UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/2 "2023-10-11T08:20:10Z")

</div>

It seems you are rebuilding the creation and annihilation matrices unnecessarily in the innermost loop. E.g. `a1t = create(cut_mb_basis, k)` only depends on the outermost variable and could be moved to

```julia
for k in 1:cut_off
  a1t = create(cut_mb_basis, k)
  [...]

```

and so on, which should save about 75%.

But actually it might be smarter to create all operators from the get-go. I suppose they are sparse matrices? Shouldn’t take up too much space. I’ve never worked with QuantumOptics, so that’s a guess.

---

<div class="post-metadata">

**Author:** ![Umut\_Can\_Turhan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/umut_can_turhan/32/49221_2.png) [@Umut\_Can\_Turhan](https://discourse.julialang.org/u/Umut_Can_Turhan)\
**Post date:** [October 11, 2023, 10:22am UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/3 "2023-10-11T10:22:01Z")

</div>

All of these matrices are sparse arrays, they help to improve code performance of course. I tried your opinion, but the results are not convincing. Because when cut\_off = 75 the code still works slowly.

---

<div class="post-metadata">

**Author:** ![kellertuer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kellertuer/32/220707_2.png) [@kellertuer](https://discourse.julialang.org/u/kellertuer)\
**Post date:** [October 11, 2023, 11:22am UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/4 "2023-10-11T11:22:49Z")

</div>

Can you post your updated code maybe? Because the one above allocates 4 matrices in every inner loop run, which looks really costly (though I am not sure what create and destroy do).

---

<div class="post-metadata">

**Author:** ![skleinbo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/skleinbo/32/36080_2.png) [@skleinbo](https://discourse.julialang.org/u/skleinbo)\
**Post date:** [October 11, 2023, 11:33am UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/5 "2023-10-11T11:33:39Z")

</div>

> [@Umut\_Can\_Turhan](#):
>
> I tried your opinion, but the results are not convincing

Did you try pre-allocating all operators outside of the loop, e.g

```julia
A = [create(cut_mb_basis, k) for k in 1:cut_off]
At = [destroy(cut_mb_basis, k) for k in 1:cut_off]

```

and then indexing into those arrays in the loop?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [October 11, 2023, 12:09pm UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/6 "2023-10-11T12:09:14Z")

</div>

Your innermost loop executes `cut_off^4` times, so anything you can do to speed it up will be critical. @skleinbo’s suggestion of pre-computing your operators looks like a good first step.

> [@Umut\_Can\_Turhan](#):
>
> `#Matrix multiplication pretty fast, no problem here`

This is still going to allocate in your inner loop. If the matrices are small, you’re going to pay a big overhead price for these allocations, which you can avoid if you pre-allocate scratch space and use in-place operations like `mul!`.

How big are your matrices? If they are a small fixed size (say 10x10 or smaller), you might use StaticArrays.jl. (This is also in the performance tips.)

---

<div class="post-metadata">

**Author:** ![simeonschaub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/simeonschaub/32/216566_2.png) [@simeonschaub](https://discourse.julialang.org/u/simeonschaub)\
**Post date:** [October 11, 2023, 3:53pm UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/7 "2023-10-11T15:53:25Z")

</div>

The first thing to think about would be whether you could take advantage of any symmetries here. Assuming your operators are all bosonic, both `coefficient` and your operator product should be symmetric under the exchanges k \leftrightarrow l and m \leftrightarrow n. Taking advantage of that (for example using `SymmetricTensor` from [GitHub - Ferrite-FEM/Tensors.jl: Efficient computations with symmetric and non-symmetric tensors with support for automatic differentiation.](https://github.com/Ferrite-FEM/Tensors.jl/)), you’d only need to calculate (n(n+1))^2/4 instead of n^4 matrix products, a 4x speedup!

You probably also have commutation relations for a\_i and a^\dagger\_i, so you might even be able get that down to {n+3 \choose 4} + n(n+1)/4 matrix products, which would be a 24x speedup

---

<div class="post-metadata">

**Author:** ![Umut\_Can\_Turhan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/umut_can_turhan/32/49221_2.png) [@Umut\_Can\_Turhan](https://discourse.julialang.org/u/Umut_Can_Turhan)\
**Post date:** [October 18, 2023, 8:12am UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/8 "2023-10-18T08:12:11Z")

</div>

Thank you! It’s a very effective way to construct any array quickly.

---

<div class="post-metadata">

**Author:** ![Umut\_Can\_Turhan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/umut_can_turhan/32/49221_2.png) [@Umut\_Can\_Turhan](https://discourse.julialang.org/u/Umut_Can_Turhan)\
**Post date:** [October 23, 2023, 8:02am UTC](https://discourse.julialang.org/t/the-optimization-problem-with-nested-loops-in-julia/104833/9 "2023-10-23T08:02:55Z")

</div>

I generally use the very big matrices/tensors that have 9000x9000 shapes.