我试图通过初等列运算找到宽矩阵的逆矩阵。该方程的形式为
A * B = C我写的代码如下
import numpy as np
A = np.random.normal(size=(3, 2))
B = np.random.normal(size=(2, 3))
result = np.matmul(A, B)
print("A: \n", A)
print("B: \n", B)
print("result: \n", result)
输出是
A:
[[-0.75475352 0.57273983]
[-0.58109354 1.60719081]
[-0.46268783 -0.94800159]]
B:
[[0.50491008 0.83252516 0.0064097 ]
[1.40104051 0.26863747 0.60867192]]
result:
[[ 0.42134905 -0.47449191 0.34377291]
[ 1.95833945 -0.05202331 0.97452729]
[-1.56180438 -0.639868 -0.57998764]]
执行操作的代码如下
rows, columns = len(result), len(result[0])
result = np.matmul(A, B)
print("B = ")
print(B)
print()
print("result = ")
print(result)
print("Start...")
for idx in range(columns):
print(f"{idx} start")
divisor = result[idx][idx]
print("Divisor: ", divisor)
B[:, idx] = B[:, idx] / divisor
result = np.matmul(A, B)
print("B = ")
print(B)
print()
print("result = ")
print(result)
for jdx in range(columns):
if idx != jdx:
subtrahend = result[idx][jdx]
print("Subtrahend: ", subtrahend)
B[:, jdx] = B[:, jdx] - (subtrahend * B[:, idx])
result = np.matmul(A, B)
print("B = ")
print(B)
print()
print("result = ")
print(result)
print(f"{idx} complete...\n")
print("-----------------------------------------------------------------")
print(result)
输出是
B =
[[0.50491008 0.83252516 0.0064097 ]
[1.40104051 0.26863747 0.60867192]]
result =
[[ 0.42134905 -0.47449191 0.34377291]
[ 1.95833945 -0.05202331 0.97452729]
[-1.56180438 -0.639868 -0.57998764]]
Start...
0 start
Divisor: 0.4213490465153189
B =
[[1.19831785 0.83252516 0.0064097 ]
[3.32513036 0.26863747 0.60867192]]
result =
[[ 1. -0.47449191 0.34377291]
[ 4.64778422 -0.05202331 0.97452729]
[-3.70667595 -0.639868 -0.57998764]]
Subtrahend: -0.47449191424025916
B =
[[1.19831785 1.40111729 0.0064097 ]
[3.32513036 1.84638494 0.60867192]]
result =
[[ 1.00000000e+00 9.01785029e-18 3.43772912e-01]
[ 4.64778422e+00 2.15331272e+00 9.74527289e-01]
[-3.70667595e+00 -2.39865577e+00 -5.79987642e-01]]
Subtrahend: 0.3437729124266499
B =
[[ 1.19831785 1.40111729 -0.40553951]
[ 3.32513036 1.84638494 -0.53441783]]
result =
[[ 1.00000000e+00 9.01785029e-18 -8.94178216e-17]
[ 4.64778422e+00 2.15331272e+00 -6.23255029e-01]
[-3.70667595e+00 -2.39865577e+00 6.94267144e-01]]
0 complete...
-----------------------------------------------------------------
1 start
Divisor: 2.1533127196651907
B =
[[ 1.19831785 0.65067989 -0.40553951]
[ 3.32513036 0.85746252 -0.53441783]]
result =
[[ 1.00000000e+00 -3.50978762e-17 -8.94178216e-17]
[ 4.64778422e+00 1.00000000e+00 -6.23255029e-01]
[-3.70667595e+00 -1.11393749e+00 6.94267144e-01]]
Subtrahend: 4.64778422285792
B =
[[-1.82590189 0.65067989 -0.40553951]
[-0.66017039 0.85746252 -0.53441783]]
result =
[[ 1.00000000e+00 -3.50978762e-17 -8.94178216e-17]
[ 1.50294027e-16 1.00000000e+00 -6.23255029e-01]
[ 1.47066515e+00 -1.11393749e+00 6.94267144e-01]]
Subtrahend: -0.6232550291421382
B =
[[-1.82590189e+00 6.50679892e-01 2.22044605e-16]
[-6.60170388e-01 8.57462516e-01 1.11022302e-16]]
result =
[[ 1.00000000e+00 -3.50978762e-17 -1.04002052e-16]
[ 1.50294027e-16 1.00000000e+00 4.94053401e-17]
[ 1.47066515e+00 -1.11393749e+00 -2.07986655e-16]]
1 complete...
-----------------------------------------------------------------
2 start
Divisor: -2.0798665491526003e-16
B =
[[-1.82590189 0.65067989 -1.06759063]
[-0.66017039 0.85746252 -0.53379532]]
result =
[[ 1.00000000e+00 -3.50978762e-17 5.00041949e-01]
[ 1.50294027e-16 1.00000000e+00 -2.37540914e-01]
[ 1.47066515e+00 -1.11393749e+00 1.00000000e+00]]
Subtrahend: 1.4706651531204138
B =
[[-0.25583354 0.65067989 -1.06759063]
[ 0.12486378 0.85746252 -0.53379532]]
result =
[[ 2.64605730e-01 -3.50978762e-17 5.00041949e-01]
[ 3.49343145e-01 1.00000000e+00 -2.37540914e-01]
[ 1.94932892e-16 -1.11393749e+00 1.00000000e+00]]
Subtrahend: -1.1139374921272234
B =
[[-0.25583354 -0.53854934 -1.06759063]
[ 0.12486378 0.2628479 -0.53379532]]
result =
[[ 2.64605730e-01 5.57015475e-01 5.00041949e-01]
[ 3.49343145e-01 7.35394270e-01 -2.37540914e-01]
[ 1.94932892e-16 -3.06368148e-17 1.00000000e+00]]
2 complete...
-----------------------------------------------------------------
[[ 2.64605730e-01 5.57015475e-01 5.00041949e-01]
[ 3.49343145e-01 7.35394270e-01 -2.37540914e-01]
[ 1.94932892e-16 -3.06368148e-17 1.00000000e+00]]
如您所见，结果在第 2 步开始时就偏离了轨道，并且不一致的情况会在以下子步骤中级联。
一些探索让我进入了 sympy 库。我尝试使用 sympy 中的一些内置函数编写上述代码。但结果是相似的。
from sympy import *
import numpy as np
init_printing(use_unicode=True)
A = np.random.normal(size=(3, 2))
B = np.random.normal(size=(2, 3))
result = np.matmul(A, B)
sympy_A = Matrix(A)
sympy_B = Matrix(B)
sympy_result = sympy_A * sympy_B
print("sympy_A = ")
pprint(sympy_A)
print()
print("sympy_B = ")
pprint(sympy_B)
print()
print("sympy_result = ")
pprint(sympy_result)
输出是
sympy_A =
⎡0.805194469349538 -0.652924642433639⎤
⎢ ⎥
⎢-0.958576236163622 -0.434740834956345⎥
⎢ ⎥
⎣0.111691146842653 -1.07586944212184 ⎦
sympy_B =
⎡ 0.29518339586785 -0.921434362266118 1.48986246464666 ⎤
⎢ ⎥
⎣-0.0431486161372306 -0.0913003438306285 -1.53001914838522⎦
sympy_result =
⎡0.265852832559516 -0.682321608015616 2.198616222001 ⎤
⎢ ⎥
⎢-0.264197323182289 0.922957070561612 -0.762984951694262⎥
⎢ ⎥
⎣0.0793916495852806 -0.00468881067914134 1.81250529492314 ⎦
以下是使用 sympy 内置函数的列操作
columns, rows = sympy_result.shape
print("Start...")
for idx in range(columns):
print(f"{idx} start")
divisor = sympy_result[idx, idx]
print("Divisor: ", divisor)
k = 1 / divisor
sympy_B = sympy_B.elementary_col_op(op='n->kn', col=idx, k=k, col1=None, col2=None)
# sympy_result = sympy_result.elementary_col_op(op='n->kn', col=idx, k=k, col1=None, col2=None)
sympy_result = sympy_A * sympy_B
print("sympy_B = ")
pprint(sympy_B)
print()
print("sympy_result = ")
pprint(sympy_result)
print()
for jdx in range(columns):
if idx != jdx:
subtrahend = sympy_result[idx, jdx]
print("Subtrahend: ", subtrahend)
k = -subtrahend
sympy_B = sympy_B.elementary_col_op(op='n->n+km', col=jdx, k=k, col1=jdx, col2=idx)
# sympy_result = sympy_result.elementary_col_op(op='n->n+km', col=jdx, k=k, col1=jdx, col2=idx)
sympy_result = sympy_A * sympy_B
print("sympy_B = ")
pprint(sympy_B)
print()
print("sympy_result = ")
pprint(sympy_result)
print()
print(f"{idx} complete...\n")
print("-----------------------------------------------------------------")
输出是
Start...
0 start
Divisor: 0.265852832559516
sympy_B =
⎡ 1.11032631484853 -0.921434362266118 1.48986246464666 ⎤
⎢ ⎥
⎣-0.162302638350002 -0.0913003438306285 -1.53001914838522⎦
sympy_result =
⎡ 1.0 -0.682321608015616 2.198616222001 ⎤
⎢ ⎥
⎢-0.99377283528903 0.922957070561612 -0.762984951694262⎥
⎢ ⎥
⎣0.298630068451527 -0.00468881067914134 1.81250529492314 ⎦
Subtrahend: -0.682321608015616
sympy_B =
⎡ 1.11032631484853 -0.163834725696617 1.48986246464666 ⎤
⎢ ⎥
⎣-0.162302638350002 -0.202042941014779 -1.53001914838522⎦
sympy_result =
⎡ 1.0 2.77555756156289e-17 2.198616222001 ⎤
⎢ ⎥
⎢-0.99377283528903 0.244884391584963 -0.762984951694262⎥
⎢ ⎥
⎣0.298630068451527 0.199072937828518 1.81250529492314 ⎦
Subtrahend: 2.19861622200100
sympy_B =
⎡ 1.11032631484853 -0.163834725696617 -0.951318982893909⎤
⎢ ⎥
⎣-0.162302638350002 -0.202042941014779 -1.17317793483535 ⎦
sympy_result =
⎡ 1.0 2.77555756156289e-17 -1.11022302462516e-16⎤
⎢ ⎥
⎢-0.99377283528903 0.244884391584963 1.42194012495613 ⎥
⎢ ⎥
⎣0.298630068451527 0.199072937828518 1.15593238204835 ⎦
0 complete...
-----------------------------------------------------------------
1 start
Divisor: 0.244884391584963
sympy_B =
⎡ 1.11032631484853 -0.669028861481253 -0.951318982893909⎤
⎢ ⎥
⎣-0.162302638350002 -0.825054384671467 -1.17317793483535 ⎦
sympy_result =
⎡ 1.0 1.11022302462516e-16 -1.11022302462516e-16⎤
⎢ ⎥
⎢-0.99377283528903 1.0 1.42194012495613 ⎥
⎢ ⎥
⎣0.298630068451527 0.812926199746991 1.15593238204835 ⎦
Subtrahend: -0.993772835289030
sympy_B =
⎡0.445463606284112 -0.669028861481253 -0.951318982893909⎤
⎢ ⎥
⎣-0.982219273472611 -0.825054384671467 -1.17317793483535 ⎦
sympy_result =
⎡ 1.0 1.11022302462516e-16 -1.11022302462516e-16⎤
⎢ ⎥
⎢5.55111512312578e-17 1.0 1.42194012495613 ⎥
⎢ ⎥
⎣ 1.10649404285483 0.812926199746991 1.15593238204835 ⎦
Subtrahend: 1.42194012495613
sympy_B =
⎡0.445463606284112 -0.669028861481253 -2.22044604925031e-16⎤
⎢ ⎥
⎣-0.982219273472611 -0.825054384671467 2.22044604925031e-16 ⎦
sympy_result =
⎡ 1.0 1.11022302462516e-16 -3.23767482109533e-16⎤
⎢ ⎥
⎢5.55111512312578e-17 1.0 1.16314824706815e-16 ⎥
⎢ ⎥
⎣ 1.10649404285483 0.812926199746991 -2.63691421801158e-16⎦
1 complete...
-----------------------------------------------------------------
2 start
Divisor: -2.63691421801158e-16
sympy_B =
⎡0.445463606284112 -0.669028861481253 0.84206229921453 ⎤
⎢ ⎥
⎣-0.982219273472611 -0.825054384671467 -0.84206229921453⎦
sympy_result =
⎡ 1.0 1.11022302462516e-16 1.22782713179679 ⎤
⎢ ⎥
⎢5.55111512312578e-17 1.0 -0.441102042350565⎥
⎢ ⎥
⎣ 1.10649404285483 0.812926199746991 1.0 ⎦
Subtrahend: 1.10649404285483
sympy_B =
⎡-0.486273311509407 -0.669028861481253 0.84206229921453 ⎤
⎢ ⎥
⎣-0.0504823556790918 -0.825054384671467 -0.84206229921453⎦
sympy_result =
⎡ -0.358583406988681 1.11022302462516e-16 1.22782713179679 ⎤
⎢ ⎥
⎢ 0.488076782152 1.0 -0.441102042350565⎥
⎢ ⎥
⎣1.31838984174237e-16 0.812926199746991 1.0 ⎦
Subtrahend: 0.812926199746991
sympy_B =
⎡-0.486273311509407 -1.35356336633193 0.84206229921453 ⎤
⎢ ⎥
⎣-0.0504823556790918 -0.140519879820786 -0.84206229921453⎦
sympy_result =
⎡ -0.358583406988681 -0.998132844197812 1.22782713179679 ⎤
⎢ ⎥
⎢ 0.488076782152 1.35858340698868 -0.441102042350565⎥
⎢ ⎥
⎣1.31838984174237e-16 8.32667268468867e-17 1.0 ⎦
2 complete...
-----------------------------------------------------------------
我知道不能像我那样使用 sympy 库。
任何有效使用 sympy 库的方法，例如矩阵的符号创建以及对这些符号矩阵执行列运算，然后替换实际值，也将非常有帮助。
重要的是，我不想使用 numpy 库已经提供的内置逆或伪逆方法，因为我打算出于研究目的修改列操作。
您可以通过设置精度然后通过fractions.Fraction 将所有数字转换为有理数来解决浮点错误。
但是，我认为不可能反转 2x3 矩阵。
从数学上讲，您有 3 个 A 行向量 A_0、A_1、A_2 和 3 个 B 列向量 B_0、B_1、B_2。实际上，您希望 A_i 和 B_j 的内积成为克罗内克增量。例如，这意味着 A_0 必须垂直于 B_1 和 B_2。但这对于 2-D 空间中的 3 个向量是无法完成的。
我可能是错的，所以如果你仍然想尝试一下，这里是代码。我已经确认它适用于方阵。
import numpy as np
from fractions import Fraction
precision = 2
wide_shape = (4, 4)
rows = wide_shape[-1]
A = np.random.normal(size=wide_shape[::-1])
B = np.random.normal(size=wide_shape)
A = np.array([[Fraction(int(a*10**precision), 10**precision) for a in row] for row in A])
B = np.array([[Fraction(int(b*10**precision), 10**precision) for b in row] for row in B])
result = np.matmul(A, B)
for i in range(rows):
divisor = result[i, i]
result[i] /= divisor
A[i] /= divisor
for j in range(i+1, rows):
multiplier = result[j,i]
result[j] -= multiplier * result[i]
A[j] -= multiplier * A[i]
for i in range(rows-1, -1, -1):
for j in range(i-1, -1, -1):
multiplier = result[j, i]
result[j] -= multiplier*result[i]
A[j] -= multiplier*A[i]