Composite Index updates for Numpy Matrices

Question

I am trying to update a set of particular rows and columns of a numpy matrix . Here is an example:

import numpy as np
A=np.zeros((8,8))
rows=[0, 1, 5]
columns=[2, 3]
#(What I am trying to achieve) The following does not update A
A[rows][:,columns]+=1
#while this just does
for i in rows:
    A[i][columns]+=1

the output I am expecting is:

In [1]:print(A)
Out[1]: 
    array([[ 0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  1.,  1.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
           [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.]])

Is there a way to perform a multiple-columnwise and -rowwise updates at the same time without looping?

Answer 1

rows needs to be a 'column' vector, e.g.

rows=[[0],[1],[5]]
cols=[2,3]
A[rows,cols]+=1

Sometimes a 2 stage indexing works, A[rows][:,cols], but not always. In particular it doesn't in this case where rows is not a slice. A[rows] is now a copy, so changing it does not change A.

There are various ways of indexing this block. plonser's use of product works, though I've rarely seen it used with numpy.

np.ix_ is a handy tool for doing this:

In [70]: np.ix_([0,1,5],[2,3])
Out[70]: 
(array([[0],
        [1],
        [5]]), array([[2, 3]]))

np.newaxis also turns a row vector into a column:

rows=np.array([0,1,5])
cols=np.array([2,3])
A[rows[:,None],cols]

This pairing of column and row vectors works because numpy broadcasts them to produce the (3,2) array of indices.

itertools.product produces the same set of indices, but as a list (or generator) of tuples

In [80]: list(itertools.product([0,1,5],[2,3]))
Out[80]: [(0, 2), (0, 3), (1, 2), (1, 3), (5, 2), (5, 3)]
In [84]: tuple(np.array(list(itertools.product([0,1,5],[2,3]))).T)
Out[84]: (array([0, 0, 1, 1, 5, 5]), array([2, 3, 2, 3, 2, 3]))

Answer 2

Indexing with arrays or lists of indices falls under the category of fancy indexing, and always generates a copy of the relevant portion of the array rather than giving you a view onto the original data.

Assignment or in-place modification of the result of a single fancy indexing operation works fine, (e.g. A[rows] += 1). However, modifying the result of chained indexing expressions, e.g. A[rows][:, columns], will not work if the first expression uses fancy indexing. The reason is that A[rows] generates a copy, so adding 1 to A[rows][:, columns] will not affect A.

You can avoid this problem if you do all of your indexing in one go. At some point you probably tried something like this:

In [38]: A[rows, columns] += 1
---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-38-7d2300f59d51> in <module>()
----> 1 A[rows, columns] += 1

ValueError: shape mismatch: objects cannot be broadcast to a single shape

Here numpy casts your lists of indices to arrays, then tries to broadcast them out to the same shape, which fails because they have incompatible dimensions. To make it succeed, you can make rows into an (n, 1) array and columns into a (1, m) array, where m and n are the number of individual row/column indices. That way they can be expanded along the size 1 dimensions to size (n, m), and used to index into A:

r = np.array(rows)
c = np.array(columns)
A[r[:, None], c[None, :]] += 1
print(A)
# [[ 0.  0.  1.  1.  0.  0.  0.  0.]
#  [ 0.  0.  1.  1.  0.  0.  0.  0.]
#  [ 0.  0.  0.  0.  0.  0.  0.  0.]
#  [ 0.  0.  0.  0.  0.  0.  0.  0.]
#  [ 0.  0.  0.  0.  0.  0.  0.  0.]
#  [ 0.  0.  1.  1.  0.  0.  0.  0.]
#  [ 0.  0.  0.  0.  0.  0.  0.  0.]
#  [ 0.  0.  0.  0.  0.  0.  0.  0.]]

Indexing with None is the same thing as using np.newaxis - it has the effect of inserting a new dimension of size 1.

There is actually a convenience function, np.ix_, for generating multidimensional indices from a sequence of 1D indices:

A = np.zeros((8, 8))
A[np.ix_(rows, columns)] += 1

Answer 3

You can do it by first getting all coordinates using itertool.product

import itertools
coord = np.array(list(itertools.product(rows,columns)))

and then using advanced indexing

A[tuple(coord.T)] += 1

As mentioned in other answers advanced indexing usually gives a copy. However, if treated correctly it still works as argued here.