2

I am trying to populate the off diagonals of a numpy array with the columns of a separate array. Is there a way to do this for off diagonals similar to numpy.fill_diagonal?

Suppose:

A = np.zeros((4,4))

array([[0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]])

B = np.tril(np.arange(1,17).reshape(4,4),-1)  

array([[ 0,  0,  0,  0],
       [ 5,  0,  0,  0],
       [ 9, 10,  0,  0],
       [13, 14, 15,  0]]))

Is there any way to populate the off diagonals of A with columns of B? Say I wanted the -1 off diagonal of A = B[1:4,0], resulting in the following array.

B[1:4,0] = array([ 5,  9, 13]) 

A = 
array([[0., 0., 0., 0.],
       [5., 0., 0., 0.],
       [0., 9., 0., 0.],
       [0., 0., 13, 0.]])

and so on until the final output for A was

A = 
array([[0., 0., 0., 0.],
       [5., 0., 0., 0.],
       [10, 9., 0., 0.],
       [15, 14, 13, 0.]])

As far as I can tell numpy.fill_diagonal provides a way to populate the main diagonal, but does not have a parameter for off diagonals. numpy.diag does have an off diagonal parameter to create an array, but it does not seem to allow more than one off diagonal per array. So does not permit this.

numpy.diag_indices also only returns the indices of the main diagonal, so i can accomplish this at the moment is something like this.

 row,col = np.diag_indices(A.shape[0])
 for i in range(1,4): 
     A[row[i:],col[:-i]]=np.trim_zeros(B[:,i-1])

But just wondering if there was a smarter way to go about it, i.e. a function that can populate the off diagonals directly or a vectorized approach for dealing with much larger arrays.

phntm
  • 511
  • 2
  • 11
  • 1
    All those functions have python code which you can examine and emulate. `fill_diagonal` for example does a `a.flat[:end:step] = val`. I think adding a `start` would fill an off-diagonal. Because each diagonal has a different length, there's an inherent messiness to 'vectorizing' the task. – hpaulj Sep 08 '20 at 18:44
  • great, thank you i'll give it a shot. – phntm Sep 08 '20 at 19:08

1 Answers1

1

If you are after convenience there are scipy.sparse.diags and scipy.sparse.spdiagswhich despite their name are capable of producing dense output.

With your specific input format spdiags works better:

scipy.sparse.spdiags(B.T[:-1],[1,2,3],4,4,format="array").T
# array([[ 0,  0,  0,  0],
#        [ 5,  0,  0,  0],
#        [10,  9,  0,  0],
#        [15, 14, 13,  0]])
Paul Panzer
  • 51,835
  • 3
  • 54
  • 99
  • awesome thanks this is incredibly helpful. Is there a method to update an existing array rather than create a brand new one? – phntm Sep 09 '20 at 05:10
  • 1
    @phntm I'm not aware of one. You could try and use `np.lib.stride_tricks` but that is for nerds only (will segfault if used wrong). – Paul Panzer Sep 09 '20 at 08:31
  • haha cool thanks for the clarification...i'm a plebe so i'll stick to the basics. – phntm Sep 09 '20 at 21:00
  • Sorry quick follow up if you're still around. Why is it that in a function spdiags(data,diag,n,m) where data is a numpy array and the off diagonal K is negative, the entire row data [i] is entered into the off diagonal. Whereas when K is positive like you have here, the kth element of the row is automatically removed as the diagonal is populated. – phntm Sep 12 '20 at 17:50
  • 1
    @phntm I think the difference is not cut-vs-not-cut, but rather cut-from-left-vs-cut-from-right. Try with inputs that have no zeros to get a better idea. – Paul Panzer Sep 12 '20 at 19:04