composing numpy array with strides

Question

I have an array A1_ijk and want to broadcast it to A2_ijmnk, but only for the case when m=n. Otherwise A2 should be filled with zeros. Of course i can create new empty array and fill it with A1 like this:

import numpy as np
A1 = np.random.rand(100, 5, 3)
A2 = np.zeros((100, 5, 2, 2, 3))
A2[..., 0, 0, :] = A1
A2[..., 1, 1, :] = A1

I feel however that this can be done in more efficient way. I tried as_strided:

as_strided = np.lib.stride_tricks.as_strided
sz = A1.itemsize
A2 = as_strided(A1, shape=(100,5,2,2,3), strides=(5*3*sz, 3*sz, 0, 0, sz))

and also broadcast_to:

broadcast_to = np.lib.stride_tricks.broadcast_to
A2=broadcast_to(A1[...,None,None,:], (100,5,2,2,3))

Unfortunately both methods fill all m, n pairs with A1 values.

My question is if this is possible to create the array of the shape i need using strides and without actual copying data?

Would `A2.shape[2],A2.shape[3]` be `(2,2)` in your actual case too? If so, for such small dimensions, your code looks okay to me in terms of readability and performance. — Divakar, Apr 22 '16 at 19:49
This can be `(2,2)`, `(3,3)`, but also more dimensions are possible, for example `m,n,o,p`, like `(2,2,2,2)` or `(3,3,3,3)`. However i have to compose `A2` quite often. I hoped that with strides i could update only `A1` and see immediately changes in `A2`. — mrkwjc, Apr 22 '16 at 19:57
Instead of constructing `A1` contiguously and trying to make `A2` a view of `A1`, can you make `A1` a view of `A2`? — user2357112, Apr 22 '16 at 20:30
Theoretically yes, but in this case updating `A1` would update only the viewed part of `A2`, isn't it? — mrkwjc, Apr 22 '16 at 20:38
Yes. If you were to make `A1` a view of `A2`, you'd have to change how your code does modifications. Updating `A1` wouldn't be enough. — user2357112, Apr 22 '16 at 21:15
So you are trying to make a diagonal matrix - but on a specific pair of dimensions (out of more than 2). You could probably simplify the problem, at least conceptually, by reshaping and transposing, so you expand a `(n,)` array to `(n,m,m)` (and then transforming it back) — hpaulj, Apr 22 '16 at 21:15

score 3 · Answer 1 · answered Apr 22 '16 at 19:48

3

My question is if this is possible to use strides to create the array of the shape i need using strides and without actual copying data?

It's not. As a simple way to see this, consider that if your array has no zeros in it, no possible view of that array will give you the off-diagonal zeros you need.

answered Apr 22 '16 at 19:48

user2357112

260,549
28
431
505

Well, maybe also zeros could be strided? – mrkwjc Apr 22 '16 at 20:24
@mrkwjc: You'd need to somehow have an array composed of views into multiple separate buffers. NumPy doesn't do that. – user2357112 Apr 22 '16 at 20:27

hpaulj · Answer 2 · 2016-04-22T22:00:22.110

I can make it simpler coding wise with np.diag_indices. I don't about efficiency relative to a strided solution (if possible). Let's see if I can simplify my development history enough

First the indices

In [2]: np.diag_indices(2)
Out[2]: (array([0, 1]), array([0, 1]))

Simpler start; we don't need 2 dimensions at the start, those can be changed with reshape. We probably don't need an ending dimension, but I'll leave that for now:

In [3]: A1=np.arange(12).reshape(4,3)

Now build the reference solution:

In [4]: A2=np.zeros((4,2,2,3),int)    
In [5]: A2[:,0,0,:]=A1
In [6]: A2[:,1,1,:]=A1

In [7]: A2
Out[7]: 
array([[[[ 0,  1,  2],
         [ 0,  0,  0]],

        [[ 0,  0,  0],    
       [[[ 3,  4,  5],
         [ 0,  0,  0]],

        [[ 0,  0,  0],
         [ 3,  4,  5]]],


       ...

       [[[ 9, 10, 11],
         [ 0,  0,  0]],

         [ 0,  1,  2]]],

          ...
        [[ 0,  0,  0],
         [ 9, 10, 11]]]])

Alternative:

In [8]: A3=np.zeros((4,2,2,3),int)

In [9]: i,j=np.diag_indices(2)

In [10]: A3[:,i,j,:]=A1
...
ValueError: shape mismatch: value array of shape (4,3) could not be broadcast to indexing result of shape (2,4,3)

Shape mismatch the first attempt

In [12]: A2[:,i,j,:]
Out[12]: 
array([[[ 0,  1,  2],
        [ 0,  1,  2]],

       [[ 3,  4,  5],
        [ 3,  4,  5]],

       [[ 6,  7,  8],
        [ 6,  7,  8]],

       [[ 9, 10, 11],
        [ 9, 10, 11]]])

In [13]: A2[:,i,j,:].shape
Out[13]: (4, 2, 3)

We need to modify A1 so it can broadcast to the destination slot.

In [14]: A1.shape
Out[14]: (4, 3)

In [15]: A3[:,i,j,:] = A1[:,None,:]

In [16]: np.allclose(A2,A3)
Out[16]: True

A2[...,i,j,:] = A1[...,None,:] should handle your example.

An even simpler version, starts with 1d array, expending to 3d

In [21]: a1=np.arange(3)

In [22]: a3=np.zeros((2,2,3),int)

In [23]: a3[...,i,j,:]=a1[...,None,:]

In [24]: a3[i,j,:]=a1   # equivalent since a1[None,:] is automatic

In [25]: a3
Out[25]: 
array([[[0, 1, 2],
        [0, 0, 0]],

       [[0, 0, 0],
        [0, 1, 2]]])

a3 doesn't have a repeated pattern of a1 values; or does it?

In [36]: a3.flatten()
Out[36]: array([0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 1, 2])

As you discovered, it's easy to fill in all the slots with as_strides, but hard to fill just the diagonals:

In [46]: ast(a1,shape=a3.shape, strides=(0,0,4))
Out[46]: 
array([[[0, 1, 2],
        [0, 1, 2]],

       [[0, 1, 2],
        [0, 1, 2]]])

score 0 · Answer 3 · answered Apr 22 '16 at 20:54

0

Youd need to pad your data with at least order n zeros; then you could set stride[n]==-stride[m] to achieve the intended effect, and avoid allocating order n*m zeros.

But something tells me there must be a more elegant solution to your problem, if you take a bigger-picture view.

answered Apr 22 '16 at 20:54

Eelco Hoogendoorn

10,459
1
44
42

I think something like this is considered here: http://stackoverflow.com/questions/18026541/make-special-diagonal-matrix-in-numpy – mrkwjc Apr 23 '16 at 11:25
yes, the second answer is along the lines of what I'm talking about; your case shouldn't be much different, apart from the extra axes – Eelco Hoogendoorn Apr 23 '16 at 13:18

composing numpy array with strides

3 Answers3