Numpy Routine for Computing Matrix Minors?

Question

Hi,

I'm interested in using numpy to compute all of the minors of a given square matrix. Is there a slick way of using array slicing to do this? I'm imagining that one can rotate the columns, delete the last column, rotate the rows of the resulting matrix and delete the last row, but I haven't found anything in the numpy documentation that indicates this is possible.

(Q: Why do this? A: I have a long sequence {M_n} of fairly large matrices, roughly 1,000,000 10,000 x 10,000 matrices, and I want to compute the determinant of each matrix. Each matrix is obtained from its predecessor by changing just one coefficient. It's going to be quite a lot faster to compute the determinant of the first matrix in the sequence, and then compute the difference det(M_{n+1}) - det(M_n), which is the product of the changed coefficient and its minor.)

Answer 1

In [34]: arr=np.random.random((4,4))

In [35]: arr
Out[35]: 
array([[ 0.00750932,  0.47917318,  0.39813503,  0.11755234],
       [ 0.30330724,  0.67527229,  0.71626247,  0.22526589],
       [ 0.5821906 ,  0.2060713 ,  0.50149411,  0.0328739 ],
       [ 0.42066294,  0.88529916,  0.09179092,  0.39389844]])

This gives the minor of arr, with the 1st row and 2nd column removed:

In [36]: arr[np.array([0,2,3])[:,np.newaxis],np.array([0,1,3])]
Out[36]: 
array([[ 0.00750932,  0.47917318,  0.11755234],
       [ 0.5821906 ,  0.2060713 ,  0.0328739 ],
       [ 0.42066294,  0.88529916,  0.39389844]])

So, you could use something like this:

def minor(arr,i,j):
    # ith row, jth column removed
    return arr[np.array(range(i)+range(i+1,arr.shape[0]))[:,np.newaxis],
               np.array(range(j)+range(j+1,arr.shape[1]))]

---

Regarding how this works:

Notice the shape of the index arrays:

In [37]: np.array([0,2,3])[:,np.newaxis].shape
Out[37]: (3, 1)

In [38]: np.array([0,1,3]).shape
Out[38]: (3,)

The use of [:,np.newaxis] was simply to give the first array the shape (3,1).

Since these are numpy arrays (instead of say, slices), numpy uses so-called "fancy" indexing. The rules for fancy indexing require the shape of the two arrays to be the same, or, when they are not the same, to use broadcasting to "pump-up" the shapes so that they do match.

In this case, the second array's shape (3,) is pumped up to (1,3). But (3,1) and (1,3) do not match, so (3,1) is pumped up to (3,3) and (1,3) is pumped up to (3,3).

Ah, at last, the two numpy arrays have (after broadcasting) the same shape, (3,3).

Numpy takes arr[<array of shape (3,3)>, <array of shape (3,3)>] and returns an array of shape (not surprisingly) (3,3).

The (i,j)-th element of the returned array will be

arr[(i,j)-th element of first array, (i,j)-th element of second array]

where the first and second arrays look (conceptually) like this:

first array:     second array:
[[0 0 0],        [[0, 1, 3],
 [2 2 2],         [0, 1, 3],
 [3 3 3]]         [0, 1, 3]]