这里不要用np.linalg.inv来求R的逆矩阵

非常感谢您的提醒，我重新修改了代码。对于一个上三角矩阵 $R$，我们可以使用回代法（back-substitution）来求解它的逆矩阵。具体而言，我们从最后一行开始，依次求解 $R$ 的每一行。对于第 $i$ 行，我们可以使用以下公式来求解：

$$R_{i,i}^{-1} = \frac{1}{R_{i,i}}$$

$$R_{i,j}^{-1} = -\frac{1}{R_{i,i}}\sum_{k=i+1}^n R_{i,k}R_{k,j}^{-1}$$

其中 $j>i$，$n$ 是 $R$ 的维度。

以下是修改后的代码：

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]
    R_inv = np.zeros((n, n))
    for i in range(n-1, -1, -1):
        R_inv[i,i] = 1 / R[i,i]
        for j in range(i+1, n):
            R_inv[i,j] = - np.sum(R[i,k] * R_inv[k,j] for k in range(i+1, n)) / R[i,i]
    
    A_inv = R_inv @ Q.T
    
    return A_inv

在修改后的代码中，我们首先计算了 $A$ 的 QR 分解得到的矩阵 $Q$ 和 $R$。然后，我们使用回代法来求解 $R$ 的逆矩阵，并使用 $R^{-1}Q^T$ 来计算 $A$ 的逆矩阵，并将其返回。