Numpy Routine for Computing Matrix Minors?

I'm interested in using numpy to compute all of the 开发者_Python百科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.)

---

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(list(range(i))+list(range(i+1,arr.shape[0])))[:,np.newaxis],
               np.array(list(range(j))+list(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]]

---

The answer provided by unutbu is already great, and to optimize the algorithm @ev-br's answer had set me up on an interesting journey.

My answer below to the question is just to make it more explict of the intent.

import numpy as np

arr = np.random.normal(0,1,(4,4))



def matrix_minor(arr, i, j):
    return np.delete(np.delete(arr,i,axis=0), j, axis=1)

# tests
arr

matrix_minor(arr, 0, 0)

matrix_minor(arr, 0, 1)

---

If you're only changing one element of a matrix at a time, you're probably better off using Sherman-Morrison type formulas, (wiki): this way, you have the complexity of N^2 instead of N^3.

---

I was thinking about this exact problem the other day and did a couple of tries and performance tests for this matter, so I'll share what I found.

Adding to the solutions provided by PaulDong and unutbu, I came up with two additional solutions. One (minor_mask()) uses fancy mask-based indexing of Numpy arrays and the other, (minor_fortran()), is a solution I came up while playing around with good ol' Fortran and slightly modified it to compile with Numba. Putting all the solutions together and performing some benchmarks:

The example code

import numpy as np
import numba


def minor_mask(A, i, j):
    """Own solution, based on NumPy fancy indexing"""
    mask = np.ones_like(A, dtype=bool)
    mask[i, :] = False
    mask[:, j] = False

    minor = A[mask].reshape(A.shape[0] - 1, A.shape[1] - 1)

    del mask

    return minor


def minor_unutbu(A, i, j):
    """Solution provided by unutbu"""
    return A[
        np.array(list(range(i)) + list(range(i + 1, A.shape[0])))[:, np.newaxis],
        np.array(list(range(j)) + list(range(j + 1, A.shape[1]))),
    ]


def minor_pauldong(A, i, j):
    """Solution by PaulDong"""
    return np.delete(np.delete(A, i, axis=0), j, axis=1)


@numba.njit()
def minor_fortran(A, i, j):
    """
    Matrix minor calculation based on a Fortran routine.
    Compiles nicely with numba.
    """

    minor = np.zeros((A.shape[0] - 1, A.shape[0] - 1))

    for m in range(A.shape[0]):
        ishift = 0
        jshift = 0

        if m > i:
            ishift = 1

        for n in range(A.shape[1]):
            if n > j:
                jshift = 1

            minor[m - ishift, n - jshift] = A[m, n]

    return minor

Performance testing

On my machine (i5 9600K, 32 GB RAM, openSUSE Leap 15.2, Python 3.8.9, Numpy 1.20.3, Numba 0.53.1, ipykernel 5.5.5), I get the following results for small and large matrices using the following code:

m_small = np.arange(1e4).reshape(100, 100)
m_large = np.arange(1e8).reshape(10000, 10000)

# First run for setup of Numba and comparing results
r1 = minor_mask(m_large, 10, 11)
r2 = minor_unutbu(m_large, 10, 11)
r3 = minor_pauldong(m_large, 10, 11)
r4 = minor_fortran(m_large, 10, 11)

print(np.all(r1 == r2))
# --> True
print(np.all(r2 == r3))
# --> True
print(np.all(r3 == r4))
# --> True

# Large matrices
%timeit minor_mask(m_large, 10, 10)
# 136 ms ± 1.95 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%timeit minor_unutbu(m_large, 10, 10)
# 247 ms ± 8.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit minor_pauldong(m_large, 10, 10)
# 217 ms ± 3.79 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
%timeit minor_fortran(m_large, 10, 10)
# 153 ms ± 1.26 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

# Small matrices
%timeit minor_mask(m_small, 10, 10)
# 12.4 µs ± 22.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit minor_unutbu(m_small, 10, 10)
# 36.7 µs ± 140 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit minor_pauldong(m_small, 10, 10)
# 14.5 µs ± 102 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
%timeit minor_fortran(m_small, 10, 10)
#10.4 µs ± 34.2 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Wrap-up

So, together, we find that the list-based approach by @unutbu does perform worst for both test cases, followed by @PaulDongs approach (although IMHO being the cleanest of all solutions).
The fancy indexing approach seems to perform consistently well for both small and large matrices, only surpassed by the compiled solution for small matrices. I anticipate the mask approach to be inefficient for very large matrices, as we need to store the boolean mask in memory to do the masking, but I did not perform any memory profiling here.

I know that the comparison is not on equal footing when including the numba solution, but I'd say it's nonetheless relevant for number crunching with Python.

---

Well, i know this post is VERY old, but I wanted to give my contribution with this question, since i didn't understood nothing of the above awnswers... Although is not that simple.

Use the prints wisely to see yourself what's going on through the code.

I hope you enjoy it!

def detMinorPrincipal(A): #function to see the minor determinants
    det = np.zeros(len(A)-1)
    for i in range(len(A)-1):
        B = np.zeros((i+1)**2).reshape(i+1,i+1) #creating a new minor matrix
        for j in range(len(B)):
            for k in range(len(B[0])):
                B[j][k] = A[j][k] #equaling the new shape with the major matrix
        #print(B)
        det[i] = determinant(B) #following to the second function
    if 0 in det: #back to this function, looking for zeros
        return False
    else:
        return True

def determinant(A): #the second funtion
    if len(A) == 1:
        return A[0][0]
    if len(A) == 2:
        return A[0][0]*A[1][1] - A[1][0]*A[0][1]
    A2 = np.hstack((A,A)) #avoiding problem with index...
    print(A2)
    positivedet_sum = 0
    for i in range(len(A2)): #this one is more 'simple'...
        positivedet = 1
        for j in range(len(A2)):
            #print(A2[j][j+i], i-j, j)
            positivedet *= A2[j][j+i]
        #print(positivedet)
        positivedet_sum += positivedet
    #print(positivedet_sum)

    n=0
    negativedet_sum = 0
    for i in [len(A2)-1 for i in range(len(A2))]: #this one is harder...
        negativedet = 1
        for j in range(len(A2)):
            #print(A2[i-j][j+n], i-j, j)
            negativedet *= A2[i-j][j+n]
        #print(negativedet)
        negativedet_sum += negativedet
        n+=1
    #print(negativedet_sum)
    print('determinant:', positivedet_sum - negativedet_sum)
    return positivedet_sum - negativedet_sum #finding the awnser

while __name__ == '__main__': #start of the program
    Tdetmin = np.array([
        [2,3,4,5],
        [1,3,5,7],
        [1,1,6,2],
        [2,4,5,2]])
    detmin = detMinorPrincipal(Tdetmin) #you can use any other square matrix!
    print('Are the minor determinants different than zero?', detmin)