使用`np.linalg.solve()计算AB'¹`

问题描述 投票:3回答:1

我需要在Python / Numpy中为两个矩阵AB⁻¹AB当然是平方)计算B

[我知道np.linalg.inv()将允许我计算B⁻¹,然后可以将其乘以A。我也知道B⁻¹A实际上是用better计算的np.linalg.solve()

[受此启发,我决定用AB⁻¹重写np.linalg.solve()。我得到了一个基于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

我的问题是:此身份总是是否正确我忽略了一些极端情况?

python numpy linear-algebra matrix-inverse
1个回答
2
投票

[通常,np.linalg.solve(B, A)等同于np.linalg.solve(B, A)。其余的只是数学。

在所有情况下,B-1A(AB)T = BTAT

在这种情况下不是必需的,但对于可逆矩阵,https://math.stackexchange.com/q/1440305/295281(AB)-1 = B-1A-1

对于可逆矩阵,https://math.stackexchange.com/q/688339/295281(A-1)T = (AT)-1也是这种情况。

由此得出https://math.stackexchange.com/q/340233/295281。只要(AB-1)T = (B-1)TAT = (BT)-1AT是可逆的,无论如何您建议的转换都不会有问题。

© www.soinside.com 2019 - 2024. All rights reserved.