计算矩阵小行列式的Numpy方法？

Question

我想用numpy来计算一个给定方阵的所有小行列式。有没有什么简单的方法可以利用数组切片来做到这一点？我在想，或许可以旋转列，删除最后一列，然后再旋转结果矩阵的行，最后删除最后一行，但我在numpy的文档中没有找到任何说明这可以实现的内容。

（问：为什么要这样做？答：我有一长串 {M_n} 的大矩阵，大约有1,000,000个 10,000 x 10,000 的矩阵，我想计算每个矩阵的行列式。每个矩阵都是通过改变前一个矩阵的一个系数得到的。计算序列中第一个矩阵的行列式会快很多，然后再计算 det(M_{n+1}) - det(M_n)，这个差值就是改变的系数和它的小行列式的乘积。）

Answer 1

如果你一次只改变矩阵中的一个元素，那么使用Sherman-Morrison类型的公式可能会更好。这样做的话，你的计算复杂度是N的平方，而不是N的立方。

Answer 2

unutbu 给出的答案已经很棒了，为了优化算法，@ev-br 的回答让我开始了一段有趣的探索之旅。

我下面的回答只是为了更清楚地表达我的意图。

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)

Answer 3

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]])

这段代码会给出数组 arr 的一个小矩阵，去掉了第一行和第二列：

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]])

所以，你可以用类似这样的方式：

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])))]

---

关于这个是怎么工作的：

注意索引数组的形状：

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,)

使用 [:,np.newaxis] 只是为了让第一个数组的形状变成 (3,1)。

因为这些是 numpy 数组（而不是切片），numpy 使用所谓的“花式”索引。花式索引的规则要求两个数组的形状要相同，或者在形状不同时，使用广播的方式来“扩展”形状，使它们能够匹配。

在这个例子中，第二个数组的形状 (3,) 被扩展成 (1,3)。但是 (3,1) 和 (1,3) 还是不匹配，所以 (3,1) 被扩展到 (3,3)，而 (1,3) 也被扩展到 (3,3)。

啊，最后，这两个 numpy 数组在经过广播后，形状都变成了 (3,3)。

numpy 会用 arr[<形状为 (3,3) 的数组>, <形状为 (3,3) 的数组>] 来返回一个形状（没意外的话）也是 (3,3) 的数组。

返回数组的 (i,j) 位置的元素将会是：

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

其中第一个和第二个数组在概念上看起来是这样的：

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