用 `np.linalg.solve()` 计算 `AB⁻¹`
Computing `AB⁻¹` with `np.linalg.solve()`
我需要在 Python / Numpy 中为两个矩阵 A
和 B
计算 AB⁻¹
(当然,B
是正方形)。
我知道 np.linalg.inv()
可以让我计算 B⁻¹
,然后我可以将其与 A
相乘。
我也知道 B⁻¹A
实际上是用 np.linalg.solve()
.
计算的
受此启发,我决定用 np.linalg.solve()
重写 AB⁻¹
。
我得到了一个基于 identity (AB)ᵀ = BᵀAᵀ
的公式,它使用 np.linalg.solve()
和 .transpose()
:
np.linalg.solve(a.transpose(), b.transpose()).transpose()
这似乎在做这项工作:
import numpy as np
n, m = 4, 2
np.random.seed(0)
a = np.random.random((n, n))
b = np.random.random((m, n))
print(np.matmul(b, np.linalg.inv(a)))
# [[ 2.87169378 -0.04207382 -1.10553758 -0.83200471]
# [-1.08733434 1.00110176 0.79683577 0.67487591]]
print(np.linalg.solve(a.transpose(), b.transpose()).transpose())
# [[ 2.87169378 -0.04207382 -1.10553758 -0.83200471]
# [-1.08733434 1.00110176 0.79683577 0.67487591]]
print(np.all(np.isclose(np.matmul(b, np.linalg.inv(a)), np.linalg.solve(a.transpose(), b.transpose()).transpose())))
# True
并且对于足够大的输入也更快:
n, m = 400, 200
np.random.seed(0)
a = np.random.random((n, n))
b = np.random.random((m, n))
print(np.all(np.isclose(np.matmul(b, np.linalg.inv(a)), np.linalg.solve(a.transpose(), b.transpose()).transpose())))
# True
%timeit np.matmul(b, np.linalg.inv(a))
# 100 loops, best of 3: 13.3 ms per loop
%timeit np.linalg.solve(a.transpose(), b.transpose()).transpose()
# 100 loops, best of 3: 7.71 ms per loop
我的问题是:这个恒等式总是是否正确或我忽略了一些极端情况?
一般来说,np.linalg.solve(B, A)
等同于B<sup>-1</sup>A
。剩下的只是数学。
在所有情况下,(AB)<sup>T</sup> = B<sup>T</sup>A<sup>T</sup>
: https://math.stackexchange.com/q/1440305/295281.
对于这种情况没有必要,但对于可逆矩阵,(AB)<sup>-1</sup> = B<sup>-1</sup> A<sup>-1</sup>
: https://math.stackexchange.com/q/688339/295281.
对于一个可逆矩阵,还有(A<sup>-1</sup>)<sup>T</sup> = (A <sup>T</sup>)<sup>-1</sup>
: https://math.stackexchange.com/q/340233/295281.
由此得出 (AB<sup>-1</sup>)<sup>T</sup> = (B<sup>- 1</sup>)<sup>T</sup>A<sup>T</sup> = (B<sup>T</sup>)<sup> -1</sup>A<sup>T</sup>
。只要 B
是可逆的,在任何情况下您提出的转换都不会有问题。
我需要在 Python / Numpy 中为两个矩阵 A
和 B
计算 AB⁻¹
(当然,B
是正方形)。
我知道 np.linalg.inv()
可以让我计算 B⁻¹
,然后我可以将其与 A
相乘。
我也知道 B⁻¹A
实际上是用 np.linalg.solve()
.
受此启发,我决定用 np.linalg.solve()
重写 AB⁻¹
。
我得到了一个基于 identity (AB)ᵀ = BᵀAᵀ
的公式,它使用 np.linalg.solve()
和 .transpose()
:
np.linalg.solve(a.transpose(), b.transpose()).transpose()
这似乎在做这项工作:
import numpy as np
n, m = 4, 2
np.random.seed(0)
a = np.random.random((n, n))
b = np.random.random((m, n))
print(np.matmul(b, np.linalg.inv(a)))
# [[ 2.87169378 -0.04207382 -1.10553758 -0.83200471]
# [-1.08733434 1.00110176 0.79683577 0.67487591]]
print(np.linalg.solve(a.transpose(), b.transpose()).transpose())
# [[ 2.87169378 -0.04207382 -1.10553758 -0.83200471]
# [-1.08733434 1.00110176 0.79683577 0.67487591]]
print(np.all(np.isclose(np.matmul(b, np.linalg.inv(a)), np.linalg.solve(a.transpose(), b.transpose()).transpose())))
# True
并且对于足够大的输入也更快:
n, m = 400, 200
np.random.seed(0)
a = np.random.random((n, n))
b = np.random.random((m, n))
print(np.all(np.isclose(np.matmul(b, np.linalg.inv(a)), np.linalg.solve(a.transpose(), b.transpose()).transpose())))
# True
%timeit np.matmul(b, np.linalg.inv(a))
# 100 loops, best of 3: 13.3 ms per loop
%timeit np.linalg.solve(a.transpose(), b.transpose()).transpose()
# 100 loops, best of 3: 7.71 ms per loop
我的问题是:这个恒等式总是是否正确或我忽略了一些极端情况?
一般来说,np.linalg.solve(B, A)
等同于B<sup>-1</sup>A
。剩下的只是数学。
在所有情况下,(AB)<sup>T</sup> = B<sup>T</sup>A<sup>T</sup>
: https://math.stackexchange.com/q/1440305/295281.
对于这种情况没有必要,但对于可逆矩阵,(AB)<sup>-1</sup> = B<sup>-1</sup> A<sup>-1</sup>
: https://math.stackexchange.com/q/688339/295281.
对于一个可逆矩阵,还有(A<sup>-1</sup>)<sup>T</sup> = (A <sup>T</sup>)<sup>-1</sup>
: https://math.stackexchange.com/q/340233/295281.
由此得出 (AB<sup>-1</sup>)<sup>T</sup> = (B<sup>- 1</sup>)<sup>T</sup>A<sup>T</sup> = (B<sup>T</sup>)<sup> -1</sup>A<sup>T</sup>
。只要 B
是可逆的,在任何情况下您提出的转换都不会有问题。