CODEWITHSUNDEEP
pythonnumpymatrix
mww
mww

Reputation: 101

Modifying block matrices in Python

I would like to take a matrix and modify blocks of it. For example, with a 4x4 matrix the {1,2},{1,2} block is to the top left quadrant ([0,1;4,5] below). The {4,1},{4,1} block is the top left quadrant if we rearrange the matrix so the 4th row/column is in position 1 and the 1st in position 2.

Let's made such a 4x4 matrix: 

a = np.arange(16).reshape(4, 4)
print(a)

## [[ 0  1  2  3]
##  [ 4  5  6  7]
##  [ 8  9 10 11]
##  [12 13 14 15]]

Now one way of selecting the block, where I specify which rows/columns I want beforehand, is as follows:

C=[3,0]
a[[[C[0],C[0]],[C[1],C[1]]],[[C[0],C[1]],[C[0],C[1]]]]

## array([[15, 12],
##        [ 3,  0]])

Here's another way:

a[C,:][:,C]

## array([[15, 12],
##        [ 3,  0]])

Yet, if I have a 2x2 array, call it b, setting

a[C,:][:,C]=b

doesn't work but

a[[[C[0],C[0]],[C[1],C[1]]],[[C[0],C[1]],[C[0],C[1]]]]=b

does.

Why is this? And is this second way the most efficient possible? Thanks!

Upvotes: 0

Views: 812

Answers (1)

hpaulj
hpaulj

Reputation: 231510

The relevant section from the numpy docs is http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html#purely-integer-array-indexing Advanced array indexing.

Adapting that example to your case:

In [213]: rows=np.array([[C[0],C[0]],[C[1],C[1]]])
In [214]: cols=np.array([[C[0],C[1]],[C[0],C[1]]])

In [215]: rows
array([[3, 3],
       [0, 0]])

In [216]: cols
array([[3, 0],
       [3, 0]])

In [217]: a[rows,cols]
array([[15, 12],
       [ 3,  0]])

due to broadcasting, you don't need to repeat duplicate indices, thus:

a[[[3],[0]],[3,0]]

does just fine. np.ix_ is a convenience function to produce just such a pair:

np.ix_(C,C) 
(array([[3],
        [0]]), 
 array([[3, 0]]))

thus a short answer is:

a[np.ix_(C,C)]

A related function is meshgrid, which constructs full indexing arrays:

a[np.meshgrid(C,C,indexing='ij')]

np.meshgrid(C,C,indexing='ij') is the same as your [rows, cols]. See the functions doc for the significance of the 'ij' parameter.

np.meshgrid(C,C,indexing='ij',sparse=True) produces the same pair of arrays as np.ix_.

I don't think there's a serious difference in computational speed. Obviously some require less typing on your part.

a[:,C][C,:] works for viewing values, but not for modifying them. The details have to do with which actions make views and which make copies. The simple answer is, use only one layer of indexing if you want to modify values.

The indexing documentation:

Thus, x[ind1,...,ind2,:] acts like x[ind1][...,ind2,:] under basic slicing.

Thus a[1][3] += 7 works. But the doc also warns

Warning The above is not true for advanced indexing.

Upvotes: 1

Related Questions