How to access lower triangle of N*M*M numpy array-CodePudding

I have a numpy array of the shape arr.shape = N,M,M.

I want to access the lower triangles for each M,M array. I tried using

arr1 = arr[:,np.tril_indices(M,-1)]
arr1 = arr[:][np.tril_indices(M,-1)]

etc, with the kernel dying in the first case, while in the second case I get an error saying that:

   ---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-23-1b36c5b12706> in <module>
----> 1 arr1 = arr[:][np.tril_indices(M,-1)]

IndexError: index 6 is out of bounds for axis 0 with size 6

Where

N=6

To clarify I want to find all the elements in the lower triangle of each M,M array(N such instances) and save the result in a new array of the shape:

arr1.shape = (N,(M*(M-1))/2)

Edit:

While np.tril(arr) works, it results in an array

arr1 = np.tril(arr)
arr1.shape

#(N,M,M)

I want the resulting array to be of the specified shape, i.e. I dont want the upper parts of the arrays

Thank you

CodePudding user response：

import numpy as np 
a = np.random.rand(2, 5, 5) 
#array([[[0.28212197, 0.29827562, 0.05151153, 0.90448236, 0.07521404],
#        [0.38938978, 0.67007919, 0.83561652, 0.5950061 , 0.73563179],
#        [0.77515285, 0.31973392, 0.91861436, 0.87386527, 0.85917542],
#        [0.12588184, 0.09173029, 0.28577701, 0.4884228 , 0.07183555],
#        [0.68656271, 0.19941039, 0.07924489, 0.15046004, 0.91011737]],
#
#       [[0.18662788, 0.45745028, 0.14557573, 0.22425571, 0.14204739],
#        [0.44502694, 0.85773626, 0.78554919, 0.07306402, 0.14608384],
#        [0.70620254, 0.81497515, 0.09397011, 0.32053184, 0.255485  ],
#        [0.50139688, 0.51539848, 0.24719375, 0.80708819, 0.39685176],
#        [0.94052069, 0.53927081, 0.39567362, 0.06065674, 0.53479994]]])

np.tril(a) 
#array([[[0.28212197, 0.        , 0.        , 0.        , 0.        ],
#        [0.38938978, 0.67007919, 0.        , 0.        , 0.        ],
#        [0.77515285, 0.31973392, 0.91861436, 0.        , 0.        ],
#        [0.12588184, 0.09173029, 0.28577701, 0.4884228 , 0.        ],
#        [0.68656271, 0.19941039, 0.07924489, 0.15046004, 0.91011737]],
#
#       [[0.18662788, 0.        , 0.        , 0.        , 0.        ],
#        [0.44502694, 0.85773626, 0.        , 0.        , 0.        ],
#        [0.70620254, 0.81497515, 0.09397011, 0.        , 0.        ],
#        [0.50139688, 0.51539848, 0.24719375, 0.80708819, 0.        ],
#        [0.94052069, 0.53927081, 0.39567362, 0.06065674, 0.53479994]]])

If you want to remove the zeros and flatten it to a (2, 15) array (note that there are 10 zeros in each lower triangle array) -

a_no_zeros = np.array([el
for mat in a_lower
for row in mat
for el in row
if el > 0
]).reshape(2, 15)

CodePudding user response：

When working with the tri... set of functions it can be useful to examine the source code. They are all python, and based on np.tri.

Make a small sample array - to illustrate and verify the answer:

In [205]: arr = np.arange(18).reshape(2,3,3)  # arange(1,19) might be better
In [206]: arr
Out[206]: 
array([[[ 0,  1,  2],
        [ 3,  4,  5],
        [ 6,  7,  8]],

       [[ 9, 10, 11],
        [12, 13, 14],
        [15, 16, 17]]])

tril sets the upper triangle values to 0. It works in this case, but application to 3d arrays is not documented.

In [207]: np.tril(arr) 
Out[207]: 
array([[[ 0,  0,  0],
        [ 3,  4,  0],
        [ 6,  7,  8]],

       [[ 9,  0,  0],
        [12, 13,  0],
        [15, 16, 17]]])

But in the code if first constructs a boolean mask from the last 2 dimensions:

In [208]: mask = np.tri(*arr.shape[-2:], dtype=bool)
In [209]: mask
Out[209]: 
array([[ True, False, False],
       [ True,  True, False],
       [ True,  True,  True]])

and uses np.where to set some values to 0. This works in the 3d case by broadcasting. mask and arr match on the last 2 dimensions, so mask can broadcast to match:

In [210]: np.where(mask, arr, 0)
Out[210]: 
array([[[ 0,  0,  0],
        [ 3,  4,  0],
        [ 6,  7,  8]],

       [[ 9,  0,  0],
        [12, 13,  0],
        [15, 16, 17]]])

Your tril_indices is just the indices of this mask:

In [217]: np.nonzero(mask)    # aka np.where
Out[217]: (array([0, 1, 1, 2, 2, 2]), array([0, 0, 1, 0, 1, 2]))
In [218]: np.tril_indices(3)
Out[218]: (array([0, 1, 1, 2, 2, 2]), array([0, 0, 1, 0, 1, 2]))

They can't be used directly to index arr:

In [220]: arr[np.tril_indices(3)].shape
Traceback (most recent call last):
  File "<ipython-input-220-e26dc1f514cc>", line 1, in <module>
    arr[np.tril_indices(3)].shape
IndexError: index 2 is out of bounds for axis 0 with size 2

In [221]: arr[:,np.tril_indices(3)].shape
Out[221]: (2, 2, 6, 3)

But unpacking the two indexing arrays:

In [222]: I,J = np.tril_indices(3)
In [223]: I,J
Out[223]: (array([0, 1, 1, 2, 2, 2]), array([0, 0, 1, 0, 1, 2]))
In [224]: arr[:,I,J]
Out[224]: 
array([[ 0,  3,  4,  6,  7,  8],
       [ 9, 12, 13, 15, 16, 17]])

The boolean mask can also be used directly:

In [226]: arr[:,mask]
Out[226]: 
array([[ 0,  3,  4,  6,  7,  8],
       [ 9, 12, 13, 15, 16, 17]])

The base np.tri works by simply doing an outer >= on indices

In [231]: m = np.greater_equal.outer(np.arange(3),np.arange(3))
In [232]: m
Out[232]: 
array([[ True, False, False],
       [ True,  True, False],
       [ True,  True,  True]])
In [234]: np.arange(3)[:,None]>=np.arange(3)
Out[234]: 
array([[ True, False, False],
       [ True,  True, False],
       [ True,  True,  True]])