描述(这是一个 hwk 问题):
我不知道从哪里开始。我计划使用拉普拉斯扩展,但我不确定如何为 nxn 矩阵实现它。任何帮助,将不胜感激。
注意:我已经有了一个为 nxn 矩阵生成随机矩阵的函数。计算的时间也不成问题。我唯一遇到的问题是如何计算行列式。
不得不删除我的班级政策的问题描述 b/c。
描述(这是一个 hwk 问题):
我不知道从哪里开始。我计划使用拉普拉斯扩展,但我不确定如何为 nxn 矩阵实现它。任何帮助,将不胜感激。
注意:我已经有了一个为 nxn 矩阵生成随机矩阵的函数。计算的时间也不成问题。我唯一遇到的问题是如何计算行列式。
不得不删除我的班级政策的问题描述 b/c。
这是用于查找矩阵行列式的 adjucate 方法的递归 python 代码。
def getMatrixMinor(m,i,j):
return [row[:j] + row[j+1:] for row in (m[:i]+m[i+1:])]
def getMatrixDeternminant(m):
#base case for 2x2 matrix
if len(m) == 2:
return m[0][0]*m[1][1]-m[0][1]*m[1][0]
determinant = 0
for c in range(len(m)):
determinant += ((-1)**c)*m[0][c]*getMatrixDeternminant(getMatrixMinor(m,0,c))
return determinant
请注意,输入是表示 nxn 矩阵的数组数组
好的,这里有一个提示。
sum
)中提琴,你有一个行列式。
另外,不要忘记,由于我们在 python 中编写列表的方式,索引会被反转。那就是如果
M = [[1, 2, 3],
[4, 5, 6],
[7, 8, 9]]
那么 m 0,1是 2 而不是 4,因为它是正常的表示法。您可以将其视为转置或使用zip
我已将Peyman Naseri用作基本思想,并希望与我的实现分享,
import numpy as np
import time
def minor(arr,i,j):
c = arr[:]
c = np.delete(c, (i),axis=0)
return [np.delete(row, (j),axis=0) for row in (c)]
def det(arr):
n = len(arr)
if n == 1 :return arr[0][0]
if n == 2 :return arr[0][0]*arr[1][1] - arr[0][1]*arr[1][0]
sum = 0
for i in range(0,n):
m = minor(arr,0,i)
sum =sum + ((-1)**i)*arr[0][i] * det(m)
return sum
matrix = np.random.randint(-5, 5, size=(10, 10)) # martix nxn with integer values in interval [-5, 5)
print(matrix)
start_time = time.time()
print("started:", start_time)
det(matrix)
end_time = time.time()
print("finished:", end_time)
print("--- %s seconds ---" % (end_time - start_time))
由于我的笔记本电脑上矩阵 10*10 的结果行列式计算大约需要 1 分钟,我知道我的代码不是最优的,但这个实现的主要原因(也许我丢失了一些东西)我只需要获得基于 Peyman Naseri解决方案的工作代码,这似乎对我来说很漂亮。
[[ 2 -1 -1 0 -2 0 4 4 3 4]
[-3 1 -1 3 0 -3 -2 0 3 -1]
[ 2 -1 -4 3 0 -2 -2 -5 3 -5]
[-2 -1 2 -2 4 -3 -5 -1 -5 3]
[ 1 -4 1 -5 -5 4 -3 -5 3 1]
[ 2 4 0 -1 -1 -5 -2 -2 -3 -5]
[ 1 4 -3 -4 -5 0 0 0 -5 -1]
[ 0 -5 -5 4 -3 -2 2 -4 2 -5]
[-3 1 -1 -4 4 -5 3 -3 -4 0]
[ 0 -2 2 -3 1 3 2 0 -1 4]]
started: 1636896388.0213237
finished: 1636896442.846928
--- 54.82560420036316 seconds ---
对于大型矩阵,不建议使用拉普拉斯方法进行行列式计算,因为它的计算量很大(递归函数)。相反,更好的方法是使用高斯消除方法将原始矩阵转换为上三角矩阵。下三角矩阵或上三角矩阵的行列式只是对角元素的乘积。
这里我们展示一个例子。
import numpy as np
from scipy import linalg
def determinant(a):
assert len(a.shape) == 2 # check if a is a two diamentional matrix
assert a.shape[0] == a.shape[1] # check if matrix is square
n = a.shape[0]
for k in range(0, n-1):
for i in range(k+1, n):
if a[i,k] != 0.0:
lam = a [i,k]/a[k,k]
a[i,k:n] = a[i,k:n] - lam*a[k,k:n]
# the matrix (a) is now in the upper triangular form
return np.prod(np.diag(a))
matrix = np.random.rand(50, 50)
print(linalg.det(matrix))
print(determinant(matrix))
结果类似于从 Scipy 行列式获得的结果!
-3032.573716363944
-3032.573716915967
然而,该函数仍然可以被表述为包括旋转。此外,通过使用 numba 库,它可以达到与 scipy 类似的效率。
def minor(array,i,j):
c = array
c = c[:i] + c[i+1:]
for k in range(0,len(c)):
c[k] = c[k][:j]+c[k][j+1:]
return c
def det(array,n):
if n == 1 :return array[0][0]
if n == 2 :return array[0][0]*array[1][1] - array[0][1]*array[1][0]
sum = 0
for i in range(0,n):
m = minor(array,0,i)
sum =sum + ((-1)**i)*array[0][i] * det(m,n-1)
return sum