# Best practices for removing singleton dimension

**URL:** <https://discourse.julialang.org/t/best-practices-for-removing-singleton-dimension/98435>\
**Category:** Performance\
**Tags:** question\
**Created:** [May 7, 2023, 11:01am UTC](https://discourse.julialang.org/t/best-practices-for-removing-singleton-dimension/98435 "2023-05-07T11:01:39Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![gaspardb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gaspardb/32/36536_2.png) [@gaspardb](https://discourse.julialang.org/u/gaspardb)\
**Post date:** [May 7, 2023, 11:01am UTC](https://discourse.julialang.org/t/best-practices-for-removing-singleton-dimension/98435/1 "2023-05-07T11:01:39Z")

</div>

Hi!  
I often use reduction like `sum` or `prod` by specifying dimensions with the `dims` keyword argument. I’m wondering now how to get rid of the singleton dimension. Typically, if `X` is an array of size `(1, n)`, we can do

- `selectdim(X, 1, 1)`, which produces a view,
- `X[1, :]`, which indexes and produces an array
- `dropdims(X; dims=1)`, which also produces an array

I thought they were equivalent, but timing indicates otherwise:

```julia
using LinearAlgebra: dot
X = randn(1, 10000)
Y = randn(1, 10000)
f1(X, Y) = dot(view(X, 1, :), view(Y, 1, :))
f2(X, Y) = dot(X[1, :], Y[1, :])
f3(X, Y) = dot(dropdims(X, dims=1), dropdims(Y, dims=1))
@btime f1($X, $Y)
# 2.471 μs (0 allocations: 0 bytes)
@btime f2($X, $Y)
# 24.621 μs (4 allocations: 156.34 KiB)
@btime f3($X, $Y)
# 2.591 μs (4 allocations: 160 bytes)

```

So I’ve discovered that indexing allocates, which is quite a surprise for me… I’ve been using this pattern (which is common with NumPy) extensively in my code since I started using Julia. Looking closely, I saw that it is mentioned [in stackoverflow](https://stackoverflow.com/questions/52505760/dropping-singleton-dimensions-in-julia) and that the behavior of `sum` and `prod` already sparked debates on [Julia’s repo](https://github.com/JuliaLang/julia/issues/16606). So my questions are the following:

1. Should it be expressed more clearly [in the documentation](https://docs.julialang.org/en/v1/manual/arrays/#man-array-indexing) that indexing allocates? Coming from NumPy, such manipulation with broadcasting are quite common, so maybe it is worth issuing a warning not to do that in Julia.
2. What should we use for the best performances? Views seem the best, but is it always the case?
3. Are there best practices on which dimension to reduce? For instance if I can design my algorithms to reduce over the first or the last dimension, is one going to be better?

For (1), maybe it’s just me 🙂  
For (3), I’m thinking for instance to batch dimension. I got that it is better to put it in the last so that each batch is contiguous in memory.

I came across [this issue](https://github.com/FluxML/NNlib.jl/issues/466) when porting my code to support GPU. Basically, it results in the following behaviour

```julia
using CUDA
using NNlib: batched_mul

CUDA.allowscalar(false)
X = cu(randn(1, 2, 12))
S = cu(randn(5, 2, 4))
Rindex = reshape(X[1, :, :, :], 2, 3, 4)
Rdrop = reshape(dropdims(X; dims=1), 2, 3, 4)
Rview = reshape(selectdim(X, 1, 1), 2, 3, 4)
batched_mul(S, Rindex) # Works
batched_mul(S, Rdrop) # Works
batched_mul(S, Rview) # Scalar indexing

```

Edit: on the previously mentioned issue (I was a bit put off by the volume of message to start reading it) [it is mentioned that](https://github.com/JuliaLang/julia/issues/16606#issuecomment-1336270235)

> `eachslice` is currently slow because it’s type-unstable. However, [#32310](https://github.com/JuliaLang/julia/pull/32310) fixes this, so `sum(eachslice(...))` should be just as fast as the other two methods starting with Julia 1.9.  
> 2.0 should probably remove `mapslices` and `dims` keyword arguments entirely, since they can be replaced by the more-consistent `map(f, eachslice(...))`.

Still, I’d be interested on feedback before version 1.9 is out; and question 1 & 3 still apply 🙂

---

<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:** [May 7, 2023, 12:00pm UTC](https://discourse.julialang.org/t/best-practices-for-removing-singleton-dimension/98435/2 "2023-05-07T12:00:20Z")

</div>

> [@gaspardb](#):
>
> if `X` is an array of size `(1, n)`,

In this case, I would tend to just suggest [`vec(X)`](https://docs.julialang.org/en/v1/base/arrays/#Base.vec), which is essentially free for an `Array` (it shares the same underlying data rather than making a copy).

You can also do `@views X[1, :]`, which returns an explicit view (a `SubArray`).

> [@gaspardb](#):
>
> Should it be expressed more clearly [in the documentation](https://docs.julialang.org/en/v1/manual/arrays/#man-array-indexing) that indexing allocates?

As always, if you are concerned about performance, you should read the Performance Tips chapter of the manual, which has a section [Consider using views for slices](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-views).

> [@gaspardb](#):
>
> Coming from NumPy, such manipulation with broadcasting are quite common, so maybe it is worth issuing a warning not to do that in Julia.

Just put `@views` in front of your code (e.g. on a line of code, a loop, or a whole function), and then do as much slicing as you want.

> [@gaspardb](#):
>
> Are there best practices on which dimension to reduce? For instance if I can design my algorithms to reduce over the first or the last dimension, is one going to be better?

This is also covered in the performance tips: [Access arrays in memory order, along columns](https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-column-major).

---

<div class="post-metadata">

**Author:** ![gaspardb](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/gaspardb/32/36536_2.png) [@gaspardb](https://discourse.julialang.org/u/gaspardb)\
**Post date:** [May 7, 2023, 12:17pm UTC](https://discourse.julialang.org/t/best-practices-for-removing-singleton-dimension/98435/3 "2023-05-07T12:17:47Z")

</div>

Indeed, I read this section a long time ago but did not fully grasp all the information! Thanks for the pointer.  
So I’ll use `vec` for arrays of `(1, n)`, and `selectdims` until Julia 1.9; for a reason which escapes me, it seems to work better with GPUs, as I mentioned in my example. If someone has an explanation for that btw I’d be glad to hear!

---

<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:** [May 7, 2023, 1:30pm UTC](https://discourse.julialang.org/t/best-practices-for-removing-singleton-dimension/98435/4 "2023-05-07T13:30:31Z")

</div>

one thing worth mentioning is that the singleton dimension can often be avoided in the first place. people coming from python and Matlab often run into extra singleton dimensions because they create matrices where they wanted vectors. (e.g. `rand(n,1)` instead of `rand(n)`)