Scilab中Householder法QR分解算法异常，输出R非三角矩阵求助

你的Householder QR分解代码存在几个关键错误，导致R无法成为上三角矩阵，Q的计算也完全错误，以下是问题点和修正后的代码：

核心错误点 #

• Householder向量与反射系数计算错误：alpha的推导公式错误，且对A(j,j)的修改不符合Householder向量的构造逻辑。
• 列更新符号错误：应用反射矩阵时，应该是减去投影项而非加上，这直接导致下方元素无法被消为0。
• Q矩阵计算错误：直接取R的逆作为Q是完全错误的，Q应该是所有Householder反射矩阵的乘积。
• 冗余的P参数：你的代码未实现列选主元QR，因此不需要返回置换矩阵P。

修正后的代码 #

clear;
clc;

function [Q,R] = QR_decompose(A)
    n = size(A,1);
    Q = eye(n);  // 初始化Q为单位矩阵
    for j = 1:n-1  // 最后一列不需要处理，已满足上三角结构
        // 提取当前列的子向量
        v = A(j:n, j);
        sigma = norm(v)^2;
        if sigma == 0 then
            disp('Singular matrix');
            return;
        end
        // 构造Householder向量的首元素修正项，保证数值稳定性
        s = -sign(v(1));
        v(1) = v(1) - s*sqrt(sigma);
        // 计算反射系数：2/(v^T v)
        beta = 2 / (v' * v);
        
        // 用反射矩阵左乘A的剩余列，消去下方元素
        for k = j:n
            tau = beta * (v' * A(j:n, k));
            A(j:n, k) = A(j:n, k) - tau * v;
        end
        
        // 累积Q矩阵：Q = Q * H（H为当前反射矩阵）
        for k = 1:n
            tau = beta * (v' * Q(j:n, k));
            Q(j:n, k) = Q(j:n, k) - tau * v;
        end
    end
    R = A;
endfunction

A = [116 80 98 113; 80 66 80 93; 98 80 98 114; 113 93 114 133];
[Q,R] = QR_decompose(A);

disp('The upper triangular matrix is');
disp(R);
disp('The orthogonal matrix is');
disp(Q);

// 验证正确性：Q*R 应近似等于原A，Q'*Q应近似为单位矩阵
disp('Verification: Q*R - A');
disp(Q*R - A);
disp('Verification: Q''*Q - I');
disp(Q'*Q - eye(size(A,1)));

代码说明 #

1. Householder向量构造：提取当前列的子向量v，修正首元素使其满足反射条件，避免数值计算中的精度损失。
2. 反射系数beta：对应公式2/(v^T v)，用于简化反射矩阵的乘法操作，减少重复计算。
3. 更新A矩阵：对A的剩余列应用反射变换，将下方元素消为0，最终得到上三角矩阵R。
4. 累积Q矩阵：每次将反射矩阵右乘到Q上，逐步累积得到正交矩阵Q。
5. 验证步骤：添加了验证代码，确认Q*R接近原矩阵A，且Q是正交矩阵（Q'*Q近似为单位矩阵）。

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