def QR(A): def householder(a): n = len(a) v = np.zeros(n) v[0] = np.linalg.norm(a) if a[0] < 0: v[0] = -v[0] v = v + a v = v / np.linalg.norm(v) H = np.eye(n) - 2 * np.outer(v, v) return H def qr_factorization(A): m, n = A.shape Q = np.eye(m) R = A.copy() for j in range(min(m, n)): a = R[j:, j] H = np.eye(m) H[j:, j:] = householder(a) Q = Q @ H.T R = H @ R return Q, R Q, R = qr_factorization(A) for i in range(10): for j in range(10): Q[i,j]=Q[j,i] #faire la transposition b = sp.Matrix(sp.symbols('a1:11')) b = Q@b x = sp.zeros(10, 1) for i in range(9, -1, -1): x[i] = b[i] for j in range(i+1, 10): x[i] -= R[i, j] * x[j] x[i] /= R[i, i] coeff_matrix = sp.Matrix(np.zeros((10, 10))) for i in range(10): for j in range(10): coeff_matrix[i, j] = x[i].coeff(sp.Symbol('a{}'.format(j + 1))) return(coeff_matrix) 我想让这个函数返回的是A的逆矩阵,请问要怎么修改?A是一个可逆的方阵
时间: 2024-02-28 11:56:37 浏览: 75
householder QR分解.rar_Householder算法QR_householder_householderQR分解
5星 · 资源好评率100%
你可以使用 QR 分解来求解逆矩阵。具体而言,对于一个可逆的方阵 A,我们可以通过 QR 分解得到 Q 和 R,然后通过求解线性方程组来得到 A 的逆矩阵。
以下是修改后的代码:
```python
import numpy as np
import sympy as sp
def QR(A):
def householder(a):
n = len(a)
v = np.zeros(n)
v[0] = np.linalg.norm(a)
if a[0] < 0:
v[0] = -v[0]
v = v + a
v = v / np.linalg.norm(v)
H = np.eye(n) - 2 * np.outer(v, v)
return H
def qr_factorization(A):
m, n = A.shape
Q = np.eye(m)
R = A.copy()
for j in range(min(m, n)):
a = R[j:, j]
H = np.eye(m)
H[j:, j:] = householder(a)
Q = Q @ H.T
R = H @ R
return Q, R
Q, R = qr_factorization(A)
n = A.shape[0]
A_inv = np.zeros((n,n))
for i in range(n):
# 求解线性方程组 Qy = e_i
e_i = np.zeros(n)
e_i[i] = 1
y = np.linalg.solve(Q, e_i)
# 求解线性方程组 Rx = y
x = np.zeros(n)
for i in range(n-1, -1, -1):
x[i] = y[i]
for j in range(i+1, n):
x[i] -= R[i, j] * x[j]
x[i] /= R[i, i]
A_inv[:,i] = x
return A_inv
```
可以看到,在修改后的代码中,我们首先得到了 A 的 QR 分解的结果 Q 和 R。然后,我们对于每个单位向量 $e_i$,都求解了线性方程组 $Qy = e_i$ 和 $Rx = y$,得到了 A 的逆矩阵的第 $i$ 列。最终,将这些列拼接起来就得到了 A 的逆矩阵。
阅读全文