# Memory best practices in Julia with Arrays vs Vectors

**URL:** <https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307>\
**Category:** New to Julia\
**Tags:** memory-allocation, arrays, staticarrays, simulations\
**Created:** [August 17, 2024, 1:59pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307 "2024-08-17T13:59:03Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![ckawa](https://avatars.discourse-cdn.com/v4/letter/c/8c91f0/32.png) [@ckawa](https://discourse.julialang.org/u/ckawa)\
**Post date:** [August 17, 2024, 1:59pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/1 "2024-08-17T13:59:04Z")

</div>

I’m transitioning from MATLAB to Julia and I’m wondering if there are any best-practices for optimizing for memory and computational time when writing simulations.

Specifically, I’m designing a simple verlet integration for the motion of N particles where `current_positions` is a N by 3 matrix where each row represents particle n and each column represents the x,y,z position respectively.

A simplified form looks like this:

```julia
# Initialize a 3D array to hold number_of_steps slices of positions
number_of_steps = # some integer value
N = # number of particles
positions = Array{Float64}(undef, N, 3, number_of_steps)
current_positions = # N by 3 matrix of initial positions
steps = [] 
for current_step = 1:number_of_steps
   
   # Update the positions
   current_positions = # some scheme of "current_positions"

   # Update the list of positions and steps
   steps = push!(steps, current_step)
   positions[:, :, current_step] = current_positions
end

```

What’s left is an array of `N` by 3 matricies where each “slice” is the data at a tilmestep. This works for what I need, as my intent is to have access to all the positions of all the particles for each time step. But I’m wondering if there is a way I could optimize the way I store the positions at each step?

The advantage to this is that I can slice along any dimension to get the positions at any step with simple indexing, but I assume that as my simulation gets large, the array can become very large and consume a significant amount of memory. This is because, as I understand it, an array is contiguous in memory.

The other options I could think of is making a vector of matrices `Vector{Matrix{Float64}}()` where each matrix in that vector is not continuous in memory. So things would prevent a single address in memory from become too large?

Thank you for your patience with educating me on this.

---

<div class="post-metadata">

**Author:** ![Tarny\_GG\_Channie](https://avatars.discourse-cdn.com/v4/letter/t/3bc359/32.png) [@Tarny\_GG\_Channie](https://discourse.julialang.org/u/Tarny_GG_Channie)\
**Post date:** [August 17, 2024, 2:27pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/2 "2024-08-17T14:27:44Z")

</div>

Julia is column-major, meaning that the memory is contiguous on the columns, so A[1,x] and A[2,x] would be adjacent to one another. With that, you could design the system as appropriate. Since the xyz are fixed, you could also make a struct or use a static array for it.

---

<div class="post-metadata">

**Author:** ![ckawa](https://avatars.discourse-cdn.com/v4/letter/c/8c91f0/32.png) [@ckawa](https://discourse.julialang.org/u/ckawa)\
**Post date:** [August 17, 2024, 2:53pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/3 "2024-08-17T14:53:48Z")

</div>

Thanks for the reply!

Since Julia is column-major, a 3x3 matrix `A` in Julia:

```julia
A = [x1 y1 z1; 
     x2 y2 z2; 
     x3 y3 z3]

```

would be stored in memory in column-major order:

```julia
Memory Layout: x1, x2, x3, y1, y2, y3, z1, z2, z3

```

Please correct my understanding if I’m wrong: if my interaction of the data is primarily calling the (x, y, z) position of particle “n”, I would want those values next to each other in memory. So I’d want to organize “A” as

```julia
A = [x1 x2 x3; 
     y1 y2 y3; 
     z1 z2 z3]

```

So the memory becomes:

```julia
Memory Layout: x1, y1, z1, x2, y2, z2, x3, y3, z3

```

So all the (x,y,z) are adjacent in memory?

---

<div class="post-metadata">

**Author:** ![mrufsvold](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mrufsvold/32/31600_2.png) [@mrufsvold](https://discourse.julialang.org/u/mrufsvold)\
**Post date:** [August 17, 2024, 3:02pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/4 "2024-08-17T15:02:08Z")

</div>

You are correct!

Also,

```julia
struct Position{T}
    x::T
    y::T 
    Z::T 
end

positions = Position{Float64}[
    Position(1.0,2.0,3.0),
# etc
]

```

This will have x,y,z next to each other in memory. You also get the benefit of indexing with `positions[n].y`.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 17, 2024, 3:32pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/5 "2024-08-17T15:32:05Z")

</div>

Also check out StaticArrays.jl

---

<div class="post-metadata">

**Author:** ![eldee](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@eldee](https://discourse.julialang.org/u/eldee)\
**Post date:** [August 17, 2024, 3:37pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/6 "2024-08-17T15:37:33Z")

</div>

> [@ckawa](#):
>
> I’m transitioning from MATLAB to Julia and I’m wondering if there are any best-practices for optimizing for memory and computational time when writing simulations.

Chances are you’re already aware of this, but if not, take a look at the [Performance Tips](https://docs.julialang.org/en/v1/manual/performance-tips/) in the documentation.

> [@Tarny\_GG\_Channie](#):
>
> Julia is column-major

Another way to think about this (which at least for me helps with higher-dimensional `Array`s), is that you want to arrange the axes of an `Array` in the order of decreasing increment frequency. For example, with `x = rand(100, 200, 300)`,

```julia
for i = axes(x, 1) # (1:100)
    for j = axes(x, 2) # (1:200)
        for k = axes(x, 3) # (1:300)
            x[i, j, k] *= 2
        end
    end
end

```

where the right-most index `k` changes most frequently, is slower than

```julia
for k = axes(x, 3) # (1:300)
    for j = axes(x, 2) # (1:200)
        for i = axes(x, 1) # (1:100)
            x[i, j, k] *= 2
        end
    end
end

```

where the left-most index `i` changes the fastest. On my machine the timings are `20.844 ms` versus `2.664 ms`.  
In your case that probably means you want `positions = Array{Float64}(undef, 3, N, number_of_steps)`. Or, as already highlighted by Tarny\_GG\_Channie, mrufsvold and DNF, you could use `Array{Position{Float64}}(undef, N, number_of_steps)` or `Array{SVector{3, Float64}}(undef, N, number_of_steps)`. All these variants have the same memory layout. Note that `Position` and `SVector` are not mutable. If mutability is desired, you can use `mutable struct Position` or `MVector`, though this might come with a performance penalty.

> [@ckawa](#):
>
> The other options I could think of is making a vector of matrices `Vector{Matrix{Float64}}()` where each matrix in that vector is not continuous in memory. So things would prevent a single address in memory from become too large?

You might have just misworded your thoughts here, but each `Matrix{Float64}` would still be contiguous in memory. Adjacent matrices (in terms of index in the `Vector`) need not be adjacent in memory though.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 17, 2024, 4:08pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/7 "2024-08-17T16:08:49Z")

</div>

> [@Tarny\_GG\_Channie](#):
>
> Julia is column-major

So is Matlab, actually.

(_Actually_, Julia `Array` is column major. There’s nothing intrinsically column-major about JuliaLang.)

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [August 17, 2024, 4:56pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/8 "2024-08-17T16:56:49Z")

</div>

In my own MD code, I store the particle positions (and velocities) as a `Vector{SVector{d, Float}}`, where d is the dimensionality. I think Molly.jl does something similar. I found this is the best compromise between ease of use and performance.

Then, every time i want the data saved (which definitely is not every time step), I save them to disk (not memory) in an HDF5 file with JLD2.jl. This circumvents the large memory use issue you mention, and, if you save infrequently enough, should also not impact performance.

---

<div class="post-metadata">

**Author:** ![eldee](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@eldee](https://discourse.julialang.org/u/eldee)\
**Post date:** [August 17, 2024, 5:04pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/9 "2024-08-17T17:04:44Z")

</div>

> [@ckawa](#):
>
> `steps = push!(steps, current_step)`

Note that you don’t need the `steps = `, and more importantly, that in your simplified code this boils down to `steps = collect(1:number_of_steps)`, i.e. at the end you get `steps = [1, 2, ..., number_of_steps] `. If in the actual code you do something more involved, presumably you could still preallocate `steps` using `number_of_steps`.

---

<div class="post-metadata">

**Author:** ![ckawa](https://avatars.discourse-cdn.com/v4/letter/c/8c91f0/32.png) [@ckawa](https://discourse.julialang.org/u/ckawa)\
**Post date:** [August 17, 2024, 5:08pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/10 "2024-08-17T17:08:52Z")

</div>

> [@eldee](#):
>
> You might have just misworded your thoughts here

You are correct!

> [@eldee](#):
>
> In your case that probably means you want `positions = Array{Float64}(undef, 3, N, number_of_steps)`.

> [@mrufsvold](#):
>
> This will have x,y,z next to each other in memory. You also get the benefit of indexing with `positions[n].y`.

This option would very nice to access the (x,y,z) coordinates once the loop (simulation) is finished. Thank you!

---

<div class="post-metadata">

**Author:** ![mkitti](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mkitti/32/12459_2.png) [@mkitti](https://discourse.julialang.org/u/mkitti)\
**Post date:** [August 17, 2024, 5:36pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/11 "2024-08-17T17:36:10Z")

</div>

> [@ckawa](#):
>
> `steps = []`

This might be a problem as it creates a `Any[]`. You might want

```julia
steps = Int[]

```

Notice the difference of types below.

```julia-repl
julia> steps = []
Any[]

julia> push!(steps, 1)
1-element Vector{Any}:
 1

julia> push!(steps, 2)
2-element Vector{Any}:
 1
 2

julia> typeof(steps)
Vector{Any} (alias for Array{Any, 1})

julia> steps = Int[]
Int64[]

julia> push!(steps, 1)
1-element Vector{Int64}:
 1

julia> push!(steps, 2)
2-element Vector{Int64}:
 1
 2

julia> typeof(steps)
Vector{Int64} (alias for Array{Int64, 1})

```

---

<div class="post-metadata">

**Author:** ![ckawa](https://avatars.discourse-cdn.com/v4/letter/c/8c91f0/32.png) [@ckawa](https://discourse.julialang.org/u/ckawa)\
**Post date:** [August 17, 2024, 7:25pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/12 "2024-08-17T19:25:09Z")

</div>

That’s interesting! So each SVector{d,Float} is the (x,y,z,vx,vy,vz) for each particle? and the “outer” vector in Vector{SVector{d,Float}} would be the vector of positions for some particle across all time?

---

<div class="post-metadata">

**Author:** ![IlianPihlajamaa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ilianpihlajamaa/32/28766_2.png) [@IlianPihlajamaa](https://discourse.julialang.org/u/IlianPihlajamaa)\
**Post date:** [August 17, 2024, 8:09pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/13 "2024-08-17T20:09:59Z")

</div>

So the vector `r` would contain the positions of all particles at the current point in time. For example,

```julia
r = [
SVector(x1, y2),
SVector(x2, y2),
SVector(x3, y3),
]

```

would represent three particles in two dimensions. I would have similar vectors `v` and F for the velocities and forces. Using `SVector`s gives you very convenient syntax. For example, the kinetic energy of particle 2 is `0.5m*sum(v[2].^2)`, and this would be a very performant implementation.

By choosing to save to disk instead of to some big trajectory variable I lose the ability to access directly the positions of particles at some earlier time during the simulation. Usually that is not necessary. Then, when you need to do some analysis afterwards, you read everything from disk and construct this big array containing the trajectories of all particles at all time frames.

Edit: of course, if you do prefer it, you could also have a vector of vectors of `SVector`s that you push these to and keep them in memory, but in that case you must be sure you have the required memory availability to store all of it.

---

<div class="post-metadata">

**Author:** ![ckawa](https://avatars.discourse-cdn.com/v4/letter/c/8c91f0/32.png) [@ckawa](https://discourse.julialang.org/u/ckawa)\
**Post date:** [August 17, 2024, 11:23pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/14 "2024-08-17T23:23:25Z")

</div>

Thanks for sharing! I must be missing something about `SVector`’s because I did a very simple bench test:

```julia
M = [
SVector(1, 2),
SVector(1.1, 2.2),
SVector(1.11, 2.22)
]

n_steps = 1E5

function bench_svector()
    for step in 1:n_steps
        delta_M = step / n_steps
        out = M .+ (SVector(delta_M, delta_M),)
    end
end

```

and I get something around `106.523 ms (1099490 allocations: 38.14 MiB)`

Then I try

```julia
M = [
    1.0 2.0;
    1.1 2.2;
    1.11 2.22
]
n_steps = 1E5

function bench_matrix()
    for step in 1:n_steps
        delta_M = step / n_steps
        out = M .+ delta_M
    end
end

```

and I get something like `13.652 ms (599490 allocations: 21.35 MiB)`

Does the benefit of `SVector`’s come from larger scales?

---

<div class="post-metadata">

**Author:** ![Oscar\_Smith](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/oscar_smith/32/25343_2.png) [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)\
**Post date:** [August 17, 2024, 11:28pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/15 "2024-08-17T23:28:31Z")

</div>

Don’t benchmark with non-const global variables.

```julia
function bench_svector(M, n_steps)
    for step in 1:n_steps
        delta_M = step / n_steps
        out = M .+ (SVector(delta_M, delta_M),)
    end
end
function bench_matrix(M, n_steps)
    for step in 1:n_steps
        delta_M = step / n_steps
        out = M .+ delta_M
    end
end

julia> @btime bench_svector(M, 10^5)
  2.988 ms (100000 allocations: 10.68 MiB)
julia> @btime bench_matrix(M, 10^5)
  3.254 ms (100000 allocations: 10.68 MiB)
```

---

<div class="post-metadata">

**Author:** ![ckawa](https://avatars.discourse-cdn.com/v4/letter/c/8c91f0/32.png) [@ckawa](https://discourse.julialang.org/u/ckawa)\
**Post date:** [August 17, 2024, 11:34pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/16 "2024-08-17T23:34:21Z")

</div>

Makes sense thank you!

---

<div class="post-metadata">

**Author:** ![PeterSimon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petersimon/32/25193_2.png) [@PeterSimon](https://discourse.julialang.org/u/PeterSimon)\
**Post date:** [August 17, 2024, 11:41pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/17 "2024-08-17T23:41:52Z")

</div>

> [@ckawa](#):
>
> `n_steps = 1E5`

Also, for future reference, note that the above produces a `Float64` quantity:

```julia
julia> n_steps = 1E5
100000.0

julia> typeof(n_steps)
Float64

```

You can use the following instead to create an `Int`:

```julia
julia> n_steps = 10^5
100000

```

---

<div class="post-metadata">

**Author:** ![lmiq](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lmiq/32/18314_2.png) [@lmiq](https://discourse.julialang.org/u/lmiq)\
**Post date:** [August 17, 2024, 11:42pm UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/18 "2024-08-17T23:42:16Z")

</div>

You may find this notebook useful depending on the type of simulation you are planning to do: [https://m3g.github.io/2021\_FortranCon/](https://m3g.github.io/2021_FortranCon/)

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [August 18, 2024, 6:20am UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/19 "2024-08-18T06:20:24Z")

</div>

> [@Oscar\_Smith](#):
>
> `out = M .+ delta_M`

Any advantage to using SVectors will completely drown in the enormous amount of array allocations going on here. (And I wouldn’t expect this particular example to benefit much from StaticArrays anyway, since a scalar is at least as fast as an SVector).

---

<div class="post-metadata">

**Author:** ![eldee](https://avatars.discourse-cdn.com/v4/letter/e/b5a626/32.png) [@eldee](https://discourse.julialang.org/u/eldee)\
**Post date:** [August 18, 2024, 8:29am UTC](https://discourse.julialang.org/t/memory-best-practices-in-julia-with-arrays-vs-vectors/118307/20 "2024-08-18T08:29:40Z")

</div>

To give a simple example where a `Vector{SVector{3, Float64}}` of length N is much faster than a 3 \times N `Matrix`, consider

```julia
using StaticArrays, BenchmarkTools

# Translate each point in in_v over the common translation vector x_v, 
# and store the result in the appropriate location of out_v
function bench_vec_svector!(out_v, in_v, x_v)
    for i = eachindex(out_v)
        out_v[i] = in_v[i] .+ x_v # (simply + also works, and is equally fast)
    end
end

# Same, but points are stored in columns of matrices in_m and out_m
function bench_matrix!(out_m, in_m, x_m)
    for i = axes(out_m, 2)
        out_m[:, i] .= @view(in_m[:, i]) .+ x_m
    end
end

N = 10^5

in_v = rand(SVector{3, Float64}, N)
out_v = similar(in_v)
x_v = rand(SVector{3, Float64})

@btime bench_vec_svector!(out_v, in_v, x_v);
# 86.700 μs (0 allocations: 0 bytes)

in_m = rand(3, N)
out_m = similar(in_m)
x_m = rand(3)

@btime bench_matrix!(out_m, in_m, x_m);
# 690.000 μs (0 allocations: 0 bytes)

```

Because the compiler knows the size of each `SVector` it can unroll the inner broadcasting loop, reducing overhead. The equivalent `Matrix` version would be

```julia
function bench_matrix_unrolled!(out_m, in_m, x_m)
    for i = axes(out_m, 2)
        out_m[1, i] = in_m[1, i] + x_m[1]
        out_m[2, i] = in_m[2, i] + x_m[2]
        out_m[3, i] = in_m[3, i] + x_m[3]
    end
end

@btime bench_matrix_unrolled!(out_m, in_m, x_m);
# 120.900 μs (0 allocations: 0 bytes)

```

I’m not sure why the `SVector` version is faster still, but hey, it sure does illustrate the point that `Vector{SVector}` can indeed be much more performant than `Matrix`.