Write vector into matrix starting from the diagonal

Question

I have a vector length n and a mxm matrix. Usually m >> n (m is way bigger than n). I need to repeatedly write the vector into the matrix, starting at the diagonal. For example:

vector v = [v_1, v_2, v_3] with a 4x4 zero-matrix results in:

v_1,  v_2,  v_3,  0
0,    v_1,  v_2,  v_3
0,    0,    v_1,  v_2
0,    0,    0,    v_1

Since I have to do this quite often, it has to be reasonably fast. Right now I am looping over every row of the matrix in raw python and writing the vector into the required position, but this is slow.

I'll try it tomorrow and come back! Sorry, was distracted with another project at work... — Leander, Jan 24 '19 at 18:13

Tarifazo · Accepted Answer · 2019-01-22T19:02:29.253

Check numpy.eye. Does this work for you?

v = [1,2,3]
N = 5
M = 10
arr = np.sum(np.eye(N, k=i, M=10) * j for i, j in enumerate(v))
arr
>>array([[1., 2., 3., 0., 0., 0., 0., 0., 0., 0.],
   [0., 1., 2., 3., 0., 0., 0., 0., 0., 0.],
   [0., 0., 1., 2., 3., 0., 0., 0., 0., 0.],
   [0., 0., 0., 1., 2., 3., 0., 0., 0., 0.],
   [0., 0., 0., 0., 1., 2., 3., 0., 0., 0.]])

Edit (thanks to hpaulj suggestion): If your matrix is very big and with lots of 0s, you can use sparse matrices

from scipy.sparse import diags
arr = diags(v,offsets=[0,1,2],shape=(N,M))
print(arr.A)
>>array([[1., 2., 3., 0., 0., 0., 0., 0., 0., 0.],
   [0., 1., 2., 3., 0., 0., 0., 0., 0., 0.],
   [0., 0., 1., 2., 3., 0., 0., 0., 0., 0.],
   [0., 0., 0., 1., 2., 3., 0., 0., 0., 0.],
   [0., 0., 0., 0., 1., 2., 3., 0., 0., 0.]])

SInce you mention `scipy.sparse`, try `sparse.diags([1,2,3],offsets=[0,1,2],shape=(4,4),dtype=int).A` — hpaulj, Jan 22 '19 at 18:03

Divakar · Answer 2 · 2019-01-22T17:48:38.060

One idea to solve it would be to pad appropriate number of zeros on either sides and get sliding windows of length m along it. We can leverage np.lib.stride_tricks.as_strided based scikit-image's view_as_windows to get sliding windows. More info on use of as_strided based view_as_windows.

from skimage.util.shape import view_as_windows

def extend2D(a, m):
    # Create zeros padded version
    p1 = np.zeros(m-1,dtype=a.dtype)
    p2 = np.zeros(m-len(a),dtype=a.dtype)    
    b = np.concatenate((p1,a,p2))

    # Get sliding windows along it of lengths `m` and finally flip rows
    return view_as_windows(b,m)[::-1]

The output would be simply sliding windowed views into a zeros padded version of input. So, if you need the output to have its own memory space, append with .copy() to the output.

Sample runs -

In [45]: a
Out[45]: array([5, 8, 6])

In [46]: extend2D(a, m=4)
Out[46]: 
array([[5, 8, 6, 0],
       [0, 5, 8, 6],
       [0, 0, 5, 8],
       [0, 0, 0, 5]])

In [47]: extend2D(a, m=5)
Out[47]: 
array([[5, 8, 6, 0, 0],
       [0, 5, 8, 6, 0],
       [0, 0, 5, 8, 6],
       [0, 0, 0, 5, 8],
       [0, 0, 0, 0, 5]])

Optimization-I

If you want to get your hands dirty with strided-indexing by sticking to NumPy using np.lib.stride_tricks.as_strided and in the process avoid that flipping at the last step in the previous approach -

def extend2D_v2(a, m):
    p1 = np.zeros(m-1,dtype=a.dtype)
    p2 = np.zeros(m-len(a),dtype=a.dtype)    
    b = np.concatenate((p1,a,p2))
    s = b.strides[0]
    return np.lib.stride_tricks.as_strided(b[m-1:],shape=(m,m),strides=(-s,s))

Optimization-II

Optimizing further, we can initialize a zeros array and then assign the input into it -

def extend2D_v3(a, m):
    b = np.zeros(2*m-1,dtype=a.dtype)
    b[m-1:m-1+len(a)] = a
    s = b.strides[0]
    return np.lib.stride_tricks.as_strided(b[m-1:],shape=(m,m),strides=(-s,s))

Timings with n=100 and m=10000 random data array -

In [97]: np.random.seed(0)
    ...: a = np.random.randint(1,9,(100))

In [98]: %timeit extend2D(a, m=10000)
    ...: %timeit extend2D_v2(a, m=10000)
    ...: %timeit extend2D_v3(a, m=10000)
10000 loops, best of 3: 51.3 µs per loop
10000 loops, best of 3: 19.6 µs per loop
100000 loops, best of 3: 12.6 µs per loop

You can also use `np.pad` instead of the concatenation or zeros + assignment. — jdehesa, Jan 22 '19 at 17:52
@jdehesa In my experience, `np.pad` seemed slower. So, avoided that. — Divakar, Jan 22 '19 at 17:52

score 0 · Answer 3 · answered Jan 22 '19 at 17:37

Here is a similar answer to Divakar's but with NumPy only. It pads the given vector with zeros and then makes a view from that buffer:

import numpy as np

def view_into_diagonals(v, m):
    # Add zeros before and after the vector
    v_pad = np.pad(v, [(m - 1, m - len(v))], mode='constant')
    # Current stride
    s, = v_pad.strides
    # Offset from which the first row starts
    offset = s * (m - 1)
    # Make ndarray
    view = np.ndarray(shape=(m, m),
                      dtype=v_pad.dtype,
                      buffer=v_pad.data,
                      offset=offset,
                      # Each column moves one forward, each row moves one backwards
                      strides=(-s, s))
    # Probably better not write to it
    view.flags.writeable = False
    return view

print(view_into_diagonals([1, 2, 3], 6))
# [[1 2 3 0 0 0]
#  [0 1 2 3 0 0]
#  [0 0 1 2 3 0]
#  [0 0 0 1 2 3]
#  [0 0 0 0 1 2]
#  [0 0 0 0 0 1]]

Write vector into matrix starting from the diagonal

3 Answers3