How to make sure the stride for discontiguous numpy array?

Question

Here is an example when I use numpy and try to understand the strides of an array

import numpy as np
m = np.arange(6, dtype=np.uint8).reshape(2, 3)
print(m.shape, m.strides)
# m.shape: (2, 3)  m.strides:(3, 1)
m = m.transpose(1, 0)
print(m.shape, m.strides)
# m.shape:(3, 2) m.strides:(1, 3)

Till now, it's all fine...the shape and the strides are all easily to understand. And array m now is discontiguous. But when I try to expand_dims for m :

x = np.expand_dims(m, 0)
print(x.shape, x.strides)
# x.shape: (1, 3, 2)  x.strides:(3, 1, 3)
y = np.expand_dims(m, 1)
print(y.shape, y.strides)
# y.shape: (3, 1, 2)  y.strides:(1, 6, 3)
z = np.expand_dims(m, 2)
print(z.shape, z.strides)
# z.shape:(3, 2, 1) z.strides:(1, 3, 3)

I just can't understand the why the expanded dims' strides of x, y, z are 3, 6, 3?

Is there a formula or document that can help me to understand the strides for discontiguous array? Thanks a lot.

Answer 1

The problem is when you do the transpose for m, the array changes from a C_CONTIGUOUS to F_CONTIGUOUS:

m = np.arange(6, dtype=np.uint8).reshape(2, 3)
print(m.flags)

flags of m

C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : False
...

For n

n = m.transpose(1, 0)
print(n.flags)

flags of n

C_CONTIGUOUS : False
F_CONTIGUOUS : True
OWNDATA : False
...

Then we look into the array_interface :

m.__array_interface__

{'data': (1840062163472, False),
 'strides': None,
 'descr': [('', '|u1')],
 'typestr': '|u1',
 'shape': (2, 3),
  'version': 3}

and for n :

n.__array_interface__

{'data': (1840062163472, False),
 'strides': (1, 3),
 'descr': [('', '|u1')],
 'typestr': '|u1',
 'shape': (3, 2),
 'version': 3}

We find that after the transpose, the key strides is assigned by a value (1, 3). Refer to this answer : How does NumPy's transpose() method permute the axes of an array?. We know the transpose swaps the shape and stride information for each axis. In this process, the key strides is assigned.

And when you use expand_dims then calculate the strides. Here I'm not sure, I think it takes the strides in the __array_interface__, then insert a stride based on the axis you called in the expand_dims. For example :

m = np.arange(6, dtype=np.uint8).reshape(2, 3)
n = m.transpose(1, 0)
y = np.expand_dims(n, 1)
y.strides
# (1, 6, 3)

So if you don't want to change the array to F_CONTIGUOUS, need to use np.ascontiguousarray after the transpose, which gives you a good result.

m = np.arange(6, dtype=np.uint8).reshape(2, 3)
a = np.ascontiguousarray(m.transpose(1, 0))
x = np.expand_dims(a, 0)
print(x.shape, x.strides)
# (1, 3, 2) (6, 2, 1)
y = np.expand_dims(a, 1)
print(y.shape, y.strides)
# (3, 1, 2) (2, 2, 1)
z = np.expand_dims(a, 2)
print(z.shape, z.strides)
# (3, 2, 1) (2, 1, 1)

Answer 2

In [25]: m = np.arange(6, dtype=np.uint8).reshape(2, 3)

As you note the strides get reversed in the transpose:

In [26]: m.strides, m.T.strides
Out[26]: ((3, 1), (1, 3))

In [27]: m
Out[27]: 
array([[0, 1, 2],
       [3, 4, 5]], dtype=uint8)

In [28]: m.T
Out[28]: 
array([[0, 3],
       [1, 4],
       [2, 5]], dtype=uint8)

We can picture that from the display. For m, to go down rows we step by 3. For m.T to go across columns we step by 3.

expand_dims uses reshape:

In [31]: m1=m.reshape(1,2,3);m1.strides,m1
Out[31]: 
((6, 3, 1),
 array([[[0, 1, 2],
         [3, 4, 5]]], dtype=uint8))

The meaning of the 6 becomes more obvious when we do a repeat:

In [32]: m1.repeat(2,0)      # (6,3,1) strides
Out[32]: 
array([[[0, 1, 2],
        [3, 4, 5]],

       [[0, 1, 2],
        [3, 4, 5]]], dtype=uint8)

To go down one plane or block, we step by 6. (m1 is a view, but the repeat is a copy, with its own data.)

Adding a dimension in the middle, we just have to "repeat" the 3:

In [36]: m2=m.reshape(2,1,3).repeat(2,axis=1);m2.strides, m2
Out[36]: 
((6, 3, 1),
 array([[[0, 1, 2],
         [0, 1, 2]],
 
        [[3, 4, 5],
         [3, 4, 5]]], dtype=uint8))

Doing the same to the transpose doesn't add any complications.

Your x case, adding a leading dimension to the transpose (with repeat for clarity):

In [38]: m.T.reshape(1,3,2).repeat(2,axis=0)
Out[38]: 
array([[[0, 3],
        [1, 4],
        [2, 5]],

       [[0, 3],
        [1, 4],
        [2, 5]]], dtype=uint8)

The new first dimension steps by a block of 6.

While the contiguity flags make some sense when dealing with 2d array, and sometimes are relevant when passing arrays to specialized compiled code, they aren't as meaningful when looking a 3d+ arrays. It's easy to make arrays, especially views, that aren't contiguous in either sense.

In recent Indexing on ndarrays result in wrong shape, I explore what happens to strides when we apply advanced indexing to columns.

edit

Your case where you get a 6:

In [61]: np.expand_dims(m.T,1)
Out[61]: 
array([[[0, 3]],

       [[1, 4]],

       [[2, 5]]], dtype=uint8)

In [62]: _.strides    # 1st dim steps by 1, last by 3
Out[62]: (1, 6, 3)

It's harder to visualize the blocks of 6 for the middle dimension, but it's doing the same thing as with the non-transpose. The strides of the repeat don't help, because it's a copy:

In [63]: np.expand_dims(m.T,1).repeat(2,1)
Out[63]: 
array([[[0, 3],
        [0, 3]],

       [[1, 4],
        [1, 4]],

       [[2, 5],
        [2, 5]]], dtype=uint8)
In [64]: _.strides
Out[64]: (4, 2, 1)
  

We could get the same [61] by first expanding, and then transposing:

In [66]: np.expand_dims(m,0).transpose(2,0,1)
Out[66]: 
array([[[0, 3]],

       [[1, 4]],

       [[2, 5]]], dtype=uint8)   
In [67]: _.shape,_.strides
Out[67]: ((3, 1, 2), (1, 6, 3))