Python实现QR分解遇问题：Q与R相乘结果与原矩阵不符

你的代码存在三个核心问题，导致Q与R相乘后无法还原原矩阵：

1. 旋转矩阵维度错误 #

调用GM函数时传入了原矩阵的列数n作为旋转矩阵的大小，但Givens旋转是对行进行变换，旋转矩阵的维度必须等于原矩阵的行数m。因为matQ是m×m的正交矩阵，只有同维度的旋转矩阵才能正确参与矩阵乘法运算。

2. 误用逐元素乘法代替矩阵乘法 #

你在验证结果时使用了*运算符，这是NumPy的逐元素乘法，而QR分解要求的是矩阵乘法，必须用@运算符。这是你看到结果不符的直接原因。

3. 未正确接收函数返回值 #

调用QRDecomposeWithGivens(A)时没有将返回的Q、R赋值给变量，且提前打印的NEW_q*NeW_R是未定义变量，会触发NameError。

修复后的完整代码 #

import numpy as np
import math 

def GM(size, p,q,cos, sin):
    if (p >= size) or (p < 0) or (q >= size) or (q < 0):
            raise ValueError("Invalid first index i or j")
    matrixRotation = np.identity(size)
    matrixRotation[p,p]=cos
    matrixRotation[q,p]=-1 * sin
    matrixRotation[p,q]=sin
    matrixRotation[q,q]=cos
    return matrixRotation

def QRDecomposeWithGivens(matIn, *, traceOutput = True):
    m, n = matIn.shape
    matR = matIn.copy()
    matQ = np.identity(m)
    for col in range(n):
        for row in range(m-1, col, -1):
            x1 = matR[row-1, col]
            x2 = matR[row,col]
            d = math.sqrt(x1**2 + x2**2)
            # 跳过已为0的元素，避免无效旋转
            if d == 0:
                continue
            cos = x1/d
            sin = x2/d
            # 传入行数m作为旋转矩阵的维度
            rotationMatrix = GM(m, row-1, row, cos, sin)
            matR = rotationMatrix @ matR
            matQ = matQ @ rotationMatrix.T
            if traceOutput:
                print("--- Next step ---")
                print(f"Rotation matrix for col {col}, row {row}; x{row-1}: {x1}, x{row}: {x2}")
                print(rotationMatrix)
                print("New R")
                print(matR)
                print("New Q")
                print(matQ)
                print("Q @ R (matrix multiplication):")
                print(matQ @ matR)
    return (matQ, matR)
    
A = np.array([[1, 2, 3],
              [4, 5, 6],
              [7, 8, 9]])

# 正确接收返回的Q和R矩阵
Q, R = QRDecomposeWithGivens(A)

print("\n--- Final Verification ---")
print("Original Matrix:")
print(A)
print("\nQ Orthogonal Matrix:")
print(Q)
print("\nUpper Triangular R Matrix:")
print(R)
print("\nQ @ R (should match original matrix):")
print(np.round(Q @ R, 6))  # 四舍五入消除浮点误差

补充说明 #

• 原矩阵A的秩为2，因此R的最后一行会全为0，这是正常现象。
• 加入了d==0的判断，避免对已经是0的元素做无效旋转。
• 最终验证时使用np.round消除浮点运算带来的微小误差，方便对比结果。

内容的提问来源于stack exchange，提问作者Parviz Pirizade