python 中的反转矩阵略有偏差
Inverting matrix in python slightly off
我正在尝试使用此 http://www.irma-international.org/viewtitle/41011/ 算法来反转 nxn 矩阵。
我运行这个矩阵上的函数
[[1.0, -0.5],
[-0.4444444444444444, 1.0]]
得到输出
[[ 1.36734694, 0.64285714]
[ 0.57142857, 1.28571429]]
正确的输出应该是
[[ 1.28571429, 0.64285714]
[ 0.57142857, 1.28571429]]
我的函数:
def inverse(m):
n = len(m)
P = -1
D = 1
mI = m
while True:
P += 1
if m[P][P] == 0:
raise Exception("Not Invertible")
else:
D = D * m[P][P]
for j in range(n):
if j != P:
mI[P][j] = m[P][j] / m[P][P]
for i in range(n):
if i != P:
mI[i][P] = -m[i][P] / m[P][P]
for i in range(n):
for j in range(n):
if i != P and j != P:
mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
mI[P][P] = 1 / m[P][P]
if P == n - 1: # All elements have been looped through
break
return mI
我哪里出错了?
https://repl.it/repls/PowerfulOriginalOpensoundsystem
输出
inverse: [[Decimal('1.285714285714285693893862813'),
Decimal('0.6428571428571428469469314065')],
[Decimal('0.5714285714285713877877256260'),
Decimal('1.285714285714285693893862813')]]
numpy: [[ 1.28571429 0.64285714] [ 0.57142857 1.28571429]]
from decimal import Decimal
import numpy as np
def inverse(m):
m = [[Decimal(n) for n in a] for a in m]
n = len(m)
P = -1
D = Decimal(1)
mI = [[Decimal(0) for n in a] for a in m]
while True:
P += 1
if m[P][P] == 0:
raise Exception("Not Invertible")
else:
D = D * m[P][P]
for j in range(n):
if j != P:
mI[P][j] = m[P][j] / m[P][P]
for i in range(n):
if i != P:
mI[i][P] = -m[i][P] / m[P][P]
for i in range(n):
for j in range(n):
if i != P and j != P:
mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
mI[P][P] = 1 / m[P][P]
m = [[Decimal(n) for n in a] for a in mI]
mI = [[Decimal(0) for n in a] for a in m]
if P == n - 1: # All elements have been looped through
break
return m
m = [[1.0, -0.5],
[-0.4444444444444444, 1.0]]
print(inverse(m))
print(np.linalg.inv(np.array(m)))
我的思考过程:
起初,我以为您可能存在潜在的浮点舍入错误。事实证明这不是真的。这就是十进制爵士乐的用途。
你的错误在这里
mI = m # this just creates a pointer that points to the SAME list as m
这里
for i in range(n):
for j in range(n):
if i != P and j != P:
mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
mI[P][P] = 1 / m[P][P]
# you are not copying mI to m for the next iteration
# you are also not zeroing mI
if P == n - 1: # All elements have been looped through
break
return mI
按照算法,每次迭代都会创建一个新的a'矩阵,不会继续修改同一个旧的a。我推断这意味着在循环不变量中,a 变为 a'。适用于您的测试用例,事实证明是正确的。
我正在尝试使用此 http://www.irma-international.org/viewtitle/41011/ 算法来反转 nxn 矩阵。
我运行这个矩阵上的函数
[[1.0, -0.5],
[-0.4444444444444444, 1.0]]
得到输出
[[ 1.36734694, 0.64285714]
[ 0.57142857, 1.28571429]]
正确的输出应该是
[[ 1.28571429, 0.64285714]
[ 0.57142857, 1.28571429]]
我的函数:
def inverse(m):
n = len(m)
P = -1
D = 1
mI = m
while True:
P += 1
if m[P][P] == 0:
raise Exception("Not Invertible")
else:
D = D * m[P][P]
for j in range(n):
if j != P:
mI[P][j] = m[P][j] / m[P][P]
for i in range(n):
if i != P:
mI[i][P] = -m[i][P] / m[P][P]
for i in range(n):
for j in range(n):
if i != P and j != P:
mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
mI[P][P] = 1 / m[P][P]
if P == n - 1: # All elements have been looped through
break
return mI
我哪里出错了?
https://repl.it/repls/PowerfulOriginalOpensoundsystem
输出
inverse: [[Decimal('1.285714285714285693893862813'), Decimal('0.6428571428571428469469314065')], [Decimal('0.5714285714285713877877256260'), Decimal('1.285714285714285693893862813')]] numpy: [[ 1.28571429 0.64285714] [ 0.57142857 1.28571429]]
from decimal import Decimal
import numpy as np
def inverse(m):
m = [[Decimal(n) for n in a] for a in m]
n = len(m)
P = -1
D = Decimal(1)
mI = [[Decimal(0) for n in a] for a in m]
while True:
P += 1
if m[P][P] == 0:
raise Exception("Not Invertible")
else:
D = D * m[P][P]
for j in range(n):
if j != P:
mI[P][j] = m[P][j] / m[P][P]
for i in range(n):
if i != P:
mI[i][P] = -m[i][P] / m[P][P]
for i in range(n):
for j in range(n):
if i != P and j != P:
mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
mI[P][P] = 1 / m[P][P]
m = [[Decimal(n) for n in a] for a in mI]
mI = [[Decimal(0) for n in a] for a in m]
if P == n - 1: # All elements have been looped through
break
return m
m = [[1.0, -0.5],
[-0.4444444444444444, 1.0]]
print(inverse(m))
print(np.linalg.inv(np.array(m)))
我的思考过程:
起初,我以为您可能存在潜在的浮点舍入错误。事实证明这不是真的。这就是十进制爵士乐的用途。
你的错误在这里
mI = m # this just creates a pointer that points to the SAME list as m
这里
for i in range(n):
for j in range(n):
if i != P and j != P:
mI[i][j] = m[i][j] + (m[P][j] * mI[i][P])
mI[P][P] = 1 / m[P][P]
# you are not copying mI to m for the next iteration
# you are also not zeroing mI
if P == n - 1: # All elements have been looped through
break
return mI
按照算法,每次迭代都会创建一个新的a'矩阵,不会继续修改同一个旧的a。我推断这意味着在循环不变量中,a 变为 a'。适用于您的测试用例,事实证明是正确的。