Inconsistent dimensions when slicing array

[Spencer Lyon]

I'm a bit confused by the following:

```julia

julia> x = reshape([1:24], 2, 3, 4);

julia> x[:, :, 2]  # 1

2x3 Array{Int64,2}:

 7   9  11

 8  10  12

julia> x[:, 2, :]  # 2

2x1x4 Array{Int64,3}:

[:, :, 1] =

 3

 4

[:, :, 2] =

  9

 10

[:, :, 3] =

 15

 16

[:, :, 4] =

 21

 22

julia> x[2, :, :]  # 3

1x3x4 Array{Int64,3}:

[:, :, 1] =

 2  4  6

[:, :, 2] =

 8  10  12

[:, :, 3] =

 14  16  18

[:, :, 4] =

 20  22  24

julia> x[:, 1, 1]  # 4

2-element Array{Int64,1}:

 1

 2

julia> y = reshape([1:12], 4, 3);

julia> y[:, 1]  # 5

4-element Array{Int64,1}:

 1

 2

 3

 4

julia> y[1, :]  # 6

1x3 Array{Int64,2}:

 1  5  9

```

Below I will use the terms "earlier" (before) and "later" (after) to refer to the location of one of the slice elements relative to the others.

It appears that when slicing a multidimensional array using ` [ ]`, if our slice expression contains a `:` after it has contained an integer, we are left with singleton  dimensions in the return value (see #2, #3, and #6). On the other hand, if we never have a `:` after an integer with the slice, the dimensionality of the return value is "squeezed" (as I would expect -- see #1, #4, and #5).

Is this the desired behavior? If so, why? It often makes setting slices of an array difficult. For example we see that the slice `x[2, :, :]` has twelve elements stored in what could feasibly be thought of as a `3x4` matrix (see in # 3 that it is a `1x3x4` array). However, because of the leading singleton dimension, we cannot directly set the slice with a `3x4` matrix:

```julia

julia> x[2, :, :] = randn(3, 4)

ERROR: argument dimensions must match

 in setindex! at array.jl:595

```

For some reason, this does work, however

```julia

julia> x[2, :] = rand(Int64, 3, 4)

3x4 Array{Int64,2}:

  9210414016196894986  5590911648019920056  -2549620222663833335  8213417336355676902

 -6399277314024803139   224843245225049497   5480293538228905038  -327610178966698264

 -4819572141754207167  1796228437444309558  -8885727529087829364  4250741295695756214

```

There is more to be said here, but I think this should illustrate the concern.

[Milan Bouchet-Valat]

This is documented in the manual:
Trailing dimensions indexed with scalars are dropped. http://docs.julialang.org/en/latest/manual/arrays/#indexing

I also found this behavior quite obscure, though I guess there are reasons why it works that way.

```julia

julia> x[2, :, :] = randn(3, 4)

ERROR: argument dimensions must match

 in setindex! at array.jl:595

```

For some reason, this does work, however

```julia

julia> x[2, :] = rand(Int64, 3, 4)

3x4 Array{Int64,2}:

  9210414016196894986  5590911648019920056  -2549620222663833335  8213417336355676902

 -6399277314024803139   224843245225049497   5480293538228905038  -327610178966698264

 -4819572141754207167  1796228437444309558  -8885727529087829364  4250741295695756214

```

There is more to be said here, but I think this should illustrate the concern.

What's the meaning of x[2, :] here? I couldn't find any documentation about that.


Regards

[Spencer Lyon]

Thanks for the response, I missed that part of the docs, I still wonder why slices work that way. I have a lot of old numpy code that relies on singleton dimensions being dropped, that is now non trivial to port to julia.

As for what x[2,:] means, I am not quite sure. I know that for some reason it takes the slice x[2,:,:], then flattens it, so it becomes a 1d array, which somehow can be filled by a 2d array. Kinda weird.

[John Lynch]

If there is likely to be a significant amount of Matlab code orphaned by this decision maybe you should raise it as an issue.  If its not changed then you are, at least, likely to get an understanding of why it shouldn't be changed.

https://github.com/JuliaLang/julia

[Stefan Karpinski]

On Sun, Nov 10, 2013 at 11:32 AM, Spencer Lyon <spence...@gmail.com> wrote:

Thanks for the response, I missed that part of the docs, I still wonder why slices work that way. I have a lot of old numpy code that relies on singleton dimensions being dropped, that is now non trivial to port to julia.

 

The fact that taking a horizontal and vertical slice of a matrix gives you the same shape of object strikes me as terribly inconvenient for linear algebra work. Even more inconvenient is that taking the transpose of such a vector gives you the same vector back.

If we did drop all singleton dimensions, then our slicing behavior would be incompatible with both Python and Matlab – that's the worst of both worlds. There is unfortunately no way to make everyone happy.

As for what x[2,:] means, I am not quite sure. I know that for some reason it takes the slice x[2,:,:], then flattens it, so it becomes a 1d array, which somehow can be filled by a 2d array. Kinda weird.

This probably shouldn't work at all. Viral, Andreas, Tim, Alessandro – is this behavior intentional?

[Stefan Karpinski]

On Sun, Nov 10, 2013 at 1:56 PM, John Lynch <mjohn...@gmail.com> wrote:

If there is likely to be a significant amount of Matlab code orphaned by this decision maybe you should raise it as an issue.  If its not changed then you are, at least, likely to get an understanding of why it shouldn't be changed.

https://github.com/JuliaLang/julia

How is this likely to orphan Matlab code? This is exactly how it works in Matlab:

>> a = randn(2,3,4);

>> size(a(2,:,:))

ans =

     1     3     4

>> size(a(:,2,:))

ans =

     2     1     4

>> size(a(:,:,2))

ans =

     2     3

[Spencer Lyon]

Maybe John Lynch misread my comment. I have a lot of old python (numpy) code that doesn’t translate directly.

My concern with not dropping singleton dimension stems from the fact that this is exactly the behavior of numpy:

In [1]: import numpy as np

In [2]: x = np.reshape(range(24), (2, 3, 4))

In [3]: x[1, :, :].shape
Out[3]: (3, 4)

In [4]: x[:, 1, :].shape
Out[4]: (2, 4)

In [5]: x[:, :, 1].shape
Out[5]: (2, 3)

[Spencer Lyon]

On Sunday, November 10, 2013 5:13:50 PM UTC-5, Stefan Karpinski wrote:

On Sun, Nov 10, 2013 at 11:32 AM, Spencer Lyon <spence...@gmail.com> wrote:

Thanks for the response, I missed that part of the docs, I still wonder why slices work that way. I have a lot of old numpy code that relies on singleton dimensions being dropped, that is now non trivial to port to julia.

 

The fact that taking a horizontal and vertical slice of a matrix gives you the same shape of object strikes me as terribly inconvenient for linear algebra work. Even more inconvenient is that taking the transpose of such a vector gives you the same vector back.

If we did drop all singleton dimensions, then our slicing behavior would be incompatible with both Python and Matlab – that's the worst of both worlds. There is unfortunately no way to make everyone happy.

I don't think this is python's behavior. See my other comment. 

[John Lynch]

Yes I did.  I only realized the extent of it when I read Stefan's reply.

[Tim Holy]

I personally think the idea of silently dropping singleton dimensions is hot
death. It seems to make sense, until you think about the case where your
indexes are generated by code with conditionals, and you didn't plan ahead for
the case where your code generated a singleton rather than a range---suddenly,
all your array dimensions mean something different than you intended. The only
way it might be tolerable (maybe) is if we don't drop dimensions in cases of
indexing with something like 3:3, which might plausibly be produced by code
that generally expects to return a range (and therefore preserves that
dimension). Bottom line is that the array dimensionality, like the rest of
well-designed Julia code, at a minimum needs to be predictable from the types
of inputs.

That said, (1) the behavior that Stefan asked about seems like a bug to me,
and (2) I largely agree with issue #4048.

--Tim

[Milan Bouchet-Valat]

Le dimanche 10 novembre 2013 à 23:58 -0600, Tim Holy a écrit :

I personally think the idea of silently dropping singleton dimensions is hot 
death. It seems to make sense, until you think about the case where your 
indexes are generated by code with conditionals, and you didn't plan ahead for 
the case where your code generated a singleton rather than a range---suddenly, 
all your array dimensions mean something different than you intended. The only 
way it might be tolerable (maybe) is if we don't drop dimensions in cases of 
indexing with something like 3:3, which might plausibly be produced by code 
that generally expects to return a range (and therefore preserves that 
dimension). Bottom line is that the array dimensionality, like the rest of 
well-designed Julia code, at a minimum needs to be predictable from the types 
of inputs.

R is plagued by the same problem: you need to do mat[,1, drop=FALSE] to avoid dropping the second dimension of a matrix.

But in Julia I think that the distinction between scalars and vectors solves the problem in an elegant way: if you want to keep dimensions, index by a vector; if you want to drop them, index by a scalar. This cannot bite you in code, since if it can contain several elements, then you will always index by a vector, and dimensions will not be dropped. Ranges like 1:1 do not appear to drop dimensions, so all is fine.

What I don't get is why only trailing dimensions are dropped. This means you cannot e.g. extract easily a vector from a row of a matrix; you need to convert the resulting array. I really wish Julia would drop all singleton dimensions indexed by a scalar (the behavior used AFAICT by NumPy). It also makes the rule much less complex to handle.

This row/column asymmetry isn't very practical in applications where your matrix is conceptually symmetric, e.g. when you have two variables without a specific direction for causation. For example, if I'm studying a table giving average wage by age range and country, I may need to extract both a vector for an age range across countries (row) and for one country across age ranges (column). Getting different results in those cases would be confusing. I'm not sure whether this is a difference of needs between "statisticians" and "data analysts".


Regards

[Alan Edelman]

My opinion:  all or nothing

I don't like treating end scalars as special (just because the underlying memory model makes it a no-op)

there is no symmetry to the situation

Either drop none or drop all. 

 My vote is drop all scalars.  ---- even if it's unfamiliar to matlab users.

It would be nice if we could some day have a new user warning system

that catered to the different classes of new user -- those coming from

matlab would get a warning when a non trailing scalar was reduced in dimension.

My reasons for my opinion in part:

 I remember when Matlab introduced multidimensional arrays

and treated  all arrays as infinite dimensional with implied 1's at the end.

... and I hated using squeeze

I agree they did not have the choice we have.

I already found out the hard way that I need to write A[i,1:1] to keep my last column

a matrix ---so I already suffered the hot death anyway.  Having A[1,:] just be a matrix

confused me. Later I realized that A[i,1:1] or A[i,[1]] was kind of nice.  It reminded

me, yes I do want the matrix.

[Jonathan Malmaud]

+1 to the idea of dropping dimensions indexed by a scalar, but keeping dimensions indexed by a length-one range. That seems harmless since it's pretty hard to imagine a function use programmatically-generated indices which are sometimes ranges and sometimes scalars, depending on the value of the data - in fact, such a function would violate Tim's "well-designed Julia code" definition, since it's indexing variables change type depending on the value and not the type of its arguments. 

On Monday, November 11, 2013 7:07:17 AM UTC-5, Milan Bouchet-Valat wrote:

[Spencer Lyon]

I also agree with this. I think it makes sense for x[:, 1, :] to be 2-dimensional, but x[:, 1:1, :] should preserve all 3 dimensions.

FWIW, this is the behavior of numpy:

In [1]: import numpy as np; x = np.random.randn(2, 3, 4)

In [2]: x[0:1, :, :].shape
Out[2]: (1, 3, 4)

[Pablo Winant]

Since I've been coding in Matlab and Python, lately in Julia and Python, I'd like to add some comments here:

1/ To me, one of the strength of Numpy's implementation of arrays is the fact that it treats all dimensions symmetrically and that there is no special treatment for two-dimensional objects. Actually there has been some natural selection here, because it is possible to use a matrix object, but the documentation strongly discourages its use and favours the ndarrray object instead (cf http://wiki.scipy.org/NumPy_for_Matlab_Users). Its main advantage is the fact that it overrides the default elementwise multiplication to replace it by matrix multiplication.

This approach is associated with the implementation of multiplication as tensor reduction operations. For instance, the following code is valid in numpy:

```
M = zeros( (4,4) )
X = zeros(4)             # a 1d  array
dot(M, X)                   # a 1d array equal to M.X
dot(X,M)                    # another 1d array
X.T                            # the transpose of a 1d array is itself
```

This is true, because `dot`, i.e. the matrix multiplication, reduces the last dimension of its first argument with the first dimension of its second argument. This multiplication is very generic and works for any combination of dimensions while keeping the output dimension predictible.

In Julia or Matlab, the 1d array behaves as a hidden column 2d matrix:

``
X.T
``
 which makes X*M invalid.

[Pablo Winant]

Since I've been coding in Matlab and Python, lately in Julia and Python, I'd like to add some comments here:

1/ To me, one of the strength of Numpy's implementation of arrays is the fact that it treats all dimensions symmetrically and that there is no special treatment for two-dimensional objects. Actually there has been some natural selection here, because it is possible to use a matrix object, but the documentation strongly discourages its use and favours the ndarrray object instead (cf http://wiki.scipy.org/NumPy_for_Matlab_Users). Its main advantage is the fact that it overrides the default elementwise multiplication to replace it by matrix multiplication.

This approach is associated with the implementation of multiplication as tensor reduction operations. For instance, the following code is valid in numpy:

```
M = zeros( (4,4) )
X = zeros(4)             # a 1d  array
dot(M, X)                   # a 1d array equal to M.X
dot(X,M)                    # another 1d array
X.T                            # the transpose of a 1d array is itself
```

This is true, because `dot`, i.e. the matrix multiplication, reduces the last dimension of its first argument with the first dimension of its second argument. This multiplication is very generic and works for any combination of dimensions while keeping the output dimension predictible.

In Julia (or Matlab), the 1d array behaves as a hidden column 2d matrix:

```
X''        # not an involution since it is a 2d array even if X is a 1d array
M*X     # a 1d array in Julia, a 2d array in Matlab
X*M     # not allowed if X is a 1d array and M a matrix
```

What I would like to argue here is simply that using 1d arrays doesn't necessary imply complications in linear algebra, since for many operations, it doesn't matter anymore whether vectors are vertical or horizontal. In some sense Julia is currently since it defines 1d arrays (`zeros(5)` is a vector not a `5x5` matrix like in Matlab), but they still behave a bit like 2d arrays.

2/ the worse thing in Matlab was the squeeze operator. For instance, suppose you have a four-dimensional array `A` and you want to get a three dimensional array, by extracting `A[i,:,:,:]`. Then you would possibly want to do `squeeze(A[i,:,:,:])` but this would lead to impredictible result. If the dimension of `A` was `2x1x1x2` you would squeeze two dimensions instead of one and get a 2d matrix instead.
Apparently, the squeeze operator in Julia behaves differently since you need to explicitly tell which dimensions you want to squeeze. But how would you replicate numpy's behaviour with it ?

```
A = zeros( 2,1,1,2 )
inds = (1,:,1,:)                                   # I would like to extract a 2d matrix A[1,:,1,:]
squeeze( A[inds...], dims=??? )    #
```
the answer is probably not very complex at all, but if it's not the default, I'd like to know which would be a simple reliable way to do it.

3/ I have never seen any disadvantage with the way numpy differenciates between scalar and range indexing: it makes it easy to guess the number of dimensions of the output. Maybe there could even be some special syntax to indicate that one doesn't want to drop a dimension, like `A[1,1-]` which would translate in `A[1,1:1]`

4/ Somebody asked about the meaning of `A[1,:]` when `A` is a three dimensional array. I would say the two operations would be expected to be associative and the result to be `A[1][:]`. In numpy that gives you the one multidimensional row (the second), or equivalently `A[1,:,:]`.  The meaning of `A[1]` is itself easy to remember, if one think of the numpy as a list of list of lists. This is also consistent with the default data ordering: if A has C-based ordering, then `A[1], or `A[1,:,:]` is a contiguous 2d array. I frankly don't know how this would translate for Julia, since it has apparently opted for Fortran order instead.

Le lundi 11 novembre 2013 17:15:46 UTC+1, Pablo Winant a écrit :

Since I've been coding in Matlab and Python, lately in Julia and Python, I'd like to add some comments here:

1/ To me, one of the strength of Numpy's implementation of arrays is the fact that it treats all dimensions symmetrically and that there is no special treatment for two-dimensional objects. Actually there has been some natural selection here, because it is possible to use a matrix object, but the documentation strongly discourages its use and favours the ndarrray object instead (cf http://wiki.scipy.org/NumPy_for_Matlab_Users). Its main advantage is the fact that it overrides the default elementwise multiplication to replace it by matrix multiplication.

This approach is associated with the implementation of multiplication as tensor reduction operations. For instance, the following code is valid in numpy:

```
M = zeros( (4,4) )
X = zeros(4)             # a 1d  array
dot(M, X)                   # a 1d array equal to M.X
dot(X,M)                    # another 1d array
X.T                            # the transpose of a 1d array is itself
```

This is true, because `dot`, i.e. the matrix multiplication, reduces the last dimension of its first argument with the first dimension of its second argument. This multiplication is very generic and works for any combination of dimensions while keeping the output dimension predictible.

In Julia or Matlab, the 1d array behaves as a hidden column 2d matrix:

``
X.T
``

[Stefan Karpinski]

That is exactly what Python does:

>>> x = np.random.randn(3,3)

>>> x[2,:].shape

(3,)

>>> x[:,2].shape

(3,)

A row slice is indistinguishable from a column slice.

[Stefan Karpinski]

On Mon, Nov 11, 2013 at 12:58 AM, Tim Holy <tim....@gmail.com> wrote:

The only way it might be tolerable (maybe) is if we don't drop dimensions in cases of indexing with something like 3:3, which might plausibly be produced by code that generally expects to return a range (and therefore preserves that dimension).

This is exactly what NumPy does, so it is defensible:

>>> x = np.random.randn(3,4)

>>> x[:,2].shape

(3,)

>>> x[:,[2]].shape

(3, 1)

>>> x[:,2:3].shape

(3, 1)

>>> x[:,2:2].shape

(3, 0)

Note that 2:2 is an empty range in Python.

One concern here is that it is currently quite easy and not uncommon to write generic code in Julia that will work with either scalars or matrices. If any of those inputs are used as indices, then this breaks down is arrays are passed in instead. Maybe that's a fake concern, maybe it isn't.

[Pablo Winant]

Le lundi 11 novembre 2013 17:59:59 UTC+1, Stefan Karpinski a écrit :

That is exactly what Python does:

>>> x = np.random.randn(3,3)

>>> x[2,:].shape

(3,)

>>> x[:,2].shape

(3,)

A row slice is indistinguishable from a column slice.

What I don't get is why do you say it is so inconvenient for linear algebra. I mean, if a 1d array is neither a column nor a row, and can be used instead of both in most applications, there is no special need to transpose it or to keep track of whether it represents a vector or a linear application. Could you elaborate a bit on that ?

 

[Spencer Lyon]

I think this conversation is important and I don't want it to get lost in the stack of new posts.

Pablo raises a good point. Allowing 1d slices of an array to just be 1d arrays, regardless of how the slice was obtained (`[:, :, 1]`  vs `[1, :, :]`), should make linear algebra easier. Stefan, what did you have in mind when you thought it might make things difficult for linear algebra programming?

Another thought I have had is that although the current convention for dropping only trailing singleton dimensions is predictable (if you can remember which side of the colon leads to dropped dims ;) ) it is not consistent. I appreciate that this is the same as Matlab, which I am guessing was a design goal so Matlab users would feel comfortable, but I would argue that it doesn't make a whole lot of sense. I am a bit biased because I prefer programming in python to Matlab, but I feel that slices using the same primitives (integers and `:`) should come out with the same number of dimensions, regardless of the position of the primitives relative to one another.

[Stefan Karpinski]

See all of the discussion here: https://github.com/JuliaLang/julia/issues/4774. There have also been various discussions along similar lines in the past.